© ESO, 2015

1. Introduction

Gamma-ray binaries are composed of a massive star in orbit with a compact object and characterised by dominant radiative output in the gamma-ray (≳MeV) range (see Dubus 2013 for a review). The compact object is widely thought to be a rotation-powered pulsar, although this remains to be proven for most systems. The interaction of the pulsar wind with the massive star wind ends up dissipating part of the pulsar’s rotation power through particle acceleration. Gamma-ray binaries offer an opportunity to study the processes involved in pulsar wind nebulae on much smaller scales, notably how an e⁺e⁻ relativistic wind is launched from a rotating magnetosphere and how its energy is released at the termination shock.

Observations of gamma-ray binaries show flux variations tied to the orbital phase in most of them. Theoretical models have focused on relating these variations to changes in the line of sight and/or the conditions at the termination shock as the pulsar follows its eccentric orbit. Amongst the different systems, LS 5039 constitutes a useful testbed because of its regular, well-documented modulations in X rays (Takahashi et al. 2009), low energy (LE, 1 − 100 MeV, Collmar & Zhang 2014), high energy (HE, 0.1 − 100 GeV, Abdo et al. 2009) and very high energy (VHE , ≥ 100 GeV, Aharonian et al. 2006) gamma rays. The X-ray, LE and VHE modulations are in phase, with a peak at inferior conjunction (when the compact object passes in front of the massive star as seen by the observer) and a minimum close to superior conjunction (compact object behind the star). In contrast, the HE gamma-ray modulation is in anti-phase, peaking at superior conjunction. LS 5039’s short orbital period of 3.9 days has made it easier to establish these modulations than in the other systems. The modulations appear stable over time, with no reported change in the orbital lightcurves.

Nearly all models invoke synchrotron and inverse Compton emission from pairs with energies up to several TeV to explain the high-energy emission from LS 5039. These processes are efficient in the sense that the energy loss timescale is usually short compared to the flow timescale for the most energetic pairs: a large fraction of the available power may thus end up as high-energy radiation (e.g. Bosch-Ramon & Khangulyan 2009).

The seed photons for inverse Compton emission are provided by the massive star. The pairs upscatter these stellar UV photons to HE and VHE gamma rays. The gamma-ray emission is at its maximum when pairs backscatter the stellar photons in the direction of the observer, at superior conjunction. However, the gamma-ray photons must propagate through the binary system before reaching the observer. VHE photons are above the threshold for pair production with stellar photons (≳ 30 GeV when interacting with the kT_⋆ ≈ 3 eV blackbody photons from the star), so a large part of the VHE flux can be lost to creating e⁺e⁻ pairs before the radiation escapes. The γγ opacity in LS 5039 is at the minimum close to inferior conjunction and maximum close to superior conjunction, explaining why HE and VHE gamma rays are anti-correlated, although both arise from inverse Compton emission of the same seed photons (see Sect. 4 in Dubus 2013 and references therein).

There are two caveats to this interpretation of the HE and VHE gamma-ray modulation. First, a straightforward application predicts that none of the VHE flux emitted in the vicinity of the compact object should make it through the system at superior conjunction, whereas observations still detect a faint source. One explanation is that there is a contribution to the emission from the electromagnetic cascade that occurs as newly created e⁺e⁻ pairs emit VHE photons that in turn create pairs, etc. (e.g. Bednarek 2006; Bosch-Ramon et al. 2008; Cerutti et al. 2010). Another, non-exclusive, explanation is that the VHE emission arises from a larger or more distant region, diluting the effets of the γγ opacity (e.g. Zabalza et al. 2013). The second caveat is that the HE and VHE spectra are very distinct, with the HE spectrum cutting off exponentially at a few GeV (Hadasch et al. 2012). Clearly, different populations of electrons must be involved in the HE and VHE domains. Their origin is uncertain. Possible sites that have been considered include the pulsar magnetosphere, the pulsar wind, various locations along the pulsar wind termination shock, and the stellar wind termination shock.

The interpretation of the X-ray and LE gamma-ray modulation also requires additional ingredients. Synchrotron emission dominates in this range. Its luminosity depends on the density of electrons and the magnetic field B, both of which can be affected by the changing distance of the termination shock to the pulsar. Indeed, with an eccentricity e = 0.35 (Casares et al. 2005), the orbital separation d in LS 5039 varies from 0.1 AU at periastron to 0.2 AU at apastron. If the winds are isotropic and have a constant velocity, the distance R_s to the termination shock is set by (e.g. Lebedev & Myasnikov 1990) $\frac{{\mathit{R}}_{\mathrm{s}}}{\mathit{d}}\mathrm{=}\frac{{\mathit{\eta }}^{\mathrm{1}\mathit{/}\mathrm{2}}}{\mathrm{1}\mathrm{+}{\mathit{\eta }}^{\mathrm{1}\mathit{/}\mathrm{2}}}\mathit{,}$(1)where $\mathit{\eta }\mathrm{=}\mathit{Ė}\mathit{/}\mathit{Ṁ}{}_{\mathrm{w}}\mathit{v}_{\mathrm{w}}\mathit{c}$, with Ė the pulsar spindown power, Ṁ_w and v_w the stellar wind mass loss rate and velocity. The termination shock distance R_s doubles from periastron to apastron, decreasing B ∝ 1 /R_s proportionally (e.g. Dubus 2006b; Takata & Taam 2009). The changing separation can also affect adiabatic cooling of the particles as they are advected away from their acceleration site (Takahashi et al. 2009). However, it is difficult to tie such changes to the observed X-ray modulation: the peaks and dips at conjunctions suggest that it is due to a geometrical effect related to the observer line of sight rather than due to intrinsic changes in the conditions experienced by the particles. One possibility related to the line of sight is Doppler boosting. If η is small, the termination shock has the appearance of a bow shock facing away from the massive star. At inferior conjunction the bow shock flows in the general direction of the observer, whereas it flows in the opposite direction at superior conjunction. If the shocked pulsar wind retains a moderately relativistic bulk motion, the synchrotron emission will be boosted at inferior conjunction and deboosted at superior conjunction. Dubus et al. (2010) have shown that the effect is significant enough to be able to explain the X-ray modulation.

A more accurate assessment of this scenario requires numerical simulations. The geometry of the termination shock, the bulk velocity of the shocked fluid at each location, and the adiabatic losses can all be derived from a relativistic hydrodynamical simulation instead of being parametrised as done in previous works. The impact of γγ absorption, Doppler boosting and particle cooling on the orbital lightcurve can then be quantified properly, tightening the constraints on the underlying particles and pulsar wind physics. Gamma-ray binaries have been simulated by several groups, focusing on the geometry of the interaction region. Relativistic hydrodynamical simulations indicate that the material re-accelerates to very high Lorentz factors in the tail of the bow shock (Bogovalov et al. 2008) and that the non-zero thermal pressure in the stellar wind leads to smaller opening angles for the bow shock than usually assumed (Lamberts et al. 2013). Pulsar wind nebulae are thought to have a low magnetisation σ at the termination shock (Gaensler & Slane 2006), in which case the magnetic field has a negligible impact on the shocked flow, and indeed Bogovalov et al. (2012) found little difference, on the orbital separation scale, between their relativistic magneto-hydrodynamical (RMHD) and relativistic hydrodynamical (RHD) simulations of colliding winds in gamma-ray binaries. Other simulation work includes Bosch-Ramon et al. (2012), who studied mixing due to orbital motion using large scale 2D relativistic hydrodynamical simulations, Paredes-Fortuny et al. (2015) who looked at the impact of a clumpy stellar wind on the shock structure, and Takata et al. (2012) who presented two 3D SPH non-relativistic simulations of the interaction of a pulsar wind with a Be disc and wind (covering very large scales with limited spatial resolution). Takata et al. (2012) computed some lightcurves and spectra from their simulations, but do not take into account particle cooling and relativistic effects. A comprehensive approach linking high-resolution simulations of the shock region with particle cooling and emission (including relativity and anisotropic effects) has not been attempted yet.

Here, we use a relativistic hydrodynamical simulation as the basis for such a comprehensive model of the emission from particles in the shocked pulsar wind, with the aim of explaining the spectral energy distribution and the flux orbital modulations observed from LS 5039. The simulation (Sect. 2) provides the spatial evolution of the density, velocity and internal energy necessary to compute the extension of the emission, the impact of Doppler boosting, the importance of adiabatic cooling. This radiative post-processing step is described in Sect. 3. The general results are presented in Sect. 4, the more detailed application to LS 5039 in Sect. 5.

2. Relativistic hydrodynamics

We perform a 3D simulation of LS 5039 using the RAMSES hydrodynamical code (Teyssier 2002). The complete description of the extension to special relativity is in Lamberts et al. (2013). Here, we recall the major aspects of the method. RAMSES uses a Cartesian grid and allows adaptive mesh refinement (AMR) so as to locally increase the resolution at a reasonable computational cost. It solves the 3D-RHD equations using an upwind second order Godunov method. These equations can be written, for an ideal fluid, as a system of conservation equations in the laboratory frame (Landau & Lifshitz 1959) $\begin{array}{ccc}& & \frac{\mathit{\partial D}}{\mathit{\partial t}}\mathrm{+}\frac{\mathit{\partial }\mathrm{\left(}\mathit{D}{\mathit{v}}_{\mathit{k}}\mathrm{\right)}}{\mathit{\partial }{\mathit{x}}_{\mathit{k}}}\mathrm{=}\mathrm{0}\\ & & \frac{\mathit{\partial }{\mathit{M}}_{\mathit{j}}}{\mathit{\partial t}}\mathrm{+}\frac{\mathit{\partial }\mathrm{(}{\mathit{M}}_{\mathit{j}}{\mathit{v}}_{\mathit{k}}\mathrm{+}\mathit{p}{\mathit{\delta }}^{\mathit{j,k}}\mathrm{)}}{\mathit{\partial }{\mathit{x}}_{\mathit{k}}}\mathrm{=}\mathrm{0}\\ & & \frac{\mathit{\partial E}}{\mathit{\partial t}}\mathrm{+}\frac{\mathit{\partial }\mathrm{(}\mathit{E}\mathrm{+}\mathit{p}\mathrm{)}{\mathit{v}}_{\mathit{k}}}{\mathit{\partial }{\mathit{x}}_{\mathit{k}}}\mathrm{=}\mathrm{0}\mathit{,}\end{array}$where the vector of conservative variables is given by ${U}\mathrm{=}\left(\begin{array}{c}\\ \mathit{D}\\ {\mathit{M}}_{\mathit{j}}\\ \mathit{E}\end{array}\right)\mathrm{=}\left(\begin{array}{c}\\ \mathrm{\Gamma }\mathit{\rho }\\ {\mathrm{\Gamma }}^{\mathrm{2}}\mathit{\rho h}{\mathit{v}}_{\mathit{j}}\\ {\mathrm{\Gamma }}^{\mathrm{2}}\mathit{\rho }{\mathit{c}}^{\mathrm{2}}\mathit{h}\mathrm{-}\mathit{p}\end{array}\right)\mathit{.}$(5)D is the density, M the momentum density and E the energy density in the frame of the laboratory. c is the speed of light, the subscripts j,k stand for the dimensions, δ^jk is the Kronecker symbol. h is the specific enthalpy, ρ is the proper mass density, v_j is the fluid three velocity, p is the gas pressure. The gravitational force is neglected since we do not treat the acceleration of the winds and the outflow speeds are well above the escape velocity of the system. The Lorentz factor is given by $\mathrm{\Gamma }\mathrm{=}\frac{\mathrm{1}}{\sqrt{\mathrm{1}\mathrm{-}{\mathit{v}}^{\mathrm{2}}\mathit{/}{\mathit{c}}^{\mathrm{2}}}}\mathrm{·}$(6)The set of Eqs. (2)–(4) is completed by an equation of state that describes the thermodynamics of the fluid. In RHD, the ideal gas approximation is inconsistent with the kinetic theory of relativistic gases (Taub 1948). The exact equation of state for relativistic fluids involves modified Bessel functions, which lead to additional computational costs. Therefore we use the approximation developed by Mignone et al. (2005), which differs from the exact solution by less than 4%. In this case, the specific enthalpy is given by $\mathit{h}\mathrm{=}\frac{\mathrm{5}}{\mathrm{2}}\frac{\mathit{p}}{\mathit{\rho }{\mathit{c}}^{\mathrm{2}}}\mathrm{+}\sqrt{\frac{\mathrm{9}}{\mathrm{4}}{\left(\frac{\mathit{p}}{\mathit{\rho }{\mathrm{c}}^{\mathrm{2}}}\right)}^{\mathrm{2}}\mathrm{+}\mathrm{1}}\mathit{.}$(7)With the specific internal energy defined as , this gives a variable equivalent adiabatic index of $\mathit{\gamma ̂}\mathrm{=}\frac{\mathit{h}\mathrm{-}\mathrm{1}}{\mathit{h}\mathrm{-}\mathrm{1}\mathrm{-}\frac{\mathit{p}}{\mathit{\rho }{\mathrm{c}}^{\mathrm{2}}}}\mathrm{·}$(8)In our simulation, we use a HLL solver with a second order reconstruction based on a minmod slope limiter. This solver limits the development of Kelvin-Helmholtz instabilities in the simulation (Lamberts et al. 2013), a deliberate choice that we justify in Sect. 4.1. The size of our cubic simulation box l_box is six times the binary separation d. With z the binary axis, the box extends along the binary axis from z = −1.5 to z = 4.5 in a coordinate system scaled by the orbital separation d. The star is at the origin and the pulsar at x = 0,y = 0,z = 1. The x-axis extends from −3 to 3 and the y-axis extends from y = −1.5 to y = 4.5, so the binary is not located at an edge of the simulation box. Our coarse grid is made of 128 cells and we use four levels of refinement. This gives an equivalent resolution of 2048³ cells. The high resolution zone, where the four levels of refinement are effectively used, covers a slab of width Δy = Δz = 2 centred on the pulsar and Δz = 1 centred on x = 0. Refinement is based on density and Lorentz factor gradients.

