# H.E.S.S. observations of the 2021 periastron passage of PSR B1259-63/LS 2883

H.E.S.S. Collaboration: F. Aharonian<sup>1,2,3</sup>, F. Ait Benkhali<sup>4</sup>, J. Aschersleben<sup>5</sup>, H. Ashkar<sup>6</sup>, M. Backes<sup>7,8</sup>, V. Barbosa Martins<sup>9</sup>, R. Batzofin<sup>10</sup>, Y. Becherini<sup>11,12</sup>, D. Berge<sup>9,13</sup>, K. Bernlöhr<sup>2</sup>, M. Böttcher<sup>8</sup>, C. Boisson<sup>14</sup>, J. Bolmont<sup>15</sup>, M. de Bony de Lavergne<sup>16</sup>, J. Borowska<sup>13</sup>, M. Bouyahiaoui<sup>2</sup>, R. Brose<sup>1</sup>, A. Brown<sup>17</sup>, F. Brun<sup>16</sup>, B. Bruno<sup>18</sup>, T. Bulik<sup>19</sup>, C. Burger-Scheidlin<sup>1</sup>, S. Caroff<sup>20</sup>, S. Casanova<sup>21</sup>, J. Celic<sup>18</sup>, M. Cerruti<sup>11</sup>, T. Chand<sup>8</sup>, S. Chandra<sup>8</sup>, A. Chen<sup>22</sup>, J. Chibueze<sup>8</sup>, O. Chibueze<sup>8</sup>, G. Cotter<sup>17</sup>, J. Damascene Mbarubucyeye<sup>9</sup>, J. Devin<sup>23</sup>, J. Djuvsland<sup>2</sup>, A. Dmytriiev<sup>8</sup>, K. Egberts<sup>10</sup>, S. Einecke<sup>24</sup>, J.-P. Ernenwein<sup>25</sup>, G. Fontaine<sup>6</sup>, S. Funk<sup>18</sup>, S. Gabici<sup>11</sup>, Y. A. Gallant<sup>23</sup>, D. Glawion<sup>18</sup>, J. F. Glicenstein<sup>16</sup>, P. Goswami<sup>11</sup>, G. Grolleron<sup>15</sup>, L. Haerer<sup>2</sup>, B. Heß<sup>26</sup>, W. Hofmann<sup>2</sup>, T. L. Holch<sup>9</sup>, M. Holler<sup>27</sup>, Zhiqiu Huang<sup>2</sup>, M. Jamrozy<sup>28</sup>, F. Jankowsky<sup>4</sup>, V. Joshi<sup>18</sup>, I. Jung-Richardt<sup>18</sup>, E. Kasai<sup>7</sup>, K. Katarzyński<sup>29</sup>, D. Khangulyan<sup>30</sup>, R. Khatoon<sup>8</sup>, B. Khélifi<sup>11</sup>, W. Kluźniak<sup>31</sup>, Nu. Komin<sup>22</sup>, K. Kosack<sup>16</sup>, D. Kostunin<sup>9</sup>, A. Kundu<sup>8</sup>, R.G. Lang<sup>18</sup>, S. Le Stum<sup>25</sup>, F. Leitl<sup>18</sup>, A. Lemière<sup>11</sup>, M. Lemoine-Goumard<sup>32</sup>, J.-P. Lenain<sup>15</sup>, F. Leuschner<sup>26</sup>, J. Mackey<sup>1</sup>, D. Malyshev<sup>26,*</sup>, G. Martí-Devesa<sup>27</sup>, R. Marx<sup>4</sup>, A. Mehta<sup>9</sup>, P. J. Meintjes<sup>33</sup>, A. Mitchell<sup>18</sup>, R. Moderski<sup>31</sup>, L. Mohrmann<sup>2</sup>, A. Montanari<sup>4</sup>, E. Moulin<sup>16</sup>, T. Murach<sup>9</sup>, M. de Naurois<sup>6</sup>, J. Niemiec<sup>21</sup>, S. Ohm<sup>9</sup>, E. de Ona Wilhelmi<sup>9</sup>, M. Ostrowski<sup>28</sup>, S. Panny<sup>27</sup>, M. Panter<sup>2</sup>, R. D. Parsons<sup>13</sup>, U. Pensec<sup>15</sup>, G. Peron<sup>11</sup>, D. A. Prokhorov<sup>34</sup>, G. Pühlhofer<sup>26,*</sup>, M. Punch<sup>11</sup>, A. Quirrenbach<sup>4</sup>, M. Regeard<sup>11</sup>, A. Reimer<sup>27</sup>, O. Reimer<sup>27</sup>, I. Reis<sup>16</sup>, H. Ren<sup>2</sup>, F. Rieger<sup>2</sup>, B. Rudak<sup>31</sup>, E. Ruiz-Velasco<sup>2</sup>, V. Sahakian<sup>3</sup>, H. Salzmann<sup>26</sup>, A. Santangelo<sup>26</sup>, M. Sasaki<sup>18</sup>, J. Schäfer<sup>18</sup>, F. Schüssler<sup>16</sup>, H. M. Schutte<sup>8</sup>, J. N. S. Shapopi<sup>7</sup>, S. Spencer<sup>18</sup>, Ł. Stawarz<sup>28</sup>, R. Steenkamp<sup>7</sup>, S. Steinmassl<sup>2</sup>, C. Steppa<sup>10</sup>, K. Streil<sup>18</sup>, I. Sushch<sup>8,*</sup>, T. Takahashi<sup>35</sup>, T. Tanaka<sup>36</sup>, A. M. Taylor<sup>9</sup>, R. Terrier<sup>11</sup>, C. Thorpe-Morgan<sup>26,*</sup>, M. Tluczykont<sup>37</sup>, T. Unbehaun<sup>18</sup>, C. van Eldik<sup>18</sup>, B. van Soelen<sup>33</sup>, M. Vecchi<sup>5</sup>, C. Venter<sup>8</sup>, J. Vink<sup>34</sup>, T. Wach<sup>18</sup>, S. J. Wagner<sup>4</sup>, F. Werner<sup>2</sup>, A. Wierzcholska<sup>21</sup>, M. Zacharias<sup>4,8</sup>, A. A. Zdziarski<sup>31</sup>, A. Zech<sup>14</sup>, and N. Żywucka<sup>8</sup>

*(Affiliations can be found after the references)*

Received 14 February 2024 / Accepted 29 April 2024

## ABSTRACT

PSR B1259-63/LS 2883 is a gamma-ray binary system that hosts a pulsar in an eccentric orbit, with a 3.4 yr period, around an O9.5Ve star (LS 2883). At orbital phases close to periastron passages, the system radiates bright and variable non-thermal emission, for which the temporal and spectral properties of this emission are, for now, poorly understood. In this regard, very high-energy (VHE) emission is especially useful to study and constrain radiation processes and particle acceleration in the system. We report on an extensive VHE observation campaign conducted with the High Energy Stereoscopic System, comprised of approximately 100 h of data taken over five months, from $t_p$ - 24 days to $t_p$ + 127 days around the system's 2021 periastron passage (where $t_p$ is the time of periastron). We also present the timing and spectral analyses of the source. The VHE light curve in 2021 is consistent overall with the stacked light curve of all previous observations. Within the light curve, we report a VHE maximum at times coincident with the third X-ray peak first detected in the 2021 X-ray light curve. In the light curve – although sparsely sampled in this time period – we see no VHE enhancement during the second disc crossing. In addition, we see no correspondence to the 2021 GeV flare in the VHE light curve. The VHE spectrum obtained from the analysis of the 2021 dataset is best described by a power law of spectral index $\Gamma = 2.65 \pm 0.04_{\text{stat}} \pm 0.04_{\text{sys}}$, a value consistent with the spectral index obtained from the analysis of data collected with H.E.S.S. during the previous observations of the source. We report spectral variability with a difference of $\Delta\Gamma = 0.56 \pm 0.18_{\text{stat}} \pm 0.10_{\text{sys}}$ at 95% confidence intervals, between sub-periods of the 2021 dataset. We also detail our investigation into X-ray/TeV and GeV/TeV flux correlations in the 2021 periastron passage. We find a linear correlation between contemporaneous flux values of X-ray and TeV datasets, detected mainly after $t_p$ + 25 days, suggesting a change in the available energy for non-thermal radiation processes. We detect no significant correlation between GeV and TeV flux points, within the uncertainties of the measurements, from $\sim t_p$ - 23 days to $\sim t_p$ + 126 days. This suggests that the GeV and TeV emission originate from different electron populations.

**Key words.** radiation mechanisms: non-thermal – binaries: general – pulsars: individual: PSR B1259-63/LS 2883

## 1. Introduction

Gamma-ray loud binaries (GRLBs) are a subclass of high-mass and intermediate-mass binary systems characterised by

their energy spectra peaking above 1 MeV, but typically at $E \gtrsim 100$ MeV, and extending to beyond 1 TeV. While hundreds of high-mass binaries have been detected in the X-ray band, the current generation of Cherenkov telescopes and gamma-ray satellites have only been able to detect about a dozen GRLB systems (see, e.g., Dubus 2013; Chernyakova et al. 2019, for recent reviews). The physical environments and mechanisms leading to

<footnote>* Corresponding authors;
e-mail: contact.hess@hess-experiment.eu</footnote>
the production of such energetic radiation in these systems are
not firmly established.

GRLB systems are comprised of a massive early-type star (spectral class O or B) and a compact object (a neutron star or a black hole). The nature of this compact object is difficult to discern in the majority of cases, in several systems however, the compact object has been identified as a non-accreting pulsar such as in the cases of PSR B1259-63/LS 2883, PSR J2032+4127 and LS I+61 303 (Johnston et al. 1992; Abdo et al. 2009a; Weng et al. 2022). Additionally, evidence of hard X-ray pulsations have been reported in the system LS 5039 (Yoneda et al. 2020), tentatively suggesting a neutron star companion as well.

