© The Authors 2023

1 Introduction

Gamma-ray binaries are binary systems that consist of a massive star and a compact object and emit the major part of their observable radiative output in the gamma-ray regime. A common approach to explaining the observed non-thermal emission from some of these systems is the wind-driven scenario (see Dubus 2013, and references therein). In this case, the compact object is assumed to be a pulsar, where the highly relativistic pulsar wind forms a wind-collision region (WCR) when interacting with the massive stellar wind from the star. At the strong shocks surrounding the WCR, particles are thought to be injected and accelerated to very high energies (Maraschi & Treves 1981; Dubus 2006).

The most studied gamma-ray binary is the LS-5039 system, for which a plethora of observations is available at different energies. Additionally, the properties of this system are comparatively well known: it contains an o-type star and a compact object in a mildly eccentric orbit around each other. For this system, many studies assume the wind-driven scenario, where a small fraction of the rotational energy of the pulsar drives a pulsar wind interacting with the wind of the o-type star. This scenario is supported by the orbital modulation observed in the non-thermal emission at different energies (see, e.g. Aharonian et al. 2005; Abdo et al. 2009; Takahashi et al. 2009). Additionally, there is tentative evidence of pulsations in X-rays from the analysis of Suzaku and NuSTAR data by Yoneda et al. (2020), which also hint at the presence of a pulsar as the compact object. While Volkov et al. (2021) have discussed the significance of the periodicity found by Yoneda et al. (2020) being rather low, they nonetheless argue from their observation of the light curve of LS 5039 that the colliding-wind scenario still is the most likely one.

In this scenario, the non-thermal emission from LS 5039 is strongly linked to the complex dynamics of the colliding pulsar and stellar winds within the WCR and beyond. Correspondingly, this dynamical interaction has already been investigated in many numerical studies. Due to the high numerical cost of a fully rel-ativistic description, early works focussed on individual aspects of such gamma-ray binary systems, where certain simplifications became feasible, such as the neglect of a relativistic description (Romero et al. 2007; Takata et al. 2012), neglect of orbital motion (Bogovalov et al. 2012; Lamberts et al. 2013; Dubus et al. 2015), or simulations with reduced dimensionality (Bosch-Ramon et al. 2012; Bogovalov et al. 2012; Paredes-Fortuny et al. 2015).

The first (three-dimensional) 3D relativistic simulations of LS 5039 were discussed in Bosch-Ramon et al. (2015), where the authors considered the evolution of instabilities and the difference between (two-dimensional) 2D and 3D simulations (see also Bosch-Ramon et al. 2017; Barkov & Bosch-Ramon 2016, for relativistic 3D simulations of other gamma-ray binaries). By using a radially increasing cell size, the authors were able to simulate a larger computational domain, while still maintaining a high resolution near the apex of the WCR. As in the previous 2D study by Bosch-Ramon et al. (2012), Bosch-Ramon et al. (2015) observed the spiral pattern caused by the orbital motion, which starts being disrupted at larger radii due to the turbulence driven in the WCR, as has also been observed in colliding-wind-binary simulations (Lamberts et al. 2012). However, due to the radially decreasing resolution, turbulence was attenuated in the outer parts of the simulations by Bosch-Ramon et al. (2015). While Huber et al. (2021b) used a homogeneous resolution through their simulated domain, those authors found that the corresponding domain was too small to contain the unshocked pulsar wind at all orbital phases.

In this work, we extend the previous studies by investigating a simulation of LS 5039 with a homogeneous high spatial resolution throughout the entire numerical domain. In particular, we increased the size of the numerical domain, as compared to Huber et al. (2021b), and further increased the spatial resolution to reduce the dissipation of turbulence in the simulation. Additionally, we simulated the dynamics of the system for three full orbits, where the first orbit is only used to allow the system to settle into a quasi steady state. With this, we are able to observe the impact of turbulence also in the outer parts of the numerical domain and not only quantify the short-term variation but also the orbit-to-orbit variability of LS 5039.

The paper is structured as follows. In Sect. 2, we introduce the specific set-up of our simulation together with the mathematical description. The results are discussed in Sect. 3, where we start by investigating the flow structure and the dynamics, ending with an analysis of the short-term and orbit-to-orbit variations. Finally, we summarise our findings in Sect. 4.

2 Physical and numerical set-up

Since we simulated the interaction of a highly relativistic pulsar wind with the massive wind of an early-type star, we used relativistic hydrodynamics (RHD) to describe the wind dynamics (for a more extensive discussion of the model, see Huber et al. 2021a). In particular, we solved the following set of partial differential equations: $$\frac{\partial D}{\partial t}+\nabla \cdot \left(\frac{1}{\gamma }Du\right)=0,$$(1) $$\frac{\partial \tau }{\partial t}+\nabla \cdot \left(\frac{1}{\gamma }\left(\tau +p\right)u\right)=S_{\tau },$$(2) $$\frac{\partial m}{\partial t}+\nabla \cdot \left(\frac{1}{\gamma }mu\right)+\nabla p=S_m,$$(3)

