Issue 
A&A
Volume 646, February 2021



Article Number  A91  
Number of page(s)  12  
Section  Astrophysical processes  
DOI  https://doi.org/10.1051/00046361/202039277  
Published online  12 February 2021 
Relativistic fluid modelling of gammaray binaries
I. The model
Institut für Astro und Teilchenphysik LeopoldFranzensUniversität Innsbruck, 6020 Innsbruck, Austria
email: david.huber@uibk.ac.at
Received:
26
August
2020
Accepted:
7
December
2020
Context. Gammaray binaries are systems that radiate the dominant part of their nonthermal emission in the gammaray band. In a winddriven scenario, these binaries are thought to consist of a pulsar orbiting a massive star, accelerating particles in the shock arising in the wind collision.
Aims. We develop a comprehensive numerical model for the nonthermal emission of shockaccelerated particles including the dynamical effects of fluid instabilities and orbital motion. We demonstrate the model on a generic binary system.
Methods. The model was built on a dedicated threedimensional particle transport simulation for the accelerated particles that were dynamically coupled to a simultaneous relativistic hydrodynamic simulation of the wind interaction. In a postprocessing step, a leptonic emission model involving synchrotron and inverseCompton emission was evaluated based on resulting particle distributions and fluid solutions, consistently accounting for relativistic boosting and γγabsorption in the stellar radiation field. The model was implemented as an extension to the CRONOS code.
Results. In the generic binary, the wind interaction leads to the formation of an extended, asymmetric windcollision region distorted by the effects of orbital motion, mixing, and turbulence. This gives rise to strong shocks terminating the pulsar wind and secondary shocks in the turbulent fluid flow. With our approach it is possible for the first time to consistently account for the dynamical shock structure in particle transport processes, which yields a complex distribution of accelerated particles. The predicted emission extends over a broad energy range, with significant orbital modulation in all bands.
Key words: radiation mechanisms: nonthermal / gamma rays: stars / methods: numerical / relativistic processes / hydrodynamics / binaries: general
© ESO 2021
1. Introduction
Gammaray binaries are composed of an earlytype massive star in orbit with a compact object, either a neutron star or a black hole. They are characterised by a dominant radiative output in the gammaray regime > 1 MeV (see Dubus 2013; Paredes & Bordas 2019, for a review). They exhibit broadband nonthermal emission ranging from radio through lowenergy (LE, 1–100 MeV) up to highenergy (HE, > 100 MeV) and very high energy (VHE, > 100 GeV) gammarays. For most systems, observations show flux variations with the orbital phase. At the time of writing, nine gammaray binaries emitting in the HE regime are confirmed: 1FGL J1018.65856 (Fermi LAT Collaboration 2012; H.E.S.S. Collaboration 2015a), 4FGL J14056119 (Corbet et al. 2019), HESS J0632+057 (Aharonian et al. 2007), HESS J1832093 (H.E.S.S. Collaboration 2015b; MartíDevesa & Reimer 2020), LMC P3 (Corbet et al. 2016; H.E.S.S. Collaboration 2018), LS 5039 (Paredes et al. 2000; Aharonian et al. 2005a), LS I +61 303 (Kniffen et al. 1997; Albert et al. 2006), PSR B125963 (Aharonian et al. 2005b), and PSR J2032+4127 (Abeysekara et al. 2018).
In the literature, two possible mechanisms are proposed to explain the nonthermal emission (see e.g. Mirabel 2006; Romero et al. 2007): A jetrelated emission scenario, where the compact object accretes matter from its stellar companion, releasing the accretion energy in the form of relativistic jets with highenergy particles (see e.g. BoschRamon & Khangulyan 2009); or a winddriven scenario, in which the compact object is commonly assumed to be a pulsar, whose rotational energy is dissipated as a relativistic pair plasma interacting with the wind of the earlytype companion star (see Maraschi & Treves 1981; Dubus 2006). In closeorbit binaries, this interaction leads to the formation of an extended windcollision region (WCR) that terminates both the pulsar and stellar wind by shocks (see e.g. Bogovalov et al. 2008; BoschRamon et al. 2012). At these sites, particles can be accelerated to ultrarelativistic energies through diffusive shock acceleration and other processes (Sironi & Spitkovsky 2011), which then produce the observed radiation. For many systems, a winddriven interpretation seems favoured. This has been investigated using a variety of different approaches ranging from simple fewzone models to more complex numerical models including different emission mechanisms (see e.g. Khangulyan et al. 2007; BoschRamon et al. 2008; Cerutti et al. 2010; Bednarek 2011; Zabalza et al. 2013; Takata et al. 2014; Dubus et al. 2015; Barkov & BoschRamon 2018; Molina & BoschRamon 2020). However, while the explored approaches succeed at modelling some observed features, they fail at reproducing others. In this work, we also focus on a winddriven scenario and aim to alleviate some of the encountered issues.
In analogy to gammaray binaries, extended WCRs can also be formed in collidingwind binaries (CWBs). These systems are very similar in terms of wind interaction, but they harbour a second earlytype star instead of the neutron star. CWBs can therefore be viewed as a nonrelativistic analogue of gammaray binaries in the winddriven scenario. The HE emission of such systems has been successfully described by Reitberger et al. (2017) employing a dynamically coupled treatment of the particle transport and the fluid. The corresponding numerical model was implemented as an extension to the CRONOS code (see Kissmann et al. 2018) solving a FokkerPlanck type particle transport equation for protons and electrons simultaneously to the equations of magnetohydrodynamics.
In this work, we extend this approach to gammaray binaries by translating the previous simulations into a relativistic framework, but we neglect the possible effect of the pulsar and stellar magnetic field on the fluid flow at this stage (Dubus et al. 2015; Bogovalov et al. 2019). This allows a production of emitting particle distributions that is consistent with the dynamics of the wind interaction; this has not been attempted in the case of gammaray binaries before. In turn, this enables the prediction of the nonthermal emission from gammaray binaries in a selfconsistent framework.
This is the first in a series of papers. With this first paper, we establish the governing equations, describe their numerical treatment, and demonstrate the approach on a generic binary system. In a second paper (Huber et al. 2020), the developed model is specifically applied to the LS 5039 system, and we compare the predicted emission to observations.
The paper is structured as follows: In Sects. 2.1 and 2.2 we introduce the governing equations for the coupled fluid and particle transport simulation, respectively, together with their numerical treatment and their implementation within CRONOS. In Sect. 2.3 we present the emission model involving synchrotron and inverseCompton emission with additional modulations by relativistic boosting and γγabsorption. In Sect. 3 we demonstrate the general properties of our model on a generic binary system and present the resulting fluid structure, particle distributions, and emission spectra and light curves. A concluding summary and outlook are given in Sect. 4.
2. Numerical model
The numerical model can be conceptually split into three distinct parts: the fluid model formulated in special relativity treating the interaction of the winds (see Sect. 2.1), the particle transport treating the evolution of the particle distributions embedded in the fluid (see Sect. 2.2), and the emission model predicting the radiation from the system (see Sect. 2.3). The first two are treated simultaneously in the simulation, whereas the radiation is computed in a postprocessing step using previously obtained particle and fluid solutions.
2.1. Fluid model
Owing to the relativistic pulsar wind, a classical treatment does no longer suffice and a special relativistic hydrodynamic framework has to be employed to accurately simulate the stellar and pulsar wind interaction. In this section we briefly summarise the necessary ingredients that are treated in an extension of the CRONOS code (Kissmann et al. 2018), which will be described in more detail in a future paper. The governing equations can be cast into a set of hyperbolic conservation equations that in Cartesian coordinates take the form
where U denotes the conserved quantities, F^{i} its respective flux, and S an additional source term. The first two can be expressed through primitive fluid variables such as the mass density ρ, the pressure p, the fourvelocity u^{μ}, and specific enthalpy h as
and
In this paper we use the normalisation c = 1 and notate the three spatial components of the fluid fourvelocity vector as u = u^{i} and its Lorentz factor as u^{0}. The presented system of five equations in six unknowns is closed by an additional equation of state. For now we employ an ideal equation of state,
where Γ = 4/3 is the constant adiabatic exponent of the fluid, as this was found to have only a negligible effect on the resulting fluid structure with respect to more sophisticate ones (BoschRamon et al. 2015).
In our simulation, we employed the HLLC solver (Mignone & Bodo 2005) alongside a secondorder subgrid reconstruction using a MINMOD slope limiter. To further increase the stability and robustness of our simulation, we directly evolved τ = E − D, which replaces the energy equation, and solved an additional conservation equation for the specific entropy s = p/ρ^{Γ} used in smooth flows,
2.1.1. Wind injection
The stellar and pulsar wind were injected in the simulation by prescribing a solution in a spherical region with radius r_{inj} of the domain. Inside each wind injection volume, we prescribed the solution as
where r denotes the vector originating at the injection centre to a given point in space, Ṁ the mass injection rate, and u the fourspeed of the injected wind. We kept the wind speed constant at its terminal velocity over the whole injection volume. The density gradient was scaled such that the prescribed massloss rate was recovered for every spherical shell inside the injection region. Because the effect of thermal pressure is neglected in our simulation, it was set to a low value scaled with the density gradient in order to keep the corresponding temperature constant.
In simulations of gammaray binaries, it is common practice to prescribe the pulsar massloss rate through the ratio η of the total momentum. The pulsar massloss rate can hence be obtained by
where quantities indicated with p and s belong to the pulsar and star, respectively.
In practice, the injection was realised by resetting the states of the respective cells after every timestep with those computed from the prescribed solution. A cell was treated as part of the star or pulsar, respectively, when its centre was within a spherical volume of given radius r_{inj} around the star or pulsar location.
2.1.2. Orbital motion
To account for the orbital motion, the wind injection volumes were relocated after every step of the simulation to their respective position on the Keplerian orbit. In standard approaches, the system is treated in the laboratory frame (see e.g. BoschRamon et al. 2012). In our approach, we treat the binary in a frame corotating with the average angular velocity of the system as described in Appendix A. This has the advantage that the vector connecting the two companions remains fixed for circular orbits instead of performing full rotations. For mildly eccentric orbits, this translates into an oscillation in a narrow angular range. Consequently, the wind injection volumes have to be relocated by a smaller amount than in the laboratory frame, which increases the stability of the simulation by avoiding unphysical conditions upon relocation.
Furthermore, because the binary motion is restricted, the orientation of the WCR is also restricted. This allows saving computational effort by reducing the simulation volume behind the star as seen from the centre of mass because we are mostly interested in the cometary taillike downstream of the WCR behind the pulsar.
2.2. Particle transport
The particle transport equation for cosmic rays has been presented in the relativistic case by Webb (1989). It takes the following form:
where is the isotropic part of the particle phase space density. All primed quantities are given in the fluid frame. The first term describes passive advection with the fluid bulk motion u^{μ} and spatial diffusion with the heat flux q^{μ}. The other terms in square brackets describe the transport in momentum space. From left to right, they correspond to adiabatic losses or gains, generic momentum losses (e.g. from radiative processes), viscous shear acceleration, the secondorder Fermi process, and losses due to noninertial frame changes. Following Vaidya et al. (2018), we neglected (for the purpose of this paper) spatial diffusion (q^{μ} = 0), viscous shearacceleration (Γ_{visc} = 0), and momentum diffusion (D_{pp} = 0). Because we assume the ultrarelativistic limit, the particle Lorentz factor can be readily obtained as γ ≈ p′/m, leading us to define , with 𝒩′ dV′ dγ′ corresponding to the number of particles within a volume dV′ in the range [γ′,γ′+ dγ′]. According to our assumptions and substitutions, we obtain the simplified transport equation
with being the total energy loss rate in the fluid frame. We now consider electrons and positrons in the particle transport because accelerated nuclei are thought to not contribute efficiently to the overall emission (e.g. BoschRamon & Khangulyan 2009).
2.2.1. Implementation
Following Reitberger et al. (2014), we implemented the transport equation Eq. (9) within CRONOS by employing a splitting scheme that separates the problem into a spatial and an energy part,
Both parts are treated independently and sequentially one after the other. The particle density was discretised in space in the same manner as the fluid variables, while its spectrum was subdivided into spectral bins . In our implementation, we allow the user to specify the limits of the simulated spectral range, that is, the first and the last edge , and the number of spectral bins N_{γ}. The range was then subdivided into logarithmic bins with
The cell centre is then obtained as the arithmetic mean of its limits.
Averaging Eq. (10) over the lth energy bin yields the finitevolume version of the spatial particle transport,
with the cell averages and . The spatial particle distribution corresponding to a given spectral bin can therefore be represented by an additional threedimensional scalar field . The spatial transport, amounting to a mere spatial convection, is solved directly by the hydrodynamic solver similar to the mass continuity equation (the first component in Eq. (1)).
Averaging Eq. (11) over a spatial cell with index i, j, k, and treating the fluid variables as constant in time yields
where is the spatial cell average of the differential particle number density at a certain energy. The spectral evolution of the particle distribution can therefore be viewed as an advection in energy with energydependent velocity .
The evolution in energy was treated independently for each spatial cell by a semiLagrangian solver in analogy to the one presented by Zerroukat et al. (2006), which we detail in the following. For better readability, we omit in the following the prime and the indices on the particle number density and the particle Lorentz factor γ ≡ γ′. We can show that the absolute number of particles in the energy range [γ_{−}(t),γ_{+}(t)] remains constant if the limits γ_{±}(t) satisfy the characteristic equation
which can be exploited to construct a conservative and explicit update scheme for the particle number densities.
We obtained the particle number density at the next timestep t^{n + 1} by first considering where the particles inside the bin were advected from, that is, in which γrange they were at time step t^{n}. This was done by performing a backward integration of the bin edges γ_{l ± 1/2} following Eq. (15) to identify the corresponding values at time t^{n} (denoting propagated edges with a superscript). Equation (15) was solved numerically using a secondorder explicit RungeKutta method, for example Ralston’s method, with the time step provided by the RHD solver.
After the interval edges were identified, the updated number density was obtained using particle conservation
where I^{n}[γ] denotes the cumulative particle count up to a value γ, for the spectrum at time t^{n} (see also Eq. (19)). To approximate the integral, we considered a secondorder, piecewise linear reconstruction in energy in each cell, in analogy to the spatial reconstruction in the hydrodynamics part,
with γ_{l − 1/2} < γ < γ_{l + 1/2}. Here, denotes the slope of the reconstruction obtained after applying the van Leer limiter (van Leer 1977),
to the right, left, and twosided finite differences, denoted by δ^{+}, δ^{−}, and δ^{C}, respectively. With this, the integral can be approximated directly as piecewise parabolas by
with γ_{l − 1/2} < γ < γ_{l + 1/2}, h_{γ} = γ − γ_{l − 1/2}, and the exact integrals at the cell edges,
where γ_{−1/2} is the lower edge of the first spectral bin l = 0.
On either end of the energy grid, we employed outflowing boundary conditions, ensuring that particles cannot be gained but only lost through the boundaries. This was realised by setting the particle number density to zero outside the energy domain, which amounts to a constant extrapolation of I^{n} in the evaluation of the righthand side of Eq. (16). To determine the slope of the first and last spectral bin according to Eq. (18), we added a ghostcell at either end of the spectrum (with l = −1 and l = N_{γ}, using Eq. (12)), and set the cell values to zero.
2.2.2. Magnetic field model
For the purpose of this work, we employed a simplified magnetic field model. We assumed the magnetic field to be amplified by plasma instabilities to a fraction ζ_{b} of the available internal energy. This fraction was kept constant throughout the simulation (see also Dubus et al. 2015; Barkov & BoschRamon 2018). The strength of the postshock magnetic field B′ in the fluid frame can therefore be expressed as
We assumed no direction on the magnetic field.
2.2.3. Radiative losses
As energyloss processes for the accelerated electronpositron pairs we considered synchrotron and inverseCompton losses. The energyloss rates are given by
with and being the magnetic and radiative energy densities in the fluid frame, respectively, and F_{KN} the KleinNishina factor.
The magnetic field in the fluid frame can be readily obtained by the transformation
The corresponding magnetic field strength is therefore given by
yielding the magnetic energy density . For disordered magnetic fields, that is, when no information on the magnetic field orientation is available, we dropped the term containing B ⋅ u in Eq. (25) (see Dubus et al. 2015).
For the inverseCompton losses, we assumed the radiation field to be monochromatic and monodirectional. The energy density of a beam of photons with direction Ω_{0} is given in the fluid frame by Doppler boosting the laboratory energy density,
with the Doppler factor 𝒟_{uΩ0} = (u^{0}−u⋅Ω_{0})^{−1}. The radiation field relevant for the inverseCompton losses emerges from the stellar binary companion,
where r_{⋆} is the distance to the star and L_{⋆} its luminosity. We employed the KleinNishina factor following Moderski et al. (2005),
assuming an isotropic seed photon and electron distribution for simplicity, with and denoting the Dopplerboosted photon energy of the stellar radiation field.
2.2.4. Particle injection
Strong shocks are commonly regarded as sites for particle acceleration, where we expect a fraction of pairs from the pulsar wind to be accelerated to high energies. However, because the converging fluid flow is the principal driver for several acceleration processes, we generally consider that strongly compressive flows participate in the acceleration. In the simulation, we identified these sites using the criterion
with a typical value of δ_{sh} = −10 c/AU to ignore weakly converging flows (see also Appendix C for a comparison to a common shockidentification criterion).
In the identified cells we injected two populations of pairs: thermal pairs from the cold pulsar wind that are isotropised, following a Maxwellian distribution; and nonthermal accelerated pairs, following a powerlaw spectrum (similar to Dubus et al. 2015), which are parametrised as
The maximum Lorentz factor for the power law is determined by the balance between acceleration and radiative losses. The acceleration timescale is given by a multiple of the gyration time of the particle τ_{acc} = 2πR_{L}/cξ_{acc} with the Larmor radius and the acceleration efficiency ξ_{acc}. We assumed synchrotron losses to be dominant for the particles at highest energies. Together with the loss timescale given by Eq. (22), this yields the maximum particle Lorentz factor expected at the current location,
Because the dominant particle acceleration mechanism in gammaray binaries could not be clearly identified, we took a generic phenomenological approach and treated both the spectral index s and the acceleration efficiency ξ_{acc} as free parameters of the model (similar to Dubus et al. 2015; Molina & BoschRamon 2020).
The average Lorentz factor of the Maxwellian γ_{t}, the lower cutoff of the powerlaw spectrum γ_{min}, and both normalisations K^{PLMW} were determined from fluid properties. This was controlled by three free parameters: the fraction of accelerated nonthermal pairs, the fraction of internal energy converted to nonthermal pairs, and ζ_{ρ} to account for the pulsar wind speed, which might be too low in our simulation and result in an overestimation of the particle number density in the pulsar wind. The respective fractions for the Maxwellian are then readily given by and . In our simulation, we kept all injection parameters constant in time and space. The injected spectra are therefore related to fluid quantities by
where 0 < ψ_{p} < 1 is a passive tracer indicating the fraction of pulsar wind material and the thermal energy density given by the fluid pressure.
For the nonthermal particles, the above expressions yield an implicit equation for , which was solved numerically by a standard rootfinding algorithm. The injection normalisation K^{PL} was then determined from Eq. (33). When the internal fluid energy was too low,
we did not inject particles in the given cell.
For the Maxwellian, Eqs. (33) and (34) can be inverted analytically in the limit γ_{t} ≫ 1 using and , yielding
2.3. Emission
The particle distributions resulting from a joint hydrodynamic and particletransport simulation (as detailed in Sects. 2.1 and 2.2) were used in a postprocessing step to predict the expected emission. In this section, we elaborate on the calculations of the emission and their implementation.
The involved radiative processes are again synchrotron (Sect. 2.3.1) and inverseCompton emission (Sect. 2.3.2). Furthermore, the produced fluxes were modified by γγabsorption on photons of the stellar radiation field depending on the location of the emitter (Sect. 2.3.3). For the sake of better readability, we omit the dependence on the location x when it is not imperative.
We formulated the emission in the fluid frame, represented by primed quantities, and then related it to the laboratory frame with a Doppler boost. The spectral emissivity j_{ϵ}, that is, the emitted power per unit photon energy per unit solid angle per unit volume, for a given emitted photon energy ϵ and direction Ω is hence given by
where 𝒟_{uΩ} = (u^{0}−Ω⋅u)^{−1} is the corresponding Doppler factor. The respective energy and direction of the emitted photons are given in the fluid frame by
Furthermore, it is convenient to express the spectral emissivity in terms of the emission produced by a single electron 𝒫′(ϵ′,Ω′,γ′), that is, the emitted power per unit photon energy per unit solid angle per electron under the assumption of an isotropic distribution of particles in the fluid frame. The spectral emissivity can then be formulated as
2.3.1. Synchrotron emission
The total power emitted by a single electron in a magnetic field B′ can be expressed as
which is peaked around the critical photon energy,
denotes the magnetic field component perpendicular to the electron motion, which gyrates with a pitch angle α around the magnetic field lines. For ultrarelativistic electrons, the emission is strongly beamed in their instantaneous direction of motion (within a cone of angle 1/γ′). Only particles whose pitch angle coincides with the angle between the magnetic field B′ and the line of sight Ω′ therefore contribute to the emission, effectively yielding = B′×Ω′. Together with the assumption of an isotropic particle distribution, this yields
In many cases, we approximated the magnetic field to be disordered, that is, we have no information on the direction of the field lines. In this case, we approximated the in Eqs. (41), (42) by its pitchangle average .
2.3.2. InverseCompton emission
Through the inverseCompton process, photons of the stellar radiation field are scattered to gammaray energies by highly energetic particles. Following the description of Moderski et al. (2005), the power emitted through this process is given by
where is the scattering rate of photons on isotropically distributed electrons. The scattering rate was presented by Aharonian & Atoyan (1981) in the case of and γ′≫1 taking the form
where n_{0}(ϵ_{0}, Ω_{0}) is the density of the seed photons with energy ϵ_{0} and direction Ω_{0} and the lowest seed photon energy that is still scattered to a given energy ϵ′ and direction Ω′
The scattering kernel f_{ic} is given by
with , w′=ϵ′(γ′m_{e}c^{2})^{−1} and the cosine of the scattering angle , which can be expressed in terms of laboratory quantities by
Together with the relativistic invariant Eq. (44) can be expressed directly in terms of the seed photon field in the laboratory frame,
Applying the factorisation of the target photon field (Appendix B) yields
which can be further simplified depending on the employed radiation field model.
2.3.3. γγabsorption
As a result of the strong stellar photon field, parts of the VHE gammaray flux produced through inverseCompton scattering are significantly attenuated by γγpair creation. This effect depends strongly on the geometry of the system and the location of the observer, yielding strongly modulated spectra with the binary orbit. For a consistent model of gammaray binary emission, this attenuation has to be taken into account for every cell in our simulation.
The absorption opacity τ_{γγ} of photons with energy ϵ originating at x in direction Ω can be expressed by a threefold integration of the differential absorption opacity given by Gould & Schréder (1967),
over the line of sight l, the solid angle Ω_{0}, and the seed photon energy ϵ_{0}. Here, μ = Ω ⋅ Ω_{0} denotes the scattering angle cosine between a VHE photon with direction Ω and a seed UV photon with direction Ω_{0}, the speed of the created e^{±} pair normalised to c and the scattering crosssection with
The lower bound of the energy integration is provided by the lowest seed photon energy to create an e^{±} pair under a given scattering angle. The stellar seed photon density n_{ϵ0} is defined in the same way as in Sect. 2.3.2.
For numerical purposes, it is more suitable to transform the occurring improper integrals into proper ones. For the lineofsight integration, we therefore used the transformation (as detailed in Dubus 2006). The integration over target photon energy ϵ_{0} was transformed into an integral over the speed of the created pairs β. Together with the radiation field factorisation (Appendix B), this yields
which can be further simplified depending on the stellar radiation field model.
The integration over solid angle Ω_{0} is approximated by a nested quadrature rule with equidistant points in azimuth and a GaussLegendre quadrature in the polar angle (see Hesse & Womersley 2012). For the integrals over ψ and β, we employed standard GaussLegendre quadratures. The necessary routines are provided by the GSL.
3. Application to a generic binary
In the following section, we apply the presented numerical model to a generic binary system. We focus on exploring the general properties of our model.
The considered binary is composed of an Otype star in a circular (e = 0) 1 d orbit with a pulsar at an orbital separation of d = 0.2 AU. We adopted a stellar massloss rate of Ṁ_{s} = 2 × 10^{−8} M_{⊙} yr^{−1} launched at a velocity of v_{s} = 2000 km s^{−1}. The pulsar wind was injected at a speed of v_{p} = 0.99 c and a massloss rate scaled with η = 0.1 in Eq. (7).
The simulation was performed on a Cartesian grid in corotating coordinates (as detailed in Sect. 2.1.2). The computational volume has the dimensions [ − 0.4, 0.4] × [−0.2, 0.6] × [−0.2, 0.2] AU^{3} and is homogeneously subdivided into 256 × 256 × 128 spatial cells. The binary orbits in the x − y plane, with the stellar companion residing at the coordinate origin.
The stellar and pulsar winds are injected as described in Sect. 2.1.1 with an injection radius of four cell widths. Initially, both injection volumes were placed in a vacuum. We ran the simulation for 0.58 d giving the slow stellar wind time to populate the computational domain.
3.1. Fluid structure
In Fig. 1 we show the resulting fluid mass, bulk fourspeed, and magnetic field strength in the fluid frame, using the magnetic field model as described in Sect. 2.2.2 with ζ_{b} = 0.5. The interaction of the stellar and pulsar wind gives rise to an extended WCR delimited by shocks on both the stellar and the pulsar side. Through the Coriolis force in the rotating system, the WCR is highly asymmetric and bends in the direction opposite to the orbital motion. This leads to the termination of the pulsar wind on the far side of the system, forming another shock at the Coriolis turnover (hereafter Coriolis shock). The location of the Coriolis shock was noted previously and is consistent with previous studies (see e.g. BoschRamon et al. 2015; Dubus et al. 2015). In the wings of the WCR, the shocked pulsar wind is reaccelerated to a Lorentz factor ∼3, again in agreement with previous simulations by Bogovalov et al. (2008).
Fig. 1. Resulting fluid structure for a generic binary system in the orbital plane. From top to bottom: fluid mass density, bulk speed, and magnetic field strength in the fluid frame. An extended WCR is formed by interaction of the pulsar (blue) and stellar (orange) wind. The reference frame is corotating with counterclockwise orbital motion. The location of the most relevant shocks is annotated. 
The strong velocity shear between the shocked winds at the contact discontinuity and the steep density gradient at the stellar wind shock leads to the development of fluid instabilities (see also BoschRamon et al. 2012). As a consequence, many secondary shocks are formed in the turbulent regions.
3.2. Particle distribution
To economise on computation time, we did not perform a particle transport simulation for the full time span used in Sect. 3.1. Instead, we restarted such a simulation on a preconverged fluid solution and ran it for a shorter time span t = 0.029 d, allowing the accelerated electronpositron pairs to populate the system, which we also refer to as electrons for brevity.
The simulation presented in this section was carried out using 50 logarithmic bins in energy ranging from γ ∈ [200, 4 × 10^{8}]. The particle injection was controlled as detailed in Sect. 2.2.4 using ζ_{ρ} = 5 × 10^{−4}, , , s = 2, and ξ_{acc} = 1, assuming the maximum acceleration timescale at the Bohm limit.
In Fig. 2 we present the resulting electron distributions at two energies. The shock structure leaves a clear imprint on the particle distribution. This is especially visible for higher energetic electrons, which remain much more confined to their acceleration site due to increased energyloss rates, as compared to lower energetic electrons. Next to the bow and Coriolis shocks, the secondary shocks arising in the turbulent flow are clearly visible, which are for the first time consistently accounted for in the particle transport.
Fig. 2. Differential particle number density for the energy bin centred at γ = 4.2 × 10^{3} (upper panel) and γ = 1.1 × 10^{7} (lower panel) in analogy to Fig. 1. The lower energetic distribution is mainly populated by Maxwellian electrons, whereas the higher energetic distribution is dominated by powerlaw electrons. For clarity we show only five orders of magnitude below the highest value. The blue regions correspond to those with lower values. 
In Fig. 3 we show the spectral energy distribution of electrons integrated over the computational volume. The contributions of the Maxwellian and the powerlaw electrons can be easily distinguished, with Maxwellian electrons dominating the lowenergy part (γ ≲ 2 × 10^{4}) and the nonthermal electrons the highenergy part. For the chosen set of parameters, Maxwellian electrons were on average injected with γ_{t} ∼ 1.0 × 10^{3} and nonthermal electrons within the limits of γ_{min} ∼ 1.1 × 10^{5} and γ_{max} = 2.4 × 10^{7}.
Fig. 3. Spectral energy distribution of accelerated electrons integrated over the computational domain for different orbital phases. 
The number and extent of particle acceleration sites identified in the turbulent Coriolis shock downstream region were determined by the chosen compression threshold in Eq. (29) (see also Appendix C). The particles accelerated in this region reach higher energies because the magnetic field is lower than in regions near the bow shock. Choosing a more restrictive threshold, that is, one with higher compression, would thus decrease the total number of particles at higher energies in the simulation.
3.3. Emission
In this section, we present the predicted emission for the generic binary. The involved computations were performed in a postprocessing step as detailed in Sect. 2.3, using the previously obtained particle distribution (see Sect. 3.2) and fluid solution (see Sect. 3.1).
We chose the observer to be located at a distance of d_{obs} = 2.5 kpc along the +y axis of the system for the orbital phase ϕ = 0. The phases ϕ = 0 and ϕ = 0.5 therefore correspond to inferior and superior conjunction, respectively. As a result of the employed corotating coordinates, the observer rotates clockwise with respect to the simulation frame. For the present purpose, the underlying particle and fluid solutions were kept constant for all orbital phases.
To compute the inverseCompton emission, we treated the stellar photon field as monochromatic, whereas for the γγ absorption, a full blackbody spectrum was considered. In both cases, the source was treated as an extended sphere with luminosity L_{⋆} = 2 × 10^{5} L_{⊙}, temperature T_{⋆} = 40 000 K, and resulting radius R_{⋆} = 9.3 R_{⊙}. The emission was computed for a system inclination of i = 30° and 20 equidistant orbital phases.
In Fig. 4 we show the resulting spectral energy distribution of emitted photons. Ranging from Xrays up to LE gammarays, the spectrum is dominated by synchrotron emission from electrons at the powerlaw tail of the spectrum. The same population of electrons is responsible for the inverseCompton scattering of stellar UV photons to the VHE gammaray regime, which are heavily attenuated by γγabsorption in the 0.1 − 10 TeV range, depending on the orbital phase. Again, the emission in the HE gammaray regime is produced through the inverseCompton process, but by the lowenergy Maxwellian electrons.
Fig. 4. Spectral energy distribution of the emission predicted by our model for a generic binary with a system inclination of 30°. The results are shown for the sampled orbital phases (grey) and their average (black). The average spectral distribution is further split into the individual radiative processes: synchrotron (green) and inverse Compton attenuated by γγ absorption (orange). Inferior (dashed red, ϕ = 0) and superior conjunction (dashed blue, ϕ = 0.5) emissions are highlighted. 
In our model, the predicted spectra can be tuned most directly by varying the parameters that regulate the emerging electron distributions, that is, the unknown injection parameters and the magnetic field model. Here, the acceleration efficiency ξ_{acc} limits the maximum energy reached γ_{max} (see Eq. (32)), hence affecting the location of the LE and VHE cutoff. The shape of the injected electron spectrum, and consequently, the main emission features, are determined by the electron specific energy density given by the ratio . In the case of the Maxwellian electrons, this can be used to shift the inverseCompton bump at HE in energy. For the nonthermal powerlaw electrons, their specific energy density affects the lowenergy cutoff of the spectrum. The overall normalisation of the spectrum is given by the amount of energy deposited in the respective populations, determined by . The index of the injected electrons translates into the spectral index of the Xray synchrotron emission and affects the spectral index at VHE gammarays. ζ_{b} is used to control the magnetic field strength in the simulation that affects the amplitude of the synchrotron emission, the VHE cutoff, and the evolution of the electron distribution itself with more subtle consequences for the emission spectra. Because the compression threshold δ_{sh} mostly affects the particle density at highest energies in the turbulent Coriolis shock downstream, this parameter can affect the synchrotron fluxes produced in the LE gammaray band, depending on the magnetic field in this region, and the VHE gammaray flux at superior conjunction.
In Fig. 5 we show the predicted light curves for different energy bands, which exhibit a significant orbital modulation in all bands. The Xray and LE gammaray light curves show a clear correlation, which is expected because both bands are dominated by synchrotron emission generated by the same particle population. Because we assumed no direction for the magnetic field, the synchrotron emissivity is isotropic in the fluid frame of our model. The driver for the orbital modulation in these bands is accordingly solely given by relativistic boosting. Because the WCR is bent by orbital motion, relativistic boosting is at its maximum at ϕ ∼ 0.1, while deboosting is most significant for ϕ ∼ 0.6.
Fig. 5. Emission light curves predicted by our model for a generic binary for different energy bands and a system inclination of 30°. Inferior and superior conjunction correspond to ϕ = 0 and ϕ = 0.5, respectively. For better visualisation, we show the data for two full orbits. 
At HE gammarays, the anisotropic inverseCompton scattering crosssection introduces a further lineofsightdependent modulation. The scattering process is most efficient for large scattering angles, which are achieved at superior conjunction (ϕ = 0.5). This is reflected in the respective light curve, which reaches its maximum shortly before superior conjunction through the onset of relativistic deboosting, yielding an anticorrelation of the HE gammaray band and the previously discussed one.
The VHE gammaray emission is again produced by inverseCompton process, but the upscattered photons are attenuated by γγpair creation with the stellar photon field. This once more introduces a lineofsightdependent modulation that reaches its maximum attenuation at superior conjunction, where the produced VHE photons have to propagate through the entire stellar radiation field. Consequently, the maximum photon flux can escape the system at inferior conjunction, with a slight shift to ϕ = 0.1 due to relativistic boosting, leaving the VHE correlated with the Xray to LE gammaray band and therefore anticorrelated to HE gammarays.
4. Summary and outlook
We presented a novel numerical model for gammaray binaries treating the transport of shockaccelerated electronpositron pairs dynamically coupled to the relativistic fluid dynamics of the stellar and pulsar wind interaction. With this, it is possible for the first time to obtain selfconsistent threedimensional particle distributions and fluid solutions while consistently accounting for dynamic fluid instabilities, mixing, and the effects of orbital motion. The proposed emission model is purely leptonic, thus we inferred the expected emission through the synchrotron and inverseCompton process acting on accelerated pairs from the pulsar wind obtained through the particle transport simulation. Together with the simultaneously obtained fluid solution, it is therefore possible to consistently account for relativistic boosting and γγabsorption. The model was implemented as a relativistic extension to the CRONOS code (Kissmann et al. 2018).
The capabilities of the presented model were demonstrated on a generic binary system. The wind interaction yields an extended asymmetric wind collision region bent by the effects of orbital motion, exhibiting a strong bowlike pulsarwind shock and a Coriolis shock behind the pulsar. The expressed features are in qualitative agreement with previous works (see e.g. BoschRamon et al. 2015).
The developing fluid instabilities give rise to the formation of secondary shocks. Our model naturally includes the reacceleration of cooled particles at these shocks, which has not been done before in the case of gammaray binaries. In our generic simulation, parts of the turbulence are damped by numerical diffusion due to our choices of spatial resolution. More highly resolved simulations might hence yield a more turbulent fluid and consequently more secondary shocks.
The emission predicted for the generic binary was extended over a broad range in energy reaching up to VHE gammarays. Our results are consistent with the simulations performed by Dubus et al. (2015), showing a strong resemblance because the investigated models are similar.
Our model predicts significant orbital modulation in every energy band. The main driver for the modulation in the Xray to LE gammaray synchrotron emission is given by relativistic boosting in the wings of the WCR. In addition to the modulation introduced by the anisotropic inverseCompton scattering producing the HE and VHE gammaray emission, the VHE band is further modulated by γγabsorption. For the investigated general binary, the model naturally predicts a correlated Xray to LE gammaray and VHE gammaray emission, while being anticorrelated to the HE band, as is seen in the observations of many gammaray binaries (e.g. LS 5039, Chang et al. 2016).
Because many of the observed gammaray binaries are known to host Be stars, a possible extension of the model foresees the inclusion of more complex, anisotropic stellar wind structures, such as adding an equatorial disc. In the forthcoming second paper of this series, we will apply the presented model specifically to the wellstudied LS 5039 system (Huber et al. 2020), for which no stellar disc has been found. The wealth of available data means that the comparison of our model predictions with observations will provide valuable insights into gammaray binaries and a calibration of our model, pointing out potential aspects that might need more modelling effort in the future.
This might, for example, consist of an inclusion of the pulsar and stellar magnetic fields in the dynamics of the wind interaction by extending the model to relativistic magnetohydrodynamics. This might be used to obtain a selfconsistent picture for the magnetic field strength and direction, which could further enable the incorporation of more sophisticated particle acceleration models in the future, for example including the obliquity of magnetised shocks (e.g. Vaidya et al. 2018, and references therein).
Another possible extension of this model could consist of easing some simplifications made in the particle transport model, for example by including spatial diffusion. Next to the obvious changes in the spatial particle distribution, an additional effect is provided by the loss term accounting for noninertial frame changes. This effect has not been studied in the context of gammaray binaries.
In future, our model can further help to shed light on dynamical phenomena that can now be addressed with a higher level of sophistication. This could involve investigations of the role of turbulence and the variability it introduces in the radiative output of gammaray binaries. Particularly, it remains to be investigated in which way the emission from these systems may vary from one orbit to the next. In a broader context, this might also yield implications for and connect to transient phenomena, such as flaring events in otherwise regular gammaray binary emissions (e.g. PSR B125963 Johnson et al. 2018).
Acknowledgments
The computational results presented in this paper have been achieved (in parts) using the research infrastructure of the Institute for Astro and Particle Physics at the University of Innsbruck, the LEO HPC infrastructure of the University of Innsbruck, the MACH2 Interuniversity Shared Memory Supercomputer and PRACE resources. We acknowledge PRACE for awarding us access to JoliotCurie at GENCI@CEA, France. This research made use of Cronos (Kissmann et al. 2018); GNU Scientific Library (GSL) (Galassi 2018); matplotlib, a Python library for publication quality graphics (Hunter 2007); and NumPy (Van Der Walt et al. 2011). We would like to thank the unknown referee/reviewer for the thoughtful comments and efforts towards improving our manuscript.
References
 Abeysekara, A. U., Benbow, W., Bird, R., et al. 2018, ApJ, 867, L19 [Google Scholar]
 Aharonian, F. A., & Atoyan, A. M. 1981, Ap&SS, 79, 321 [NASA ADS] [Google Scholar]
 Aharonian, F., Akhperjanian, A. G., Aye, K. M., et al. 2005a, Science, 309, 746 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Aharonian, F., Akhperjanian, A. G., Aye, K. M., et al. 2005b, A&A, 442, 1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Aharonian, F. A., Akhperjanian, A. G., BazerBachi, A. R., et al. 2007, A&A, 469, L1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Albert, J., Aliu, E., Anderhub, H., et al. 2006, Science, 312, 1771 [Google Scholar]
 Banyuls, F., Font, J. A., Ibáñez, J. M., Martí, J. M., & Miralles, J. A. 1997, ApJ, 476, 221 [Google Scholar]
 Barkov, M. V., & BoschRamon, V. 2018, MNRAS, 479, 1320 [Google Scholar]
 Bednarek, W. 2011, MNRAS, 418, L49 [NASA ADS] [Google Scholar]
 Bogovalov, S. V., Khangulyan, D. V., Koldoba, A. V., Ustyugova, G. V., & Aharonian, F. A. 2008, MNRAS, 387, 63 [Google Scholar]
 Bogovalov, S. V., Khangulyan, D., Koldoba, A., Ustyugova, G. V., & Aharonian, F. 2019, MNRAS, 490, 3601 [Google Scholar]
 BoschRamon, V., & Khangulyan, D. 2009, Int. J. Mod. Phys. D, 18, 347 [Google Scholar]
 BoschRamon, V., Khangulyan, D., & Aharonian, F. A. 2008, A&A, 482, 397 [Google Scholar]
 BoschRamon, V., Barkov, M. V., Khangulyan, D., & Perucho, M. 2012, A&A, 544, A59 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 BoschRamon, V., Barkov, M. V., & Perucho, M. 2015, A&A, 577, A89 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Cerutti, B., Malzac, J., Dubus, G., & Henri, G. 2010, A&A, 519, A81 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Chang, Z., Zhang, S., Ji, L., et al. 2016, MNRAS, 463, 495 [Google Scholar]
 Corbet, R. H. D., Chomiuk, L., Coe, M. J., et al. 2016, ApJ, 829, 105 [Google Scholar]
 Corbet, R. H. D., Chomiuk, L., Coe, M. J., et al. 2019, ApJ, 884, 93 [Google Scholar]
 Dubus, G. 2006, A&A, 451, 9 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Dubus, G. 2013, A&ARv, 21, 64 [Google Scholar]
 Dubus, G., Lamberts, A., & Fromang, S. 2015, A&A, 581, A27 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Fermi LAT Collaboration (Ackermann, M., et al.) 2012, Science, 335, 189 [Google Scholar]
 Galassi, M. 2018, GNU Scientific Library Reference Manual [Google Scholar]
 Gould, R. J., & Schréder, G. P. 1967, Phys. Rev., 155, 1404 [Google Scholar]
 H.E.S.S. Collaboration (Abramowski, A., et al.) 2015a, A&A, 577, A131 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 H.E.S.S. Collaboration (Abramowski, A., et al.) 2015b, MNRAS, 446, 1163 [NASA ADS] [Google Scholar]
 H.E.S.S. Collaboration (Abdalla, H., et al.) 2018, A&A, 610, L17 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Hesse, K., & Womersley, R. S. 2012, Adv. Comput. Math., 36, 451 [Google Scholar]
 Huber, D., Kissmann, R., Reimer, A., & Reimer, O. 2020, A&A, submitted [arXiv:2012.04975] [Google Scholar]
 Hunter, J. D. 2007, Comput. Sci. Eng., 9, 90 [Google Scholar]
 Johnson, T. J., Wood, K. S., Kerr, M., et al. 2018, ApJ, 863, 27 [Google Scholar]
 Khangulyan, D., Hnatic, S., Aharonian, F., & Bogovalov, S. 2007, MNRAS, 380, 320 [Google Scholar]
 Kissmann, R., Kleimann, J., Krebl, B., & Wiengarten, T. 2018, ApJS, 236, 53 [Google Scholar]
 Kniffen, D. A., Alberts, W. C. K., Bertsch, D. L., et al. 1997, ApJ, 486, 126 [Google Scholar]
 Maraschi, L., & Treves, A. 1981, MNRAS, 194, 1P [Google Scholar]
 MartíDevesa, G., & Reimer, O. 2020, A&A, 637, A23 [EDP Sciences] [Google Scholar]
 Mignone, A., & Bodo, G. 2005, MNRAS, 364, 126 [Google Scholar]
 Mignone, A., & McKinney, J. C. 2007, MNRAS, 378, 1118 [Google Scholar]
 Mirabel, I. F. 2006, Science, 312, 1759 [Google Scholar]
 Misner, C. W., Thorne, K. S., & Wheeler, J. A. 1973, Gravitation (San Francisco: W. H. Freeman) [Google Scholar]
 Moderski, R., Sikora, M., Coppi, P. S., & Aharonian, F. 2005, MNRAS, 363, 954 [Google Scholar]
 Molina, E., & BoschRamon, V. 2020, A&A, 641, A84 [CrossRef] [EDP Sciences] [Google Scholar]
 Papadopoulos, P., & Font, J. A. 2000, Phys. Rev. D, 61, 024015 [Google Scholar]
 Paredes, J. M., & Bordas, P. 2019, ArXiv eprints [arXiv:1901.03624] [Google Scholar]
 Paredes, J. M., Martí, J., Ribó, M., & Massi, M. 2000, Science, 288, 2340 [Google Scholar]
 Pons, J. A., Font, J. A., Ibanez, J. M., Marti, J. M., & Miralles, J. A. 1998, A&A, 339, 638 [NASA ADS] [Google Scholar]
 Reitberger, K., Kissmann, R., Reimer, A., Reimer, O., & Dubus, G. 2014, ApJ, 782, 96 [Google Scholar]
 Reitberger, K., Kissmann, R., Reimer, A., & Reimer, O. 2017, ApJ, 847, 40 [Google Scholar]
 Romero, G. E., Okazaki, A. T., Orellana, M., & Owocki, S. P. 2007, A&A, 474, 15 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Sironi, L., & Spitkovsky, A. 2011, ApJ, 741, 39 [Google Scholar]
 Takata, J., Leung, G. C. K., Tam, P. H. T., et al. 2014, ApJ, 790, 18 [Google Scholar]
 Vaidya, B., Mignone, A., Bodo, G., Rossi, P., & Massaglia, S. 2018, ApJ, 865, 144 [Google Scholar]
 Van Der Walt, S., Colbert, S. C., & Varoquaux, G. 2011, Comput. Sci. Eng., 13, 22 [Google Scholar]
 van Leer, B. 1977, J. Comput. Phys., 23, 276 [Google Scholar]
 Webb, G. M. 1989, ApJ, 340, 1112 [Google Scholar]
 Zabalza, V., BoschRamon, V., Aharonian, F., & Khangulyan, D. 2013, A&A, 551, A17 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Zerroukat, M., Wood, N., & Staniforth, A. 2006, Int. J. Numer. Methods Fluids, 51, 1297 [Google Scholar]
