© ESO, 2013

1. Introduction

Composed of a massive star and a compact object, gamma-ray binaries emit most of their power at GeV and TeV energies. Up to now a handful of such systems have been detected: PSR B1259-63 (Aharonian et al. 2005b), LS 5039 (Aharonian et al. 2005a), LSI +61°303 (Albert et al. 2006), and more recently 1FGL J1018.6-5856 (Fermi LAT Collaboration et al. 2012) and HESS J0632+057 (Bongiorno et al. 2011).

Pulsed radio emission in PSR B1259-63 (Johnston et al. 1992) indicates the compact object is a fast-rotating pulsar, while the nature of the compact object in the other systems is less certain. In PSR B1259-63 the high-energy emission probably arises from the interaction between the relativistic wind from the pulsar and the wind from the companion star (Tavani et al. 1994). It displays extended radio emission showing a structure looking like a cometary tail whose orientation changes with orbital phases (Moldón et al. 2011a). It is interpreted as the evolution of a shock as the pulsar orbits around the massive star. The similarities in the variable high-energy emission and the extended radio emission between PSR B1259-63 and the other detected gamma-ray binaries suggest the wind collision scenario is at work in all the systems (Dubus 2006).

These binaries share a common structure with colliding wind binaries composed of two massive stars. Colliding stellar winds have been extensively studied in the past few years using numerical simulations (see e.g. Parkin et al. 2011) that have highlighted the importance of various instabilities in the colliding wind region. The Kelvin-Helmholtz instability (KHI) can have a strong impact, especially when the velocity difference between the winds is important. Lamberts et al. (2012, hereafter Paper II) show that it may destroy the expected large scale spiral structure. However, strong density gradients have a stabilizing effect. For strongly cooling winds, the shocked region is very narrow and the non-linear thin shell instability (Vishniac 1994) can strongly distort it (Lamberts et al. 2011, hereafter Paper I). For highly eccentric systems, the instability may grow only at periastron, where the higher density in the collision region leads to enhanced cooling (Pittard 2009). A very high resolution is required to model this instability, and up to now large scale studies have been limited by their numerical cost.

Both these instabilities create a turbulent colliding wind region, and may be responsible for the variability observed in stellar colliding wind binaries. They can distort the positions of the shocks, thereby shifting the location of particle acceleration. Strong instabilities also lead to a significant amount of mixing between the winds (Paper II). This may enhance thermalization of non-thermal particles and decrease their high-energy emission. All these effects may be at work in gamma-ray binaries to affect their structure and emission properties (Bosch-Ramon & Barkov 2011; Bosch-Ramon et al. 2012).

The first numerical studies of the interaction between pulsar winds and their environment were performed in the hydrodynamical limit (van der Swaluw et al. 2003; Bucciantini 2002). While providing good insight into the overall structure, a relativistic temperature or bulk motion of the plasma can have subtle effects. For instance, the normal and transverse velocities are coupled via the Lorentz factor even for a normal shock. The specific enthalpy of a relativistic plasma is also greater than for a cold plasma. A relativistic treatment is also necessary for estimating the impact of Doppler boosting on emission properties of the flow. Relativistic simulations have studied the morphology of the interaction between pulsar winds and supernova remnants (Bucciantini et al. 2003), the interstellar medium (Bucciantini et al. 2005), or the wind of a companion star (Bogovalov et al. 2008). Studies have focussed on the impact of magnetization and anisotropy in the pulsar wind or orbital motion (Vigelius et al. 2007; Bogovalov et al. 2012). They also provide models for non-thermal emission from pulsar wind nebulae (Del Zanna et al. 2006; Volpi et al. 2008; Komissarov & Lyubarsky 2004).

One of our goals is to determine the importance of relativistic effects both on the structure and stability of gamma-ray binaries. We perform a set of 2D simulations of gamma-ray binaries and study the impact of the momentum flux ratio and the Lorentz factor of the pulsar wind on the colliding wind region. We also determine the impact of the Kelvin-Helmholtz instability on the colliding wind region (Sect. 3). We therefore extend the hydrodynamics code RAMSES to special relativistic hydrodynamics (Sect. 2). We explain the implementation within the AMR scheme in particular and present some tests (Appendix A) that validate this new code for relativistic hydrodynamics. We discuss how our results can be extrapolated to more realistic values for the Lorentz factor of the pulsar wind (Sect. 4) and conclude (Sect. 5).

2. Relativistic hydrodynamics with RAMSES

RAMSES (Teyssier 2002) is a numerical method for multidimensional astrophysical hydrodynamics and magnetohydrodynamics (Fromang et al. 2006). RAMSES uses an upwind second-order Godunov method. The vector of conserved variables U is averaged over the volume of the cells and fluxes F at cell interfaces are averaged in time. The updates in time are then given by $\frac{{{U}}_{\mathit{i}}^{\mathit{n}\hspace{0.17em}\mathrm{+}\hspace{0.17em}\mathrm{1}}\mathrm{-}{{U}}_{\mathit{i}}^{\mathit{n}}}{\mathrm{\Delta }\mathit{t}}\mathrm{+}\frac{{{F}}_{\mathit{i}\hspace{0.17em}\mathrm{+}\hspace{0.17em}\mathrm{1}\mathit{/}\mathrm{2}}^{\mathit{n}\hspace{0.17em}\mathrm{+}\hspace{0.17em}\mathrm{1}\mathit{/}\mathrm{2}}\mathrm{-}{{F}}_{\mathit{i}\hspace{0.17em}\mathrm{-}\hspace{0.17em}\mathrm{1}\mathit{/}\mathrm{2}}^{\mathit{n}\hspace{0.17em}\mathrm{+}\hspace{0.17em}\mathrm{1}\mathit{/}\mathrm{2}}}{\mathrm{\Delta }\mathit{x}}\mathrm{=}\mathrm{0}\mathit{,}$(1)where the subscript i stands for the index of the cell and n refers to the time step.

RAMSES uses a Cartesian grid and allows the use of adaptive mesh refinement (AMR), which enables the resolution to be locally increased at a reasonable computational cost.

2.1. Equations of relativistic hydrodynamics (RHD)

Throughout this paper the speed of light is c ≡ 1. In the frame of the laboratory, the 3D-RHD equations for an ideal fluid can be written as a system of conservation equations (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}\\ & & \end{array}$These equations can be expressed in compact form $\frac{\mathit{\partial }{U}}{\mathit{\partial t}}\mathrm{+}\sum _{\mathit{k}}^{\mathit{N}}\frac{\mathit{\partial }{{F}}_{\mathit{k}}}{\mathit{\partial }{\mathit{x}}_{\mathit{k}}}\mathrm{=}\mathrm{0}\mathit{,}$(5)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 h}\mathrm{-}\mathit{P}\end{array}\right)\mathrm{·}$(6)The fluxes F_k along each of the N directions are given by ${{F}}_{\mathit{k}}\mathrm{=}\left(\begin{array}{c}\\ \mathit{\rho }\mathrm{\Gamma }{\mathit{v}}_{\mathit{k}}\\ \mathit{\rho h}{\mathrm{\Gamma }}^{\mathrm{2}}{\mathit{v}}_{\mathit{j}}{\mathit{v}}_{\mathit{k}}\mathrm{+}\mathit{P}{\mathit{\delta }}^{\mathit{jk}}\\ \mathit{\rho h}{\mathrm{\Gamma }}^{\mathrm{2}}{\mathit{v}}_{\mathit{k}}\end{array}\right)\mathit{,}$(7)where D is the density, M the momentum density, and E the energy density in the frame of the laboratory. The subscripts j,k stand for the dimensions, δ^jk is the Kronecker symbol, and h is the specific enthalpy given by $\mathit{h}\mathrm{=}\mathit{h}\mathrm{\left(}\mathit{P,\rho }\mathrm{\right)}\mathrm{=}\mathrm{1}\mathrm{+}\frac{\mathit{\gamma }}{\mathit{\gamma }\mathrm{-}\mathrm{1}}\frac{\mathit{P}}{\mathit{\rho }}\mathrm{=}\frac{\mathit{e}\mathrm{+}\mathit{P}}{\mathit{\rho }}\mathit{,}$(8)and ρ is the proper mass density, v_j the fluid three-velocity, P the gas pressure, and e is the sum of the internal energy and rest mass energy of the fluid. The Lorentz factor is given by $\mathrm{\Gamma }\mathrm{=}\frac{\mathrm{1}}{\sqrt{\mathrm{1}\mathrm{-}{\mathit{v}}^{\mathrm{2}}}}\mathrm{·}$(9)A passive scalar s can be included in the simulations using S = sρΓ as the conserved variable and F = ρsvΓ to compute its flux.

An equation of state closes the system of Eqs. (2)–(4). The most commonly used is the so-called classical equation of state, where the rest mass energy is removed from the total internal energy e: $\mathit{P}\mathrm{=}\mathrm{(}\mathit{\gamma }\mathrm{-}\mathrm{1}\mathrm{)}\mathrm{(}\mathit{e}\mathrm{-}\mathit{\rho }\mathrm{)}\mathit{,}$(10)where γ is the ratio of specific heats (or adiabatic index), which is constant and should not be confused with the Lorentz factor Γ. In the non-relativistic limit γ = 5/3, in the ultrarelativistic limit γ = 4/3. The sound speed is given by ${\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}\mathrm{\equiv }{\left(\frac{\mathit{\partial P}}{\mathit{\partial e}}\right)}_{\mathit{s}}\mathrm{=}\mathit{\gamma }\frac{\mathit{P}}{\mathit{h\rho }}\mathrm{·}$(11)This leads to c_s < 1/3 in the ultrarelativistic limit and c_s < 2/3 in the non-relativistic limit. The kinetic theory of relativistic gases (Synge 1957) shows that the ratio of specific heats cannot be kept constant when passing from non-relativistic to highly relativistic temperatures and provides an equation of state that is valid for all temperatures. Owing to the complexity of the expressions relating different thermodynamical quantities, approximations have been developed for the implementation in numerical schemes (Falle & Komissarov 1996; Mignone et al. 2005; Ryu et al. 2006; Mignone & McKinney 2007). This has not been implemented yet in the relativistic extension of RAMSES that we present.

The RHD equations have a similar structure to the Euler equations and can be solved following the same numerical method and making only localized changes. Still, the strong coupling of the equations through the presence of the Lorentz factor and the enthalpy makes the resolution of the equations more complex. Moreover, in relativistic dynamics, velocities are bounded by the speed of light, resulting in an important numerical constraint. These changes particularly affect the determination of the primitive variables, the development of a second-order scheme (Sect. 2.2), and the AMR framework (Sect. 2.3).

2.2. Second-order numerical scheme

Several steps in the algorithm require the use of the primitive variables $\begin{array}{ccc}{q}\mathrm{=}\left(\begin{array}{c}\\ \\ \\ \end{array}\right)\mathit{.}& & \end{array}$(12)In classical hydrodynamics, analytic relations allow a straightforward conversion from conservative to primitive variables. In RHD, there are no such relations, and several methods have been developed to recover the primitive variables (see e.g. Aloy & Martí 1999; Noble et al. 2006; Ryu et al. 2006). We use the method described in Mignone & McKinney (2007), which uses a Newton-Raphson algorithm to solve an equation on W′ = W − D = ρhΓ² − ρΓ. The details of the computation are given in Appendix A.1.

The determination of the fluxes F cells in Eq. (1) involves solving the Riemann problem at the interface between cells. Different solvers have been developed, and some have been extended to RHD. We have implemented the relativistic HLL and HLLC solver in RAMSES. For the HLL Riemann solver, the expression of the Godunov flux is identical to the non-relativistic expression (Schneider et al. 1993). However, the computation of the maximal wavespeeds propagating towards the left and right is different. In Newtonian hydrodynamics, the wavespeed is the sum of the sound speed and the advection speed of the flow. The relativistic composition of velocities couples the velocity of the flow parallel and perpendicular to the direction of spatial derivation and all components of the velocity need to be taken into (Del Zanna & Bucciantini 2002). An HLLC solver was developed by Mignone & Bodo (2005), whose method we closely follow.

To reach high enough accuracy in scientific applications, second-order schemes are, however, necessary. In such cases, fluxes are determined half a time step ahead of the current time step. We have implemented the monotonic upstream scheme for conservation laws (MUSCL) and piecewise linear (PLM) methods. In both schemes, the update is performed on the primitive variables, which are determined by a Taylor expansion $\begin{array}{ccc}{{q}}_{\mathit{i}\mathrm{\pm }\mathrm{1}\mathit{/}\mathrm{2}\mathit{,}\mathrm{L}\mathit{,}\mathrm{R}}^{\mathit{n}\hspace{0.17em}\mathrm{+}\hspace{0.17em}\mathrm{1}\mathit{/}\mathrm{2}}\mathrm{=}{{q}}_{\mathit{i}}^{\mathit{n}}\mathrm{\pm }\mathrm{\Delta }{{q}}_{\mathit{i}}^{\mathit{n}}\mathrm{+}\frac{\mathrm{d}{{q}}_{\mathit{i}}^{\mathit{n}}}{\mathrm{d}\mathit{t}}\frac{\mathrm{\Delta }\mathit{t}}{\mathrm{2}}\mathrm{·}& & \end{array}$(13)The subscripts L and R respectively stand for the left-and right-hand sides of a cell boundary, $\mathrm{\Delta }{{q}}_{\mathit{i}}^{\mathit{n}}$ are the slopes of the variables in the cell, computed using a slope limiter, which is not affected by special relativity, $\mathit{d}{{q}}_{\mathit{i}}^{\mathit{n}}$ is the temporal variation, given by $\begin{array}{ccc}\frac{\mathit{\partial }{q}}{\mathit{\partial t}}\mathrm{+}\sum _{\mathit{k}}^{\mathit{N}}{A}\frac{\mathit{\partial }{q}}{\mathit{\partial x}}\mathrm{=}\mathrm{0}\mathit{,}& & \end{array}$(14)where A = ∂F(q)/∂U(q) is the Jacobian matrix of the system. In classical hydrodynamics, flows in a given direction are not affected by motions perpendicular to that direction. In RHD all spatial directions are coupled through the Lorentz transformation: A is 5 × 5 matrix even for 1D flows (Pons et al. 2000). The value of A and its equivalents along the y and z directions are given in Appendix A.3.

