The radio to GeV picture of PSR B1259-63 during the 2021 periastron passage.

PSR B1259-63 is a gamma-ray binary system with a radio pulsar orbiting an O9.5Ve star, LS 2883, with a period of ∼ 3 . 4 yr. Close to the periastron the system is detected at all wavelengths, from radio to the TeV energies. The emission in this time period is believed to originate from the interaction of LS 2883 and pulsar’s outflows. The observations of 4 periastra passages taken in 2010-2021 show strong correlation of the radio and X-ray light curves with two peaks just before and after the periastron. The observations of the latest 2021 periastron passage reveal the presence of the 3rd X-ray peak and subsequent disappearance of radio / X-ray flux correlation. In this paper we present the results of our optical, radio and X-ray observational campaigns on PSR B1259-63 performed in 2021 accompanied with the analysis of the publicly available GeV Fermi / LAT data. We compare the properties of di ff erent periastron passages, discuss the obtained results and show that they can be explained in terms of the 2-zone emission model proposed by us previously.


INTRODUCTION
PSR B1259-63 is a gamma-ray binary with a 47.76 ms radio pulsar orbiting a O9.5Ve star (LS 2883) with a period of ∼ 1236.7 days in a highly eccentric orbit (e ∼ 0.87) (Johnston et al. 1992b;Negueruela et al. 2001;Shannon, Johnston & Manchester 2014).Based on the parallax data in the Gaia DR2 Archive (Gaia Collaboration et al. 2018) the distance to the system is 2.39 ± 0.19 kpc, which is consistent with the value of 2.6 +0.4  −0.3 kpc reported by Miller-Jones et al. (2018).
PSR B1259-63 was discovered in 1992 during the Parkes Galactic plane survey (Johnston et al. 1992b,a).Further observations revealed that the pulsed radio emission is completely absorbed from around 20 days before to 20 days after periastron, being obscured by the circumstellar disc of the Be star which is inclined to the orbital plane (Wang, Johnston & Manchester 2004). 1he interaction of the pulsar wind with the equatorial wind of the Be star leads to the appearance of unpulsed radio emission about 20 days before the periastron, quickly rising to a flux exceeding the pulsed flux by a factor of several tens, and lasts for at least 100 hundred days after periastron.The light curve of this unpulsed emission shows two peaks associated with the time of pulsar crossing the plane of the circumstellar disk, before and after the periastron (Connors et al. 2002;Johnston et al. 2005).A similar two-peak structure is also observed in X-rays (see e.g.Chernyakova et al. 2006Chernyakova et al. , 2009)).
During the periastron passage, PSR B1259-63 is also detected at high and very high energies (Chernyakova et al. 2021; H. E. S. S. Collaboration et al. 2020).Current H.E.S.S. observations indicate that the TeV light curve may also have a two peak structure, similar to what is observed at radio and X-ray energies (H.E. S. S. Collaboration et al. 2020), but more sensitive observations are needed to confirm this.Hopefully CTA will address this issue in the very near future (Chernyakova et al. 2019).The GeV emission, however, shows a very different behaviour and is characterised by a strong flare with an average luminosity approaching that of the pulsar spin-down luminosity, occurring at least 10 days after the second X-ray peak with a position varying from periastron to periastron with no obvious counterpart at other wavelength (e.g.Chernyakova et al. 2021).
Several models have been proposed so far to explain the observed multiwavelength emission from this system.The broadband emission has been proposed to originate from electrons accelerated in the shock that forms between the pulsar wind and stellar wind and is produced by synchrotron and IC radiation, e.g.Dubus (2006); Bogovalov et al. (2008); Khangulyan et al. (2011); Chernyakova et al. (2014); Chen et al. (2019).
In Chernyakova et al. (2015) the GeV flare was explained as a result of synchrotron cooling of monoenergetic relativistic electrons injected into the system during the disruption of the stellar disc.The comptonization of unshocked pulsar wind particles (Khangulyan et al. 2012) as well as Doppler boosting (Dubus, Cerutti & Henri 2010;Kong, Cheng & Huang 2012) were also suggested to play a role in producing the GeV flares.Yi & Cheng (2017) proposed that the GeV flare can be a result of the transition between the ejector, propeller, and accretor phases.In this model the compact object is working as an ejector all along its orbit and being powered by the propeller effect when it is close to periastron, in a so-called flip-flop state.
Despite the variety of the models mentioned above, none of them can explain or predict the ultra-luminous short-time subflares discovered during 2017 periastron passage.It turned out that the GeV flare consists of a large number of bright sub-flares which were as short as 15 minutes and had luminosities exceeding the pulsar spin-down luminosity by more than an order of magnitude (Johnson et al. 2018;Tam et al. 2018;Chang et al. 2018;Chernyakova et al. 2020Chernyakova et al. , 2021)).To explain these observations Chernyakova et al. (2020) propose a model in which the observed X-ray and TeV emission is generated by the highly relativistic electrons of the pulsar wind strongly accelerated at the apex of the shock formed between the pulsar and stellar winds.The GeV emission, in turn, is produced by the IC emission of weakly accelerated electrons flowing from the system along the shock, with a possible addition of bremsstrahlung emission on the clumps of the stellar wind material which penetrated beyond the shock cone.In this case one can expect to see a magnification of the GeV emission when the cone is looking in the direction of the observer.
To test this model we have organised an intensive multiwavelength campaign to observe the 2021 periastron passage (Chernyakova et al. 2021).This campaign revealed a number of unexpected features in the system, like a significant delay of the GeV flare, the presence of a third X-ray peak and disappearance of the correlation of the radio and X-ray data observed during the second X-ray peak.
In the current paper we present more details on the radio behaviour of the system including data at 9 GHz in addition to the 5.5 GHz presented earlier, the evolution of the pulsed flux and the polarisation properties (Section 2.1).In Section 2.2 we present 2021 SALT optical observations and compare them with the reanalysed 2014 and 2017 data.In Section 2.3 we present the analysis of newly available XMM-Newton data.In addition to that we have reanalysed all available Fermi/LAT data in order to compare the main properties of the GeV flare (mainly spectrum and timing structure) from different periastron passages (Section 2.4).Broad band spectral modeling of the new data is presented in Section 3, and Section 4 is devoted to the discussion of the relation between radio and X-ray components in the source.Finally the discussion and conclusions are presented in Section 5.

