© The Authors 2024

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 ≳ 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 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)¹. 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_sd = 8.2 × 10³⁵ erg s⁻¹ (Johnston et al. 1994), with a characteristic age of 330 kyr (Johnston et al. 2005).

The companion star LS 2883 has a bolometric luminosity of L_* = 2.3 × 10³⁸ erg s⁻¹ (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 10 R_⊙ (0.05 AU) (Negueruela et al. 2011), and its mass is ∼24 M_⊙ (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 satellite (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 periastron passage, the X-ray light curve exhibited a third maximum between ∼t_p + 30 and t_p + 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 asymmetric double peak profile (H.E.S.S. Collaboration 2020). Maxima derived from a Bayesian block analysis of stacked data from previous 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σ and 39.8σ, 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 methodology 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 passage. In Sect. 4 these results are discussed in the context of previous periastron passages and in the context of the unique findings at other wavelengths in 2021. We also present some theoretical interpretation of the findings of this study. Finally, Sect. 5 contains 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 operated under dark conditions. Because of this, H.E.S.S. is not operated 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 unbiased 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 analysis. The analysis presented in this paper used the reflected regions background method, for light curve production and spectral analysis, as well as the ring background method for the creation of maps (see Berge et al. 2007, for further details of these anaylsis methods). Observations were performed using pointing offsets from the source position, all offsets were exclusively 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 contained 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 periastron 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, representing 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 uncertainties 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.

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°.

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 α = 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 system during the periastron passage was undertaken and several spectra were derived. The total spectrum of the available periastron passage data was calculated, and the spectra of key intervals 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 measured (here defined as t_p + 25 days to t_p + 36 days). Additionally, 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. observational 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 = ϕ₀(E/E₀)^−Γ, where Γ represents the photon index of the power law with a normalisation ϕ₀ and a decorrelation energy of E₀. 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σ. 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.

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.

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σ.

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.

We find a spectral index for the total periastron passage of Γ = 2.78  ±  0.05_stat  ±  0.10_sys (for the fixed energy range spectrum) which is consistent with the spectral index of previous years, Γ = 2.76  ±  0.03_stat  ±  0.10_sys (H.E.S.S. Collaboration 2020).

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 ΔΓ = 0.56  ±  0.18_stat  ±  0.10_sys, implying a sub-orbital spectral variation at a c.i. of greater than 95%.

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 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.

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 tp − 16 days and tp + 16 days correspond to the time at which the disc crossing is thought to occur, the dot-dashed red line at tp = 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.

Individual flux points and their uncertainties were calculated using a reference spectrum, in this work this was a a power-law model in the energy range (0.4 − 100.0) TeV. An index of Γ = 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². 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 Γ = 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 Γ = 2.42 (from the spectrum of the emission from dataset D) and Γ = 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.

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.

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 Lsd = 8.2 × 1035 erg s−1 (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 ∼tp  ±  16 days represent the assumed disc crossing times and the line at tp = 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 tp + 16 days, the purple region represents the time frame tp + 16 to tp + 30 days, the cyan region is the time frame tp + 30 to tp + 50 days (coincident with the third X-ray peak reported in 2021 X-ray data), and the orange region corresponds to any time after tp + 50 days.

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 ΔΓ = 0.56  ±  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 GehrelsSwift (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³) 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 η² (where η² is a linear combination of χ² tests, see appendix A for details of the η² 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.38}^{+0.41}$ and $b=0.{50}_{-0.13}^{+0.12}\times {10}^{-11}$ erg cm⁻² s⁻¹ (where F_X and F_TeV are the X-ray and TeV fluxes, respectively). These values of the fit gave a total η² = 105.08 for 26 correlated pairs, resulting in ${\overline{\eta }}^2=2.10$.

Fig. 6. Linear correlations between X-ray-TeV and GeV-TeV data sets. Left panel: linear correlation (FX = aFTeV + 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−11 erg cm−2 s−1; this outlier originates from rapid X-ray variability in the period from periastron up to ∼tp + 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).

To estimate the statistical uncertainties we performed numerical simulations, considering N = 10⁶ 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 ${\overline{\eta }}^2=2.10$ by comparing the number of trials that provided better η² values than the original data, to the number of trials with a worse η². Making this comparison we found a chance probability of 2.27 × 10⁻³.

We therefore conclude that there is a positive correlation between the X-ray and TeV flux during the time periods of the 2021 perisatron 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σ 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.

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 η² = 57.77 that, for 22 correlated pairs, results in a ${\overline{\eta }}^2=1.38$ (compared to a value of ${\overline{\eta }}^2=2.10$ for the unrestricted dataset) and gave a chance probability of 8.3 × 10⁻⁴.

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 ${\overline{\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 ${\overline{\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σ. 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 decretion disc.

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

The termination of the pulsar wind occurs at a distance of R_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:

$$\begin{array}{cc}\hfill B\lesssim B_{\max }& =\sqrt{\frac{L_{\mathrm{sd}}}{c{R}_{\mathrm{ts}}^2}}\hfill \\ \hfill & \approx 3{\left(\frac{L_{\mathrm{sd}}}{8.2\times {10}^{35}\mathrm{erg}{\mathrm{s}}^{-1}}\right)}^{1/2}{\left(\frac{R_{\mathrm{ts}}}{0.1\mathrm{AU}}\right)}^{-1}\mathrm{G},\hfill \end{array}$$(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

$$\begin{array}{cc}\hfill w_{\mathrm{ph}}& =\frac{L_{\ast }}{4\pi c{R}^2}\hfill \\ \hfill & \approx 3\left(\frac{L_{\ast }}{2.3\times {10}^{38}\mathrm{erg}{\mathrm{s}}^{-1}}\right){\left(\frac{R}{1\mathrm{AU}}\right)}^{-2}{\mathrm{erg}\mathrm{cm}}^{-3},\hfill \end{array}$$(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 ∝ 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 γ–γ 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 ≈ 154° (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°. In this case, an emission of energy 1 TeV may be generated by electrons with energy E_TeV ≈ 1.6 TeV; here we adopt a photon field temperature of T_* ≈ 3 × 10⁴ K (Negueruela et al. 2011), and use the approximation from Khangulyan et al. (2014).

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).

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 ≈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⁻² 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:

$$\begin{array}{c}\hfill L_X\propto L_{\mathrm{TeV}}{R}^{1/2}.\end{array}$$(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 (∝ R⁻²) on the orbital phase. On the other hand, the rate of adiabatic losses scale differently, ∝ R⁻¹ (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 Γ = 2.7 ± 0.1_stat ± 0.1_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., Γ = 2.42 ± 0.1_stat ± 0.1_sys for the dataset C). A complicating factor is that during the pre-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 perisatron 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 Γ = 2.75  ±  0.05_stat  ±  0.1_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_C^{95\%}=27.1$ TeV. We also present, for the first time, evidence of spectral variability on a sub-orbital scale. A difference of ΔΓ = 0.56  ±  0.18_stat  ±  0.10_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 negligible and the change of the IC scattering angle is small, the revealed hardening indicates on a change of the energy distribution 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.

Acknowledgments

The support of the Namibian authorities and of the University 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 Education and Research (BMBF), the Max Planck Society, the Helmholtz Association, the French Ministry of Higher Education, Research and Innovation, the Centre 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 Education 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 Council (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 construction and operation of the equipment. This work benefited from services provided by the H.E.S.S. Virtual Organisation, supported by the national resource providers of the EGI Federation.

References

Appendix A: The η² parameter

In order to to include the uncertainties of both X-ray and TeV data, we utilise a linear combination of χ² tests (see e.g. Bausch 2013), denoted here as η², and defined as:

$$\begin{array}{c}\hfill {\eta }^2(f,\mathrm{\Omega })=\sum _i\frac{{(X_i-f(T_i,\mathrm{\Omega }))}^2}{\delta X_i^2}+\sum _i\frac{{(T_i-f^{-1}(X_i,\mathrm{\Omega }))}^2}{\delta T_i^2}\end{array}$$(A.1)

where X_i and T_i are the i-th X-ray and TeV flux values from the time-correlated dataset. Accordingly δX_i and δ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, Ω) with Ω standing for the variable parameter(s) of the function f. The inverse function f⁻¹ is given by: T = f⁻¹(X, Ω). For an accurate comparison between η² values, we also implement a method of reduced η², named ${\overline{\eta }}^2$. After summing the reduction of the two constituent χ² values, we reduce η² by applying

$$\begin{array}{c}\hfill {\overline{\eta }}^2=\frac{{\eta }^2}{2(N-1)}\end{array}$$

where N is the number of correlated pairs.

We note that by design the η² test is symmetric with respect to the interchange of the X ↔ T datasets.

All Tables

All Figures

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°.

In the text

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.

In the text

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.

In the text

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 tp − 16 days and tp + 16 days correspond to the time at which the disc crossing is thought to occur, the dot-dashed red line at tp = 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.

In the text

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 Lsd = 8.2 × 1035 erg s−1 (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 ∼tp  ±  16 days represent the assumed disc crossing times and the line at tp = 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 tp + 16 days, the purple region represents the time frame tp + 16 to tp + 30 days, the cyan region is the time frame tp + 30 to tp + 50 days (coincident with the third X-ray peak reported in 2021 X-ray data), and the orange region corresponds to any time after tp + 50 days.

In the text

Fig. 6. Linear correlations between X-ray-TeV and GeV-TeV data sets. Left panel: linear correlation (FX = aFTeV + 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−11 erg cm−2 s−1; this outlier originates from rapid X-ray variability in the period from periastron up to ∼tp + 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).

In the text

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).

In the text