We initialise the winds in a spherical region with a continuous outflow. The spherical regions both have a radius of 0.15 (in units of d). This is large enough to yield spherical symmetry, but small enough to avoid an impact on the formation of the shock. The winds are updated in the spherical region at every timestep. We assume the winds have reached terminal velocity at launch. A constant low-density medium initially fills the simulation box. This medium is cleared out as the winds propagate out and interact. After a transition phase corresponding to about 10 t_dyn (where we define t_dyn as the time for the pulsar wind to reach the back shock located at x ≈ 3), the simulation converges to a laminar, stationary state where the position of the shocks does not evolve with time anymore i.e. time drops out of Eq. (2)−(4). Our radiative calculation is based on a snapshot in the (y,z) plane of the hydrodynamics in this state. The geometry of the interaction depends only on the ratio of wind momentum fluxes η (Sect. 1). For LS 5039, η is thought to be ≃ 0.1, implying a reasonable value for the pulsar spindown power Ė ≈ 4 × 10³⁶ erg s^-1 assuming Ṁ_w ≈ 10^-7M_⊙ yr^-1, v_w ≈ 2000 km s^-1 (Szostek & Dubus 2011; Zabalza et al. 2011).

In the simulation, the stellar wind has a mass loss rate of Ṁ_w = 2 × 10^-8M_⊙ yr^-1 and a velocity of 2000 km s^-1. It has a Mach number ℳ = 30 at the location of the pulsar. Conventional models of pulsar winds invoke ultra-relativistic Lorentz factors up to Γ_p = 10⁶. Such high values are well beyond the range that can be achieved by standard fluid methods. We set the velocity of the pulsar wind to v_p = 0.99c, or Γ_p = 7.08. This is high enough to capture the relativistic effects in the shocked flow, especially near the apex where the shock is perpendicular and the post-shock flow speed tends to when Γ_p ≫ 1. In the wings, the shocked flow gradually accelerates to a velocity close to the initial pre-shock velocity. The spatial scale on which this occurs increases with the initial Lorentz factor of the wind (Lamberts et al. 2013). The pulsar mass loss rate is scaled to have η = 0.1, so that $\dot{{\mathit{M}}_{\mathrm{p}}}\mathrm{=}\mathit{\eta }\frac{\mathit{Ṁ}\mathrm{w}{\mathit{v}}_{\mathrm{w}}}{{\mathrm{\Gamma }}_{\mathrm{p}}{\mathit{v}}_{\mathrm{p}}}\mathrm{·}$(9)The flow dynamics and geometry will be correct even if the density must be scaled. For reference, the pulsar wind power corresponding to our choice of Ṁ_w = 2 × 10^-8M_⊙ yr^-1 and η = 0.1 in the simulation is $\mathit{Ė}\mathrm{=}{\mathrm{\Gamma }}_{\mathrm{p}}\mathit{Ṁ}{}_{\mathrm{p}}\mathit{c}^{\mathrm{2}}\mathrm{\approx }\mathrm{7.6}\mathrm{\times }{\mathrm{10}}^{\mathrm{35}}\hspace{0.17em}\mathrm{erg}\hspace{0.17em}{\mathrm{s}}^{-1}$. However, note that we can scale the spindown power up or down without changing the hydrodynamical simulation as long as Ṁ_w is scaled in proportion to keep η = 0.1. The classical Mach number of the pulsar wind is set to 20, which gives a relativistic Mach number ℳ_rel = ℳΓ_p ≃ 140. Orbital motion turns the shocked structure into a spiral of stepsize of order S ≳ v_sP_orb ≃ 4 AU ≃ 20 d (Lamberts et al. 2012). Our simulation covers a smaller size (≃ 5 d) and we make the assumption that we can neglect orbital motion on this scale¹. Even if previous simulations indicate that this is a reasonable approximation when the weaker wind is also the fastest (see e.g. Fig. 6 in Lamberts et al. 2012 or Fig. 1 in Bosch-Ramon et al. 2012), a more realistic model should include the orbital motion. We discuss in Sect. 6.2 how this might change our results.

The ratio η does not change along the orbit since the winds are isotropic and are at their terminal velocity. The dynamical timescale corresponds to t_dyn ≈ 300 s for the typical orbital separation of 0.2 AU appropriate to LS 5039 (see above). The time taken to form the interaction region, on the scales that we consider, is very short compared to the orbital timescale of 3.9 days for LS 5039. The structure has ample time to reach a stationary state at each orbital phase. Hence, since the stationary structure only depends on η, which is constant with orbital phase, one simulation is enough for all orbital phases even if the orbital separation changes. The velocities, Ṁ and wind temperature (Mach number) remain constant. However, for each orbital phase we rescale the distances to the actual orbital separation i.e. by a factor d. The density and pressure in a given pixel are thus rescaled by a factor 1 /d².

By doing a RHD simulation instead of a RMHD simulation, we implicitly assume that the magnetic field has no dynamical role in the interaction on the spatial scales probed by the simulation. This assumption is supported by the study of Bogovalov et al. (2012), who carried out both relativistic hydrodynamical and relativistic MHD simulations of interacting winds, on spatial scales comparable to those studied here, and found little difference between the two when the pulsar magnetisation σ ≤ 0.1. Conventional models of pulsar wind nebulae assume σ ≪ 1 at the termination shock to ensure the efficient conversion of the wind kinetic energy to particle energy at the shock (Kennel & Coroniti 1984; Gaensler & Slane 2006; Bucciantini et al. 2011). We come back to this issue in Sect. 6.2.

3. High-energy radiation

Fig. 1 Model geometry showing our choices of axis and angles. The interaction has cylindrical symmetry around the z axis. The two thick solid lines in the (x,z) plane represent the pulsar wind termination shock and the contact discontinuity with the shocked stellar wind. Dashed lines represent the same in the (y,z) plane. Particles injected at the termination shock follow streamlines in the region in between these two surfaces.

The interaction of stellar winds classically leads to a double shock structure separated by a contact discontinuity (Fig. 1). Non-thermal emission may be expected from the shock associated with both winds (Bednarek 2011). There is plenty of evidence for efficient high-energy non-thermal emission from pulsar wind termination shocks, but rather less from stellar wind termination shocks with only Eta Carinae detected so far in gamma rays (Werner et al. 2013). Here, we have assumed that the non-thermal emission results exclusively from particles accelerated at the pulsar wind termination shock. We have not taken into account thermal bremsstrahlung from the shocked stellar wind, since this would be best treated by dedicated simulations à la Stevens et al. (1992), and because the X-ray emission from LS 5039 appears entirely non-thermal with, as yet, no evidence for thermal emission (Zabalza et al. 2011).

Particles are randomised and/or accelerated to very high energies at the termination shock of the pulsar. We wish to follow the evolution of these particles once they are injected in the shocked pulsar wind flow. Once their energy distribution at each location is known, we can compute the radiation seen at a given line of sight to the binary.

We adopted some basic assumptions to make the problem tractable. First, the evolution of the particles is decoupled from the hydrodynamical simulation and treated in the test particle limit as a post-processing step i.e. radiative cooling does not impact the dynamics of the flow. We come back to this in Sects. 4.2 and 5.1. Second, the high-energy particles are assumed to follow the flow. Spatial diffusion is expected to remain negligible if the Larmor radius remains small compared to the flow spatial scales, which is the case here. Third, plasma processes are assumed to keep the particle momentum distribution isotropic. Last, the flow is assumed to be stationary, simplifying the problem to that of following the evolution of particles along streamlines.

We start by explaining how we calculate the total emission based on the particle evolution along streamlines (Sect. 3.1). We then describe how we choose the initial distribution of the injected particles (Sect. 3.2), how we compute the particle evolution and streamline emission (Sect. 3.3), how we estimate the magnetic field (Sect. 3.4), and how we deal with the changing aspect due to orbital motion (Sect. 3.5).

3.1. Total flux

The total flux from the emission region, measured in the laboratory frame at a frequency ν, is given by $\mathit{F}\mathrm{\left(}\mathit{\nu }\mathrm{\right)}\mathrm{=}\frac{\mathrm{1}}{{\mathit{D}}^{\mathrm{2}}}{\mathrm{\int }}_{\mathit{V}}{\mathrm{𝒟}}_{\mathrm{obs}}^{\mathrm{2}}\mathit{j}\mathrm{(}\mathit{\nu }\mathit{/}{{\mathrm{𝒟}}_{\mathrm{obs}}}^{\mathrm{)}}\exp \mathrm{(}\mathrm{-}{\mathit{\tau }}_{\mathit{\nu }}\mathrm{)}\mathrm{d}\mathit{V,}$(10)where the integral is over the volume V of the region (measured in the laboratory frame), D is the distance to the source, j is the local particle emissivity in the co-moving frame (j is the sum of the synchrotron j_sync and inverse Compton j_ic emissivities), and τ_ν is the opacity at frequency ν due to pair production as the photons emitted in the volume dV travel along the line of sight to the observer (see Sect. 3.3). is the Doppler boost to the observer ${\mathrm{𝒟}}_{\mathrm{obs}}\mathrm{=}{\left[\mathrm{\Gamma }\mathrm{(}\mathrm{1}\mathrm{-}{v}\mathit{.}{{{e}}_{\mathrm{obs}}}^{\mathrm{)}}\right]}^{-1}\mathit{,}$(11)with v the flow velocity vector and e_obs the unit vector in the direction of the observer (Sect. 3.5).

The simulated flow is laminar. To compute this integral, we extract from a snapshot streamlines in the emission region. Each streamline starts at the position of the shock and ends where it leaves the computational domain. These streamlines subdivide the emission region into streamtubes². The integral over the volume can be written as an integral over time of the particles flowing along each streamtube: $\mathit{F}\mathrm{\left(}\mathit{\nu }\mathrm{\right)}\mathrm{=}\frac{\mathrm{1}}{{\mathit{D}}^{\mathrm{2}}}\sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathrm{𝒩}}{\mathrm{\int }}_{\mathrm{0}}^{{\mathit{t}}_{\mathit{i}}}{\mathrm{𝒟}}_{\mathrm{obs}}^{\mathrm{2}}\mathit{j}\mathrm{(}\mathit{\nu }\mathit{/}{\mathrm{𝒟}}_{\mathrm{obs}}\mathrm{)}\exp \mathrm{(}\mathrm{-}{\mathit{\tau }}_{\mathit{\nu }}\mathrm{)}{v}\mathit{.}{S}\mathrm{d}\mathit{t,}$(12)where S is a cross section of the streamtube i and t = t_i is the time taken by particles flowing along streamline i to leave the computational domain, with t = 0 at their injection at the shock. Since the flow is stationary, the particle flux along each streamtube $\mathit{Ṅ}\mathit{i}\mathrm{\equiv }{\mathrm{\Gamma }\mathit{n}{v}\mathit{.}{S}|}_{\mathit{i}}$(13)is constant (n = ρ/m_e in the pulsar wind composed of e⁺e⁻ pairs). The total flux from the shocked region becomes $\mathit{F}\mathrm{\left(}\mathit{\nu }\mathrm{\right)}\mathrm{=}\frac{\mathrm{1}}{{\mathit{D}}^{\mathrm{2}}}\sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathrm{𝒩}}{\mathit{f}}_{\mathit{i}}\mathrm{\left(}\mathit{\nu }\mathrm{\right)}\mathit{Ṅ}\mathit{i}\mathit{,}$(14)where we have defined the fluence per particle f_i of the streamline i as ${\mathit{f}}_{\mathit{i}}\mathrm{\left(}\mathit{\nu }\mathrm{\right)}\mathrm{\equiv }{\mathrm{\int }}_{\mathrm{0}}^{{\mathit{t}}_{\mathit{i}}}\frac{{\mathrm{𝒟}}_{\mathrm{obs}}^{\mathrm{2}}\mathit{j}\mathrm{(}\mathit{\nu }\mathit{/}{{\mathrm{𝒟}}_{\mathrm{obs}}}^{\mathrm{)}}\exp \mathrm{(}\mathrm{-}{\mathit{\tau }}_{\mathit{\nu }}\mathrm{)}}{\mathrm{\Gamma }\mathit{n}}\mathrm{d}\mathit{t}\mathit{.}$(15)This is an integral over the emission per particle as they flow along their streamline. The integral can also be recast as an integral over position along the streamline since the curvilinear distance l is related to time by dl = vdt. We obtain Ṅ_i numerically by measuring the particle flux across the shock surface for each streamtube in the RHD simulation (Eq. (13)).