Total radio emission from PSR B1259-63
The radio data were obtained with the Australia Telescope Compact Array (ATCA) using the 4 cm observing band, starting from ∼10 days after the 2021 periastron passage.For the initial calibration, we followed the procedure as described in Chernyakova et al. (2021), using the miriad software package, with B1934−638 as the bandpass and flux calibrator and J1322−6532 as the phase calibrator, which we also used for correcting for the instrumental polarization following standard techniques (Sault, Teuben & Wright 1995).The miriad task uvflux was then used to determine the total flux density and the degree of polarization of PSRB1259−63, using op-tions=uvpol and an inner uv-cut of 20 kilo-wavelengths.The inner uv-cut was employed to minimise contamination of the flux from other sources in the field of view and any diffuse Galactic emission.This provided measurements at both 5.5 GHz and 9 GHz (using the full 2 GHz of bandwidth in both cases), as shown in Table 1 and Fig. 1 (left panel).The uncertainties on individual values were obtain from adding in quadrature the rms noise of the visibilties and the nominal absolute flux density uncertainty of 10%.Some total emission values at 9 GHz are missing due to the failure of the calibration on those days (59286.88, 59315.77, 59317.77).In the case of the pulsar flux, it was not detectable during the first three observations, while for some proceeding days, reliable flux densities could not be extracted, either at 5.5 GHz (59290.85),or both 5.5 and 9 GHz (59296.85,59347.78).The first detection of the pulsed emission occurred ∼ 16 days after the periastron.
The radio spectral index of the total flux α (EdN/dE ∝ E α ) between 5.5 and 9 GHz (see middle panel of Fig. 1) was calculated assuming a power law spectrum interpolation between the observed data points: where F 5 (F 9 ) is the total flux density at 5.5 GHz (9 GHz).The statistical uncertainties α err of the derived indexes were obtained from the Gaussian uncertainties propagation: The average value of the spectral index is ∼ −0.76 with the characteristic uncertainty2 of ±0.25 across all days for which it could be reliably measured, confirming the nature of the radio emission as optically thin synchrotron emission.
For optically thin synchrotron emission, the degree of linear polarization (p) depends on how uniform the magnetic field distribution is across the emission region, with a maximum possible value of ∼72%.From Fig. 1 (left panel), we can see that the polarization's degree at 5.5 GHz p 5 increases substantially when the pulsar emerges from behind the disk, going from ≲5% to ∼25%.It is also notable that p 5 at 5.5 GHz is systematically lower than p 9 at 9 GHz, see right panel of Fig. 1.This decrease in the degree of polarization towards lower frequencies can be ascribed to the effect In blue we show photon indexes that would explain the overall radio to X-ray spectrum (using data points separated by no more than 1.5 days).Right: Ratio of degrees of polarization at 9 GHz and 5.5 GHz.One outlier happened 34 days after the periastron is not shown at a value of 10.7. of Faraday rotation on scales smaller than angular resolution of the observations (Faraday rotation rotates the linear polarization angle proportional to the square of the wavelength, which means that the effect is stronger at lower frequencies).In general, this is known as Faraday depolarization, and provides a way to probe variations in the free electron density and magnetic field structure on scales smaller than the size of the emission region (Burn 1966).
One way to estimate the density and magnetic field variations is by using a model of external Faraday dispersion (Burn 1966), described by the equation where p 0 is the intrinsic degree of polarization (i.e. at zero wavelength), and σ RM is the standard deviation of the Faraday rotation measure within the telescope beam, due to the turbulent magnetic field in the ionized gas (Knuettel et al. 2019) where n e is the electron number density in cm −3 , B the magnetic field strength in G, L the path length through the Faraday rotating medium, and d the scale of fluctuations.The ratio of the degree of polarization between 9 and 5.5 GHz, varies from ∼1 (no depolarization) to ∼10 (strong depolarization), with an average value of ∼ 1.7 (see Fig. 1, right panel).Considering the external Faraday dispersion model, a depolarization factor of 1.7 corresponds to a value of σ RM ∼ 135 rad m −2 (Eq.4), which can be explained by n e B of order 0.5 cm −3 G.Such low values of the density and field are unlikely to be produced close to the stellar disk, which in turn suggests that the observed depolarization is produced by an extended region, where the field and density are significantly lower.However, the time-variable depolarization indicates that the effect is not related to the ISM on large-scales, but instead due to an inhomogeneous medium associated with the system, possibly a clumpy, extended wind.We do not expect the pulsar emission to dominate the polarized signal, since depolarization in pulsars is generally not expected to be as large as observed here, due to the compact nature of the pulsar emission region (Sobey et al. 2019).