The PSR B1259-63/LS 2883 system was discovered during a high-frequency radio survey intending to search for nearby pulsars (Johnston et al. 1992). Subsequent radio and optical observations resulted in the identification of the compact object in the system as a young radio pulsar (spin period ~48 ms), in a highly eccentric (e = 0.87) 3.4-yr (1236.724526 ± 6 × 10⁻⁶ day) orbit around the O9.5Ve star LS 2883 (Johnston et al. 1992, 1994; Negueruela et al. 2011; Shannon et al. 2014; Miller-Jones et al. 2018)<fnref n="1" />. The system is located at a distance of 2.39 ± 0.19 kpc from Earth (Gaia Collaboration 2018), and recent measurements of the inclination angle suggest that the binary orbit is observed at an angle of 154° to the line of sight (Miller-Jones et al. 2018). The projected semi-major axis is *a sin i* = 1296.27448 ± 0.00014 lt-s (Miller-Jones et al. 2018), which for the pulsar's orbital eccentricity corresponds to apastron and periastron separations of 11 AU and 0.8 AU, respectively. Additionally, the spin-down luminosity of the pulsar was estimated to be *L*<sub>sd</sub> = 8.2 × 10<sup>35</sup> erg s<sup>-1</sup> (Johnston et al. 1994), with a characteristic age of 330 kyr (Johnston et al. 2005).

The companion star LS 2883 has a bolometric luminos-
ity of $L_* = 2.3 \times 10^{38}$ erg s$^{-1}$ (Negueruela et al. 2011)
and hosts a decretion disc that extends up to at least
20 stellar radii (Johnston et al. 1992; Negueruela et al. 2011;
Chernyakova et al. 2014) from the star (0.56 AU). The radius of
LS 2883 is about $10R_\odot$ (0.05 AU) (Negueruela et al. 2011), and
its mass is ~24 $M_\odot$ (Shannon et al. 2014). The disappearance of
pulsed radio emission at ~$t_p$ – 16 days (where $t_p$ is the time of
periastron), and its reappearance at ~$t_p$ + 16 days (Johnston et al.
2005), as well as observations of the dispersion measure along
the periastron passage, both suggest that the stellar decretion
disc is inclined with respect to the orbital plane (Johnston et al.
1996). Measurements of this inclination angle between the plane
of the pulsar’s orbit and the circumstellar disc suggest an angle
of ~35° (Johnston et al. 1994; Shannon et al. 2014).

Following its optical and radio detection, the system was
later detected in the X-ray band with the ROSAT satel-
lite (Cominsky et al. 1994). In the X-ray regime, PSR B1259-
63/LS 2883 is detected during its entire orbit with a non-thermal,
non-pulsed spectrum (Marino et al. 2023). While the X-ray flux
level is minimal around apastron, close to the periastron passage
the keV light curve is typically characterised by two maxima
roughly coinciding with the times of the disappearance and re-
appearance of pulsed radio emission (see e.g. Chernyakova et al.
2015). These peaks are usually interpreted as being connected to
the pulsar crossing the Oe stellar disc. During the 2021 perias-
tron passage, the X-ray light curve exhibited a third maximum
between ~t<sub>p</sub> + 30 and t<sub>p</sub> + 50 days (Chernyakova et al. 2021)
(henceforth referred to as the third X-ray peak), in addition to

the two X-ray peaks at ~$t_p$ ± 16 days detected in all observed
periastron passages.

PSR B1259-63/LS 2883 was detected in the GeV band with Fermi-LAT (Abdo et al. 2011; Tam et al. 2011). At these energies, the system is characterised by a relatively low flux level in the period between $t_p$ – 30 and $t_p$ + 30 days. It later enters a high flux state (coined as the “GeV flare”) that has been detected following all periastron passages observed with Fermi-LAT to date (2010–2021 Abdo et al. 2011; van Soelen & Meintjes 2015; Caliandro et al. 2015; Wood et al. 2018; Chang et al. 2021). However, for all periastron passages to date during which very high-energy (VHE; ≥100 GeV) observations were taken contemporaneously with the corresponding GeV flare, no clear counterpart at very high energies has been seen (H.E.S.S. Collaboration 2013, 2020). The GeV flare in 2017 began after a noticeable delay, starting at up to ~50 days after the periastron passage (Chang et al. 2021). The light curve of the GeV flare obtained from the 2017 periastron passage also showed a number of extremely strong and rapid sub-flares on timescales as short as ~10 min. The observed luminosity of these sub-flares reached values of 30 times the spin-down luminosity of the pulsar (Johnson et al. 2018).

In the VHE band, PSR B1259-63/LS 2883 was detected with the High Energy Stereoscopic System (H.E.S.S.) for the first time in 2004 (Aharonian et al. 2005a), after which the array regularly observed the system at orbital phases close to its periastron passages. VHE observations of PSR B1259-63/LS 2883 are summarised in H.E.S.S. Collaboration (2020) which reports on five (2004–2017) periastron passages observed with H.E.S.S. (see also Aharonian et al. 2005a, 2009; H.E.S.S. Collaboration 2013, 2020, for individual analyses of the 2004–2017 periastron passages, respectively). See Table 1 for specific periastron passage dates and a summary of each passage’s VHE observation campaign.

The VHE light curve obtained from the stacked analysis
of the orbital-period folded data collected during the previous
observations of the system indicate the presence of an asymmet-
ric double peak profile (H.E.S.S. Collaboration 2020). Maxima
derived from a Bayesian block analysis of stacked data from pre-
vious periastron passages were reported between $t_p$ – 32 and
$t_p$ – 26 days (with a hint of a sub-peak at around $t_p$ – 15 days)
and between $t_p$ + 16 and $t_p$ + 57 days, with significances of 12.1$\sigma$
and 39.8$\sigma$, respectively.

In this work we present the results of the most recent H.E.S.S. observational campaign on PSR B1259-63/LS 2883, performed around the 2021 periastron passage. Extensive coverage of the system during this observation campaign has allowed an unprecedented amount of observational data to be taken post-periastron passage. In particular, observations extended up to the largest post-periastron orbital phase interval in the TeV band to date (29 days more than the previous longest in 2004), see Table 1 for details.

Following this introduction, Sect. 2 outlines the method-
ology and details of the H.E.S.S. array and its data pipeline.
Moreover, this section covers specific details of prior H.E.S.S.
observation campaigns and data analysis of the source during
periastron passages up to and including 2021. In Sect. 3 the
results of the analysis are presented, including studies of the flux
behaviour and light curve trends, as well as spectral analysis of
the source with a search for spectral variability. In this section
we also present our investigation into a correlation between the
X-ray/TeV flux and the GeV/TeV flux in the 2021 periastron pas-
sage. In Sect. 4 these results are discussed in the context of previ-
ous periastron passages and in the context of the unique findings

<footnote n="1">¹ In the following, we assume that the 2021 periastron of PSR B1259-63/LS 2883 occurred at $t_p = 59254.867359$ MJD.</footnote>
at other wavelengths in 2021. We also present some theoretical interpretation of the findings of this study. Finally, Sect. 5 con- tains our concluding remarks.

## 2. Method

H.E.S.S. is an array of five imaging atmospheric Cherenkov telescopes (IACTs), where each telescope is abbreviated and numbered CT1-5, and is located in the Khomas Highlands of Namibia (see Aharonian et al. 2006; H.E.S.S. Collaboration 2020, for detailed descriptions of the H.E.S.S. array).

In order to detect Cherenkov light, H.E.S.S. can only be oper- ated under dark conditions. Because of this, H.E.S.S. is not oper- ated during periods of bright moonlight (defined as above ~40% illumination). This results in a cycle-wise data taking period of 28 days. The fundamental data-taking unit of the H.E.S.S. array is an observational run, defined as a period of data acquisition lasting ~28 mins.

The VHE data presented in this paper are exclusively taken from runs where a minimum of three telescopes from CT1-4 were present (stereo mode). We use CT1-4 data to allow unbi- ased direct comparison to the majority of the other periastron passages covered by H.E.S.S., in which only CT1-4 data was available. For this reason CT5 data were not used in the anal- ysis. The analysis presented in this paper used the reflected regions background method, for light curve production and spec- tral analysis, as well as the ring background method for the creation of maps (see Berge et al. 2007, for further details of these analysis methods). Observations were performed using pointing offsets from the source position, all offsets were exclu- sively performed along right ascension due to the presence of the nearby bright source HESS J1303-631 at an angular separation of 0.75° North from PSR B1259-63/LS 2883. The dataset con- tained almost exclusively 0.5° telescope offsets with two runs of 221 (the total run number after data quality selection cuts had been applied) at an offset of 0.7°.

Prior to 2021, H.E.S.S. observed PSR B1259-63/LS 2883 covering a number of orbital phases close to previous peria- stron passages. These include the 2004 (Aharonian et al. 2005a), 2007 (Aharonian et al. 2009), 2011 (H.E.S.S. Collaboration 2013), 2014, and most recently the 2017 periastron passages (see H.E.S.S. Collaboration 2020, for both the 2014 and 2017 periastron passages). See Table 1 for details of observations in previous periastron passages.

During the 2021 periastron passage H.E.S.S. attained a total of 100 h of observations in the stereo configuration after data quality selection cuts had been applied. See Table 1 for details of the 2021 observations.

The results presented in this paper were produced using the HAP (H.E.S.S. Analysis Package)/ImPACT (Image Pixel-wise fit for Atmospheric Cherenkov Telescopes) template-based method chain (Parsons & Hinton 2014). Results have been cross-checked using the Paris Analysis chain (de Naurois & Rolland 2009).

All light curves and spectra in this work were produced from data that had passed the spectral quality selection cuts, represent- ing the strictest cut criterion for H.E.S.S. data (Aharonian et al. 2006). The data were also subject to a maximum event offset of 2.5°.

To estimate the systematic uncertainties we adopt the values outlined in Aharonian et al. (2006) for stereo analyses, as well as compare the reconstructed fluxes and spectral indices between the two major H.E.S.S. analysis chains. This study indicated a systematic uncertainty in the flux at an estimated level of 20%

and an uncertainty in spectral indices of 0.1. Statistical uncer- tainties on values/figures in this work (with the exception of spectral parameters that are reported at 95% confidence interval –c.i.–) are given at 68% c.i., unless explicitly stated otherwise.

In calculating the spectra in this study we utilise the forward- folding method (see Piron et al. 2001, for further information on the forward folding method).

## 3. Results