We make the implicit assumption in Eq. (14) that the emissivity j and the opacity τ_ν do not vary substantially across the streamtube for a given position. In practice, taking a sufficient number of streamlines ensures this is achieved. We used streamlines, evenly spread along the shock surface, and we checked that this is more than enough for numerical convergence of Eq. (14).

Alternatively, one could divide up the particle distribution into energy bins, recast the evolution equation of the particles in each energy bin in conservative form and add them to the set solved by RAMSES (e.g. Reitberger et al. 2014). The advantage is that this can deal with non-stationary flows and mixing. However, a major drawback is that in our case the particle cooling time at very high energies can be very short: the particles cool on small spatial scales, requiring a very high spatial resolution (see Sect. 4.2 below). Another difficulty is that radiative cooling does not scale with the orbital separation d as adiabatic cooling (d^-2 and d^-1 respectively, see Eq. (23) below). Each orbital phase would thus require a full simulation to account for this difference in cooling spatial length. Each new choice of parameters for radiative cooling would also require a new simulation run. Treating the particle evolution in a post-processing step makes the problem much more tractable, allowing for a wider exploration of the parameter space involved in the flow emission.

We explain our choices for particle injection at the shock first before turning to the computation of the streamline fluence (Eq. (15)) in Sect. 3.3.

3.2. Particle injection at the shock

We distribute the particles as a power-law function of their Lorentz factor γ: ${\frac{\mathrm{d}\mathit{n}}{\mathrm{d}\mathit{\gamma }}|}_{\mathit{t}\hspace{0.17em}\mathrm{=}\hspace{0.17em}\mathrm{0}}\mathrm{=}\mathit{K}{\mathit{\gamma }}^{\mathrm{-}\mathit{s}}\mathrm{with}{\mathit{\gamma }}_{\min }\mathrm{\le }\mathit{\gamma }\mathrm{\le }{\mathit{\gamma }}_{\max }\mathit{.}$(16)The particle density at the shock sets the normalisation K of the distribution. The two other parameters are γ_min and γ_max.

The maximum Lorentz factor γ_max is set by the balance between acceleration and radiative losses, since the gyroradius will be much smaller than the characteristic size of the acceleration region in our case (Dubus 2006b). We assume that the acceleration timescale in the comoving frame τ_acc is some multiple of the Böhm limit: ${\mathit{\tau }}_{\mathrm{acc}}\mathrm{=}\mathrm{2}\mathit{\pi \xi }{\mathit{R}}_{\mathrm{L}}\mathit{/}\mathit{c,}$(17)with R_L the gyroradius. We expect ξ ≥ 1 for diffusive shock acceleration, with ξ ≤ 10 corresponding to “extreme acceleration” (Khangulyan et al. 2008); ξ< 1 may be possible for acceleration at magnetic reconnection sites, as proposed to explain the gamma-ray flares from the Crab nebula (Cerutti et al. 2012). Synchrotron losses dominate over inverse Compton losses at very high energies (Dubus 2006b). The synchrotron loss timescale is ${\mathit{\tau }}_{\mathrm{sync}}\mathrm{\equiv }\mathit{\gamma }{\left(\frac{\mathrm{d}\mathit{\gamma }}{\mathrm{d}{\mathit{t}}^{\mathrm{\prime }}}\right)}^{-1}\mathrm{=}\frac{\mathit{\gamma }{\mathit{m}}_{\mathrm{e}}{\mathit{c}}^{\mathrm{2}}}{\frac{\mathrm{4}}{\mathrm{3}}{\mathit{\sigma }}_{\mathit{T}}\mathit{c}\mathrm{(}\mathit{\beta \gamma }{\mathrm{)}}^{\mathrm{2}}{\mathit{u}}_{\mathrm{b}}}\mathrm{\approx }\mathrm{77}\hspace{0.17em}{\left(\frac{\mathrm{1}\hspace{0.17em}\mathrm{G}}{\mathit{b}}\right)}^{\mathrm{2}}\left(\frac{{\mathrm{10}}^{\mathrm{7}}}{\mathit{\gamma }}\right)\mathrm{s}\mathit{,}$(18)where b and t′ are the magnetic field and time in the comoving frame. Setting τ_acc ≤ τ_sync gives ${\mathit{\gamma }}_{\max }\mathrm{=}{\left(\frac{\mathrm{3}\mathit{e}}{\mathit{\xi }{\mathit{\sigma }}_{\mathit{T}}\mathit{b}}\right)}^{\mathrm{1}\mathit{/}\mathrm{2}}\mathrm{\approx }\mathrm{5}\mathrm{\times }{\mathrm{10}}^{\mathrm{7}}{\mathit{\xi }}^{\mathrm{-}\mathrm{1}\mathit{/}\mathrm{2}}{\left(\frac{\mathrm{1}\hspace{0.17em}\mathrm{G}}{\mathit{b}}\right)}^{\mathrm{1}\mathit{/}\mathrm{2}}\mathrm{·}$(19)With γ_max and K known, we derive γ_min by writing that the energy in the particle distribution is a fraction ζ_p of the downstream specific energy ϵ: ${\mathit{\zeta }}_{\mathrm{p}}\mathit{ϵ}\mathrm{\equiv }{\mathit{ϵ}}_{\mathrm{nt}}\mathrm{=}\left({\mathrm{\int }}_{{\mathit{\gamma }}_{\min }}^{{\mathit{\gamma }}_{\max }}\mathit{\gamma }{\mathit{m}}_{\mathrm{e}}{\mathit{c}}^{\mathrm{2}}\frac{\mathrm{d}\mathit{n}}{\mathrm{d}\mathit{\gamma }}\mathrm{d}\mathit{\gamma }\right){}^{\mathrm{/}}{\left({\mathrm{\int }}_{{\mathit{\gamma }}_{\min }}^{{\mathit{\gamma }}_{\max }}\frac{\mathrm{d}\mathit{n}}{\mathrm{d}\mathit{\gamma }}\mathrm{d}\mathit{\gamma }\right)}^{}\mathrm{·}$(20)The numerator is the total energy in the distribution and the denominator is the particle density. This equation implicitly sets γ_min, which can be found numerically with a Newton-Raphson scheme. If γ_max ≫ γ_min and s> 2 then the integrals can be simplified and solved for γ_min, such that ${\mathit{\gamma }}_{\min }\mathrm{\approx }{\mathit{\zeta }}_{\mathrm{p}}\mathit{ϵ}\left(\frac{\mathit{s}\mathrm{-}\mathrm{2}}{\mathit{s}\mathrm{-}\mathrm{1}}\right)\mathrm{=}\frac{{\mathit{\zeta }}_{\mathrm{p}}}{\mathrm{(}\mathit{\gamma ̂}\mathrm{-}\mathrm{1}\mathrm{)}}\left(\frac{\mathit{s}\mathrm{-}\mathrm{2}}{\mathit{s}\mathrm{-}\mathrm{1}}\right)\frac{\mathit{p}}{\mathit{n}{\mathit{m}}_{\mathrm{e}}{\mathit{c}}^{\mathrm{2}}}\mathrm{·}$(21)ζ_p enables us to scale the emission to realistic values since the (lab-frame) specific enthalpy of the cold pulsar wind at the shock is Γ_p ≈ 7 in the simulation, which is too low to reproduce the very high values of ϵ (equivalently Γ_p) required to fit the observations. Raising ζ_p implies that the effective mass loss rate from the pulsar wind is decreased by ${\mathit{\zeta }}_{\mathrm{p}}^{-1}$ because Ė/c = ζ_pΓ_pṀ_pc is fixed by η for given Ṁ_w, v_w, Γ_p (in other words, the same total available energy is distributed amongst fewer particles).

We have also experimented with a relativistic Maxwellian distribution because Fermi-LAT observations of gamma-ray binaries in HE gamma rays require an additional population of particles with a narrow range in energy (Sect. 1). Moreover, particle-in-cell simulations of relativistically-shocked pair plasmas typically show a prominent Maxwellian distribution of shock-heated particles together with the power-law distribution of accelerated particles (Sironi & Spitkovsky 2009, 2011). The relativistic Maxwellian distribution is ${\frac{\mathrm{d}\mathit{n}}{\mathrm{d}\mathit{\gamma }}|}_{\mathit{t}\hspace{0.17em}\mathrm{=}\hspace{0.17em}\mathrm{0}}\mathrm{\equiv }\mathit{K}{\mathit{\gamma }}^{\mathrm{2}}\exp \mathrm{(}\mathrm{-}\mathit{\gamma }\mathit{/}{\mathit{\gamma }}_{\mathrm{t}}\mathrm{)}\mathit{.}$(22)Again, K is derived by imposing that the integral of the distribution scales with the particle density at the shock. The mean Lorentz factor of the distribution γ_t is derived from Eq. (20), which gives γ_t ≈ e_nt/ 3 when γ_min ≪ γ_t ≪ γ_max.

We assumed that ζ_p and ξ do not vary along the shock for lack of strong justifications for more sophisticated assumptions.

3.3. Streamline emission

For a stationary flow, the evolution of a particle injected at a given location and followed along the associated streamline in the comoving frame, is uniquely set by the evolution equation (e.g. Del Zanna et al. 2006; Mimica et al. 2009; Porth et al. 2014, for applications to pulsar wind nebulae or AGN jets) $\frac{\mathrm{1}}{\mathit{\gamma }}\frac{\mathrm{d}\mathit{\gamma }}{\mathrm{d}{\mathit{t}}^{\mathrm{\prime }}}\mathrm{=}\frac{\mathrm{d}\ln \mathit{ϵ}}{\mathrm{d}{\mathit{t}}^{\mathrm{\prime }}}\mathrm{-}\frac{\mathrm{1}}{{\mathit{\tau }}_{\mathrm{sync}}}\mathrm{-}\frac{\mathrm{1}}{{\mathit{\tau }}_{\mathrm{ic}}}\mathit{,}$(23)where the terms on the right hand side represent adiabatic, synchrotron and inverse Compton losses – the only relevant cooling processes here. The elapsed time in the laboratory frame dt is related to the proper time in the comoving frame by dt = Γdt′. The bulk Lorentz factor Γ and the elapsed time dt are derived from the position and velocity along the streamline as calculated with RAMSES. The adiabatic loss term is also directly derived from the simulation.

The evolution of the magnetic field along the streamline must be known to compute the synchrotron losses τ_sync (Eq. (18)) and the associated emissivity j_sync. Our choices for the magnetic field are explained below in Sect. 3.4. The synchrotron emissivity j_sync is computed using the usual formula involving Bessel functions (Rybicki & Lightman 1979).

For the inverse Compton losses τ_ic, we take the massive star as the only source of seed photons and compute the electron energy losses τ_ic using the Jones (1968) scattering kernel. The star is modelled as a blackbody of temperature T_⋆ = 39 000 K and radius R_⋆ = 9.3 R_⊙. The density n_⋆ of photons with an energy ϵ_⋆ (in units of m_ec²) seen in the comoving frame by pairs at a distance d_⋆ from the star is, in photons per cm³ per unit energy, ${\mathit{n}}_{\mathit{\star }}\mathrm{=}\mathrm{2}\mathit{\pi }{\left(\frac{\mathit{mc}}{\mathit{h}}\right)}^{\mathrm{3}}{\left(\frac{{\mathit{R}}_{\mathit{\star }}}{{\mathit{d}}_{\mathit{\star }}}\right)}^{\mathrm{2}}\frac{{\mathit{ϵ}}_{\mathit{\star }}^{\mathrm{2}}\mathit{/}{\mathrm{𝒟}}_{\mathit{\star }}^{\mathrm{2}}}{\exp \left(\frac{{\mathit{ϵ}}_{\mathit{\star }}{\mathit{m}}_{\mathrm{e}}{\mathit{c}}^{\mathrm{2}}}{{\mathrm{𝒟}}_{\mathit{\star }}\mathit{k}{\mathit{T}}_{\mathit{\star }}}\right)\mathrm{-}\mathrm{1}}\mathrm{·}$(24) is the Doppler boost required to transform the stellar radiation field into the comoving frame, ${\mathrm{𝒟}}_{\mathit{\star }}\mathrm{=}{\left[\mathrm{\Gamma }\left(\mathrm{1}\mathrm{-}{v}\mathit{.}{{e}}_{{\star }}\right)\right]}^{-1}\mathit{,}$(25)with e_⋆ the unit vector giving the direction from the star to the flow element containing the pairs.

The evolution of the particle distribution can be calculated semi-analytically when inverse Compton losses are in the Thomson regime (Begelman & Li 1992). A numerical solution is required in our case since stellar photons are up-scattered to gamma-ray energies in the Klein-Nishina regime (Dubus et al. 2008). Following Bošnjak et al. (2009), the particle distribution of each streamline is discretised in “Lagrangian" bins ${\mathit{n}}_{\mathit{k}}\mathrm{=}{\mathrm{\int }}_{{\mathit{\gamma }}_{\mathit{k}}}^{{\mathit{\gamma }}_{\mathit{k}\mathrm{+}\mathrm{1}}}\frac{\mathrm{d}\mathit{n}}{\mathrm{d}\mathit{\gamma }}\mathrm{d}\mathit{\gamma }\mathrm{with}\sum _{\mathit{k}}{\mathit{n}}_{\mathit{k}}\mathrm{=}\mathit{n}\mathit{.}$(26)The relative number of particles n_k/n in each energy bin is conserved but the bin boundaries [γ_k,γ_{k + 1}] vary along the streamline according to the energy loss equation (Eq. (23)). We use 400 energy bins, initially logarithmically distributed between γ_min and γ_max. For each γ_k, Eq. (23) is integrated in time steps representing a small fraction (5%) of the minimum energy loss timescale. To ease the computational burden, we stop following the energy losses once γ becomes lower than 10³: we verified that these particles do not contribute emission in the frequency range that we are interested in.