where we have the vector U of conserved quantities $$U=\left(\begin{array}{c}D\\ \tau \\ m\end{array}\right)=\left(\gamma \rho \underset{Dhu}{\overset{\gamma \rho }{\left(\gamma h-1\right)}}-p\right).$$(4)

Here, ρ is the mass density, u is the spatial vector of the rel-ativistic four velocity with иⁱ = γυⁱ, h is the specific enthalpy, and p is the thermal pressure. Additionally, τ is related to E via E = τ + D. With our normalisation c = 1, υ is given in units of the speed of light and we find for the relativistic Lorentz factor: $$\gamma =\frac{1}{\sqrt{1-{\upsilon }^2}}=\sqrt{1+u^2},$$(5)

namely, the Lorentz factor can be directly computed from the spatial components of the relativistic four velocity.

The system of equations given in Eqs. (1)–(2) is closed by an ideal equation of state: $$h\left(\rho ,p\right)=1+\frac{\text{Γ}}{\text{Γ}-1}\frac{p}{\rho },$$(6)

where we use a constant adiabatic exponents Г = 4/3 (see Huber et al. 2021a, for a related discussion). Here, we did not solve for E as is often done in other simulations (see, e.g. Mignone et al. 2007), since D and E can become very similar in the low-pressure, non-relativistic regime. With both D and E being conserved quantities, τ is also conserved, which is more suitable in our case especially for the treatment of the stellar wind. Additionally, we solved a supplementary conservation equation for the specific entropy s = p/ρ^Г, which is only conserved in smooth flows where it can be used to overcome possible problems related to negative pressure (see Huber et al. 2021a).

These equations were solved using the CRONOS code (Kissmann et al. 2018), which was recently extended to allow simulations of relativistic hydrodynamics (see Huber & Kissmann 2021). In this work, we applied piecewise linear spatial reconstruction using the minmod limiter (van Leer 1979; Harten et al. 1983) together with an HLLC RHD Riemann solver (Mignone & Bodo 2005). The simulation analysed in this work was run on the Joliot-Curie system at GENCI@CEA, France.

In our model of LS 5039, we simulated the wind dynamics of the system, consisting of a pulsar emitting a relativistic pulsar wind in orbit with a massive O-type star, for three full orbits. Here, we used the same physical parameters for the system as given in Huber et al. (2021b): we adopted the orbital parameters from Casares et al. (2005) using a period of P_orbit = 3.9 d and eccentricity of e = 0.35, along with masses of M_star = 23 M_⊙ and M_pulsar = 1.4 M_⊙ (see also Dubus 2013).

As discussed in the introduction, we increased both the extent of the numerical domain and also the spatial resolution in comparison to our previous simulations of LS 5039. With the centre of mass at the coordinate origin, we used a numerical domain with dimensions [−2, 2] × [−0.5, 2.5] × [−1, 1] AU³ that is homogeneously filled with 2048 × 1536 × 1024 cubic cells. In total, this leads to a resolution of ~0.002 AU or ~0.42 R_⊙ in each spatial dimension, namely, about the same resolution as used in the inner region by Bosch-Ramon et al. (2015) or a doubling of spatial resolution together with a significant increase of the simulation volume in comparison to the previous study by Huber et al. (2021b). In the present study, we used homogeneous resolution throughout the numerical domain. Adaptive mesh refinement was not considered, since small-scale fluctuations filled most of the computational volume (see also Bosch-Ramon et al. 2015).

In our simulation, the z-axis is perpendicular to the orbital plane. As before, the simulation was performed in a reference frame co-rotating with the average angular velocity of the system, which allowed us to use a smaller non-symmetric domain than in a non-rotating set-up that contains the full unshocked pulsar wind at all times of the simulation. By using the so-called 3+1 Valencia formulation of general relativistic hydrodynamics (Banyuls et al. 1997), this requires using the source term: $$S_m={\left(\text{Ω}{m}_y,-\text{Ω}{m}_x,0\right)}^{\top },$$(7)

where Ω is the average angular velocity of the binary (see Huber et al. 2021a, for further details). Since we are dealing with an eccentric orbit in the case of LS 5039, the two stars still show residual motion in the system rotating with the average angular velocity. This is apparent by the axis connecting the binaries not being aligned with the coordinate axes for most of the phases depicted, for instance, in Fig. 1. Periastron is marked by phase Φ = 0, where the axis connecting the binaries is aligned with the y-axis of our simulation volume.

In the present case, the simulation was performed for 3 full orbits of the system, where the analysis focusses on the second and the third orbit (the first orbit was used to allow the system to change from the initial homogeneous low-density medium to the fully turbulent state). Here, we allowed more time for this initialisation phase than in our previous study because of the larger spatial domain. In the present simulation, the undisturbed wind of the massive star reaches the farthest corner of the numerical domain at t_crossing ≃ 0.75 P_orbit ≃ 2.9 d. This is actually a more conservative estimate than in Bosch-Ramon et al. (2015), who initialised most of the domain with the expanding stellar wind, allowing ≃0.5 d for initialisation.