The PSR B1259-63/LS 2883 system is located at the J2000 coordinates RA = 13h02m47.65s, Dec = -63°50'8.6" and is situated in the Galactic plane (Johnston et al. 1992). It is near to by the bright source HESS J1303-631, a pulsar wind nebula that is spatially coincident with the pulsar PSR J1301- 6305 (H.E.S.S. Collaboration 2012). The significance map of the source and its surrounding region are shown in Fig. 1 using Li and Ma significances (Li & Ma 1983) and were created by utilising all H.E.S.S. data passing spectral cuts from the 2021 periastron passage ($t_p$ – 23 days to $t_p$ + 127 days). HESS J1303- 631 is known to have an energy-dependent morphology with a large spillover at GeV energies (H.E.S.S. Collaboration 2012; Acero et al. 2013). This spillover corresponds to an extended and energy dependent emission profile of the source, to a degree that it has the potential to contaminate the emission of nearby sources such as PSR B1259-63/LS 2883. This required us to ensure that the effect of spillover was non-existent or negligible at very- high energies by measuring the effect of the spillover in runs far from the periastron passage (using combined data in the period of ~$t_p$ + 100 days to $t_p$ + 500 days) where VHE emission from PSR B1259-63/LS 2883 was consistent with zero. No evidence of contaminant emission at VHE energies was found.

During the analysis of PSR B1259-63/LS 2883 a standard angular distance cut for point sources of 0.005 deg² was applied (defined as the angular distance between a reconstructed event and the expected source position).

The background acceptance ratio between the ON and OFF region had a value of $\alpha$ = 0.07, resulting in a total excess of 1668.40 events. In total, for an acceptance corrected live time of 100.02 h, we obtain a Li and Ma significance of 36.0 from PSR B1259-63/LS 2883.

### 3.1. Spectral analysis

A full investigation into the VHE spectral properties of the sys- tem during the periastron passage was undertaken and several spectra were derived. The total spectrum of the available peri- astron passage data was calculated, and the spectra of key inter- vals were created to investigate spectral variability. The first of these time frames included the two H.E.S.S. observational cycles (from $t_p$ – 3.9 days to $t_p$ + 15.3 days) that occurred concurrently with the periastron passage. Secondly, we created a spectrum for the period in which the peak levels of VHE flux were mea- sured (here defined as $t_p$ + 25 days to $t_p$ + 36 days). Addition- ally, we created a spectrum from the data contemporaneous with the 2021 GeV flare (here referring to the period $t_p$ + 55 days to $t_p$ + 108 days as defined in Chernyakova et al. 2021). Finally, we created a spectrum of the data from the final two H.E.S.S. obser- vational cycles from ~$t_p$ + 81 days to ~$t_p$ + 127 days (from now on referred to as the “TeV low flux” period). These datasets will henceforth be referred to as A, B, C, and D, respectively (please refer to Table 2). Each spectrum of these periods was fit with a power-law model, $dN/dE = \phi_0(E/E_0)^{-\Gamma}$, where $\Gamma$ represents the photon index of the power law with a normalisation $\phi_0$ and
<table><caption>**Table 1.** Summary of analysed H.E.S.S. observations of PSR B1259-63/LS 2883 from 2004 to 2021.</caption><thead><tr><th></th><th>2004</th><th>2007</th><th>2011</th><th>2014</th><th>2017</th><th>2021</th></tr></thead><tbody><tr><td>Start date</td><td>Feb. 27</td><td>Apr. 09</td><td>Jan. 10</td><td>Mar. 07</td><td>Aug. 10</td><td>Jan. 16</td></tr><tr><td>End date</td><td>Jun. 15</td><td>Aug. 08</td><td>Jan. 16</td><td>Jul. 21</td><td>Aug. 20</td><td>Jun. 16</td></tr><tr><td>Time - t<sub>p</sub> [days]</td><td>-7 to +98</td><td>-110 to +11</td><td>+26 to +32</td><td>-39 to +78</td><td>-42 to -37</td><td>-23 to +127</td></tr><tr><td>CT5 mono</td><td>N<sub>runs</sub></td><td>-</td><td>-</td><td>-</td><td>141</td><td>12</td><td>(216)</td></tr><tr><td></td><td>t<sub>L</sub> [h]</td><td>-</td><td>-</td><td>-</td><td>62.2</td><td>6</td><td>(99.4)</td></tr><tr><td></td><td>θ [°]</td><td>-</td><td>-</td><td>-</td><td>41.8</td><td>57</td><td>(42.8)</td></tr><tr><td>CT1-4 stereo</td><td>N<sub>runs</sub></td><td>138</td><td>213</td><td>11</td><td>163</td><td>-</td><td>221</td></tr><tr><td></td><td>t<sub>L</sub> [h]</td><td>57.1</td><td>93.9</td><td>4.8</td><td>68.1</td><td>-</td><td>100.0</td></tr><tr><td></td><td>θ [°]</td><td>42.5</td><td>45.1</td><td>47.6</td><td>41.9</td><td>-</td><td>42.8</td></tr></tbody></table>

**Notes.** The number of runs passing spectral quality selection cuts in a given periastron passage is defined as $N_{\text{runs}}$ (see text for further details), $t_L$ refers to a periastron passage's total acceptance-corrected observation time, and finally $\theta$ indicates the mean zenith angle of the periastron passage's observations. Values for years prior to 2021 are taken from H.E.S.S. Collaboration (2020) and references therein. CT1-4 stereo data are those runs in which third or more of the CT1-4 telescopes were active. CT5 mono data corresponds to the data obtained from only the central (and largest) telescope CT5. The 2021 CT5 mono observations are displayed in brackets as they are not presented or utilised in this study.

<figure><img src="image_1.png" /></figure>

**Fig. 1.** Significance and excluded significance maps. Left panel: VHE significance map displaying PSR B1259-63/LS 2883 (here labelled as HESS J1302-638) and the surrounding region. PSR B1259-63/LS 2883 itself is the central object in the image, where also of note is the nearby pulsar wind nebula HESS J1303-631 (Aharonian et al. 2005a) directly to the north (see text for further details on this source). Also shown in blue in the lower left of the image, is the 68% c.i of the point spread function for these observations. Right panel: the significance map of PSR B1259-63/LS 2883 and its surrounding region after masking the two sources. The bins in both maps are correlated within a circle of radius 0.1°.

a decorrelation energy of $E_0$. The best-fit parameters of these models are presented at 95% c.i. unless otherwise stated.

To define the energy range for the spectral analysis, two different approaches were used. For the total 2021 spectrum the lower energy bound, at 0.27 TeV, was defined using Monte Carlo simulations which ensure that the energy reconstruction bias is less than 10% of the energy (Aharonian et al. 2006). The upper bound of the energy range was defined by the highest energy bin that could be fit with a significance of $2\sigma$. Henceforth we refer to this energy range as unfixed.

However, for the total spectrum used for comparison to the sub-periods, and for the sub-periods themselves, a fixed energy range was applied allowing accurate comparison of the different spectra. Thus, we apply an energy fitting range of (0.4–10.0) TeV for these sub-periods. The lower bound was chosen such that it supersedes the safe energy threshold for any of the data subsets. The higher energy threshold was chosen to ensure sufficient statistics up to the cut energy for all subsets.

Figure 2 shows the total 2021 periastron passage spectrum where no pre-fixed energy range for the fit was applied. This spectrum includes all the data taken with H.E.S.S. during

the 2021 periastron passage and spans the energy range (0.30–39.6) TeV (the centres of the lowest and highest energy bins, respectively). Figure 2, shows that the data largely follow a power law, however, there are hints that the spectrum may contain substructure. These substructures could be a result of systematic effects, though an investigation into these effects is beyond the scope of this paper.

Figure 3 displays two of the three sub-spectra (datasets B and D) created to investigate the spectral behaviour of the system over the course of a single periastron passage. The inclusion of dataset D (the TeV low flux period) allows direct comparison between two unique flux states of the system to search for spectral variability. Additionally the total spectrum of the periastron passage, calculated with a fixed energy range, is shown in this figure. The data points in both Fig. 2 and for the sub-spectra in Fig. 3 were binned to ensure that every flux point has a statistical significance of at least $2\sigma$.