We compute the streamline fluence f_i(ν) (Eq. (15)) once the evolution of the particle distribution along the streamline is known. The fluence depends on the location in the flow but also on the line of sight to the observer because of the relativistic boost associated with the bulk motion of the flow (Eq. (11)). Unlike the particle evolution calculation, which is symmetric around the binary axis and requires only a 2D integration (see Sect. 3.5), the emission calculation requires a full 3D integration since e_obs also varies with the azimuth θ around the binary axis (Fig. 1).

Synchrotron radiation is isotropic in the comoving frame for a random orientation of the magnetic field and an isotropic distribution of particles, both reasonable assumptions at this stage. However, the inverse Compton emission is clearly not isotropic in the comoving frame, even for an isotropic distribution, because the seed photons from the star are anisotropic³. We follow Dubus et al. (2010) to take this effect into account when computing the up-scattered emissivity j_ic towards the observer line of sight. Finally, we also calculate the γγ absorption of the VHE flux due to pair production with the stellar photons as in Dubus (2006a). The γγ opacity τ_γγ depends on the path of the VHE photon from its emission location to the observer. We neglected the finite size of the star to ease the computational requirements. This approximation appears reasonable here since the extension of the VHE emission region would probably reduce and smooth out any effect of the finite size on both the inverse Compton emission and the γγ opacity. Knowing τ_γ, j_ic, and j_sync along each streamline enables us to evaluate Eq. (15).

3.4. Magnetic field

Following the standard description of pulsar winds, the laboratory frame magnetic field B_p at a distance d_p from the pulsar in the unshocked wind is given by $\frac{{\mathit{B}}_{\mathrm{p}}^{\mathrm{2}}}{\mathrm{4}\mathit{\pi }}\left(\frac{\mathrm{1}\mathrm{+}\mathit{\sigma }}{\mathit{\sigma }}\right)\mathrm{=}\frac{\mathit{Ė}}{\mathrm{4}\mathit{\pi }{\mathit{d}}_{\mathrm{p}}^{\mathrm{2}}\mathit{c}}\mathit{,}$(27)where σ is the wind magnetisation $\mathit{\sigma }\mathrm{\equiv }\frac{{\mathit{B}}_{\mathrm{p}}^{\mathrm{2}}}{\mathrm{4}\mathit{\pi }{\mathrm{\Gamma }}_{\mathrm{p}}^{\mathrm{2}}{\mathit{n}}_{\mathrm{p}}{\mathit{m}}_{\mathrm{e}}{\mathit{v}}_{\mathrm{p}}^{\mathrm{2}}}\mathit{,}$(28)and the subscript p identifies a value in the unshocked wind. Pulsar winds are thought to have a low σ so we neglect the influence of the magnetic field on the flow dynamics (Sect. 2).

The magnetic field is purely toroidal at large distances from the magnetosphere of the pulsar, so the shock is perpendicular and the field is amplified at the shock by the compression ratio χ = B/B_p = N/N_p, where N = Γn is the laboratory frame number density. The compression ratio depends on the adiabatic index , which can change along the shock (see Eq. (7)). Even if is constant, the compression ratio changes in a relativistic shock when the transverse speed is non-negligible. For a ultra-high wind Lorentz factor, the shock behaves like a normal shock with except very far out in the wings, where the speed normal to the shock becomes non-relativistic and . We have taken everywhere along the shock, so that $\mathit{B}\mathrm{=}\frac{\mathrm{1}}{\mathit{\gamma ̂}\mathrm{-}\mathrm{1}}{\left[\frac{\mathit{Ė}}{{\mathit{d}}_{\mathrm{p}}^{\mathrm{2}}\mathit{c}}\left(\frac{\mathit{\sigma }}{\mathrm{1}\mathrm{+}\mathit{\sigma }}\right)\right]}^{\mathrm{1}\mathit{/}\mathrm{2}}$(29)In practice we take B = ζ_b/d_p, where ζ_b is a constant, since Ė and σ are not fixed but will be determined by the comparison to data.

Alternatively, the magnetic field could be amplified by plasma instabilities to a value that represents a fraction of the available internal energy $\frac{{\mathit{b}}^{\mathrm{2}}}{\mathrm{8}\mathit{\pi }}\mathrm{=}{\mathit{\zeta }}_{\mathrm{b}}\frac{\mathit{p}}{\mathit{\gamma ̂}\mathrm{-}\mathrm{1}}\mathit{,}$(30)where b = B/ Γ is the magnetic field in the comoving frame, again assumed to be toroidal, and (again) ζ_b is a constant. From the conservation of Γh across the shock, the pressure is related to the upstream conditions by

$\mathit{p}\mathrm{=}\frac{\mathit{\gamma ̂}\mathrm{-}\mathrm{1}}{\mathit{\gamma ̂}}{\left(\frac{{\mathrm{\Gamma }}_{\mathrm{p}}}{\mathrm{\Gamma }}\right)}^{\mathrm{2}}\frac{{\mathit{v}}_{\mathrm{p}}}{\mathit{v}}\left(\mathrm{1}\mathrm{-}\frac{\mathrm{\Gamma }}{{\mathrm{\Gamma }}_{\mathrm{p}}}\right){\mathit{n}}_{\mathrm{p}}{\mathit{m}}_{\mathrm{e}}{\mathit{c}}^{\mathrm{2}}$when the upstream pressure is negligible⁴. This equation shows that p ∝ (Γd_p)^-2, hence, that B ∝ 1 /d_p. Thus, there is little difference in practice between the magnetic field B as given by Eqs. (29) and (30), the two being exactly proportional when the pulsar wind Lorentz factor is high (Γ_p ≫ 1).

Since we assume no influence of the magnetic field on the flow, the evolution of B beyond the shock is set solely by the induction equation. For a stationary flow, with cylindrical symmetry around the binary axis, the induction equation for the purely toroidal B field becomes $\frac{\mathit{\partial }{B}}{\mathit{\partial t}}\mathrm{-}\mathrm{\nabla }\mathrm{\times }\mathrm{(}{v}\mathrm{\times }{B}\mathrm{)}\mathrm{=}\frac{\mathit{\partial }{\mathit{v}}_{\mathit{r}}\mathit{B}}{\mathit{\partial r}}\mathrm{+}\frac{\mathit{\partial }{\mathit{v}}_{\mathit{z}}\mathit{B}}{\mathit{\partial z}}\mathrm{=}\mathrm{0}\mathit{,}$(31)with z the coordinate along the symmetry axis and r the radial cylindrical coordinate (Fig. 1). The induction equation is identical to the continuity equation if B is replaced by Γnr (Micono et al. 1999; Bogovalov et al. 2012). Hence, the evolution of B along the flow streamline simply follows B ∝ Γnr, with the proportionality constant set by the initial conditions at the shock.

3.5. Geometrical factors and orbital motion

Taking the origin of the coordinate system at the centre of the massive star, z along the binary axis and r the radial coordinate perpendicular to the binary axis then the flow element is at (rcosθ,rsinθ,z) where θ is the azimuthal angle around the binary axis of symmetry (see Fig. 1). The unit vector e_⋆ is ${{e}}_{\mathit{\star }}\mathrm{=}\mathrm{\left(}\sin \mathit{\alpha }\cos \mathit{\theta ,}\sin \mathit{\alpha }\sin \mathit{\theta ,}\cos \mathit{\alpha }\mathrm{\right)}\mathit{,}$(32)with tanα = r/z. Taking advantage of the symmetry around the binary axis, the speed in the laboratory frame is v = (v_rcosθ,v_rsinθ,v_z). It is straightforward to see that the boost (Eq. (25)) that applies to the stellar emission seen in the comoving frame does not depend on θ: the evolution can be computed using a set of streamlines taken in a plane including the binary axis. The unit vector to the observer is ${{e}}_{\mathrm{obs}}\mathrm{=}\mathrm{\left(}\cos \mathit{i,}\sin \mathit{\omega }\sin \mathit{i,}\cos \mathit{\omega }\sin \mathit{i}\mathrm{\right)}\mathit{,}$(33)with i is the inclination of the system and ω the true anomaly of the orbit.

For each orbital phase, the simulation is scaled with the orbital separation d: $\mathit{d}\mathrm{=}\frac{\mathit{a}\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{e}}^{\mathrm{2}}\mathrm{)}}{\mathrm{1}\mathrm{-}\mathit{e}\sin \mathrm{(}\mathit{\omega }\mathrm{-}{\mathit{\omega }}_{\mathrm{p}}\mathrm{)}}\mathrm{·}$(34)The orbital parameters are those of LS 5039 (McSwain et al. 2004; Casares et al. 2005, 2011) i.e. the semi-major axis $\mathit{a}\mathrm{=}\mathrm{(}\mathit{GM}{\mathit{P}}_{\mathrm{orb}}^{\mathrm{2}}\mathit{/}\mathrm{4}{\mathit{\pi }}^{\mathrm{2}}{\mathrm{)}}^{\mathrm{1}\mathit{/}\mathrm{3}}$ with P_orb = 3.9 days, M = 1.4 M_⊙ + 23 M_⊙ the total mass, e = 0.35 the eccentricity, ω the true anomaly, and ω_p = 212° the binary angle at periastron. We divide the orbit into 30 phases, each of which requires a (2D) calculation of the particle evolution and a (3D) calculation of the observed emission. We use 20 cells for θ, 400 for the electron Lorentz factor γ, 20 for the stellar photon energy ϵ_⋆. The calculations were parallelised using OpenMP.

4. Results

4.1. The structure of the shocked flow

The numerical simulation shows the expected double shock structure. Numerical diffusivity, induced by our choice of Riemann solver, stabilises the structure despite the presence of a strong velocity shear at the interface between the shocked pulsar and stellar winds. The diffusivity leads to gradual mixing between the winds i.e. numerical spreading of the contact discontinuity, quenching the development of the Kelvin-Helmholtz (KH) instability (Lamberts et al. 2011). Strong KH mixing could impact the emission of the region, for instance by reducing the Lorentz factor of the flow, and by generating strong turbulence. The fluctuation timescale of the interface would be short since the flow is relativistic. However, the strong velocity shear is accompanied by a strong density contrast between the dense stellar wind and the tenuous pulsar wind. The ratio of the KH growth timescale to the advection timescale is ∝ (ρ₂/ρ₁)^1/2Δv for two fluids of density ρ₁ and ρ₂ ≪ ρ₁, sheared by a velocity difference Δv (see appendix in Lamberts et al. 2012 and Bodo et al. 2004 for the growth rate in the relativistic regime). Hence, the KH growth is dampened for high density contrasts, such as that expected between the tenuous highly relativistic pulsar wind and the dense stellar wind, making it debatable whether KH-induced mixing is dynamically important in gamma-ray binaries on the scales that we consider here. Bosch-Ramon et al. (2012, 2015) find mixing occurs mostly on larger scales in their simulations and attribute it rather to instabilities triggered by orbital motion. Large, dense clumps in the stellar wind could also affect the shock structure and variability (Paredes-Fortuny et al. 2015), although it is unclear whether there is enough time for the clumps to grow before reaching the termination shock in LS 5039 (located within 1−2 stellar radii of the star). Our simulation implicitly assumes that any mixing is limited and roughly captured – in a time-averaged sense – by the numerical diffusivity.

The basic structure of the shocked pulsar wind is illustrated in Fig. 2, where we show only part of the full simulation domain. The stellar wind (shocked and unshocked) and the unshocked pulsar wind have been edited out of this map as well as subsequent ones since our focus is entirely on the shocked pulsar wind. The head of the pulsar wind is shocked in a bow-shaped region with asymptotic angles ≈ 40° (termination shock) to 50° (contact discontinuity), measured from the z axis. This is larger than the 30° to 45° found by Bogovalov et al. (2008) for η = 0.1, but our simulation domain is smaller and our angles may not yet have reached their true asymptotic values.

Fig. 2 RHD flow in the shocked pulsar wind. The star is at the origin (0,0) and the pulsar is at (1,0) in units of the orbital separation d. A few selected streamlines have been plotted. Three different regions have been coloured. They correspond to the bow shock (blue), the reflected shock (light blue), and the back shock (dark blue) regions.

Besides the bow-shaped shock, our simulation shows that the pulsar is also terminated at the back instead of propagating freely. The structure is the one expected for Mach reflection on the binary axis as described in Bogovalov et al. (2008) in the context of gamma-ray binaries and as observed in e.g. simulations of pulsar bow-shock nebulae (Gaensler et al. 2004). Material flowing in the bow shock region abruptly changes direction when it crosses into the light blue region (black streamlines in Fig. 2). The change is due to a reflection shock that appears in order to accommodate the back shock region. This reflected shock region is separated by a contact discontinuity from the back shock region (boundary between medium and dark blue regions in Fig. 2).

