The major role of eccentricity in the evolution of colliding pulsar-stellar winds

Binary systems that host a massive star and a non-accreting pulsar can be powerful non-thermal emitters. The relativistic pulsar wind and the non-relativistic stellar outflows interact along the orbit, producing ultrarelativistic particles that radiate from radio to gamma rays. To properly characterize the physics of these sources, and better understand their emission and impact on the environment, careful modelling of the outflow interactions, spanning a broad range of spatial and temporal scales, is needed. Full 3-dimensional approaches are very computationally expensive, but simpler approximate approaches, while still realistic at the semi-quantitative level, are available. We present here the results of calculations done with a quasi 3-dimensional scheme to compute the evolution of the interacting flows in a region spanning in size up to a thousand times the size of the binary. In particular, we analyze for the first time the role of different eccentricities in the large scale evolution of the shocked flows. We find that the higher the eccentricity, the closer the flows behave like a one-side outflow, which becomes rather collimated for eccentricity values $\gtrsim 0.75$. The simulations also unveil that the pulsar and the stellar winds become fully mixed within the grid for low eccentricity systems, presenting a more stochastic behavior at large scales than in the highly eccentric systems.


Introduction
Binary systems hosting a massive star and a non-accreting pulsar, or pulsar highmass binaries (PHMB), can be powerful sources of gamma rays. The objects of this kind capable of gamma-ray emission pertain to the wider class of gamma-ray binaries, in which most of the non-stellar radiation is released in the gamma-ray energy range (see, e.g., Dubus, 2013;Paredes and Bordas, 2019a,b, for these and related sources). The radiation is produced through the interaction of a relativistic pulsar wind and the outflows ejected by the star: Massive stars produce strong non-relativistic winds of supersonic nature and, in cases of very fast stellar rotation, quasi-Keplerian equatorial disks (or decretion disks) that flow outwards at subsonic speeds. These outflows interact with the pulsar wind and later on interstellar matter in a process in which ultrarelativistic particles are accelerated and produce emission from radio to gamma rays (see, e.g., Tavani and Arons, 1997;Sierpowska and Bednarek, 2005;Dubus, 2006;Neronov and Chernyakova, 2007;Khangulyan et al., 2007;Kong et al., 2012;Zabalza et al., 2013;Dubus et al., 2015;Molina and Bosch-Ramon, 2020;Huber et al., 2021;Lyutikov et al., 2020;Khangulyan et al., 2021).
The pulsar wind-stellar outflow interactions are complex, and different regions that can influence each other are relevant when trying to understand the evolution of the shocked flows. Even the associated radiation and its reprocessing can feedback on the flow dynamics, making the whole physical system highly non-linear. At small scales, there is the region right between the star and the pulsar in which flows are stopped and shocked by colliding against each other. After the collision, the flows become subsonic and hot, and start moving symmetrically sideways while pressure gradients lead to their reacceleration, getting supersonic again. This picture is roughly similar in both sides of the contact discontinuity separating the stellar and the pulsar shocked flows. Later on, the respective evolution of the flows largely differs due, for instance, to very different wind momentum rates and initial velocities, plus orbital effects.
Due to the large momentum rate of stellar winds, the stellar wind confines the pulsar wind, which bends over the pulsar. Shocked winds form an approximately axisymmetric curved structure that becomes conical further away from the binary. On those larger scales, if orbital motion were neglected, the shocked winds would move ballistically and form a conical shell made of shocked pulsar wind, surrounded by another shell of shocked stellar wind. The half-opening angle of the conical contact discontinuity would converge to a value that can be derived from the pulsar-to-stellar wind momentum rate ratio (Bogovalov et al., 2008): where L sd is pulsar spin-down luminosity,Ṁ w and v w are the stellar massloss rate and wind speed, respectively. Following Eichler and Usov (1993) and Bogovalov et al. (2008), the contact discontinuity of the cone-like structure has an approximate opening angle of However, orbital motion is to be included in the colliding-wind picture, which makes a Coriolis force appear, a force that affects differently the pulsar and the stellar outflows due to the large relative velocity and density contrast. This differential Coriolis effect makes the stellar wind push on the shocked pulsar wind against the orbital rotation sense, creating a strong deflection of the interaction structure in that direction, and triggering a strong lateral shock in the shocked pulsar wind. As the two shocked flows have very different densities and velocities, they are prone to the occurrence of strong instabilities, such as Rayleigh-Taylor, Kelvin-Helmholtz, and Richtmyer-Meshkov, in the contact discontinuity (Bosch-Ramon et al., 2015). Thus, as the flows move, the outflow contact surface gets partially disrupted, stellar wind mixes with the shocked pulsar wind, and the latter develops strong turbulence and decelerates. The result is that the shocked flow structure shape becomes a one-arm spiral that fills much of the volume and is expected to disrupt after a few turns.
Between the apex of the interaction structure, located at the two-shocked flow stagnation point (Bogovalov et al., 2012), and the starting point of the Coriolis shock, on the leading edge of the interaction structure, the shocked pulsar wind gets compressed and thus heated by the Coriolis force-related lateral pressure of the stellar wind. This can weaken the mentioned shocked flow reacceleration caused by pressure gradients. On the other hand, in the trailing edge of that interaction structure, the shocked pulsar wind quickly expands and accelerates through rarefaction waves.
The presence of decretion disks can significantly alter the geometry of the interaction structure, which must develop now embedded in a much more complex circumstellar environment. Nevertheless, the disk is rather massive and marginally bounded to the star, so part of the material may not even escape the binary. In addition, the accumulated disk mass can be in fact just comparable to that of the stellar wind. Thus, on large scales, the shocked flow dynamics is likely dominated by the pulsar and the stellar wind, the disk and radiation processes can be important for flow dynamics on small and middle scales. Assuming then that the disk is mostly relevant closer to the binary, and neglecting the role of the magnetic field, the system eccentricity may turn out to be as important as η to describe the evolution of the shocked flows on large scales. In particular, Bosch-Ramon et al. (2017); Barkov and Bosch-Ramon (2018) show that for very high eccentricities the shocked pulsar wind becomes strongly focused along the periastron-apastron direction, as it gets deflected by the stellar wind in that direction for most of the orbit. To date, however, an exploration of how different eccentricity values affect the large-scale shocked flow structure is missing, mainly, at which eccentricities one-sided outflows form.
In this work, we perform a numerical study of how the shocked flows from PHMB evolve, and propagate, up to large distances from the binary for different eccentricities. Our major goal is to find for which orbit eccentricity the mentioned one-sided outflow forms. The study is carried out using the quasi 3-dimensional (3D) calculation scheme developed by Barkov and Bosch-Ramon (2016), used to study PSR B1259−63, and HESS J0632+057 (Bosch-Ramon et al., 2017;Barkov and Bosch-Ramon, 2018). The advantage of this method, which employs spherical coordinates, is that it focuses on the orbital plane, but sacrifices resolution for zenital angles far from that plane. Furthermore, the shocked flow geometry just outside the binary turns out to be amenable to be simplified such that the colliding-wind apex region does not need to be modelled, which allows a computationally much cheaper resolution. All this largely reduces the cost of the simulations, allowing one to probe a large region surrounding the binary. The accuracy of the method is appropriate at a semiquantitative level, as shown by comparison with results obtained using full 3D calculations encompassing overlapping regions (Bosch-Ramon et al., 2015).
The article is organized as follows: In Sect. 2, the simulations are described, and in Sect. 3, their results presented. Then, in Sect 4, the simulations results are summarized and discussed.

