PeV photon and neutrino flares from galactic gamma-ray binaries

The high-energy radiation from short period binaries containing a massive star with a compact relativistic companion was detected from radio to TeV gamma rays. We show here that PeV regime protons can be efficiently accelerated in the regions of collision of relativistic outflows of a compact object with stellar winds in these systems. The accelerated proton spectra in the presented Monte Carlo model have an upturn in the PeV regime and can provide very hard spectra of sub-PeV photons and neutrinos by photo-meson processes in the stellar radiation field. Recent report of a possible sub-PeV gamma-ray flare in coincidence with a high-energy neutrino can be understood in the frame of this model. The gamma-ray binaries may contribute substantially to the Galactic component of the detected high-energy neutrino flux.


INTRODUCTION
The cosmic accelerators of petaelectronovolt (PeV) energy particles revealed themselves through the measured fluxes of the Galactic cosmic rays as well as by the high-energy neutrinos first detected by IceCube (Aartsen et al. 2013;IceCube Collaboration 2013). Ground based γ-ray telescopes found a population of sources with power-law spectra extending above 10 TeV without a clear signature of the spectral cut-off (Aharonian et al. 2019;HESS Collaboration et al. 2016) which can be associated with cosmic pevatrons. Strong attenuation of γ-ray fluxes above 100 TeV due to photon-photon pair production limits the possibility of the γ-ray pevatron detection mainly to a population of nearby Galactic sources at distance D < ∼ 10 kpc (Nikishov 1962;Dermer & Menon 2009). The Large High Altitude Air Shower Observatory (LHAASO) collaboration reported a significant detection of 12 γ-ray sources above 100 TeV and up to 1.4 PeV highlighting the problem of the origin of these PeV accelerators (Cao et al. 2021). The Carpet-2 experiment team (Dzhappuev et al. 2021) recently reported a detection of a 3.1σ excess of γ-ray flux above 300 TeV from the Cygnus region associated with a 150-TeV neutrino event detected by IceCube (Lagunas Gualda et al. 2020) and most likely consistent with a flare of a few months duration. The γ-ray energy flux during the flare is 10 −9 erg cm −2 s −1 , an order of magnitude higher than the 95% CL upper limit on the steady-state flux, 1.2·10 −10 erg cm −2 s −1 , obtained for the same source by Carpet-2. The flare flux is well above the fluxes in the TeV regime detected from the population of γ-ray sources in the region (see e.g. Amenomori et al. 2021). It also highly exceeds the fluxes from known γ-ray sources including the γ-ray pulsars, supernova remnants (e.g. Tibet ASγ Collaboration et al. 2021), superbubbles (e.g. Abeysekara et al. 2021), γ-ray binaries (e.g. Aharonian et al. 2006) as well as from yet unidentified PeV candidate sources (e.g. The LHAASO collaboration 2021; Abdalla et al. 2021).
Most of high-energy neutrinos are probably extragalactic. Combined IceCube and ANTARES data (Albert et al. 2018) limit the Galactic contribution to 15%, and first indications to the presence of this component have been found (Aartsen et al. 2019). Some Galactic contribution helps to explain naturally (Palladino & Vissani 2016;Neronov & Semikoz 2016) the tension between different measurements of the neutrino spectrum, see e.g. Abbasi et al. (2021). For a review of particular classes of Galactic sources see Kheirandish (2020). Gamma-ray binaries were proposed as high-energy neutrino sources long ago (see e.g. Levinson & Waxman 2001;Distefano et al. 2002;Bednarek 2005;Neronov & Ribordy 2009;Zdziarski et al. 2010;Sahakyan et al. 2014).
The angular resolution of both Carpet-2 and IceCube does not allow one to identify a particular source of the photon flare and of the contemporaneous neutrino in the densely populated Cygnus region. The time variability on a few months timescale and the high observed γ-ray flux make any associations of the source with extended supernova remnants or superbubbles very unlikely. A number of compact Galactic sources of high-energy radiation are located in the Cygnus region, including gamma-ray binaries Cyg X-3 and PSR J2032+4127.
The maximum energies of protons accelerated by outflows with frozen-in magnetic fields of a kinetic/magnetic luminosity L K can be estimated from the equation: where the dimensionless velocity of the flow is β f = u f /c, c is the speed of light, Γ f = 1/ 1 − β 2 f and Ω is the opening angle of the outflow (see e.g. Lemoine & Waxman 2009;Bykov et al. 2012, and references therein).
It follows from the equation that for a given L K , the higher values of E max can be achieved for the mildly relativistic flows with β f Γ f ∼ 1 with relatively narrow opening angle Ω < 1. Then L K > 3 × 10 35 erg s −1 and Ω < 1/3 are needed to reach the energy of the accelerated proton ∼ 10 PeV. Eq. (1) does not account for radiative losses, which can reduce E max substantially.
The observed flare flux transforms to the γ-ray luminosity above 300 TeV of L γ ∼ 4 × 10 35 (D/1.5 kpc) 2 erg s −1 . Thus, if the source is indeed located in the Cygnus star-forming region then, according to Eq. (1), L K required to produce PeV regime particles radiating 300-TeV photons is comparable to L γ .
One can conclude that a very hard spectrum of accelerated particles and a fast cooling of PeV particles are needed in order the required kinetic/magnetic luminosity to be consistent with that available in compact relativistic sources of stellar masses. We present below a model of the compact γ-ray sources that can convert a substantial fraction of their kinetic luminosity (provided by the magnetic braking of pulsar, magnetic field reconnection in magnetars, or accretion onto black hole) to the PeV regime γ-rays and neutrinos by photo-meson production mechanism in proton-photon collisions. Gamma-ray binary sources (LS 5039, PSR B1259-63, LSI +61°303, PSR J2032+4127, and others, see e.g. Dubus (2013)) can be considered as possible candidates. Well-known powerful microquasar Cyg X-3 demonstrated giant γ-ray flares (Corbel et al. 2012). In all of these sources, which have been subjects of extensive modeling for a long time (see e.g., Tavani & Arons 1997;Chernyakova et al. 2020b), the compact object has a massive early type star companion.
The model of a PeV source we discuss here suggests the interaction of a fast outflow from a compact relativistic object with the stellar wind (SW) of a bright massive star. The colliding magnetized flows provide a plausible way to a fast Fermi-type acceleration to the PeV regime of TeV-energy particles pre-accelerated in the vicinity of the compact object. The acceleration mechanism at the colliding wind flows (CWFs) may form very hard spectrum of particles in the TeV-PeV energy range. In addition to γ-ray radiation produced in proton-proton collisions, the accelerated PeV protons interact efficiently with the optical photons of the luminous massive star by photo-meson mechanism providing fast cooling in PeV γ-rays and neutrinos. In the next section we will discuss the PeV γ-ray and neutrino production in the generic source PSR J2032+4127, though a similar model could be applied to other binary sources including either pulsar winds (PWs) or jets of the accreting black holes.