The back shock structure in our simulation is very similar to the back shock structure in the 2D and 3D relativistic simulations of Bosch-Ramon et al. (2012, 2015), who use a different code (PLUTO) but similar values for η and Γ_p (η = 0.1 and Γ_p = 2 for their 3D simulation). All their simulations include orbital motion and they interpret the presence of this structure as an effect of orbital motion. Since this cannot be the case in our simulation, we suspect that the confinement depends on a subtle combination of 3D + relativistic + pressure (Mach number) effects (Lamberts et al. 2011, 2012, 2013). We defer a resolution of this possible issue to future studies. In the present case, as we shall see, the back shock and reflected shock only play a minor role in shaping the high-energy emission.

Fig. 3 Maps of various quantities in the shocked pulsar wind. The particle density n, pressure p, and magnetic field are displayed on a logarithmic scale ranging down to 10-4 of the maximum value (top colour bar). The magnetic field b1 (resp. b2) is calculated using Eq. (29) (resp. Eq. (30)) at the shock. The adiabatic index and Lorentz factor Γ are displayed on a linear scale (bottom colour bars). The map spatial scale is in units of the orbital separation d with the pulsar at (1,0) and the star at (0,0).

Figure 3 shows maps of the various flow quantities in the shocked pulsar wind. The jumps in density single out the reflected shock region. The jump in pressure identifies the interface with the bow shock flow as a shock while the matching pressures identifies the interface with the back flow as a contact discontinuity. The bow shock flow is re-energised by adiabatic compression when it crosses the reflected shock. The magnetic field distribution is identical regardless of the assumption adopted for B at the pulsar termination shock (Sect. 3.4). The highest magnetic field intensities are found at the contact discontinuity with the stellar wind, where streamlines from the bow shock head pile up. The magnetic field increases with the density in the reflected shock region (b ∝ nr). The last two panels show the adiabatic index and the Lorentz factor. The adiabatic index is that of a relativistic gas () when the shock is perpendicular and decreases towards its non-relativistic value () as material flows in the bow shock region due to adiabatic expansion. The bow shock flow accelerates back up to a fraction of the initial Lorentz factor of the free pulsar wind before being slowed down again by the reflected shock. The properties in the back flow region vary little on the scales examined here: the flow speed remains close to v = c/ 3 with , and a slowly varying pressure and density.

4.2. Particle cooling

Fig. 4 Left: evolution of particle energy as a function of the distance l along the streamline (in units of the orbital separation). Right: evolution of the particle energy distribution along the streamline (lighter colours for later times i.e. increasing distance l). The top panels correspond to the streamline starting at o = 48° and the bottom panels to the streamline starting at o = 115° (both are identified by an initial dot in Fig. 2). The distribution in the bottom right panel “jumps” when the particles cross the reflected shock and are re-energised. The evolution is calculated at periastron for our reference simulation.

Our reference model for the emission has ξ = 1 (acceleration timescale), s = 2 (power-law slope of electron spectrum, no Maxwellian), and B calculated using Eq. (29). The parameters ζ_p and ζ_b are set to 1. At the apex of the termination shock at periastron (d = 0.098 AU), ζ_p = 1 corresponds to γ_min ≈ 8 × 10⁴ and ζ_b = 1 corresponds to B ≈ 40 G. Figure 4 illustrates the typical evolution of particles along two streamlines, at periastron φ = 0. We label streamlines by the angle o that their starting point makes with the binary axis, with o = 0° corresponding to the bow shock head on the binary axis, o = 90° perpendicular to the binary axis, 130° ≤ o ≤ 180° to the back shock (see Fig. 2, where the black dots identify the two streamlines for which the particle evolution is shown).

For the first streamline (o = 48°), the initial particle distribution ranges from γ_min = 6.8 × 10⁴ to γ_max = 8.8 × 10⁶ (top panels of Fig. 4). The highest energy electrons radiate away half of their energy on a scale l ≲ 0.003 (in units of the orbital separation, here d = 0.092 AU). For comparison, this is comparable to the spatial resolution at maximum grid refinement in our simulation. Properly resolving the cooling spatial scale within the RHD simulation would require an additional 1−2 levels of refinement, at large computational cost as mentioned in Sect. 3.1. Even the lowest energy electrons cool on a small scale l ≲ 0.05 compared to the length of the streamline. Synchrotron and inverse Compton burnoff at high energy is seen in the evolution of the particle distribution (right panel). The late evolution of the distribution is set by adiabatic losses, which do not modify the shape of the distribution.

For the second streamline (o = 115°), the initial distribution ranges from γ_min = 1.2 × 10⁴ to γ_max = 3.1 × 10⁷. The higher γ_max is due to the lower magnetic field at the shock (Eq. (19)) while the lower γ_min is due to the lower pressure (Eq. (20)). Radiative cooling is weakened by the distance and by the higher Lorentz factor (decreasing the comoving density of stellar photons). There is a moderate evolution of the particle distribution before the particles are re-energised by passing through the reflection shock. Compression at the shock heats the particles and enhances the magnetic field. Radiative cooling is much more important in the subsequent evolution of the highest energy particles.

Fig. 5 Fraction of the mean particle energy remaining after radiation losses along each streamline, assuming the conditions at periastron (φ = 0). The streamlines are ordered by angle o (Fig. 2). Dashed line shows the same at apastron (φ = 0.5).

Without radiative losses, the particles lose ≈ 43% of their energy adiabatically within the box. The particle energy losses are enhanced by radiation. Figure 5 shows the fraction of the total energy losses that are due to radiative losses depending on streamline. For each streamline, the ratio of the specific energy in non-thermal particles ϵ_nt to the thermal energy ϵ is compared at the beginning and at the end of the streamline. $\frac{{\mathit{ϵ}}_{\mathrm{nt}}}{\mathit{ϵ}}\mathrm{\propto }\frac{\mathrm{(}\mathit{\gamma ̂}\mathrm{-}\mathrm{1}\mathrm{)}}{\mathit{p}}{\mathrm{\int }}_{{\mathit{\gamma }}_{\min }}^{{\mathit{\gamma }}_{\max }}\mathit{\gamma }{\mathit{m}}_{\mathrm{e}}{\mathit{c}}^{\mathrm{2}}\frac{\mathrm{d}\mathit{n}}{\mathrm{d}\mathit{\gamma }}\mathrm{d}\mathit{\gamma }\mathit{.}$(35)The value of ϵ_nt/ϵ will be the same at the beginning and at the end of the streamline (ϵ_nt/ϵ)_in = (ϵ_nt/ϵ)_out = ζ_p if the particle losses (or gains) are only due to the adiabatic term i.e. the ratio (ϵ_nt/ϵ)_in/ (ϵ_nt/ϵ)_out plotted in Fig. 5 will be 1 if there are no radiative losses. The figure shows this is not the case in our reference simulation: radiative losses dominate the overall energy losses as the particles propagate along the streamlines. The radiative losses are most important for the streamlines that start close to the apex (o ≲ 45°) and in the reflected shock region. Integrating ϵ_ntṄ over the whole flow, we find that ≈ 70% of the power given to non-thermal particles is lost to radiation within our simulation box. Adiabatic losses have a minor influence on the spectrum and lightcurve: nearly identical results are obtained when our baseline calculation is run without taking adiabatic losses into account.

We can only speculate on the feedback that radiative cooling could have on the flow dynamics, since our simulation does not take it into account. The shock region width is likely to decrease as the plasma loses pressure support, raising the density. Since the magnetic field intensity is tied to the density, this may cause particles to radiate even faster and at a higher synchrotron frequency than in our computation. Thin shell instabilities may also disrupt the interaction region. We leave this for future investigations.

Fig. 6 Emission maps of the shocked pulsar wind at various frequencies (top to bottom) and for orbital phases φ corresponding to the conjunctions (left and right columns). The emission is displayed on a logarithmic scale ranging from 1 down to 10-4, with 1 corresponding to the maximum value of the emission along the orbit at this frequency (i.e. each map is normalised to its maximum value over any pixel and any φ). The VHE maps do not take the pair production opacity into account. The case shown here corresponds to i = 25° in Fig. 8. The spatial scale is in units of the orbital separation d. See Figs. A.1−A.4 for associated movies of the evolution with orbital phase.

4.3. Emission maps

Emission maps were built for the baseline case with i = 25° (Fig. 6). The maps represent the unabsorbed emissivity integrated over azimuth θ (Eq. (10)) $\begin{array}{ccc}\mathrm{\int }{\mathrm{𝒟}}_{\mathrm{obs}}^{\mathrm{2}}\mathit{j}\mathrm{(}\mathit{\nu }\mathit{/}{{\mathrm{𝒟}}_{\mathrm{obs}}}^{\mathrm{)}}\mathit{r}\mathrm{d}\mathit{\theta }\mathit{.}& & \end{array}$This quantity was sampled along each streamline and binned to form maps at 1 keV, 1 MeV, 100 GeV, and 1 TeV. Figure 6 compares the maps at superior (φ = 0.08) and inferior conjunctions (φ = 0.77). Animations showing the evolution at all orbital phases are available online (See Figs. A.1−A.4). Fast cooling concentrates emission at the highest frequencies to thin layers close to the pulsar termination shock (e.g., compare the synchrotron 1 MeV and 1 keV maps). The emission is more concentrated at φ = 0.08 than at φ = 0.77 because the orbital separation is smaller (d ≈ 0.10 AU compared to d ≈ 0.15 AU), leading to stronger radiative losses. The simulation box covers well the emission zones at 1 TeV and 1 MeV, but misses some of the 100 GeV and 1 keV, especially in the back region. Note, however, that the map flux scale is logarithmic so the impact on the overall lightcurve is negligible.

Reheating in the reflection shock region is easily seen in the maps, especially at 1 keV where a significant fraction of the flux may come from this region (and hence escape X-ray absorption by the stellar wind, see Szostek & Dubus 2011). The bow shock emission is concentrated towards the head while the back shock emission covers a much wider area. The VHE emission from the back region suffers less from γγ absorption and actually contributes nearly all the TeV flux when the opacity is highest (around φ = 0.08, see Fig. 8).

Figure 7 presents a different way of looking at where particles cool. We have plotted the integrated contribution of each streamline to the total emission at different frequencies. Streamlines that start at o ≤ 90° correspond to the head of the bow shock, streamlines with o ≥ 130° correspond to the back shock (Fig. 2). The TeV inverse Compton emission originates mostly at streamline angles o larger than for the 100 GeV emission, and the same applies when comparing the MeV and keV synchrotron emission. This is because the streamline initial magnetic field decreases with o, allowing for higher initial particle energies (γ_max). The back shock contribution clearly dominates the absorbed VHE flux at φ = 0.08. At φ = 0.77 the emission is much more concentrated in the streamlines with o ≈ 100°, which pass through the reflected shock and strongly benefit from the relativistic Doppler boost since the flow in the back region is then aligned with the observer line of sight (i = 25°).

Fig. 7 Contribution to the total flux by each streamline for our baseline model with i = 25°. The streamlines are ordered by angle o (Fig. 2). Top panel is for orbital phase φ = 0.08 (superior conjunction), bottom panel for φ = 0.77 (inferior conjunction). The flux fractions in each panel correspond to the four frequencies mapped in Fig. 6. Solid lines correspond to the flux unabsorbed by pair production, dashed lines are for the absorbed flux (affecting the fractions only at 100 GeV and 1 TeV).

4.4. Spectra and lightcurves

Spectral energy distributions and lightcurves were computed for several system inclinations using our baseline model. The results are displayed in Fig. 8. The spectra and lightcurves are normalised by a coefficient $\mathrm{𝒦}\mathrm{=}\mathrm{5}\mathrm{\times }{\mathrm{10}}^{-10}\hspace{0.17em}{\left(\frac{\mathrm{2.5}\hspace{0.17em}\mathrm{kpc}}{\mathit{D}}\right)}^{\mathrm{2}}\left(\frac{\mathit{Ė}\mathrm{p}}{\mathrm{7.6}\mathrm{\times }{\mathrm{10}}^{\mathrm{35}}\hspace{0.17em}\mathrm{erg}\hspace{0.17em}{\mathrm{s}}^{-1}}\right)\hspace{0.17em}\mathrm{erg}\hspace{0.17em}\mathrm{c}{\mathrm{m}}^{-2}\hspace{0.17em}{\mathrm{s}}^{-1}\mathit{,}$(36)where we have made explicit the dependence on the injected power in particles Ė_p = 7.6 × 10³⁵ erg s^-1. Ė_p is related to the pulsar spindown power by Ė = Ė_p(1 + σ). The results can be scaled with Ė_p, or σ, as long as Ṁ_w and Ė change in parallel to keep η = 0.1 (Sect. 2). Figure 8 shows the spectral energy distribution sampled at various orbital phases to highlight the spectral variability, the average spectral energy distribution (thick black line) and (in blue) the average spectrum corresponding to phases 0.45 <φ ≤ 0.9 (INFC) and φ ≤ 0.45 or φ> 0.9 (SUPC), allowing a comparison with the HESS spectral analysis in Aharonian et al. (2006).