For the stellar wind, we used a stellar mass-loss rate of Ṁ_s = 2 × 10⁻⁸ M_⊙ yr⁻¹ together with a terminal velocity of v_s = 2000 km s⁻¹ (Dubus et al. 2015). The stellar wind was injected as an isotropic outflow in a sphere around the position of the star with radius r_inj = 0.012 AU, that is, with a radius of six computational cells, where the star moves in the co-rotating frame due to its eccentric orbit. The pulsar wind was similarly injected with a speed of v_p = 0.99c, with c as the speed of light. For the pulsar’s mass-loss rate we used ${\dot{M}}_{\text{p}}=\eta \frac{{\dot{M}}_{\text{s}}{\upsilon }_{\text{s}}}{u_{\text{p}}}$ with η = 0.1 as also used in Bosch-Ramon et al. (2015). With this, we assume a pulsar spin-down luminosity of L_SD = 7.55 × 10²⁸ W. The thermal pressure at the outer radius of the injection sphere was set to be a fraction of 10⁻⁹ and 10⁻⁷ of the local rest-energy density for the star and the pulsar, respectively.

For the analysis, we stored 20 output files evenly distributed over orbital phase for each orbit. For a little more than half of these outputs after the first orbit, we produced a series of six successive output datasets with only about 8.3 min between the individual datasets. These can be used to investigate short-term fluctuations and compute statistics of relevant quantities.

Fig. 1 Mass density in the orbital plane during the second full orbit. Corresponding phases are given in the individual sublots together with an arrow towards the observer location. The centre of gravity is at x = 0, y = 0. For the central plot at the top, the black box indicates the region for which we performed the turbulence analysis.

3 Results

In Fig. 1, we show the density in the orbital plane for six different phases of the second simulated orbit. In these figures, the dark and mostly featureless region is the radially expanding supersonic pulsar wind, where the bright dot marks the local density peak at the location of the pulsar. The homogeneous region to the left and lower-left contains the unshocked, radially expanding stellar wind, where the star is visible as the large bright circle near (x, y) = (0, 0). Between the star and the pulsar, we can see the wind-collision region (WCR), where both winds interact. on both sides, this WCR is bounded by a shock, where the wind-component parallel to the shock normal becomes subsonic. For the stellar wind, this termination shock shows a spiral pattern owing to the orbital motion of the stars. This becomes even more apparent at larger scales as investigated in Bosch-Ramon et al. (2015). The termination shock of the pulsar wind can be split into a u-shaped bow shock and a so-called Coriolis shock formed due to the orbital motion of the system (Bosch-Ramon & Barkov 2011; Bosch-Ramon et al. 2012), where the latter is also a possible site for particle acceleration (Zabalza et al. 2013; Huber et al. 2021b).

In the downstream regions of the WCR, Fig. 1 shows space-filling turbulence that is driven by different instabilities within the WCR. On the one hand, the shear-flow between the shocked stellar and pulsar winds at the contact discontinuity within the WCR triggers Kelvin-Helmholtz (KH) instabilities (Bosch-Ramon et al. 2012, 2015; Lamberts et al. 2013). Additionally, the Richtmyer-Meshkov (RM; Richtmyer 1960; Meshkov 1972) and Rayleigh-Taylor (RT; Rayleigh 1882; Taylor 1950) instabilities are active in this highly dynamical environment, where Bosch-Ramon et al. (2015) discuss the fact that they act together as an important driver of the turbulence in these systems. As is also visible in our Fig. 1, turbulence is much stronger at the leading edge of the WCR (to the left of the cavity containing the unshocked pulsar wind).

This is a strong indication that the RT instability is acting as one of the important drivers of turbulence in these systems, because acceleration (by the Coriolis force) and density gradient only point in the opposite direction at the leading edge of the WCR (see also Bosch-Ramon et al. 2015). While Guzdar et al. (1982) pointed out that a velocity shear will modify the RT instability, they only find a significant reduction of the growth rate for wavelengths shorter than the scale of the density gradient, which in our case is the grid scale. At large scales, where the relevant driving occurs in our simulations, the RT instability is expected to be unaffected by the velocity shear at the contact discontinuity.

The presence of the RT and the RM instabilities is particularly important, since our simulation cannot capture ultra-relativistic Lorentz factors as conventionally expected for pulsar winds (Aharonian et al. 2012; Khangulyan et al. 2012). Since larger Lorentz factors lead to a correspondingly smaller mass-loss rate Ṁ_p for the pulsar and thus a larger density contrast, we expect longer growth timescales for the KH instability (Ferrari et al. 1980), which is also observed in the classical limit (see, e.g. Chandrasekhar 1961). In contrast, the growth rate of the RT instability depends on the difference of the densities of the stellar and the pulsar wind (see, e.g. Eq. (7) in Allen & Hughes 1984), where the density of the pulsar wind is already negligible in our set-up – and even more so for higher Lorentz factors. Thus, we expect little change with respect to higher Lorentz factors for the RT instability.