MODEL
PSR J2032+4127 is located at distance ∼ 1.4 kpc and orbits around a massive Be star MT91 213 (B0Vp) with a long period ∼ 50 yr (Lyne et al. 2015;Ho et al. 2017). Its spin-down powerĖ may reach 3×10 35 erg s −1 (Camilo et al. 2009). This is close to the γ-ray luminosity ∼ 2.7×10 35 erg s −1 at the lower limit of Carpet-2 measurement uncertainty band. This suggests that most of the CWFs acceleration source power should be converted into PeV range protons energy. The total available power of the CWF is a sum of the pulsar spindown power and the fraction of the SW power released in the CWFs. Fermi I type acceleration in the CWFs can produce particle energy distribution f (E) ∝ E −s with s ∼ 1, where energy is mainly accumulated by the most energetic particles (see Bykov et al. 2017). Efficient acceleration at CWFs may be caused by passage of the pulsar through the equatorial region of the Be star SW. In November 2020 the pulsar was at ∼ 20 AU from MT91 213 that is about ∼ 400 stellar radii allowing interaction of the PW with the Be star equatorial disk (Klement et al. 2017).
The model of PeV flares of PSR J2032+4127 we discuss here uses only the hadronic emission processes. Severe synchrotron radiation losses do not allow acceleration of PW leptons to PeV energies at considered orbital phase. The toroidal SW magnetic field component scaling ∝ r −1 gives ∼ 1 G at distance r ∼ 20 AU, if the surface dipolar stellar field is ∼ 1 kG. For a fast rotating star with mass M * ∼ 15M and radius R * ∼ 10R , even higher field ∼ 2 kG may be expected (see, e.g., Shultz et al. 2019).
We simulate acceleration at the CWFs using the kinetic Monte Carlo model described in Section 4 of Bykov et al. (2017) adapted for this problem. This model allows simulating the diffusive particle propagation in the region where the PW collides with the ambient matter flow, launching the PW termination shock and possibly the bow shock. The model plausibly catches the spatial structure of CWFs system represented by spatial zones with reasonable parameters of diffusion and magnetic fields and includes a simple model of flows velocities distribution (see Fig. 1).
Protons are injected in the CWFs region with a wide soft power-law spectrum f inj (E) ∝ E −s with s ∼ 2.2 (e.g. Bykov et al. 2012;Sironi et al. 2015;Amato 2020). It is produced by Fermi I type acceleration at the termination shock formed in the collision of the pairdominated pulsar wind with magnetized plasma of stellar wind. Both pairs and protons have to be accelerated in this system. In the simulation we generate a population of particles with the energy distribution function f inj (E) and inject them in the CWFs at the contact discontinuity (see Fig. 1).
Particles then propagate through the CWFs according to the adopted diffusion model. The diffusion coefficients are chosen in the Bohm limit, i.e., defined by the particle energy and local magnetic field. Each particle propagates to a distance defined by its mean free path (mfp) given by the diffusion coefficient at its location, then it is scattered isotropically in the local plasma rest frame, then again propagates to a mfp distance -and so on, until it leaves the simulation. All the generated particles are propagated through the simulation area one by one. Multiple subsequent scatterings and crossings of the contact discontinuity in the zone of almost head-on winds collision allows particle to gain energy. Particles can be accelerated until either their mfp exceeds the size of accelerator allowing them to escape through the border with free escape boundary condition, or they lose too much energy due to adiabatic and radiative losses and therefore join the low-energy particle pool.
Relativistic protons radiate mainly due to interactions with stellar photons (photo-meson process, also dubbed as 'pγ') and SW protons (pp process). The radiative losses (including synchrotron and inverse Compton) are accounted for at each scattering by subtraction of the energy radiated by a particle in all the radiation processes during the preceding free propagation flight. We calculate the radiation loss rates using the approaches of Kelner et al. (2006); Kafexhiu et al. (2014) for the pp process and Kelner & Aharonian (2008) for the pγ process.
Detection of momenta of all particles leaving the simulation box with the Monte Carlo techniques allows us to get the spectral energy distribution of particles accelerated by the CWFs, averaged over the simulation volume. Then, the emission produced in pγ and pp processes is calculated using parametrizations of Kelner et al. (2006); Kelner & Aharonian (2008); Kafexhiu et al. (2014). We take into account the finite size of the simulation box, scaling the photon field ∝ r −2 with distance from the star r. To do so, we divide the box into a number of spatial bins, calculate the emission flux produced by each bin separately and finally summarize the results.
For a thermal distribution of stellar photons with T ≈ 30000 K, the pγ process works above a threshold proton Lorentz-factor Γ p > ∼ 10 7 (e.g. Dermer & Menon 2009). Thus, to enable the pγ process for the peak of CWFproduced hard energy distribution, the peak should lie at > ∼ 10 PeV.
The structure of the outer disks of Be stars was studied by Klement et al. (2017). They estimated the outer disk extensions to be up to ∼ 400 equatorial radii of the star and found that the disk height H scales as H ∝ r 3/2 . This means that the disk could extend to the distance of about 20 AU with the width well above 1 AU.
Maximum standoff distance R sd = Ė /4πρu 2 c ∼ 10 AU defining the CWFs acceleration region size is limited by the orbital separation distance. Here ρ and u are the SW mass density and velocity at the winds interaction region. Confinement of a particle with a gyroradius r g in the CWFs, r g < ∼ R sd , then requires the magnetic field > ∼ 1 G. For the surface dipolar field of the star ∼ 1 kG the toroidal SW field allows ∼ 2 G at r ∼ 2 × 10 14 cm. At r R * radial component of SW velocity is about its terminal velocity v ∞ ∼ 2GM * /R * ≈ 800 km s −1 , where G is the gravitational constant. This is much faster than the orbital velocity v K ∼ 30 km s −1 and the wind toroidal velocity ∼ few km s −1 . Thus, we assume the SW velocity u ≈ v ∞ and SW protons number density ∼ 3000 cm −3 corresponding to the SW disk density at stellar surface ∼ 10 −12 g cm −3 (Klement et al. 2017).
The timescale of the acceleration in the CWFs is much shorter than flare duration ∼ 10 7 s. According to Eq. (21) from Bykov et al. (2017) it takes ∼ 10 6 s to produce a hard energy distribution peaking at tens of PeV. The CWFs acceleration may provide the high enough value of γ-ray flux during the pulsar passage through the disk that also may take H/u ∼ 10 7 s. Figure 2 presents the simulated spectrum of accelerated protons. The red curve shows the injection spectrum, the purple curve -particle spectrum after the CWFs acceleration. The latter is much harder and peaks at a few tens of PeV, where most of particle energy is accumulated.