In a MUSCL-Hancock scheme (Van Leer 1979), one uses $\mathit{d}{q}\mathrm{=}\mathrm{-}{A}\frac{\mathit{\partial }{q}}{\mathit{\partial x}}\mathrm{d}\mathit{t}$(15)directly to reconstruct the variables. In RAMSES, we have found that this reconstruction gives satisfactory results with the minmod (Roe 1986) slope limiter but fails with the moncen limiter (van Leer 1977). Therefore we implemented the PLM method, which works with both limiters.

In the Piecewise Linear Method (Collela 1990), one rewrites Eq. (14) $\frac{\mathit{\partial }{q}}{\mathit{\partial t}}\mathrm{+}\sum _{\mathit{k}}^{\mathit{N}}{{L}}^{\mathit{\alpha }}{\mathit{\lambda }}^{\mathit{\alpha }}{{R}}^{\mathit{\alpha }}\frac{\mathit{\partial }{q}}{\mathit{\partial x}}\mathrm{=}\mathrm{0}$(16)with λ^α the eigenvalues, and (L^α,R^α) the eigenvectors of the Jacobian matrix. In this case, the slopes are projected on the characteristic variables, using the left-hand eigenvectors. The waves that cannot reach the interface between cells in the given time step are filtered out and do not contribute to the interface states. One then recovers the primitive variables, multiplying by the right-hand eigenvectors. The final interface states are then computed by adding the contributions of all the selected waves: $\begin{array}{ccc}{{q}}_{\mathit{i,}\mathrm{R}}^{\mathit{n}\mathrm{+}\mathrm{1}\mathit{/}\mathrm{2}}& \mathrm{=}& {{q}}_{\mathit{i}}^{\mathit{n}}\mathrm{+}\sum _{{\mathit{\lambda }}_{\mathit{i}}^{\mathit{\alpha }}\hspace{0.17em}\mathit{<}\hspace{0.17em}\mathrm{0}}\left(\frac{\mathrm{1}}{\mathrm{2}}\left(\mathrm{-}\mathrm{1}\mathrm{-}{\mathit{\lambda }}_{\mathit{i}}^{\mathit{\alpha }}\frac{\mathrm{d}\mathit{t}}{\mathrm{d}\mathit{x}}\right){{L}}_{\mathit{i}}^{\mathit{\alpha }}\mathit{\partial }{{q}}_{\mathit{i}}\right)\mathrm{·}{{R}}_{\mathit{i}}^{\mathit{\alpha }}\\ {{q}}_{\mathit{i,}\mathrm{L}}^{\mathit{n}\mathrm{+}\mathrm{1}\mathit{/}\mathrm{2}}& \mathrm{=}& {{q}}_{\mathit{i}}^{\mathit{n}}\mathrm{+}\sum _{{\mathit{\lambda }}_{\mathit{i}}^{\mathit{\alpha }}\hspace{0.17em}\mathit{>}\hspace{0.17em}\mathrm{0}}\left(\frac{\mathrm{1}}{\mathrm{2}}\left(\mathrm{1}\mathrm{-}{\mathit{\lambda }}_{\mathit{i}}^{\mathit{\alpha }}\frac{\mathrm{d}\mathit{t}}{\mathrm{d}\mathit{x}}\right){{L}}_{\mathit{i}}^{\mathit{\alpha }}\mathit{\partial }{{q}}_{\mathit{i}}\right)\mathrm{·}{{R}}_{\mathit{i}}^{\mathit{\alpha }}\mathit{.}\end{array}$The full set of eigenvalues and eigenvectors is given in Appendix A.3.

Since all directions are reconstructed separately, although each component is subluminal, nothing guarantees that the norm of the total velocity remains subluminal. When a superluminal velocity is obtained, one possibility is to reconstruct Γv and renormalize the velocity when necessary (Aloy & Martí 1999) or to reconstruct the Lorentz factor independently (Wang et al. 2008). When superluminous velocities are obtained, we choose to switch to first-order reconstruction (q_{i    ±    1/2,L,R} = q_i) in all directions, for all variables. No switch occurred in our test simulation of a relativistic jet (see Sect. B.3). It typically occurs once every 10⁷ updates in our gamma-ray binary simulations with a Lorentz factor of 7, probably due to the much larger region covered by the high-velocity flow and its very low pressure (see Sect. 3.3).

Figure 1 shows the result of a shock tube test for the two different methods. Both methods give comparable results and show that the code works well on a uniform grid. In this test (Martí & Müller 2003), ρ_L = 10, ρ_R = 1, P_L = 13.3 and P_R = 10^-6 and there is no initial velocity. The flow is only mildly relativistic in the dynamical sense but e ≫ ρ in the right-hand state, resulting in a thermodynamically highly relativistic flow. Figure 2 illustrates the differences between a non-relativistic adiabatic index (γ = 5/3) and the ultrarelativistic adiabatic index (γ = 4/3). The compression ratio in the shocked region is higher than for Newtonian hydrodynamics, especially when the gas is thermodynamically ultrarelativistic. The results compare very well with Fig. 2 in Ryu et al. (2006).

Fig. 1 Density, velocity, and pressure in the laboratory frame at t = 0.45. The thin solid line shows the PLM method, the dashed line the MUSCL metod, which gives very similar results. The analytic solution (Giacomazzo & Rezzolla 2006) is given in blue. The resolution is nx = 256.

Fig. 2 Density, velocity, and pressure in the laboratory frame at t = 0.45 for tests with different adiabatic indices. The black solid line shows γ = 5/3, the dashed green line shows γ = 4/3, the resolution is nx = 256.

2.3. Implementation within the AMR structure

Mesh refinement is an important tool when modelling a wide range of spatial structures. Several methods exist and have been successfully associated with RHD schemes: static mesh refinement (Beckwith & Stone 2011), block-based AMR (Hughes et al. 2002; Zhang & MacFadyen 2006; Mignone et al. 2012) or tree-based AMR as in Keppens et al. (2012). Moving meshes offer another possibility to enable locally increased resolution (Duffell & MacFadyen 2011).

In RAMSES, cells are related in a recursive tree structure (Kravtsov et al. 1997) and gathered together in octs of 2^N cells (N is the number of dimensions), which share the same parent cell. When creating new refined cells or computing the fluxes at the interface between two levels, the conserved variables need to be interpolated from level l to level l + 1. The interpolation is done at second order, using a minmod limiter for the linear reconstruction. Conversely, since hydrodynamical updates are only performed at the highest level of refinement, variables need to be determined at lower levels by computing the average over the whole oct.

In RHD, the restriction step may lead to failures, when the energy is strongly dominated by the kinetic energy. In RHD, positive pressure and subluminal velocity are guaranteed when E² > M² + D² (Mignone & Bodo 2005, see Appendix A for the demonstration). Although cells at a given level satisfy this condition, nothing guarantees that ${\mathit{E}}_{\mathrm{oct}}^{\mathrm{2}}\mathit{>}{\mathit{M}}_{\mathrm{oct}}^{\mathrm{2}}\mathrm{+}{\mathit{D}}_{\mathrm{oct}}^{\mathrm{2}}$ where the subscript oct means variables are summed over an oct. The resulting state can be non-physical. We found that this problem can be bypassed by performing the reconstruction on the specific internal energy ϵ (i.e. the temperature) rather than on the total energy using $\mathit{ϵ}\mathrm{=}\frac{\mathit{P}}{\mathit{\rho }}\frac{\mathrm{1}}{\mathit{\gamma }\mathrm{-}\mathrm{1}}\mathit{,}$(19)where P and ρ are computed with the Newton-Raphson scheme used to computed the primitive variables. After the restriction, one recovers the total energy using $\mathit{E}\mathrm{=}{\mathrm{\Gamma }}^{\mathrm{2}}\mathit{\rho h}\mathrm{-}\mathit{P}\mathrm{=}\mathit{Dh}\mathrm{\Gamma }\mathrm{-}\mathrm{(}\mathit{\gamma }\mathrm{-}\mathrm{1}\mathrm{)}\frac{\mathit{D}}{\mathrm{\Gamma }}\mathit{ϵ,}$(20)where $\mathit{h}\mathrm{=}\mathrm{1}\mathrm{+}\mathit{\gamma ϵ} \mathrm{and} \mathrm{\Gamma }\mathrm{=}\sqrt{\frac{{\mathit{M}}^{\mathrm{2}}}{{\mathit{D}}^{\mathrm{2}}{\mathit{h}}^{\mathrm{2}}}\mathrm{+}\mathrm{1}}\mathit{.}$(21)This method makes the numerical scheme non-conservative but guarantees that the pressure is positive and the speed is subluminal.

Refinement is often based on gradients in the flow variables: either the density, pressure or velocity. In highly relativistic flows, the velocity does not vary significantly and the variations in the Lorentz factor may not be captured properly, thus strongly altering the dynamics of highly relativistic flows. To avoid this, we have implemented refinement based on the gradients of the Lorentz factor. Figure 3 shows two simulations of a shock tube test. We reproduced the test by Pons et al. (2000), with ρ_L = ρ_R = 1, P_L = 10³, and P_R = 10^-2, v_x,L = v_x,R = 0, v_y,L = v_y,R = 0.99 initially. This is a very stringent test, with a maximum Lorentz factor of 120. In the first case, refinements are only based on pressure gradients and do not properly capture the contact discontinuity and positions of the shock. In the other simulation, refinement is also based on the Lorentz factor. In this case, the simulation is in very good agreement with the analytic result. Similar accuracy can also be found by refining according to both the pressure and density gradients. However, it increases the computational cost by ≃15% in this test. This test validates the implementation within the AMR framework, which is well adapted for shocks in highly relativistic flows. We estimated the level in non-conservation by measuring the variation of the total density, momentum, and energy (in the laboratory frame) at the end of the test. We found relative errors of 3 × 10^-9,   9 × 10^-5, and − 6 × 10^-8, respectively.

Fig. 3 Density, parallel velocity, transverse velocity, and pressure in the laboratory frame for two different simulations. In both cases, the coarse grid is set by nx = 64 with 11 levels of refinement based on pressure gradients. The solid line shows a simulation where additional refinements occur based on Lorentz factor gradients. The result is superimposed to the analytic solution (Giacomazzo & Rezzolla 2006) given in blue. The red dashed line in the upper panel shows the refinement levels, based on the Lorentz factor gradients.

This method is currently implemented in RAMSES and succesfully passes all the commonly used numerical tests. A few of them are detailed in Appendix B.

3. 2D simulations of gamma-ray binaries

Pulsar winds have a Lorentz factor of about 10⁴ − 10⁶ (see e.g. Kirk et al. 2009, for a review), which is far beyond the reach of current computer power and numerical schemes. Most present multidimensional simulations model Lorentz factors of at most a few 10. Modelling higher Lorentz factors, especially in low-density or low-pressure environments, require highly sophisticated numerical methods, particularly in 2D and 3D, since spatial directions are strongly coupled. Up to now, there has been no extensive study of the impact of relativistic effects regarding the interaction of a pulsar wind with its surroundings.

The following simulations are demanding, not only because they need high resolution, but also because of the different timescales involved. The longest dynamical time is set by the speed of the stellar wind, while the time step in the simulation is limited by the speed of the pulsar wind. Because v_p/v_∗ ≃ 100, simulating a gamma-ray binaries takes roughly a hundred times longer than simulating colliding stellar winds. The simulations in this section typically took between 10 000 and 15 000 CPU hours to complete.

3.1. Collision with supersonic relativistic winds

We have compared the interaction between a pulsar wind and a stellar wind to the collision between two stellar winds descried in Paper I. Similarly to the non-relativistic case, one can determine the position of the contact discontinuity between the winds. The standoff point, which is situated at the intersection between the line-of-centres and the contact discontinuity, is given by the equation of momentum fluxes ${\mathit{\rho }}_{\mathit{p}}{\mathit{h}}_{\mathit{p}}{\mathrm{\Gamma }}_{\mathit{p}}^{\mathrm{2}}{\mathit{v}}_{\mathit{p}}^{\mathrm{2}}\mathrm{+}{\mathit{P}}_{\mathit{p}}\mathrm{=}{\mathit{\rho }}_{\mathit{s}}{\mathit{v}}_{\mathit{s}}^{\mathrm{2}}\mathrm{+}{\mathit{P}}_{\mathit{s}}\mathit{,}$(22)where the subscript s represents the variables in the stellar wind, and the variables with subscript p refer to the pulsar. We neglect thermal pressure in the winds, which gives ${\mathit{\rho }}_{\mathit{p}}{\mathrm{\Gamma }}_{\mathit{p}}^{\mathrm{2}}{\mathit{v}}_{\mathit{p}}^{\mathrm{2}}\mathrm{=}{\mathit{\rho }}_{\mathit{s}}{\mathit{v}}_{\mathit{s}}^{\mathrm{2}}\mathit{.}$(23)The mass loss rate Ṁ ≡ 4πr²Γρv is constant so the location of standoff point, where Eq. (22) is realized, is set by the dimensionless ratio ${\mathit{\eta }}_{\mathrm{rel}}\mathrm{=}\frac{\mathit{Ṁ}\mathit{p}{\mathrm{\Gamma }}_{\mathit{p}}{\mathit{v}}_{\mathit{p}}}{\dot{{\mathit{M}}_{\mathit{s}}}{\mathit{v}}_{\mathit{s}}}\mathrm{=}{\mathit{\eta }}_{\mathrm{cl}}\mathrm{\Gamma }$(24)where η_cl (Stevens et al. 1992) is the usual definition of the momentum flux ratio of the winds in the non-relativistic limit. These equations suggest that we can expect a similar structure to the one in colliding stellar winds, provided we define the momentum flux ratio of the winds using Eq. (24).

