Properties of a hypothetical cold pulsar wind in LS~5039

LS~5039 is a powerful gamma-ray binary that probably hosts a non-accreting pulsar. Despite the wealth of data available, the power source of the non-thermal emitter is still unknown. We use a dynamical-radiative numerical model and multiwavelength data to constrain the properties of a pulsar wind that may power the non-thermal emitter in LS~5039. We ran simulations of an ultrarelativistic (low-$B$) cold $e^\pm$-wind that Compton scatters stellar photons and that dynamically interacts with the stellar wind. The effects of energy losses on the unshocked $e^\pm$-wind dynamics, and the geometry of the two-wind contact discontinuity, are computed for different wind models. The predicted unshocked $e^\pm$-wind radiation at periastron, when expected to be highest, is compared to LS~5039 data. The minimum possible radiation from an isotropic cold $e^\pm$-wind overpredicts the X-ray to gamma-ray fluxes at periastron by a factor of $\sim 3$. In the anisotropic wind case X-ray and $\gtrsim 100$ MeV data are not violated by wind radiation if the wind axis is at $\lesssim 20-40^\circ$ from the line of sight (probability of $\lesssim 6-24$\%), depending on the anisotropic wind model, or if the wind Lorentz factor $\in 10^2-10^3$, in which case the wind power can be higher, but it requires $e^\pm$-multiplicities of $\sim 10^6$ and $10^9$ for a $10^{-2}$~s and 10~s pulsar period, respectively. The studied model predicts that a low-$B$ cold pulsar $e^\pm$-wind in LS~5039 should be strongly anisotropic, with either a wind Lorentz factor $\in 10^2-10^3$ and very high multiplicities or with a fine-tuned wind orientation. A low-$B$, cold baryon-dominated wind would be possible, but then the multiplicities should be rather low, while the baryon-to-$e^\pm$ energy transfer should be very efficient at wind termination. A strongly magnetized cold wind seems to be the most favorable (least constrained) option.