In this work, we feature a larger Lorentz factor than Bosch-Ramon et al. (2015), where the turbulence still shows the same general structure in our simulation: we also find ubiquitous turbulence and corresponding mixing downstream of the WCR. Also, in our simulation, we find that the instabilities within the WCR drive the mass loading, namely, dense matter from the stellar wind is mixed with the dilute, high-velocity pulsar-wind material (see also Bosch-Ramon et al. 2012). This is visible in the dense clumps of shocked stellar-wind material embedded in this region, where in our case the mixing already seems stronger near the apex of the WCR than observed by Bosch-Ramon et al. (2015). Additionally, the fluctuations of the WCR are much stronger than in Bosch-Ramon et al. (2015), where (in particular) the size of the region containing the unshocked pulsar wind varies more strongly. This difference can be attributed at least to two effects: first, we use a slightly larger eccentricity in our simulation and, secondly, due to the homogeneous resolution in our simulation, turbulence levels in the outer parts of the numerical domain are not damped and can therefore lead to a disruption of the structure of the supersonic pulsar wind volume.

Fig. 2 Different snapshots of mass density in the orbital plane during a short period in the second full orbit as indicated in the plots. As given by the indicated orbital phases, these snapshots cover a period of slightly more than 40 min.

3.1 Coriolis-shock region

The position and structure of the Coriolis shock in our simulation is subject to dynamical changes, as visualised in Fig. 2. There, we show a zoom into the region containing the unshocked pulsar wind over a period of slightly more the 40 min. Apparently, the distance of the Coriolis shock from the apex of the WCR changes significantly, where it is absent for parts of that period, namely, the flanks of the WCR converge into a point instead of connecting to both ends of the Coriolis shock. This will also have implications for particle acceleration at this shock, where a change in distance from the star will also change the energy loss rates of particles at the shock.

Further peculiarities of the shocks in the wind-collision region are visible in Fig. 3, where we show the absolute value of the spatial component of the relativistic four velocity of the fluid и = vγ. Apparently, sometimes the Coriolis shock marks the transition to a region with the highly turbulent medium containing a mixture of stellar and pulsar wind, while at other phases it is embedded within the relativistic fluid from the wings of the WCR propagating between the unshocked pulsar wind and the contact discontinuity, or in a mixed state (right plot in Fig. 3).

This rather laminar fluid in the wings of the WCR, contains streamlines connecting to the apex of the WCR, where the fluid is re-accelerated by the large pressure gradient (Bogovalov et al. 2008), and streamlines connecting to the flank of the pulsar-wind termination shock, where the velocity component perpendicular to the shock normal is still highly relativistic – in total we find Lorentz factors up to ∼3 in this region. In all cases, this flow also terminates at an extended shock, when entering the region of highly turbulent, mixed matter downstream of the WCR. The material in this region of mostly laminar flow correspondingly still moves supersonically, as can be seen in Fig. 4, where we plot the Mach number in the orbital plane during the second full orbit. The highly turbulent matter beyond the extended shock finally features large fractions of subsonic gas. Also beyond the curved termination shock of the stellar wind, we see large regions of turbulent low-Mach-number flow. Only the trailing-edge region beyond the stellar-wind termination shock features supersonic flow that later terminates in an extended shock. The turbulence beyond these extended shocks is also obvious via the strong fluctuation in the Mach number.

At some times, the region of laminar shocked material is also permeated by large-scale discontinuities, as visible at orbital phases Φ = 1.310 and Φ = 1.811 in Fig. 3. As an example, we show mass density and thermal pressure for Φ = 1.310 in Fig. 5, which feature the same large-scale discontinuity as the flow velocity. Apparently, these discontinuities are shock waves connected to eddies in the contact discontinuity separating the shocked stellar wind from the shocked pulsar wind. By visual inspection it seems that the fluctuations increase beyond the point, where the large-scale shocks traverse the contact discontinuity. Since we do not have outputs at sufficiently short time intervals available to analyse the short-timescale evolution of these fluctuations, this can only be viewed as a hint that the RM instability might be responsible at least for parts of the fluctuations visible in our simulations (see also Bosch-Ramon et al. 2015). Apart from that, it will be interesting to see how relevant these shock structures are for particle acceleration and the subsequent non-thermal emission, which has not been studied so far.

Fig. 3 Absolute value of the spatial component of the fluid’s four velocity in the orbital plane at selected phases as indicated in the plots. Superimposed over the images, we indicate the direction of the flow via streamlines. Here, we only show results for a part of the orbital plane, focussing on the region around the Coriolis shock. We note the different x-axis in the right-handed plot.

Fig. 4 Mach number in the orbital plane during the second full orbit, for orbital phases as indicated in the subplots. The arrows indicate the direction towards the observer.

Fig. 5 Pressure (left) and mass density (right) in the region surrounding the unshocked pulsar wind for phase Φ = 1.310.

Fig. 6 Distribution function of the relativistic gamma factor fγ within the orbital domain at different phases as indicated in the plots. Here, we show results for the second orbit (left), for the third orbit (middle), and also for the same short time interval as visualised in Fig. 2 (right).