Numerical model
Quasi-3D simulations of PHMB wind-wind collisions with different orbit eccentricities were performed using the PLUTO code 1 (Mignone et al., 2007). PLUTO is a modular Godunov-type code entirely written in C and intended mainly for astrophysical applications and high Mach number flows in multiple spatial dimensions. Spatial parabolic interpolation, a 3rd order Runge-Kutta approximation in time, and an HLLC Riemann solver were used (Li, 2005). The simulations were performed on the CFCA XC30 cluster of the National Astronomical Observatory of Japan (NAOJ). To reduce computing costs, the flow was approximated by a simple equation of state enough for our purposes: that of an ideal relativistic gas with adiabatic index 4/3. We adopted spherical coordinates (R, θ, φ), with 768 cells in both the radial and the azimuthal directions. To reduce the computation costs, we took only 3 cells in the zenital direction. The domain size was taken to be R ∈ [1, 500]R min , θ ∈ [π/4, 3π/4] and φ ∈ [0, 2π]. We set R min = 2(1 − e 2 )a, with a being the semi-major axis of the orbit, and e its eccentricity; thus, the scales captured by the simulations are larger than a few times the orbital separation distance at periastron. The computational grid was made logarithmic in the radial direction, that is, cells grow with a constant aspect ratio. Our treatment of the θ-direction allows a reasonably realistic characterization of the orbital plane physics on scales beyond the pulsar, where the interaction structure expansion becomes approximately linear with distance, although a more quantitative account would require a complete 3D treatment (see Barkov and Bosch-Ramon, 2016;Bosch-Ramon et al., 2017;Barkov and Bosch-Ramon, 2018).
In the quasi-3D calculation scheme adopted here the injected pulsar wind has a half-opening angle φ c , which depends on the momentum rate relation, which is set to η = 0.1, an intermediate value for this parameter, to restrict the degrees of freedom of the problem. This η-value corresponds to φ c = 0.87 radians, so the computational domain was divided in two non-equal parts. The first one, in which φ ∈ (φ c , 2π − φ c ), was filled by a radial stellar wind with velocity value v w = 2400 km s −1 . The second part, in which φ ∈ [−φ c , φ c ], was filled by a radial pulsar wind with Lorentz factor Γ = 1.9. Despite this Lorentz factor being just moderately relativistic, the calculations are enough relativistic for our purposes because internal energy density already plays a significant inertial role (see Bosch-Ramon et al., 2012. The winds were assumed to be highly supersonic at injection, with Mach numbers M w = 6.9 and M j = γ p v p /γ s,p c s,p = 14 for the stellar and the pulsar wind, respectively, where "s" refers to the sound speed. For simplicity, any wind azimuthal velocity component, coming for instance from object rotation and angular momentum conservation, was neglected in our calculations as its value would be well below v w and c. Also, as discussed in Sect.1, we neglected the role of a decretion disk, although more quantitative studies should include it. The initial pulsar position is at the left of the computational domain, which means that the simulated pulsar wind cone is also initially directed to the left, which corresponds to the periastron-apastron direction. The adopted orbital period is T orb = 16.6 days, and the stellar masses are M M S = 31 M and M P SR = 1.44 M for the star and the pulsar, respectively, so from Kepler's third law the corresponding orbital semi-major axis is a = 6 × 10 12 cm. During the simulation, the φ-intervals within which the pulsar and the stellar winds are injected rotate along the orbit with the pace and sense of the corresponding orbital velocity. The studied cases have orbit eccentricities e = 0, 0.25, 0.5 and 0.75. Simulations with even larger e-values can be found in Barkov and Bosch-Ramon (2016); Bosch-Ramon et al. (2017); Barkov and Bosch-Ramon (2018). These e-values are characteristic of the different known high-mass gamma-ray binaries, although we note that a pulsar has been confirmed to be present only in PSR B1259−63 and PSR J2032+4127 (e.g. Aharonian et al., 2005;Lyne et al., 2015). On the other hand, the simulated orbital period is significantly shorter than those of HESS J0632+057, PSR B1259−63 and PSR J2032+4127. This was done as previous studies already explored cases with long periods (Barkov and Bosch-Ramon, 2016;Bosch-Ramon et al., 2017;Barkov and Bosch-Ramon, 2018). It is worth noting that the asymmetry of the interaction structure on large scales should not depend significantly on T orb , because what characterizes this asymmetry is the relative change with φ of the shocked pulsar wind energy rate along the orbit, which is independent of T orb . We note that a potential v φ component of the stellar wind would affect the angular distribution of the flow, but since v φ ∝ 1/R this component will become negligible at R R min .

Results
Maps of the density distributions and the velocity vector fields in the orbital plane are presented in Fig. 1 for the different eccentricities studied. Regardless of the eccentricity, the spiral structure starts to disrupt after one orbital turn, and after 2-3 orbital turns, spiral structure disruption leaves a more or less uniform medium with randomly located density and velocity irregularities. In particular, the velocity field shows spiral motion in the first turns and, farther away from the system, it mostly shows a radial outflow. Asymmetry of the interaction structure on the orbital plane starts to become significant for e 0.5. In the case of e = 0.75, the effect becomes extreme, with the density (velocity) in the periastronapastron direction being much smaller (larger) than in the other directions. For e = 0.5, this effect is also present, but not so strong. A fluid property used to track fluid motion -the tracer behavior, shown in Fig. 2, is similar to that of the density, and one can see there as well that the spiral structure disappears after 2-3 orbital turns regardless of eccentricity. The higher the e-value, the more prominent the pulsar wind becomes in the periastron-apastron direction. The pressure spatial distribution is presented in Fig. 3, showing the same trends as density and tracer. In addition, pressure shows a smooth drop after spiral structure disruption. As for density and tracer, the main difference between cases with low and high eccentricities is that pressure falls anisotropically, its decrease being slightly shallower in the periastron-apastron direction when e is large enough. We averaged, weighting in mass, the radial velocity and the Mach number in three φ-ranges, or sectors: s1, with φ s1 ∈ [150 o ,   averaged over a certain R-or φ-range, was calculated as and respectively, where ρ is density, q the value being averaged, and sX the corresponding sector. The radial distributions of v r averaged over s1, s2, and s3, at simulation time T = 1.11 × 10 5 are presented in Fig. 4. We note that time in this work is given in simulation units, which are (a/c) = 200 s and are implicit. The figure includes the case for s1 at T = 1.09 × 10 5 as well to illustrate short-term variability on top of the longer term behavior. Independently of the eccentricity, there are strong spatial variations in the radial velocity up to R ∼ 300 a that are related to the spiral interaction structure. Beyond R ∼ 300 a, there is a gradual acceleration of the flow, which is at that point already made of a mixture of pulsar and stellar winds. This acceleration is equally prominent for all sectors when e = 0, although the higher the eccentricity, the stronger the acceleration becomes for s1, around the periastron-apastron direction, and weaker for s2 and s3. This is also seen in Fig. 5, which displays the dependence of the sectoraveraged v r with e at R = 400 a. In general, the acceleration slows down at R ∼ 400 − 600 a depending on e, which is expected as the flow becomes highly supersonic. In Fig. 6, the radial distribution of averaged Mach number is also shown. At smaller radii, R 300 a, large variations of this quantity are related to the spiral structure, whereas at larger radii the flow averaged Mach number approaches ∼ 6 (being still strongly variable for high e). Despite the shallower drop of pressure in s1 shown above (see Fig. 5), v r grows faster in that direction because more energy is invested in the flow motion whereas the average density is lower.
To illustrate the temporal evolution of the system, color maps of sectoraveraged radial velocity and Mach number in the R (y-axis) versus T (x-axis) plane are presented in Figs. 7 and 8, respectively. These maps show two different colour regions that indicate the transition at R < 300 a from a spiral structure to a more homogeneous, mixed, outflow at larger radii. This effect is very prominent in all directions for low eccentricities, both in averaged v r and Mach number, whereas for high eccentricities the same effect is much more prominent in s1 (around the periastron-apastron direction) than in s2 and s3 for v r , whereas for Mach number the opposite happens, although less dramatically. The slower increase in the s1-averaged Mach number for high e-values is related to the presence of shocks that reheat the flow. The contrast of shocked flow behavior depending on e is also illustrated in Fig. 9, which shows the whole φ-averaged in sector s1 v r versus T at R = 400 a for different eccentricities. The velocity jumps grows with eccentricity and became more pronounced for e ≥ 0.5.
The radial velocity, averaged over R in the interval [350 a, 435 a], versus φ is shown in Fig. 10. The radial velocity is also color mapped on the φ (y-axis) versus T (x-axis) plane in Fig. 11. Consistently with previous figures, one sees that the low e cases present a quasi-isotropic flow with significant stochastic behavior on top of the longer term behavior, whereas high e-cases show faster flows concentrate around the periastron-apastron direction (φ = π).

Summary and discussion
In this article, we have presented a detailed analysis of the results obtained by the simulation of colliding pulsar and stellar winds. For the first time, cases with different eccentricities are systematically explored, so the impact of this prediction can be assessed. The further development of these models should impact the studies of very powerful sources such as HMPB or microquasars (Dubus et al., 2010;Zdziarski et al., 2018;Sinitsyna and Sinitsyna, 2021;Massi et al., 2020). We note that the first simulation of a jet-stellar wind interaction in a microquasar along a full orbit has been already done in Barkov and Bosch-Ramon (2022).
The evolution of the shocked pulsar and stellar winds was studied analytically in Bosch-Ramon and Barkov (2011). The mixed-wind eventual velocity       away (expel) from the binary can be estimated as: where η −1 = η/10 −1 , and our numerical results, v exp = 0.40±0.03 c (see Fig. 4), confirm this prediction. In our numerical calculations, we find that the shocked flow structure on large scales is a slowly-accelerating, supersonic mixture of shocked stellar and pulsar winds, with an approximately isotropic propagation in φ. We note that the simulations could not properly probe the expansion in θ due to low resolution, although 3D simulations by (Bosch-Ramon et al., 2015) suggest that on scales a the shocked structure can become wider in that direction than conical expansion predicts because of internal energy confinement. The present work also shows that this mixed supersonic wind is very clumpy in density and velocity, so particle acceleration may easily occur in such an environment.
Our results are fully consistent with those obtained by Barkov and Bosch-Ramon (2016); Bosch-Ramon et al. (2017), who showed for the first time the dramatic impact that eccentricity can have on the evolution of the shocked flows on large scales (non-thermal processes were discussed in Barkov and Bosch-Ramon 2018). In addition to that, the present work also characterizes the eccentricity of a transition between an approximately isotropic supersonic wind, made of shocked stellar and pulsar wind (for small e), and a sort of two-component structure, one fast, light and collimated, directed along the periastron-apastron direction, and the other slow, dense and broad, directed elsewhere (for high e). In particular, we find that the transition is somewhere between e = 0.5 and 0.75, probably close to the latter.
As mentioned in Sect. 2, our results should not be sensitive to the orbital period on large scales. This implies that at a semi-quantitative level, our predictions can be extrapolated to wider systems. However, higher accuracy in the estimate of the eccentricity associated to a structure geometry transition requires fully 3D calculations. On the other hand, the magnetic field could also play an important role in the evolution of the shocked flows, and should be included in future stages of this research. Moreover, future numerical work should tackle the issue of how the shocked mixed flows interact with the interstellar medium, both in the low e (see Bosch-Ramon, 2011, for analytical predictions) and the high e regimes. Finally, the consequences of the eccentricity dependence of the shocked flow evolution with respect to non-thermal emission should be studied in more detail than what has been done so far.

Data availability
The original data and its analysis can be requested by email.