3.2. Numerical setup

We performed 2D numerical simulations of the interaction in a region extending up to four times the binary separation a. Our 2D setup is cylindrical, as described in Paper I. We neglected orbital motion in our simulations to enable comparisons with the analytic estimates and clarify the differences with our previous non-relativistic results (Paper I). Orbital motion can be neglected without affecting the dynamics on scales much smaller than the spiral stepsize S ≳ v_sP_orb (Paper II). For LS 5039, the gamma-ray binary with the shortest known period, the step of the spiral is at least S ≃ 4 AU, which about 20 times the binary separation (Moldón et al. 2012).

Since we want to study the impact of relativistic effects, we do not set the pulsar wind velocity to ≃c but keep it as a free parameter. We set η_rel and v_p and derive the pulsar’s mass loss rate $\begin{array}{ccc}\dot{{\mathit{M}}_{\mathit{p}}}\mathrm{=}{\mathit{\eta }}_{\mathrm{rel}}\frac{\mathit{Ṁ}\mathit{s}{\mathit{v}}_{\mathit{s}}}{{\mathit{v}}_{\mathit{p}}{\mathrm{\Gamma }}_{\mathit{p}}}\mathrm{·}& & \end{array}$(25)As for the classical case, the winds are initialized in masks following the method described in Lemaster et al. (2007). Density is determined by mass conservation in the winds, and the velocity is set to the terminal velocity. As in stellar wind simulations, we set the Mach number ℳ = v/c_s = 30 at the distance r = a. We introduce a passive scalar, which is set to one in the stellar wind and set to nought in the pulsar wind.

To ease comparisons between the different simulations, we set the adiabatic index γ to 5/3 and kept the parameters of the stellar wind identical in all simulations. We set Ṁ_s = 10^-7  M_⊙ yr^-1 and v_s = 3000  km s^-1, corresponding to typical values for winds from early type stars (Puls et al. 2008). At this stage, we use the HLL Riemann solver to numerically quench the development of the Kelvin-Helmholtz instability at the contact discontinuity between the winds. We use the MUSCL scheme, combined with the minmod slope limiter. In most simulations, the binary is at the centre of a region of size l_box = 8 a, the coarse grid is set by n_x = 64, and the refinement is based on density and Lorentz factor gradients.

We performed test simulations for η_rel = 1 at different resolutions to determine the level of refinement required to reach numerical accuracy. The Sod test problem (Sect. B.1) suggests that the higher the Lorentz factor of the flow, the higher the resolution needed to obtain a satisfactory simulation. We consider a simulation provides satisfactory results when the positions of the discontinuities do not vary with increasing resolution. The position of the contact discontinuity is set by finding where the passive scalar equals 0.5. The positions of the two shocks are given by the minimum and maximum value of the derivative of the pressure. Figure 4 shows the positions of the discontinuities for different levels of resolution, for β_p = v_p/c = 0.99. This corresponds to Γ_p = 7.08. For l_max ≥ 12, the positions do not vary anymore. In comparison, l_max = 9 is sufficient for the non-relativistic wind in Fig. 4. A much higher resolution is required to properly model the interaction with a relativistic flow.

Fig. 4 Position of both shocks and the contact discontinuity for ηrel = 1, βp = 0.99 for increasing resolution lmax = 10 (thin dotted line), 11 (thick dotted line), 12 (dashed line), 13 (solid line). The last two resolutions give the same result.

3.3. Geometry of the colliding wind region

Figure 5 shows the density map and the Lorentz factor map for simulations with η_rel = 0.1,   1,   10, and β_p = 0.99. We set the maximal resolution to l_max = 12 and checked that, for η_rel = 0.1 and 10, the positions of the discontinuities do not vary with higher resolution (see Sect. 3.2). The stellar wind is the dense flow on the right, the pulsar wind is located on the left. The analytic solution for the contact discontinuity is obtained by combining the relativistic definition of the momentum flux ratio η_rel (Eq. (24)) with the position of the contact discontinuity derived by Antokhin et al. (2004) for colliding stellar winds. The solution assumes a thin shell geometry where both shocks and contact discontinuity are merged in one single layer. The overall structure appears similar to what is found for non-relativistic flows (see e.g Fig. 5 in Paper I). However, relativistic effects have a non-negligible influence on the structure. For instance, the symmetry is broken between the cases η_rel = 0.1 and η_rel = 10, the latter showing a reconfinement shock while the former does not. The direction of the velocity in the winds is indicated on the map of the Lorentz factor. Along the line-of-centres, the winds collide head on, and when getting farther from the binary, the velocity is mostly parallel to the direction of the shocks. This corresponds to the shock tube test simulation (Sect. B.1), which required very high resolution because of the high velocity perpendicular to the shock normal (v_⊥ = 0.99). The importance of the velocity transverse to the shock normal is the reason high resolution is necessary in simulations with increasing Lorentz factors in the pulsar wind.

Fig. 5 Density (upper row) and Lorentz factor (lower row) maps for simulations with ηrel = 0.1,   1,   10 (from left to right) and βp = 0.99. The pulsar is located at x =  −0.5, the star at x = 0.5. The density is given in g cm-2, the arrows show the velocity field.

To quantify the impact of relativistic effects, we carried out simulations with η_rel = 1 for increasing values of the pulsar wind velocity β_p =  { 0.01,0.5,0.99,0.998 }. The corresponding Lorentz factors are Γ =  { 1.00005,1.15,7.08,   15.8 }. Figure 6 shows the positions of both shocks and the contact discontinuity in the different simulations. This plot shows that, on the line-of-centres, the location of the contact discontinuity is identical in all simulations and is set by η_rel. At the edges of the box, the contact discontinuity is not exactly midplane between the winds and is closer to the pulsar. We checked that, for Mach numbers above ten, the contact discontinuity stays midplane in non-relativistic simulations with η = 1 and a different Mach number for each wind. We interpret the tilt of the CD as the impact of thermal pressure in the winds. For a given momentum flux ratio, when the pulsar wind speed is increased, its pressure is decreased, because we keep a fixed Mach number, while the stellar wind is unchanged. Thus, the higher the Lorentz factor of the pulsar wind, the more the thermal pressure from the stellar wind dominates the pulsar wind pressure. The tilt of the CD is related to our definition of the pressure in the pulsar wind, which decreases with the Lorentz factor. When we set the pressure in the pulsar wind equal to the pressure in the stellar wind, both still being small compared to the ram pressure on the line-of-centres, we find that the CD is exactly on the axis of symmetry. We consider this case to be unrealistic because the pulsar wind pressure is expected to be negligible due to adiabatic expansion (Kirk et al. 2009). The shape of the shock in the pulsar wind changes with its Lorentz factor, especially far from the binary axis. In these zones, the velocity is mostly transverse to the shock normal. Unlike the non-relativistic case, the transverse and normal velocities are tied in the jump conditions via the Lorentz factor and this has an impact on the shock location (Taub 1948). The position of the shock in the stellar wind remains unaffected by the speed of the pulsar wind. We conclude that the shape of the interaction region is affected by relativistic effects, the reason being that the pressure and transverse speed influence the relativistic jump conditions through h and Γ. The numerical tests in Appendix B.2 highlight that very high resolution is needed to properly model the shock propagation speed, when there is a high velocity perpendicular to the shock propagation direction. Although we carefully determined an adequate resolution for our simulations, a numerical effect proper to relativistic simulations cannot be excluded.

Fig. 6 Position of both shocks and the contact discontinuity in a simulation with ηrel = 1, in simulations with different values for the velocity of the pulsar wind : βp= 0.01 (dotted line), 0.5 (dashed line), 0.99 (solid line) and 0.998 (thick solid line). For β = 0.998 we used a smaller computational box.

To evaluate the impact of the Lorentz factor on the hydrodynamics of the shocked region, we followed the hydrodynamical quantities along streamlines in the shocked pulsar wind. We find that, for all the values of β_p, starting from the shock, the flow reaccelerates up to its initial velocity. The acceleration tends to a spherical, adiabatic expansion, with 2πsv_pΓ_p ≡ kṀ_p ≃ constant (our simulations are 2D), where s is the curvilinear distance to the shock. Figure 7 shows the value of k for different speeds of the pulsar wind, for η_rel = 1. In all cases the shocked flow reaches k ≈ 3 on a scale comparable to the binary separation (increasing with β_p). For η_rel = 0.1, we find the same re-acceleration, although with a higher value of k ≃ 8. The higher corresponding mass flow is due to the smaller opening angle of the shocked pulsar wind.

Fig. 7 Mass flow (Ṁp = 2πsvpΓpρp) along a streamline in the shocked pulsar wind for βp = 0.01 (dotted line), 0.5 (dashed line), and 0.99 (solid line). The mass flow is normalized to the initial mass loss rate of the pulsar, we use ηrel = 1. The simulations with βp = 0.01 and 0.5 give the same results.

3.4. Kelvin-Helmholtz instability in gamma-ray binaries

The stellar wind has a velocity of 3000 km s^-1, while the pulsar wind is almost two orders of magnitude faster. Similarly to what happens for classical flows, at the contact discontinuity between the winds, the KHI is likely to modify the structure of the flow (Blandford & Pringle 1976; Turland & Scheuer 1976). Bodo et al. (2004) find analytic solutions to the dispersion relation in RHD and show that, in the frame of the laboratory, the stability criteria are the same as in the classical case, provided one uses the relativistic definition of the Mach number (ℳ_rel = ℳΓ/Γ_cs, where Γ_cs is the Lorentz factor of the sound speed). In 2D simulations, we thus expect the interface between the winds to be unstable provided ${\mathrm{ℳ}}_{\mathrm{rel}}\mathit{>}\sqrt{\mathrm{2}}$.

Fig. 8 Density (upper panel) and mixing (lower panel) maps for a simulation with ηrel = 1 and βp = 0.99.

To verify the impact of the KHI on gamma-ray binaries, we performed a test simulation with η_rel = 1 and β_p = 0.99. Radice & Rezzolla (2012) highlight the importance of an adapted Riemann solver and a low-diffusivity scheme in the study of the KHI. We used the HLLC Riemann solver to enable the development of the KHI at the contact discontinuity between the winds. This simulation displays a velocity difference (v_p/v_s ≃ 100). In Paper II, we found that in the classical limit, such a velocity difference leads to a significant amount of mixing between the wind and destroys the expected large scale spiral structure.

Figure 8 shows the density map and mixing between the stellar wind and the pulsar wind. The Kelvin-Helmholtz instability is clearly present in the whole colliding wind region. This is consistent with Fig.7 by Bosch-Ramon et al. (2012). In their simulation, Γ_p=10, η_rel = .3 and the stellar wind parameters are very close to ours, making comparison possible. Although orbital motion is included in their work, it has no impact at the location we are considering, close to the star. At a distance of about twice the binary separation, the re-acceleration of the pulsar wind means that the conditions in the shocked flow are incompatible with the analytic criterion for the developing the KHI. This suggests that the eddies present far away from the binary are due to the advection of eddies formed closer in. The eddies affect the position of the contact discontinuity but also modify the positions of both shocks. We find important differences in the position of the shock in the pulsar wind between two successive outputs, which indicates variation on a timescale of at most ten days. The KHI induces some mixing between the winds, although it occurs mostly in low-density parts of the shocked region.

The overall structure is similar to what we found for colliding non-relativistic stellar winds, for η = 1,v₁/v₂ = 20 (see Fig. 5 in Paper I) and in a test simulation with η = 1,v₁/v₂ = 100 (with the same resolution as in the relativistic simulation). Still in the classical limit, the mixed region covers a larger domain, with parts of the faster wind crossing the contact discontinuity and being totally embedded in the slower, denser wind. The mixed region seems smoother and less broken up in the classical limit. Comparison between both simulations is not straightforward because the density jump between both winds is about two orders of magnitude higher in the relativistic simulation, which probably decreases the growth of the intsability.

Bosch-Ramon et al. (2012) have simulated the large scale structure of gamma-ray binaries including orbital motion in 2D. Their simulations of the short-period binary LS 5039 show strong instabilities far from the binary, which they partly attribute to the accelerated transonic flow being deflected by the slower stellar wind and producing an unstable shocked structure. Whether this effect holds for higher values of the Lorentz factor and systems with longer orbital periods (such as PSR B1259-63) is worth investigating further.

4. Discussion: scaling to realistic pulsar winds

The simulations we performed have highlighted some similarities and differences between the interaction of classical stellar winds and the interaction between a stellar wind and a relativistic pulsar wind. This raises two questions:

• To what extent can simulations in the classical limit provideinformation on the structure of gamma-ray binaries?

• How can the results of relativistic simulations be extrapolated to more realistic values for the pulsar wind velocity (Γ_p ≃ 10^4 − 6) ?