3.2 Downstream flow and turbulence

Beyond the extended shock structures we observe a highly turbulent downstream region embedded on the left in the unshocked stellar wind. Here, the mixing of stellar- and pulsar-wind material is obvious both in Figs. 1 and 3, where we observe high density clumps from the stellar wind or regions with a high Lorentz factor from the pulsar wind. In our case, this mixing shows up already close to the apex of the WCR, showing the efficiency of the instabilities acting to produce the turbulence. Corresponding statistics is visible in Fig. 6, where we show the distribution function of the Lorentz-factor within the entire numerical domain for different orbital phases.

Given the high resolution together with the homogeneous grid in our simulation, we also investigated the longitudinal structure functions as a measure for the turbulence in our computational domain. Since only the downstream medium shows strong fluctuations, we restricted the corresponding analysis to this region. In particular, we extracted the spatial part of the relativistic four-velocity field from a region with an extent of 1 AU in all spatial dimensions, giving a 512 × 512 × 512 data cube. This subdomain is centred around the orbital plane and located at x = 0… 1 AU, y = 1.36… 2.36 AU, as indicated via the black box in Fig. 1. The longitudinal structure functions were computed according to (see, e.g. Zrake & MacFadyen 2013): $$S_{\text{p}}^{\parallel }\left(l\right)=\langle {\left|\frac{1}{l}\delta u^{\mu }\delta x_{\mu }\right|}^p\rangle ,$$(8)

where $\delta u^{\mu }=u_1^{\mu }-u_1^{\mu }$ and l = (δx^µδx_µ)^1/2 with $\delta x^{\mu }=x_2^{\mu }-x_1^{\mu }$ the separation four vector between pairs of points chosen randomly. Here, we selected points that were simultaneous in the co-rotating frame.

Examples for corresponding structure functions are shown in Fig. 7. We observe only a tiny inertial range (if any) at scales around 0.07 AU, where S3 in most phases shows a linear dependence on 1. Especially for phases just around apastron, the structure functions often show rather erratic behaviour, where we give one of the more well-behaved examples in Fig. 7.

In contrast to homogeneous turbulence simulations, where structure functions are often analysed and show typical behaviour given by the dissipative structures (She & Leveque 1994), the driving in our simulations is not constrained to the largest spatial scales. Instead it follows from the different instabilities, which drive the fluctuations at a range of spatial scales. Therefore, we do not have an extended inertial range in our simulations, where only the effect of the turbulent cascade dominates. Even though most of this driving occurs in the vicinity of the contact discontinuity, the corresponding fluctuations are advected into our analysis regions. This becomes particularly important for orbital phases around and after apastron, where the region filled by unshocked pulsar wind becomes rather large, thus shifting the region with the turbulent driving closer to the turbulence-analysis region. In some of these phases additionally, parts of the unshocked pulsar wind extend into the volume, where we analyse the fluctuations. Despite all this, the structure functions are related to the amplitude of fluctuations within the orbital domain, where we observe strong differences between the different phases depicted in Fig. 7. This strong variability also observed in the other previous plots, motivates an analysis of the short- and long-term variability of different quantities in this system.

3.3 Short-term and orbit-to-orbit variations

An important property of our new set of simulation results is the long timescale considered, here. Apart from the initialisation phase in the first orbit, we considered the evolution of the system for two further full orbits. Therefore, we are in a position to quantify dynamical changes during a single orbit as well as orbit-to-orbit variations.

As a first relevant variable that can be easily obtained from the simulation results, we investigated the size of the volume filled with unshocked pulsar wind. This was obtained by adding the volume of all cells with a four velocity fulfilling γ > 7. As can also be seen in Fig. 3, outside the unshocked pulsar wind, such speeds are not achieved despite the re-acceleration of shocked pulsar wind in the flanks of the WCR.

Figure 8 shows some remarkable features regarding the stability of the unshocked pulsar-wind structure. First, its size apparently is very stable right after periastron. After that, the size expectedly grows, but the fluctuations become very large, especially after apastron. As a measure of the fluctuations, we computed the standard deviation of volumes extracted from six consecutive output files. A visual inspection of Fig. 8 shows that this estimate sometimes underestimates the actual variation in our simulations, where strong changes occur at times, for which we have no output data available.

Considering that this is volume and not linear distance, length scales will show correspondingly smaller variations, but it can still be expected that emission from this system will show significantly higher fluctuations around and in particular after apastron. From the given averaged values we find the smallest volume of the unshocked pulsarwind to be 0.0015 AU³ and the largest, the peak in the third orbit, 0.19 AU³. Converting this to a linear distance, this corresponds to a factor of more than 5 in spatial extent, which is nicely correlated with the expectation, when considering Fig. 1.