Introduction
LS 5039 is a strong galactic source of variable and periodic gamma rays (Paredes et al. 2000;Aharonian et al. 2005a;Abdo et al. 2009;Hadasch et al. 2012;Collmar & Zhang 2014;Chang et al. 2016). This source displays non-thermal radio emission with milliarcsecond to subarcsecond components (Paredes et al. 2000(Paredes et al. , 2002 and X-rays of likely non-thermal origin (see, e.g., Bosch-Ramon et al. 2007), presenting orbital changes in morphology (Ribó et al. 2008;Moldón et al. 2012a) and flux (see, e.g., Bosch-Ramon et al. 2005;Takahashi et al. 2009), respectively. The source is a compact binary with period P ≈ 3.91 d, semi-major axis a ≈ 2.1×10 12 cm, and eccentricity e ≈ 0.35, and hosts a O6.5V star and an undetermined compact object (CO) of a few solar masses (Casares et al. 2005). The source was initially thought to be a radio-loud X-ray binary (Marti et al. 1998), and a gamma-ray emitting microquasar later on (Paredes et al. 2000). The proposed microquasar nature and gamma-ray association led to the development of microquasar jet models for the non-thermal emitter (e.g., Paredes et al. 2000;Bosch-Ramon & Paredes 2004;Paredes et al. 2006;Bednarek 2006;Dermer & Böttcher 2006;Khangulyan et al. 2008), although the lack of accretion signatures in X-rays, among other source features, and the gamma-ray detection of PSR B1259−63/LS2883 (Aharonian et al. 2005b, a similar source hosting a non-accreting pulsar), led to the suggestion that LS 5039 hosted a young non-accreting pulsar (Martocchia et al. 2005;Dubus 2006b).
shorter period was not possible for Yoneda et al. (2020) as the data statistics only allowed the search of periods > 1 s.
Regardless of the specific pulsar nature, if a pulsar wind transports the energy from the CO to the non-thermal emitter, the wind can radiate much of its energy even before being shocked by the stellar wind. In particular, inverse Compton (IC) scattering off stellar photons by an ultrarelativistic cold pulsar e ± -wind (Bogovalov & Aharonian 2000) is very efficient in a compact binary (see, e.g., Ball & Kirk 2000;Sierpowska & Bednarek 2005;Khangulyan et al. 2007;Sierpowska-Bartosik & Torres 2007;Cerutti et al. 2008;Khangulyan et al. 2011;Hu et al. 2020), and in LS 5039 the expected IC fluxes may be higher than the observed ones (Cerutti et al. 2008). The total non-thermal luminosity of LS 5039 is consistently high all along the orbit, at a level L NT ≈ 10 36 erg s −1 (mostly released in the MeV-GeV range, Abdo et al. 2009;Collmar & Zhang 2014;at 2.1 kpc, Gaia Collaboration et al. 2018;Luri et al. 2018), which is hardly compatible with Doppler boosting (Db) models (see below), meaning that the pulsar wind luminosity L p must indeed be > L NT .
In this work we revisit the (weakly magnetized) cold pulsar e ± -wind model in the context of LS 5039 to constrain the hypothetical wind properties, using for the first time the most recent data in X-rays and soft and hard gamma rays. We model how such a wind transports the energy from the pulsar to its termination shock, which is produced by the interaction with the stellar wind. Unlike previous works, to compute the geometry of the stellar-pulsar wind interaction region we include the pulsar wind momentum-flux losses and orbital motion. As in Cerutti et al. (2008), we consider both isotropic and anisotropic wind models, and assume that the e ± -wind does not heat up due to IC braking. In Sects. 2 and 3 we present our wind model and its results, respectively, adopting parameter values that minimize the e ± -wind IC emission at periastron, the orbital phase when this emission should be the brightest (i.e., the most constraining choice). We conclude and discuss the results in Sect. 4.

Model
The cold pulsar e ± -wind propagation, radiation, and termination location were computed including the orbit effect on the geometry of the contact discontinuity (CD) between the pulsar and the stellar wind. A description of the model and the approach adopted for these calculations are presented in what follows.
To derive the most stringent constraints on the wind radiation towards the observer, we focused on periastron, when this radiation is expected to be highest (see Cerutti et al. 2008; LS 5039 periastron orbital phase, φ = 0, is close to superior conjunction of the CO, φ ≈ 0.058, see Casares et al. 2005). At the same time, we were interested in a parameter choice predicting the minimum possible IC emission for that orbital phase. Thus, regarding the inclination, a value of i ≈ 45 • was chosen, which also avoids black hole masses for the CO or stellar eclipse (Casares et al. 2005). Moreover, we fixed L p = 2 × 10 36 erg s −1 , which is L p ≈ 2L NT ; in other words, the non-thermal emitter must emit ≈ 50% of the available power. As L NT ∼ 10 36 erg s −1 all through the orbit, Db cannot be behind the high (but then apparent) L NT unless the emitting flow is always roughly pointing to the observer, which is implausible in a colliding wind system (although Db can still significantly affect the lightcurve; see, e.g., Zabalza et al. 2013;Dubus et al. 2015;Molina & Bosch-Ramon 2020). For comparison, the luminosities injected in the emitting e ± of the shocked wind in Zabalza et al. (2013), Dubus et al. (2015), and Molina & Bosch-Ramon (2020), which must be < L p , were L rel−e = 4×10 35 , 10 35 , and 6×10 35 erg s −1 , respectively. All three works accounted for Db in the shocked-wind emitter, and the last two severely underpredicted the flux in MeV (not considered by Zabalza et al. 2013).
We computed for all directions the energy evolution of the wind e ± , which moved from the CO along straight trajectories as they IC scattered stellar photons. We used an iterative scheme to follow the e ± energy, with IC energy losses and emission determined by IC interactions in the anisotropic blackbody stellar field (Khangulyan et al. 2014). Descriptions of how to calculate the e ± propagation, energy loss, and emission can be found in Cerutti et al. (2008) and Khangulyan et al. (2011), among others. For simplicity, a point-like star was considered when computing the IC emission in different directions, which is reasonable as i is such that gamma rays are not eclipsed. When computing the anisotropic IC losses for all directions the interaction angle was set to be at least the angle subtended by the star and at most its supplementary, as seen from the emitting point, accounting in this way for the actual star finite size.
The wind magnetic field (B) was taken to be irrelevant both dynamically and, as a cold e ± wind does not produce synchrotron radiation, radiatively. The explored initial Lorentz factors of the wind were Γ p,0 ∈ 10−10 5 . Gamma-ray absorption and reprocessing were not included in the calculations because Γ p,0 10 5 is needed for efficient pair creation and cascading in the stellar photon field. Not accounting for such a high Γ p,0 is reasonable as in LS 5039 the reprocessed gamma-ray luminosity of an e ± -wind should be still a significant fraction of L p , peaking around multi-GeV energies (see Sierpowska & Bednarek 2005, for a similar system), which is at odds with the rather low flux of LS 5039 at ∼ 10 GeV (see the overall gamma-ray spectrum in, e.g., Hadasch et al. 2012;Collmar & Zhang 2014). Not accounting for gamma-ray reprocessing simplifies the calculations greatly as otherwise one would need a detailed knowledge of the background medium (for these process complexities in LS 5049, see, e.g., Aharonian et al. 2006 The IC braking effect on Γ p was computed for all the possible trajectories up to three times the orbital separation (i.e., 3 × R orb ), and the shape of the two-wind CD was calculated including the orbit effect. Five approximations were adopted for the initial e ±energy and for the energy and momentum flux angular dependence: (case A) an isotropic e ± -wind; (cases B and B') a phenomenological axisymmetric e ± -wind with constant Γ p,0 and for B and for B', being θ p the angle with the pulsar wind symmetry axis; and (cases C and C') a more physical axisymmetric wind with dL p /dΩ ∝ Γ p,0 (θ p ) = (10 +Γ p,0 sin 2 θ p ) for C and dL p /dΩ ∝ Γ p,0 (θ p ) = (10 +Γ p,0 sin 4 θ p ) for C' (see, e.g., Cerutti et al. 2008;Tchekhovskoy et al. 2016, andKhangoulyan 2002 in the context of Crab).
The explored values forΓ p,0 were ∈ 10 − 10 5 (as they were for Γ p,0 ). The line-of-sight angle of the e ± -wind axis in cases B and Article number, page 2 of 6 V. Bosch-Ramon: On the hypothetical pulsar wind of LS 5039 C was set to ψ = 20 • , and in cases B' and C' to ψ = 40 • , not to violate the observed fluxes. Several e ± -wind orientations were tested fulfilling ψ ≈ 20 • and 40 • , and to be conservative we chose those with the lowest IC fluxes. The probability of having a smaller ψ-value was P 6% for cases B and C, and P 24% for cases B' and C'. The stellar wind was assumed to be spherically symmetric, although a prescription for its momentum flux ∝ sin 2 θ w was also tested, with the plane defined by θ w = π/2 passing through the pulsar center (the most conservative case); both prescriptions yielded very similar results.
The stellar wind was assumed to be cold, with momentum flux oḟ and velocity of in the non-inertial frame rotating with the orbit angular velocity (Ω), where θ is the angle to the orbit normal, R * the stellar radius, and d the stellar distance. To normalizeṗ w , the maximum possible mass-loss rate for this wind,Ṁ = 7×10 −7 M yr −1 , was adopted (Casares et al. 2005), minimizing the distance covered by the unshocked e ± -wind, and thus its IC emission. Nevertheless, theṀ-dependence of the e ± -wind IC luminosity is weak, as the pulsar wind facing the observer at periastron typically loses most of its energy (see below). The stellar wind velocity at infinity v ∞ was fixed to 2.44 × 10 8 cm s −1 , and β to 1, while the minimum radial velocity was set to the escape velocity, ≈ 10 8 cm s −1 (McSwain et al. 2004).
The e ± -wind, little affected by Coriolis forces due to its high speed, was assumed with zero pressure, radial from the pulsar, and with momentum flux oḟ and velocity of The e ± -wind anisotropy was introduced via ρ p (B and B') and Γ p,0 (C and C'). Losses affected Γ p . The CD was computed using an approach based on the thin shell axisymmetric approximation (see, e.g., Antokhin et al. 2004), but modified to account for the lack of axisymmetry due to the orbit-related Coriolis force and the wind anisotropy.
First, the CD stagnation point (SP), at whicḣ is found. We looked for the SP assuming that it was on the orbital plane, in a direction from the pulsar tilted from the pulsar-star direction clockwise by an angle where ω orb is the orbit angular velocity and η = (L p /Ṁ|v w |c) is the momentum rate ratio of the pulsar and the stellar wind. This approximation works well for small τ, so we adopted this method in our calculations as a first-order approximation. Once found, we took the SP position as the initial one and computed the paths shaping the CD in all directions through an iterative process. The iterative steps in each path provide segments that characterize the CD surface step by step. To find the directions of the new segments, we used the relations |ṗ w | sin 2 α w = |ṗ p | sin 2 α p (11) anḋ p s =ṗ w sin α w +ṗ p sin α p (12) at each step, where α w and α p are the complementary angles between the stellar and the pulsar wind directions and the CD normal, respectively, andṗ s is the momentum flux remaining after wind collision, at the segment location on the CD. The thin shell approximation allows the cancellation of the opposing momentum components. Once a direction on the CD surface from the SP is chosen, the next point is defined by the local direction ofṗ s and a sufficiently small distance in that direction. Even if α w and α p are not known, the direction ofṗ s can be computed by eliminating one of the α parameters using Eq. (11), and renormalizinġ p s to eliminate the other one from Eq. (12).
To account for the shocked region of the e ± -wind between the termination shock and the CD, the e ± -wind was terminated at ≈ 2/3 of the distance from the pulsar to the CD (see Bogovalov et al. 2008 for cases with η 0.01, although if B were strong it would largely modify the overall interaction structure, see Bogovalov et al. 2019). In the presence of orbital motion, analytical and numerical calculations show that a strong shock should terminate the pulsar wind due to Coriolis forces in the directions away from the star (see, e.g., Bosch-Ramon & Barkov 2011;Bosch-Ramon et al. 2012Huber et al. 2020). Therefore, the unshocked pulsar wind, and the CD itself, were simply assumed to stop in these directions at the distance where such a Coriolis shock is expected to form. Since we focus on periastron, our results are not affected by this assumption as the relevant line of sight does cross the computed CD at that phase.
The paths characterizing the CD surface for case A are shown in the top panel of Fig. 1 for Γ p,0 = 10 and 10 5 , with and without orbital motion (the CD shapes for B and C, and B' and C' were very similar, although a detailed discussion is beyond the scope of this work). The CD shrinks by a factor of ∼ 2 due to stronger IC losses for Γ p,0 = 10 5 , as the cooling time is ∝ Γ −1 p,0 in the explored range. Orbital motion clearly affects the CD geometry on the scales of the orbital separation for the parameter values adopted, although the computed e ± -wind IC fluxes without orbital motion (see the dashed line in the top panel of Fig. 1) were almost the same. The bottom panel of Figure 1 shows two maps of the fraction of the energy remaining in the e ± after IC losses for case A, with Γ p,0 = 10 (right) and 10 5 (left), together with the corresponding orbital-plane CD shapes. As expected, for Γ p,0 = 10 e ± lose ∼ 1% of their energy or less when reaching the CD, whereas for Γ p,0 = 10 5 they lose ∼ 80% towards the star, and ∼ 20% away from it.
The e ± -wind spectral energy distributions are shown in Fig. 2 for Γ p,0 ∈ 10 − 10 5 , and for cases A, B, and C. The spectral Orbital plane projection of the computed paths on the CD surface, at periastron, for case A and Γ p,0 = 10 (green; more extended) and 10 5 (blue; more compact), including orbital motion. The red dashed line shows the profile of the structure for Γ p,0 = 10 5 when orbital motion is not included. The projected line of sight points along −ŷ and the star is at (0, 0). Bottom panel: Color maps of the orbital plane of the remaining energy fraction, after IC losses, for particles propagating from the pulsar in all directions, for Γ p,0 = 10 (left) and 10 5 (right). The white solid lines are the profiles of the projected CD on the orbital plane at periastron for case A including orbital motion. The pulsar is at (0, 0). Both panels: x-and y-axis units are the orbital distance. energy distributions for cases B' and C' (not shown) are very similar to those of B and C, as expected given that the ψ-value choice is determined by the same observational constraints. An isotropic e ± -wind (A) overpredicts the observed fluxes by at least a factor of ≈ 3. The anisotropic e ± -wind of cases B and C, with ψ 20 • between the wind symmetry axis and the line of sight, does not overpredict the observed X-ray or 100 MeV fluxes if Γ p,0 10 2 − 10 3 . The IC flux of cases B and C is also lower than the observed ∼ 0.1 − 30 MeV fluxes by a factor of ≈ 3 − 5 for Γ p,0 ∈ 10 2 − 10 3 , which allows the relaxation of the constraint on ψ from 20 • to 45 • (P 30%) for L p 2 × 10 36 erg s −1 . In the anisotropic e ± -wind model of cases B' and C', the observed X-ray and 100 MeV fluxes constrain ψ to 40 • (P 24%) if Γ p,0 10 2 −10 3 . On the other hand, the observed ∼ 0.1−30 MeV fluxes do not constrain ψ at all in cases B' and C' for Γ p,0 ∈ 10 2 − 10 3 if L p 2 × 10 36 erg s −1 . The X-ray data, having high statistics, constrains a narrow spectral component to a level of a ∼ 10% of the observed fluxes (see, e.g., Bosch-Ramon et al. 2007). Thus, Cases B and B' are also constrained by X-ray data even though the predicted flux is ∼ 3 times lower than observed. Fig. 2. Computed spectral energy distributions of the IC cold pulsar wind emission, at periastron, for 11 different Γ p,0 values from Γ p,0 = 10 to 10 5 , for cases A (dotted lines), B (dashed lines), and C (solid lines) (cases B' and C' are very similar; see main text). A schematic spectral energy distribution based on the observed X-ray to gamma-ray emission around superior conjunction (adapted from fig. 9 in Collmar & Zhang 2014; see references therein), which is close to periastron, is also shown.

Weakly magnetized cold wind model
If LS 5039 has a weakly magnetized cold pulsar e ± -wind carrying the energy to the non-thermal emitter, this wind cannot be isotropic as this is inconsistent with observations, even for ratios of the total broadband non-thermal radiation luminosity to the pulsar wind power of L NT /L p ∼ 1.
An anisotropic wind with angular dependence ∝ sin 2 θ p does not violate the observational constraints, but requires either a very high L NT /L p ∼ 1, even for the most favorable Γ p,0 values, or an unlikely orientation. For cases B and C and Γ p,0 10 2 − 10 3 , even for an unrealistic L NT /L p ≈ 0.5 a rather improbable ψ 20 • (P 6%) is required for the predictions to be consistent with X-ray and gamma-ray data. For the same cases but where Γ p,0 ∈ 10 2 − 10 3 , ψ is less constrained, 45 • (P 30%), if L NT /L p ≈ 0.5 is kept. Smaller ψ-values relax the energetic constraints when Γ p,0 ∈ 10 2 − 10 3 , but reduce the associated probability: for L NT /L p ≈ 0.1, ψ 20 • and thus P 6%. We note that L NT /L p ≈ 0.1 is still quite high, typically values of L rel−e /L p ∼ 0.1 are adopted in the modeling of the whole LS 5039 non-thermal emission, resulting in an even lower L NT /L p value. We recall that invoking Db can reduce the energetic requirements in some orbital phases, but not everywhere in the orbit.
If the wind is even more anisotropic, with angular dependence ∝ sin 4 θ p , the situation becomes less dramatic, although quite high L NT /L p values are still needed in general. For cases B and C and Γ p,0 10 2 − 10 3 , setting L NT /L p ≈ 0.5 implies that ψ 40 • (P 24%) does not violate the observed fluxes, whereas for Γ p,0 ∈ 10 2 − 10 3 the angle ψ is not constrained at all. Adopting ψ-values significantly smaller than ≈ 40 • , L NT /L p becomes virtually unconstrained for Γ p,0 ∈ 10 − 10 5 in B', and for Γ p,0 ∈ 10 2 − 10 5 in case C' (in C' taking Γ p,0 10 2 still overpredicts the X-ray fluxes). Therefore, in cases B' and C' the data are less restrictive than B and C of the e ± -wind properties, but some fine-tuning in ψ or Γ p,0 is still needed if a high L NT /L p is to be avoided.
Despite all this discussion on wind anisotropy, it is worth noting that fine-tuning the wind orientation is less effective in relaxing the energetic requirements on L p if the pulsar precesses (Stairs et al. 2000). In addition, the wind anisotropy is likely to be more complex than just described (see, e.g., Philippov & Spitkovsky 2018), which can smear its effects and thus make wind orientation fine-tuning less effective in relaxing the observational constraints.
To explore the e ± -multiplicity κ in the studied wind model, one can take for instance Γ p,0 = 3 × 10 2 and L p ∼ 10 37 erg s −1 , values allowed by observations for ψ 20 • in B and C, and in a broader range of ψ-values in B' and C'. With this choice of parameters, κ should be ∼ 10 6 (∼ 10 9 ) for a 10 km radius pulsar with 10 12 G (10 15 G) surface B and 10 −2 s (10 s) period.
A weakly magnetized cold baryonic wind is an alternative to a pure e ± -wind model. Such a wind would require that baryons reach Γ p,0 5 × 10 4 ( 5 × 10 7 -10 s-) for L p 2 × 10 36 erg s −1 , while κ 10 ( 10 4 -10 s-), so as not to violate the 10 GeV ( 1 TeV -10 s-) fluxes for an isotropic wind, and 10 2 ( 10 5 -10 s-) for an anisotropic wind with ψ ≈ 20 • (B and C) or ≈ 40 • (B' and C'). We note that this scenario requires very efficient proton-to-e ± transfer at wind termination, as multiplicity limits become tighter for higher proton wind luminosity budgets. We consider here protons solely for energy transport as a nonthermal proton emitter in LS 5039 is likely to be radiatively inefficient due to not enough ambient photon energy, and radiation and wind densities (Bosch-Ramon & Khangulyan 2009).

Magnetized winds
A strongly magnetized flow produced by a young magnetar in LS 5039 was proposed by Yoneda et al. (2020) prompted by the evidence of a ≈ 9 s period in X-ray data. As indicated in Sect. 1, the young magnetar hypothesis put forward in that work is at odds with the lack of evidence for a SNR, and the physics of such a flow would be more uncertain than that of more standard pulsar wind models. Interestingly, it is also possibile to envision an exotic explanation for the ≈ 9 s period unrelated to pulsar rotation. LS 5039 might instead host not one CO, but an extremely close CO binary of semi-major axis ≈ 10 9 cm. Such a binary, given the mass function of the system (Casares et al. 2005), would most likely be formed by two neutron stars. The neutron stars might have rotation periods much shorter than ≈ 9 s, for instance in the millisecond range. Now, whether such a configuration is plausible from the point of view of stellar evolution is to be studied elsewhere (see, e.g., Abbott et al. 2020, in the context of gravitational wave detections of related systems).
Nevertheless, for the single-CO scenario a pulsar wind with a dynamically significant magnetic field would be in agreement with recent research on pulsar winds, which proposes average magnetization parameters σ ∼ 1 instead of ∼ 10 −2 −10 −3 (Amato 2020), where withṁ being the wind mass rate. Unfortunately, the complexity of a magnetized wind (as exemplified for instance in the simulations by Philippov & Spitkovsky 2018) does not allow the use of a simple, but still physically consistent, prescription for the wind like the one employed here for the cold, matter-dominated case. However, we can already show that a strongly magnetized wind may be, given the observational constraints, the most favored scenario in LS 5039: As in the case of a cold baryonic wind, a magnetized wind would not produce prominent narrow X-ray and gamma-ray spectral features if the wind e ± did not get energized by Bdissipation until wind termination (and e ± -energy redistribution). However, for a strongly magnetized wind, the constraints on κ and Γ p,0 are looser than in the baryonic wind case, as now Γ p,0 is also not constrained. For Γ p,0 10 4 , we obtain κΓ p,0 10 7 ( 10 10 -10 s-) for case A, implying σ 10. For Γ p,0 10 4 , σ must be higher due to tighter observational constraints above the GeV range. For cases B, C, B', and C' the observational constraints on a weakly magnetized wind are looser for specific orientations, and thus less demanding on sigma, but high sigma values render wind orientation fine-tuning (which has its caveats, as already mentioned) unnecessary. Therefore, we conclude that a high-sigma wind is less constrained and thus seems more favorable than the other two scenarios.

Hot winds
It is possible that, relatively early in its propagation, the pulsar wind turns a significant amount of energy (magnetic or kinetic) into non-thermal energy, meaning that the non-thermal emitter would start much closer to the pulsar than the wind termination shock (see, e.g., Sierpowska-Bartosik & Torres 2007;Pétri & Dubus 2011;Derishev & Aharonian 2012, for the presence of non-thermal particles in the wind of LS 5039 due, e.g., to magnetic dissipation, IC cascades). This scenario is perhaps the hardest to model, and may not be consistent with evidence of fast e ± -wind acceleration in the Crab pulsar (see Aharonian et al. 2012), although in a high-mass binary system the unshocked pulsar wind has a very different environment than in an isolated pulsar. On the other hand, it is worth noting that observations of LS 5039 strongly suggest that significant particle acceleration is still required far from the CO as the X-ray and particularly the very high-energy emission seem to originate in the binary outskirts (see Sect. 1). This peripheric accelerator must be very efficient as acceleration rates close to the electrodynamical limit have been inferred from the very high-energy data (see Khangulyan et al. 2008). Nevertheless, it cannot be discounted that the bulk of the ∼ 0.1 − 30 MeV emission may still be produced in a hot wind, relatively close to the pulsar (as proposed in Derishev & Aharonian 2012, effectively similar to the cold e ∓ -wind case with Γ p,0 ∈ 10 2 − 10 3 discussed here) 1 .

Similar sources
We finish by noting that, in addition to LS 5039, four other highmass gamma-ray binaries, LS I +61 303, 1FGL J1018.6−5856, LMC P3, and 4FGL J1405−6119, are approximately as compact and at least as powerful as LS 5039, with a relatively similar phenomenology and unknown CO (see Sect. 1 in Molina & Bosch-Ramon 2020, and references therein). Although all these sources deserve specific studies of their own, an approach such as the one presented here can be helpful to constrain the properties of a hypothetical pulsar wind powering their non-thermal emitter.