We find a spectral index for the total periastron passage of $\Gamma = 2.78 \pm 0.05_{\text{stat}} \pm 0.10_{\text{sys}}$ (for the fixed energy range spectrum) which is consistent with the spectral index of previous years, $\Gamma = 2.76 \pm 0.03_{\text{stat}} \pm 0.10_{\text{sys}}$ (H.E.S.S. Collaboration 2020).
<table><caption>**Table 2.** Comparison of the best-fitting spectral parameters when a power-law model is fit to different datasets within the 2021 periastron passage of PSR B1259-63/LS 2883.</caption><thead><tr><th rowspan="2">Dataset</th><th rowspan="2">Time - t<sub>p</sub><br>days</th><th rowspan="2">Spectral index (Γ)</th><th>Normalisation (φ<sub>0</sub>)</th><th>Normalisation (φ<sub>0</sub>)</th></tr><tr><th>decorrelation<br>energy (E<sub>0</sub>)<br>10<sup>-12</sup> TeV<sup>-1</sup> cm<sup>-2</sup> s<sup>-1</sup></th><th>(1 TeV)<br>10<sup>-12</sup> TeV<sup>-1</sup> cm<sup>-2</sup> s<sup>-1</sup></th></tr></thead><tbody><tr><td>Total 2021<br>(0.27-33.6 TeV)</td><td>-23.58 to +127.26</td><td>2.65 ± 0.04<sub>stat</sub> ± 0.10<sub>sys</sub></td><td>1.28 ± 0.05<sub>stat</sub> ± 0.26<sub>sys</sub><br>(0.95 TeV)</td><td>1.13 ± 0.04<sub>stat</sub> ± 0.23<sub>sys</sub></td></tr><tr><td>Total 2021<br>(0.4-10 TeV)</td><td>-23.58 to +127.26</td><td>2.78 ± 0.05<sub>stat</sub> ± 0.10<sub>sys</sub></td><td>1.31 ± 0.05<sub>stat</sub> ± 0.26<sub>sys</sub><br>(0.95 TeV)</td><td>1.15 ± 0.04<sub>stat</sub> ± 0.23<sub>sys</sub></td></tr><tr><td>2021 Periastron period (Dataset A)<br>(0.4-10 TeV)</td><td>-3.9 to +15.3</td><td>2.75 ± 0.11<sub>stat</sub> ± 0.10<sub>sys</sub></td><td>1.52 ± 0.12<sub>stat</sub> ± 0.30<sub>sys</sub><br>(0.95 TeV)</td><td>1.34 ± 0.10<sub>stat</sub> ± 0.27<sub>sys</sub></td></tr><tr><td>2021 Peak TeV flux (Dataset B)<br>(0.4-10 TeV)</td><td>+25 to +36</td><td>2.98 ± 0.07<sub>stat</sub> ± 0.10<sub>sys</sub></td><td>5.00 ± 0.22<sub>stat</sub> ± 1.00<sub>sys</sub><br>(0.79 TeV)</td><td>2.45 ± 0.12<sub>stat</sub> ± 0.49<sub>sys</sub></td></tr><tr><td>2021 GeV flare period (Dataset C)<br>(0.4-10 TeV)</td><td>+55 to +108</td><td>2.42 ± 0.10<sub>stat</sub> ± 0.10<sub>sys</sub></td><td>0.38 ± 0.03<sub>stat</sub> ± 0.08<sub>sys</sub><br>(1.27 TeV)</td><td>0.68 ± 0.06<sub>stat</sub> ± 0.14<sub>sys</sub></td></tr><tr><td>TeV Low flux period (Dataset D)<br>(0.4-10 TeV)</td><td>+81 to +127</td><td>2.42 ± 0.17<sub>stat</sub> ± 0.10<sub>sys</sub></td><td>0.17 ± 0.02<sub>stat</sub> ± 0.03<sub>sys</sub><br>(1.40 TeV)</td><td>0.37 ± 0.05<sub>stat</sub> ± 0.07<sub>sys</sub></td></tr><tr><td>Average of previous periastron passages</td><td>-106 to +98</td><td>2.76 ± 0.09<sub>stat</sub> ± 0.10<sub>sys</sub></td><td>-</td><td>1.55 ± 0.12<sub>stat</sub> ± 0.31<sub>sys</sub></td></tr></tbody></table>

**Notes.** Γ represents the best-fit photon index of the power law $dN/dE = \phi_0(E/E_0)^{-\Gamma}$, with $\phi_0$ denoting the best-fit normalisation level, and $E_0$ being the decorrelation energy (Abdo et al. 2009b). Presented here are the spectral properties of the total dataset of the periastron passage (for both fixed and unfixed energy ranges), as well as four additional sub-orbital periods from the 2021 dataset. All 2021 spectra in this table (with the exception of the total periastron passage spectrum marked “unfixed energy range”) are calculated using an energy range of (0.4–10.0) TeV to enable comparison. The average spectral properties of H.E.S.S. stereo analysis from previous periastron passages (2004, 2007, 2011 and 2014, taken from H.E.S.S. Collaboration 2020) are also included for comparison, from an averaging of the values reported in these papers. The spectra of previous periastron passages, however, do not use the same fixed energy values as the 2021 data. The two errors associated with each value are representative of, first, the magnitude of the statistical uncertainty in the value at a 95% confidence level, and secondly the systematic error in the measurements as adopted from Aharonian et al. (2006).

<figure><img src="image_2.png" /><figcaption>**Fig. 2.** Total spectrum of PSR B1259-63/LS 2883’s 2021 periastron passage. The spectrum was produced from all data taken with H.E.S.S. during the 2021 periastron passage observations of PSR B1259-63/LS 2883, using an unfixed energy range. The spectrum has been fit with a power-law model. For details on the properties of the spectra displayed see Table 2. The red band indicates the 68% c.i. of the statistical error for the fitted model.</figcaption></figure>

We note a statistically significant difference between the spectral indices obtained through the power-law model describing dataset B and D (see Table 2). Accounting for the uncertainties, the spectral index change between these two datasets is $\Delta\Gamma = 0.56 \pm 0.18_{\text{stat}} \pm 0.10_{\text{sys}}$, implying a sub-orbital spectral variation at a c.i. of greater than 95%.

<figure><img src="image_3.png" /><figcaption>**Fig. 3.** Comparison of spectra from PSR B1259-63/LS 2883’s 2021 periastron passage. Shown is the total spectrum of the 2021 periastron passage compared to two sub-spectra. Each model was calculated with the forward-folding method, using a power law in an energy range of (0.4–10.0) TeV to allow a comparison between them. Comparing spectra allowed us to search for VHE spectral variability on the scale of a single periastron passage. Displayed in red (black circles) is the spectrum of the total 2021 periastron passage, in cyan (cyan diamonds) the spectrum from the dataset D and in magenta (magenta triangles) the spectrum of dataset B. For details on the properties of the spectra displayed see Table 2. Shaded regions indicate the 68% c.i. for the fitted model.</figcaption></figure>

For the total unfixed spectrum, we attempted to fit an exponentially cut-off power-law model in order to determine if the data shows a preference for a high-energy cut-off. This revealed
<figure><img src="image_4.png" /><figcaption>**Fig. 4.** VHE light curve of PSR B1259-63/LS 2883's 2021 periastron passage. VHE Data is shown in different binnings: night-wise binned fluxes (blue triangles), cycle-wise binned fluxes (red dots), and stacked data from H.E.S.S. observations of previous periastron passages (green diamonds). Horizontal error bars correspond to times of the earliest and latest runs that were merged to make the data point. Cycle-wise data contain the merged runs within one observational cycle (~28 days). Stacked data points have a weekly binning and correspond to data from the 2004, 2007, 2011, 2014 and 2017 periastron passages; these points were taken from H.E.S.S. Collaboration (2020). The two dashed black lines at $t_p$ – 16 days and $t_p$ + 16 days correspond to the time at which the disc crossing is thought to occur, the dot-dashed red line at $t_p$ = 0 days represents the point of periastron in the system. Finally, the shaded areas with arrows are displayed as a visual representation of the time periods from which sub datasets were taken. These regions are also marked with the letters of the dataset they represent. The grey region represents the time period encompassed by dataset A, the magenta region corresponds to the time period of dataset B, and the yellow region and the cyan regions indicate the time periods of datasets C and D, respectively; see Table 2 for the full details of these sub-periods.</figcaption></figure>

that a model with a cutoff is not preferred, with a lower limit on the cut-off energy of $E_C^{95\%} = 27.1$ TeV.

## 3.2. Flux analysis

For the 2021 periastron passage dataset, light curves were produced for the H.E.S.S. data in two different integration timescales, see Fig. 4. These were: night-wise binning and cycle-wise binning – grouping runs to individual observational cycles of ~28 days.

Individual flux points and their uncertainties were calculated using a reference spectrum, in this work this was a power-law model in the energy range (0.4–100.0) TeV. An index of $\Gamma = 2.65$ was used, corresponding to the total 2021 dataset spectral index value.

We performed a search for variability in PSR B1259-63/LS 2883 over a number of different timescales, however, statistical uncertainties prevented us from establishing the presence of variability on run-to-run timescales<fnref n="2" />. Our analysis of the VHE data at 25–35 days after periastron (the time period of fastest VHE flux increase) indicates that a model with a linearly increasing

flux is a better fit to the data in this period than a constant flux model. This was determined by comparing the chi-squared values of a linear increase model, and a constant flux model, in this time period. The comparison of the fits of these methods showed that a linear flux increase is preferred at a ~4σ level. During this time period the flux increased by a factor of two. Other than this increase in a period of ~10 days, we did not find significant evidence for a linear flux increase at shorter timescales. Thus, we see variability on timescales of down to ~10 days. It is possible that there exists variability on shorter timescales, however, we are unable to probe this due to statistics.

We investigated the impact of using an assumed spectral index of $\Gamma = 2.65$ to calculate the night-wise fluxes, given the discovery of sub-orbital spectral index variation. We investigated this by calculating binned fluxes using the two extreme values of the spectral index $\Gamma = 2.42$ (from the spectrum of the emission from dataset D) and $\Gamma = 2.98$ (dataset B). We then evaluated the difference between the nightwise fluxes of the two light curves that these indices produced. The percentage difference between the flux of the two new light curves yielded a maximum systematic error in the flux of ±10%, a comparable value to that of Aharonian et al. (2006) from which the systematic error values of this study were taken (see Sect. 2 for details). This value represents an additional systematic flux error in the light curves exclusively, and does not have an impact on any scientific conclusions in the paper.

<footnote n="2">² Most commonly, subsequent runs were taken during the same or the following night, with the exception of several breaks due to moonlight or bad weather. Thus, run-to-run timescales range from a few hours to a few days.</footnote>
Although the 2021 VHE light curve presented in Fig. 4 shows an overall trend similar to the light curve obtained from the stacked analysis of the orbital-period-folded data collected during previous observations (H.E.S.S. Collaboration 2020), we argue that a detailed comparison of the system's flux behaviour is complicated by the different coverage of the H.E.S.S. datasets. Despite observing at orbital phases close to the second disc crossing in the 2021 dataset, we do not see a VHE flux enhancement around this time. However, we report a VHE maximum occurring between $t_p + 20$ and $t_p + 50$ days (seen during the period of the maximum reported in previous years, H.E.S.S. Collaboration 2020).

### 3.3. 2021 GeV flare

The 2021 GeV flare (shown in Fig. 5) differed in considerable ways from those of previous periastron passages (although the GeV behaviour appears inherently variable between periastron passages). As in 2017, the 2021 GeV flare started at ~$t_p + 55$ days with GeV activity extending to ~$t_p + 108$ days, see Chernyakova et al. (2021). The system underwent numerous rapid and energetic sub flares on very short timescales (in some cases as short as ~10 min) reaching up to 30 times the spin-down luminosity of the pulsar.

We see no correspondence to the 2021 GeV flare, from ~$t_p + 55$ days to ~$t_p + 108$ days, in the VHE light curve (see Sect. 3.5 for further investigation into this). We do, however, note that our ability to monitor this is somewhat complicated by a large gap in our observations during the 2021 GeV flare period, as no observations were performed in the time frame of ~$t_p + 65$ days to ~$t_p + 81$ days.