In Fig. 9, we show the orbital variation of the distance of bow and Coriolis shocks from the pulsar. These were computed by following the direction connecting star and pulsar from the pulsar until the Lorentz-factor is smaller than 7. As expected, the distance of the bow shock from the pulsar varies smoothly over the orbit – fluctuations of this distance are small and the variation corresponds to the variation of the orbit. At all phases, the distance of the bow shock is at a similar fraction of the approximate distance of the contact discontinuity: $$d_{\text{CD}}=\frac{\sqrt{\eta }}{1+\sqrt{\eta }},$$(9)

as given by Eichler & Usov (1993). In contrast, the distance of the Coriolis shock varies significantly and not always in accordance with the orbit. Again, strong variations around apastron are present, but the variation can also be quite large around periastron. This again reflects the sudden displacements of the Coriolis shock as shown in Fig. 2 and discussed in Sect. 3. The implications for gamma-ray emission from the particle population related to the Coriolis shock will be investigated in a future study.

The distribution of the Lorentz factor in Fig. 6 also shows a distinct orbital variation related to the variation of turbulence. In this plot, we can also see the variation of the volume of the unshocked pulsar wind, as discussed above, from the changing magnitude of the peak near γ = 7. The varying fraction of unshocked stellar wind at low gamma factors simply follows through the motion of the stars within the co-rotating domain, which leads to different volumes of unshocked stellar wind, as also visible in Fig. 1. The imprint of turbulence can rather be seen at intermediate gamma factors. Also, the changing volume due to the moving stars has some impact here, but apart from that (especially during the second orbit), the 4 < γ < 7 regime shows changes related to the dynamics of turbulence, mixing, and dynamical changes of the region around the pulsar-wind termination shock. Apparently, for phases just after apastron, the fraction of shocked pulsar wind in the computational domain is much larger, leading to higher contributions of the intermediate gamma values. Especially at these phases, the unshocked pulsar-wind volume is rather large with a surrounding region of sometimes smooth flow (see Sect. 3.1 and the right-handed figure of Fig. 3), where reacceleration and shocks with large angles between shock normal and pulsar-wind velocity can lead to high gamma factors in this post-shock region.

This behaviour is qualitatively also visible in the third orbit (see middle plot in Fig. 6), but it is less pronounced in this case. This hints at strong orbit-to-orbit variability of the turbulence within this system. Apart from that, we illustrate the short-time variability of the velocity in the right-handed plot of Fig. 6. Apparently, distinct changes of the distribution function of the gamma factor occur even on timescales of a few minutes. In this case, the distribution decreases with time in the region 4 < γ < 7. This is in accordance with the observations in Fig. 2, where we see that during this time, a large region filled with a smooth high-gamma-factor wind is disrupted and filled with slower, more turbulent matter.

From the distribution of the Mach number in the computational domain, we additionally computed the fraction of supersonic flow in the downstream region. For this, the unshocked stellar wind was identified by its small speed of sound and the unshocked pulsar wind by its large Lorentz factor. Both unshocked winds were excluded from this analysis since the flow is highly supersonic in each case. As can be seen in Fig. 10, the fraction of downstream medium that moves supersonically in our simulation depends strongly on orbital phase, with a clear peak after apastron. This behaviour is also visible in Fig. 4, where phases near and after apastron feature large regions with high Mach-number flow, especially relating to the smooth flow surrounding the unshocked pulsar wind.

Since neither the mere amplitude of the velocity as given by the distribution of the gamma factor nor the Mach number of the flow do directly give the turbulence amplitude, we additionally computed an estimate for the mean rate of inertial energy dissipation from the third-order longitudinal structure functions S₃ Using Kolmogorov’s four-fifth law: $$S_3=-\frac{4}{5}ϵl,$$(10)

ϵ was computed by integrating S₃ over the range of scales, where S₃ follows the expected S₃ ∝ l dependence (see Fig. 7), and then solving for ϵ. Its orbital dependence shown in Fig. 11 is similar to the previously discussed ones, with peaks around apastron and strong variation around the same time. In particular around apastron, we observed several phases, where the region in which we evaluate the structure functions is permeated by parts of the unshocked pulsar wind. This leads on the one hand to non-turbulent regions within the analysis volume. On the other hand the contact discontinuity is close by the analysis region and is part of it in some cases (more precisely the phases near Φ = 1.4− 1.6, 2.55−2.6 are affected). This also means that the instability driving the turbulence is active within the analysis volume at scales, which are investigated for fully developed turbulence. As a result, the structure functions are not representative for fully developed turbulence. Nonetheless, in nearly all cases, we find S₃ ∝ l in the regime from which we compute ϵ. Therefore, we still used the result as a measure of the fluctuations inside the relevant volume. We observe strong short-term fluctuations in ϵ, especially for the phases around apastron.

Apart from the short-time variation discussed so far, we also observed the orbit-to-orbit variations. In Figs. 8–11, we find that the phase dependence around periastron is very similar in the second and the third orbit. Around apastron, however, the short-term fluctuations are stronger than the orbital ones leading to different phase dependencies in the two investigated orbits. This can also be seen when comparing the mass density in the orbital plane in the third orbit shown in Fig. 12 to the one in the second orbit shown in Fig. 1. It is especially for the shape of the unshocked pulsar-wind region and the immediate surroundings that we find strong differences after apastron (see the density distributions at Φ = 1.6 vs. Φ = 2.6 or at Φ = 1.75 vs. Φ = 2.75). This corresponds to the times when we observe large fluctuations in the volume of the unshocked pulsar wind in Fig. 8.