Appendix A: Corotating frame
We considered a reference frame that corotates with the mean angular velocity of the binary Ω in the x − y plane. The corotation velocity is in this case given by the orbital period of the system as Ω = 2π/P_{orb}. The set of new coordinates, indicated by a bar, can be readily expressed as
Transforming the flat Cartesian Minkowski metric into the corotating system yields
Generic metric tensors can be accounted for in the broader framework of general relativistic hydrodynamics. Here, the equations for energy and momentum conservation are given by ∇_{μ}T^{μν} = 0 with the energymomentum tensor T^{μν} = ρhu^{μ}u^{ν} + pg^{μν}. Because of the nonvanishing, offdiagonal terms in the metric, the SRHD schemes can no longer be applied without modification.
In a 3+1 decomposition (Misner et al. 1973), the line element is expressed as ds^{2} = −(α^{2} − β_{i}β^{i}) dt^{2} + 2β_{i} dx^{i} dt + γ_{ij} dx^{i} dx^{j}, with the lapse, the shift, and the spatial metric. These gauge functions can be directly obtained from the metric tensor Eqs. (A.5) and (A.6). For the presented corotating system, they are given by
In our approach, we relied on the socalled 3+1 Valencia formulation of general relativistic hydrodynamics (Banyuls et al. 1997), and chose the conserved variables
consisting of the restmass density D, the momentum density in the jdirection m_{j}, and the total energy density E, measured by a family of Eulerian observers. In their respective frames, the fluid threevelocity is given by w^{i}, which is related to the fourvelocity in the rotating system by . The corresponding Lorentz factor is given by W = (1−w_{i}w^{i})^{−1/2} with w_{i} = γ_{ij}w^{j}.
The chosen variables closely resemble those in SRHD Eq. (2) up to the usage of covariant momentum densities instead of contravariant ones, which are equivalent and interchangeable in our case because . This allows an effective reuse of the variable inversion scheme developed for the special relativistic case that is already implemented in CRONOS.
The evolution equations take the generic form
with γ = det(γ_{ij}) and the fluxes
The system of equations again closely resemble those in special relativity Eq. (1) up to the additional fluxes introduced by β^{i} ≠ 0. We treated the equations using the scheme presented by Pons et al. (1998), describing how special relativistic Riemann solvers can be used with general metrics. The required transformations between the Eulerian and the local inertial frames in our case take the especially simple form
The source terms in Eq. (A.11) are given by (Banyuls et al. 1997)
and reduce to
in our case. This result may seem counterintuitive at first glance because the terms formally do not correspond to their Newtonian analogue for the centrifugal and Coriolis force in a rotating system. However, we recall that the conserved quantities are measured by the Eulerian observers instead of in the rotating frame.
Finally, our approach is equivalent to others such as Papadopoulos & Font (2000), who chose T^{0ν} as conserved variables and evolved the quantities directly in the rotating frame, which recovers the wellknown terms for the fictitious forces. This formulation, however, is less efficient in practice because the involved fluxes are more complex, and the whole metric has to be considered for contractions, for example in the computation of Lorentz factors, instead of just the spatial metric in the 3+1 Valencia formulation, which is unity in our case. Furthermore, the chosen conserved variables formally differ from Eq. (2), prohibiting the reuse of already developed SRHD variable inversion schemes (e.g. Mignone & McKinney 2007).
Ultimately, we are interested in quantities measured in the stationary laboratory frame to compute the radiative output of the system. This amounts to a Lorentz transformation of the corotating velocity field , which is obtained from the velocities measured by the Eulerian observers to the stationary u^{μ}. For our case, we find after some algebra that the Eulerian velocities directly correspond to those in the stationary system up to a rotation in space. This rotation is accounted for in the computation of relativistic boosting described in Sect. 2.3.
Appendix B: Radiation fields
Depending on the radiation field at hand, the integrals in the computation of inverseCompton emission or γγ absorption can be simplified. For example, for monochromatic seed photon fields, the energy integration reduces to a trivial substitution. To be able to map these properties over to our implementation, we employed a factorisation of the radiation field. In general, a radiation field can be described as the number of photons n per unit volume, unit energy, and unit solid angle, depending on the location x, the photon energy ϵ, and the direction Ω. For our purposes, however, the spectral shape is direction independent, which allows the following factorisation of the photon density
with a factor g_{n} describing the total photon density, a factor g_{ϵ} describing the energy dependence, and a factor g_{Ω} describing the directional profile,
We note that g_{ϵ} and g_{Ω} are normalised with respect to ϵ and Ω, respectively.
When we consider spherically expanding photons injected at a location x_{⋆} with a given rate Ṅ_{⋆}, the photon number density is given by
where r = ‖x−x_{⋆}‖.
A monodirectional or beamed radiation field in a given direction Ω_{beam} is in general described by
The radiation field of a pointlike source can therefore be viewed as a special case of a field beamed in the radial direction Ω_{beam} = e_{r} = r/r. An extended spherical source is approximated as a disc with constant brightness and angular extent limited by the source radius R,
where μ_{0} is the cosine of the angle enclosed between Ω and e_{r} and μ_{0, min} = (1−R^{2}/r^{2})^{1/2}.
The radiation field of the companion star can be approximated as a blackbody spectrum with temperature T,
This can be further approximated as a monochromatic spectrum with the average photon energy ,
In our simulation, the radiation field of the star was parametrised by its luminosity L_{⋆} and its temperature T_{⋆}. The extent of the source is therefore given by and its photon injection rate by .
Appendix C: Comparison to a common shockidentification criterion
In the literature, particle acceleration is thought to occur at strong shocks, which are usually identified by a diverging fluid flow and a jump in pressure (see e.g. Vaidya et al. 2018). As described in Sect. 2.2.4, we employed a less restrictive criterion on the fluid compression alone because a jump in pressure might not yet have developed for a given instance in time.
For the generic binary system considered in Sect. 3, we found both sets of criteria to identify the same main shock features, that is, the bow and Coriolis shock, with reflected shocks at their encounter and secondary shocks. In comparison to our criterion, the usual approach identifies fewer and less extended acceleration sites in the turbulent Coriolis shock downstream region because of the additional pressure constraint. This introduces a slight decrease in the emission at all wavelengths through the reduced particle injection, which could be compensated for by a different set of injection parameters. Because the magnetic field in the Coriolis shock downstream region is lower, the effect becomes more relevant for particles at higher energies, leading to a decrease in the synchrotron emission in the LE gammaray band. The temporal behaviour remains unaffected in all bands. For a compression threshold of δ_{sh} ∼ −50 c/AU, we found that our criterion identified the same acceleration sites as the classical approach. We treated the threshold as a free parameter of the simulation, which might be constrained for a given system compared with observations.
All Figures
Fig. 1. Resulting fluid structure for a generic binary system in the orbital plane. From top to bottom: fluid mass density, bulk speed, and magnetic field strength in the fluid frame. An extended WCR is formed by interaction of the pulsar (blue) and stellar (orange) wind. The reference frame is corotating with counterclockwise orbital motion. The location of the most relevant shocks is annotated. 

In the text 
Fig. 2. Differential particle number density for the energy bin centred at γ = 4.2 × 10^{3} (upper panel) and γ = 1.1 × 10^{7} (lower panel) in analogy to Fig. 1. The lower energetic distribution is mainly populated by Maxwellian electrons, whereas the higher energetic distribution is dominated by powerlaw electrons. For clarity we show only five orders of magnitude below the highest value. The blue regions correspond to those with lower values. 

In the text 
Fig. 3. Spectral energy distribution of accelerated electrons integrated over the computational domain for different orbital phases. 

In the text 
Fig. 4. Spectral energy distribution of the emission predicted by our model for a generic binary with a system inclination of 30°. The results are shown for the sampled orbital phases (grey) and their average (black). The average spectral distribution is further split into the individual radiative processes: synchrotron (green) and inverse Compton attenuated by γγ absorption (orange). Inferior (dashed red, ϕ = 0) and superior conjunction (dashed blue, ϕ = 0.5) emissions are highlighted. 

In the text 
Fig. 5. Emission light curves predicted by our model for a generic binary for different energy bands and a system inclination of 30°. Inferior and superior conjunction correspond to ϕ = 0 and ϕ = 0.5, respectively. For better visualisation, we show the data for two full orbits. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.