The spectrum of dataset C (derived during the 2021 GeV flare period) has a spectral index notably similar to that of dataset D. There is, therefore, also a discrepancy in the spectral index between datasets B and C at $\Delta\Gamma = 0.56 \pm 0.18$, the same level as the previously discussed discrepancy between dataset B and D.

### 3.4. X-ray-TeV correlation

We investigated a potential correlation between VHE and X-ray flux in the 2021 periastron passage data. In this study we utilise the results reported in Chernyakova et al. (2021), from the *Neil Gehrels Swift (Swift)* X-ray telescope (XRT) and the Neutron Star Interior Composition Explorer NICER (both in the 0.3–10 keV range), covering the time period January, 19th, 2021 to May, 24th, 2021 ($t_p - 21$ to $t_p + 103$ days).

In order to perform the correlation study on timescales relevant to the system's behaviour, we binned the data from H.E.S.S. on a nightly basis, resulting in a total of 57 TeV points, compared to a total of 96 X-ray points binned by observations (32 of which from NICER, 64 from Swift). X-ray points had an average separation in time of 1.24 days (excluding time gaps of greater than a week) over the whole dataset.

A sub-selection of observations from both the TeV and X-ray datasets was made, ensuring the selection of only X-ray points occurring within one day of a TeV point. This correlation timescale of a day was selected because this is the shortest timescale in which the available statistics could confirm a lack of variability of PSR B1259-63/LS 2883 in X-rays. To make this sub-selection, we iterated over the data in steps of correlation timescale, where any TeV and X-ray points within this time became a correlated pair. Instances where two or more points of

the same data type (X-ray or TeV<fnref n="3" />) were found to be within a single correlation timescale were handled by averaging the flux, time and uncertainties of the respective points. Any data points that did not contain a counterpart within one day of the point were not considered in the correlation study.

After this selection a total of 26 correlated pairs were found. This selection of correlated pairs and their distribution across the periastron passage can be seen in Fig. 5. The first correlated pair is at a time of $t_p - 1.53$ days, extending up to the time of the final pair at $t_p + 97.47$ days. The majority of pairs occurred at times later than $t_p + 50$ days.

Figure 6 shows the results of the correlation investigation between the two datasets. By minimising $\eta^2$ (where $\eta^2$ is a linear combination of $\chi^2$ tests, see appendix A for details of the $\eta^2$ test) we obtained the following best-fitting values of the linear fit parameters for the model $F_X = aF_{TeV} + b$: $a = 2.62^{+0.41}_{-0.38}$ and $b = 0.50^{+0.12}_{-0.13} \times 10^{-11}$ erg cm$^{-2}$ s$^{-1}$ (where $F_X$ and $F_{TeV}$ are the X-ray and TeV fluxes, respectively). These values of the fit gave a total $\eta^2 = 105.08$ for 26 correlated pairs, resulting in $\bar{\eta}^2 = 2.10$.

To estimate the statistical uncertainties we performed numerical simulations, considering $N = 10^6$ random trial datasets. The integrated X-ray and TeV photon flux for each trial dataset were simulated from the original data, assuming a Gaussian distribution of uncertainties. The quoted errors for each parameter correspond to a 68% c.i. of all best-fit values obtained during random trial datasets when fit with the same model. We estimated the chance probability of finding a correlation at $\bar{\eta}^2 = 2.10$ by comparing the number of trials that provided better $\eta^2$ values than the original data, to the number of trials with a worse $\eta^2$. Making this comparison we found a chance probability of $2.27 \times 10^{-3}$.

We therefore conclude that there is a positive correlation between the X-ray and TeV flux during the time periods of the 2021 periastron that were probed by the study. While all points follow this linear trend, there are two notable outliers from the correlation. These pairs represent X-ray and TeV flux points that were measured shortly before $t_p + 16$ days (see Figs. 5 and 6). With regards to the conclusion of a linear correlation, it is important to note that this initial study, which placed no restrictions on the time of the flux points used in the correlation, consisted mostly of flux points that were measured at times greater than $t_p + 25$ days (see Fig. 5). This uneven sampling across the periastron passage prevents us from establishing the presence of such a correlation before this time period.

We separately assessed the correspondence of the VHE flux level to either the third X-ray peak or to the gradual decay seen in the X-ray flux profiles of previous years. To achieve this, we performed model fitting on the full dataset of 2021 VHE flux data, separately fitting both a negative exponential function and a negative exponential function combined with a Gaussian. Here, the negative exponential model is representative of the X-ray flux behaviour (corresponding well to the behaviour seen in previous years) and the Gaussian represents the flux profile of the third X-ray peak in 2021. In this comparison, the VHE data was better fit by the negative exponential function summed with a Gaussian, at a $5.5\sigma$ level. However, this fit is based on the available VHE data points that occur only immediately before the peak (see Fig. 5). Therefore, the limited number of points immediately after the peak makes it difficult to draw more robust claims.

<footnote n="3">³ Formally, due to the TeV data being binned on a nightly basis it was not possible for two runs to occur within the same correlation timescale of a day. Thus this step only affected X-ray data points, which frequently occurred within a day of other X-ray points</footnote>
<figure><img src="image_5.png" /><figcaption>**Fig. 5.** Comparison of the 2021 GeV light curve to the correlated X-ray and VHE light curves. Upper panel: the GeV light curve of PSR B1259-63/LS 2883's 2021 periastron passage observed with Fermi-LAT (0.1–10 GeV). Figure adapted from Chernyakova et al. (2021). Data points represent daily binnings, with green points indicating 95% upper limits and red points denoting detections (TS > 4). The horizontal black dashed line marks the flux value corresponding to the spin-down luminosity of PSR B1259-63/LS 2883's pulsar $L_{sd} = 8.2 \times 10^{35}$ erg s⁻¹ (see Chernyakova et al. 2021). Lower panel: light curve of PSR B1259-63/LS 2883's 2021 periastron passage displaying both VHE flux points from H.E.S.S. and X-ray data from the Swift and NICER observatories. Points with no transparency are points that were selected by the time correlation step and fell within a day of an alternate type of observation (see text for further details). Translucent points are data from each respective dataset that did not pass the time correlation step. The lines at ~$t_p \pm 16$ days represent the assumed disc crossing times and the line at $t_p = 0$ days marks the time of periastron. In both panels, coloured regions correspond to the orbital phases of correlation points in Fig. 6 that share the same colour. Here the grey region represents orbital phases before $t_p + 16$ days, the purple region represents the time frame $t_p + 16$ to $t_p + 30$ days, the cyan region is the time frame $t_p + 30$ to $t_p + 50$ days (coincident with the third X-ray peak reported in 2021 X-ray data), and the orange region corresponds to any time after $t_p + 50$ days.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>**Fig. 6.** Linear correlations between X-ray-TeV and GeV-TeV data sets. Left panel: linear correlation ($F_X = aF_{\text{TeV}} + b$) between the flux of the TeV and X-ray datasets, where exclusively integrated energy flux points measured within a day of one another were utilised. Shown in red is the line of best-fit resulting from a linear fit with an offset. A large outlier can be seen at TeV flux ~0.3 × 10⁻¹¹ erg cm⁻² s⁻¹; this outlier originates from rapid X-ray variability in the period from periastron up to ~$t_p$ + 16 days. The dashed lines mark the zero points of each axis. As can be seen, there is a non-zero X-ray flux value at zero TeV flux in the fit. Right panel: correlation plot of the flux of the TeV and GeV datasets, where exclusively flux points measured within a day of one another were used. Note that transparent points correspond to upper limits taken from the *Fermi*-LAT daily light curve, where each upper limit has been assumed to have a flux of half the upper limit value and a flux error corresponding to a 95% c.i. on the *Fermi* upper limit. GeV flux points have been adapted from Chernyakova et al. (2021). In both figures the orbital phases from which correlation points are taken are denoted by their colour (see Fig. 5).</figcaption></figure>

In addition, we undertook a separate correlation study utilising only data occurring after the time of the second X-ray peak in 2021 ($t_p$ + 16 days). The data from before this point, notably, were the cause of the blue coloured outliers at greater than 68% c.i. seen in the linear correlation in Fig. 6. We apply exactly the same method as described previously. This results in a goodness-of-fit value of $\eta^2 = 57.77$ that, for 22 correlated pairs, results in a $\bar{\eta}^2 = 1.38$ (compared to a value of $\bar{\eta}^2 = 2.10$ for the unrestricted dataset) and gave a chance probability of $8.3 \times 10^{-4}$.

### 3.5. GeV-TeV correlation

We also conducted a study into the correlation between the 2021 TeV and GeV datasets to further quantify the apparent absence of a TeV counterpart to the 2021 GeV flare. We utilise daily-binned GeV flux data from Chernyakova et al. (2021) for comparison with the nightly-binned TeV data. The method of correlation used was identical to that of the previous study between X-ray and TeV data. Despite known sub-flares at GeV energies, sometimes on timescales of ten minutes, we opted to utilise a timescale of one day. This correlation length matched the binning of the two light curves, and is also consistent with the X-ray-TeV correlation results.

Following the correlated pair selection, a total of 56 correlated pairs (from 57 TeV points and 187 GeV points) were selected. The first of these is at a time of $t_p - 23.5$ days, extending up to $t_p + 126.5$ days. Because all but one TeV points are time-correlated with a GeV point, in the GeV/TeV correlation these pairs are relatively evenly sampled across the time range of the periastron passage, reflecting the distribution of the TeV points. The majority of the correlated GeV points, however, are upper limits from the *Fermi* analysis (144 of 187 points). We therefore tested for a correlation using several approaches. Firstly, we simply omitted the upper limits from the correlation study. This resulted in an $\bar{\eta}^2 = 11.56$. However such a selection introduces a