The size of the simulation domain limits how far we can follow particle cooling, as the maps of Fig. 6 illustrate. The emission is thus necessarily incomplete below some energy, which we estimate to be ≲ 100 eV (synchrotron emission component) and ≲ 1 GeV (inverse Compton component) – based on comparing spectra obtained with a reduced domain size.

The spectra produce broad band X-ray to TeV emission but, as could be expected (see Dubus 2013), cannot reproduce the peaked GeV emission observed with the Fermi-LAT. This emission component requires a completely different population of electrons, with a narrow distribution in energy. Zabalza et al. (2013) speculated this could arise from the back shock but we find no obvious difference between bow and back shocks. The average spectrum from the back region is shown as a dashed line in the top panels. Emission from the back region can dominate near superior conjunction, when γγ absorption is important (see dashed lightcurve in the bottom panels), but its contribution to the average spectrum remains minor. The spectra of the bow and back region are similar; they would need to have very different acceleration parameters to produce significantly different spectra (Zabalza et al. 2013; Takata et al. 2014). We come back to the question of the origin of the HE gamma-ray emission in Sect. 5.

We show lightcurves for the X-ray (1−10 keV) and VHE (> 100 GeV) gamma-ray bands, where the spectra and orbital modulations are well-known from Suzaku and HESS observations (Takahashi et al. 2009; Aharonian et al. 2006). The lightcurves were computed by integrating F(ν) over the relevant energy range. When the system is seen face-on (i = 0°), the VHE modulation is directly related to the varying stellar photon density which increases both inverse Compton emission and pair production. The synchrotron emission varies little. The synchrotron loss timescale increases at larger orbital separations since τ_sync ~ b^{− 3/2}ν^1/2 ~ d^3/2 at a given frequency. However, the actual size of the computational domain increases as d, while the particle distribution slope s = 2 ensures equal power per particle energy, so the emission should vary roughly as d^0.5, a factor 1.4 from periastron to apastron, a bit more than what the full calculation gives (bottom left panel of Fig. 8).

The emission received by the observer changes dramatically with the inclination angle of the system. The synchrotron emission is changed by relativistic boosting as changes with orbital phase. The bow shock region creates a hollow cone of high-velocity material, surrounding a filled cone of lower-velocity material flowing away from the back shock. Increasing the inclination boosts the X-ray emission around inferior conjunction φ_inf ≈ 0.77, when the shocked flow is oriented towards the observer, and de-boosts the X-ray emission around superior conjunction φ_sup ≈ 0.08, when the flow is directed away. At higher inclination the observer line of sight starts crossing the emission cone of the bow shock along its full length. This result in maximum boost at φ_inf and 40° ≲ i ≲ 50° when the shocked flow is going directly in the direction of the observer. For i ≳ 50° the line of sight crosses first one edge of the hollow cone then the other, resulting in double-peaked emission with a minimum at φ_inf. At i = 90°, the line of sight crosses the cone first at orbital phases 0.54 ≲ φ ≲ 0.60 and then at 0.87 ≲ φ ≲ 0.89, in agreement with the position of the two peaks. Although the cone is symmetric, the second peak is narrower because of the faster orbital motion during the second crossing (nearer to periastron passage).

The VHE emission is influenced by the anisotropy of the inverse Compton and pair production cross-sections. Inverse Compton emission is enhanced around φ_sup, when stellar photons are backscattered towards the observer, and diminished around φ_inf, when the stellar photons are forward-scattered (Fig. 9). For the same reasons, the pair production opacity is important around φ_sup and lower at φ_inf. The latter can be verified by comparing the dark (absorbed) and grey (unabsorbed) lines in the middle panels of Fig. 8. However, the effects of Doppler boosting dominate the VHE lightcurve. Figure 9 shows the expected lightcurve at i = 25° but with the relativistic effects turned off when computing the emission (the electron populations are identical). The X-ray synchrotron lightcurve is nearly constant in the absence of Doppler boosting effects. Comparing with the full calculation (Fig. 8), relativistic effects displace the peak VHE and X-ray emission towards φ_inf and then lead to a double-peaked structure at higher i. Relativistic boosting also hardens the VHE spectrum around φ_inf. The intrinsic anisotropic inverse Compton emission is harder at these orbital phases because scattering is increasingly within the Thomson regime when the stellar photons are closer to being forward-scattered (Dubus et al. 2008). This effect is amplified by the bulk Doppler boosting. The strong spectral evolution with orbital phase at high i can be followed in the top panels of Fig. 8. The grey lines show the spectral energy distribution at different phases. The inverse Compton spectrum pivots around 100 GeV for i = 50°. Above this energy, maximum emission occurs around φ_inf (the INFC spectrum is brighter) whereas, below this pivot energy, the gamma-ray emission peaks around φ_sup (the SUPC spectrum is brighter). The inverse Compton lightcurve at different frequencies can thus behave in anti-phase because of the subtle hardening effects brought about by scattering angle and bulk Doppler boost.

Fig. 8 Spectral energy distributions and lightcurves depending on inclination for our reference model. Top panels: thick black line is the average spectrum, thick coloured lines are the average INFC (dark blue) and SUPC (light blue) spectra (Sect. 4.4), thick dashed line shows the contribution of the back shock region to the average spectrum, thin grey lines show the spectral evolution with orbital phase. Bottom panels: VHE gamma-ray and X-ray lightcurves (thick solid lines). Again, the dashed line shows the back shock contribution. The thin grey line in the middle panels is the unabsorbed VHE flux.

Fig. 9 Same as Fig. 8 for i = 25° and 50° except that relativistic aberrations has not been taken into account when computing the emission.

4.5. Exploring parameter space

We carried out a limited exploration of the parameter space around our reference model. The dependence on inclination i has already been shown in Fig. 8. The dependence on the other parameters, namely ζ_p, ζ_b, ξ, and s, is shown in Fig. 10. The normalisation of the spectra is the same as for our reference model (Eq. (36)). We remind that ζ_p controls the mean energy of the particles (Eq. (20)), ζ_b controls the magnetic field intensity at the shock (Eq. (29)), ξ controls the maximum particle energy at the shock (Eq. (17)), and s is the slope of the injected power-law distribution of electrons.

We successively changed the value of each parameter, keeping the others to their reference value. A higher ζ_b increases synchrotron losses, leading to a pronounced νF_ν ~ ν^{(2 − s)/2} ~ ν⁰ spectrum of cooled particles and lowering the inverse Compton emission component. Conversely, a lower ζ_b increases the inverse Compton component relative to the synchrotron component and, in our case, leads to a hard synchrotron spectrum because the inverse Compton energy losses are in the Klein-Nishina regime. A higher ζ_p increases γ_min and thus narrows down the energy range of the injected power law. The low-energy slope of the synchrotron component for ζ_p = 10 corresponds to the νF_ν ~ ν^4/3 expected for a tail of emission from electrons at a high γ_min whereas, in the ζ_p = 0.1 case, the slope is the νF_ν ~ ν^{(3 − s)/2} ~ ν^1/2 slope expected from uncooled electrons emitting in our frequency range with (smaller) Lorentz factors ≥ γ_min. Changing the acceleration timescale ξ directly impacts the maximum synchrotron frequency, but has no influence on the inverse Compton emission because, in our case, γ_max is always high enough for the interaction with stellar photons to occur in the inefficient Klein-Nishina regime. Finally, changing the slope s of the injected distribution, unsurprisingly, produces a harder synchrotron spectrum when the electron distribution is harder (smaller s). The ratio of synchrotron to inverse Compton emission is also higher because a harder distribution implies more very high energy electrons that radiate more efficiently synchrotron emission compared to inverse Compton emission in the Klein-Nishina regime.

We do not show how the TeV and X-ray modulations were affected by the changes in parameters in Fig. 10. The reason is that the modulations did not change much compared to the reference model. These lightcurves are predominantly shaped by the inclination rather than by changes in the other parameters. All cases are also comparably radiatively-efficient: the bolometric luminosities of the different models vary, in normalised units (Eq. (36)), between ≈13 (ξ = 100) and ≈ 35 (ξ = 0.01). As in the reference case, most of the pulsar power is converted into radiation.

5. Application to LS 5039

These results guided us towards a model reproducing the emission from LS 5039 based on our RHD simulation. This model is compared with the observed spectral energy distribution of LS 5039 in Fig. 11 and the X-ray, MeV, GeV, and TeV lightcurves in Fig. 12. The cost of the calculations (several hours per model) does not allow an extensive exploration of parameter space. The parameter combination given here is only indicative of what seems to work i.e. this is not a best-fit model.

5.1. Model parameters

The main drivers in deriving this model were (1) reproducing the VHE spectral variations; (2) accounting for the comparable X-ray and VHE gamma-ray fluxes; and (3) understanding the origin of the HE gamma-ray emission. We start with the latter.

As the results from the previous section should make clear, the HE gamma-ray emission observed with the Fermi-LAT requires an additional component. In principle, a low value of ξ pushing the synchrotron component to GeV energies might account for the Fermi-LAT spectrum with an exponential cutoff (at the price of supposing a faster-than-Bohm acceleration timescale). However, the HE modulation would then be in phase with the X-ray modulation, which is ruled out by the observations. The HE modulation is actually consistent with expectations for inverse Compton scattering of stellar photons (Abdo et al. 2009).

Fig. 10 Dependence of the spectrum on the model parameters. The spectra should be compared to the baseline model with ζb = ζp = ξ = 1, s = 2 and i = 25° (second column of Fig. 8). For line labels, see caption to Fig. 8.

We explored the possibility that the HE emission could be due to the inverse Compton emission from a narrow Maxwellian distribution of electrons, as we had for the case of PSR B1259-63 in Dubus & Cerutti (2013). We assumed that a fraction of the pulsar wind particles injected at the shock are accelerated to a power law, accounting for the broad band X-ray to TeV emission, while the rest are randomised to this Maxwellian distribution. Adjusting the HE spectrum with the inverse Compton emission from the Maxwellian fixes its mean Lorentz factor to γ_t ≈ 5000. This also fixes ζ_p to ≈ 0.017 through Eq. (20). The available specific internal energy is a priori identical for the Maxwellian and power-law distributions so ζ_p is thus fixed for both populations of electrons. The relative contribution of each population is adjusted by the fraction of the particle density going to each (or, equivalently, total energy since ϵ_nt is the same). The contributions to the average spectrum from each population of particles are highlighted in the top panel of Fig. 11.

The HESS INFC spectrum is hard, best described as a power law of photon index of 1.8 combined with an exponential cutoff at 8.7 TeV. The SUPC spectrum is a steeper power law with an index of 2.5 (Fig. 11). The two spectra pivot at an energy ≈ 200 GeV. Reproducing both the hard INFC spectrum and the comparable levels of X-ray and VHE emission turned out to be difficult. The models explored in Sect. 4 all cut off around 100 GeV because the high-energy particles responsible for this emission are strongly cooled by synchrotron losses. Lowering ζ_b decreases the mean magnetic field, hence increases the VHE cutoff, but also lowers the X-ray flux relative to the VHE gamma-ray flux. There is a trade-off between having enough synchrotron emission to account for the X-ray flux, and lowering ζ_b to enable the highest-energy particles to radiate enough VHE photons. As discussed in Sect. 4.5, changing ξ has little influence on the VHE spectrum so we kept ξ = 1. Limited exploration showed the best agreement was obtained by slightly lowering ζ_b to 0.5 and by taking an injection slope s = 1.5 instead of 2.

The inclination has an effect on both lightcurves and spectral variations. The VHE spectral variations are more pronounced with higher i although too high an inclination results in a INFC spectrum with a flux that is too low (Fig. 8). A high inclination also results in a pronounced dip of emission at φ_inf = 0.77. The HESS lightcurve appears double-peaked with a shallow minimum at φ_inf, favouring a model with 25° <i< 50° (Fig. 8). The X-ray modulation is single-peaked at φ_inf, favouring models towards the low end of this range of i. Using i = 35° turned out to be a good compromise.

The “best-adjusted” model shown in Figs. 11, 12 has ζ_p ≈ 0.017, ζ_b = 0.5, ξ = 1, s = 1.5, i = 35° and an injected population of particles consisting of a Maxwellian plus a power law. An additional parameter is the value of η = 0.1 that we fixed all throughout this study. The hard injection spectral slope s = 1.5 hints at reconnection rather than Fermi acceleration. Adjusting the model to the observations requires that the injected power in particles is Ė_p ≈ 10³⁵ erg s^-1. The majority of the particles or available power (88%) goes to the Maxwellian at the termination shock. Only 12% goes to the particles distributed as a power law. However, these power-law particles are more radiatively efficient than those in the Maxwellian: about 80% of the power injected as a power law ends up radiated away within the simulation domain (like the models shown in Sect. 4) compared with only 25% of the power injected as a Maxwellian. Hence, most of the radiation comes from the injected power law but most of the power is in the Maxwellian. As a consequence, most of the particles evolve adiabatically in this model.