Pulsar winds are characterized by the spindown power Ė_p carried by the wind, which can be estimated from timing of the pulse period. The corresponding momentum flux ratio can be expressed as a function of Ė_p by noticing that $\mathit{Ė}{}_{\mathit{p}}\mathrm{=}\mathrm{\Gamma }\dot{{\mathit{M}}_{\mathit{p}}}{\mathit{c}}^{\mathrm{2}}$ and v_p ≈ c so ${\mathit{\eta }}_{\mathrm{rel}}\mathrm{=}\frac{\mathit{Ė}\mathit{p}}{\mathit{Ṁ}\mathit{s}{\mathit{v}}_{\mathit{s}}\mathit{c}}\frac{{\mathit{v}}_{\mathit{p}}}{\mathit{c}}\mathrm{\approx }\frac{\mathit{Ė}\mathit{p}}{\mathit{Ṁ}\mathit{s}{\mathit{v}}_{\mathit{s}}\mathit{c}}\mathrm{·}$(26)Using this η_rel provides a good indication for the position of the contact discontinuity between the winds and may be used to rescale simulations of non-relativistic colliding stellar winds. Nonetheless, we found that the pulsar wind is more collimated than expected when we carried out relativistic simulations where the pulsar wind pressure P_p is lower than the stellar wind pressure P_s, even though both pressures are negligible compared to the ram pressure term ρΓ²v² on the line-of-centers (Eq. (22)). The difference probably occurs because the relativistic jump conditions are coupled by Γ and h, with the strongest impact in the shock wings where the speed is mostly transverse to the shock. In contrast, Bogovalov et al. (2008) find that the location of the contact discontinuity in their RHD simulations is identical to the classical case, with η calculated using Eq. (26). We suspect that they obtain this result because they set the pressure exactly to zero in both winds. We recover the classical result if we set the pressure to the same value in both winds.

Since a realistic pulsar wind has P_p ≈ 0 ≪ P_s ≠ 0, our results suggest that pulsar winds are more collimated than assumed by extrapolating from classical results. This may have a significant impact on the interpretation of the radio VLBI maps of gamma-ray binaries, where the collimated radio emission has been modelled as synchrotron emission from the shocked pulsar wind (Dhawan et al. 2006; Moldón et al. 2011a,b, 2012). In particular, Romero et al. (2007) notes that the collimation angles in the binary LSI +61°303 were too narrow compared to the ones expected in the classical limit. Simulations at higher values of the Lorentz factor are highly desirable to confirm this and fully explore the impact of the wind pressure on shock location.

To easy comparison with colliding stellar winds, we assumed a constant adiabatic index γ = 5/3 in our simulation. While this is realistic for the stellar wind and (cold) unshocked pulsar wind, the adiabatic index of the (warm) shocked pulsar wind is probably close to four thirds. Performing a test simulation with η_rel = 1, β_p = .99, and γ = 4/3 in both winds, we find that the shocked region in the pulsar wind is about 60% denser than when γ = 5/3, and slightly narrower. These results are consistent with what we observed in Fig. 2. However, the exact impact on gamma-ray binaries is not straightforward, since the adiabatic index transitions from the non-relativistic to the ultrarelativistic limit. As can be seen in shock tests in Figs. 5–7 of Ryu et al. (2006), when a relativistic equation of state is used, although the density jump is confined by the two limiting cases, the exact hydrodynamical structure is more complex than a simple intermediate between the cases γ = 4/3 and 5/3.

Simulations with moderate values for the Lorentz factors provide information on the density and velocity structure of the shocked flow that is useful for calculating emission properties of the system. We found that the shocked pulsar wind reaccelerates to its initial Lorentz factor and reaches an asymptotic regime with $\mathrm{4}\mathit{\pi }{\mathit{s}}^{\mathrm{2}}{\mathit{\rho }}_{\mathit{p}}^{\mathrm{\ast }}{\mathit{v}}_{\mathit{p}}^{\mathrm{\ast }}{\mathrm{\Gamma }}_{\mathit{p}}^{\mathrm{\ast }}\mathrm{=}\mathit{k}\dot{{\mathit{M}}_{\mathit{p}}}\mathrm{=}$ constant (or equivalently for our cylindrical 2D geometry $\mathrm{2}\mathit{\pi s}{\mathit{\rho }}_{\mathit{p}}^{\mathrm{\ast }}{\mathit{v}}_{\mathit{p}}^{\mathrm{\ast }}{\mathrm{\Gamma }}_{\mathit{p}}^{\mathrm{\ast }}\mathrm{=}\mathit{k}\dot{{\mathit{M}}_{\mathit{p}}}\mathrm{=}$ constant, with stars for quantities in the shocked region). We find that k is of order unity and is related to the opening angle of the shocked region and seems independent of the Lorentz factor of the pre-shock flow. Since the shocked wind accelerates up to its initial (i.e. pre-shock) velocity on a distance of a few times the binary separation, we have ${\mathit{v}}_{\mathit{p}}^{\mathrm{\ast }}\mathrm{=}{\mathit{v}}_{\mathit{p}}$, and this implies that, asymptotically, the density in the shocked wind is given by ${\mathit{\rho }}^{\mathrm{\ast }}\mathrm{\simeq }\frac{\mathit{Ė}\mathit{p}}{\mathrm{4}\mathit{\pi }{\mathit{s}}^{\mathrm{2}}\mathit{c}}\frac{\mathit{k}}{{\mathit{v}}_{\mathit{p}}^{\mathrm{2}}{\mathrm{\Gamma }}_{\mathit{p}}^{\mathrm{2}}}\mathit{,}$(27)where we have used $\mathit{Ė}{}_{\mathit{p}}\mathit{/}\mathit{c}\mathrm{\equiv }{\mathrm{\Gamma }}_{\mathit{p}}\mathit{Ṁ}{}_{\mathit{p}}\mathit{v}_{\mathit{p}}$.

We have verified this relation by measuring the density ratio in the shocked pulsar wind for the simulation with β_p = 0.5 and β_p = 0.99. Following Eq. (27), we expect ρ_0.5/ρ_0.99 ≃ 150. At the edge of our simulation domain, where s ≃ 4a, we find that ρ_0.5/ρ_0.99 varies between 140 and 155. This method assumes the flow has reached its terminal velocity. It can be used to estimate the density in the shocked pulsar wind for any value of its Lorentz factor, even using simulations in the non-relativistic limit. However, it is unable to precisely indicate the properties of the flow close to the binary. Similar caveats were found by Komissarov & Falle (1996), who proposed scaling relations to model non-relativistic equivalents of relativistic jets using the mass, energy or momentum flux. Equating the momentum flux leads to the most satisfactory results on the global morphology of the jet but still underestimates the size of the region with the highest Lorentz factors (Rosen et al. 1999).

The above ideal picture for the shocked flow is likely to be affected by the development of instabilities. Analytic work and simulations indicate that highly relativistic flows are less subject to the Kelvin-Helmholtz instability. Martí et al. (1997); Perucho et al. (2004b,a) indicate that cold jets with high Lorentz factors are more stable than slow, warm jets. In fast flows the timescale for the growth of the KHI is too long with respect to the dynamical time, while in cold flows the sound speed is low and perturbations propagate too slowly. In those particular systems, the specific axisymmetric geometry plays a role in the development of particular modes of the instability and makes it hard to directly remap the results to gamma-ray binaries.

Our simulations show that the KHI affects the positions of the shocks and induces some mixing between the winds. The most favourable conditions for developing the KHI occur close to the binary, where the shocked pulsar wind is subsonic and only mildly relativistic. Since the downstream velocity cannot be higher than c/3 for ultrarelativistic flows, the conditions for the developing the KHI are likely to be met in real systems (having Γ_p ≃ 10^4 − 6), at least close to the binary. The presence of the KHI further may be due to advection alone. Simulations at higher Lorentz factors are necessary to determine whether the unstable region is large enough to enable important growth before the advection phase begins. This is a key aspect since strong instabilities will induce mixing between the stellar and pulsar wind, which is likely to thermalize the high-energy particles of the pulsar wind and thus decrease their high-energy emission.

5. Conclusion and perspectives

We have developed a special relativistic extension to the adaptive-mesh refinement code RAMSES. We implemented two different second-order schemes and Riemann solvers. So far, only the so-called classical equation of state has been implemented, but the code was written in order to make the implementation of more complex equations of state straightforward. The code is fully three-dimensional and allows using adaptive mesh refinement. To model high-velocity flows, we included the possibility of refining cells according to the gradient of the Lorentz factor. Our code passes all of the common numerical tests in 1D, 2D, and 3D, including the most stringent tests with high Lorentz factors and velocities perpendicular to the flow.

Using this new relativistic code, we performed a 2D study of the geometry of the gamma-ray binaries, with Lorentz factors up to 16. The relativistic nature of the pulsar wind changes the momentum balance between the winds. The position of the shock in the pulsar wind is affected by the impact of the pre-shock wind pressure and transverse velocity on the shock. We found this makes the shocked pulsar wind more collimated than expected using classical hydrodynamics. This relativistic effect is strongest where the flow is essentially along the shock, far from the stagnation point. Comparing simulations with increasing values for the Lorentz factor of the pulsar wind, we find that the shocked pulsar wind re-accelerates up to its initial velocity on a scale of a few times the binary separation. We determine a prescription for the density in the shocked pulsar wind, which can be used for any value of the Lorentz factor. We find that the Kelvin-Helmholtz instability develops at the contact discontinuity affects the structure of the shocked region. Still, the exact understanding of the instability at realistic values of the pulsar wind Lorentz factor, in the particular geometry of gamma-ray binaries still needs to be confirmed. Our simulations provide clues do understanding the physics of gamma-ray binaries. However, detailed comparison with observations requires including more complex physical effects.

This code will be part of the next public release of RAMSES. It is well suited to study many astrophysical flows including gamma-ray bursts, pulsar wind nebulae, and jets in galactic binaries or active galactic nuclei. As an example, the Lorentz factor inferred from the prompt emission of gamma-ray bursts is of order 100, while the afterglow emission is related to a bulk motion of Γ ≃ 10 (see Gehrels et al. 2009, for a review) for a review). Following the fireball model, such systems display a very wide range of lengthscales, from the progenitor to decelaration radius of the forward shock, including the reverse shock and the structure of the internal shocks. Although the fully consistent study of such systems is still out of reach, the possibility for adaptive mesh refinement makes RAMSES a useful tool to model gamma-ray bursts. Our work on gamma-ray binaries can be easily adapted to simulating pulsar winds interacting with the interstellar medium or a supernova remnant. In such cases, the shock occurs further from the pulsar than in gamma-ray binaries, and the typical size of the interaction region is a few pc. This increase in length scales makes the use of AMR even more fundamental when modelling the structure and instabilities in pulsar wind nebulae.

Online material

Appendix A: Numerical scheme for RAMSES-RHD

This Appendix provides the Jacobian matrices and their eigenstructure for reconstructing the primitive variables in the 3D-RHD case. It also explains the recovery of the primitive variables.

Appendix A.1: Recovering the primitive variables

To determine the primitive variables (ρ,v_j,P)^T, we simplify the method by Mignone & McKinney (2007), initially developed for relativistic magnetohydrodynamics. One rewrites the total energy $\mathit{E}\mathrm{=}{\mathit{W}}^{\mathrm{\prime }}\mathrm{+}\mathit{D}\mathrm{-}\mathit{P}$(A.1)and solves an equation on W′ = W − D = ρhΓ² − ρΓ.

To avoid numerical problems in the non-relativistic or ultrarelativistic limits, one introduces ${\mathit{u}}^{\mathrm{2}}\mathrm{=}\frac{{\mathit{M}}^{\mathrm{2}}}{\mathrm{(}{\mathit{W}}^{\mathrm{\prime }}\mathrm{+}\mathit{D}{\mathrm{)}}^{\mathrm{2}}\mathrm{-}{\mathit{M}}^{\mathrm{2}}}\mathrm{=}{\mathrm{\Gamma }}^{\mathrm{2}}\mathit{v}2$(A.2)such that the Lorentz factor is given by $\mathrm{\Gamma }\mathrm{=}\mathrm{(}\mathrm{1}\mathrm{+}{\mathit{u}}^{\mathrm{2}}{\mathrm{)}}^{\mathrm{1}\mathit{/}\mathrm{2}}\mathit{.}$(A.3)This gives ${\mathit{W}}^{\mathrm{\prime }}\mathrm{=}\frac{\mathit{D}{\mathit{u}}^{\mathrm{2}}}{\mathrm{\Gamma }\mathrm{+}\mathrm{1}}\mathrm{+}{\mathrm{\Gamma }}^{\mathrm{2}}\frac{\mathit{\gamma }}{\mathit{\gamma }\mathrm{-}\mathrm{1}}\mathit{P,}$(A.4)which can be used to replace P in Eq. (A.1). Once W′ is found using, say, a Newton-Raphson method, one recovers the Lorentz factor using Eqs. (A.2) and (A.3), then the density ρ = D/Γ and finally pressure using Eq. (A.4). We determine W′ with a precision of 10^-10.

The method can be easily adapted to different equations of state, by simply changing Eq. (A.4). It is well suited to codes with AMR since an initial guess of the value of W₀, which guarantees the positivity of pressure, can be found analytically (Mignone & McKinney 2007). This avoids storing values of W′ between time steps, which can be cumbersome on a changing grid. One can determine W₀, which then gives ${\mathit{W}}_{\mathrm{0}}^{\mathrm{\prime }}\mathrm{=}{\mathit{W}}_{\mathrm{0}}\mathrm{-}\mathit{D}$, rewriting the pressure as $\mathit{P}\mathrm{=}\mathit{W}\mathrm{-}\mathit{E}\mathrm{=}\frac{{\mathit{M}}^{\mathrm{2}}\mathrm{-}{\mathit{v}}^{\mathrm{2}}{\mathit{W}}^{\mathrm{2}}\mathrm{+}\mathrm{4}{\mathit{W}}^{\mathrm{2}}\mathrm{-}\mathrm{4}\mathit{EW}}{\mathrm{4}\mathit{W}}$(A.5)since M² ≡ W²v². Pressure is thus guaranteed to be positive when ${\mathit{M}}^{\mathrm{2}}\mathrm{-}{\mathit{v}}^{\mathrm{2}}{\mathit{W}}^{\mathrm{2}}\mathrm{+}\mathrm{4}{\mathit{W}}^{\mathrm{2}}\mathrm{-}\mathrm{4}\mathit{EW}\mathit{>}\mathrm{0.}$(A.6)This happens for any W > W₊, with W₊ the largest of the two roots of the quadratic equation Eq. (A.6). As v < 1, solving Eq. (A.6) for v = 1, implies that the corresponding root W_+ ,1 > W_+ ,v, which guarantees positive pressure. One thus uses W_+ ,1 as the initial guess to start the Newton-Raphson scheme.