This might also have implications for the variability of non-thermal radiation observed from the system. From our investigations of the dynamics of the system, we expect stable gamma-ray emission around periastron with possibly strong variations around apastron. Huber et al. (2021b) found that the gamma-ray emission at the time shortly after apastron is particularly affected by relativistic beaming related to the emission of energetic particles within the leading edge of the shocked pulsar wind. From Figs. 1 and 12, we see that the flow direction in this region is aligned with the direction to the observer sometime around relative phase Φ = 0.6, where in the current simulation this situation seems to occur a little later than in the simulations by Huber et al. (2021b). However, as also discussed in Huber et al. (2021b) this phase was not well represented in their simulations due to the limited size of their numerical domain.

This relativistic beaming occurs just in the region, where we have the rather laminar flow within the leading edge of the shocked pulsar wind, which is also visible as the dark yellow region in Fig. 3. That same figure shows that at some phases, there even appear two distinct such regions with different overall flow directions. Regarding the impact on gamma-ray emission, only the relative phases around Φ = 0.6 when the flow is aligned with the observer direction will be important. Due to the rather different structure of this region at phases Φ = 1.6 and Φ = 2.6, however, we can expect to find a different impact of relativistic beaming for the gamma-ray emission in the second and the third orbit. This will be studied in more detail in a forthcoming publication, where we will investigate the propagation of energetic particles within the simulated stellar winds, discussed here. There, it will be particularly interesting if we find short-term and orbit-to-orbit variations on a similar level as we observe in the present study.

Fig. 7 Structure functions at four different orbital phases as indicated in the figures. For comparison, we show a linear function by the black dotted line. By our normalisation, the structure functions are given in units of cp.

Fig. 8 Volume of unshocked pulsar wind Vpuisai as a function of orbital phase Φ. Data are shown for the full second and third orbit. The error bars indicate the corresponding standard deviation, that was computed from multiple output files at similar orbital phases, where we averaged over six phases spread over slightly more than 40 min.

Fig. 9 Distance of the Coriolis shock (blue plus signs) and bow shock (red crosses) from the location of the pulsar as functions of orbital phase, where distances are given along the line connecting star and pulsar. Average values and error bars were determined as stated in Fig. 8. The y-ахis uses logarithmic units to allow showing both shock distances simultaneously. Additionally, the dashed orange line shows the approximate distance of the contact discontinuity from the pulsar in the direction of the star using the formula by Eichler & Usov (1993).

Fig. 10 Fraction of supersonic downstream medium as a function of orbital phase. Average values and error bars were determined as stated in Fig. 8.

Fig. 11 Variation of the mean rate of inertial energy dissipation per unit mass ϵ as a function of orbital phase Φ over the second and the third orbit. As in Fig. 8, the error bars indicate the standard deviation computed from six consecutive output flies, ϵ is given in units of с3 AU−1.

Fig. 12 Same as Fig. 1 but for the third full orbit.

4 Summary and discussion

In this study, we investigated the dynamical interaction of the pulsar wind and the stellar wind in the gamma-ray binary system LS 5039 via RHD simulations, where we assumed a wind-driven scenario to explain the observed gamma-ray emission. Here, we did not take the magnetic field in the winds into account, as previous 2D simulations led to the expectation of a rather low impact on the wind dynamics (see the discussion in Bogovalov et al. 2012; Bosch-Ramon et al. 2015). Recent 3D relativistic magnetohydrodynamics (RMHD) simulations (Barkov et al. 2022), however, have shown that for future simulations it will also be interesting to include the effects of the magnetic field. The investigated simulation was done with unprecedented resolution over three full orbits to allow for a detailed analysis of the turbulence driven by the wind-wind interaction in the WCR, together with the short-term and long-term variability of the wind dynamics.

In the simulation, we observed strong turbulence in the downstream region of the WCR, where we did not see a clear, broad inertial range due to our driving force by the different instabilities being active at a range of spatial scales. Here, we did not take clumping for the stellar wind into account, which can be shown to further increase fluctuations (Paredes-Fortuny et al. 2015; Kefala & Bosch-Ramon 2023). Thus, short-term and orbit-to-orbit variation might actually be even stronger than found in our model. However, it is important to be aware that the Lorentz factor of the pulsar wind, chosen in the present and similar simulations, was significantly smaller than expected in a realistic pulsar wind. While this is expected to diminish growth rates for the KH instability (Ferrari et al. 1980), the RM and RT instabilities, which seem to be responsible for a large part of the observed turbulence should be much less affected by a larger Lorentz factor and a corresponding change in density contrast (Allen & Hughes 1984; Bosch-Ramon et al. 2015).