Pulsed radio emission from PSR B1259−63
The ATCA correlator was configured such that we were able to simultaneously observe in a pulsar-binning mode at 5.5 and 9 GHz, providing an ability to measure the pulsed flux density (i.e. from the pulsar itself), when the pulsar was visible.The initial calibration for the pulsar data followed the same procedure as the continuum.The calibrated data were then dedispersed with the miriad task psrfix, using a dispersion measure of 146.73 cm −3 pc and period 3 of 47.7625 ms, since the ATCA data were not of sufficient time and frequency resolution to derive these values independently.The pulse profile of flux across the 32 bins at each frequency band was then plotted using the miriad task psrplt, in order to determine if the pulsed emission was visible for a particular observation (i.e.on/off).The ATCA pulse profile is consistent with previous observations from the Parkes telescope at higher time resolution. 4o obtain the pulsed flux density, we first subtracted out the continuum baseline (i.e.emission not from the pulsar) using the miriad task psrbl, selecting only the bins in which the pulsar was not visible.After this, we used uvflux on the output file to obtain the pulsed flux density and its uncertainty (in the same manner as was done for the total flux density described above).The flux density of the pulsed emission is shown in Table 1 and Fig. 1 (left panel) for all days for which a reliable measurement could be obtained.Initially the pulsar flux density is less than 10% of the total flux, but beyond MJD 59300 it is a significant fraction, ranging between ∼30 and 50%.We do not have measurements of the pulsar polarized flux, but it is known to be highly linearly polarized (∼70 to 90%) from the study of Connors et al. (2002).This means that beyond MJD 59300 the pulsar polarized flux possibly contributes up to about half of the total polarized flux.We know that it cannot be all of the polarized emission also because we observed polarized flux even when the pulsar is not visible.

Optical Data
Optical observations were obtained with SALT (Buckley, Swart & Meiring 2006) during the 2021 periastron passage, using both the RSS and HRS instruments (Burgh et al. 2003;Bramall et al. 2010).As we noted in Chernyakova et al. (2021), there is a hint that the equivalent width (EW) of the H α emission line was, on average, slightly weaker than it was during the 2014 periastron passage.For completeness we have, therefore, self-consistently remeasured the EW for the previously reported 2014, 2017 and 2021 periastron passages (Chernyakova et al. 2015;van Soelen et al. 2016;Chernyakova et al. 2020Chernyakova et al. , 2021) ) to minimize an systematic differences that may have arisen from difference in the continuum correction or method of measurement used.A few additional points are also included post 80 days after the 2021 periastron which were not included in Chernyakova et al. (2021).All extracted spectra were continuum corrected using a low-order polynomial in the wavelength region around the H α emission line, and the different continuum corrections were checked against the 2021 RSS observations, for consistency.The EW was then measured between λ = 6528 − 6607 Å to avoid other lines in the continuum, and a linear continuum fit was used between these points to remove any residual continuum.The EW was determined by summing over the data points (using a trapezium method) and the uncertainty was estimated following the method in Vollmann & Eversberg (2006), where the signal-to-noise ratio for each spectrum was determined following the method in Stoehr et al. (2008).
In general, the results are consistent with what has previously been reported, however, there are two changes to the 2014 observations.First, the H α EW for the third observation taken with the SAAO 1.9-m (≈ 33 d after periaston) was found to be slightly lower.Second, was noted that for a few of the RSS observations, the peak of the emission line was slightly too bright and the CCD response was starting to become non-linear.This was realized as shorter, 20 s exposure had been taken at the same time, and a comparison showed the core of the emission line was slightly stronger in the short exposures.However, since the first few observations did not show a difference in the line profile, this effect had not previously been noticed.For correctness, we checked each night and if the profile of the normalized spectrum of the short exposure was more than 2.5 per cent larger than the longer exposure (in the line core) we measured the EW from the shorter exposure.
Last, it was noted that the Modified Julian Date was recorded incorrectly by the telescope software during the 2017 observations.The time has been recalculated from the correct date and time, which shifts the time of observations by 0.5 d but does not change any interpretations or conclusions that were previously drawn.
The results are shown in Fig. 2 (bottom panel).This consistent re-analysis shows a remarkable similarity between the three periastron passages.The observations in 2021, show the EW of the H α line is slightly smaller, but the difference is not large.The increase of the line strength also starts at approximately the same time, and begins to decrease around 35 days after periastron.While this time was co-incident with the start of the Fermi-LAT flare in 2014, the flare occurred much later, and no large change in line strength was observed in 2021 around the increase in GeV activity (Fig. 2).