Provided the vector of conservative variables is physical, this method converges to a physical primitive state. The positivity of the density can be verified by D > 0. Subluminal velocity is guaranteed by E² > M² + D². Indeed, using Eq. (6) one has ${\mathit{E}}^{\mathrm{2}}\mathit{>}{\mathit{M}}^{\mathrm{2}}\mathrm{+}{\mathit{D}}^{\mathrm{2}}\mathrm{\Rightarrow }{\left(\frac{\mathit{M}}{\mathit{v}}\mathrm{-}\mathit{P}\right)}^{\mathrm{2}}\mathit{>}{\mathit{M}}^{\mathrm{2}}\mathrm{+}{\mathit{D}}^{\mathrm{2}}\mathit{.}$(A.7)This necessarily implies positive pressure and v < 1. Whenever a non-physical conserved state occurs, we floor the density to ρ = 10^-10 and the pressure to P = 10^-20. This occurs roughly once every 10⁶ updates in our gamma-ray binary simulations. It occurs in the unshocked pulsar wind, where the density and pressure are lowest and the Lorentz factor the highest. It mostly occurs along the line of centers of the binary, since the Cartesian grid is unable to perfectly model the spherical symmetry of the wind. The number of failures strongly reduces with resolution. Beckwith & Stone (2011) note a few failures when pressure is very low and suggest using the entropy instead of the energy as conserved variable in those cases. Reconstruction by using the neighbouring cells is very cumbersome owing to the AMR structure in RAMSES, so we have not considered this option.

Appendix A.2: Computing the time step

The time step is determined by the Courant condition for unsplit schemes $\mathrm{\Delta }\mathit{t}\mathrm{=}{\mathit{C}}_{\mathrm{CFL}}\mathrm{\times }\begin{array}{c}\min \end{array}\left(\frac{\mathrm{\Delta }\mathit{x}}{{\mathit{\lambda }}_{\mathit{x}}\mathrm{+}{\mathit{\lambda }}_{\mathit{y}}\mathrm{+}{\mathit{\lambda }}_{\mathit{z}}}\right)$(A.8)where C_CFL is the Courant number and λ_k is the propagation speed along direction k determined by ${\mathit{\lambda }}_{\mathit{k}}\mathrm{=}\begin{array}{c}\max \end{array}\mathrm{\left(}\mathrm{\right|}{\mathit{\lambda }}_{\mathrm{+}}\mathrm{\left|}\mathit{,}\mathrm{\right|}{\mathit{\lambda }}_{\mathrm{-}}\mathrm{\left|}\mathrm{\right)}$(A.9)with ${\mathit{\lambda }}_{\mathrm{\pm }}\mathrm{=}\frac{{\mathit{v}}_{\mathrm{\parallel }}\mathrm{(}\mathrm{1}\mathrm{-}{{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}}^{\mathrm{)}}\mathrm{+}\sqrt{\left(\mathrm{1}\mathrm{-}{\mathit{v}}^{\mathrm{2}}\right)\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{v}}_{\mathrm{\parallel }}^{\mathrm{2}}\mathrm{-}{\mathit{v}}_{\mathrm{\perp }}{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}}\mathrm{)}}{\mathrm{1}\mathrm{-}{\mathit{v}}^{\mathrm{2}}{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}}$(A.10)where v_∥ and v_⊥ are the velocity parallel and perpendicular to the direction k. All the simulations in this paper were performed using C_CFL = .8.

Our determination of the time step is more restrictive than the more commonly used (see e.g. Mignone et al. 2012) $\mathrm{\Delta }\mathit{t}\mathrm{=}{\mathit{C}}_{\mathrm{CFL}}\mathrm{\times }\begin{array}{c}\min \end{array}\left(\frac{\mathrm{\Delta }\mathit{x}}{{\mathit{\lambda }}_{\mathit{x}}}\mathit{,}\frac{\mathrm{\Delta }\mathit{x}}{{\mathit{\lambda }}_{\mathit{y}}}\mathit{,}\frac{\mathrm{\Delta }\mathit{x}}{{\mathit{\lambda }}_{\mathit{z}}}\right)\mathrm{·}$(A.11)Using the above expression for the time step, we find that the simulation of the 3D jet (Appendix B.3) shows no failure for C_CFL ≤ .5, which is comparable to other codes (Zhang & MacFadyen 2006).

Appendix A.3: Jacobian matrices of the 3D-RHD equations

This section will only be available online. One has $\begin{array}{ccc}\mathit{\partial }{q}& \mathrm{=}& \mathrm{-}{{A}}_{\mathit{x}}\frac{\mathit{\partial }{q}}{\mathit{\partial x}}\mathrm{d}\mathit{t}\mathrm{-}{{A}}_{\mathit{y}}\frac{\mathit{\partial }{q}}{\mathit{\partial y}}\mathrm{d}\mathit{t}\mathrm{-}{{A}}_{\mathit{z}}\frac{\mathit{\partial }{q}}{\mathit{\partial z}}\mathrm{d}\mathit{t}\mathit{.}\end{array}$(A.12)Following Font et al. (1994) for a 1D flow with transverse velocity and the general definition of the enthalpy $\begin{array}{ccc}\frac{\mathit{\partial h}}{\mathit{\partial x}}& \mathrm{=}& \frac{\mathit{\partial h}}{\mathit{\partial \rho }}\frac{\mathit{\partial \rho }}{\mathit{\partial x}}\mathrm{+}\frac{\mathit{\partial h}}{\mathit{\partial P}}\frac{\mathit{\partial P}}{\mathit{\partial x}}\mathrm{\equiv }\mathit{\chi }\frac{\mathit{\partial \rho }}{\mathit{\partial x}}\mathrm{+}\mathit{\kappa }\frac{\mathit{\partial P}}{\mathit{\partial x}}\mathrm{·}\end{array}$(A.13)With these definitions, the sound speed can be written as $\begin{array}{ccc}{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}\mathrm{=}\frac{\mathit{\rho \chi }}{\mathit{h}\left(\mathrm{1}\mathrm{-}\mathit{\rho \kappa }\right)}\mathrm{·}& & \end{array}$(A.14)For the classical equation of state, one has $\begin{array}{ccc}\mathit{\chi }& \mathrm{=}& \mathrm{-}\frac{\mathit{\gamma }}{\mathit{\gamma }\mathrm{-}\mathrm{1}}\frac{\mathit{P}}{{\mathit{\rho }}^{\mathrm{2}}}\\ \mathit{\kappa }& \mathrm{=}& \frac{\mathit{\gamma }}{\mathit{\gamma }\mathrm{-}\mathrm{1}}\frac{\mathrm{1}}{\mathit{\rho }}\mathrm{·}\end{array}$The complete Jacobian matrix is given by ${{A}}_{\mathit{x}}\mathrm{=}\left(\begin{array}{c}\\ {\mathit{c}}_{\mathrm{11}}& {\mathit{c}}_{\mathrm{12}}& {\mathit{c}}_{\mathrm{13}}& {\mathit{c}}_{\mathrm{14}}& {\mathit{c}}_{\mathrm{15}}\\ {\mathit{c}}_{\mathrm{21}}& {\mathit{c}}_{\mathrm{22}}& {\mathit{c}}_{\mathrm{23}}& {\mathit{c}}_{\mathrm{24}}& {\mathit{c}}_{\mathrm{25}}\\ {\mathit{c}}_{\mathrm{31}}& {\mathit{c}}_{\mathrm{32}}& {\mathit{c}}_{\mathrm{33}}& {\mathit{c}}_{\mathrm{34}}& {\mathit{c}}_{\mathrm{35}}\\ {\mathit{c}}_{\mathrm{41}}& {\mathit{c}}_{\mathrm{42}}& {\mathit{c}}_{\mathrm{43}}& {\mathit{c}}_{\mathrm{44}}& {\mathit{c}}_{\mathrm{45}}\\ {\mathit{c}}_{\mathrm{51}}& {\mathit{c}}_{\mathrm{52}}& {\mathit{c}}_{\mathrm{53}}& {\mathit{c}}_{\mathrm{54}}& {\mathit{c}}_{\mathrm{55}}\end{array}\right)\mathit{.}$(A.17)$\begin{array}{ccc}\mathit{N}{\mathit{c}}_{\mathrm{11}}& \mathrm{=}& {\mathit{v}}_{\mathit{x}}\mathit{N}\\ \mathit{N}{\mathit{c}}_{\mathrm{12}}& \mathrm{=}& \mathit{\rho h}{\mathrm{\Gamma }}^{\mathrm{2}}\mathrm{(}\mathit{\rho \kappa }\mathrm{-}\mathrm{1}\mathrm{)}\\ \mathit{N}{\mathit{c}}_{\mathrm{13}}& \mathrm{=}& \mathrm{0}\\ \mathit{N}{\mathit{c}}_{\mathrm{14}}& \mathrm{=}& \mathrm{0}\\ \mathit{N}{\mathit{c}}_{\mathrm{15}}& \mathrm{=}& {\mathit{v}}_{\mathit{x}}\mathrm{(}\mathrm{1}\mathrm{-}\mathit{\rho \kappa }\mathrm{)}\\ \mathit{N}{\mathit{c}}_{\mathrm{21}}& \mathrm{=}& \mathrm{0}\\ \mathit{N}{\mathit{c}}_{\mathrm{22}}& \mathrm{=}& {\mathit{v}}_{\mathit{x}}{\mathrm{\Gamma }}^{\mathrm{2}}\mathrm{(}\mathit{\rho h\kappa }\mathrm{-}\mathit{h}\mathrm{+}\mathit{\rho \chi }\mathrm{)}\\ \mathit{N}{\mathit{c}}_{\mathrm{23}}& \mathrm{=}& \mathrm{0}\\ \mathit{N}{\mathit{c}}_{\mathrm{24}}& \mathrm{=}& \mathrm{0}\\ \mathit{N}{\mathit{c}}_{\mathrm{25}}& \mathrm{=}& \frac{\mathrm{1}}{\mathit{\rho h}{\mathrm{\Gamma }}^{\mathrm{2}}}\mathrm{(}\mathit{\rho h}{\mathrm{\Gamma }}^{\mathrm{2}}\mathit{\kappa }\mathrm{(}\mathrm{1}\mathrm{-}{{\mathit{v}}_{\mathit{x}}^{\mathrm{2}}}^{\mathrm{)}}\mathrm{-}\mathit{h}\mathrm{-}\mathit{h}{\mathrm{\Gamma }}^{\mathrm{2}}\mathrm{(}{\mathit{v}}_{\mathit{y}}^{\mathrm{2}}\mathrm{+}{{\mathit{v}}_{\mathit{z}}^{\mathrm{2}}}^{\mathrm{)}}\mathrm{+}\mathit{\rho \chi }{\mathrm{\Gamma }}^{\mathrm{2}}{\mathrm{(}{\mathit{v}}_{\mathit{y}}^{\mathrm{2}}\mathrm{+}{{\mathit{v}}_{\mathit{z}}^{\mathrm{2}}}^{\mathrm{)}}}^{\mathrm{)}}\\ \mathit{N}{\mathit{c}}_{\mathrm{31}}& \mathrm{=}& \mathrm{0}\\ \mathit{N}{\mathit{c}}_{\mathrm{32}}& \mathrm{=}& {\mathit{v}}_{\mathit{y}}\mathit{\rho \chi }\\ \mathit{N}{\mathit{c}}_{\mathrm{33}}& \mathrm{=}& {\mathit{v}}_{\mathit{x}}\mathit{N}\\ \mathit{N}{\mathit{c}}_{\mathrm{34}}& \mathrm{=}& \mathrm{0}\\ \mathit{N}{\mathit{c}}_{\mathrm{35}}& \mathrm{=}& \mathrm{-}\frac{{\mathit{v}}_{\mathit{x}}{\mathit{v}}_{\mathit{y}}\mathrm{(}\mathrm{-}\mathit{h}\mathrm{+}\mathit{\rho \chi }\mathrm{+}\mathit{\rho h\kappa }\mathrm{)}}{\mathit{\rho h}}\\ \mathit{N}{\mathit{c}}_{\mathrm{41}}& \mathrm{=}& \mathrm{0}\\ \mathit{N}{\mathit{c}}_{\mathrm{42}}& \mathrm{=}& {\mathit{v}}_{\mathit{z}}\mathit{\rho \chi }\\ \mathit{N}{\mathit{c}}_{\mathrm{43}}& \mathrm{=}& \mathrm{0}\\ \mathit{N}{\mathit{c}}_{\mathrm{44}}& \mathrm{=}& {\mathit{v}}_{\mathit{x}}\mathit{N}\\ \mathit{N}{\mathit{c}}_{\mathrm{45}}& \mathrm{=}& \mathrm{-}\frac{{\mathit{v}}_{\mathit{x}}{\mathit{v}}_{\mathit{z}}\mathrm{(}\mathrm{-}\mathit{h}\mathrm{+}\mathit{\rho \chi }\mathrm{+}\mathit{\rho h\kappa }\mathrm{)}}{\mathit{\rho h}}\\ \mathit{N}{\mathit{c}}_{\mathrm{51}}& \mathrm{=}& \mathrm{0}\\ \mathit{N}{\mathit{c}}_{\mathrm{52}}& \mathrm{=}& \mathrm{-}{\mathit{\rho }}^{\mathrm{2}}\mathit{h}{\mathrm{\Gamma }}^{\mathrm{2}}\mathit{\chi }\\ \mathit{N}{\mathit{c}}_{\mathrm{53}}& \mathrm{=}& \mathrm{0}\\ \mathit{N}{\mathit{c}}_{\mathrm{54}}& \mathrm{=}& \mathrm{0}\\ \mathit{N}{\mathit{c}}_{\mathrm{55}}& \mathrm{=}& {\mathit{v}}_{\mathit{x}}{\mathrm{\Gamma }}^{\mathrm{2}}\mathrm{(}\mathit{\rho \chi }\mathrm{-}\mathit{h}\mathrm{+}\mathit{h\rho \kappa }\mathrm{)}\end{array}$with N = ⁽ − h + hρκ + ρχv²⁾Γ².