Fig. 11 Our best-adjusted model to the observed spectral energy distribution of LS 5039. Top: the black curve is the average model spectrum; the synchrotron (~ 1 eV) and inverse Compton (~ 1 GeV) contributions from the Maxwellian population of electrons are in dark blue; the contributions from the power-law population are in light blue. Bottom: the thick, dark blue curve (resp. light blue curve) represents the INFC (resp. SUPC) spectrum, the thin grey curves represent the evolution of the SED in orbital phase steps of 1/30. From left to right, the data shown are: the orbital maximum and minimum X-ray bowties from Suzaku (Takahashi et al. 2009), the BATSE data points and INTEGRAL bowtie (Harmon et al. 2004; Hoffmann et al. 2009), the average COMPTEL data points with bowtie (Collmar & Zhang 2014), the average Fermi-LAT spectrum with data points from 100 MeV to 50 GeV and the best-fit power law with exponential cutoff (Hadasch et al. 2012), the HESS INFC (dark blue bowtie) and SUPC spectra (light blue bowtie) with the associated data points from 100 GeV to 30 TeV (Aharonian et al. 2006).

5.2. Comparison to the observations

The (rough) adjustment provides a reasonable, albeit imperfect, description of the data. The X-ray flux is too low by a factor ~ 2 in the present model. X-ray emission arising beyond our computational domain might account for this mismatch (Sect. 4.4 and Fig. 6). A higher ζ_b would raise the X-ray flux, but lower the VHE cutoff in the INFC spectrum to values that would not be consistent with the observations. The largest discrepancy is the COMPTEL data at MeV energies, recently associated with LS 5039 (Collmar & Zhang 2014). Although the spectral slope of the model, as well as its evolution with orbital phase (Fig. 12), are compatible with the observations (the COMPTEL spectrum is harder at SUPC phases), the flux is clearly underestimated by an order of magnitude. A higher ζ_b would raise the synchrotron luminosity but a cooling break in the synchrotron spectrum is hard to avoid and this would also be at the expense of the hard VHE INFC spectrum. If ζ_b is low then the synchrotron spectrum can extend from X rays to MeV, as in Takahashi et al. (2009) where adiabatic cooling is assumed to dominate over radiative cooling, but our model shows the inverse Compton component would then be too narrow and luminous (see Fig. 10). Our model is unlikely to be able to account for the observed MeV flux without additional ingredients (discussed in Sect. 6), unless the flux is contaminated by diffuse emission or the flux from other sources due to the poor angular resolution at these energies. Progress is much needed in this difficult observational band.

The orbital modulations from the model are compared with the observations in Fig. 12. The fluxes were calculated in each band in the units given by the observations. Because of the mismatch in X-ray and MeV fluxes, we had to multiply the model flux by a factor 2 and 10 (respectively) to obtain a level comparable to the observations. The lightcurves are reasonably well reproduced. A slightly larger inclination would deepen the VHE minimum at φ ≈ 0.7. The Maxwellian component dominates in the HE band, with a modulation in anti-phase with the other wavelengths. The simulation output shows that the HE emission is concentrated at the head of the bow shock and is not affected much by relativistic boosting. The HE modulation is dominated by the variation with phase of the stellar photon density and scattering angle, resulting in peak emission at ≈ φ_sup. The Maxwellian component also contributes some flux in the 10−30 MeV range. A small peak at φ ≈ 0.75 is visible in the TeV, MeV, and X-ray lightcurves. This appears to be due to the observer line of sight grazing the top of the emission cone at this orbital phase when i = 35°. The small-scale, stable, structures observed in the X-ray modulation lightcurve (Kishishita et al. 2009) might thus be directly related to small structures in the shocked flow that are probed when our line of sight passes through. Emission from the back shock region provides some residual TeV flux around φ_sup, when emission from the bow shock region is strongly absorbed by pair production. This is still insufficient to explain the VHE detections. Emission from the cascade, initiated when the newly created e⁺e⁻ pairs are able to radiate VHE gamma rays, is very likely to be responsible for the residual flux at φ_sup. Cerutti et al. (2010) found from their study of cascade emission that a good fit required an inclination i ≈ 40°, consistent with the present model.

The synchrotron emission from the Maxwellian component peaks in the visible band, where it will be difficult to detect against the bright V = 11.2 companion O star. This emission is modulated, with a lightcurve shape (not shown) similar to the MeV modulation. The peak-to-peak amplitude of the V band modulation is 0.25 mJy. This translates into a 1.3 mmag modulation, below the current upper limit of 2 mmag (Sarty et al. 2011).

Fig. 12 Comparison between the LS 5039 model (Fig. 11) and the observed lightcurves. The model 10−30 MeV and 1−10 keV lightcurves were multiplied by a factor 10 and 2, respectively, to match the observations in Takahashi et al. (2009), Collmar & Zhang (2014). The grey lightcurve in the 0.1−10 GeV and 10−30 MeV panels represents the contribution from the relativistic Maxwellian component. The VHE and HE gamma-ray lightcurves are taken from Aharonian et al. (2006) and Abdo et al. (2009), respectively.

6. Discussion

6.1. Influence of the hydrodynamics on the flow emission

Our motivation for developing this radiative code, based on a relativistic hydrodynamical simulation, was to obtain a more realistic and coherent treatment of the emission geometry, adiabatic losses, and Doppler boosting. We discuss these points in turn below.

The simulation shows a complex shock structure, fully containing the pulsar wind, with a reflected shock that re-energises the shocked pulsar wind. In our models, most of the highest-energy emission remains concentrated towards the head of the bow shock where the electrons cool quickly. Hence, both the spectral energy distribution and the modulations are predominantly shaped by the head of the bow shock region. This is fortunate as the back shock and reflected shock structures are likely to have some dependence on our choice of Mach number and could change in the presence of orbital motion, strong mixing, dynamically important radiative losses or magnetic fields, etc. It is primarily at lower frequencies, notably in X rays, that these structures contribute significantly, when the particles cool more slowly (on larger spatial scales) and/or re-heated to mild energies. These conclusions depend on the assumptions that we made on what happens at the various shocks, namely that there is no difference in particle acceleration between the bow and back shock, and that the reflected shock only compresses particles. The first assumption is likely to be incorrect at some level because pulsar winds are not isotropic. Some latitude dependence is expected in the pulsar wind due, for instance, to differences between the propagation of the high-latitude regions and the equatorial region (defined by the pulsar rotation axis) where the pulsar wind is striped and prone to reconnection. This may be manifest in a latitude dependence of the Lorentz factor of the pulsar wind, as seems to be required by models of the Crab pulsar wind nebula (Bogovalov & Khangoulian 2002; Porth et al. 2014), and/or a dependence of the pulsar wind magnetisation σ with distance (Zabalza et al. 2013; Takata et al. 2014). A latitude-dependent ϵ or B could have important observational signatures, even if there is no dramatic change in the structure of the flow (Vigelius et al. 2007; Bogovalov et al. 2012). On the second assumption, having the reflected shock only adiabatically compress the particles has some justification because it is usually found to be difficult to accelerate particles at shock in pair plasmas except in special circumstances (very low magnetisation, shock-induced reconnection of the striped wind). Yet, we cannot exclude that particles are re-accelerated to a power law, with some influence on the overall emission. We leave the exploration of these possibilities to future work.

The models we have explored are radiatively very efficient despite the fast flow timescale, the size of the emission region, and the decreasing magnetic field strength with distance. Adiabatic losses play a minor role in the high-energy emission, excluding them as the main driver of the X-ray and VHE modulation in LS 5039 as proposed by Takahashi et al. (2009). In principle, strong radiative losses should be taken into account in the dynamics of the flow region. These would result in a narrower and denser shocked flow, experiencing a higher magnetic field (Sect. 4.2). However, in order to reproduce the spectral component seen with the Fermi-LAT, we have proposed that most of the particles in the shocked flow are actually injected in the form of a Maxwellian component rather than accelerated to a power-law distribution. If the particles in the Maxwellian were only randomised at the shock, their mean Lorentz factor γ_t corresponds to the Lorentz factor Γ_p of the pulsar wind so Γ_p ≈ γ_t = 5000. These particles do not cool efficiently in the shocked flow and actually dominate the energy budget⁵. Hence, perhaps counter-intuitively, the flow remains essentially adiabatic and the assumption of the simulation is verified.

Doppler boosting has a very strong effect in shaping the modulation lightcurves. The geometry is basically a rotating cone whose emission is boosted at the phases where its wings cross the observer line of sight. The result is a double-peaked modulation at high inclinations, affecting the synchrotron emission and the VHE inverse Compton emission. In our models, the modulation due to the anisotropy of the inverse Compton cross-section is more important than the Doppler modulation only at lower energies, when the scattering occurs in the Thomson regime. The dependence of the lightcurve shape on inclination allows us to constrain i to ≈ 35°. A side effect of Doppler boosting is that it contributes to steepening the VHE emission near superior conjunction, erasing the dip around a few 100 GeV in the SUPC spectrum due to γγ absorption and typically seen in previous models (e.g. Dubus et al. 2008; Khangulyan et al. 2008; Yamaguchi & Takahara 2012). Our simulation assumed a pulsar wind Lorentz factor Γ_p = 7, not quite high enough to obtain ultra-relativistic shock conditions along most of the bow shock. A higher Γ_p could enhance the contribution from the wings, although we consider this unlikely given the sharp decrease in flux as the shock becomes mostly transverse (Fig. 6). Even if a simulation with higher Γ_p is desirable, we do not expect our results to change much.

6.2. Towards more realistic models

Adjusting to the observations of LS 5039 highlights both the difficulties and the progress to be expected from the present approach. Parameters that would be well-suited to the observations in one energy band are easily discarded as a result of the comprehensive approach taken here, because they fail to reproduce the observations in another band.

We have discussed in Sect. 5 the difficulties in reproducing both the level of the X-ray flux and the hard VHE spectrum. Part of the difficulty may be alleviated if cascade emission is taken into account. It would contribute to the VHE flux at all phases, not only at superior conjunction, and also to the X-ray emission via synchrotron radiation from the pairs (Bednarek 2007; Bosch-Ramon et al. 2008; Cerutti et al. 2010). Combining the present model with a 3D cascade model represents a daunting task.

The COMPTEL flux level presents a similar challenge to models. The modulation, in phase with the X-rays and in anti-phase with the HE gamma rays, excludes that the COMPTEL emission arises from the same electrons responsible for the HE gamma-ray emission. It is more natural to attribute it to the extension of the synchrotron emission (Takahashi et al. 2009), yet it appears very difficult to achieve this without a high magnetic field. As explained above in Sect. 6.1, more complex dependencies of ϵ or B, motivated by the physics of pulsar winds and their termination shock, might be able to simultaneously explain the VHE emission and the strong synchrotron emission.

Our results are based on a single simulation with a given η. They should hold qualitatively for other values of η. The strongest impact would certainly be on the value of the inclination required to adjust the observations since the cone opening angle depends directly on η. More subtle effects may appear if η changes along the orbit. This is to be expected at some level in LS 5039 because the stellar wind is still in its acceleration phase when it encounters the pulsar wind at a distance of one stellar radius from the surface of the star. Using a beta law for the stellar wind velocity, we find that this leads to relative changes of 20% in η. The opening angle of the cone would be slightly higher at periastron than at apastron. Much stronger effects are expected in the case of gamma-ray binaries like LS I+61°303 and PSR B1259-63 where the pulsar wind interacts with a dense equatorial outflow from the companion star. Modelling these systems requires orbit-dependent simulations. Besides the orbital phase dependency of η, such simulations will also be able to address the impact of orbital motion on the shape of the interaction cone. We expect that the leading arm (the part of the shocked pulsar wind that moves into the stellar wind due to orbital motion) will be compressed and the trailing arm will expand (see Lamberts et al. 2012; Bosch-Ramon et al. 2015 and references therein). The impact on the lightcurves should be limited since the emission arises mostly from the innermost, less-affected regions. The back shock in our simulation looks similar to the “Coriolis shock” identified by Bosch-Ramon & Barkov (2011) and Bosch-Ramon et al. (2012, 2015) in simulations including orbital motion. We speculate that the presence of a back shock is not related to the Coriolis force, and note that some of our previous simulations and those of Bogovalov et al. (2008) indeed show full confinement without orbital motion, albeit with a different back shock geometry (Sect. 4.1). Dedicated 2.5D (cylindrical) relativistic simulations would be useful to clearly define the conditions for full confinement, the shape of the structure, and resolve the issue.

The parameter ζ_b imposes the value of the magnetic field at the apex of the bow shock: B ≈ 20 G at a distance ≈ 3 × 10¹¹cm from the pulsar and at periastron passage. Taking into account the toroidal (beyond the light cylinder radius r_LC = cP/ 2π) and dipolar (within light cylinder) nature of the magnetic field, the intensity at the pulsar surface is B₀ ≈ 20 (3 × 10¹¹ cm /r_LC)(r_LC/r_ns)³ ≈ 1.4 × 10¹² (P/ 0.1 s)²G where P is the pulsar spin period and for a neutron star radius r_ns ≈ 10⁶cm. This value of B₀, for a given P = 0.1 s, is standard for rotation-powered gamma-ray pulsars (Abdo et al. 2013).