bias towards high GeV fluxes and could mask an existing correlation. We therefore examined two approaches including upper limits, these methods corresponded to adopting a Gaussian distribution as the probability density function (PDF) for the upper limits (see Kelly 2007, for further details and alternative treatments of the PDF). The first of these approaches was to utilise half the upper limit value as the flux value and a dispersion corresponding to 95% c.i. of the *Fermi* upper limits. The second method utilised a zero value as the flux estimate (and once again a dispersion corresponding to 95% c.i. of the *Fermi* upper limits). The resulting $\bar{\eta}^2$ values were 15.02 and 10.67 for the Gaussian centered on half the upper limit value and the Gaussian centered at zero, respectively. All of these tests excluded a correlation at levels greater than $5\sigma$. We conclude that (within the uncertainties of the measurements) we detect no significant GeV-TeV correlation throughout the entire probed time period.

## 4. Discussion

Understanding the broadband emission from GRLBs is a complex problem, which still awaits a definite solution. Despite the difficulties, some progress, however, has been made in the modeling of emission from these systems. PSR B1259-63/LS 2883 contains a non-accreting pulsar, thus, in what follows we discuss the properties of the emission in the framework of a binary pulsar model. This scenario implies that the relativistic outflow from a rotation powered pulsar interacts with the stellar wind which, in the case of PSR B1259-63/LS 2883, consists of a radiation-driven polar wind, and a significantly more dense Keplerian-like decrease disc.

### 4.1. Orbital dependence of the X-ray and TeV emission

The termination of the pulsar wind occurs at a distance of $R_{\text{ts}}$ (from the pulsar) where the ram pressure of the stellar and
pulsar winds are equal. Given the much higher anticipated speed of the pulsar wind, the energy injection of non-thermal particles into the interaction region is dominated by contributions from the pulsar. Therefore, one expects that the radiation processes in binary pulsar systems are similar to those taking place in pulsar wind nebulae (Tavani et al. 1994; Tavani & Arons 1997), however, with the caveat of some important modifications.

Firstly, as the magnetic field is provided by the pulsar wind, a smaller termination distance necessarily implies a stronger magnetic field:

$$B \lesssim B_{\text{max}} = \sqrt{\frac{L_{\text{sd}}}{cR_{\text{ts}}^2}} \\ \approx 3 \left( \frac{L_{\text{sd}}}{8.2 \times 10^{35} \text{ erg s}^{-1}} \right)^{1/2} \left( \frac{R_{\text{ts}}}{0.1 \text{ AU}} \right)^{-1} \text{ G}, \quad (1)$$

where $L_{sd}$ is the pulsar's spin down luminosity. The second important difference is that the photon field is dominated by contributions from the optical companion. This provides an intense photon field with an energy density of

$$w_{\text{ph}} = \frac{L_*}{4\pi c R^2} \\ \approx 3 \left( \frac{L_*}{2.3 \times 10^{38} \text{ erg s}^{-1}} \right) \left( \frac{R}{1 \text{ AU}} \right)^{-2} \text{ erg cm}^{-3}, \quad (2)$$

where $R$ is the separation between the star and the pulsar (system separation). For simplicity, we assume that the production region is located close to the pulsar.

For a Gauss-strength magnetic field, VHE electrons generate synchrotron emission in the hard X-ray band. Binary pulsar systems were predicted, therefore, to be TeV sources, provided that the energy density of the stellar photon field is comparable to the expected energy density of the magnetic field (Kirk et al. 1999). For an accurate calculation of the expected TeV flux level, one needs to account for a number of effects including the Klein-Nishina cutoff, IC scattering in the anisotropic regime, and gamma-gamma attenuation (Kirk et al. 1999).

Under the assumption of isotropic winds, the pulsar wind termination distance is proportional to the system separation distance, $R_{ts} \propto R$. Thus, the ratio of the photon to magnetic field energy density does not depend on the orbital phase, and one may expect quite similar X-ray and TeV light curves unless $\gamma-\gamma$ attenuation is significant (Dubus 2006; Khangulyan et al. 2007; Sushch & van Soelen 2017, 2023). However, one needs to take into account that some physical parameters can change their values depending on the orbital phase. For example, the magnetic field strength, which is expected to be proportional to the distance to the termination shock, may undergo a significant change with orbital phase. Consequently, this may induce a change of the cooling regime and/or of the synchrotron component (Khangulyan et al. 2012; Dubus & Cerutti 2013).

Unless the stellar disc or locally generated fields provide significant targets for IC scattering, the temperature of the target photons does not vary with orbital phase. However, one needs to account for the change of the scattering angle. For an orbital inclination angle of $i \approx 154^\circ$ (Miller-Jones et al. 2018), and for a production region in the orbital plane during the epochs of the H.E.S.S. observations (see Fig. 7), the IC scattering angle is approximately $65^\circ$. In this case, an emission of energy 1 TeV may be generated by electrons with energy $E_{\text{TeV}} \approx 1.6$ TeV; here we adopt a photon field temperature of $T_* \approx 3 \times 10^4$ K (Negueruela et al. 2011), and use the approximation from Khangulyan et al. (2014).

<figure><img src="image_7.png" /><figcaption>Fig. 7. Comparison of the scattering angle of IC processes to separation of the two objects comprising PSR B1259-63/LS 2883. Left axis, green curve: The angle between the line-of-sight and the direction from the optical star at the pulsar's location, as a function of time to periastron passage. If the production region is close to the pulsar, this angle is equal to the scattering angle for IC processes. Right axis, blue curve: The ratio of the separation distance, $R$, to the periastron separation. The shaded regions in this plot represent the periods of the sub spectra and are defined as in Fig. 4 (see Table 2 for the full details of these sub-periods).</figcaption></figure>

Because of the eccentric orbit of the pulsar in the system, the system separation changes by a factor of four during the period relevant for the H.E.S.S. observations. This will therefore induce a proportional change of the magnetic field strength in the production region, meaning that the energy of the X-ray emitting electrons may change by a factor of $\approx 2$. For a typical X-ray spectrum slope of 1.5 (the average value obtained from Swift observations by Chernyakova et al. 2021), X-ray emitting electrons have an $E^{-2}$ energy distribution and so, even in an idealised case of isotropic winds, the relationship between X-ray and TeV luminosity should depend on separation as:

$$L_X \propto L_{\text{TeV}} R^{1/2}. \quad (3)$$

In Fig. 6, the linear fit is mostly constrained by pairs of correlated X-ray and TeV runs occurring at $t > t_p + 50$ days, i.e., when $R$ is large and changes more slowly with time. For smaller separation distances, using Eq. (3) one should expect that the linear fit overestimates the X-ray flux level. However, from Fig. 6 one can see that for certain pairs the measured X-ray flux is significantly higher than the value given by the linear fit. This could be considered as a hint of the wind interaction in a non-isotropic regime (e.g a Keplerian decretion disc, a non-isotropic pulsar wind, or changes in the scattering angle between relativistic electrons and soft photons). Indeed, the points with high relative X-ray flux correspond to orbital phases close to $t_p + 16$ days, where the pulsar may interact with the stellar disc. Providing a significantly dense stellar disc, the pulsar wind will terminate significantly closer to the pulsar, enhancing the magnetic field strength without a proportional enhancement of the photon field (Khangulyan et al. 2007; Takata et al. 2012; Sushch & Böttcher 2014). This results in increasing X-ray flux and a softer X-ray spectrum during the disc crossing, consistent with available observations (see e.g., Chernyakova et al. 2014).

Analysis of the X-ray–TeV emission correlation reveals that there is a contribution to the X-ray flux that depends weakly on the TeV flux level (the contribution due to the constant term in the linear relationship). This could indicate the presence of two
or more zones that generate non-thermal emission. This is also supported by the absence of a strong correlation between the X-ray and radio emission. Up to about 30 days after periastron, radio and X-ray emission show a very good correlation, but following this the X-ray flux starts to increase while the radio flux continues decreasing (see in Chernyakova et al. 2021).

This third X-ray peak, occuring 30–50 days after periastron (Chernyakova et al. 2021), has not been reported in previous periastron passages. Although we lack TeV observations in 2021 during a larger fraction of this period, there is good evidence (from a significant TeV flux rise around $t_p$ + 30 days, a good correlation with X-ray data in this period, and from the light curve template fitting discussed in Sect. 3.2) that there is a correspondence of the third X-ray peak in the 2021 TeV light curve. However, because of the gap in TeV coverage after $t_p$ + 35 days, the time at which the maximum occurs in the TeV light curve is not well constrained and could be shifted with respect to the X-ray peak.

Summarising the multiwavelength data from Chernyakova et al. (2021) during the period 30–60 days post-periastron, we see that the radio emission decreases, X-ray emission increases and then decreases, GeV emission stays in a low-emission state. From the H.E.S.S. data, we see that the TeV emission increases and then decreases. Immediately following this period, the GeV emission increases strongly with no corresponding increase in any other band. Given the variation seen in emission profiles, it is difficult to reconcile all of these observational trends with a simple one-zone model for the post-periastron time-evolution of the non-thermal emission. A multi-zone configuration can be produced by the complex geometry and dynamics of the interaction between the pulsar and stellar winds (Bogovalov et al. 2008; Bosch-Ramon et al. 2012; Dubus et al. 2015; Huber et al. 2021), and it appears likely that such models are required to explain the data. The correlation of the TeV and X-ray light curves 30 days after periastron, suggests either that the electron population responsible for the third X-ray peak also emits in the TeV regime. Alternatively, given that the observed TeV light curve is compatible with previous periastron passages, it is possible that the X-ray emission accompanying the TeV peak was suppressed for some reason during previous periastron passages.

The nature of the GeV flare, which is not accompanied by an increase in emission at any other waveband, remains puzzling. This scenario could potentially be connected to a complex evolution of the wind termination shock, i.e. strong confinement of the pulsar wind due to either the eccentricity of the orbit (Barkov & Bosch-Ramon 2016), or the interaction with the circumstellar disc (Khangulyan et al. 2012) followed by a rapid expansion of the pulsar wind bubble later on.

## 4.2. Spectral variability

Another important finding in the H.E.S.S. 2021 dataset is spectral variability of the VHE emission. While in other GRLBs spectral variability is an established feature of the TeV emission (most notably in LS 5039, Aharonian et al. 2005b), the previously reported VHE spectra of PSR B1259-63/LS 2883 have a power-law shape with statistically indistinguishable photon indexes (H.E.S.S. Collaboration 2020).