The matrix A_y is given by $\begin{array}{ccc}\mathit{N}{\mathit{c}}_{\mathrm{11}}& \mathrm{=}& {\mathit{v}}_{\mathit{y}}\mathit{N}\\ \mathit{N}{\mathit{c}}_{\mathrm{12}}& \mathrm{=}& \mathrm{0}\\ \mathit{N}{\mathit{c}}_{\mathrm{13}}& \mathrm{=}& \mathit{\rho h}{\mathrm{\Gamma }}^{\mathrm{2}}\mathrm{(}\mathit{\rho \kappa }\mathrm{-}\mathrm{1}\mathrm{)}\\ \mathit{N}{\mathit{c}}_{\mathrm{14}}& \mathrm{=}& \mathrm{0}\\ \mathit{N}{\mathit{c}}_{\mathrm{15}}& \mathrm{=}& {\mathit{v}}_{\mathit{y}}\mathrm{(}\mathrm{-}\mathit{\rho \kappa }{\mathrm{\Gamma }}^{\mathrm{2}}\mathrm{+}\mathit{\rho \kappa }{\mathrm{\Gamma }}^{\mathrm{2}}{\mathit{v}}^{\mathrm{2}}\mathrm{+}{\mathrm{1}}^{\mathrm{)}}\\ \mathit{N}{\mathit{c}}_{\mathrm{21}}& \mathrm{=}& \mathrm{0}\\ \mathit{N}{\mathit{c}}_{\mathrm{22}}& \mathrm{=}& {\mathit{v}}_{\mathit{y}}\mathit{N}\\ \mathit{N}{\mathit{c}}_{\mathrm{23}}& \mathrm{=}& {\mathit{v}}_{\mathit{x}}\mathit{\rho \chi }\\ \mathit{N}{\mathit{c}}_{\mathrm{24}}& \mathrm{=}& \mathrm{0}\end{array}$$\begin{array}{ccc}& & \\ \mathit{N}{\mathit{c}}_{\mathrm{25}}& \mathrm{=}& \frac{\mathrm{-}\mathrm{\left(}{\mathit{v}}_{\mathit{x}}{\mathit{v}}_{\mathit{y}}\mathrm{\right(}\mathrm{-}\mathit{h}\mathrm{+}\mathit{\rho \chi }\mathrm{+}\mathit{h\rho \kappa }\mathrm{\left)}\mathrm{\right)}}{\mathit{\rho h}}\\ \mathit{N}{\mathit{c}}_{\mathrm{31}}& \mathrm{=}& \mathrm{0}\\ \mathit{N}{\mathit{c}}_{\mathrm{32}}& \mathrm{=}& \mathrm{0}\\ \mathit{N}{\mathit{c}}_{\mathrm{33}}& \mathrm{=}& {\mathit{v}}_{\mathit{y}}{\mathrm{\Gamma }}^{\mathrm{2}}\mathrm{(}\mathit{\rho \chi }\mathrm{-}\mathit{h}\mathrm{+}\mathit{h\rho \kappa }\mathrm{)}\\ \mathit{N}{\mathit{c}}_{\mathrm{34}}& \mathrm{=}& \mathrm{0}\\ \mathit{N}{\mathit{c}}_{\mathrm{35}}& \mathrm{=}& \frac{\mathrm{1}}{\mathit{\rho h}{\mathrm{\Gamma }}^{\mathrm{2}}}\mathrm{(}\mathit{\rho h}{\mathrm{\Gamma }}^{\mathrm{2}}\mathit{\kappa }\mathrm{(}\mathrm{1}\mathrm{-}{{\mathit{v}}_{\mathit{y}}^{\mathrm{2}}}^{\mathrm{)}}\mathrm{-}\mathit{h}\mathrm{-}\mathit{h}{\mathrm{\Gamma }}^{\mathrm{2}}\mathrm{(}{\mathit{v}}_{\mathit{x}}^{\mathrm{2}}\mathrm{+}{{\mathit{v}}_{\mathit{z}}^{\mathrm{2}}}^{\mathrm{)}}\mathrm{+}\mathit{\rho \chi }{\mathrm{\Gamma }}^{\mathrm{2}}{\mathrm{(}{\mathit{v}}_{\mathit{x}}^{\mathrm{2}}\mathrm{+}{{\mathit{v}}_{\mathit{z}}^{\mathrm{2}}}^{\mathrm{)}}}^{\mathrm{)}}\\ \mathit{N}{\mathit{c}}_{\mathrm{41}}& \mathrm{=}& \mathrm{0}\\ \mathit{N}{\mathit{c}}_{\mathrm{42}}& \mathrm{=}& \mathrm{0}\\ \mathit{N}{\mathit{c}}_{\mathrm{43}}& \mathrm{=}& {\mathit{v}}_{\mathit{z}}\mathit{\rho \chi }\\ \mathit{N}{\mathit{c}}_{\mathrm{44}}& \mathrm{=}& {\mathit{v}}_{\mathit{y}}\mathit{N}\\ \mathit{N}{\mathit{c}}_{\mathrm{45}}& \mathrm{=}& \frac{\mathrm{-}{\mathit{v}}_{\mathit{x}}{\mathit{v}}_{\mathit{z}}\mathrm{(}\mathrm{-}\mathit{h}\mathrm{+}\mathit{\rho \chi }\mathrm{+}\mathit{\rho h\kappa }\mathrm{)}}{\mathit{\rho h}}\\ \mathit{N}{\mathit{c}}_{\mathrm{51}}& \mathrm{=}& \mathrm{0}\\ \mathit{N}{\mathit{c}}_{\mathrm{52}}& \mathrm{=}& \mathrm{0}\\ \mathit{N}{\mathit{c}}_{\mathrm{53}}& \mathrm{=}& \mathrm{-}{\mathit{\rho }}^{\mathrm{2}}\mathit{h}{\mathrm{\Gamma }}^{\mathrm{2}}\mathit{\chi }\\ \mathit{N}{\mathit{c}}_{\mathrm{54}}& \mathrm{=}& \mathrm{0}\\ \mathit{N}{\mathit{c}}_{\mathrm{55}}& \mathrm{=}& {\mathit{v}}_{\mathit{y}}{\mathrm{\Gamma }}^{\mathrm{2}}\mathrm{(}\mathit{\rho \chi }\mathrm{-}\mathit{h}\mathrm{+}\mathit{h\rho \kappa }\mathrm{)}\mathit{.}\end{array}$Along z one has A_z given by $\begin{array}{ccc}\mathit{N}{\mathit{c}}_{\mathrm{11}}& \mathrm{=}& {\mathit{v}}_{\mathit{z}}\mathit{N}\\ \mathit{N}{\mathit{c}}_{\mathrm{12}}& \mathrm{=}& \mathrm{0}\\ \mathit{N}{\mathit{c}}_{\mathrm{13}}& \mathrm{=}& \mathrm{0}\\ \mathit{N}{\mathit{c}}_{\mathrm{14}}& \mathrm{=}& \mathit{h\rho }{\mathrm{\Gamma }}^{\mathrm{2}}\mathrm{(}\mathit{\rho \kappa }\mathrm{-}\mathrm{1}\mathrm{)}\\ \mathit{N}{\mathit{c}}_{\mathrm{15}}& \mathrm{=}& {\mathit{v}}_{\mathit{z}}\mathrm{(}\mathrm{-}\mathit{\rho }{\mathrm{\Gamma }}^{\mathrm{2}}\mathit{\kappa }\mathrm{+}\mathit{\rho }{\mathit{v}}^{\mathrm{2}}{\mathrm{\Gamma }}^{\mathrm{2}}\mathit{\kappa }\mathrm{+}{\mathrm{1}}^{\mathrm{)}}\\ \mathit{N}{\mathit{c}}_{\mathrm{21}}& \mathrm{=}& \mathrm{0}\\ \mathit{N}{\mathit{c}}_{\mathrm{22}}& \mathrm{=}& {\mathit{v}}_{\mathit{z}}\mathit{N}\\ \mathit{N}{\mathit{c}}_{\mathrm{23}}& \mathrm{=}& \mathrm{0}\\ \mathit{N}{\mathit{c}}_{\mathrm{24}}& \mathrm{=}& {\mathit{v}}_{\mathit{z}}\mathit{\rho \chi }\\ \mathit{N}{\mathit{c}}_{\mathrm{25}}& \mathrm{=}& \frac{\mathrm{-}{\mathit{v}}_{\mathit{x}}{\mathit{v}}_{\mathit{z}}\mathrm{(}\mathrm{-}\mathit{h}\mathrm{+}\mathit{\rho \chi }\mathrm{+}\mathit{\rho h\kappa }\mathrm{)}}{\mathit{\rho h}}\\ \mathit{N}{\mathit{c}}_{\mathrm{31}}& \mathrm{=}& \mathrm{0}\\ \mathit{N}{\mathit{c}}_{\mathrm{32}}& \mathrm{=}& \mathrm{0}\\ \mathit{N}{\mathit{c}}_{\mathrm{33}}& \mathrm{=}& {\mathit{v}}_{\mathit{z}}\mathit{N}\\ \mathit{N}{\mathit{c}}_{\mathrm{34}}& \mathrm{=}& {\mathit{v}}_{\mathit{y}}\mathit{\rho \chi }\\ \mathit{N}{\mathit{c}}_{\mathrm{35}}& \mathrm{=}& \frac{\mathrm{-}{\mathit{v}}_{\mathit{y}}{\mathit{v}}_{\mathit{z}}\mathrm{(}\mathrm{-}\mathit{h}\mathrm{+}\mathit{\rho \chi }\mathrm{+}\mathit{\rho h\kappa }\mathrm{)}}{\mathit{\rho h}}\\ \mathit{N}{\mathit{c}}_{\mathrm{41}}& \mathrm{=}& \mathrm{0}\\ \mathit{N}{\mathit{c}}_{\mathrm{42}}& \mathrm{=}& \mathrm{0}\\ \mathit{N}{\mathit{c}}_{\mathrm{43}}& \mathrm{=}& \mathrm{0}\\ \mathit{N}{\mathit{c}}_{\mathrm{44}}& \mathrm{=}& {\mathit{v}}_{\mathit{z}}{\mathrm{\Gamma }}^{\mathrm{2}}\mathrm{(}\mathit{\rho \chi }\mathrm{-}\mathit{h}\mathrm{+}\mathit{h\rho \kappa }\mathrm{)}\\ \mathit{N}{\mathit{c}}_{\mathrm{45}}& \mathrm{=}& \frac{\mathrm{1}}{\mathit{\rho h}{\mathrm{\Gamma }}^{\mathrm{2}}}\mathrm{(}\mathit{\rho h}{\mathrm{\Gamma }}^{\mathrm{2}}\mathit{\kappa }\mathrm{(}\mathrm{1}\mathrm{-}{{\mathit{v}}_{\mathit{z}}^{\mathrm{2}}}^{\mathrm{)}}\mathrm{-}\mathit{h}\mathrm{-}\mathit{h}{\mathrm{\Gamma }}^{\mathrm{2}}\mathrm{(}{\mathit{v}}_{\mathit{x}}^{\mathrm{2}}\mathrm{+}{{\mathit{v}}_{\mathit{y}}^{\mathrm{2}}}^{\mathrm{)}}\mathrm{+}\mathit{\rho \chi }{\mathrm{\Gamma }}^{\mathrm{2}}{\mathrm{(}{\mathit{v}}_{\mathit{x}}^{\mathrm{2}}\mathrm{+}{{\mathit{v}}_{\mathit{y}}^{\mathrm{2}}}^{\mathrm{)}}}^{\mathrm{)}}\\ \mathit{N}{\mathit{c}}_{\mathrm{51}}& \mathrm{=}& \mathrm{0}\\ \mathit{N}{\mathit{c}}_{\mathrm{52}}& \mathrm{=}& \mathrm{0}\\ \mathit{N}{\mathit{c}}_{\mathrm{53}}& \mathrm{=}& \mathrm{0}\\ \mathit{N}{\mathit{c}}_{\mathrm{54}}& \mathrm{=}& \mathrm{-}{\mathit{\rho }}^{\mathrm{2}}\mathit{h}{\mathrm{\Gamma }}^{\mathrm{2}}\mathit{\chi }\\ \mathit{N}{\mathit{c}}_{\mathrm{55}}& \mathrm{=}& {\mathit{v}}_{\mathit{z}}{\mathrm{\Gamma }}^{\mathrm{2}}\mathrm{(}\mathit{\rho \chi }\mathrm{-}\mathit{h}\mathrm{+}\mathit{h\rho \kappa }\mathrm{)}\mathit{.}\end{array}$

