© ESO 2021

1. Introduction

Gamma-ray binaries are composed of an early-type 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 gamma-ray regime > 1 MeV (see Dubus 2013; Paredes & Bordas 2019, for a review). They exhibit broadband non-thermal emission ranging from radio through low-energy (LE, 1–100 MeV) up to high-energy (HE, > 100 MeV) and very high energy (VHE, > 100 GeV) gamma-rays. For most systems, observations show flux variations with the orbital phase. At the time of writing, nine gamma-ray binaries emitting in the HE regime are confirmed: 1FGL J1018.6-5856 (Fermi LAT Collaboration 2012; H.E.S.S. Collaboration 2015a), 4FGL J1405-6119 (Corbet et al. 2019), HESS J0632+057 (Aharonian et al. 2007), HESS J1832-093 (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 B1259-63 (Aharonian et al. 2005b), and PSR J2032+4127 (Abeysekara et al. 2018).

In the literature, two possible mechanisms are proposed to explain the non-thermal emission (see e.g. Mirabel 2006; Romero et al. 2007): A jet-related emission scenario, where the compact object accretes matter from its stellar companion, releasing the accretion energy in the form of relativistic jets with high-energy particles (see e.g. Bosch-Ramon & Khangulyan 2009); or a wind-driven 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 early-type companion star (see Maraschi & Treves 1981; Dubus 2006). In close-orbit binaries, this interaction leads to the formation of an extended wind-collision region (WCR) that terminates both the pulsar and stellar wind by shocks (see e.g. Bogovalov et al. 2008; Bosch-Ramon et al. 2012). At these sites, particles can be accelerated to ultra-relativistic energies through diffusive shock acceleration and other processes (Sironi & Spitkovsky 2011), which then produce the observed radiation. For many systems, a wind-driven interpretation seems favoured. This has been investigated using a variety of different approaches ranging from simple few-zone models to more complex numerical models including different emission mechanisms (see e.g. Khangulyan et al. 2007; Bosch-Ramon et al. 2008; Cerutti et al. 2010; Bednarek 2011; Zabalza et al. 2013; Takata et al. 2014; Dubus et al. 2015; Barkov & Bosch-Ramon 2018; Molina & Bosch-Ramon 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 wind-driven scenario and aim to alleviate some of the encountered issues.

In analogy to gamma-ray binaries, extended WCRs can also be formed in colliding-wind binaries (CWBs). These systems are very similar in terms of wind interaction, but they harbour a second early-type star instead of the neutron star. CWBs can therefore be viewed as a non-relativistic analogue of gamma-ray binaries in the wind-driven 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 Fokker-Planck type particle transport equation for protons and electrons simultaneously to the equations of magneto-hydrodynamics.

In this work, we extend this approach to gamma-ray 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 gamma-ray binaries before. In turn, this enables the prediction of the non-thermal emission from gamma-ray binaries in a self-consistent 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 inverse-Compton 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 post-processing 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

$$\begin{array}{c}\hfill {\partial }_t\mathit{U}+{\partial }_i{\mathit{F}}^i\left(\mathit{U}\right)=\mathit{S},\end{array}$$(1)

where U denotes the conserved quantities, Fⁱ 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 four-velocity u^μ, and specific enthalpy h as

$$\begin{array}{c}\hfill \mathit{U}=\left(\begin{array}{c}D\\ E\\ m^j\end{array}\right)=\left(\begin{array}{c}\rho u^0\\ Dh{u}^0-p\\ Dh{u}^j\end{array}\right)\end{array}$$(2)

and

$$\begin{array}{c}\hfill {\mathit{F}}^i(\mathit{U})=\left(\begin{array}{c}D{u}^i/u^0\\ m^i\\ m^j{u}^i/u^0+p{\delta }^{\mathit{ij}}\end{array}\right).\end{array}$$(3)

In this paper we use the normalisation c = 1 and notate the three spatial components of the fluid four-velocity vector as u = uⁱ and its Lorentz factor as u⁰. 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,

$$\begin{array}{c}\hfill h(\rho ,p)=1+\frac{\mathrm{\Gamma }}{\mathrm{\Gamma }-1}\frac{p}{\rho },\end{array}$$(4)

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 (Bosch-Ramon et al. 2015).

In our simulation, we employed the HLLC solver (Mignone & Bodo 2005) alongside a second-order sub-grid 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,

$$\begin{array}{c}\hfill {\partial }_t(sD)+{\partial }_i(sD{u}^i/u^0)=0.\end{array}$$(5)

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

$$\begin{array}{c}\hfill \rho \left(\mathit{r}\right)=\frac{\dot{M}}{4\pi r^{2}u},\mathit{u}(\mathit{r})=u{\mathit{e}}_r,p(\mathit{r})={10}^{-6}\rho \left(\mathit{r}\right),\end{array}$$(6)

where r denotes the vector originating at the injection centre to a given point in space, Ṁ the mass injection rate, and u the four-speed 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 mass-loss 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 gamma-ray binaries, it is common practice to prescribe the pulsar mass-loss rate through the ratio η of the total momentum. The pulsar mass-loss rate can hence be obtained by

$$\begin{array}{c}\hfill {\dot{M}}_{\mathrm{p}}=\eta \frac{{\dot{M}}_{\mathrm{s}}{u}_{\mathrm{s}}}{u_{\mathrm{p}}},\end{array}$$(7)

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 time-step 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. Bosch-Ramon 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 tail-like 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:

$$\begin{array}{cc}\hfill 0& ={\mathrm{\nabla }}_{\mu }\left(u^{\mu }{f^{\prime }}_0+q^{\mu }\right)+\frac{1}{p^2}\frac{\partial }{\partial p^{\prime }}[-\frac{{p^{\prime }}^3}{3}{\mathrm{\nabla }}_{\mu }{u}^{\mu }{f^{\prime }}_0+{\dot{p^{\prime }}}_{\mathrm{loss}}{p^{\prime }}^2{f^{\prime }}_0\hfill \\ \hfill & -{\mathrm{\Gamma }}_{\mathrm{visc}}{p^{\prime }}^4\tau \frac{\partial f{\prime }_0}{\partial p\prime }-p^{\prime 2}{D}_{\mathit{pp}}\frac{\partial f{\prime }_0}{\partial p\prime }-p\prime {(p^{\prime 0})}^2{\dot{u}}_{\mu }{q}^{\mu }],\hfill \end{array}$$(8)

where $f_0^{\prime }(x^{\mu },p^{\prime })$ 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 second-order Fermi process, and losses due to non-inertial frame changes. Following Vaidya et al. (2018), we neglected (for the purpose of this paper) spatial diffusion (q^μ = 0), viscous shear-acceleration (Γ_visc = 0), and momentum diffusion (D_pp = 0). Because we assume the ultra-relativistic limit, the particle Lorentz factor can be readily obtained as γ ≈ p′/m, leading us to define ${\mathcal{N}}^{\prime }(x^{\mu },{\gamma }^{\prime })=4\pi p^{\prime 2}{f}_0^{\prime }m$, 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

$$\begin{array}{c}\hfill {\mathrm{\nabla }}_{\mu }(u^{\mu }\mathcal{N}\prime )+\frac{\partial }{\partial \gamma \prime }[\underset{=⟨\frac{\mathrm{d}\gamma \prime }{\mathrm{d}t\prime }⟩}{\underbrace{(-\frac{\gamma \prime }{3}{\mathrm{\nabla }}_{\mu }{u}^{\mu }+{\dot{{\gamma }^{\prime }}}_{\mathrm{rad}})}}\mathcal{N}\prime ]=0,\end{array}$$(9)

with $⟨\frac{\text{d}\gamma \prime }{\text{d}t\prime }⟩$ 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. Bosch-Ramon & 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,

$$\begin{array}{cc}\hfill {\partial }_t(u^0\mathcal{N}\prime )+{\mathrm{\nabla }}_i(u^i\mathcal{N}\prime )& =0\hfill \end{array}$$(10)

$$\begin{array}{cc}\hfill {\partial }_t(u^0\mathcal{N}\prime )+{\partial }_{\gamma \prime }(\langle \frac{\mathrm{d}\gamma \prime }{\mathrm{d}t\prime }\rangle \mathcal{N}\prime )& =0.\hfill \end{array}$$(11)

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 $[{\gamma }_{l-1/2}^{\prime },{\gamma }_{l+1/2}^{\prime }]$. In our implementation, we allow the user to specify the limits of the simulated spectral range, that is, the first and the last edge ${\gamma }_{-1/2}^{\prime },{\gamma }_{N_{\gamma }-1/2}^{\prime }$, and the number of spectral bins N_γ. The range was then subdivided into logarithmic bins with

$$\begin{array}{c}\hfill \gamma {\prime }_{l-1/2}=\gamma {\prime }_{-1/2}{\left(\frac{\gamma {\prime }_{N_{\gamma }-1/2}}{\gamma {\prime }_{-1/2}}\right)}^{l/N_{\gamma }},\mathrm{for}l\in \}0,\cdots ,N_{\gamma }\{.\end{array}$$(12)

The cell centre ${\gamma }_l^{\prime }$ is then obtained as the arithmetic mean of its limits.

Averaging Eq. (10) over the l-th energy bin yields the finite-volume version of the spatial particle transport,

$$\begin{array}{c}\hfill {\mathrm{\nabla }}_{\mu }\left(u^{\mu }\mathcal{N}{\prime }_l\right)=0\end{array}$$(13)

with the cell averages ${\mathcal{N}}_l^{\prime }=\Delta {\gamma }^{\prime }{}_l^{-1}{\int }_{{\gamma }_{l-1/2}^{\prime }}^{{\gamma }_{l+1/2}^{\prime }}{\mathcal{N}}^{\prime }\text{d}{\gamma }^{\prime }$ and $\mathrm{\Delta }{\gamma }_l^{\prime }={\gamma }_{l+1/2}^{\prime }-{\gamma }_{l-1/2}^{\prime }$. The spatial particle distribution corresponding to a given spectral bin can therefore be represented by an additional three-dimensional scalar field ${\mathcal{N}}_l^{\prime }$. 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

$$\begin{array}{c}\hfill {\partial }_t\mathcal{N}{\prime }_{i,j,k}+{\partial }_{\gamma \prime }\left(\frac{1}{u^0}⟨\frac{\mathrm{d}\gamma \prime }{\mathrm{d}t\prime }⟩{\mathcal{N}\prime }_{i,j,k}\right)=0,\end{array}$$(14)

where ${\mathcal{N}}^{\prime }{}_{i,j,k}$ 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 energy-dependent velocity $\dot{\gamma }\prime =⟨\frac{\text{d}\gamma \prime }{\text{d}t\prime }⟩/u^0$.

The evolution in energy was treated independently for each spatial cell by a semi-Lagrangian 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 $\mathcal{N}\equiv {\mathcal{N}}_{i,j,k}^{\prime }$ 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

$$\begin{array}{c}\hfill \frac{\mathrm{d}}{\mathrm{d}t}\gamma =\dot{\gamma }(\gamma ),\end{array}$$(15)

which can be exploited to construct a conservative and explicit update scheme for the particle number densities.

We obtained the particle number density ${\mathcal{N}}_l^{n+1}$ 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ⁿ. This was done by performing a backward integration of the bin edges γ_l ± 1/2 following Eq. (15) to identify the corresponding values ${\gamma }_{l\pm 1/2}^n$ at time tⁿ (denoting propagated edges with a superscript). Equation (15) was solved numerically using a second-order explicit Runge-Kutta method, for example Ralston’s method, with the time step provided by the RHD solver.

After the interval edges ${\gamma }_{l\pm 1/2}^n$ were identified, the updated number density was obtained using particle conservation

$$\begin{array}{c}\hfill \mathrm{\Delta }{\gamma }_l{\mathcal{N}}_l^{n+1}={\displaystyle {\int }_{{\gamma }_{l-1/2}^n}^{{\gamma }_{l+1/2}^n}{\mathcal{N}}^n\mathrm{d}\gamma =I^n\left[{\gamma }_{l+1/2}^n\right]-I^n\left[{\gamma }_{l-1/2}^n\right],}\end{array}$$(16)

where Iⁿ[γ] denotes the cumulative particle count up to a value γ, for the spectrum at time tⁿ (see also Eq. (19)). To approximate the integral, we considered a second-order, piece-wise linear reconstruction in energy in each cell, in analogy to the spatial reconstruction in the hydrodynamics part,

$$\begin{array}{c}\hfill {\mathcal{N}}^n(\gamma )={\mathcal{N}}_l^n+{\delta }_l^n(\gamma -{\gamma }_l),\end{array}$$(17)

with γ_{l − 1/2} <  γ <  γ_l + 1/2. Here, ${\delta }_l^n$ denotes the slope of the reconstruction obtained after applying the van Leer limiter (van Leer 1977),

$$\begin{array}{c}\hfill {\delta }_{\mathrm{van}\mathrm{Leer}}=\frac{max({\delta }^{+}{\delta }^{-},0)}{{\delta }^C},\end{array}$$(18)

to the right-, left-, and two-sided finite differences, denoted by δ⁺, δ⁻, and δ^C, respectively. With this, the integral can be approximated directly as piecewise parabolas by

$$\begin{array}{cc}\hfill I^n\left[\gamma \right]& ={\displaystyle {\int }_{{\gamma }_{-1/2}}^{\gamma }{\mathcal{N}}^n(\stackrel{\sim }{\gamma })\mathrm{d}\stackrel{\sim }{\gamma }}\hfill \\ \hfill & \approx I^n\left[{\gamma }_{l-1/2}\right]+({\mathcal{N}}_l^n-{\delta }_l^n{\gamma }_l)h_{\gamma }+\frac{{\delta }_l^n}{2}{h}_{\gamma }^2,\hfill \end{array}$$(19)

with γ_{l − 1/2} <  γ <  γ_l + 1/2, h_γ = γ − γ_{l − 1/2}, and the exact integrals at the cell edges,

$$\begin{array}{c}\hfill I^n\left[{\gamma }_{l-1/2}\right]={\displaystyle {\int }_{{\gamma }_{-1/2}}^{{\gamma }_{l-1/2}}{\mathcal{N}}^n\mathrm{d}\gamma =\sum _{j=0}^{l-1}\mathrm{\Delta }{\gamma }_j{\mathcal{N}}_j^n,}\end{array}$$(20)

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ⁿ in the evaluation of the right-hand side of Eq. (16). To determine the slope of the first and last spectral bin according to Eq. (18), we added a ghost-cell 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 & Bosch-Ramon 2018). The strength of the post-shock magnetic field B′ in the fluid frame can therefore be expressed as

$$\begin{array}{c}\hfill \frac{B^{\prime 2}}{2{\mu }_0}={\zeta }_b\frac{p}{\mathrm{\Gamma }-1}.\end{array}$$(21)

We assumed no direction on the magnetic field.

2.2.3. Radiative losses

As energy-loss processes for the accelerated electron-positron pairs we considered synchrotron and inverse-Compton losses. The energy-loss rates are given by

$$\begin{array}{c}\hfill \dot{\gamma }{\prime }_{\mathrm{syn}}=\frac{4}{3}\frac{{\sigma }_T}{m_{e}c}{w}_B\prime {\gamma }^{\prime 2}\end{array}$$(22)

$$\begin{array}{c}\hfill \dot{\gamma }{\prime }_{\mathrm{ic}}=\frac{4}{3}\frac{{\sigma }_T}{m_{e}c}{w}_{}\text{rad}\prime {\gamma }^{\prime 2}{F}_{\mathrm{KN}}({ϵ}_0\prime ,\gamma \prime ),\end{array}$$(23)

with $w_B{}^{\prime }$ and $w_{\text{rad}}^{\prime }$ being the magnetic and radiative energy densities in the fluid frame, respectively, and F_KN the Klein-Nishina factor.

The magnetic field in the fluid frame can be readily obtained by the transformation

$$\begin{array}{c}\hfill \mathit{B}\prime =\frac{1}{u^0}[\mathit{B}+\frac{\mathit{B}·\mathit{u}}{u^0+1}\mathit{u}].\end{array}$$(24)

The corresponding magnetic field strength is therefore given by

$$\begin{array}{c}\hfill B\prime =\frac{1}{u^0}\sqrt{B^2+{(\mathit{B}·\mathit{u})}^2},\end{array}$$(25)

yielding the magnetic energy density ${\mathit{w}}_B\prime =\frac{{B\prime }^2}{2{\mu }_0}$. 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 inverse-Compton losses, we assumed the radiation field to be monochromatic and mono-directional. The energy density of a beam of photons with direction Ω₀ is given in the fluid frame by Doppler boosting the laboratory energy density,

$$\begin{array}{c}\hfill w_{}\text{rad}\prime ={\mathcal{D}}_{u{\mathrm{\Omega }}_0}^{-2}{w}_{}\text{rad},\end{array}$$(26)

with the Doppler factor 𝒟_uΩ0 = (u⁰−u⋅Ω₀)⁻¹. The radiation field relevant for the inverse-Compton losses emerges from the stellar binary companion,

$$\begin{array}{c}\hfill w_{}\text{rad}=\frac{L_{\mathrm{star}}}{4\pi r_{\mathrm{star}}^{2}c},\end{array}$$(27)

where r_⋆ is the distance to the star and L_⋆ its luminosity. We employed the Klein-Nishina factor following Moderski et al. (2005),

$$\begin{array}{cc}\hfill F_{\mathrm{KN}}(\stackrel{\sim }{b}\prime )& =\frac{9}{{\stackrel{\sim }{b}}^{\prime 3}}[(\frac{1}{2}\stackrel{\sim }{b\prime }+6+\frac{6}{\stackrel{\sim }{b\prime }})log(1+\stackrel{\sim }{b\prime })+2{\mathrm{Li}}_2(-\stackrel{\sim }{b\prime })\hfill \\ \hfill & -(\frac{11}{12}{\stackrel{\sim }{b}}^{\prime 3}+6{\stackrel{\sim }{b}}^{\prime 2}+9\stackrel{\sim }{b\prime }+4)·{(1+\stackrel{\sim }{b\prime })}^{-2}-2],\hfill \end{array}$$(28)

assuming an isotropic seed photon and electron distribution for simplicity, with $\stackrel{\sim }{b}\prime =4{ϵ\prime }_0\gamma \prime {(m_e{c}^2)}^{-1}$ and ${ϵ}_0^{\prime }={\mathcal{D}}_{u{\mathrm{\Omega }}_0}^{-1}{ϵ}_0$ denoting the Doppler-boosted 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

$$\begin{array}{c}\hfill {\mathrm{\nabla }}_{\mu }{u}^{\mu }<{\delta }_{\mathrm{sh}},\end{array}$$(29)

with a typical value of δ_sh = −10 c/AU to ignore weakly converging flows (see also Appendix C for a comparison to a common shock-identification 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 non-thermal accelerated pairs, following a power-law spectrum (similar to Dubus et al. 2015), which are parametrised as

$$\begin{array}{c}\hfill {\mathcal{N}}^{\mathrm{MW}}=K^{\mathrm{MW}}{\gamma }^{2}exp(-\gamma /{\gamma }_t)\end{array}$$(30)

$$\begin{array}{c}\hfill {\mathcal{N}}^{\mathrm{PL}}=K^{\mathrm{PL}}{\gamma }^{-s}\mathrm{for}{\gamma }_{\min }^{\mathrm{PL}}<\gamma <{\gamma }_{\max }^{\mathrm{PL}}{\mathcal{N}}^{\mathrm{MW}}\end{array}$$(31)

The maximum Lorentz factor ${\gamma }_{\text{max}}^{\text{PL}}$ 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 $R_L=\frac{\gamma m}{eB\prime }$ and the acceleration efficiency ξ_acc. We assumed synchrotron losses to be dominant for the particles at highest energies. Together with the loss timescale ${\tau }_{\mathrm{loss}}=\gamma \prime /\dot{\gamma }\prime$ given by Eq. (22), this yields the maximum particle Lorentz factor expected at the current location,

$$\begin{array}{c}\hfill {\gamma }_{\max }^{\mathrm{PL}}=\sqrt{\frac{3}{4\pi }\frac{{\mu }_{0}ec}{{\xi }_{\mathrm{acc}}{\sigma }_{T}B\prime }}\approx 4.65\times {10}^7\sqrt{\frac{1\mathrm{G}}{{\xi }_{\mathrm{acc}}B\prime }}.\end{array}$$(32)

Because the dominant particle acceleration mechanism in gamma-ray 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 & Bosch-Ramon 2020).

The average Lorentz factor of the Maxwellian γ_t, the lower cutoff of the power-law spectrum γ_min, and both normalisations K^PLMW were determined from fluid properties. This was controlled by three free parameters: ${\zeta }_n^{\text{PL}}$ the fraction of accelerated non-thermal pairs, ${\zeta }_e^{\text{PL}}$ the fraction of internal energy converted to non-thermal 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 ${\zeta }_n^{\text{MW}}=1-{\zeta }_n^{\text{PL}}$ and ${\zeta }_e^{\text{MW}}=1-{\zeta }_e^{\text{PL}}$. In our simulation, we kept all injection parameters constant in time and space. The injected spectra are therefore related to fluid quantities by

$$\begin{array}{c}\hfill {\displaystyle {\int }_{{\gamma }_{\min }^{\mathrm{PL},\mathrm{MW}}}^{{\gamma }_{\max }^{\mathrm{PL},\mathrm{MW}}}{\mathcal{N}}^{\mathrm{PL},\mathrm{MW}}(\gamma )\mathrm{d}\gamma =\underset{=:n^{\mathrm{PL},\mathrm{MW}}}{\underbrace{{\zeta }_n^{\mathrm{PL},\mathrm{MW}}{\zeta }_{\rho }{\psi }_{\mathrm{p}}\rho /m_e}}}\end{array}$$(33)

$$\begin{array}{c}\hfill {\displaystyle {\int }_{{\gamma }_{\min }^{\mathrm{PL},\mathrm{MW}}}^{{\gamma }_{\max }^{\mathrm{PL},\mathrm{MW}}}\gamma m_e{c}^2{\mathcal{N}}^{\mathrm{PL},\mathrm{MW}}(\gamma )\mathrm{d}\gamma =\underset{=:e^{\mathrm{PL},\mathrm{MW}}}{\underbrace{{\zeta }_e^{\mathrm{PL},\mathrm{MW}}{\psi }_{\mathrm{p}}{e}_{\mathrm{thermal}}}},}\end{array}$$(34)

where 0 <  ψ_p <  1 is a passive tracer indicating the fraction of pulsar wind material and the thermal energy density $e_{\mathrm{thermal}}=\frac{p}{\mathrm{\Gamma }-1}$ given by the fluid pressure.

For the non-thermal particles, the above expressions yield an implicit equation for ${\gamma }_{\text{min}}^{\text{PL}}$, which was solved numerically by a standard root-finding algorithm. The injection normalisation K^PL was then determined from Eq. (33). When the internal fluid energy was too low,

$$\begin{array}{c}\hfill {\displaystyle {\int }_{{\gamma }_{\min }^{\mathrm{PL}}}^{{\gamma }_{\max }^{\mathrm{PL}}}\gamma m_e{c}^2{\mathcal{N}}^{\mathrm{PL}}(\gamma )\mathrm{d}\gamma >e^{\mathrm{PL}},}\end{array}$$(35)

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 ${\gamma }_{\text{min}}^{\text{MW}}=1$ and ${\gamma }_{\text{max}}^{\text{MW}}=\infty$, yielding

$$\begin{array}{c}\hfill {\gamma }_t=\frac{1}{3}\frac{e^{\mathrm{MW}}}{n^{\mathrm{MW}}{m}_e{c}^2},K^{\mathrm{MW}}=\frac{n^{\mathrm{MW}}}{2{\gamma }_t^3}.\end{array}$$(36)

2.3. Emission

The particle distributions resulting from a joint hydrodynamic and particle-transport simulation (as detailed in Sects. 2.1 and 2.2) were used in a post-processing 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 inverse-Compton 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

$$\begin{array}{c}\hfill j_{ϵ}(ϵ,\mathbf{\Omega })={\mathcal{D}}_{u\mathrm{\Omega }}^{2}j{\prime }_{ϵ}(ϵ\prime ,\mathbf{\Omega }\prime ),\end{array}$$(37)

where 𝒟_uΩ = (u⁰−Ω⋅u)⁻¹ is the corresponding Doppler factor. The respective energy and direction of the emitted photons are given in the fluid frame by

$$\begin{array}{c}\hfill ϵ\prime ={\mathcal{D}}_{u\mathrm{\Omega }}^{-1}ϵ,\mathbf{\Omega }\prime ={\mathcal{D}}_{u\mathrm{\Omega }}[\mathbf{\Omega }+(\frac{\mathrm{\Omega }·\mathit{u}}{u^0+1}-1)\mathit{u}].\end{array}$$(38)

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

$$\begin{array}{c}\hfill j{\prime }_{ϵ}(ϵ\prime ,\mathbf{\Omega }\prime )={\displaystyle \int }\mathcal{P}\prime (ϵ\prime ,\mathbf{\Omega }\prime ,\gamma \prime )\mathcal{N}\prime \mathrm{d}\gamma \prime .\end{array}$$(39)

2.3.1. Synchrotron emission

The total power emitted by a single electron in a magnetic field B′ can be expressed as

$$\begin{array}{c}\hfill {{\mathcal{P}}^{\prime }}^{}{\text{tot}}_{}\text{syn}({ϵ}^{\prime },{\gamma }^{\prime })=\frac{\sqrt{3}{e}^3}{8{\pi }^2{\epsilon }_0{m}_{e}cħ}{B^{\prime }}_{\perp }F\left(\frac{{ϵ}^{\prime }}{{ϵ}_{c^{\prime }}}\right),\end{array}$$(40)

which is peaked around the critical photon energy,

$$\begin{array}{c}\hfill {ϵ}_c\prime =\frac{3eħ}{2m_e}B{\prime }_{\perp }{\gamma }^{\prime 2}=1.74\times {10}^{-8}\text{eV}\frac{B{\prime }_{\perp }}{1\text{G}}{\gamma }^{\prime 2}.\end{array}$$(41)

$B_{\perp }^{\prime }=B^{\prime }sin(\alpha )$ denotes the magnetic field component perpendicular to the electron motion, which gyrates with a pitch angle α around the magnetic field lines. For ultra-relativistic 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_{\perp }^{\prime }$ = |B′×Ω′|. Together with the assumption of an isotropic particle distribution, this yields

$$\begin{array}{c}\hfill \mathcal{P}{\prime }_{}\text{syn}(ϵ\prime ,\mathbf{\Omega }\prime ,\gamma \prime )=\frac{\sqrt{3}{e}^3}{32{\pi }^3{\epsilon }_0{m}_{e}cħ}B{\prime }_{\perp }F\left(\frac{ϵ\prime }{{ϵ}_c\prime }\right).\end{array}$$(42)

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 $B_{\perp }^{\prime }$ in Eqs. (41), (42) by its pitch-angle average ${\overline{B}\prime }_{\perp }=\frac{\pi }{4}B\prime$.

2.3.2. Inverse-Compton emission

Through the inverse-Compton process, photons of the stellar radiation field are scattered to gamma-ray energies by highly energetic particles. Following the description of Moderski et al. (2005), the power emitted through this process is given by

$$\begin{array}{c}\hfill \mathcal{P}{\prime }_{}\text{ic}(ϵ\prime ,\mathbf{\Omega }\prime ,\gamma \prime )=ϵ\prime \frac{\partial \dot{N}{\prime }_{}\text{sc}}{\partial ϵ\prime \partial \mathrm{\Omega }\prime },\end{array}$$(43)

where $\frac{\partial {\dot{N}\prime }_{}\text{sc}}{\partial ϵ\prime \partial \mathrm{\Omega }\prime }$ is the scattering rate of photons on isotropically distributed electrons. The scattering rate was presented by Aharonian & Atoyan (1981) in the case of ${ϵ}^{\prime }\gg {ϵ}_0^{\prime }$ and γ′≫1 taking the form

$$\begin{array}{cc}\hfill \frac{\partial \dot{N}{\prime }_{}\text{sc}}{\partial ϵ\prime \partial \mathrm{\Omega }\prime }& =\frac{3}{16\pi }c{\sigma }_T\frac{1}{{\gamma }^{\prime 2}}\hfill \\ \hfill & {\displaystyle \int {\int }_{ϵ{\prime }_{0,\text{min}}}^{\infty }}{f}_{}\text{ic}(ϵ\prime ,ϵ{\prime }_0,\gamma \prime ,\mu \prime )\frac{n_{{ϵ}_0}\prime ({ϵ}_0\prime ,{\mathbf{\Omega }}_0\prime )}{{ϵ}_0\prime }\mathrm{d}ϵ{\prime }_0{\mathrm{d}}^2{\mathrm{\Omega }}_0\prime ,\hfill \end{array}$$(44)

where n₀(ϵ₀, Ω₀) is the density of the seed photons with energy ϵ₀ and direction Ω₀ and ${ϵ}_{0,\text{min}}^{\prime }$ the lowest seed photon energy that is still scattered to a given energy ϵ′ and direction Ω′

$$\begin{array}{c}\hfill {ϵ}_{0,\text{min}}\prime =ϵ\prime {\left[2(1-\mu \prime )(1-\frac{ϵ\prime }{\gamma \prime m_e{c}^2})\right]}^{-1}.\end{array}$$(45)

The scattering kernel f_ic is given by

$$\begin{array}{c}\hfill f_{}\text{ic}(ϵ\prime ,ϵ{\prime }_0,\gamma \prime ,\mu \prime )=1+\frac{w^{\prime 2}}{2(1-w\prime )}-\frac{2w\prime }{b{\prime }_{\mu }(1-w\prime )}+\frac{2w^{\prime 2}}{{b{\prime }_{\mu }}^2{(1-w\prime )}^2},\end{array}$$(46)

with $b_{\mu }^{\prime }=2(1-{\mu }^{\prime }){ϵ}_0^{\prime }{\gamma }^{\prime }{(m_e{c}^2)}^{-1}$, w′=ϵ′(γ′m_ec²)⁻¹ and the cosine of the scattering angle ${\mu }^{\prime }={\mathbf{\Omega }}^{\prime }\cdot {\mathbf{\Omega }}_0^{\prime }$, which can be expressed in terms of laboratory quantities by

$$\begin{array}{c}\hfill \mu \prime =1+{\mathcal{D}}_{u{\mathrm{\Omega }}_0}{\mathcal{D}}_{u\mathrm{\Omega }}(\mu -1).\end{array}$$(47)

Together with the relativistic invariant $\frac{n_{{ϵ}_0}\prime }{{ϵ}_{0\prime }}\text{d}{ϵ\prime }_0{\text{d}}^2{\mathrm{\Omega }}_0\prime =\frac{n_{{ϵ}_0}}{{ϵ}_0}\text{d}{ϵ}_0{\text{d}}^2{\mathrm{\Omega }}_0$ Eq. (44) can be expressed directly in terms of the seed photon field in the laboratory frame,

$$\begin{array}{cc}& \frac{\partial \dot{N}{\prime }_{}\text{sc}}{\partial ϵ\prime \partial \mathrm{\Omega }\prime }=\frac{3}{16\pi }c{\sigma }_T\frac{1}{{\gamma }^{\prime 2}}{\displaystyle \int {\int }_{{ϵ}_{0,\text{min}}}^{\infty }}f(ϵ\prime ,ϵ{\prime }_0,\gamma \prime ,\mu \prime )\frac{n_{{ϵ}_0}({ϵ}_0,{\mathbf{\Omega }}_{\mathbf{0}})}{{ϵ}_0}\mathrm{d}{ϵ}_0{\mathrm{d}}^2{\mathrm{\Omega }}_0.\hfill \end{array}$$(48)

Applying the factorisation of the target photon field (Appendix B) yields

$$\begin{array}{cc}\hfill \frac{\partial \dot{N}{\prime }_{}\text{sc}}{\partial ϵ\prime \partial \mathrm{\Omega }\prime }& =\frac{3}{16\pi }c{\sigma }_T\frac{1}{{\gamma }^{\prime 2}}{g}_n(\mathit{x}){\displaystyle \int }{g}_{\mathrm{\Omega }}({\mathbf{\Omega }}_0)\hfill \\ \hfill & {\displaystyle {\int }_{{ϵ}_{0,\text{min}}}^{\infty }}f(ϵ\prime ,ϵ{\prime }_0,\gamma \prime ,\mu \prime )\frac{g_{ϵ}({ϵ}_0)}{{ϵ}_0}\mathrm{d}{ϵ}_0{\mathrm{d}}^2{\mathrm{\Omega }}_0,\hfill \end{array}$$(49)

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 gamma-ray flux produced through inverse-Compton 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 gamma-ray 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 three-fold integration of the differential absorption opacity given by Gould & Schréder (1967),

$$\begin{array}{c}\hfill \mathrm{d}{\tau }_{\gamma \gamma }=(1-\mu ){\sigma }_{\gamma \gamma }(\beta )n_{{ϵ}_0}(\mathit{x},{ϵ}_0,{\mathbf{\Omega }}_0)\mathrm{d}l{\mathrm{d}}^2{\mathrm{\Omega }}_0\mathrm{d}{ϵ}_0,\end{array}$$(50)

over the line of sight l, the solid angle Ω₀, and the seed photon energy ϵ₀. Here, μ = Ω ⋅ Ω₀ denotes the scattering angle cosine between a VHE photon with direction Ω and a seed UV photon with direction Ω₀, ${\beta }^2=1-\frac{m_e^2{c}^4}{ϵ{ϵ}_0}\frac{2}{1-\mu }$ the speed of the created e^± pair normalised to c and ${\sigma }_{\gamma \gamma }=\frac{3}{16}{\sigma }_T{f}_{\gamma \gamma }(\beta )$ the scattering cross-section with

$$\begin{array}{c}\hfill f_{\gamma \gamma }(\beta )=(1-{\beta }^2)[(3-{\beta }^4)log\left(\frac{1+\beta }{1-\beta }\right)-2\beta (2-{\beta }^2)].\end{array}$$(51)

The lower bound of the energy integration is provided by the lowest seed photon energy ${ϵ}_{0,\text{min}}=\frac{m_e^2{c}^4}{ϵ}\frac{2}{1-\mu }$ 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 line-of-sight integration, we therefore used the transformation $l=d_0\frac{sin(\psi -{\psi }_0)}{sin(\psi )}$ (as detailed in Dubus 2006). The integration over target photon energy ϵ₀ was transformed into an integral over the speed of the created pairs β. Together with the radiation field factorisation (Appendix B), this yields

$$\begin{array}{cc}\hfill {\tau }_{\gamma \gamma }& =\frac{3}{16}{\sigma }_T{\displaystyle {\int }_{{\psi }_0}^{\pi }\frac{\mathrm{d}l}{\mathrm{d}\psi }{g}_n(\mathit{x}\prime )\int (1-\mu )g_{\mathrm{\Omega }}(\mathit{x}\prime ,{\mathbf{\Omega }}_0)}\hfill \\ \hfill & {\displaystyle {\int }_0^1}\frac{\mathrm{d}{ϵ}_0}{\mathrm{d}\beta }{g}_{ϵ}(\mathit{x}\prime ,{ϵ}_0)f_{\gamma \gamma }(\beta )\mathrm{d}\beta {\mathrm{d}}^2{\mathrm{\Omega }}_0\mathrm{d}\psi ,\hfill \end{array}$$(52)

which can be further simplified depending on the stellar radiation field model.

The integration over solid angle Ω₀ is approximated by a nested quadrature rule with equidistant points in azimuth and a Gauss-Legendre quadrature in the polar angle (see Hesse & Womersley 2012). For the integrals over ψ and β, we employed standard Gauss-Legendre 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 O-type 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 mass-loss rate of Ṁ_s = 2 × 10⁻⁸ M_⊙ yr⁻¹ launched at a velocity of v_s = 2000 km s⁻¹. The pulsar wind was injected at a speed of v_p = 0.99 c and a mass-loss 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³ 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 four-speed, 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. Bosch-Ramon 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 counter-clockwise 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 Bosch-Ramon 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 electron-positron 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⁸]. The particle injection was controlled as detailed in Sect. 2.2.4 using ζ_ρ = 5 × 10⁻⁴, ${\zeta }_n^{\text{PL}}=5\times {10}^{-3}$, ${\zeta }_e^{\text{PL}}=0.5$, 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 energy-loss 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 × 103 (upper panel) and γ = 1.1 × 107 (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 power-law 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 power-law electrons can be easily distinguished, with Maxwellian electrons dominating the low-energy part (γ ≲ 2 × 10⁴) and the non-thermal electrons the high-energy part. For the chosen set of parameters, Maxwellian electrons were on average injected with γ_t ∼ 1.0 × 10³ and non-thermal electrons within the limits of γ_min ∼ 1.1 × 10⁵ and γ_max = 2.4 × 10⁷.

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 post-processing 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 inverse-Compton 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⁵ 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 X-rays up to LE gamma-rays, the spectrum is dominated by synchrotron emission from electrons at the power-law tail of the spectrum. The same population of electrons is responsible for the inverse-Compton scattering of stellar UV photons to the VHE gamma-ray 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 gamma-ray regime is produced through the inverse-Compton process, but by the low-energy 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 ${\zeta }_e^{\text{PL,MW}}/({\zeta }_{\rho }{\zeta }_n^{\text{PL,MW}})$. In the case of the Maxwellian electrons, this can be used to shift the inverse-Compton bump at HE in energy. For the non-thermal power-law electrons, their specific energy density affects the low-energy cutoff of the spectrum. The overall normalisation of the spectrum is given by the amount of energy deposited in the respective populations, determined by ${\zeta }_e^{\text{PL,MW}}$. The index of the injected electrons translates into the spectral index of the X-ray synchrotron emission and affects the spectral index at VHE gamma-rays. ζ_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 gamma-ray band, depending on the magnetic field in this region, and the VHE gamma-ray 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 X-ray and LE gamma-ray 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 gamma-rays, the anisotropic inverse-Compton scattering cross-section introduces a further line-of-sight-dependent 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 anti-correlation of the HE gamma-ray band and the previously discussed one.

The VHE gamma-ray emission is again produced by inverse-Compton process, but the upscattered photons are attenuated by γγ-pair creation with the stellar photon field. This once more introduces a line-of-sight-dependent 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 X-ray to LE gamma-ray band and therefore anti-correlated to HE gamma-rays.

4. Summary and outlook

We presented a novel numerical model for gamma-ray binaries treating the transport of shock-accelerated electron-positron 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 self-consistent three-dimensional 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 inverse-Compton 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 bow-like pulsar-wind shock and a Coriolis shock behind the pulsar. The expressed features are in qualitative agreement with previous works (see e.g. Bosch-Ramon 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 gamma-ray 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 gamma-rays. 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 X-ray to LE gamma-ray synchrotron emission is given by relativistic boosting in the wings of the WCR. In addition to the modulation introduced by the anisotropic inverse-Compton scattering producing the HE and VHE gamma-ray emission, the VHE band is further modulated by γγ-absorption. For the investigated general binary, the model naturally predicts a correlated X-ray to LE gamma-ray and VHE gamma-ray emission, while being anti-correlated to the HE band, as is seen in the observations of many gamma-ray binaries (e.g. LS 5039, Chang et al. 2016).

Because many of the observed gamma-ray 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 well-studied 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 gamma-ray 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 self-consistent 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 non-inertial frame changes. This effect has not been studied in the context of gamma-ray 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 gamma-ray 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 gamma-ray binary emissions (e.g. PSR B1259-63 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 Joliot-Curie 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

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

$$\begin{array}{cc}\hfill \overline{t}& =t\hfill \end{array}$$(A.1)

$$\begin{array}{cc}\hfill \overline{x}& =xcos\left(\mathrm{\Omega }t\right)+ysin\left(\mathrm{\Omega }t\right)\hfill \end{array}$$(A.2)

$$\begin{array}{cc}\hfill \overline{y}& =ycos\left(\mathrm{\Omega }t\right)-xsin\left(\mathrm{\Omega }t\right)\hfill \end{array}$$(A.3)

$$\begin{array}{cc}\hfill \overline{z}& =z.\hfill \end{array}$$(A.4)

Transforming the flat Cartesian Minkowski metric into the corotating system yields

$$\begin{array}{cc}& g_{\overline{\mu }\overline{\nu }}=0.4\left(\begin{array}{cccc}-1+{\mathrm{\Omega }}^2({\overline{x}}^2+{\overline{y}}^2)& -\mathrm{\Omega }\overline{y}& \mathrm{\Omega }\overline{x}& 0\\ -\mathrm{\Omega }\overline{y}& 1& 0& 0\\ \mathrm{\Omega }\overline{x}& 0& 1& 0\\ 0& 0& 0& 1\end{array}\right)\hfill \end{array}$$(A.5)

$$\begin{array}{cc}\hfill & g^{\overline{\mu }\overline{\nu }}=0.4\left(\begin{array}{cccc}-1& -\mathrm{\Omega }\overline{y}& \mathrm{\Omega }\overline{x}& 0\\ -\mathrm{\Omega }\overline{y}& 1-{\mathrm{\Omega }}^2{\overline{y}}^2& {\mathrm{\Omega }}^2\overline{x}\overline{y}& 0\\ \mathrm{\Omega }\overline{x}& {\mathrm{\Omega }}^2\overline{x}\overline{y}& 1-{\mathrm{\Omega }}^2{\overline{x}}^2& 0\\ 0& 0& 0& 1\end{array}\right).\hfill \end{array}$$(A.6)

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 energy-momentum tensor T^μν = ρhu^μu^ν + pg^μν. Because of the non-vanishing, off-diagonal 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² = −(α² − β_iβⁱ) dt² + 2β_i dxⁱ dt + γ_ij dxⁱ dx^j, with $\overline{\alpha }$ the lapse, ${\overline{\beta }}^i$ the shift, and ${\overline{\gamma }}_{\mathit{ij}}$ 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

$$\begin{array}{cc}\hfill \overline{\alpha }& =1\hfill \end{array}$$(A.7)

$$\begin{array}{cc}\hfill {\overline{\beta }}^i& ={\overline{\beta }}_i={(-\mathrm{\Omega }\overline{y},\mathrm{\Omega }\overline{x},0)}_i\hfill \end{array}$$(A.8)

$$\begin{array}{cc}\hfill {\overline{\gamma }}_{\mathit{ij}}& =\mathbb{1}.\hfill \end{array}$$(A.9)

In our approach, we relied on the so-called 3+1 Valencia formulation of general relativistic hydrodynamics (Banyuls et al. 1997), and chose the conserved variables

$$\begin{array}{c}\hfill \mathit{U}=\left(\begin{array}{c}D\\ E\\ m_j\end{array}\right)=\left(\begin{array}{c}\rho W\\ DhW-p\\ DhW{w}_i\end{array}\right),\end{array}$$(A.10)

consisting of the rest-mass density D, the momentum density in the j-direction m_j, and the total energy density E, measured by a family of Eulerian observers. In their respective frames, the fluid three-velocity is given by wⁱ, which is related to the four-velocity in the rotating system by ${\mathit{w}}^i=(u^{\overline{i}}/u^{\overline{t}}+{\beta }^{\overline{i}})/\alpha$. The corresponding Lorentz factor is given by W = (1−w_iwⁱ)^−1/2 with w_i = γ_ijw^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 ${\gamma }_{\mathit{ij}}=\mathbb{1}$. 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

$$\begin{array}{c}\hfill {\partial }_t(\sqrt{\gamma }\mathit{U})+{\partial }_i(\sqrt{\gamma }{\mathit{F}}^i)=\sqrt{\gamma }\mathit{S},\end{array}$$(A.11)

with γ = det(γ_ij) and the fluxes

$$\begin{array}{c}\hfill {\mathit{F}}^i=\left(\begin{array}{c}\alpha w^{i}D-{\beta }^{i}D\\ \alpha m^i-{\beta }^{i}E\\ \alpha (m_j{w}^i+p{\delta }_j^i)-{\beta }^i{m}_j\end{array}\right).\end{array}$$(A.12)

The system of equations again closely resemble those in special relativity Eq. (1) up to the additional fluxes introduced by βⁱ ≠ 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

$$\begin{array}{c}\hfill {M^{\widehat{i}}}_i{|}_{{\mathrm{\Sigma }}_{\overline{x}}}=\mathbb{1},{M^{\widehat{i}}}_i{{|}_{{\mathrm{\Sigma }}_{\overline{y}}}=\mathbb{1},{M^{\widehat{i}}}_i|}_{{\mathrm{\Sigma }}_{\overline{z}}}=\mathbb{1}.\end{array}$$(A.13)

The source terms in Eq. (A.11) are given by (Banyuls et al. 1997)

$$\begin{array}{c}\hfill \mathit{S}=\left(\begin{array}{c}0\\ {\alpha }^2(T^{\mu t}{\partial }_{\mu }ln\alpha -{\mathrm{\Gamma }}_{\mu \lambda }^t{T}^{\mu \lambda })\\ \alpha T^{\mu \nu }({\partial }_{\mu }{g}_{\nu j}-{\mathrm{\Gamma }}_{\nu \mu }^{\lambda }{g}_{\lambda j})\end{array}\right),\end{array}$$(A.14)

and reduce to

$$\begin{array}{c}\hfill \mathit{S}={(0,0,\mathrm{\Omega }{m}_{\overline{y}},-\mathrm{\Omega }{m}_{\overline{x}},0)}^{\top }\end{array}$$(A.15)

in our case. This result may seem counter-intuitive 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 well-known 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 $u^{\overline{\mu }}$, 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 inverse-Compton 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

$$\begin{array}{c}\hfill n(\mathit{x},ϵ,\mathbf{\Omega })=\frac{\mathrm{d}N}{\mathrm{d}V\mathrm{d}ϵ\mathrm{d}\mathrm{\Omega }}=g_n(\mathit{x})·g_{ϵ}(\mathit{x},ϵ)·g_{\mathrm{\Omega }}(\mathit{x},\mathbf{\Omega }),\end{array}$$(B.1)

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,

$$\begin{array}{c}\hfill g_n(\mathit{x}):=\frac{\mathrm{d}N}{\mathrm{d}V}\end{array}$$(B.2)

$$\begin{array}{c}\hfill g_{ϵ}(\mathit{x},ϵ):={\left(\frac{\mathrm{d}N}{\mathrm{d}V}\right)}^{-1}\frac{\mathrm{d}N}{\mathrm{d}V\mathrm{d}ϵ}\end{array}$$(B.3)

$$\begin{array}{c}\hfill g_{\mathrm{\Omega }}(\mathit{x},\mathbf{\Omega }):={\left(\frac{\mathrm{d}N}{\mathrm{d}V}\right)}^{-1}\frac{\mathrm{d}N}{\mathrm{d}V\mathrm{d}\mathrm{\Omega }}.\end{array}$$(B.4)

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

$$\begin{array}{c}\hfill g_n^{\text{spherical}}(\mathit{x})=\frac{1}{r^2}\frac{{\dot{N}}_{\star }}{4\pi c},\end{array}$$(B.5)

where r = ‖x−x_⋆‖.

A mono-directional or beamed radiation field in a given direction Ω_beam is in general described by

$$\begin{array}{c}\hfill g_{\mathrm{\Omega }}^{\text{beamed}}(\mathit{x},\mathbf{\Omega })=\delta (\mathbf{\Omega }-{\mathbf{\Omega }}_{}\text{beam}).\end{array}$$(B.6)

The radiation field of a point-like 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,

$$\begin{array}{c}\hfill g_{\mathrm{\Omega }}^{\text{extended}}(\mathit{x},\mathbf{\Omega })=\{\begin{array}{cc}{(2\pi (1-{\mu }_{0,\text{min}}))}^{-1}\hfill & {\mu }_0>{\mu }_{0,\text{min}}\hfill \\ 0\hfill & \text{else}\hfill \end{array},\end{array}$$(B.7)

where μ₀ is the cosine of the angle enclosed between Ω and e_r and μ_0, min = (1−R²/r²)^1/2.

The radiation field of the companion star can be approximated as a blackbody spectrum with temperature T,

$$\begin{array}{c}\hfill g_{ϵ}^{\text{blackbody}}(\mathit{x},ϵ)=\frac{1}{2\zeta (3)k_{\mathrm{b}}T}\frac{{(ϵ/k_{\mathrm{b}}T)}^2}{exp(ϵ/k_{\mathrm{b}}T)-1}.\end{array}$$(B.8)

This can be further approximated as a monochromatic spectrum with the average photon energy $\overline{ϵ}=k_{\mathrm{B}}T\frac{{\pi }^4}{30\zeta (3)}$,

$$\begin{array}{c}\hfill g_{ϵ}^{\text{mono}}(\mathit{x},ϵ)=\delta (ϵ-\overline{ϵ}).\end{array}$$(B.9)

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 $R_{\star }^2=L_{\star }/(\frac{4\pi }{3}\sigma T_{\star }^4)$ and its photon injection rate by ${\dot{N}}_{\star }=L_{\star }/{\overline{ϵ}}_{\star }$.

Appendix C: Comparison to a common shock-identification 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 gamma-ray 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 counter-clockwise 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 × 103 (upper panel) and γ = 1.1 × 107 (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 power-law 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