The total injected power in particles Ė_p ≈ 10³⁵ erg s^-1 of our best model is also standard for pulsars detected in gamma rays (Abdo et al. 2013). However, combining Eq. (29) for the magnetic field with Ė = Ė_p(1 + σ), we find that our model requires σ ≈ 1. The pulsar spindown power is equally spread between magnetic and kinetic energy. Such a value is much higher than has been usually assumed in pulsar wind nebulae, starting with the work of Kennel & Coroniti (1984), though it is not necessarily surprising since pulsar winds are thought to be launched with very high values of σ and to convert the magnetic energy to kinetic energy as they propagate (the “σ problem”, see e.g. the reviews by Kirk et al. 2009 and Arons 2011). The shock is much closer to the pulsar in LS 5039 (≈ 3 × 10¹¹cm) than in pulsar wind nebulae (0.1 pc in the Crab nebula) so a higher value of σ is not problematic per se. A more worrying issue is that with σ ≈ 1 the assumption of hydrodynamics breaks down. The higher magnetic field at the termination shock means that less of the pulsar wind energy will be transferred to the particles. A full RMHD simulation should be carried out to investigate whether a substantial magnetisation alleviates some of the difficulties we have encountered in reproducing the observations.

Finally, a pulsar spindown power of Ė = 2 × 10³⁵ erg s^-1 implies a stellar wind mass loss rate Ṁ_w = 5 × 10^-9M_⊙ yr^-1 given η = 0.1 and v_w = 2000 km s^-1. This is at the low end of the range of estimated Ṁ_w, even if we take into account that the spindown power would have to be increased by a factor at least 3 to account for the 1−30 MeV luminosity of ≈ 6 × 10³⁵ (D/ 2.5 kpc)² erg s^-1 (assuming we could reproduce the peculiar spectral shape). The mass loss rate estimated from Hα measurements give values in the range 2 − 75 × 10^-8 M_⊙ yr^-1 (McSwain et al. 2004; Casares et al. 2005; Sarty et al. 2011). However, wind clumping is known to bias this estimator, leading to mass loss rates that can be overestimated by a factor ≈ 20 for O6.5 stars like LS 5039 (Fullerton et al. 2006). Indeed, the lack of signatures from X-ray thermal emission or absorption in the stellar wind favours the lower end of the range of estimated Ṁ_w (Szostek & Dubus 2011; Zabalza et al. 2011). Alternatively, η may be smaller than the value we assumed.

7. Conclusions

We have developed a post-processing radiative code to investigate high-energy non-thermal emission based on relativistic hydrodynamical simulations. Our code includes synchrotron emission, anisotropic inverse Compton emission, the opacity due to pair production at VHE, and takes into account relativistic effects using the velocity field from the simulation. The particle energy distribution is evolved according to the adiabatic losses derived from the simulation and radiative losses. Our goal was to provide a coherent model of the spectral modulations observed from X rays to VHE gamma rays in the gamma-ray binary LS 5039.

• (i)

  The simulation shows a complex shock structure even whenorbital motion is neglected. The pulsar flow is fully confined by abow shock and a back shock. The presence of the back shockinduces a reflected shock in the bow shock region. The back shockand reflected shock have a limited impact on the overall emission.

• (ii)

  The VHE emission remains very concentrated towards the apex of the bow shock and strongly absorbed by pair production at superior conjunction. The back shock contribution dominates the VHE flux at superior conjunction but its flux is insufficient to explain the HESS detection without emission from the pair cascade. The back shock thus has a very minor influence on the gamma-ray emission from the system. The X-ray emission region is much larger, which will help smooth out X-ray absorption signatures from the stellar wind (Szostek & Dubus 2011).

• (iii)

  Doppler boosting plays the major role in modulating the X-ray and VHE emission with orbital phase. Its impact is predominantly set by the inclination of the system i, with double-peaked lightcurves expected at high i. We constrain the inclination of LS 5039 to i ≈ 35°.

• (iv)

  There is tension between the hard VHE spectrum and the level of X-ray emission as they require differing intensities of the magnetic field. This issue is aggravated by the recent COMPTEL detection that, if fully attributed to LS 5039, implies an even stronger synchrotron component (hence higher B) and a sharp cutoff between 10 and 100 MeV. These observations cannot be accommodated in our current model. Possible options that may ease the issue include missing X-ray emission from the simulation box, a more intense magnetic field in the regions where radiative cooling is strong (leading to a denser flow and a more compressed B), contributions from the pair cascade triggered by the absorption of VHE gamma rays, and a latitude or distance-dependent magnetisation σ or wind Lorentz factor Γ_p.

• (v)

  We attribute the Fermi-LAT emission component to particles randomised to a Maxwellian distribution at the shock, as shown by simulations of particle acceleration at pair-dominated shocks. This implies that the Lorentz factor of the wind is Γ_p ≈ 5000. We find that these particles represent the bulk of the power injected in particles, with only 12% going to the electrons accelerated to a power law. Synchrotron emission from the Maxwellian population produces a weak (~1 mmag) orbital modulation of the optical flux.

• (vi)

  The power-law electrons radiate very efficiently, while the “thermal” particles lose energy primarily through adiabatic losses. A modest overall injected power of a few 10³⁵ erg s^-1 is sufficient to account for the broad band X-ray and TeV emission. For our choice of η = 0.1 this implies a stellar wind mass loss rate of the order of 10^-8M_⊙ yr^-1 at the low end of currently estimated values.

• (vii)

  This power, combined with the magnetic field intensity required by our best model, implies a pulsar magnetisation σ ≈ 1. Such a high value supports the picture that has pulsar winds launched with high σ, but fails our assumption of hydrodynamics. Relativistic MHD simulations will be required to further investigate the issue and, perhaps, resolve some of the difficulties encountered in reproducing the observations.

While gamma-ray binaries may be expected to shed light into the processes involved in propagation and termination of pulsar winds, we believe that robust conclusions will require the type of coherent approach linking dynamical and radiative aspects that we have explored here.

Online material

Movie of Fig. A.1 Access here

Movie of Fig. A.2 Access here

Movie of Fig. A.3 Access here

Movie of Fig. A.4 Access here

Acknowledgments

We thank Geoffroy Lesur for his advice and for allowing computations on his private machine. This work was partly supported by the European Community via contract ERC-StG-200911, by the French “Programme National Hautes Énergies”, and by the Centre National d’Études Spatiales. A.L. is supported by the UWM Research Growth Initiative, the NASA ATP program through NASA grant NNX13AH43G, and NSF grant AST-1255469. The RHD simulations have been performed using HPC resources from GENCI- [CINES] (Grant 2013046391) and Texas Advanced Computing Center (TACC) at the University of Texas at Austin (Grant TG-AST-130004).

References

Appendix A: Movies

Fig. A.1 Movie showing the evolution with orbital phase of the emission at 1 keV for our reference model (see Fig. 6). The movie is available online.

Fig. A.2 Movie showing the evolution with orbital phase of the emission at 1 MeV for our reference model (see Fig. 6). The movie is available online.

Fig. A.3 Movie showing the evolution with orbital phase of the emission at 1 TeV for our reference model (see Fig. 6). The movie is available online.

Fig. A.4 Same as Fig. A.3 but with the 1 TeV flux after γγ absorption (averaged over the azimuth θ). The movie is available online.

All Figures

Fig. 1 Model geometry showing our choices of axis and angles. The interaction has cylindrical symmetry around the z axis. The two thick solid lines in the (x,z) plane represent the pulsar wind termination shock and the contact discontinuity with the shocked stellar wind. Dashed lines represent the same in the (y,z) plane. Particles injected at the termination shock follow streamlines in the region in between these two surfaces.

In the text

	Fig. 2 RHD flow in the shocked pulsar wind. The star is at the origin (0,0) and the pulsar is at (1,0) in units of the orbital separation d. A few selected streamlines have been plotted. Three different regions have been coloured. They correspond to the bow shock (blue), the reflected shock (light blue), and the back shock (dark blue) regions.
In the text

Fig. 3 Maps of various quantities in the shocked pulsar wind. The particle density n, pressure p, and magnetic field are displayed on a logarithmic scale ranging down to 10-4 of the maximum value (top colour bar). The magnetic field b1 (resp. b2) is calculated using Eq. (29) (resp. Eq. (30)) at the shock. The adiabatic index and Lorentz factor Γ are displayed on a linear scale (bottom colour bars). The map spatial scale is in units of the orbital separation d with the pulsar at (1,0) and the star at (0,0).

In the text

Fig. 4 Left: evolution of particle energy as a function of the distance l along the streamline (in units of the orbital separation). Right: evolution of the particle energy distribution along the streamline (lighter colours for later times i.e. increasing distance l). The top panels correspond to the streamline starting at o = 48° and the bottom panels to the streamline starting at o = 115° (both are identified by an initial dot in Fig. 2). The distribution in the bottom right panel “jumps” when the particles cross the reflected shock and are re-energised. The evolution is calculated at periastron for our reference simulation.

In the text

	Fig. 5 Fraction of the mean particle energy remaining after radiation losses along each streamline, assuming the conditions at periastron (φ = 0). The streamlines are ordered by angle o (Fig. 2). Dashed line shows the same at apastron (φ = 0.5).
In the text

Fig. 6 Emission maps of the shocked pulsar wind at various frequencies (top to bottom) and for orbital phases φ corresponding to the conjunctions (left and right columns). The emission is displayed on a logarithmic scale ranging from 1 down to 10-4, with 1 corresponding to the maximum value of the emission along the orbit at this frequency (i.e. each map is normalised to its maximum value over any pixel and any φ). The VHE maps do not take the pair production opacity into account. The case shown here corresponds to i = 25° in Fig. 8. The spatial scale is in units of the orbital separation d. See Figs. A.1−A.4 for associated movies of the evolution with orbital phase.

In the text

Fig. 7 Contribution to the total flux by each streamline for our baseline model with i = 25°. The streamlines are ordered by angle o (Fig. 2). Top panel is for orbital phase φ = 0.08 (superior conjunction), bottom panel for φ = 0.77 (inferior conjunction). The flux fractions in each panel correspond to the four frequencies mapped in Fig. 6. Solid lines correspond to the flux unabsorbed by pair production, dashed lines are for the absorbed flux (affecting the fractions only at 100 GeV and 1 TeV).

In the text

Fig. 8 Spectral energy distributions and lightcurves depending on inclination for our reference model. Top panels: thick black line is the average spectrum, thick coloured lines are the average INFC (dark blue) and SUPC (light blue) spectra (Sect. 4.4), thick dashed line shows the contribution of the back shock region to the average spectrum, thin grey lines show the spectral evolution with orbital phase. Bottom panels: VHE gamma-ray and X-ray lightcurves (thick solid lines). Again, the dashed line shows the back shock contribution. The thin grey line in the middle panels is the unabsorbed VHE flux.

In the text

	Fig. 9 Same as Fig. 8 for i = 25° and 50° except that relativistic aberrations has not been taken into account when computing the emission.
In the text

	Fig. 10 Dependence of the spectrum on the model parameters. The spectra should be compared to the baseline model with ζb = ζp = ξ = 1, s = 2 and i = 25° (second column of Fig. 8). For line labels, see caption to Fig. 8.
In the text

Fig. 11 Our best-adjusted model to the observed spectral energy distribution of LS 5039. Top: the black curve is the average model spectrum; the synchrotron (~ 1 eV) and inverse Compton (~ 1 GeV) contributions from the Maxwellian population of electrons are in dark blue; the contributions from the power-law population are in light blue. Bottom: the thick, dark blue curve (resp. light blue curve) represents the INFC (resp. SUPC) spectrum, the thin grey curves represent the evolution of the SED in orbital phase steps of 1/30. From left to right, the data shown are: the orbital maximum and minimum X-ray bowties from Suzaku (Takahashi et al. 2009), the BATSE data points and INTEGRAL bowtie (Harmon et al. 2004; Hoffmann et al. 2009), the average COMPTEL data points with bowtie (Collmar & Zhang 2014), the average Fermi-LAT spectrum with data points from 100 MeV to 50 GeV and the best-fit power law with exponential cutoff (Hadasch et al. 2012), the HESS INFC (dark blue bowtie) and SUPC spectra (light blue bowtie) with the associated data points from 100 GeV to 30 TeV (Aharonian et al. 2006).

In the text

Fig. 12 Comparison between the LS 5039 model (Fig. 11) and the observed lightcurves. The model 10−30 MeV and 1−10 keV lightcurves were multiplied by a factor 10 and 2, respectively, to match the observations in Takahashi et al. (2009), Collmar & Zhang (2014). The grey lightcurve in the 0.1−10 GeV and 10−30 MeV panels represents the contribution from the relativistic Maxwellian component. The VHE and HE gamma-ray lightcurves are taken from Aharonian et al. (2006) and Abdo et al. (2009), respectively.

In the text

	Fig. A.1 Movie showing the evolution with orbital phase of the emission at 1 keV for our reference model (see Fig. 6). The movie is available online.
In the text

	Fig. A.2 Movie showing the evolution with orbital phase of the emission at 1 MeV for our reference model (see Fig. 6). The movie is available online.
In the text

	Fig. A.3 Movie showing the evolution with orbital phase of the emission at 1 TeV for our reference model (see Fig. 6). The movie is available online.
In the text

	Fig. A.4 Same as Fig. A.3 but with the 1 TeV flux after γγ absorption (averaged over the azimuth θ). The movie is available online.
In the text