Appendix A.4: Eigenstructure

Fig. B.1 Density, AMR levels (red dashed line), parallel and transverse velocity and pressure in the frame of the laboratory for the shock test with Γmax = 120. Left panels: uniform grid nx = 512; middle panels: uniform grid nx = 217 = 131072; right panels: nx = 64, lmax = 17. The black solid line represents the MUSCL method, while the black dashed line represents the PLM method. At high resolution, they give the same result, which is very close to the analytic solution.

The eigenvalues of A_x are (see e.g. Falle & Komissarov 1996) $\begin{array}{ccc}{\mathit{\lambda }}_{\mathrm{-}}& \mathrm{=}& {\mathit{v}}_{\mathit{x}}\frac{\mathrm{(}\mathrm{1}\mathrm{-}{{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}}^{\mathrm{)}}\mathrm{-}{\mathit{c}}_{\mathrm{s}}{\mathrm{\Gamma }}^{-1}\mathit{\omega }}{\mathrm{1}\mathrm{-}{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}{\mathit{v}}^{\mathrm{2}}}\\ {\mathit{\lambda }}_{\mathrm{0}\mathit{x}}& \mathrm{=}& {\mathit{v}}_{\mathit{x}}\\ {\mathit{\lambda }}_{\mathrm{0}\mathit{y}}& \mathrm{=}& {\mathit{v}}_{\mathit{x}}\\ {\mathit{\lambda }}_{\mathrm{0}\mathit{z}}& \mathrm{=}& {\mathit{v}}_{\mathit{x}}\\ {\mathit{\lambda }}_{\mathrm{+}}& \mathrm{=}& {\mathit{v}}_{\mathit{x}}\frac{\mathrm{(}\mathrm{1}\mathrm{-}{{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}}^{\mathrm{)}}\mathrm{+}{\mathit{c}}_{\mathrm{s}}{\mathrm{\Gamma }}^{-1}\mathit{\omega }}{\mathrm{1}\mathrm{-}{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}{\mathit{v}}^{\mathrm{2}}}\mathit{,}\end{array}$with $\mathit{\omega }\mathrm{=}\sqrt{\mathrm{1}\mathrm{-}{\mathit{v}}_{\mathit{x}}^{\mathrm{2}}\mathrm{-}{\mathit{c}}_{\mathrm{s}}\mathrm{(}{\mathit{v}}_{\mathit{y}}^{\mathrm{2}}\mathrm{+}{{\mathit{v}}_{\mathit{z}}^{\mathrm{2}}}^{\mathrm{)}}}\mathit{.}$(A.23)The left eigenvectors are given by $\begin{array}{ccc}{\mathit{L}}_{\mathrm{-}}& \mathrm{=}& \left(\mathrm{0}\mathit{,}\mathrm{-}\frac{\mathit{\rho }\mathrm{\Gamma }}{\mathrm{2}{\mathit{c}}_{\mathrm{s}}\mathit{\omega }}\mathit{,}\mathrm{0}\mathit{,}\mathrm{0}\mathit{,}\frac{\mathrm{1}}{\mathrm{2}{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}\mathit{h}}\right)\\ {\mathit{L}}_{\mathrm{0}\mathit{x}}& \mathrm{=}& \left(\mathrm{1}\mathit{,}\mathrm{0}\mathit{,}\mathrm{0}\mathit{,}\mathrm{0}\mathit{,}\mathrm{-}\frac{\mathrm{1}}{{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}\mathit{h}}\right)\\ {\mathit{L}}_{\mathrm{0}\mathit{y}}& \mathrm{=}& \left(\mathrm{0}\mathit{,}\frac{{\mathit{v}}_{\mathit{x}}{\mathit{v}}_{\mathit{y}}}{\mathrm{1}\mathrm{-}{\mathit{v}}_{\mathit{x}}^{\mathrm{2}}}\mathit{,}\mathrm{1}\mathit{,}\mathrm{0}\mathit{,}\frac{{\mathit{v}}_{\mathit{y}}}{\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{v}}_{\mathit{x}}^{\mathrm{2}}\mathrm{)}{\mathrm{\Gamma }}^{\mathrm{2}}\mathit{\rho h}}\right)\\ {\mathit{L}}_{\mathrm{0}\mathit{z}}& \mathrm{=}& \left(\mathrm{0}\mathit{,}\frac{{\mathit{v}}_{\mathit{x}}{\mathit{v}}_{\mathit{z}}}{\mathrm{1}\mathrm{-}{\mathit{v}}_{\mathit{x}}^{\mathrm{2}}}\mathit{,}\mathrm{0}\mathit{,}\mathrm{1}\mathit{,}\frac{{\mathit{v}}_{\mathit{z}}}{\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{v}}_{\mathit{x}}^{\mathrm{2}}\mathrm{)}{\mathrm{\Gamma }}^{\mathrm{2}}\mathit{\rho h}}\right)\\ {\mathit{L}}_{\mathrm{+}}& \mathrm{=}& \left(\mathrm{0}\mathit{,}\frac{\mathit{\rho }\mathrm{\Gamma }\mathit{\omega }}{\mathrm{2}{\mathit{c}}_{\mathrm{s}}}\mathit{,}\mathrm{0}\mathit{,}\mathrm{0}\mathit{,}\frac{\mathrm{1}}{\mathrm{2}{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}\mathit{h}}\right)\mathrm{·}\end{array}$The right eigenvectors are given by $\begin{array}{ccc}{\mathit{R}}_{\mathrm{-}}& \mathrm{=}& {\left(\mathrm{1}\mathit{,}\mathrm{-}\frac{{\mathit{c}}_{\mathrm{s}}\mathit{\omega }}{\mathit{\rho }\mathrm{\Gamma }}\mathit{,}\frac{\mathrm{-}{\mathit{v}}_{\mathit{y}}{\mathit{c}}_{\mathrm{s}}\mathrm{(}\mathrm{\Gamma }\mathit{\omega }{\mathit{v}}_{\mathit{x}}\mathrm{+}{\mathit{c}}_{\mathrm{s}}\mathrm{)}}{\mathit{\rho }{\mathrm{\Gamma }}^{\mathrm{2}}\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{v}}_{\mathit{x}}^{\mathrm{2}}\mathrm{)}}\mathit{,}\frac{\mathrm{-}{\mathit{v}}_{\mathit{z}}{\mathit{c}}_{\mathrm{s}}\mathrm{(}\mathrm{\Gamma }\mathit{\omega }{\mathit{v}}_{\mathit{x}}\mathrm{+}{\mathit{c}}_{\mathrm{s}}\mathrm{)}}{\mathit{\rho }{\mathrm{\Gamma }}^{\mathrm{2}}\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{v}}_{\mathit{x}}^{\mathrm{2}}\mathrm{)}}\mathit{,}{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}\mathit{h}\right)}^{\mathrm{T}}\\ {\mathit{R}}_{\mathrm{0}\mathit{x}}& \mathrm{=}& \mathrm{(}\mathrm{1}\mathit{,}\mathrm{0}\mathit{,}\mathrm{0}\mathit{,}\mathrm{0}\mathit{,}\mathrm{0}{\mathrm{)}}^{\mathrm{T}}\\ {\mathit{R}}_{\mathrm{0}\mathit{y}}& \mathrm{=}& \mathrm{(}\mathrm{0}\mathit{,}\mathrm{0}\mathit{,}\mathrm{1}\mathit{,}\mathrm{0}\mathit{,}\mathrm{0}{\mathrm{)}}^{\mathrm{T}}\\ {\mathit{R}}_{\mathrm{0}\mathit{z}}& \mathrm{=}& \mathrm{(}\mathrm{0}\mathit{,}\mathrm{0}\mathit{,}\mathrm{0}\mathit{,}\mathrm{1}\mathit{,}\mathrm{0}{\mathrm{)}}^{\mathrm{T}}\\ {\mathit{R}}_{\mathrm{+}}& \mathrm{=}& {\left(\mathrm{1}\mathit{,}\frac{{\mathit{c}}_{\mathrm{s}}\mathit{\omega }}{\mathit{\rho }\mathrm{\Gamma }}\mathit{,}\frac{{\mathit{v}}_{\mathit{y}}{\mathit{c}}_{\mathrm{s}}\mathrm{(}\mathrm{\Gamma }\mathit{\omega }{\mathit{v}}_{\mathit{x}}\mathrm{-}{\mathit{c}}_{\mathrm{s}}\mathrm{)}}{\mathit{\rho }{\mathrm{\Gamma }}^{\mathrm{2}}\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{v}}_{\mathit{x}}^{\mathrm{2}}\mathrm{)}}\mathit{,}\frac{{\mathit{v}}_{\mathit{z}}{\mathit{c}}_{\mathrm{s}}\mathrm{(}\mathrm{\Gamma }\mathit{\omega }{\mathit{v}}_{\mathit{x}}\mathrm{-}{\mathit{c}}_{\mathrm{s}}\mathrm{)}}{\mathit{\rho }{\mathrm{\Gamma }}^{\mathrm{2}}\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{v}}_{\mathit{x}}^{\mathrm{2}}\mathrm{)}}\mathit{,}{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}\mathit{h}\right)}^{\mathrm{T}}\mathrm{·}\end{array}$Similarly for A_y one has $\begin{array}{ccc}{\mathit{\lambda }}_{\mathrm{-}}& \mathrm{=}& {\mathit{v}}_{\mathit{y}}\frac{\mathrm{(}\mathrm{1}\mathrm{-}{{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}}^{\mathrm{)}}\mathrm{-}{\mathit{c}}_{\mathrm{s}}{\mathrm{\Gamma }}^{-1}\mathit{\omega }}{\mathrm{1}\mathrm{-}{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}{\mathit{v}}^{\mathrm{2}}}\\ {\mathit{\lambda }}_{\mathrm{0}\mathit{x}}& \mathrm{=}& {\mathit{v}}_{\mathit{y}}\\ {\mathit{\lambda }}_{\mathrm{0}\mathit{y}}& \mathrm{=}& {\mathit{v}}_{\mathit{y}}\\ {\mathit{\lambda }}_{\mathrm{0}\mathit{z}}& \mathrm{=}& {\mathit{v}}_{\mathit{y}}\\ {\mathit{\lambda }}_{\mathrm{+}}& \mathrm{=}& {\mathit{v}}_{\mathit{y}}\frac{\mathrm{(}\mathrm{1}\mathrm{-}{{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}}^{\mathrm{)}}\mathrm{+}{\mathit{c}}_{\mathrm{s}}{\mathrm{\Gamma }}^{-1}\mathit{\omega }}{\mathrm{1}\mathrm{-}{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}{\mathit{v}}^{\mathrm{2}}}\mathrm{·}\end{array}$with $\mathit{\omega }\mathrm{=}\sqrt{\mathrm{1}\mathrm{-}{\mathit{v}}_{\mathit{y}}^{\mathrm{2}}\mathrm{-}{\mathit{c}}_{\mathrm{s}}\mathrm{(}{\mathit{v}}_{\mathit{x}}^{\mathrm{2}}\mathrm{+}{{\mathit{v}}_{\mathit{z}}^{\mathrm{2}}}^{\mathrm{)}}}\mathit{.}$(A.39)The left eigenvectors are given by $\begin{array}{ccc}{\mathit{L}}_{\mathrm{-}}& \mathrm{=}& \left(\mathrm{0}\mathit{,}\mathrm{0}\mathit{,}\mathrm{-}\frac{\mathit{\rho }\mathrm{\Gamma }}{\mathrm{2}{\mathit{c}}_{\mathrm{s}}\mathit{\omega }}\mathit{,}\mathrm{0}\mathit{,}\frac{\mathrm{1}}{\mathrm{2}{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}\mathit{h}}\right)\\ {\mathit{L}}_{\mathrm{0}\mathit{x}}& \mathrm{=}& \left(\mathrm{0}\mathit{,}\mathrm{1}\mathit{,}\frac{{\mathit{v}}_{\mathit{y}}{\mathit{v}}_{\mathit{x}}}{\mathrm{1}\mathrm{-}{\mathit{v}}_{\mathit{y}}^{\mathrm{2}}}\mathit{,}\mathrm{0}\mathit{,}\frac{{\mathit{v}}_{\mathit{x}}}{\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{v}}_{\mathit{y}}^{\mathrm{2}}\mathrm{)}{\mathrm{\Gamma }}^{\mathrm{2}}\mathit{\rho h}}\right)\end{array}$$\begin{array}{ccc}& & \\ {\mathit{L}}_{\mathrm{0}\mathit{y}}& \mathrm{=}& \left(\mathrm{1}\mathit{,}\mathrm{0}\mathit{,}\mathrm{0}\mathit{,}\mathrm{0}\mathit{,}\mathrm{-}\frac{\mathrm{1}}{{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}\mathit{h}}\right)\\ {\mathit{L}}_{\mathrm{0}\mathit{z}}& \mathrm{=}& \left(\mathrm{0}\mathit{,}\mathrm{0}\mathit{,}\frac{{\mathit{v}}_{\mathit{y}}{\mathit{v}}_{\mathit{z}}}{\mathrm{1}\mathrm{-}{\mathit{v}}_{\mathit{y}}^{\mathrm{2}}}\mathit{,}\mathrm{1}\mathit{,}\frac{{\mathit{v}}_{\mathit{z}}}{\mathrm{(}\mathrm{1}\mathrm{-}{{\mathit{v}}_{\mathit{y}}^{\mathrm{2}}}^{\mathrm{)}}{\mathrm{\Gamma }}^{\mathrm{2}}\mathit{\rho h}}\right)\\ {\mathit{L}}_{\mathrm{+}}& \mathrm{=}& \left(\mathrm{0}\mathit{,}\mathrm{0}\mathit{,}\frac{\mathit{\rho }\mathrm{\Gamma }}{\mathrm{2}{\mathit{c}}_{\mathrm{s}}\mathit{\omega }}\mathit{,}\mathrm{0}\mathit{,}\frac{\mathrm{1}}{\mathrm{2}{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}\mathit{h}}\right)\mathrm{·}\end{array}$The right eigenvectors are given by $\begin{array}{ccc}{\mathit{R}}_{\mathrm{-}}& \mathrm{=}& {\left(\mathrm{1}\mathit{,}\frac{\mathrm{-}{\mathit{v}}_{\mathit{x}}{\mathit{c}}_{\mathrm{s}}\mathrm{(}\mathrm{\Gamma }\mathit{\omega }{\mathit{v}}_{\mathit{y}}\mathrm{+}{\mathit{c}}_{\mathrm{s}}\mathrm{)}}{\mathit{\rho }{\mathrm{\Gamma }}^{\mathrm{2}}\mathrm{(}\mathrm{1}\mathrm{-}{{\mathit{v}}_{\mathit{y}}^{\mathrm{2}}}^{\mathrm{)}}}\mathit{,}\mathrm{-}\frac{{\mathit{c}}_{\mathrm{s}}\mathit{\omega }}{\mathit{\rho }\mathrm{\Gamma }}\mathit{,}\frac{\mathrm{-}{\mathit{v}}_{\mathit{z}}{\mathit{c}}_{\mathrm{s}}\mathrm{(}\mathrm{\Gamma }\mathit{\omega }{\mathit{v}}_{\mathit{y}}\mathrm{+}{\mathit{c}}_{\mathrm{s}}\mathrm{)}}{\mathit{\rho }{\mathrm{\Gamma }}^{\mathrm{2}}\mathrm{(}\mathrm{1}\mathrm{-}{{\mathit{v}}_{\mathit{y}}^{\mathrm{2}}}^{\mathrm{)}}}\mathit{,}{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}\mathit{h}\right)}^{\mathrm{T}}\\ {\mathit{R}}_{\mathrm{0}\mathit{x}}& \mathrm{=}& {\mathrm{(}\mathrm{0}\mathit{,}\mathrm{1}\mathit{,}\mathrm{0}\mathit{,}\mathrm{0}\mathit{,}{\mathrm{0}}^{\mathrm{)}}}^{\mathrm{T}}\\ {\mathit{R}}_{\mathrm{0}\mathit{y}}& \mathrm{=}& {\mathrm{(}\mathrm{1}\mathit{,}\mathrm{0}\mathit{,}\mathrm{0}\mathit{,}\mathrm{0}\mathit{,}{\mathrm{0}}^{\mathrm{)}}}^{\mathrm{T}}\\ {\mathit{R}}_{\mathrm{0}\mathit{z}}& \mathrm{=}& {\mathrm{(}\mathrm{0}\mathit{,}\mathrm{0}\mathit{,}\mathrm{0}\mathit{,}\mathrm{1}\mathit{,}{\mathrm{0}}^{\mathrm{)}}}^{\mathrm{T}}\\ {\mathit{R}}_{\mathrm{+}}& \mathrm{=}& {\left(\mathrm{1}\mathit{,}\frac{{\mathit{v}}_{\mathit{x}}{\mathit{c}}_{\mathrm{s}}\mathrm{(}\mathrm{\Gamma }\mathit{\omega }{\mathit{v}}_{\mathit{y}}\mathrm{-}{\mathit{c}}_{\mathrm{s}}\mathrm{)}}{\mathit{\rho }{\mathrm{\Gamma }}^{\mathrm{2}}\mathrm{(}\mathrm{1}\mathrm{-}{{\mathit{v}}_{\mathit{y}}^{\mathrm{2}}}^{\mathrm{)}}}\mathit{,}\mathrm{+}\frac{{\mathit{c}}_{\mathrm{s}}\mathit{\omega }}{\mathit{\rho }\mathrm{\Gamma }}\mathit{,}\frac{{\mathit{v}}_{\mathit{z}}{\mathit{c}}_{\mathrm{s}}\mathrm{(}\mathrm{\Gamma }\mathit{\omega }{\mathit{v}}_{\mathit{y}}\mathrm{-}{\mathit{c}}_{\mathrm{s}}\mathrm{)}}{\mathit{\rho }{\mathrm{\Gamma }}^{\mathrm{2}}\mathrm{(}\mathrm{1}\mathrm{-}{{\mathit{v}}_{\mathit{y}}^{\mathrm{2}}}^{\mathrm{)}}}\mathit{,}{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}\mathit{h}\right)}^{\mathrm{T}}\mathrm{·}\end{array}$For A_z one has $\begin{array}{ccc}{\mathit{\lambda }}_{\mathrm{-}}& \mathrm{=}& {\mathit{v}}_{\mathit{z}}\frac{\mathrm{(}\mathrm{1}\mathrm{-}{{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}}^{\mathrm{)}}\mathrm{-}{\mathit{c}}_{\mathrm{s}}{\mathrm{\Gamma }}^{-1}\mathit{\omega }}{\mathrm{1}\mathrm{-}{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}{\mathit{v}}^{\mathrm{2}}}\\ {\mathit{\lambda }}_{\mathrm{0}\mathit{x}}& \mathrm{=}& {\mathit{v}}_{\mathit{z}}\\ {\mathit{\lambda }}_{\mathrm{0}\mathit{y}}& \mathrm{=}& {\mathit{v}}_{\mathit{z}}\\ {\mathit{\lambda }}_{\mathrm{0}\mathit{z}}& \mathrm{=}& {\mathit{v}}_{\mathit{z}}\\ {\mathit{\lambda }}_{\mathrm{+}}& \mathrm{=}& {\mathit{v}}_{\mathit{z}}\frac{\mathrm{(}\mathrm{1}\mathrm{-}{{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}}^{\mathrm{)}}\mathrm{+}{\mathit{c}}_{\mathrm{s}}{\mathrm{\Gamma }}^{-1}\mathit{\omega }}{\mathrm{1}\mathrm{-}{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}{\mathit{v}}^{\mathrm{2}}}\mathit{,}\end{array}$with $\mathit{\omega }\mathrm{=}\sqrt{\mathrm{1}\mathrm{-}{\mathit{v}}_{\mathit{z}}^{\mathrm{2}}\mathrm{-}{\mathit{c}}_{\mathrm{s}}\mathrm{(}{\mathit{v}}_{\mathit{x}}^{\mathrm{2}}\mathrm{+}{{\mathit{v}}_{\mathit{y}}^{\mathrm{2}}}^{\mathrm{)}}}\mathit{.}$(A.55)The left eigenvectors are given by $\begin{array}{ccc}{\mathit{L}}_{\mathrm{-}}& \mathrm{=}& \left(\mathrm{0}\mathit{,}\mathrm{0}\mathit{,}\mathrm{0}\mathit{,}\mathrm{-}\frac{\mathit{\rho }\mathrm{\Gamma }}{\mathrm{2}{\mathit{c}}_{\mathrm{s}}\mathit{\omega }}\mathit{,}\frac{\mathrm{1}}{\mathrm{2}{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}\mathit{h}}\right)\\ {\mathit{L}}_{\mathrm{0}\mathit{x}}& \mathrm{=}& \left(\mathrm{0}\mathit{,}\mathrm{1}\mathit{,}\mathrm{0}\mathit{,}\frac{{\mathit{v}}_{\mathit{x}}{\mathit{v}}_{\mathit{z}}}{\mathrm{1}\mathrm{-}{\mathit{v}}_{\mathit{z}}^{\mathrm{2}}}\mathit{,}\frac{{\mathit{v}}_{\mathit{x}}}{\mathrm{(}\mathrm{1}\mathrm{-}{{\mathit{v}}_{\mathit{z}}^{\mathrm{2}}}^{\mathrm{)}}{\mathrm{\Gamma }}^{\mathrm{2}}\mathrm{)}\mathit{\rho h}}\right)\\ {\mathit{L}}_{\mathrm{0}\mathit{y}}& \mathrm{=}& \left(\mathrm{0}\mathit{,}\mathrm{0}\mathit{,}\frac{{\mathit{v}}_{\mathit{y}}}{\mathrm{(}\mathrm{1}\mathrm{-}{{\mathit{v}}_{\mathit{z}}^{\mathrm{2}}}^{\mathrm{)}}{\mathrm{\Gamma }}^{\mathrm{2}}\mathit{\rho h}}\mathit{,}\frac{{\mathit{v}}_{\mathit{y}}{\mathit{v}}_{\mathit{z}}}{\mathrm{1}\mathrm{-}{\mathit{v}}_{\mathit{z}}^{\mathrm{2}}}\mathit{,}\mathrm{0}\mathit{,}\right)\\ {\mathit{L}}_{\mathrm{0}\mathit{z}}& \mathrm{=}& \left(\mathrm{1}\mathit{,}\mathrm{0}\mathit{,}\mathrm{0}\mathit{,}\mathrm{0}\mathit{,}\mathrm{-}\frac{\mathrm{1}}{{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}\mathit{h}}\right)\\ {\mathit{L}}_{\mathrm{+}}& \mathrm{=}& \left(\mathrm{0}\mathit{,}\mathrm{0}\mathit{,}\mathrm{0}\mathit{,}\frac{\mathit{\rho }\mathrm{\Gamma }}{\mathrm{2}{\mathit{c}}_{\mathrm{s}}\mathit{\omega }}\mathit{,}\frac{\mathrm{1}}{\mathrm{2}{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}\mathit{h}}\right)\mathrm{·}\end{array}$The right eigenvectors are given by $\begin{array}{ccc}{\mathit{R}}_{\mathrm{-}}& \mathrm{=}& {\left(\mathrm{1}\mathit{,}\frac{\mathrm{-}{\mathit{v}}_{\mathit{x}}{\mathit{c}}_{\mathrm{s}}\mathrm{(}\mathrm{\Gamma }\mathit{\omega }{\mathit{v}}_{\mathit{z}}\mathrm{+}{\mathit{c}}_{\mathrm{s}}\mathrm{)}}{\mathit{\rho }{\mathrm{\Gamma }}^{\mathrm{2}}\mathrm{(}\mathrm{1}\mathrm{-}{{\mathit{v}}_{\mathit{z}}^{\mathrm{2}}}^{\mathrm{)}}}\mathit{,}\frac{\mathrm{-}{\mathit{v}}_{\mathit{y}}{\mathit{c}}_{\mathrm{s}}\mathrm{(}\mathrm{\Gamma }\mathit{\omega }{\mathit{v}}_{\mathit{z}}\mathrm{+}{\mathit{c}}_{\mathrm{s}}\mathrm{)}}{\mathit{\rho }{\mathrm{\Gamma }}^{\mathrm{2}}\mathrm{(}\mathrm{1}\mathrm{-}{{\mathit{v}}_{\mathit{z}}^{\mathrm{2}}}^{\mathrm{)}}}\mathit{,}\mathrm{-}\frac{{\mathit{c}}_{\mathrm{s}}\mathit{\omega }}{\mathit{\rho }\mathrm{\Gamma }}\mathit{,}{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}\mathit{h}\right)}^{\mathrm{T}}\\ {\mathit{R}}_{\mathrm{0}\mathit{x}}& \mathrm{=}& {\mathrm{(}\mathrm{0}\mathit{,}\mathrm{1}\mathit{,}\mathrm{0}\mathit{,}\mathrm{0}\mathit{,}{\mathrm{0}}^{\mathrm{)}}}^{\mathrm{T}}\\ {\mathit{R}}_{\mathrm{0}\mathit{y}}& \mathrm{=}& {\mathrm{(}\mathrm{0}\mathit{,}\mathrm{0}\mathit{,}\mathrm{1}\mathit{,}\mathrm{0}\mathit{,}{\mathrm{0}}^{\mathrm{)}}}^{\mathrm{T}}\\ {\mathit{R}}_{\mathrm{0}\mathit{z}}& \mathrm{=}& {\mathrm{(}\mathrm{1}\mathit{,}\mathrm{0}\mathit{,}\mathrm{0}\mathit{,}\mathrm{0}\mathit{,}{\mathrm{0}}^{\mathrm{)}}}^{\mathrm{T}}\\ {\mathit{R}}_{\mathrm{+}}& \mathrm{=}& {\left(\mathrm{1}\mathit{,}\frac{{\mathit{v}}_{\mathit{x}}{\mathit{c}}_{\mathrm{s}}\mathrm{(}\mathrm{\Gamma }\mathit{\omega }{\mathit{v}}_{\mathit{z}}\mathrm{-}{\mathit{c}}_{\mathrm{s}}\mathrm{)}}{\mathit{\rho }{\mathrm{\Gamma }}^{\mathrm{2}}\mathrm{(}\mathrm{1}\mathrm{-}{{\mathit{v}}_{\mathit{y}}^{\mathrm{2}}}^{\mathrm{)}}}\mathit{,}\frac{{\mathit{v}}_{\mathit{x}}{\mathit{c}}_{\mathrm{s}}\mathrm{(}\mathrm{\Gamma }\mathit{\omega }{\mathit{v}}_{\mathit{z}}\mathrm{-}{\mathit{c}}_{\mathrm{s}}\mathrm{)}}{\mathit{\rho }{\mathrm{\Gamma }}^{\mathrm{2}}\mathrm{(}\mathrm{1}\mathrm{-}{{\mathit{v}}_{\mathit{y}}^{\mathrm{2}}}^{\mathrm{)}}}\mathit{,}{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}\mathit{h}\right)}^{\mathrm{T}}\mathit{.}\end{array}$