In our configuration, we find strong variability due to the turbulence driven within the WCR. While different parameters turn out to be stable around periastron, we observe particularly strong fluctuations both on short and on orbit-to-orbit timescales around apastron. These fluctuations seem to be stronger than found in previous simulations as in Bosch-Ramon et al. (2015), for instance. This might be partly attributed to the higher resolution, especially in the outer parts of the numerical domain, but also to the somewhat higher eccentricity used in our simulation.

In a future analysis, we will further investigate how the strong dynamics observed here will impact energetic particle transport and ensuing emission of non-thermal radiation from this system. Again the phase shortly after apastron will be particularly interesting, because this is the time when we observed the largest fluctuations and it is also when relativistic beaming in the direction of the observer can enhance emission from within the WCR.

Movie

Movie provided by the authors Access here

Acknowledgements

We thankfully acknowledge PRACE for granting us access to Joliot-Curie at GENCI@CEA, France. We thankfully acknowledge the access to the research infrastructure of the Institute for Astro- and Particle Physics at the University of Innsbruck (Server Quanton AS-220tt-trt8n16-g11 x8). This research made use of Cronos (Kissmann et al. 2018); GNU Scientific Library (GSL) (Galassi et al. 2018); FFTW3 (Frigo & Johnson 2005), matplotlib, a Python library for publication quality graphics (Hunter 2007); Scipy (Virtanen et al. 2020); and NumPy (Harris et al. 2020).

Appendix A Supplementary material

For a further illustration of the dynamics of the system in our simulation, we supply a video of the evolution of the mass density in the orbital plane of the LS-5039 system. The video was produced from the series of six successive outputs for different orbital phases during the full second and third orbit. Since these outputs are not homogeneously distributed over all phases of the orbit, the video continuously jumps from phase to phase, where it shows a brief moment of the short-time dynamics of the system.

This video is useful in showing several of the effects discussed in this paper: the vast difference in speed of the stellar and the pulsar wind becomes directly obvious from the very different short-time motion of the material. Additionally, the highly dynamical changes of the shape and the size of the volume filled by unshocked pulsar wind become very apparent. Finally, this video alone can clearly show the large-scale motion of the material over longer timescales.

References

All Figures

	Fig. 1 Mass density in the orbital plane during the second full orbit. Corresponding phases are given in the individual sublots together with an arrow towards the observer location. The centre of gravity is at x = 0, y = 0. For the central plot at the top, the black box indicates the region for which we performed the turbulence analysis.
In the text

	Fig. 2 Different snapshots of mass density in the orbital plane during a short period in the second full orbit as indicated in the plots. As given by the indicated orbital phases, these snapshots cover a period of slightly more than 40 min.
In the text

	Fig. 3 Absolute value of the spatial component of the fluid’s four velocity in the orbital plane at selected phases as indicated in the plots. Superimposed over the images, we indicate the direction of the flow via streamlines. Here, we only show results for a part of the orbital plane, focussing on the region around the Coriolis shock. We note the different x-axis in the right-handed plot.
In the text

	Fig. 4 Mach number in the orbital plane during the second full orbit, for orbital phases as indicated in the subplots. The arrows indicate the direction towards the observer.
In the text

	Fig. 5 Pressure (left) and mass density (right) in the region surrounding the unshocked pulsar wind for phase Φ = 1.310.
In the text

	Fig. 6 Distribution function of the relativistic gamma factor fγ within the orbital domain at different phases as indicated in the plots. Here, we show results for the second orbit (left), for the third orbit (middle), and also for the same short time interval as visualised in Fig. 2 (right).
In the text

	Fig. 7 Structure functions at four different orbital phases as indicated in the figures. For comparison, we show a linear function by the black dotted line. By our normalisation, the structure functions are given in units of cp.
In the text

	Fig. 8 Volume of unshocked pulsar wind Vpuisai as a function of orbital phase Φ. Data are shown for the full second and third orbit. The error bars indicate the corresponding standard deviation, that was computed from multiple output files at similar orbital phases, where we averaged over six phases spread over slightly more than 40 min.
In the text

Fig. 9 Distance of the Coriolis shock (blue plus signs) and bow shock (red crosses) from the location of the pulsar as functions of orbital phase, where distances are given along the line connecting star and pulsar. Average values and error bars were determined as stated in Fig. 8. The y-ахis uses logarithmic units to allow showing both shock distances simultaneously. Additionally, the dashed orange line shows the approximate distance of the contact discontinuity from the pulsar in the direction of the star using the formula by Eichler & Usov (1993).

In the text

	Fig. 10 Fraction of supersonic downstream medium as a function of orbital phase. Average values and error bars were determined as stated in Fig. 8.
In the text

	Fig. 11 Variation of the mean rate of inertial energy dissipation per unit mass ϵ as a function of orbital phase Φ over the second and the third orbit. As in Fig. 8, the error bars indicate the standard deviation computed from six consecutive output flies, ϵ is given in units of с3 AU−1.
In the text

	Fig. 12 Same as Fig. 1 but for the third full orbit.
In the text