XMM-Newton
Since the work of Chernyakova et al. (2021), four new XMM-Newton observations become available.We have analysed them for completeness, see Table 2.The analysis was performed with the latest XMM-Newton Science Analysis software v.21.0.0.Known hot pixels and electronic noise were removed, and data were fil- Table 2. XMM-Newton observations of PSR B1259-63.The table summarizes the Observational Id, time since periastron and best-fit parameters of the spectrum with an absorbed powerlaw model ("cflux*tbabs*po") -the total flux in 0.3-10 keV range (not corrected for the absorption), the best-fit slope and n H . p-value shows the null hypothesis probability of the spectral fit.χ 2 LC /ndf shows the χ 2 of the fit of the lightcurve of the corresponding observation with a constant.The moment of periastron was assumed to be t p = 59254.867MJD.tered to exclude soft proton flares episodes (following the standard procedures).The spectrum was extracted from a 40 ′′ -radius circle centred at the position of PSR B1259-63 and the background was extracted from a nearby source-free region of 80 ′′ radius.The RMFs and ARFs were extracted using the RMFGEN and ARFGEN tools, respectively.The obtained spectra of MOS and PN cameras were fit simultaneously with the absorbed powerlaw model with the help of XSPEC software v.12.13.0c (cflux*Tbabs*po model in terms of XSPEC with "wilms" abundances (Wilms, Allen & Mc-Cray 2000)).The best-fit neutral hydrogen density n H , the powerlaw photon index5 Γ (dN/dE ∝ E −Γ ) and 0.3 − 10 keV are summarized in Table 2.We note that off-periastron observations at ∼ 753 d are characterized by the softer slope and lower n H values comparing to the rest of observations taken around the first entrance to the disk.These changes could be connected to additional acceleration of the electrons as the pulsar approaches the disk of Be star and the additional absorption on the neutral hydrogen present in the disk.The obtained results are inline with previous observations reported in Chernyakova et al. (2009).

ObsId
The exposure corrected count rates detected by the PN XMM-Newton camera are shown as a function of time in Fig. 5 for each of the considered observations.The shown light curves are binned with 600 s per point.None of the light curves demonstrate strong variability and are all consistent with a constant level of flux (varying from observation to observation).The χ 2 of the best-fit with a constant for each observation are summarized in the last column of Tab. 2.

Swift/XRT
The Swift/XRT data presented in this work were adopted from Chernyakova et al. (2021).These data include the spectral indexes in 1-10 keV energy band as well as the flux in the same energy band.In order to determine radio-to-X-ray spectral slope we explicitly assumed a power-law spectrum continuing between these energy bands.The radio-to-X-ray index was then determined based on the flux measurements in corresponding energy bands for measurements separated by no more then 1.5 days.The X-ray, radio and radio-to-X-ray spectral indexes are shown in Fig. 1 (middle panel) with red, green and blue points, respectively.The shown errorbars correspond to 1σ uncertainties.The Swift/XRT data similarly to XMM-Newton were fitted with an absorbed powerlaw spectral model (cflux*tbabs*po in XSpec terms) with the n H absorption fixed to 0.7 × 10 22 cm −2 -the value derived from the joint fit of all Swift/XRT observations.Please note, that the lower n H values seen in XMM-Newton observations, see Tab. 2 would correspond to the even harder powerlaw indexes due to the strong correlation of these parameters.