In the context of GRLBs there are three major factors that cause changes of the VHE spectral slope: γ-γ attenuation, anisotropic IC scattering, and changes in the distribution of emitting particles due to the orbital phase. In the specific case of PSR B1259-63/LS 2883, γ-γ attenuation might be relevant only at points close to the periastron passage which, most likely,

has no significant impact during the orbital phases relevant for H.E.S.S. observations in 2021. Similarly, there is no significant change of the scattering angle during this period (see Fig. 7). With these aforementioned factors accounted for, the H.E.S.S. spectral variation measurement implies a hardening of the electron distribution that could, for example, be caused by a change of the cooling regime. If one assumes that the winds interact in an isotropic regime, then the rate of IC and synchrotron losses have a similar dependence (proportional to $R^{-2}$) on the orbital phase. On the other hand, the rate of adiabatic losses scale differently, proportional to $R^{-1}$ (see, e.g., Khangulyan et al. 2008a). Hence, one expects that at large system separations, the transition to an adiabatic loss-dominated cooling regime occurs at higher energies.

If this process indeed defines the hardening of the VHE spectrum, then one should also expect an analogous hardening during similar epochs prior to the periastron passage. The stacked analysis of the H.E.S.S. data collected in 2004, 2007, 2011, 2014, and 2017 indicates that VHE emission during the interval $t_p - 109$ days to $t_p - 47$ days has a photon index of $\Gamma = 2.7 \pm 0.1_{\text{stat}} \pm 0.1_{\text{sys}}$ (H.E.S.S. Collaboration 2020) which is, in fact, significantly softer than the value obtained from “symmetric” orbital phases in 2021 (e.g., $\Gamma = 2.42 \pm 0.1_{\text{stat}} \pm 0.1_{\text{sys}}$ for the dataset C). A complicating factor is that during the periastron passage period the IC scattering angle is larger (see in Fig. 7) and the resulting VHE spectrum is expected to be softer (see, e.g., Khangulyan et al. 2008b).

In summary, it appears that the observed spectral change can be explained in the context of a hardening of the electron spectrum. This is, in turn, driven by changes in the scaling of cooling timescales as a result of varying orbital separation. A detailed numerical model is required to quantitatively test the viability of such a scenario and is beyond the scope of this paper. The possible important role of adiabatic losses supports the general conclusion that in binary pulsar systems, (magneto)hydrodynamic processes play an essential role for non-thermal radiation formation (Bogovalov et al. 2008, 2012, 2019; Bosch-Ramon et al. 2012). Hydrodynamic processes may also lead to the formation of several distinct production regions (Zabalza et al. 2013; Dubus et al. 2015; Huber et al. 2021), and the existence of these seem to be supported by observational evidence. In particular we note the lack of a firm correlation between X-ray and radio emission (Chernyakova et al. 2021), and the very different properties of the GeV and TeV emission detected from the system. A similar absence of correlation between GeV and TeV emission was seen during previous periastron passages (H.E.S.S. Collaboration 2020).

## 5. Conclusions

This work summarises the results from the H.E.S.S. observations and analysis of the 2021 periastron of the PSR B1259-63/LS 2883 system in the VHE band. As displayed in Table 2, our spectral studies reveal that the periastron-averaged spectrum can be described by a power-law model, with a spectral index of $\Gamma = 2.75 \pm 0.05_{\text{stat}} \pm 0.1_{\text{sys}}$. This value is consistent with the average value reported in previous periastron passages (H.E.S.S. Collaboration 2020). We find that the fit has no preference to a power-law containing a cut-off component, with a lower limit on the cut-off energy of $E_{\gamma}^{95\%} = 27.1$ TeV. We also present, for the first time, evidence of spectral variability on a sub-orbital scale. A difference of $\Delta\Gamma = 0.56 \pm 0.18_{\text{stat}} \pm 0.10_{\text{sys}}$ (at greater than 95% c.i.) is seen between the spectral slopes of datasets B and D, see Table 2. Since during the epochs corresponding to the datasets B and D, the γ-γ absorption is neg-
ligible and the change of the IC scattering angle is small, the
revealed hardening indicates on a change of the energy distribu-
tion of the emitting particles, which can be caused by a change
of the cooling regime.

The study of contemporaneous X-ray and TeV fluxes allowed the establishment of a linear correlation between the two energy bands. While the majority of the dataset is fitted relatively well by the applied linear fit (see in Fig. 5), two data pairs show significantly higher X-ray flux levels. The two outliers correspond to orbital phases when the pulsar likely interacts with the disc, therefore the structure of the flow deviates considerably from an axially-symmetric configuration. During this period, it is expected that the pulsar wind terminates at a significantly smaller distance, thereby strongly enhancing the magnetic field.

Regarding the TeV data taken during the time period of the third X-ray peak, we argue that there is good evidence for a correspondence of this TeV data to the third X-ray peak, in the 2021 TeV light curve. However, the time of the maximum in the TeV light curve is not well constrained because of a lack of data 35–55 days post-periastron. Nevertheless, this feature is very interesting and requires further investigation.

The correlation obtained contains a significant constant term, which implies a presence of X-ray emitting electrons with no proportional TeV component. This supports the existence of a multiple emission zone geometry within the system. The evidence for a multi-zone setup can also be obtained from the uncorrelated radiation in the GeV and TeV energy bands, as well as from the absence of a strong X-ray-radio correlation. The formation of a multi-zone setup can originate as a result of the complexity of the hydrodynamics within the pulsar and stellar wind interaction. The detection of spectral hardening at TeV energies after the 2021 periastron passage, together with the measured X-ray-TeV correlation, provide new constraints that will contribute to building a consistent physical model for the multiwavelength emission from PSR B1259-63/LS 2883.

Acknowledgements. The support of the Namibian authorities and of the Uni-
versity of Namibia in facilitating the construction and operation of H.E.S.S. is
gratefully acknowledged, as is the support by the German Ministry for Educa-
tion and Research (BMBF), the Max Planck Society, the Helmholtz Association,
the French Ministry of Higher Education, Research and Innovation, the Cen-
tre National de la Recherche Scientifique (CNRS/IN2P3 and CNRS/INSU), the
Commissariat à l'énergie atomique et aux énergies alternatives (CEA), the U.K.
Science and Technology Facilities Council (STFC), the Irish Research Council
(IRC) and the Science Foundation Ireland (SFI), the Polish Ministry of Educa-
tion and Science, agreement no. 2021/WK/06, the South African Department
of Science and Innovation and National Research Foundation, the University
of Namibia, the National Commission on Research, Science & Technology of
Namibia (NCRST), the Austrian Federal Ministry of Education, Science and
Research and the Austrian Science Fund (FWF), the Australian Research Coun-
cil (ARC), the Japan Society for the Promotion of Science, the University of
Amsterdam and the Science Committee of Armenia grant 21AG-1C085. We
appreciate the excellent work of the technical support staff in Berlin, Zeuthen,
Heidelberg, Palaiseau, Paris, Saclay, Tübingen and in Namibia in the construc-
tion and operation of the equipment. This work benefited from services pro-
vided by the H.E.S.S. Virtual Organisation, supported by the national resource
providers of the EGI Federation.

References

Abdo, A. A., Ackermann, M., Ajello, M., et al. 2009a, Science, 325, 840

Abdo, A. A., Ackermann, M., Ajello, M., et al. 2009b, ApJ, 707, 1310

Abdo, A. A., Ackermann, M., Ajello, M., et al. 2011, ApJ, 736, L11

Acero, F., Ackermann, M., Ajello, M., et al. 2013, ApJ, 773, 77

Aharonian, F., Akhperjanian, A. G., Aye, K. M., et al. 2005a, A&A, 442, 1

Aharonian, F., Akhperjanian, A. G., Aye, K. M., et al. 2005b, Science, 309, 746

Aharonian, F., Akhperjanian, A. G., Bazer-Bachi, A. R., et al. 2006, ApJ, 636, 777

Aharonian, F., Akhperjanian, A. G., Anton, G., et al. 2009, A&A, 507, 389

Barkov, M. V., & Bosch-Ramon, V. 2016, MNRAS, 456, L64

Bausch, J. 2013, J. Phys. A: Math. Theor., 46, 505202

Berge, D., Funk, S., & Hinton, J. 2007, A&A, 466, 1219

Bогоvalov, S. V., Khangulyan, D. V., Koldoba, A. V., Ustyugova, G. V., & Aharonian, F. A. 2008, MNRAS, 387, 63

Bогоvalov, S. V., Khangulyan, D., Koldoba, A. V., Ustyugova, G. V., & Aharonian, F. A. 2012, MNRAS, 419, 3426

Bогоvalov, S. V., Khangulyan, D., Koldoba, A., Ustyugova, G. V., & Aharonian, F. 2019, MNRAS, 490, 3601

Bosch-Ramon, V., Barkov, M. V., Khangulyan, D., & Perucho, M. 2012, A&A, 544, A59

Caliandro, G. A., Cheung, C. C., Li, J., et al. 2015, ApJ, 811, 68

Chang, Z., Zhang, S., Chen, Y.-P., et al. 2021, Universe, 7, 472

Chernyakova, M., Abdo, A. A., Neronov, A., et al. 2014, MNRAS, 439, 432

Chernyakova, M., Neronov, A., van Soelen, B., et al. 2015, MNRAS, 454, 1358

Chernyakova, M., Malyshev, D., Paizis, A., et al. 2019, A&A, 631, A177

Chernyakova, M., Malyshev, D., van Soelen, B., et al. 2021, Universe, 7, 242

Cominsky, L., Roberts, M., & Johnston, S. 1994, ApJ, 427, 978

de Neurois, M., & Rolland, L. 2009, Astropart. Phys., 32, 231

Dubus, G. 2006, A&A, 451, 9

Dubus, G. 2013, A&ARv, 21, 64

Dubus, G., & Cerutti, B. 2013, A&A, 557, A127

Dubus, G., Lambert, A., & Fromang, S. 2015, A&A, 581, A27

Gaia Collaboration (Brown, A. G. A., et al.) 2018, A&A, 616, A1

H.E.S.S. Collaboration (Abramowski, A., et al.) 2012, A&A, 548, A46

H.E.S.S. Collaboration (Abramowski, A., et al.) 2013, A&A, 551, A94