RESULTS AND DISCUSSION
In Figure 3 we show the simulated emission spectra produced by the modeled particle distribution. The pp gamma-ray flux (dotted curve) is well below the Carpet-2 measured flux, whereas pγ flux (dashed curve) produced by particles accelerated up to tens PeV allows us to explain the results of Carpet-2 flare measurements. At TeV energies the total flux produced in both hadronic processes does not exceed even the steady-state fluxes observed by VERITAS and MAGIC (Abeysekara et al. 2018). The Fermi LAT 1-10 GeV lightcurve of PSR J2032+4127 is stable along the whole orbit with the flux ∼ 10 −10 erg cm −2 s −1 (Chernyakova et al. 2020a). It is likely produced by the pulsar magnetospheric emission and the flux is well above that from hadronic interactions described by our model, see Fig. 3. The GeV-TeV spectra are successfully modelled within the leptonic scenario (see e.g. Chernyakova et al. 2020a;Takata et al. 2017). As it was mentioned above, severe radiation losses exclude the leptonic origin of PeV gamma-ray emission. Interestingly, the indications of leptons spectral break from s ∼ 2 to a harder s ∼ 1, similar to shown in Fig. 2 for protons spectrum, were obtained in the hard X-ray spectrum of LS 5039 detected with IBIS INTEGRAL (Falanga et al. 2021).
The γ-ray binaries may represent a possible population of cosmic ray sources with hard spectra. This can be important, e.g., to resolve the issues with energetics of high-energy radiation (see Murase & Fukugita 2019).
The effects of the absorption of GeV-TeV γ-rays in close γ-ray binaries containing luminous massive stars were discussed in the case of LS 5039, LSI +61°303 and Cyg X-3 (see, e.g., Bednarek 2006;Cerutti et al. 2011). However, in the source considered above, the attenuation of PeV photons of interest due to e ± pair production on the intensive stellar radiation is not essential.
The predicted neutrino flux is compared in Figure 3 with typical estimates of the flare and steady-state neutrino emission from IceCube data. Like in γ-rays, the flare neutrino flux is order of magnitude higher than the steady one. The overall contribution of this and similar sources to the IceCube astrophysical neutrino flux is limited by short duty cycles and by a small number of sources. We estimate that the time-averaged flux from all gamma-ray binaries can reach at most ∼ (10 − 13)% of the full-sky astrophysical neutrino flux, in agreement with constraints on the diffuse Galactic component (Albert et al. 2018); studies of point-source contributions (e.g., Liu & Kheirandish 2021) are less constraining yet.
In the PeV source model we considered a particular geometry of the mildly relativistic colliding flows of PW -massive Be star wind shown in Fig. 1. The same mechanism can accelerate the protons rapidly up to PeV energies in other γ-ray binaries, e.g. in case of the black hole jet colliding with the massive star wind also illustrated in Fig. 1. The realistic models of interaction of the fast outflow from a compact object with the wind of a massive star are still under debate. The models of LS 5039 with anisotropic PW (Bosch-Ramon 2021) and microquasars like Cyg X-1 and Cyg X-3 with relativistic jets from the accreting stellar mass black holes (see, e.g., Romero et al. 2017, and the references therein) represent a few of these. Apart from the microquasars and pulsar wind nebulae (e.g., Arons 2012; Amato 2020) the powerful outflows of the high enough kinetic luminosities (consistent with that given by Eq.1) in supernova remnants (Cristofari 2021), compact stellar clusters (Bykov et al. 2015) and superbubbles can accelerate protons above PeV and produce (sub-)PeV γ-rays and neutrino. The distinctive feature of the model of the γ-ray binary PSR J2032+4127 -MT91 213 discussed above is the combination of both the fast production of hard energy spectra of protons peaking in the PeV regime and the fast efficient photo-meson cooling in the vicinity of the B0Vp star which provides the highest PeV regime luminosity.