GeV Data
A detailed description of our analysis of the Fermi/LAT data of the PSR B1259-63 during its 2021 periastron passage is presented in Chernyakova et al. (2021).In order to perform a proper comparison of the GeV flares during 2010, 2014, 2017 and 2021 periastron passages we reanalyzed all the data using the same version of the software (Fermitools version 2.2.0), using the latest Pass 8 reprocessed data (P8R3) from the SOURCE event class for all event types (FRONT+BACK) with the latest IRFs (V3).The binned likelihood analysis was performed for photons within the energy range 0.1-100 GeV (summarised by > 0.1 GeV) that arrived around 2010, 2014, 2017 and 2021 periastron passages within a 20 • -radius region around PSR B1259-63 position, based on the analysis procedure recommended by the Fermi collaboration6 (see also Malyshev & Mohrmann 2023).The selected maximum zenith angle was 90 • .The source model for the likelihood analysis included all sources within a 25 • -radius region around PSR B1259-63 using the 4FGL-DR4 catalogue (Abdollahi et al. 2020).For the likelihood fit of the total data for any periastron, the spectral parameters of all sources outside a 5 • -radius were fixed.All other sources only had the normalization value as a free parameter, except PSR B1259-63 which had a free normalization and index.For subsequent light curve and spectral points built with a likelihood fit, the index of PSR B1259-63 was fixed to the previously calculated value.
Initially, light curves were made for each periastron period using aperture photometry.The binning was done so that each time bin accommodates 9 photons with energies above 0.1 GeV in a 1 •radius circle around the position of PSR B1259-63.The resulting time bins have a duration of 5 min to 7 hours with an average duration of about 2 h.In Fig. 4 we show a comparison of the GeV light curves during the flares of different years with the same adaptive binning built through the likelihood analysis of Fermitools to verify the results of the aperture photometry.Upper limits value for both light curves and spectra were calculated with a 95% confidence level using the IntegralUpperLimits python module available in Fermitools.It is clearly seen from these figures that short GeV sub-flares with luminosities exceeding the spin-down luminosity are present after each periastron.The number and luminosity of these sub-flares, however, differs from periastron to periastron.
In order to compare the spectral properties of the sub-flares to the rest of the flare, we built combined spectra from all data bins with a peak flux above the spin-down luminosity and TS > 16 (pkfl), and from all data bins below the spin-down luminosity (lowfl).These spectral parameters, for each periatron passage, are given in Table 3.Also shown are the spectral parameters of the average spectrum over the flare period (avfl), and the average over the pre-flare period calculated for both 20 days before (prfl1) and 20 days after (prfl2) periastron.For ease of comparison we fit all the data with a power-law spectrum.Note however, that at least in 2017 the prfl2 data set was better described with a cut-off power law model (Chernyakova et al. 2020).
Unfortunately, the quality of the data does not allow us to unambiguously interpret the obtained results.The Tab. 4 summarizes the results of the fit of the spectral indexes (see Tab. 3) with a constant for different time periods discussed above.The periods include corresponding periastra averaged (labeled 2010-2021) and specific periods averaged over all observed periastra (labeled pfl1avfl).The other columns indicate the best-fit χ 2 /d.o.f. of the fit of the spectral index with a constant; significance σ of the consistency of the data with constant; the best-fit (average) spectral index and its uncertainty.Indeed, the significance of the variability of the average photon index Γ, during each individual periastron passage is Since the quality of the data are significantly better during the flare state, we have also tried to analyze the average spectral variability during different phases, by determining how well the slopes are fit by a constant value (Table 4).We note that during the preflare (prfl1 and prfl2), and low-flare (lowfl) periods, the spectral indices are consistent with a constant value.However, the peakflare (pkfl) and average flare (avfl) data sets are not consistent with a constant spectral index at a ∼ 3σ and a ∼ 4σ level, respectively (Tab.4).One can also conclude that, on average, the spectrum of the source during the period of GeV flare is softer than before the flare.

MODELLING
In our analysis we have used the same model as was used to describe the data from 2017 periastra passages in Chernyakova et al. (2020).In this model we assume that the observed emission is com-ing from two main populations of relativistic electrons: (i) electrons of the unshocked and weakly shocked pulsar wind (producing GeV emission via IC and/or bremsstrahlung mechanisms) and (ii) strongly shocked electrons accelerated near the apex of the shock (producing X-ray and TeV emission via synchrotron/IC mechanisms), see Fig. 5 of Chernyakova et al. (2020).In this model one can expect to see an enhancement of the GeV emission when the bow shock cone is pointing in the direction of the observer, see Appendix A for more details.Please note that in the discussed 2-zone toy model several potentially important effects (anisotropy of the IC emission; exact geometry of the system) are not explicitly considered.The model is not designed to calculate the exact variability of the emission along the orbit, but rather to illustrate that the observed enhancement of the emission can be understood in general.
In our simplified model we have calculated all the IC emission on the soft photons of the Be star only (T = 27 000 K), and assumed the same photon density in both emission zones.
The spectrum of unshocked electrons was selected to be a power law with a slope of Γ e = −2 in the energy range E e = 0.2 − 1 GeV.A small fraction of electrons are additionally accelerated at the strong shock near the apex to E e ∼ 500 TeV energies with the similar slope Γ e = −2.The rest of the electrons flying into the shock direction are reverted to flow along the shock at the surface of stellar-pulsar wind interaction cone far from the Table 5.Details of the model parameters, see section 3. D is a distance from the Be star to the emission region, Γ 2 is a slope of accelerated electrons above 1 GeV, and t esc is time the electrons stay in the emitting region before the escape.Effective luminosity L of the pulsar wind electrons (without considering beaming effects) is measured in units of spin-down luminosity L sd = 8.2 × 10 35 erg/s.apex, and could be additionally mildly accelerated on a weak shock.This leads to a power law tail in the spectrum of diverted electrons with a softer slope Γ 2 which continues above 1 GeV to at least E e ∼ 5 GeV.A slope which is softer than 2 is a characteristic of particle acceleration on weak shocks (see e.g.Bell 1978;Blandford & Eichler 1987).
After the injection, the spectra of both populations are modified by radiative (Inverse Compton, synchrotron, or bremsstrahlung) and non-radiative (adiabatic or escape) losses operating in the system.In our calculations, the resulting electron spectrum was determined by numerically calculating the radiative losses of a continuously injected spectrum of electrons, until a steady solution has been obtained.The time that electrons spend in the emitting region is about a few thousand seconds, and is long enough to substantially modify the injected spectrum.
In the analysis of the broadband data of 2021 periastron passage we considered separately the spectra during the 20 days before (prfl1) and after (prfl2) the periastron passage, and during the high (pkfl) and low (lowfl) flaring states.We found that the model describes the data reasonably well, and the contribution of bremsstrahlung emission is important only during the flare period, see Table 5 and Fig. 6.We note that for the modelling of the flaring state we have assumed the escape time to be rather short, t esc = 1000 s in order to match to the duration of the shortest subflares discussed in section 2.4.We would like to note also a strong correlation between modeled clumps' density n clump and the escape time from the system t esc .In order to maintain the total energy in the relativistic electrons present in the system to be roughly the same, the longer t esc requires lower values of n clump .
Note that all results discussed above are strongly model dependent and are used here as general illustration of the model.Moreover, the parameters of the model (magnetic filed strength, soft photon field density, electrons spectrum, escape time) are strongly correlated and form a degenerate parameter space.Thus the observations can be explained via somewhat different model parameters.More accurate modelling taking into account the full geometry of the system is needed for a better reconstruction of system parameters.