H.E.S.S. Collaboration (Abdalla, H., et al.) 2020, A&A, 633, A102

Huber, D., Kissmann, R., Reimer, A., & Reimer, O. 2021, A&A, 646, A91

Johnston, S., Manchester, R. N., Lyne, A. G., et al. 1992, ApJ, 387, L37

Johnston, S., Manchester, R. N., Lyne, A. G., Nicastro, L., & Spyromilio, J. 1994, MNRAS, 268, 430

Johnston, S., Manchester, R. N., Lyne, A. G., et al. 1996, MNRAS, 279, 1026

Johnston, S., Ball, L., Wang, N., & Manchester, R. N. 2005, MNRAS, 358, 1069

Johnson, T. J., Wood, K. S., Kerr, M., et al. 2018, ApJ, 863, 27

Kelly, B. C. 2007, ApJ, 665, 1489

Khangulyan, D., Hnatic, S., Aharonian, F., & Bогоvalov, S. 2007, MNRAS, 380, 320

Khangulyan, D. V., Aharonian, F. A., Bогоvalov, S. V., Koldoba, A. V., & Ustyugova, G. V. 2008a, IJMPD, 17, 1909

Khangulyan, D., Aharonian, F., & Bosch-Ramon, V. 2008b, MNRAS, 383, 467

Khangulyan, D., Aharonian, F. A., Bогоvalov, S. V., & Ribó, M. 2012, ApJ, 752, L17

Khangulyan, D., Aharonian, F. A., & Kelner, S. R. 2014, ApJ, 783, 100

Kirk, J. G., Ball, L., & Skjæråasen, O. 1999, Astropart. Phys., 10, 31

Li, T. P., & Ma, Y. Q. 1983, ApJ, 272, 317

Marino, A., Driessen, L., Lenc, E., et al. 2023, ATel, 15942, 1

Miller-Jones, J. C. A., Deller, A. T., Shannon, R. M., et al. 2018, MNRAS, 479, 4849

Negueruela, I., Ribó, M., Herrero, A., et al. 2011, ApJ, 732, L11

Parsons, R., & Hinton, J. 2014, Astropart. Phys., 56, 26

Piron, F., Djannati-Atai, A., Punch, M., et al. 2001, A&A, 374, 895

Shannon, R. M., Johnston, S., & Manchester, R. N. 2014, MNRAS, 437, 3255

Sushch, I., & Böttcher, M. 2014, J. High Energy Astrophys., 3, 18

Sushch, I., & van Soelen, B. 2017, ApJ, 837, 175

Sushch, I., & van Soelen, B. 2023, ApJ, 959, 30

Takata, J., Okazaki, A. T., Nagataki, S., et al. 2012, ApJ, 750, 70

Tam, P. H. T., Huang, R. H. H., Takata, J., et al. 2011, ApJ, 736, L10

Tavani, M., & Arons, J. 1997, ApJ, 477, 439

Tavani, M., Arons, J., & Kaspi, V. M. 1994, ApJ, 433, L37

van Soelen, B., & Meintjes, P. J. 2015, Mem. Soc. Astron. It., 86, 123

Weng, S.-S., Qian, L., Wang, B.-J., et al. 2022, Nat. Astron., 6, 698

Wood, K. S., Johnson, T., Ray, P. S., et al. 2018, Am. Astron. Soc. Meet. Abstr., 231, 233,04

Yoneda, H., Makishima, K., Enoto, T., et al. 2020, Phys. Rev. Lett., 125, 111103

Zabalza, V., Bosch-Ramon, V., Aharonian, F., & Khangulyan, D. 2013, A&A, 551, A17

<footnote n="1">¹ Dublin Institute for Advanced Studies, 31 Fitzwilliam Place, Dublin 2, Ireland</footnote>

<footnote n="2">² Max-Planck-Institut für Kernphysik, PO Box 103980, 69029 Heidelberg, Germany</footnote>

<footnote n="3">³ Yerevan State University, 1 Alek Manukyan St, Yerevan 0025, Armenia</footnote>

<footnote n="4">⁴ Landessternwarte Universität Heidelberg, Königstuhl, 69117 Heidelberg, Germany</footnote>
⁵ Kapteyn Astronomical Institute, University of Groningen, Landleven 12, 9747 AD Groningen, The Netherlands

⁶ Laboratoire Leprince-Ringuet, École Polytechnique, CNRS, Institut Polytechnique de Paris, 91128 Palaiseau, France

⁷ University of Namibia, Department of Physics, Private Bag 13301, Windhoek 10005, Namibia

⁸ Centre for Space Research, North-West University, Potchefstroom 2520, South Africa

⁹ Deutsches Elektronen-Synchrotron DESY, Platanenallee 6, 15738 Zeuthen, Germany

¹⁰ Institut für Physik und Astronomie, Universität Potsdam, Karl-Liebknecht-Strasse 24/25, 14476 Potsdam, Germany

¹¹ Université de Paris, CNRS, Astroparticule et Cosmologie, 75013 Paris, France

¹² Department of Physics and Electrical Engineering, Linnaeus University, 351 95 Växjö, Sweden

¹³ Institut für Physik, Humboldt-Universität zu Berlin, Newtonstr. 15, 12489 Berlin, Germany

¹⁴ Laboratoire Univers et Théories, Observatoire de Paris, Université PSL, CNRS, Université Paris Cité, 5 Pl. Jules Janssen, 92190 Meudon, France

¹⁵ Sorbonne Université, CNRS/IN2P3, Laboratoire de Physique Nucléaire et de Hautes Energies, LPNHE, 4 Place Jussieu, 75005 Paris, France

¹⁶ IRFU, CEA, Université Paris-Saclay, 91191 Gif-sur-Yvette, France

¹⁷ University of Oxford, Department of Physics, Denys Wilkinson Building, Keble Road, Oxford OX1 3RH, UK

¹⁸ Friedrich-Alexander-Universität Erlangen-Nürnberg, Erlangen Centre for Astroparticle Physics, Nikolaus-Fiebiger-Str. 2, 91058 Erlangen, Germany

¹⁹ Astronomical Observatory, The University of Warsaw, Al. Ujazdowskie 4, 00-478 Warsaw, Poland

²⁰ Université Savoie Mont Blanc, CNRS, Laboratoire d'Annecy de Physique des Particules – IN2P3, 74000 Annecy, France

²¹ Instytut Fizyki Jądrowej PAN, ul. Radzikowskiego 152, 31-342 Kraków, Poland

²² School of Physics, University of the Witwatersrand, 1 Jan Smuts Avenue, Braamfontein, Johannesburg 2050, South Africa

²³ Laboratoire Univers et Particules de Montpellier, Université Montpellier, CNRS/IN2P3, CC 72, Place Eugène Bataillon, 34095 Montpellier Cedex 5, France

²⁴ School of Physical Sciences, University of Adelaide, Adelaide 5005, Australia

²⁵ Aix Marseille Université, CNRS/IN2P3, CPPM, Marseille, France

²⁶ Institut für Astronomie und Astrophysik, Universität Tübingen, Sand 1, 72076 Tübingen, Germany

²⁷ Universität Innsbruck, Institut für Astro- und Teilchenphysik, Technikerstraße 25, 6020 Innsbruck, Austria

²⁸ Obserwatorium Astronomiczne, Uniwersytet Jagielloński, ul. Orla 171, 30-244 Kraków, Poland

²⁹ Institute of Astronomy, Faculty of Physics, Astronomy and Informatics, Nicolaus Copernicus University, Grudziadzka 5, 87-100 Torun, Poland

³⁰ Department of Physics, Rikkyo University, 3-34-1 Nishi-Ikebukuro, Toshima-ku, Tokyo 171-8501, Japan

³¹ Nicolaus Copernicus Astronomical Center, Polish Academy of Sciences, ul. Bartycka 18, 00-716 Warsaw, Poland

³² Université Bordeaux, CNRS, LP2I Bordeaux, UMR 5797, 33170 Gradignan, France

³³ Department of Physics, University of the Free State, PO Box 339, Bloemfontein 9300, South Africa

³⁴ GRAPPA, Anton Pannekoek Institute for Astronomy, University of Amsterdam, Science Park 904, 1098 XH Amsterdam, The Netherlands

³⁵ Kavli Institute for the Physics and Mathematics of the Universe (WPI), The University of Tokyo Institutes for Advanced Study (UTIAS), The University of Tokyo, 5-1-5 Kashiwa-no-Ha, Kashiwa, Chiba 277-8583, Japan

³⁶ Department of Physics, Konan University, 8-9-1 Okamoto, Higashinada, Kobe, Hyogo 658-8501, Japan

³⁷ Universität Hamburg, Institut für Experimentalphysik, Luruper Chaussee 149, 22761 Hamburg, Germany
## Appendix A: The $\eta^2$ parameter

In order to include the uncertainties of both X-ray and TeV data, we utilise a linear combination of $\chi^2$ tests (see e.g. Bausch 2013), denoted here as $\eta^2$, and defined as:

$$ \eta^2(f, \Omega) = \sum_i \frac{(X_i - f(T_i, \Omega))^2}{\delta X_i^2} + \sum_i \frac{(T_i - f^{-1}(X_i, \Omega))^2}{\delta T_i^2} \quad (A.1) $$

where $X_i$ and $T_i$ are the i-th X-ray and TeV flux values from the time-correlated dataset. Accordingly $\delta X_i$ and $\delta T_i$ are the uncertainties of these values. The dependency between X-ray and TeV fluxes was assumed to have the functional form $X = f(T, \Omega)$ with

$\Omega$ standing for the variable parameter(s) of the function $f$. The inverse function $f^{-1}$ is given by: $T = f^{-1}(X, \Omega)$. For an accurate comparison between $\eta^2$ values, we also implement a method of reduced $\eta^2$, named $\bar{\eta}^2$. After summing the reduction of the two constituent $\chi^2$ values, we reduce $\eta^2$ by applying

$$ \bar{\eta}^2 = \frac{\eta^2}{2(N-1)} $$

where $N$ is the number of correlated pairs.

We note that by design the $\eta^2$ test is symmetric with respect to the interchange of the $X \leftrightarrow T$ datasets.