CONCLUSIONS
The model we discuss in the paper allows converting a sizeable fraction of the kinetic power of the outflows of the compact objects in close binary systems to PeV regime radiation due to the high efficiency of Fermi mechanism in colliding flows. Recent modeling by Pittard et al. (2021) of TeV regime particle acceleration in colliding wind binary with wind velocities ∼ few × 10 3 km s −1 and ∼ mG magnetic fields in the acceleration region demonstrated that ∼ 30% of the wind power was transferred to non-thermal particles. In our model with the short binary separation allowing for > ∼ G magnitude magnetic fields and mildly relativistic flows produced by shocked relativistic outflow of a compact object one can reach a similarly high efficiency but for proton acceleration to PeV energy regime. The rapid photo-meson cooling of the PeV protons converts a significant part of the available kinetic power to (sub-)PeV photons and neutrinos.
Our model generally predicts a transient character of the very bright PeV regime radiation. The main reason of that is the presence of the variations of the magnetic field, particle and seed photon densities along the orbit of the compact object due to the disks (for a Be companion star) or anisotropic stellar winds. This provides variability of the very-high-energy radiation along the compact object orbit. The timescale to cross the Be star disk could be a few months for the orbital parameters of PSR J2032+4127 which may explain the estimated duration of the γ-ray flare detected by Carpet-2.
An important factor is also the threshold character of the photo-meson radiation mechanism. Indeed, the energy of a seed photon in the rest frame of the accelerated proton must exceed ∼ 200 MeV to start the mechanism, which requires protons of energies above 10 PeV to interact with the intense optical radiation of the massive star.
A transient activation of the efficient and fast Fermi I type acceleration up to tens of PeV in the colliding flows during passage of the compact companion through the equatorial disk and rapid photo-meson cooling allows the interpretation of the very bright PeV photon flare detected by Carpet-2. Moreover, our mechanism unavoidably provides high neutrino flux, which could be detected by IceCube and other observatories. PeV neutrinos produced by the gamma-ray binaries in starforming galaxies may contribute to the observed population of the very high energy neutrinos.
20020. M.E.K. was supported by RFBR grant 20-32-90156. S.V.T. was supported by the Ministry of science and higher education of Russian Federation under the contract 075-15-2020-778. Some of the modeling was performed at the Joint Supercomputer Center JSCC RAS and at the "Tornado" subsystem of the St. Petersburg Polytechnic University supercomputing center.