RADIO TO X-RAY BEHAVIOR
The relation between the radio and X-ray components in PSR B1259-63 has a non-trivial behavior.As it was reported in Chernyakova et al. ( 2021) at first the radio observations corresponding to the time of the second X-ray peak (from ∼ 11 to 27 days after the 2021 periastron passage) show a clear correlation with the X-ray flux.After this time the radio to X-ray correlation stops, see Fig. 3. Using the data obtained during the correlation period we performed a search for a possible delay between the X-ray and radio emission (expected, for example, in the case when the emission in these bands is produced in spatially separated regions).For this purpose we selected the time interval of 11-27 days after the periastron passage (white region in Fig. 3) and re-normalized the radio-band light curve to match the average X-ray flux during this time period.We modelled both the X-ray and re-normalized radio light curves with a linear increase -constant -linear decrease of the flux model, see dotted black curve in Fig. 3.In addition, we suggested that exactly the same model describes both X-ray and radio data sets except that the model describing the radio data is delayed in time in relation to the X-ray data.We fit this model to both the X-ray and radio data sets simultaneously, assuming the flux increase/decrease rates, the start/end times of the constant flux level and the radio-X-ray flux delay to be free parameters.The best-fit of the model to the data is shown with dotted black line in Fig. 3.We did not find any significant delay between X-ray and radio emissions with a best-fit delay value of 0 ± 0.20 d.
In terms of the spectral-energy distribution, the radio data in general lie on the continuation of the X-ray spectrum (Figures 1  & 6), though the spectral index in the radio band is systematically softer than the index in the X-ray band and the index that would be required to simultaneously fit the radio and X-ray data.
The fact that observations do not show a single power-law from radio to X-rays is not surprising.Within our proposed model the energetics constraints do not allow the population of electrons to produce a single, continuous power-law spectrum of photons from the radio to the X-ray band.
The radio emission could be partially explained as synchrotron emission from the same population of electrons which produces the GeV emission via IC/bremsstrahlung mechanisms.However, the predicted spectral slope in this case significantly differs from the observed one that does not allow us to explain all the observed radio emission in this way.The dashed magenta curve in the left panel of Fig. 6 illustrates the emission in a characteristic magnetic field of ∼ 0.2 G around the periastron.Note that any higher value of the magnetic field in the radio emission region will produce emission at a higher level than the observed one and is thus disfavored by our model.During the GeV flare a two orders of magnitude lower magnetic field is required.Note that the exact value is strongly model dependent and can be higher if one takes into account possible absorption of the radio emission or/and assume higher clump density and attributing most of the GeV emission to bremsstrahlung emission.
We argue that the majority of observed radio emission is produced either by the high-energy electrons of the pulsar wind that have cooled down during their movement along the shock front, or by a different population of electrons, e.g.electrons in the Be-star's wind which are accelerated at the shock.
In the case that the radio emission is produced by the cooled down population of the relativistic electrons, the characteristic properties of the population and the magnetic field can be estimated as follow.The characteristic frequency of the photons emitted by electrons with an energy E e in a magnetic field B is: GHz.
(5) 10 5 10 2 10 1 10 4 10 7 10 10 10 13 E, eV Dashed lines represent the contribution of synchrotron emission, dotted lines the IC emission, and dash-dotted lines the bremsstrahlung emission.Note that the contribution of the IC component in the GeV region is the same for both the low and high states of the flare (dotted blue curve).The radio points F5 and F9 in the left panel correspond to radio emission at 5 and 9.5 GHz 20 days after the periastron, and 50 days after periastron in the right panel.The blue, green and magenta butterflies in the keV range correspond to the range of the observed X-ray emission by SWIFT/XRT 20 days before and after the periastron, and during the period of GeV flare, respectively.The TeV points are taken from (Aharonian et al. 2005), and the shaded regions at TeV energies represent the range of multiyear H.E.S.S. measurements (see Fig. The corresponding cooling time for the same electrons is: The available radio-band observations indicate variability on time scales as short as ∼ 1 day.If the radio emission is produced by the cooled down population of relativistic electrons, this variability time scale should not exceed the cooling time, i.e. t synch ≲ 1 d.Together with this, the requirement that ν synch is equal to the radio frequency of the observations (ν synch ∼ 10 GHz) provides two equations for the two unknowns (E e and B).The solution given by the equations results in unreasonably high values of the magnetic field B ∼ 24 G.Such high values of the magnetic field would indicate that the emission region is located much closer to the star (to keep the ratio of the energy densities of the magnetic and photon fields at the same level), but will not strongly affect the spectral shape of the electron population.The requirement of the high magnetic field could be overcome in the case when the escape time of the electrons is significantly shorter than t synch .In this case, however, the cooling of the electrons is very inefficient and it is difficult to explain that the full energy release in the radio band is roughly consistent with the continuation of the X-ray spectrum.
The introduction of a separate population of electrons does not require high magnetic fields, but instead requires a fine tuning to explain the similarities of the X-ray and radio spectral indexes as well as the similarities of the radio flux and the power-law continuation of the X-ray spectrum.The origin of the initial correlation between radio and X-ray emission and its further disappearance (Chernyakova et al. 2021) is also not clear in this case.

DISCUSSION AND CONCLUSIONS
In the current paper we present the results of the analysis of the radio-to-GeV data for the recent periastron passage of 2021 and systematic re-analysis of all previous GeV (Fermi/LAT) and optical (SAAO/SALT) data.
We present here for the first time the behaviour of PSR B1259-63 at 9 GHz, the radio spectral index, and degree of polarization at both 5.5 and 9 GHz during the 2021 periastron passage (Table 1).The total radio flux is initially well correlated with the X-ray emission (i.e. the second X-ray peak), however the third X-ray peak is uncorrelated with the radio emission, which continues decreasing at both 5 and 9 GHz.The spectral index of the radio emission is consistent with optically thin synchrotron emission, with no evidence for free-free absorption during the observed period.In general, the radio spectrum is inconsistent with a single population of electrons explaining both X-ray and radio emission.The degree of polarization of the radio emission increases substantially (from ≲5 to ∼25%) as the pulsar emerges from behind the disk, as it experiences less depolarization overall and/or the magnetic field in the emission region becomes more ordered.Time-variable depolarization is also observed, which is most likely due to a clumpy, extended wind associated with the system but not close to the stellar disk (where the field and density are too high to explain the observations).
The consistent re-analysis shows a remarkable similarity between the behaviour of the EW of the H α lines during last three periastron passages.Many Be X-ray binary systems show superorbital periods, linked to variations in the circumstellar disc (e.g.Rajoelimanana, Charles & Udalski 2011), and this is found for the gamma-ray binary LS I +61 303 (e.g.Paredes-Fortuny et al. 2015).The lack of large changes in the equivalent width seems to suggest a similar outflow from the star, and a consistent behaviour for the disc over the last few periastron passages.The observed differences in optical ranges are much smaller than the difference in the GeV behavior.While the decrease in line strength was co-incident with the start of the Fermi/LAT flare in 2014, the flare in 2021 occurred much later, with no large change in line strength around the start of GeV activity, see Fig. 2. Similar to the previous report of Chernyakova et al. (2021) we do not find any counterparts of the 2021 GeV flare at any other wavelengths.We note, however, that the binary separation at the time of the GeV flare is ≈ 3 − 5 AU.
Such distances significantly exceed the size of the region where the majority of the H α emission is expected to arise.This suggests that the GeV flare could manifest itself at longer wavelengths (e.g.Klement et al. 2017).Such observations can be more sensitive to the changes in the colder, outer regions of the Be star disk and would be of extreme importance for the identification of GeV flare counterpart at other wavelengths.
In this paper we also present analysis of all recent (after 2021) XMM-Newton observations, see Table 2 and Fig. 5.We note that, similar to previous periastron passages, the spectral slope of the system is significantly harder shortly before the periastron (Γ ∼ 1.3 − 1.4) than near the apastron (Γ ∼ 1.7).We do not see any significant variability on short (10 minutes) time scale, see Fig. 5.At the same time we found that the spectral index in the radio band is systematically softer than the one in the X-ray band, or the index derived from the interpolation between radio and X-ray energies.We also do not find any evidence for a delay between X-ray and radio emission with the best-fit delay to be 0 ± 0.20 d.
A systematic comparison of timing and spectral properties of GeV flares over the last 4 periastron passages demonstrates that all flares differ in their total duration, overall short-timescale structure and spectral properties, see Tables 3 and 4.
Our model describes the X-ray-to-TeV data reasonably well with parameters similar to presented before (Chernyakova et al. 2020), see Table 5.At the same time this model under predicts the radio flux (see Figure 6).The radio data seems to be located on the powerlaw continuation of the X-ray data, although the X-rayto-radio index is systematically harder than radio one, see Figure 1.To explain this one may consider either effective cooling (requires very high magnetic field B∼ 20 G) or a second population of electrons (e.g.electrons of the Be-star's wind accelerated at the shock).The introduction of a second electron population would need a fine tuning to explain the initial radio/X-ray correlation and its subsequent disappearance as well as the fact that the observed radio emission is lying roughly on the continuation of the X-ray spectrum.The assumption of a very high magnetic field would mean that the magnetic field in the emitting region is, at least at some orbital phases, dominated by the magnetic field of the Be star (as assumed in e.g.Melatos, Johnston & Melrose 1995), which can lead to interesting consequences, e.g.correlation of X-ray/TeV emission (as both magnetic and soft photon fields are proportional to the distance from the star).i.e.Ω cl ∼ 1/30π for the typical values of the assumed parameters (α ∼ 30, γ ∼ 3).We note also that the discussed condition Ω emi > Ω cl indeed holds as Ω emi ≈ π/γ 2 > 1/30π = Ω cl .
We note also, that in case if the clump penetrates through both shock waves and interacts with the electrons of the unshocked wind, then the flux of electrons intersected by the clump from the isotropic central source would be F cl,el = F iso,el Ω cl 4π The observed photon flux from the clump (under the assumption that in this case Ω cl > Ω emi ) is equal to This means that in this case the flux of the photons detected by the observer cannot exceed the spin-down luminosity.

Figure 1 .
Figure1.Left: Evolution of the total flux, pulsed flux and degree of polarization at 5.5 GHz.Middle: Evolution of X-ray (red) and Radio (green, 1 − α) photon indexes around the periastron.In blue we show photon indexes that would explain the overall radio to X-ray spectrum (using data points separated by no more than 1.5 days).Right: Ratio of degrees of polarization at 9 GHz and 5.5 GHz.One outlier happened 34 days after the periastron is not shown at a value of 10.7.

Figure 2 .
Figure2.Evolution of GeV emission (top panel) and of the equivalent width of the H α emission (bottom panel) as measured over the last three periastron passages.Fermi/LAT flux measurements are given in the E > 100 MeV energy range with a weekly bin size taken from(Chernyakova et al. 2021).Flux is given in 10 −6 cm −2 s −1 .

Figure 3 .
Figure3.The X-ray (red) and radio (green) light curves of PSR B1259-63 for 11 -40 days after the 2021 periastron passage.The light-grey shaded region shows the time when the radio/X-ray flux correlation starts to disappear.X-ray points are taken fromChernyakova et al. (2021).The dotted black curve illustrates the model used to describe the variations of the flux during the 2nd X-ray peak, see text for the details.

Figure 4 .Figure 5 .
Figure 4.The GeV flares during the 2010, 2014, 2017 and 2021 periastron passages with adaptive time binning (see text).Each light curve covers the period when the highest GeV flux was detected for each periastron passage.Data points are shown in blue circles while upper limits are shown in orange squares.The horizontal dashed line indicates the isotropic flux level corresponding to the spin-down luminosity (L sd = 8.2 × 10 35 erg/s).

Figure 6 .
Figure6.The GeV emission of PSR B1259-63 during (left: ) twenty days before (blue points, prfl1) and after (red points, prfl2) the periastron, and (right:) during the 2021 flaring period.Magenta points correspond to the high state of the flare (2021-pkfl) and blue points to the low state during the flare (2021-lowfl).Dashed lines represent the contribution of synchrotron emission, dotted lines the IC emission, and dash-dotted lines the bremsstrahlung emission.Note that the contribution of the IC component in the GeV region is the same for both the low and high states of the flare (dotted blue curve).The radio points F5 and F9 in the left panel correspond to radio emission at 5 and 9.5 GHz 20 days after the periastron, and 50 days after periastron in the right panel.The blue, green and magenta butterflies in the keV range correspond to the range of the observed X-ray emission by SWIFT/XRT 20 days before and after the periastron, and during the period of GeV flare, respectively.The TeV points are taken from(Aharonian et al. 2005), and the shaded regions at TeV energies represent the range of multiyear H.E.S.S. measurements (see Fig.2of H. E. S. S.Collaboration et al. 2020).
Figure6.The GeV emission of PSR B1259-63 during (left: ) twenty days before (blue points, prfl1) and after (red points, prfl2) the periastron, and (right:) during the 2021 flaring period.Magenta points correspond to the high state of the flare (2021-pkfl) and blue points to the low state during the flare (2021-lowfl).Dashed lines represent the contribution of synchrotron emission, dotted lines the IC emission, and dash-dotted lines the bremsstrahlung emission.Note that the contribution of the IC component in the GeV region is the same for both the low and high states of the flare (dotted blue curve).The radio points F5 and F9 in the left panel correspond to radio emission at 5 and 9.5 GHz 20 days after the periastron, and 50 days after periastron in the right panel.The blue, green and magenta butterflies in the keV range correspond to the range of the observed X-ray emission by SWIFT/XRT 20 days before and after the periastron, and during the period of GeV flare, respectively.The TeV points are taken from(Aharonian et al. 2005), and the shaded regions at TeV energies represent the range of multiyear H.E.S.S. measurements (see Fig.2of H. E. S. S.Collaboration et al. 2020).

Table 1 .
ATCA observations of PSR B1259-63 in 2021.The table summarizes the modified Julian date of the observation (MJD),

Table 3 .
Spectral parameters around 2010Spectral parameters around  , 2014Spectral parameters around  , 2017Spectral parameters around   and 2021periastron passages when fitted with a power-law model (dN/dE ∝ E −Γ ).The second column gives either the dates of the observations relative to periastron or the total length of the state if it is shown in brackets.The last column gives values of test statistics of the source detection during the given period.For the definition of different periods see the text.

Table 4 .
The χ 2 value for each period if the Fermi indexes are fit with a constant value.The significance of the data not being described by a constant, is given by σ.