Issue 
A&A
Volume 661, May 2022



Article Number  A117  
Number of page(s)  28  
Section  Stellar structure and evolution  
DOI  https://doi.org/10.1051/00046361/202142968  
Published online  13 May 2022 
3D RMHD simulations of jetwind interactions in highmass Xray binaries
^{1}
Departament d’Astronomia i Astrofísica, Universitat de València, Dr. Moliner 50, 46100 Burjassot, València, Spain
email: jose.lopezmiralles@uv.es
^{2}
Observatori Astronòmic, Universitat de València, C/ Catedràtic José Beltrán 2, 46980 Paterna, València, Spain
^{3}
Aurora Technology for the European Space Agency, ESAC/ESA, Camino Bajo del Castillo s/n, Urb. Villafranca del Castillo, 28691 Villanueva de la Cañada, Madrid, Spain
^{4}
Institut de Ciències del Cosmos (ICC), Universitat de Barcelona (IEECUB), Martí i Franquès 1, 08028 Barcelona, Spain
Received:
21
December
2021
Accepted:
21
February
2022
Context. Relativistic jets are ubiquitous in the Universe. In microquasars, especially in highmass Xray binaries, the interaction of jets with the strong winds driven by the massive and hot companion star in the vicinity of the compact object is fundamental for understanding the jet dynamics, nonthermal emission, and longterm stability. However, the role of the jet magnetic field in this process is unclear. In particular, it is still debated whether the magnetic field favors jet collimation or triggers more instabilities that can jeopardize the jet evolution outside the binary.
Aims. We study the dynamical role of weak and moderate to strong toroidal magnetic fields during the first several hundred seconds of jet propagation through the stellar wind, focusing on the magnetized flow dynamics and the mechanisms of energy conversion.
Methods. We developed the code Lóstrego v1.0, a new 3D relativistic magnetohydrodynamics code to simulate astrophysical plasmas in Cartesian coordinates. Using this tool, we performed the first 3D relativistic magnetohydrodynamics numerical simulations of relativistic magnetized jets propagating through the clumpy stellar wind in a highmass Xray binary. To highlight the effect of the magnetic field in the jet dynamics, we compared the results of our analysis with those of previous hydrodynamical simulations.
Results. The overall morphology and dynamics of weakly magnetized jet models is similar to previous hydrodynamical simulations, where the jet head generates a strong shock in the ambient medium and the initial overpressure with respect to the stellar wind drives one or more recollimation shocks. On the timescales of our simulations (i.e., t < 200 s), these jets are ballistic and seem to be more stable against internal instabilities than jets with the same power in the absence of fields. However, moderate to strong toroidal magnetic fields favor the development of currentdriven instabilities and the disruption of the jet within the binary. A detailed analysis of the energy distribution in the relativistic outflow and the ambient medium reveals that magnetic and internal energies can both contribute to the effective acceleration of the jet. Moreover, we verified that the jet feedback into the ambient medium is highly dependent on the jet energy distribution at injection, where hotter, more diluted and/or more magnetized jets are more efficient. This was anticipated by feedback studies in the case of jets in active galaxies.
Key words: magnetohydrodynamics (MHD) / Xrays: binaries / ISM: jets and outflows / stars: winds, outflows / relativistic processes
© ESO 2022
1. Introduction
Jets are the most spectacular and powerful consequences of accretion onto compact objects; they have been observed in systems containing white dwarfs, neutron stars (NSs), and black holes (BHs) of all mass scales, from stellar mass in Xray binaries (XRBs) to supermassive in active galactic nuclei (AGN). These outflows are relativistic and extract a large and possibly dominant fraction of the total available accretion energy (Ghisellini et al. 2014), and they may even tap the power of the BH spin (Blandford & Znajek 1977). Studies of plasma dynamics and radiative processes in relativistic jets, and especially the comparison of different classes of compact objects (Migliari & Fender 2006), allows us to investigate the physics of stronggravity and curved spacetime; the presence of a stellar surface in NSs, or the existence of an event horizon in BHs; the role of magnetic fields in jet formation, collimation, and evolution; and the shockacceleration mechanisms of plasma particles, to name some representative examples (see, e.g., Event Horizon Telescope Collaboration 2019; Kim et al. 2020; Janssen et al. 2021). Moreover, powerful relativistic jets are one of the main ways in which accreting black holes provide kinetic feedback to their surroundings (see, e.g., Bordas et al. 2009; McNamara & Nulsen 2007, in the case of microquasars and AGN, respectively). However, most of the fundamental questions regarding jet formation mechanisms, composition, acceleration, collimation, and interaction with the interstellar medium (ISM) are still debated.
XRBs, also called microquasars, are binary systems hosting a compact object (i.e., a stellar mass BH or NS) and a companion star that has not collapsed and supplies matter to the compact object through Rochelobe overflow (lowmass XRBs) or through the capture of stellar winds (highmass XRBs; HMXBs). This matter accumulates around the central object in the form of an accretion disk. Most XRBs are transient objects (Belloni & Motta 2016), but in some spectral states, these systems produce powerful bipolar relativistic jets launched by magnetocentrifugal forces (Blandford & Znajek 1977; Blandford & Payne 1982) that emit nonthermal synchroton radiation (Mirabel & Rodríguez 1999). On smaller scales, these processes mimic most of the phenomena observed in quasars and AGNs, but in XRBs, accretion varies on much faster humanly accessible timescales than in supermassive BHs. This allows us to observe and follow the evolution of the system and to investigate the link between jet formation and the different accretion states. Therefore, numerical simulations that study the physical conditions that may reproduce existing observations are strongly needed.
In order to understand how XRBs produce relativistic ejections and how these jets affect the binary system and the surrounding ISM, we may distinguish three different regions of interest: (1) the innermost scale close to the compact object, where the jet is launched by magnetocentrifugal forces from the inner accretion disk (Migliari et al. 2003, 2004; Fender et al. 2004; Kylafis et al. 2012; Marino et al. 2020), (2) the binary scale, where the collimated outflow interacts with the wind of the companion star and may presumably be disrupted (Perucho & BoschRamon 2008, 2012; Perucho et al. 2010a; BoschRamon & Barkov 2016; de la Cita et al. 2017; Barkov & BoschRamon 2022), and (3) the interaction of the jet with the ISM, where nonthermal particles could produce a significant amount of radiation at different wavelengths (Bordas et al. 2009; BoschRamon et al. 2011; Yoon et al. 2011; MonceauBaroux et al. 2014).
For HMXBs, jetwind interactions are particularly relevant because of the strong winds that the massive and hot companion star drives in the vicinity of the compact object. As mentioned above, this wind can have a strong impact on jet dynamics by jeopardizing its stability and eventually preventing a detection of the jet (Perucho & BoschRamon 2008; Perucho et al. 2010a). Moreover, the onesided impact of the wind on the relativistic flow may produce strong collisionless shocks that lead to efficient particle acceleration (Rieger et al. 2007), resulting in significant nonthermal radiation from synchroton, inverse Compton (IC), or even protonproton collisions (see, e.g., Romero et al. 2003; Molina et al. 2019). These processes may be at the origin of the gammaray emission detected in some XRBs such as LS 5039 (e.g., Paredes et al. 2000; Aharonian et al. 2005), LS I+61303 (e.g., Tavani et al. 1998; Albert et al. 2006) (if they really host accreting BHs), Cygnus X1 (e.g., Albert et al. 2007; Zanin et al. 2016), and Cygnus X3 (e.g., Albert et al. 2007; Zanin et al. 2016).
Several attempts have been made in the past to describe jetwind interactions from the point of view of theory and numerical simulations. Perucho & BoschRamon (2008) performed numerical twodimensional (2D) simulations of jets crossing the stellar wind, reporting the formation of strong recollimation shocks that can accelerate particles efficiently and produce nonthermal radiation. Similar results, but with threedimensional (3D) hydrodynamical simulations, were found by Perucho et al. (2010a). This work also confirmed that jets with total power L_{j} ≲ 10^{36}Ṁ_{w,−6} erg s^{−1} (Ṁ_{w,−6} ≡ Ṁ_{w}/10^{−6} M_{⊙} yr^{−1}: the wind massloss rate) may be disrupted by the impact of the stellar wind, conditioning the radioband detection of this population of XRBs. Perucho & BoschRamon (2012) analyzed the effect of the wind clumpiness on the jet dynamics and found significant differences in the flow stability and disruption degree even with jets with luminosities L_{j} ∼ 10^{37} Ṁ_{w,−6} erg s^{−1}. Numerical calculations of the highenergy emission from one clumpjet interaction were performed by de la Cita et al. (2017), who predicted that these interactions may collectively dominate the nonthermal radiation. Yoon & Heinz (2015) developed numerical simulations to derive an analytic formula for the asymptotic jet bending and applied it to two wellknown XRBs, Cygnus X1 and Cygnus X3. This work was then extended by Yoon et al. (2016) and BoschRamon & Barkov (2016). More recently, Charlet et al. (2022) analyzed the dynamical and structural effects of radiative losses in the same two fiducial cases. They performed largescale 3D RHD simulations of jet outbreak and early propagation, finding that radiative cooling effects are more relevant for Cygnus X3 than for Cygnus X1. Furthermore, Barkov & BoschRamon (2022) used 3D hydrodynamical simulations to study the combined effect of stellar winds and orbital motion on the scales of the binary and beyond, finding that jets with power L_{j} ∼ 10^{37} Ṁ_{w,−6} erg s^{−1} can be disrupted on scales ∼1 AU.
However, all of these numerical works did not regard the dynamical effect of one fundamental ingredient: the existence of magnetic fields. Toroidal magnetic fields may have a relevant role in jet evolution by either favoring collimation through magnetic hoopstress or by acting as a destabilising agent via the growth of currentdriven instabilities (see, e.g., Martí et al. 2016; Perucho 2019).
Jets are thought to be magnetically dominated close to the compact object, but as the jet propagates away from the binary center, different magnetohydrodynamical processes convert the magnetic energy into kinetic energy. This accelerates the flow, which may be also wind massloaded through mixing produced by smallscale instabilities. The scale of the binary is the main focus of this paper. We expect the kinetic energy at this scale to dominate the total jet power, but it is still unclear to which extent even a relatively weak magnetic field can affect the jet dynamical evolution by favoring jet collimation or triggering more instabilities. Moreover, even moderately weak magnetic fields can be locally reinforced and thus play an important role in shaping the nonthermal emission of the jet by affecting particle acceleration, cooling, and triggering synchrotron radiation. On the other hand, it is also possible that magnetic dissipation may not be as efficient as generally expected, and thus jet magnetic and kinetic or thermal powers are still near equipartition on the scale of the binary.
Previously, relativistic magnetohydrodynamics (RMHD) simulations have concentrated on the formation of jets from MHD mechanisms (McKinney & Blandford 2009; Tchekhovskoy et al. 2011; McKinney et al. 2012; Porth 2013; CruzOsorio et al. 2022) and the morphological characterization of kiloparsec scale magnetized jets (see, e.g., Martí 2019, and references therein). Different authors focused on the properties of stationary onedimensional and 2D axisymmetric relativistic magnetized jet models (Komissarov 1999a; Leismann et al. 2005; Mignone et al. 2005; Keppens et al. 2008; Komissarov et al. 2015; Martí 2015a; MoyaTorregrosa et al. 2021). Mizuno et al. (2007) studied the longterm stability of 3D magnetized spinesheath relativistic jets. The first highresolution numerical simulations of the propagation in three dimensions of an RMHD jet were performed by Mignone et al. (2010). Using similar methods, Guan et al. (2014) studied the propagation of Poyntingfluxdriven jets in AGNs and provided a detailed analysis of energy conversion. Porth et al. (2014) also performed the first 3D RMHD simulations of pulsar wind nebulae with parameters most suitable for the Crab nebula.
In this work, we have performed the first 3D RMHD numerical simulations of magnetized microquasar jets interacting with a moderate to strong stellar wind in an HMXB. Our main goal is to analyze the role of a toroidal magnetic field configuration in the dynamics of a collection of relativistic jets with different parameters. To do this, we compared the results with those obtained in Perucho et al. (2010a) and Perucho & BoschRamon (2012). These simulations can be used as a numerical benchmark upon which to build new highresolution simulations that can be directly compared with observations by radiative postprocess calculations. We plan to address this work in the near future.
The paper is organized as follows: in Sect. 2 we describe the physical scenario that we consider in this work and the numerical setup of the collection of simulations we performed. We also present Lóstrego v1.0, our new 3D RMHD code, and provide the main technical information about the code configuration we use in this paper. In Sect. 3 we present the results of the simulations. In Sect. 4 we discuss the results by comparing our outcome with previous studies of HMXBs jetwind interactions and analyze different mechanisms of energy conversion. In Sect. 5 we summarize the main conclusions.
2. Simulations
For our numerical simulations, we developed a new proprietary code named Lóstrego v1.0, a program that solves the conservative equations of special RMHD in three dimensions using highresolution shockcapturing (HRSC) methods in Cartesian coordinates. This code is fully parallelized with a hybrid scheme with parallel processes (message passing interface, MPI) and parallel threads (OpenMP, OMP), based on the same configuration as the hydro code Ratpenat (Perucho et al. 2010b). This hybrid MPI+OMP scheme exploits the architecture of modern supercomputers and keeps all the cores busy inside each node with OMP instructions. A description of the numerical methods we employed to solve the RMHD equations and the performance of the new code solving classical tests in RMHD can be found in Appendix A. All simulations were performed in Tirant, the supercomputer at the Servei d’Informàtica de la Universitat de València, using up to 1024 cores.
2.1. Physical scenario
We studied the scenario of a relativistic magnetized jet propagating through the clumpy wind of the companion star in an HMXB, similar to the winds of typical XRBs such as Cygnus X1 or LS 5039 (if it really hosts an accreting BH). For simplicity, we located the star in the plane perpendicular to the direction of propagation of the jet, at R_{orb} = 3 × 10^{12} cm from the position of the compact object, and no orbital motion was considered because the characteristic simulation time is much shorter than the typical orbital period of the system. The jet was injected with an initial radius R_{j} = 6 × 10^{9} cm at a distance to the compact object y_{0} = 6 × 10^{10} cm, propagating along the yaxis with an initially mildly relativistic velocity v_{j} = 0.55c. The distance between the injection point and the compact object was enough to guarantee that general relativistic effects (e.g., the curvature of spacetime) can be disregarded from these simulations. For the ambient medium, we considered an inhomogeneous stellar wind with a radial velocity field v_{0} = 2 × 10^{8} cm s^{−1}, a massloss rate ∼10^{−6} M_{⊙} yr^{−1}, a temperature T_{w} ∼ 10^{4} K, thermal pressure p_{w} ∼ 1, 6 × 10^{−3} erg cm^{−3}, and a mean density ρ_{w} ∼ 3 × 10^{−15} g cm^{−3} (we made the assumption that the mean wind density is roughly constant at the scales of jet outbreak within the binary, i.e., y < R_{orb}). These are typical values for a moderate to strong stellar wind from a primary OBtype star (see, e.g., Perucho & BoschRamon 2008; Muijres et al. 2012; Krtička 2014). The stellar wind parameters are based on a simplified model. For example, the wind may still be accelerating when it interacts with the jets, but on the other hand, we did not take the beaming toward the compact object due to the ionization of the accelerating stellar wind close to the star into account (Molina et al. 2019; Vilhu et al. 2021). The characterization of these nontrivial effects, among other complexities of the wind structure, requires dedicated simulations that are beyond the scope of this paper. The stellar wind filled the computational box from the beginning, and we did not account for radiation pressure from the star or other related physical effects. We also considered that the wind was dominated by the kinetic component at the scales of the binary, so we neglected the magnetization of the ambient medium. We modeled the plasma as an ideal gas with Γ = 5/3, and we neglected both thermal and nonthermal cooling.
2.2. Numerical setup
We simulated three different magnetized jets in order to study the effect of the magnetic field in the jet evolution inside the binary. The physical parameters of each simulation are listed in Table 1. The total luminosities of jet A and jet B are ∼10^{35} and ∼10^{37} erg s^{−1}, respectively, which is similar to the luminosities studied in Perucho et al. (2010a). Jet A and jet B are cold (h ≈ 1) and kinetically dominated, such that the internal and magnetic energy fluxes that are injected in the grid represent a small fraction of the total jet luminosity (i.e., L_{B} ≈ 0.004 × L_{h}, where L_{h} is the sum of the kinetic and internal energy powers, and L_{B} is the magnetic power). The total luminosity of jet C is similar to that of jet B, but in this case, the jet is set in equipartition between kinetic or thermal and magnetic energy fluxes. For all of our models, the β ratio (i.e., the ratio of the average magnetic pressure, , to the average gas pressure, ), remains close to unity, although magnetic pressure dominates gas pressure by a factor of 1.5 in jet C. The initial densities were ρ_{A} = 0.088 ρ_{w} for jet A, ρ_{B} = 8.8 ρ_{w} for jet B, and ρ_{C} = 0.88 ρ_{w} for jet C. Because jet C is more dilute and the total pressure (gas+magnetic) is higher than in jet B, the flow is initially hotter (h ≈ 1.66).
Summary of the main parameters of the three jets.
We performed the simulations in a numerical grid box with physical dimensions of 80 × 240 × 80 R_{j} (in units of the jet radius at the injection plane) with an initial resolution of 6 cells/R_{j}, so that the box involves 480 × 1440 × 480 cells. The resolution of the grid improves the resolution used in Perucho et al. (2010a) and Perucho & BoschRamon (2012) by a factor of 1.5, but we only considered 75% of the longitudinal length in Perucho et al. (2010a) to keep the computational requirements within reasonable limits. No grid extensions or mesh refinement were employed in our simulations. As the jets expanded in the numerical box, the effective resolution of the simulation was increased because the transversal section of the jet involves more computational cells.
Timeresolved spectroscopic monitoring of hot star emission lines supports the idea that the stellar wind is formed by clumps. These clumps are density inhomogenities that are created as a result of thermal instabilities during the wind acceleration process (Runacres & Owocki 2002; Cassinelli et al. 2008; Moffat et al. 2008). Moreover, the similarity of the scaling laws between these clumps and the molecular clouds favors the interpretation that the clumpy structure of the wind is also a result of turbulence, with scaling laws similar to those of the fractalized ISM (Moffat et al. 2008). Thus, the initial inhomogenous density distribution of the stellar wind was constructed using the publicly available PyFC code^{1} (Wagner & Bicknell 2011), which is a Python package that is useful to represent a dense, inhomogenous component embedded in a smooth background. This method is different from the approach adopted in Perucho & BoschRamon (2012), where clouds were modeled as Gaussianshape spheres that are randomly distributed in the computational box. The distribution of inhomogeneities was established by an iterative process following the original work on terrestrial cloud models by Lewis & Austin (2002). The 3D density field followed a lognormal singlepoint statistics in real space, while the fractal structure of the system was achieved by first Fourier transforming and then multiplying the computational box by a power law with a Kolmogorov spectral index β = −5/3. The mean (μ) of the lognormal parent distribution was μ = 1.0 and the variance σ^{2} = 0.6, such that the minimum density of the wind is ρ_{min} ≈ 0.1 ρ_{w} and the maximum density is ρ_{max} ≈ 10 ρ_{w}. Based on the dimensions of the computational box, we set a minimum sampling wavenumber k_{min} = 3, so that the radius of the largest fractal structure in the cube (i.e., the maximum cloud size) was approximately R_{c, max} = 6.5 R_{j} (R_{c, max} < 2 R_{⊙}, ∼20% of the stellar radius).
In each simulation, the jet was injected into the y = y_{0} plane (in the system of coordinates with origin in the compact object) as a cylindrical nozzle with (where length units are normalized to the jet injection radius). In this plane, density (ρ) and velocity (v) profiles are given by
where θ, ϕ are the characteristic angles of spherical coordinates with respect to the location of the star,
and the subscript j = 1 indicates that the values of density and velocity in the yboundary cells outside the nozzle r < 1 are a copy of the values in the first plane inside the box. Beyond the cylindrical nozzle (i.e., r > 1), we distinguish two different inlet boundaries: the jet cocoon, a dynamical region in which boundaries reflect to simulate the presence of a counterjet, and the stellar wind, following the radial velocity field with origin in the location of the star (see Fig. 1). To distinguish between these two regions during the simulation, the jet cocoon was defined as the set of cells for which the gas pressure p_{g} > p_{c}, where we chose the threshold p_{c} = 10^{−8} (in units of ρ_{0}c^{2}). The azimuthal component of the jet axisymmetric magnetic field in the laboratory frame has the form (Lind et al. 1989; Komissarov 1999b; Leismann et al. 2005; Martí 2015a)
Fig. 1. Scheme of the simulation setup in the XY plane. The radial velocity field that represents the stellar wind fills the box from the beginning. We also show a top view of the injection plane at Y = 0, where we distinguish the cylindrical nozzle (blue), the shear layer (red), and the jet cocoon (orange), where boundaries are reflecting. 
where R_{Bϕ, m} is the magnetization radius normalized to the jet radius, and is the maximum value of the magnetic field. Thus, the magnetic field grows linearly for r < R_{Bϕ, m}, reaches a maximum at r = R_{Bϕ, m}, and decreases as 1/r for r > R_{Bϕ, m}. Hereinafter, we fixed R_{Bϕ, m} = 0.37, and we chose for each simulation according to the magnetic power of the jet (see Table 1),
In order to avoid the appearance of nonphysical magnetic monopoles, boundaries were always free with respect to the magnetic field vector components. For jets without rotation (v^{ϕ} = 0), the gas pressure profile can be derived from a single ordinary differential equation for the transversal equilibrium across the beam,
where W is the Lorentz factor. When we integrate Eq. (8) by separation of variables, the gas pressure profile p(r) yields
where
In Eq. (10), represents the value of the toroidal magnetic field evaluated at the jet radius, and we assumed that B^{r} = B^{y} = 0. This assumption is supported by two main arguments: (1) an axial component contained in the boundary zone would lead to open field lines and thus to a violation of the divergencefree condition, and (2) the injection scale likely exceeds the jet acceleration region, and so the magnetic field is expected to be mainly toroidal. The total pressure in r = R_{j}, p_{a}, can be derived from the hydrodynamic (kinetic+internal) jet luminosity equation,
where h(r) is the specific enthalpy. The stellar wind enters the grid from the XMAX boundary, while in all the remaining boundaries of the computational box (i.e., XMIN, YMAX, ZMIN, ZMAX), we considered freeflow conditions.
The simulations presented in Sect. 3 were performed with Lóstrego v1.0, a new conservative, finitevolume 3D(S)RMHD code in Cartesian coordinates. The configuration of the code is based on our testing benchmark results (see Appendix A). We employed a secondorder Godunovtype scheme with the HLLD Riemann solver (Mignone et al. 2009) and the piecewise linear method (PLM) for cell reconstruction, with the VAN LEER slope limiter (van Leer 1974; Mignone & Bodo 2006). The limiter was degraded to MINMOD (Roe 1986) when strong shocks are detected in order to avoid spurious numerical oscillations around shocks (Mignone & Bodo 2006). When shockflattening is applied to cell reconstruction, the HLLD Riemann solver is also degraded to the simpler and more diffusive HLL solver. The advance in time was performed using the thirdorder TVDpreserving RungeKutta (Shu & Osher 1989) with CFL = 0.2. The relativistic correction algorithm CA2 of Martí (2015b) was used to correct the conserved variables after each time iteration. In our test simulations, this scheme was tested to be robust even when the magnetic pressure dominates the gas pressure by more than two orders of magnitude. The magnetic field divergencefree constraint is preserved with the constrained transport (CT) method (Evans & Hawley 1988; Ryu et al. 1998; Balsara & Spicer 1999), where electromotive forces at cell corners are interpolated following Gardiner & Stone (2005). An additional equation that describes the advection of a tracer function f was included in the system of equations to indicate the composition of the fluid in every cell of the box as a function of time. This tracer was used to distinguish between jet material (f = 1), wind material (f = 0), and the regions in which jet has been mixed with the environment (0 < f < 1).
In RMHD jet simulations, the plane of injection is one of the most challenging regions of the grid and is critical for the longterm survival of the whole simulation. This might become especially relevant when jets are injected as axisymmetric tophat distributions with low numerical resolution in Cartesian coordinates. For example, we have found from analyzing the first seconds of jet propagation that the axial velocity is slightly perturbed along the jet surface. This perturbation is translated into the toroidal magnetic field, which starts to develop an axial component that grows with time and may lead to local nonphysical solutions. Porth (2013) also described the introduction of a significant amount of noise caused by the Cartesian discretization and quadrantal symmetry of the grid, which led to the pump of multiples of the m = 4 mode in all flow quantities (see Sect. 3.2 in Porth 2013). These pathologies are controlled in our simulations first by using inhomogenous stellar winds that naturally break the original Cartesian symmetry, and second, by smoothing the initial tophat field distribution. In order to smooth the initial profile to avoid the growth of random perturbations at the jet base, we replaced the discontinuous functions that we described above by smooth functions of the form (Bodo et al. 1994)
To guarantee that the injected toroidal field has zero divergence up to machine accuracy, we fixed one of the field components according to its analytical expression (smoothed with the shear layer of Eq. (12)) and calculated the other component numerically using the solenoidal condition (i.e., ∇ ⋅ B = 0), where spatial derivatives are approximated with finite differences. For the sake of completeness, a schematic representation of the whole setup is shown in Fig. 1.
Units are used in which the light speed (c), the mean density of the stellar wind (ρ_{w}), and the jet radius (R_{j}) are set to unity. A factor of is absorbed in the definition of the magnetic field, so that the actual field is smaller by .
3. Results
To illustrate the early evolution of a fiducial jet model of our three simulations, Fig. 2 shows the restmass density distribution of jet B within a small computational box soon after injection (after the first 20 R_{j} of propagation through the ambient medium). In the plot, we show a 3D cut of the logarithmic density of the ambient medium, a volume render of the tracer function and a gas pressure contour to show the position of the jet bow shock. When the jet starts to propagate, the highly supersonic gas generates a strong forward shock that pushes the ambient medium and a reverse shock that decelerates the jet gas at the terminal region (i.e., the hot spot), deflecting backward off the jet flow and forming a hot and light cavity surrounding the jet called cocoon. In the following, we refer equally to cocoon, shocked cavity, or jet cavity. Surrounding the cocoon is a dense shell of shocked gas that is separated from the cocoon by a contact discontinuity and is isolated from the unshocked ambient medium by the bow shock. In these simulations, the cocoon and the shear layer are both subject to the development of irregular structures through the interaction of the bow shock with the clumps of the stellar wind. In Fig. 3 we show that the toroidal magnetic field is advected following the jet gas through the clumps of the stellar wind, represented with gold isosurfaces at ρ = 3 ρ_{w}. Field vectors of Fig. 3 were integrated to represent a collection of magnetic field lines in the hot spot of Fig. 2 (solid black lines), where the density of the lines and the restmass density reach their maximum.
Fig. 2. Propagation of a fiducial jet (jet B) during the first seconds after injection. The dimensions of the box are , with an effective resolution of 6 cells/R_{j}. The density distribution of the ambient medium is represented in logarithmic scale, where we limited the maximum density to ρ_{max} = 10 ρ_{w} to highlight the wind clumps. Magnetic field lines are represented in the head of the jet, where the density is highest. The 3D jet render is constructed using the jet tracer function and is colored with the same colour scale as the density distribution (although we do not show the tracer legend to avoid redundancy), where jet particles (f = 1) are represented in light blue and white and the ambient medium (f = 0) in yellow and white. A gas pressure contour (faint yellow) is also included to show the position of the jet bow shock. 
Fig. 3. Toroidal magnetic field vectors in the same time frame as in Fig. 2. Stellar wind clumps are represented with gold isosurfaces at ρ = 3 ρ_{w}. 
3.1. Jet A
Jet A propagates up to y ≤ 220 R_{j} at t = 1020 (R_{j}/c), which corresponds to y ≤ 1.3 × 10^{12} cm at t ≈ 200 s in physical units. This means that the propagation velocity of the jet head through the stellar wind is v_{j} ≈ 0.21 c, slightly higher than the onedimensional relativistic approximation of Martí et al. (1997) The relativistic velocity of the jet head can be estimated based on a onedimensional hydrodynamical momentum balance throughout the working surface, v_{1D} ≈ 0.15 c. This small difference in the velocity of propagation can be explained by considering the ballistic nature of jet A, as described below.
In the vertical panels of Fig. 4, we show the logarithms of restmass density, gas pressure, and magnetic pressure at t = 1020 (R_{j}/c). A collection of ten tracer contours is overplotted together with the logarithm of gas pressure to show the position of the jet core and the composition of the jet cavity. In the horizontal panels of this figure, we show the distribution of axial velocity at three different time frames: t = 310, t = 610, and t = 1020 (R_{j}/c). In the top panel of Fig. 5, we show a 3D render of the jet morphology and the clumpy wind distribution after jet propagation through the numerical domain, using the jet tracer and the restmass density, respectively. We also include a gas pressure contour (in faint yellow) to show the position of the jet bow shock at the last time frame. The overall structure of the jet is highly ballistic in the timescales of the simulation, but the interaction with the denser clumps of the stellar wind deforms the bow shock and produces an irregular shell. Near the base of the jet, the cavity seems to be inflated because of the effect of reflecting boundary conditions inside the shocked gas, producing a small equatorial bulge. The jet central spine shows several pinches due to weak reconfinement shocks, but jet collimation is preserved during the whole evolution. The jet bends slightly to the left (in the XY plane of Fig. 4) due to the lateral impact of the stellar wind and the subsequent difference in the total pressure on either side of the jet. Instabilities caused by this asymmetry trigger some degree of mixing between the jet material and the ambient gas at the jet boundaries close to the head. In this simulation, the jet spine maintains the injection velocity until the terminal region (i.e., v_{j} ≃ 0.55 c). In the bottom panel of Fig. 5, we show a collection of density isocontours that represent the jet cocoon (ρ = 0.1 ρ_{w}, red), the unperturbed lowdensity clumps of the stellar wind (ρ = 2.5 ρ_{w}, gold), and the clumps compressed by the bow shock at two density levels: ρ = 5 ρ_{w} (dark blue) and ρ = 7.5 ρ_{w} (light blue). The toroidal magnetic field lines remain anchored to the jet central spine, whereas the gas that fills the cocoon drives a largescale mildly entangled toroidal to helical field structure. At y ≤ 80 R_{j}, magnetic field lines near the jet surface are subject to shear effects, which lead to the development of a poloidal component. This effect, which explains the melted appearance of the toroidal field in this model, might be produced by the combined effect of numerical and possibly nonphysical instabilities and the irregular structure of the poloidal velocity profile.
Fig. 4. Axial cuts of jet A simulation. Vertical panels: Logarithmic restmass density, logarithmic gas pressure, and logarithmic magnetic pressure at t = 1020 (R_{j}/c). Ten tracer contours are overplotted together with the gas pressure, from f ≈ 0 (white to blue) to f = 1 (yellow to white). Horizontal panels: evolution of the jet velocity at t = 310 (top), t = 610 (middle), and t = 1020 (bottom) (R_{j}/c). 
Fig. 5. 3D model of jet A simulation. Top panel: 3D render of tracer at t = 1010 (R_{j}/c) and stellar wind clumps. A gas pressure contour (faint yellow) is included to show the position of the jet bow shock. Bottom panel: density isocontours at ρ = 0.1 ρ_{w} (red), ρ = 2.5 ρ_{w} (gold), ρ = 5.0 ρ_{w} (dark blue), and ρ = 7.5 ρ_{w} (light blue) together with a collection of magnetic field lines. 
3.2. Jet B
The propagation of jet A and jet B is similar from a qualitative point of view, although important dynamical differences are triggered by the asymmetric impact of the stellar wind and the reconfinement shocks produced within the jet. Figure 6 is similar to Fig. 4, but it shows the solution of jet B at t = 560 (R_{j}/c) and three frames of the jet axial velocity along its evolution, at t = 110, t = 360, and t = 560 (R_{j}/c). This jet propagates rather ballistically up to y ≤ 220 R_{j} at t = 560 (R_{j}/c), which corresponds to y ≤ 1.3 × 10^{12} cm at t ≈ 110 s, in physical units. This means that the propagation velocity of the jet head through the stellar wind is v_{j} ≈ 0.39 c, which agrees well with the onedimensional relativistic estimate, v_{1D} ≈ 0.43 c.
Fig. 6. Axial cuts of jet B simulation. Vertical panels: Logarithmic restmass density, logarithmic gas pressure, and logarithmic magnetic pressure at t = 560 (R_{j}/c). Ten tracer contours are overplotted together with the gas pressure, from f ≈ 0 (white to blue) to f = 1 (yellow to white). Horizontal panels: Evolution of the jet velocity at t = 110 (top), t = 360 (middle), and t = 560 (bottom) (R_{j}/c). 
Because the total power of jet B is two orders of magnitude larger than for jet A, the effect of the impact of the lateral wind on the longterm evolution of the jet is clearly smaller and the symmetry of the cocoon is almost preserved. The interaction of the backflow and the shocked ambient medium in the shear layer leads to the development of KelvinHelmholtz instabilities and turbulence, but the jet core is unmixed (f ≈ 1) and maintains the velocity of the injection point with good accuracy until the end of the simulation (i.e., v_{j} ≃ 0.55 c). The cocoon is more elongated than in jet A because the jet head advances faster. In this case, the jet is initially denser than the average ambient medium (i.e., ρ_{j} = 8.8 ρ_{w}), and the total pressure is also higher. This produces an initial lateral expansion near the base with a corresponding drop in total pressure, which finally leads to underpressure with respect to the cocoon and to the generation of a recollimation shock far from the injection plane. The position of the shock changes from y ≈ 120 R_{j} at t = 360 (R_{j}/c) to y ≈ 160 R_{j} at t = 560 (R_{j}/c), meaning that it moves with v_{shock} ≈ 0.2 c. At the terminal shock, the magnetic pressure also increases because the field lines are compressed. Near the reconfinement shock region, jet material cannot flow freely and accumulates in the head of the jet, leading to backflow blobs of jet gas that drag the magnetic field filling the cocoon (Fig. 7, bottom panel). As the shock moves downstream, the stretched nozzle deposits more plasma in the cocoon and creates a filamentous structure, where f < 1 because of mixing with the ambient medium gas.
Fig. 7. 3D model of jet B simulation. Top panel: 3D render of tracer at t = 560 (R_{j}/c) and stellar wind clumps. A gas pressure contour (faint yellow) is included to show the position of the jet bowshock. Bottom panel: density isosurfaces at ρ = 0.1 ρ_{w} (red), ρ = 2.5 ρ_{w} (gold), ρ = 5.0 ρ_{w} (dark blue), and ρ = 7.5 ρ_{w} (light blue) together with a collection of magnetic field lines. 
3.3. Jet C
Even though the total luminosity of jet C is the same as in jet B, the jet dynamics are vastly different, as are evolution and overall morphology. This is triggered by the moderate to strong magnetic fields within the jet. Figure 8 is similar to Figs. 4 and 6, but it shows the solution of jet C at t = 720 (R_{j}/c), and three panels of the axial velocity distribution at times t = 210, t = 490, and t = 720 (R_{j}/c) to illustrate the jet evolution. In this simulation, the jet propagates up to y ≤ 200 R_{j} at t = 720 (R_{j}/c), which corresponds to y ≤ 1.2 × 10^{12} cm at t ≈ 142 s, in physical units. This means that the propagation velocity of the jet head through the stellar wind is v_{j} ≈ 0.26 c, which is lower than the onedimensional relativistic estimate of Martí et al. (1997), v_{1D} ≈ 0.33 c. It is important to note, however, that this simulation cannot be prolonged in time because the bow shock has touched the lateral walls of the grid at t < 720 (R_{j}/c). Because the magnetic field is initially larger than in the two previous models (as it is required to achieve power equipartition with the energy flux associated with the particles), the initial overpressure of the jet nozzle with respect to the medium is even higher than in jet B, leading to a pronounced lateral expansion and a quick, strong recollimation. During this quick process, the jet flow increases its bulk velocity from the mildly relativistic value at injection (i.e., v_{j} ≃ 0.55 c) to a maximum of v_{j} ≃ 0.95 c at the jet core. The difference with respect to the onedimensional estimate is even more evident if we calculate the velocity of propagation considering the maximum velocity after the first quick expansion (i.e., v_{j} ≃ 0.95c), yielding v_{1D} ≈ 0.76 c (a factor ∼3 higher than the actual velocity of propagation).
Fig. 8. Axial cuts of jet C simulation. Vertical panels: Logarithmic restmass density, logarithmic gas pressure, and logarithmic magnetic pressure at t = 720 (R_{j}/c). Ten tracer contours are overplotted together with the gas pressure, from f ≈ 0 (white to blue) to f = 1 (yellow to white). Horizontal panels: Evolution of the jet velocity at t = 210 (top), t = 490 (middle), and t = 720 (bottom) (R_{j}/c). 
The effect of the toroidal field on the jet dynamics is visible from the early jet evolution by the development of currentdriven instabilities that destabilize the jet beyond the first reconfinement shock. The consequent increase in crosssection in which the head deposits its momentum slows down the velocity of propagation of the head with respect to jet B. The morphology and dynamics of the cocoon are also different than in the other two simulations (Fig. 9, top panel) for the same reason. The first reconfinement shock drives a chain of shocks that are evident throughout the entire jet spine, where magnetic pressure reaches maximum within the jet due to the compression of toroidal field lines. At t = 720 (R_{j}/c), at least four reconfinement shocks have appeared after jet evolution, at y ≈ 35 R_{j}, y ≈ 70 R_{j}, y ≈ 100 R_{j}, and y ≈ 130 R_{j}. This means that at the end of the simulation, there is a shock at approximately every y ≈ 30 R_{j}. The backflow interacts strongly with the shocked ambient medium, developing KelvinHelmholtz instabilities that lead to turbulence and mixing with the jet gas within the cocoon (f < 1). The morphology of the magnetic field lines is also different in the lower body of the cocoon and in the jet head (Fig. 9, bottom panel): at y ≤ 100 R_{j}, the magnetic field remains roughly ordered and the field lines preserve the toroidal morphology. However, at y > 100 R_{j}, the field lines are highly tousled, forming an entangled lobe structure around the head of the central axis. This structure of the magnetic field is translated into the morphology of the cocoon material, where gas exhibits toroidal filaments at y ≤ 100 R_{j}, but it appears heavily disordered downstream from the chain of shocks (Fig. 9, top panel).
Fig. 9. 3D model of jet C simulation. Top panel: 3D render of tracer at t = 720 (R_{j}/c) and stellar wind clumps. A gas pressure contour (faint yellow) is included to show the position of the jet bow shock. Bottom panel: density isosurfaces at ρ = 0.1 ρ_{w} (red), ρ = 2.5 ρ_{w} (gold), ρ = 5.0 ρ_{w} (dark blue), and ρ = 7.5 ρ_{w} (light blue) together with a collection of magnetic field lines. 
4. Discussion
4.1. Jet power and jet morphology
In our RMHD simulations, the total jet power was chosen to be comparable to those of the RHD simulations performed in Perucho et al. (2010a). Jet A and jet B have a power similar to that of jet 1 and jet 2, respectively. We focus first on the case of jet A compared to jet 1. Although the magnetized jet propagates up to almost the same distance (y ≃ 1.3 × 10^{12} vs. y ≃ 1.5 × 10^{12} cm), there is a remarkable difference in the amount of time that these two models take to achieve this distance: whereas jet A needs approximately 200 s, jet 1 takes 1250 s (i.e., a difference of a factor ∼6). The reason for this difference is related to the destabilisation and disruption of jet 1 beyond a moderate to strong recollimation shock at y ≃ 5 × 10^{11} cm. The jet flow is strongly decelerated following that shock, and it is thus prone to the fast development of helical instabilities triggered by the impact of the stellar wind. In jet A, magnetic tension of the toroidal field acts against expansion (furthermore, the gas pressure profile is set up initially to ensure transversal equilibrium with the magnetic pressure), whereas in jet 1, the jet expands only as determined by the jettococoon total pressure ratio.
Although jet B is also slightly faster than jet 2 in Perucho et al. (2010a) (the first reaches y ≃ 1.3 × 10^{12} cm in t ≃ 110 s, and the second takes almost twice as long to reach y ≃ 2 × 10^{12} cm), the difference is reduced precisely because jet 2 is more collimated and develops a very oblique recollimation shock, which does not decelerate the flow enough to favor its disruption, in contrast to jet 1. Both jets are fairly similar in their evolution, with the minor difference in both advance velocity and location of the first recollimation shock (which is slightly farther downstream in jet 2), probably due to the longer simulation time in the RHD simulation. Nevertheless, the hydrodynamical jet shows a lower degree of collimation than the RMHD one (see Fig. 7). Again, although in this case the magnetic tension does not seem to play such an obvious role as in the case of jet A, its effect can be crucial for jet evolution and longterm stability. Whether this difference is related to the development of a poloidal sheared component around the jet flow (see the bottom panel in Fig. 7) should be a matter of study, but it is beyond the scope of this paper.
Overall, we find that the presence of a toroidal field, as long as it is nondominant from an energetic point of view, might contribute to stabilizing the jet evolution in microquasars. Nevertheless, a stronger magnetic field (as is the case for jet C) causes the jet to become prone to the development of fast recollimation shocks and currentdriven instabilities, which do not develop in the case of jet A and jet B. Therefore, future work should address the threshold from which the magnetic field can become a destabilizing factor for the jet structure at the scales of the binary. This comment is of remarkable relevance in terms of the longterm stability of relativistic jets not only for microquasars, but also for AGNs: a nondominant toroidal field, as is expected to be the case beyond collimation and acceleration scales in both microquasar and AGN jets, could have a stabilizing role with respect to purely hydrodynamical jets. The reasons are magnetic tension to avoid rapid expansion and largeangle recollimation shocks with the consequent flow deceleration on the one hand, and the possible generation of a sheared poloidal component in the backflow that shields the jet against the development of instabilities on the other hand (although KelvinHelmholtz instabilities, mainly helical, can still develop).
Interestingly, jet C is more similar to jet 1 (which is two orders of magnitude less luminous) than to jet 2 in Perucho et al. (2010a) or jet B in this paper, which stresses the relevance of the development (or lack of it) of instabilities in outflows in terms of energy dissipation and the subsequent kinematics and morphology. This jet resembles the structures observed in axisymmetric simulations (see, e.g., Martí et al. 2016; MoyaTorregrosa et al. 2021) for the case of jets in which the Poynting flux is relevant or dominates the total energy flux. In this type of models, the presence of a chain of recollimation shocks forced by the toroidal field favors the dissipation of kinetic energy. However, it is important to note that this process takes place after the jet has reached higher velocities than are expected in microquasar jets, which is precisely a consequence of the acceleration provided by the presence of strong fields and high internal energies at injection. Jets with higher kinetic energies or inertia are expected to be less prone to developing disrupting instabilities, at least in the nonmagnetized case (Perucho et al. 2005). Despite the large inertia achieved by the jet after acceleration, kinetic energy is efficiently dissipated at the series of recollimation shocks. These shocks therefore destabilize the flow not only by themselves, but also by making the jet sensitive to the effect of the lateral wind, as reported in Perucho et al. (2010a) for jet 1. Whether microquasar jets eventually reach such velocities and are later decelerated at larger distances (beyond the binary system) should be tested by comparison of the synthetic radiative output from simulations and observational features from known sources.
Finally, we highlight that jet 1 and jet 2 in Perucho et al. (2010a) faced a homogeneous medium, but this is not the case of jet A and jet B in this paper. However, our results show that the differences caused by clumps are likely minor at the timescales of our simulations. In contrast, in Perucho & BoschRamon (2012) the authors adopted a setup in which the jet was established in equilibrium with the ambient medium, a simplified version of the clumpy stellar wind. The stage of direct windjet interaction is only plausible when the forward bow shock has propagated far enough from the simulated binary region, that is, at longer timescales than those studied here.
4.2. Energy channels
The analysis of the energy distribution and the evolution of the energy fluxes along the jet can shed light on the interplay among the different energy channels (i.e., kinetic, internal, and magnetic) and, in particular, the possibility of jet acceleration via extraction of magnetic and/or internal energies. Furthermore, it is well understood that the relativistic nature of jets favors the very efficient exchange of a good fraction of the kinetic and internal energy fluxes with the ambient medium, mediated by strong shocks. This effect was previously described based on 2D axisymmetric models (Perucho et al. 2011, 2014) and was extended to longterm 3D hydrodynamical simulations in Perucho et al. (2019) for AGN jet evolution. However, this effect has not been explored so far with magnetized jet models or in the particular context of microquasars.
The accumulated energies in each of the three channels described above can be computed numerically considering the relativistic version of the energy equations and summation over all cells in the 3D numerical box,
where τ_{k}(t) is the kinetic energy, τ_{ϵ}(t) is the internal energy, τ_{B}(t) is the magnetic energy (all expressed in the observer’s reference frame), and ξ represents the jet tracer, ξ = f, to calculate the energy of jet particles, or ξ = 1 − f, to calculate the energy transferred to the ambient medium. Finally, Δx, Δy, and Δz are the cell sizes for each spatial dimension. Considering that the wind is stationary and the bow shock has not crossed the outer boundaries of the numerical grid, the sum of the different energies should be equal to the total energy injected by the jet nozzle at a given time, t (plus the internal and kinetic energy of the stellar wind).
For the purposes of this section, we restrict our analysis to simulations of jet B and jet C because these two models have an equivalent total jet power, but the energy flux is distributed in a different way among the three energy channels; while in jet B the kinetic energy flux dominates the internal energy flux and the magnetic energy flux by almost two orders of magnitude, in jet C the magnetic flux is initially in equipartition with the sum of internal and kinetic energy fluxes (whereas the internal energy flux is also a factor ∼5 higher than the kinetic energy flux). Figure 10 shows the time evolution of the logarithm of the different types of energies for jet B (left panel) and jet C (right panel) at t = 110, 360, 560 (R_{j}/c) and t = 10, 110, 210, 490, 720 (R_{j}/c), respectively. In the case of jet C, we included two additional time frames to those shown in Fig. 8 in order to analyze energy conversion during the early moments of jet propagation, because some of these processes occur very quickly after the injection of the jet. In Fig. 10 the solid lines (without marks) show the logarithm of the energy injected in the computational box as a function of time. We recall that the inclusion of a shear layer to smooth the initial tophat axial velocity and magnetic field profiles as described in Sect. 2.2 implies that the eventually injected power is lower than the theoretical value given in Table 1. This is translated into a power reduction of a ∼20% for the simulation of jet B and ∼10% for jet C. The difference in the percentages of the two models can be explained by the higher sensitivity of kinetic energy flux to our shear layer than the magnetic and internal energy fluxes.
Fig. 10. Time evolution of the logarithm of the energy in the simulation for jet B (left) and jet C (right) at the three time frames of Figs. 6 and 8, respectively. In the right panel (jet C), we have included two additional time frames at t = 10 (R_{j}/c) and t = 110 (R_{j}/c). Energy stored by the jet plasma (f ≠ 0) is represented with dots and dashed lines, ambient medium energies (f = 0) with dots and dashdotted lines, and total energies (jet+ambient medium) with squares and solid lines. Continuous solid lines (without marks) are the theoretical curves (i.e., the injected values) for each type of energy: kinetic (red), internal (black), magnetic (blue), and total (green). Energy units are . 
As said before, the sum of the jet and ambient medium energies accumulated in the numerical grid must be equal to the injected energy plus the original internal and kinetic energies of the ambient medium. The small difference (a few percent; see the next paragraphs) comes from the use of a correction algorithm for the conserved variables (see the appendix) and, more importantly, the low resolution of the numerical representation of the circular injection nozzle by a Cartesian grid. On the other hand, assuming that no energy conversion between the different channels (internal, magnetic, and kinetic) takes place within the jet, and that the only transfer of energy within each of those channels is that between the jet and the ambient medium, the sum of the jet and the ambient energies (solid lines with square marks) should be close to the theoretical injected energy (solid lines) for each energy channel. Thus, the difference between the numerical values for the total energy in each channel and the theoretical curves implies that different physical processes of energy conversion are operating during jet evolution. We describe these processes for the two jet simulations in the remaining paragraphs of the section.
Jet B. In the left panel of Fig. 10, we show that total energy (solid green line with squares) is roughly conserved during the evolution from t = 110 to t = 560 (R_{j}/c) because the three points we considered for this analysis overlap (with an error smaller than 1%) with the total injected energy (solid green line without marks). Although the total kinetic energy (solid red line with squares), total internal energy (solid black line with squares) and total magnetic energy (solid line with squares) do not coincide with their respective injected energies (because of effective energy conversion), the initial global energy distribution is preserved: kinetic energy dominates the whole evolution, and magnetic energy remains residual (below 1% of the total energy) even when the jet is fully developed. The jet magnetic energy (dashed blue line with dots) is lower than the injected value even at early t = 110 (R_{j}/c) (we also confirmed that this is true for earlier times). Discarding a conversion of magnetic energy into kinetic energy (because the jet is kinetically dominated), the origin of the difference between the total magnetic energy and the injected energy can be attributed to the numerical resistivity of the code, which can produce a continuous conversion of magnetic into internal energy.
On the other hand, the blue points connected with dashdotted lines shown in the left panel of Fig. 10 represent the total magnetic energy accumulated in the ambient medium. Optimally, this energy would have to be zero because in ideal MHD, there is no physical mechanism to transfer magnetic energy from the jet fluid into the nonmagnetized ambient medium. However, it is important to note that the magnetic energy assigned to the ambient medium comes mainly from numerical cells with a mixture of jet and ambient fluids, where according to Eq. (13) the magnetic energy is distributed arbitrarily proportional to the corresponding mass fraction. Then, the amount of magnetic energy that is accumulated in the ambient medium (∼10% of the injected magnetic energy at the end of the simulation) may be interpreted as an ample upper bound of the magnetic diffusion of our numerical code.
In this simulation, the main manifestation of energy transfer concerns the conversion of jet kinetic energy (dashed red line with dots) into kinetic and internal energy of the ambient medium (dashdotted red and black lines with dots, respectively), mediated by the propagation of the bow shock, and into internal energy of jet plasma in the shocked cavity (dashed black line with dots). At the end of the simulation, the jet plasma preserves more than 70% of the injected kinetic energy. At the same time, the transferred portion (∼30%) has been used to accelerate (∼10%) and heat (∼20%) the shocked stellar wind, and to increase the total internal energy of the jet plasma by a factor ∼ 2.
Jet C. In the right panel of Fig. 10, we show the energy distribution for jet C simulation between t = 10 and t = 720 (R_{j}/c), where we included two additional time frames at t = 10 and t = 110 (R_{j}/c) with respect to those shown in the horizontal panels of Fig. 8. The total energy is conserved during the evolution with an error smaller than 1%^{2}, whereas the magnetic energy assigned to the ambient medium by the end of the simulation is lower than 5% of the total injected magnetic energy.
The evolution of the energy distribution in this simulation is strongly conditioned by the physical processes operating at very early times. At t = 10 (R_{j}/c), magnetic or internal and kinetic energies are already lower and higher than their injected values, respectively, with τ_{B, jet}/τ_{B, inj} ≈ 0.40, τ_{ϵ, jet}/τ_{ϵ, inj} ≈ 0.60 and τ_{k, jet}/τ_{k, inj} ≈ 3, meaning that about 60% of the jet magnetic energy and 40% of the internal energy have been converted into jet kinetic energy by that time. Thermal acceleration occurs in hot, expanding flows by the Bernoulli process. Magnetic acceleration, in contrast, needs a particular expanding flow geometry (differential collimation; see, e.g., Vlahakis 2004; Vlahakis & Königl 2004; Komissarov 2011, and references therein). According to the results of Komissarov et al. (2007) for axisymmetric models, magnetic acceleration is very efficient with as much as ∼75% of the Poynting flux converted into kinetic energy at long enough spatial scales. The fast expansion experienced by get C close to the injection nozzle (see Fig. 8) could provide the correct dynamical and geometric conditions to trigger the magnetic acceleration process, which would then have been observed in microquasar jet simulations for the first time and at much shorter spatial scales.
The plot of the evolution of the energy accumulated in the different channels, with the curves overtaking each other several times, reveals the complex nature of jet propagation in model C. As an example, we note that the internal energy of the jet dominates the magnetic and kinetic energies between t = 10 (R_{j}/c) and t = 210 (R_{j}/c), but at t = 490 (R_{j}/c), it has exchanged its position with the kinetic energy. Additionally, the figure shows that the magnetic channel loses energy steadily with time. The decrease of magnetic energy can take place within the jet at the magnetic acceleration sites after the reconfinement shocks (where the flow expansion could admit the differential collimation), but also at the magnetically dominated regions of the turbulent shocked cavity (see Fig. 11).
Fig. 11. Fraction of magnetic energy density with respect to total energy density for jet C in the XY (top) and ZY (bottom) plane. 
Finally, we observe a relevant difference between models B and C in terms of kinetic and internal energy distribution: in the case of jet B, most of the energy in the grid is hosted by the jet matter, whereas the opposite seems to be the case for jet C. It has been shown in a series of papers (Perucho et al. 2011, 2014, 2019) that relativistic jets are very efficient in transferring energy to the ambient medium. In Perucho et al. (2017), the authors attributed this efficiency to the relativistic nature of jets and showed that massive, cold jets are less efficient in the process of energetic transfer because a large fraction of the energy flux is in the form of kinetic energy, dominated by the mass of the particles. This agrees with the results we derived here, where jet B is significantly denser and colder than jet C, which makes it highly inefficient with respect to jet C in terms of energy transfer to the ambient medium. This shows the relevance of the relativistic nature of outflows in the feedback process.
4.3. Radiative processes in the magnetized jet
As discussed, for instance, in Rieger et al. (2007) and BoschRamon & Rieger (2012), particle acceleration is expected to be more efficient in strong shocks in microquasar jets, with electrons and positrons typically cooling much more efficiently than protons and nuclei (see, e.g., BoschRamon & Khangulyan 2009). In these systems, strong magnetic field dissipation through magnetic reconnection could also be relevant for particle acceleration (BoschRamon 2012; BoschRamon & Rieger 2012).
With respect to hydrodynamical shocks, the strong recollimation shock(s) in relativistic jets can be a perfect candidate for generating a significant nonthermal population of particles. The nonthermal emission related to the jetwind interaction has been discussed in previous works (see, e.g., Molina et al. 2019, and references therein), although in all but one case (i.e., de la Cita et al. 2017, where the one clumpjet interaction was simulated using relativistic hydrodynamics), the jet dynamics was semianalytical. Nevertheless, in these works, the magnetic field was not informed using devoted RMHD simulations, but assuming either a fixed ratio of the magnetic to internal energy density, or establishing frozenin conditions in the hydrodynamical flow. Therefore, in this context, the simulations presented here provide with new insights into the nonthermal processes of HMXB jets.
The most straightforward conclusion that we can derive from the results of this work is that synchrotron emission and inverse Compton (IC) may be anticorrelated spacewise in the jet. Nonthermal electrons in highly magnetized regions may emit more synchrotron and less IC radiation, whereas in regions that are not strongly magnetized, it may be the opposite (see, e.g., Khangulyan & Aharonian 2005, for subtle effects of synchrotron and IC losses on the radiation spectra). All this applies in particular to electrons emitting synchrotron Xrays or very high energy IC photons, as they are more sensitive to radiative cooling versus advection in the jet.
On the other hand, our results for jet C show more shocks than in jet B. In these shocks, kinetic energy turns into internal energy, favoring plausible particle acceleration within the scales of the binary. Magnetic energy is also transferred to matter rather efficiently through magnetic acceleration at favorable expansions. Magnetic acceleration introduces greater complexity to radiation coming from the region where this occurs through a complex pattern of Doppler boosting. Another nonthermal process linked to the magnetic field, which is difficult to capture here but might be relevant in highly magnetized jets as well, is magnetic reconnection. Magnetic energy can be injected into particles by reconnecting magnetic lines, which can trigger (multiple) shocks in the region in which this occurs. These nonthermal processes and others that might be relevant (e.g., acceleration driven by turbulence or shear motion or gammaray reprocessing in the complex jetwind environment), are very difficult to quantify unless specific detailed studies are carried out, which is beyond the scope of this work.
As shown, jets might be more stable if they are moderately magnetized. If this were the case, they could reduce the amount of energy transferred to nonthermal particles within the binary system. However, as suggested in BoschRamon & Barkov (2016) and shown by Barkov & BoschRamon (2022), the impact of orbital motion n longer timescales than simulated here can further destabilize the jet, and so most of the energy dissipation may not take place within the system, but farther out, where orbital motion strongly shapes the jet evolution (see Molina & BoschRamon 2018, for an example of such a situation).
If the jet kinetic and magnetic powers are roughly similar, or if the latter dominates the former, the jet may become severely disrupted on the scales of the binary. This may seem incompatible with milliarcsecond observations of the jet of Cygnus X1 (Stirling et al. 2001), which show that the jet is still rather collimated on scales of tens of orbital separation distances. However, the combined effects of the stellar wind and orbital motion may cause the binary system to act as a nozzle for the highly magnetized, disordered jet, allowing the jet to reaccelerate and potentially recollimate on relatively short timescales given the fast drop in density and pressure outside this region. In addition, the magnetic field of the jet may act as a protective shield against mixing with the stellar wind, which would allow the jet to again reach relativistic speeds if the environment and/or the magnetic field were appropriate (e.g., Komissarov 2011). No such effect is expected from a hydrodynamical jet, which could hardly become relativistic again due to heavy wind entrainment (BoschRamon & Barkov 2016; Barkov & BoschRamon 2022).
Depending on the magnetic field strength, the jet will present different emission properties. Specifying these differences for comparison with observations requires specific simulations covering regions that are larger than the orbital separation distance. Simulations of quasisteady state jets would be also needed, as they are complementary to the present ones because the timescales involved in observations are often long enough to capture this jet regime. Moreover, the magnetic flux in the jet may be different depending on the particular accretionejection state of HMXBs. This different magnetization in the ejections, which would be continuous or discrete depending on the state of the source, would be a major factor (in addition to the flow energetics) determining the properties of the nonthermal emission. This possibility requires a detailed analysis of the observed outflow properties and it is left for future work.
Nevertheless, recent largescale 3D hydrodynamical simulations (Charlet et al. 2022) show that for two fiducial HMXBs (namely, Cygnus X1 and Cygnus X3), thermal cooling (i.e., freefree losses) could dominate the emission compared to nonthermal mechanisms in the beam and inner cocoon. In particular, the authors analyzed the jet outbreak and early propagation beyond the scale of the binary and concluded that these radiative losses do not play a dynamical or structural role for the space of parameters in Cygnus X1, although they can be relevant for the evolution of Cygnus X3 due to the denser stellar winds of the WolfRayet primary. Thus, according to our models and the timescales of our simulations, we do not expect that (non)thermal cooling produces a significant effect on the relativistic jet dynamics, but the characterization of the system emission would require a dedicated analysis. This is beyond the scope of this paper.
5. Conclusions
We have performed three different numerical simulations of microquasar jets propagating through an inhomogeneous stellar wind using the new code Lóstrego v1.0 for computations in relativistic astrophysics. Our numerical setup was based on the space of parameters described in previous 3D hydrodynamical jetwind simulations, including a toroidal magnetic field in transversal equilibrium with the gas to analyze its role for jet dynamics and longterm stability. We present the jet evolution throughout the initial minutes after the trigger of jet formation, which allowed us to keep a steady injection point in the grid and neglect the orbital motion of the binary.
Our simulations show that the magnetic field could play a stabilizing role in jet evolution as long as the flux of magnetic energy is lower than the kinetic energy flux. This stabilizing effect provides the simulated, magnetized jets, with additional collimation along the grids, as compared to similar, nonmagnetized jets. In contrast, a nonnegligible magnetic energy flux translates into a destabilization of the flow, with the consequence of increased energy dissipation and a slower jet head propagation velocity. The study of the threshold at which the magnetic field changes its role should be addressed in future publications.
The evaluation of the probability of jet propagation beyond the binary region has thus to be revisited in the presence of magnetic fields. In Perucho & BoschRamon (2008); Perucho et al. (2010a) it was concluded that this depended basically on the jet power, whereas in this work, we see that lowpower jets threaded by relatively weak fields could propagate longer distances without being destabilized. However, more powerful jets, in contrast, might be disrupted by currentdriven instabilities if the magnetic field is still dynamically relevant at the scales of the binary. Therefore, we conclude that a more thorough characterization of the magnetic field role in microquasar jets at these scales is needed in order to set a proper limit on the minimum energy flux to allow for propagation beyond the progenitor system.
In RMHD jet simulations, the presence of recollimation shocks within the binary region is more probable than in the case of RHD jets due to the role of magnetic tension, which acts against expansion and contributes to jet collimation (Martí et al. 2016; Fuentes et al. 2018), independently of drops in the ambient pressure (in contrast to RHD flows, Perucho & BoschRamon 2008; Perucho et al. 2010a), but the pressure jumps at these shocks and its consequent radiative properties still depend on the obliqueness with respect to the jet stream lines.
A detailed analysis of the energy distribution in the jet shows that magnetic and internal energies can both contribute very efficiently to jet acceleration in expansion regions. At shock waves (recollimation, standing, and terminal shocks), both magnetic and internal energies may grow at the expense of the kinetic energy budget. In the turbulent backflow, these magnetically dominated regions can convert magnetic energy into kinetic and internal energy of the flow (the latter by magnetic diffusion), while if the kinetic energy is dominant, this can be used as a reservoir to locally reinforce the magnetic field and to increase the internal energy. Moreover, we have also contributed further evidence to the relevant role of the relativistic nature of outflows in the study of jet energy transfer to the ambient medium, as predicted before by Perucho et al. (2017) for AGNs, and as analyzed here in the context of microquasar jet simulations.
Acknowledgments
The project that gave rise to these results received the support of a fellowship from “la Caixa” Foundation (ID 100010434). The fellowship code is LCF/BQ/DR19/11740030. J.L.M. acknowledges additional support from the Spanish Ministerio de Ciencia through grant PID2019105510GBC31. M.P. and J.MM. acknowledge support by the Spanish Ministerio de Ciencia through grants PID2019107427GBC33 and PGC2018095984BI00, and from the Generalitat Valenciana through grant PROMETEU/2019/071. M.P. acknowledges additional support from the Spanish Ministerio de Ciencia through grant PID2019105510GBC31. V.B.R. acknowledges financial support from the State Agency for Research of the Spanish Ministry of Science and Innovation under grant PID2019105510GBC31 and through the ”Unit of Excellence María de Maeztu 20202023” award to the Institute of Cosmos Sciences (CEX2019000918M), and by the Catalan DEC grant 2017 SGR 643. V.BR. is Correspondent Researcher of CONICET, Argentina, at the IAR. Computer simulations have been carried out in the Servei d’Informàtica de la Universitat de València (Tirant supercomputer). We thank the referee for all the constructive comments and suggestions that helped to improve the quality of the manuscript.
References
 Aharonian, F., Akhperjanian, A. G., Aye, K. M., et al. 2005, Science, 309, 746 [NASA ADS] [CrossRef] [Google Scholar]
 Albert, J., Aliu, E., Anderhub, H., et al. 2006, Science, 312, 1771 [NASA ADS] [CrossRef] [Google Scholar]
 Albert, J., Aliu, E., Anderhub, H., et al. 2007, ApJ, 665, L51 [NASA ADS] [CrossRef] [Google Scholar]
 Anile, A. M. 1989, Relativistic Fluids and Magnetofluids: with Applications in Astrophysics and Plasma Physics (Cambridge; New York: Cambridge University Press) [Google Scholar]
 Antón, L., Miralles, J. A., Martí, J. M., et al. 2010, ApJS, 188, 1 [Google Scholar]
 Balsara, D. 2001, ApJS, 132, 83 [NASA ADS] [CrossRef] [Google Scholar]
 Balsara, D. S., & Spicer, D. S. 1999, J. Comput. Phys., 149, 270 [NASA ADS] [CrossRef] [Google Scholar]
 Barkov, M. V., & BoschRamon, V. 2022, MNRAS, 510, 3479 [NASA ADS] [CrossRef] [Google Scholar]
 Beckwith, K., & Stone, J. M. 2011, ApJS, 193, 6 [NASA ADS] [CrossRef] [Google Scholar]
 Belloni, T. M., & Motta, S. E. 2016, in Transient Black Hole Binaries, ed. C. Bambi, 440, 61 [Google Scholar]
 Blandford, R. D., & Payne, D. G. 1982, MNRAS, 199, 883 [CrossRef] [Google Scholar]
 Blandford, R. D., & Znajek, R. L. 1977, MNRAS, 179, 433 [NASA ADS] [CrossRef] [Google Scholar]
 Bodo, G., Massaglia, S., Ferrari, A., & Trussoni, E. 1994, A&A, 283, 655 [NASA ADS] [Google Scholar]
 Bordas, P., BoschRamon, V., Paredes, J. M., & Perucho, M. 2009, A&A, 497, 325 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 BoschRamon, V. 2012, A&A, 542, A125 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 BoschRamon, V., & Barkov, M. V. 2016, A&A, 590, A119 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 BoschRamon, V., & Khangulyan, D. 2009, Int. J. Mod. Phys. D, 18, 347 [Google Scholar]
 BoschRamon, V., & Rieger, F. M. 2012, in Astroparticle, Particle, Space Physics and Detectors For Physics Applications, Proceedings of the 13th ICATPP Conference, eds. G. Simone (World Scientific Publishing Co. Pte. Ltd.), 219 [CrossRef] [Google Scholar]
 BoschRamon, V., Perucho, M., & Bordas, P. 2011, A&A, 528, A89 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Brackbill, J. U., & Barnes, D. C. 1980, J. Comput. Phys., 35, 426 [NASA ADS] [CrossRef] [Google Scholar]
 Brio, M., & Wu, C. C. 1988, J. Comput. Phys., 75, 400 [NASA ADS] [CrossRef] [Google Scholar]
 Cassinelli, J. P., Ignace, R., Waldron, W. L., et al. 2008, ApJ, 683, 1052 [NASA ADS] [CrossRef] [Google Scholar]
 Castro, M. J., Gallardo, J. M., & Marquina, A. 2017, Comput. Phys. Commun., 219, 108 [NASA ADS] [CrossRef] [Google Scholar]
 Charlet, A., Walder, R., Marcowith, A., et al. 2022, A&A, 658, A100 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 CruzOsorio, A., Fromm, C. M., Mizuno, Y., et al. 2022, Nat. Astron., 6, 103 [NASA ADS] [CrossRef] [Google Scholar]
 de la Cita, V. M., del Palacio, S., BoschRamon, V., et al. 2017, A&A, 604, A39 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Dedner, A., Kemm, F., Kröner, D., et al. 2002, J. Comput. Phys., 175, 645 [Google Scholar]
 Del Zanna, L., & Bucciantini, N. 2002, A&A, 390, 1177 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Del Zanna, L., Bucciantini, N., & Londrillo, P. 2003, A&A, 400, 397 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Dolezal, A., & Wong, S. S. M. 1995, J. Comput. Phys., 120, 266 [NASA ADS] [CrossRef] [Google Scholar]
 Dubal, M. R. 1991, Comput. Phys. Comm., 64, 221 [NASA ADS] [CrossRef] [Google Scholar]
 Evans, C. R., & Hawley, J. F. 1988, ApJ, 332, 659 [Google Scholar]
 Event Horizon Telescope Collaboration (Akiyama, K., et al.) 2019, ApJ, 875, L1 [Google Scholar]
 Fender, R. P., Belloni, T. M., & Gallo, E. 2004, MNRAS, 355, 1105 [NASA ADS] [CrossRef] [Google Scholar]
 Fuentes, A., Gómez, J. L., Martí, J. M., & Perucho, M. 2018, ApJ, 860, 121 [Google Scholar]
 Gardiner, T. A., & Stone, J. M. 2005, J. Comput. Phys., 205, 509 [NASA ADS] [CrossRef] [Google Scholar]
 Ghisellini, G., Tavecchio, F., Maraschi, L., Celotti, A., & Sbarrato, T. 2014, Nature, 515, 376 [NASA ADS] [CrossRef] [Google Scholar]
 Giacomazzo, B., & Rezzolla, L. 2006, J. Fluid Mech., 562, 223 [NASA ADS] [CrossRef] [Google Scholar]
 Giacomazzo, B., & Rezzolla, L. 2007, Class. Quant. Grav., 24, S235 [NASA ADS] [CrossRef] [Google Scholar]
 Guan, X., Li, H., & Li, S. 2014, ApJ, 781, 48 [NASA ADS] [CrossRef] [Google Scholar]
 Harten, A., Lax, P., & van Leer, B. 1983, SIAM Rv., 25, 35 [CrossRef] [Google Scholar]
 Janssen, M., Falcke, H., Kadler, M., et al. 2021, Nat. Astron., 5, 1017 [NASA ADS] [CrossRef] [Google Scholar]
 Keppens, R., Meliani, Z., van der Holst, B., & Casse, F. 2008, A&A, 486, 663 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Khangulyan, D., & Aharonian, F. 2005, in High Energy GammaRay Astronomy, eds. F. A. Aharonian, H. J. Völk, & D. Horns, Am. Inst. Phys. Conf. Ser., 475, 359 [NASA ADS] [CrossRef] [Google Scholar]
 Kim, J.Y., Krichbaum, T. P., & Broderick, A. E. 2020, A&A, 640, A69 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Komissarov, S. S. 1999a, MNRAS, 308, 1069 [Google Scholar]
 Komissarov, S. S. 1999b, MNRAS, 303, 343 [NASA ADS] [CrossRef] [Google Scholar]
 Komissarov, S. S. 2011, Mem. Soc. Astron. It., 82, 95 [NASA ADS] [Google Scholar]
 Komissarov, S. S., Barkov, M. V., Vlahakis, N., & Königl, A. 2007, MNRAS, 380, 51 [Google Scholar]
 Komissarov, S. S., Porth, O., & Lyutikov, M. 2015, Comput. Astrophys. Cosmol., 2, 9 [Google Scholar]
 Krtička, J. 2014, A&A, 564, A70 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Kylafis, N. D., Contopoulos, I., Kazanas, D., & Christodoulou, D. M. 2012, A&A, 538, A5 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Leismann, T., Antón, L., Aloy, M. A., et al. 2005, A&A, 436, 503 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Lewis, G. M., & Austin, P. H. 2002, in 11th Conference on Atmospheric Radiation,American Meteorological Society Conference Series, eds. G. Smith, & J. Brodie, 123 [Google Scholar]
 Lind, K. R., Payne, D. G., Meier, D. L., & Blandford, R. D. 1989, ApJ, 344, 89 [Google Scholar]
 Londrillo, P., & Del Zanna, L. 2000, ApJ, 530, 508 [NASA ADS] [CrossRef] [Google Scholar]
 Marino, A., Malzac, J., Del Santo, M., et al. 2020, MNRAS, 498, 3351 [CrossRef] [Google Scholar]
 Martí, J.M. 2015a, MNRAS, 452, 3106 [Google Scholar]
 Martí, J.M. 2015b, Comput. Phys. Comm., 191, 100 [CrossRef] [Google Scholar]
 Martí, J.M. 2019, Galaxies, 7, 24 [Google Scholar]
 Martí, J. M., & Müller, E. 2015, Liv. Rev. Comput. Astrophys., 1, 3 [Google Scholar]
 Martí, J. M., Müller, E., Font, J. A., Ibáñez, J. M. Z., & Marquina, A. 1997, ApJ, 479, 151 [Google Scholar]
 Martí, J. M., Perucho, M., & Gómez, J. L. 2016, ApJ, 831, 163 [Google Scholar]
 McKinney, J. C., & Blandford, R. D. 2009, MNRAS, 394, L126 [Google Scholar]
 McKinney, J. C., Tchekhovskoy, A., & Blandford, R. D. 2012, MNRAS, 423, 3083 [Google Scholar]
 McNamara, B. R., & Nulsen, P. E. J. 2007, ARA&A, 45, 117 [NASA ADS] [CrossRef] [Google Scholar]
 Migliari, S., & Fender, R. P. 2006, MNRAS, 366, 79 [CrossRef] [Google Scholar]
 Migliari, S., Fender, R. P., Rupen, M., et al. 2003, MNRAS, 342, L67 [CrossRef] [Google Scholar]
 Migliari, S., Fender, R. P., Rupen, M., et al. 2004, Nucl. Phys. B Proc. Suppl., 132, 628 [NASA ADS] [CrossRef] [Google Scholar]
 Mignone, A., & Bodo, G. 2006, MNRAS, 368, 1040 [Google Scholar]
 Mignone, A., & Del Zanna, L. 2021, J. Comput. Phys., 424, 109748 [NASA ADS] [CrossRef] [Google Scholar]
 Mignone, A., & McKinney, J. C. 2007, MNRAS, 378, 1118 [Google Scholar]
 Mignone, A., Massaglia, S., & Bodo, G. 2005, Space Sci. Rev., 121, 21 [NASA ADS] [CrossRef] [Google Scholar]
 Mignone, A., Ugliano, M., & Bodo, G. 2009, MNRAS, 393, 1141 [NASA ADS] [CrossRef] [Google Scholar]
 Mignone, A., Rossi, P., Bodo, G., Ferrari, A., & Massaglia, S. 2010, MNRAS, 402, 7 [NASA ADS] [CrossRef] [Google Scholar]
 Mignone, A., Zanni, C., Tzeferacos, P., et al. 2012, ApJS, 198, 7 [Google Scholar]
 Mirabel, I. F., & Rodríguez, L. F. 1999, ARA&A, 37, 409 [Google Scholar]
 Mizuno, Y., Hardee, P., & Nishikawa, K.I. 2007, ApJ, 662, 835 [Google Scholar]
 Moffat, A. F. J. 2008, in Clumping in HotStar Winds, eds. W. R. Hamann, A. Feldmeier, & L. M. Oskinova, 17 [Google Scholar]
 Molina, E., & BoschRamon, V. 2018, A&A, 618, A146 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Molina, E., del Palacio, S., & BoschRamon, V. 2019, A&A, 629, A129 [EDP Sciences] [Google Scholar]
 MonceauBaroux, R., Porth, O., Meliani, Z., & Keppens, R. 2014, A&A, 561, A30 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 MoyaTorregrosa, I., Fuentes, A., Martí, J. M., Gómez, J. L., & Perucho, M. 2021, A&A, 650, A60 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Muijres, L. E., Vink, J. S., de Koter, A., Müller, P. E., & Langer, N. 2012, A&A, 537, A37 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Noh, W. 1987, J. Comput. Phys., 72, 78 [NASA ADS] [CrossRef] [Google Scholar]
 Orszag, S. A., & Tang, C. M. 1979, J. Fluid Mech., 90, 129 [NASA ADS] [CrossRef] [Google Scholar]
 Paredes, J. M., Martí, J., Ribó, M., & Massi, M. 2000, Science, 288, 2340 [Google Scholar]
 Perucho, M. 2019, Galaxies, 7, 70 [Google Scholar]
 Perucho, M., & BoschRamon, V. 2008, A&A, 482, 917 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Perucho, M., & BoschRamon, V. 2012, A&A, 539, A57 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Perucho, M., Martí, J. M., & Hanasz, M. 2005, A&A, 443, 863 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Perucho, M., BoschRamon, V., & Khangulyan, D. 2010a, A&A, 512, L4 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Perucho, M., Martí, J. M., Cela, J. M., et al. 2010b, A&A, 519, A41 [CrossRef] [EDP Sciences] [Google Scholar]
 Perucho, M., Quilis, V., & Martí, J.M. 2011, ApJ, 743, 42 [NASA ADS] [CrossRef] [Google Scholar]
 Perucho, M., Martí, J.M., Quilis, V., & Ricciardelli, E. 2014, MNRAS, 445, 1462 [NASA ADS] [CrossRef] [Google Scholar]
 Perucho, M., Martí, J.M., Quilis, V., & BorjaLloret, M. 2017, MNRAS, 471, L120 [NASA ADS] [CrossRef] [Google Scholar]
 Perucho, M., Martí, J.M., & Quilis, V. 2019, MNRAS, 482, 3718 [NASA ADS] [CrossRef] [Google Scholar]
 Porth, O. 2013, MNRAS, 429, 2482 [NASA ADS] [CrossRef] [Google Scholar]
 Porth, O., Komissarov, S. S., & Keppens, R. 2014, MNRAS, 438, 278 [Google Scholar]
 Powell, K. G. 1994, Approximate Riemann Solver for Magnetohydrodynamics (that works in more than one dimension) [Google Scholar]
 Quirk, J. J. 1994, Int. J. Numer. Methods Fluids, 18, 555 [NASA ADS] [CrossRef] [Google Scholar]
 Rieger, F. M., BoschRamon, V., & Duffy, P. 2007, Ap&SS, 309, 119 [NASA ADS] [CrossRef] [Google Scholar]
 Roe, P. L. 1986, Ann. Rev. Fluid Mech., 18, 337 [Google Scholar]
 Romero, G. E., Torres, D. F., Kaufman Bernadó, M. M., & Mirabel, I. F. 2003, A&A, 410, L1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Runacres, M. C., & Owocki, S. P. 2002, A&A, 381, 1015 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Ryu, D., Miniati, F., Jones, T. W., & Frank, A. 1998, ApJ, 509, 244 [NASA ADS] [CrossRef] [Google Scholar]
 Shu, C.W., & Osher, S. 1989, J. Comput. Phys., 83, 32 [NASA ADS] [CrossRef] [Google Scholar]
 Stirling, A. M., Spencer, R. E., de la Force, C. J., et al. 2001, MNRAS, 327, 1273 [Google Scholar]
 Suresh, A., & Huynh, H. T. 1997, J. Comput. Phys., 136, 83 [Google Scholar]
 Tavani, M., Kniffen, D., Mattox, J. R., Paredes, J. M., & Foster, R. S. 1998, ApJ, 497, L89 [Google Scholar]
 Tchekhovskoy, A., Narayan, R., & McKinney, J. C. 2011, MNRAS, 418, L79 [NASA ADS] [CrossRef] [Google Scholar]
 Tóth, G. 2000, J. Comput. Phys., 161, 605 [Google Scholar]
 van Leer, B. 1974, J. Comput. Phys., 14, 361 [NASA ADS] [CrossRef] [Google Scholar]
 van Leer, B. 1977, J. Comput. Phys., 23, 276 [Google Scholar]
 van Putten, M. H. P. M. 1993, J. Comput. Phys., 105, 339 [NASA ADS] [CrossRef] [Google Scholar]
 Vilhu, O., Kallman, T. R., Koljonen, K. I. I., & Hannikainen, D. C. 2021, A&A, 649, A176 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Vlahakis, N. 2004, Ap&SS, 293, 67 [NASA ADS] [CrossRef] [Google Scholar]
 Vlahakis, N., & Königl, A. 2004, ApJ, 605, 656 [Google Scholar]
 Wagner, A. Y., & Bicknell, G. V. 2011, ApJ, 728, 29 [NASA ADS] [CrossRef] [Google Scholar]
 Wang, P., Abel, T., & Zhang, W. 2008, ApJS, 176, 467 [NASA ADS] [CrossRef] [Google Scholar]
 Wright, A. J., & Hawke, I. 2019, ApJS, 240, 8 [NASA ADS] [CrossRef] [Google Scholar]
 Yoon, D., & Heinz, S. 2015, ApJ, 801, 55 [NASA ADS] [CrossRef] [Google Scholar]
 Yoon, D., Morsony, B., Heinz, S., et al. 2011, ApJ, 742, 25 [NASA ADS] [CrossRef] [Google Scholar]
 Yoon, D., Zdziarski, A. A., & Heinz, S. 2016, MNRAS, 456, 3638 [NASA ADS] [CrossRef] [Google Scholar]
 Zanin, R., FernándezBarral, A., de Oña Wilhelmi, E., et al. 2016, A&A, 596, A55 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Appendix A: Lóstrego v1.0. A new tool for computational relativistic astrophysics
In this appendix, we present Lóstrego v1.0, a new computational tool for simulating astrophysical relativistic plasmas in 1D, 2D, and 3D Cartesian coordinates. The code is entirely written in the FORTRAN programming language and solves the conservative equations of special RMHD in finite volumes. It is especially designed to run in multiple cores to exploit the capacity of modern parallel architectures. Lóstrego v1.0 is based on the high resolution shockcapturing (HRSC) methods that have been probed to be robust and accurate in other similar codes (see, e.g., Martí & Müller 2015, and references therein), extending the techniques implemented in Martí (2015a) to three dimensions. Essentially, the algorithm follows the reconstructsolveaverage (RSA) strategy: first, cellaverage primitive variables are reconstructed to the cell interfaces, where a Riemann problem is then solved numerically by any approximate Riemann solver. Once fluxes are known at each cell face, the conserved variables are evolved explicitly in time using total variation diminishing (TVD) RungeKutta algorithms. Finally, an inversion scheme is applied to recover the primitive variables from the timeevolved solution. Several onedimensional and multidimensional tests are presented in this appendix to probe the ability of the new code to solve classical problems in the field of RMHD.
A.1. Basic equations
The RMHD system of partial differential equations in the Minkowski metric^{3} and Cartesian coordinates can be written as a system of conservation laws (Komissarov 1999b),
where U = {D, S^{j}, τ, B^{j}} is a vector of conserved variables, and F^{i} are the vector of fluxes for each spatial direction. Both the vector of conserved variables U and the vector of fluxes F^{i} can be expressed in terms of a set of primitive variables V = {ρ, v^{j}, p, B^{j}} using the following relations:
where D is the relativistic rest mass density, S^{j} is the momentum density of the magnetized fluid, τ is the energy density measured in the Eulerian frame, and δ^{ij} is the Kronecker delta. The set of primitive variables are the fluid restmass density ρ, fluid threevelocity v_{j}, and gas pressure p. The magnetic field vector B^{j} is at the same time primitive and conserved variable, so the conversion from one set to the other is trivial. In Eqs. A.2 and A.3 we introduced the hydromagnetic specific enthalpy h^{*},
where ϵ is the specific internal energy. The total pressure p^{*} of the magnetic fluid is given by
In terms of the vectors v and B in the laboratory frame, the fourvelocity u^{α} and magnetic field b^{α} covariant vectors are
with the normalizations
where is the Lorentz factor, and summation over repeated indices is assumed. The time evolution of the system of equations of classical and relativistic MHD must also preserve the magnetic field divergencefree constraint,
Finally, the system of RMHD equations is closed with an equation of state (EoS) of the form
which relates the thermodynamic primitive variables with each other. For the particular case of ideal gases, this relation becomes
where Γ is the adiabatic index. This index should be set as Γ = 5/3 for nonrelativistic or mildly relativistic flows and as Γ = 4/3 in the ultrarelativistic limit.
A.2. Numerical methods
A.2.1. Spatial and temporal discretization
We have implemented multidimensional HRSC methods following the finite volumes (FV) strategy and dimensional splitting, taking the integral form of the conservation laws (Eq. A.1) and cellaveraged values. This means that operators that involve spatial derivatives are applied dimension by dimension, permuting the spatial coordinates cyclically. Eq. A.1 admits a conservative semidiscretization in a rectangular grid of size L_{x} × L_{y} × L_{z} and resolution N_{x} × N_{y} × N_{z},
where {i, j, k} is the position of each cell center, and {i ± 1/2, j ± 1/2, k ± 1/2} is the position of the right and left cell interface, respectively. The elements Δx = L_{x}/N_{x}, Δy = L_{y}/N_{y}, Δz = L_{z}/N_{z} are the cell sizes for each spatial dimension. As mentioned before, in the FV formalism, the conserved variables represent an approximation to the average over the cell volume,
while the numerical fluxes represent an approximation to the average in the surface of the cell,
A.2.2. Spatial cell reconstruction
To advance the set of conserved variables in time, numerical fluxes must be known at cell interfaces, the latter being functions of the primitive variables. These are reconstructed from cellcentered values by means of interpolation routines that are designed to reconstruct a piecewise polynomial approximation inside each cell preserving the monotonicity and TVstability properties of the algorithm. Consistently with our dimensional splitting, reconstruction is performed in one dimension following each spatial coordinate. The current version of Lóstrego admits either Godunov (first order) or linear (second order) spatial reconstruction. The piecewise linear method (PLM) in the xdirection follows
where is the linear slope at cell i. Similar expressions hold for the y, z spatial directions. The most popular slope limiters are MINMOD (Roe 1986), MC (van Leer 1977), and VAN LEER (van Leer 1974). These limiters avoid the generation of spurious extrema at cell interfaces and reduce the accuracy of the method to first order at extrema with a vanishing slope. Although conserved variables can be directly reconstructed to cell interfaces, we interpolate the primitive variables since experience has shown it to be more robust than the former approach. If the PLM reconstruction leads to nonphysical situations in which v^{2} > 1, we reduce the method to first order (i.e., Godunov reconstruction). Highorder schemes like the MP reconstruction (Suresh & Huynh 1997) are already available in the code, and we plan to test them in future publications.
A.2.3. Numerical fluxes and the Riemann problem
When the cellcentered primitive variables are spatially reconstructed to obtain (i.e., the values of V at the left and right sides of the interface i + 1/2), computation of numerical fluxes requires solving a Riemann problem at each zone edge of the type
Although exact Riemann solvers are available in the literature for RMHD (Giacomazzo & Rezzolla 2006), computation of numerical fluxes usually involves upwind strategies based on approximate solutions. In the current version of Lóstrego, three different approximate solvers are available, all of the HartenLaxvan Leer (HLL) family: HLL (Harten et al. 1983), HLLC (Mignone & Bodo 2006), and HLLD (Mignone et al. 2009). In this collection of nonlinear solvers, the solution is approximated by N < 7 waves that travel with characteristic speeds λ_{k + 1} > λ_{k}, k = 1, ..., N − 1, separated by N + 1 intermediate states. The characteristic speed of the outermost waves in the Riemann fan (i.e., fast magnetosonic waves), which are the maximum and minimum eigenvalues of the Jacobian matrix of the left and right states, require finding the roots of a quartic equation (see, e.g., Anile 1989; Antón et al. 2010), which must be solved by any rootfinding numerical method. As the number of waves and intermediate states increases, the diffusion of the method is reduced, and thus we obtain more accurate representations of the real solution. Figure A.1 shows a typical Riemann fan representing all different states in the three approximate solvers contained in Lóstrego. In HLL (N = 2), which is the simplest method of the family, the initial discontinuity is decomposed into two fast magnetosonic waves with characteristic speeds λ_{L}, λ_{R}, such that the internal flux in the fan can be derived from the RankineHugoniot jump conditions across the magnetosonic waves. This method is robust and attractive for RMHD applications because it preserves positivity of pressure and v< 1, but it may suffer from high levels of numerical diffusion compared to more sophisticated solvers such as HLLC or HLLD. HLLC (N = 3) involves two internal states separated by a contact wave. This solver increases the accuracy of HLL significantly, but it requires a separate treatment of the hydrodynamical limit (B = 0) and may show different pathologies when B → 0, especially in 3D applications. For this reason, we have also included the fivewave Riemann solver HLLD (N = 5), where in addition to the contact mode and the two fast magnetosonic waves, it also considers two rotational discontinuities. In this case, the solution of the system requires applying the RankineHugoniot jump conditions across all the waves in the fan of Fig. A.1. Although this solver is more complex and thus more expensive from the point of view of computational efficiency, we find it very robust and accurate for multidimensional applications. However, less dissipative solvers such as HLLC and HLLD might develop undesirable artifacts (see, e.g., Quirk 1994; Wang et al. 2008) that require a careful treatment, for example, degrading the solver to HLL in the proximity of strong shocks.
Fig. A.1. Riemann fan for the HLLfamily of Riemann solvers: HLL (blue), HLLC (blue+black), and HLLD (blue+black+red). 
A.2.4. Time integration
When numerical fluxes are resolved at cell interfaces, conserved variables can be advanced one time step following Eq. A.14 and predictorcorrector methods. Numerical time integration is based on second (RK2) and thirdorder (RK3) total variation dimishing (TVD) RungeKutta time integrators described by Shu & Osher (1989). The secondorder method (RK2) follows
while the thirdorder method (RK3) is given by
where L(U^{(n)}) is the upwind differencing operator,
A.2.5. Divergencefree constraint
In general, HRSC methods do not preserve the divergencefree condition given by Eq. A.11. This can produce numerical artificial forces that can lead to incorrect solutions and the eventual failure of the code. Thus, different approaches exist to control the solenoidal condition in RMHD codes. The most popular approaches are (see, e.g., Tóth 2000; Martí & Müller 2015, and references therein) the eightwave method (Powell 1994), the hyperbolicparabolic divergence cleaning (Dedner et al. 2002), the projection scheme (Brackbill & Barnes 1980), and the constrained transport method (Evans & Hawley 1988; Ryu et al. 1998; Balsara & Spicer 1999, CT). The latter has been probed to be robust and accurate because it guarantees that the divergencefree constraint is fulfilled up to machine roundoff errors by definition, preserving the divergence of the initial setup during the simulation. Thus, it is the only scheme we have implemented in Lóstrego v1.0. The method was originally proposed by Evans & Hawley (1988) and adapted to Godunovtype Riemann solvers by Ryu et al. (1998), Balsara & Spicer (1999). The main characteristic of this approach is that it requires the introduction of a staggered representation of the magnetic field at cell interfaces, which is evolved according to a semidiscrete version of the induction equation,
where are the staggered field components, which represent an average of the magnetic field at cell surface,
The electromotive forces (emf) Ω^{i} are discrete representations of the components of Ω = v × B defined at cell corners. To compute them, we follow the fluxCT formalism (Balsara & Spicer 1999; Giacomazzo & Rezzolla 2007), taking the simple average of the neighboring upwind numerical fluxes,
where is the μcomponent of the flux along the νdirection.
Arithmetic averaging is the simplest and most often suggested CT scheme, but this approach may include insufficient dissipation and does not reduce to the planeparallel algorithm for gridalign flows (Gardiner & Stone 2005; Mignone & Del Zanna 2021). Thus, we modified the definitions of Eqs. A.30 to A.32 following the CTcontact scheme of Gardiner & Stone (2005), where the CT algorithm is presented as a spatial integration procedure. This scheme has been probed to be robust and accurate in practical applications, and it does reduce to the base upwind method for gridaligned planar flows. In the CTcontact approximation, the average electromotive forces are modified according to
where similar expressions hold for Ω^{x}, Ω^{y} by cyclic permutation of the spatial indices. Gardiner & Stone (2005) proposed several ways to compute the derivatives of Eq. A.33. Here, we follow the CT algorithm, where derivatives are selected in an upwind way according to the sign of the velocity of the contact mode,
To obtain the derivatives that appear in Eq. A.34, Gardiner & Stone (2005) proposed to subtract the face centered emf computed with the upwind fluxes and the emf evaluated at the cell center,
The cellcenter magnetic field can be recovered from the staggered representation using a simple linear interpolation at the end of the CT step,
The interpolations performed to obtain the required fluxes at cell edges in the CTcontact formalism, and to recover the cellcentered magnetic field from the staggered solution, limit the accuracy of the algorithm to second order. Other strategies such as the UCTHLL(D) method of Mignone & Del Zanna (2021) were followed to extend the CT method to higher than second order.
A.2.6. Recovery of primitive variables
Finally, in order to start a new time iteration, we must recover the primitive variables from the already evolved conserved solution. In RHD and RMHD codes, this task requires solving a highly nonlinear algebraic system of equations, which constitutes one of the most challenging tasks from the point of view of computational efficiency. The recovery of primitive variables is usually considered as one of the bottlenecks for the speedup of the code (Wright & Hawke 2019). Of the different strategies that exist in the literature (see, e.g, Martí & Müller 2015, and references therein), we chose to recover the set of primitive variables following the inversion scheme of Mignone & McKinney (2007). This algorithm can be extended to general EoS, and it has been probed to avoid numerical problems due to loss of precision in the nonrelativistic and ultrarelativistic limit. For example, this might be relevant for microquasar jet simulations, where we showed that other inversion schemes such as the one used by Leismann et al. (2005); Martí (2015a) failed to recover the low pressure of the nonrelativistic stellar wind with acceptable accuracy. The relevant equations of this inversion scheme are
where E′ = E − D, Z′ = Z − D, Z = DhW, S_{B} = S ⋅ B and u^{2} = W^{2}v^{2}. Equation A.39 can be solved for Z′ by means of a onedimensional NewtonRaphson iterative method using Eq. A.40 to express u^{2} as a function of Z′. The method involves the computation of the derivative
where dp/dZ′ depends on the EoS of the code. To avoid catastrophic losses of accuracy in the nonrelativistic and ultrarelativistic limits, Mignone & McKinney (2007) chose p = p(χ, ρ) with χ = ρϵ + p to obtain
For ideal gases, ∂p/∂χ and ∂p/∂ρ become
Moreover, the variable χ can be written as a function of Z′ as
and its derivative with respect to Z′ as
where
On the other hand,
When Z′ and p have been found with some desirable degree of accuracy, the inversion scheme is completed by recovering the threevelocity and the density of the gas through
We refer to Mignone & McKinney (2007) for further details of the full algorithm. If the inversion step fails to recover a proper physical solution (i.e., p < 0), we follow the same strategy as in the PLUTO code and fix the pressure to a positive threshold before solving for v^{2} using a bisections root finder with the following function:
Using this approach, Z has to be recomputed after each iteration using Eqs. A.39 to A.46. Then, the total energy needs to be corrected for according to the new solution.
A.2.7. Correction of conserved variables
The RMHD codes based on constrained transport algorithms are subject to several pathologies that may lead to nonphysical solutions, especially for highly magnetized flows under the influence of strong shocks. The decoupled evolution of the staggered and cellcentered magnetic field makes the algorithm inconsistent because the conserved variables have been evolved according to the fluxes computed from the cellcentered field, while the staggered solution has been computed with the CT method. This inconsistency relies on the basis to develop correction algorithms of the conserved variables, which may be especially important when the magnetic pressure dominates the gas pressure by more than two orders of magnitude (Martí 2015b). Lóstrego v1.0 admits both the nonrelativistic energy correction of Mignone & Bodo (2006)
and the fully relativistic correction of both energy and momentum of Martí (2015b). In this last approach, the nonrelativistic correction of the energy, now , is used to obtain a first approximation of the primitive variables and then use the new flow velocity v^{(1)} to complete the relativistic correction of the momentum and energy,
Although these algorithms can be applied before or after recovering the primitive variables, in Lóstrego, we decided to perform the correction before. Following the notation of Martí (2015b), these algorithms are called CA1 (Eq. A.51) and CA2 (Eq. A.52 and Eq. A.53), respectively.
A.3. Parallelization
Lóstrego v1.0 is parallelized with an hybrid scheme using OpenMP (OMP) and MPI library instructions in order to exploit the computational power of distributed and shared memory. This parallelization configuration is based on the parallel architecture of the code Ratpenat (Perucho et al. 2010b). The hybrid scheme is only available in the code for fully 3D simulations, while (2D) computations like the tests described in the next section might benefit from shared memory parallelization with OMP directives. For the large 3D simulations for which the hybrid scheme is available, the computational box is initially decomposed into several subdomains such that each one is assigned to an MPI node that performs all the calculations independently using sharedmemory OMP threads. At the beginning of each time step, all MPI blocks that share internal boundaries must interchange a collection of ghost cells with the neighbors in the three spatial directions. These cells are used as boundary conditions for the blocks that do not limit the physical boundaries of the grid. The latter are imposed according to the physical conditions on the box boundaries, which are particular for each numerical simulation. Moreover, the implementation of the CTcontact algorithm (see Sec. A.2) requires one more MPI communication step before the time evolution of the staggered components of the magnetic field. This CT algorithm requires not only the directionwise boundaries in the three spatial directions, but also the corners of each subbox. Nevertheless, this information can be also transmitted directionwise because the corners are already available in the neighboring boundaries after the first MPI comunnication.
A.4. Testing benchmark
We provide an extensive collection of onedimensional and multidimensional numerical tests to probe the performance of the new Lóstrego v1.0 code. Unless otherwise stated, all tests shown in this section were performed with the HLLD Riemann solver and the secondorder piecewise linear method (PLM) with the VAN LEER slope limiter for cell reconstruction. A weak form of flattening is introduced that degrades the slope limiter to MINMOD whenever a strong shock is detected (Mignone & Bodo 2006), but no degradation of the HLLD Riemann solver is applied in the tests. We used a thirdorder TVDpreserving RungeKutta algorithm for time integration with CFL = 0.3. The magnetic field divergencefree constraint is preserved with the CT method, where electromotive forces were averaged according to the CTcontact formalism. The relativistic correction scheme CA2 was used to correct the conserved variables after each time integration.
A.4.1. One dimension
Sinusoidal perturbation As a first 1D application, we show the purely hydrodynamical problem initially formulated by Dolezal & Wong (1995) with a nuclear EoS and adapted to ideal gases in Del Zanna & Bucciantini (2002). A onedimensional unitary shock tube with a resolution of 2000 computational zones is perturbed in the right zone (0.5 < x < 1.0) with a density sinusoidal profile such that U_{L} = (ρ, v, p) = (5, 0, 50) and U_{R} = (ρ, v, p) = (2 + 0.3sin(50x), 0, 5). Figure A.2 shows the solution of the problem at t = 0.35. The original jump in pressure generates a blast wave that interacts with the density perturbation. Qualitatively, our results are similar to those presented in Del Zanna & Bucciantini (2002).
Fig. A.2. Gas density (top), pressure (middle), and xvelocity (bottom) at t = 0.35 for the hydrodynamical sinusoidal perturbation test. The shock tube contains 2000 grid points in the xdirection. 
RMHD shock tubes A collection of 1D RMHD test problems were initially proposed by Dubal (1991), van Putten (1993) and extended in Komissarov (1999b); Balsara (2001). These tests constitute a set of Riemann problems that have been established as a benchmark to test the robustness, accuracy, stability, and diffusion of any RMHD code. We consider the five shock tube problems of Balsara (2001) (hereinafter, BA15), although we only show the solution of two selected tests. The initial conditions of thesefive tests can be found in the aforementioned reference, so we do not reproduce them here. For all of these tests, we considered a onedimensional unit grid with a resolution of 1600 cells. The adiabatic coefficient was set to Γ = 5/3 for all tests except for BA1, where Γ = 2.0. The BA1 problem, which was originally proposed by Brio & Wu (1988) and extended to relativistic MHD by van Putten (1993), describes the formation of a leftgoing fast rarefaction wave, a leftgoing compound wave (that only appears in the numerical solution), a contact discontinuity, a rightgoing slow shock, and a rightgoing fast rarefaction wave in a moderate relativistic flow. The right part of the shock tube is magnetically dominated. The solution of this test is not shown in the paper. BA2 and BA3 tests are both blast waves problems, although the initial pressure jump is moderate for BA2 and strong for BA3. The BA2 problem, which is not shown in this paper either, is wellresolved with our code, and we clearly identify a leftgoing fast rarefaction wave, a leftgoing slow rarefaction, a contact discontinuity, a rightgoing slow shock, and a rightgoing fast shock. Figure A.3 shows the solution of the BA3 test at t = 0.4 (red dots) and the analytical solution provided by Giacomazzo & Rezzolla (2006) (blue line). The collection of waves that develops as a consequence of the initial pressure discontinuity is similar to those obtained in BA2, but due to the jump in pressure of several orders of magnitude, it develops a strong relativistic flow. As a direct consequence of the large Lorentz factors achieved, the contact discontinuity and the rightgoing shocks are underresolved, even with the less diffusive HLLD Riemann solver. This was not the case of BA2, where a weaker blast wave leads to lower Lorentz factors and thus to mildly relativistic flows. The solution of the BA4 test at t = 0.4 is shown in Fig. A.4 (red dots), overplotted with the analytical solution provided by Giacomazzo & Rezzolla (2006) (blue line) for this test. This problem was originally proposed in Noh (1987), and it describes the interaction of two streams moving in opposite direction with high Lorentz factor, so the problem is strongly relativistic. As a byproduct of this collision, four shocks are generated: two strong external fast shocks and two internal slow shocks, one of each leftgoing and rightgoing, respectively. Finally, BA5 is the only generic RMHD Riemann problem that allows all seven waves to appear after the decay of the initial discontinuity: a leftgoing fast shock, a leftgoing Alfvén wave, a leftgoing slow rarefaction, a contact discontinuity, a rightgoing slow shock, a rightgoing Alfvén wave, and a rightgoing fast shock. The initial setup includes transverse components of the velocity and magnetic field vectors. The results of this test, which are not shown in this section, show that the code is able to capture the whole set of waves predicted by the analytical solution.
Fig. A.3. From left to right (top): Density, total pressure, xvelocity, and yvelocity. From left to right (bottom): zvelocity, B_{y}, B_{z}, and logarithmic Lorentz factor at t = 0.4 for the test BA3 (red dots). The analytic solution of the Riemann problem is overplotted (solid blue line). The shock tube contains 1600 grid points in the xdirection. 
Fig. A.4. From left to right (top): Density, total pressure, xvelocity, and yvelocity. From left to right (bottom): zvelocity, B_{y}, B_{z}, and logarithmic Lorentz factor at t = 0.4 for the test BA4 (red dots). The analytic solution of the Riemann problem is overplotted (solid blue line). The shock tube contains 1600 grid points in the xdirection. 
A.4.2. Two dimensions
Cylindrical magnetized blast wave The cylindrical magnetized blast wave is a classical problem in RMHD to test the performance of the code handling MHD wave degeneracies parallel and perpendicular to the field orientation (see, e.g., Martí & Müller 2015, and references therein). Our initial setup follows the version proposed by Beckwith & Stone (2011) and adapted from Leismann et al. (2005). We considered the domain [ − 6.0, 6.0]^{2} with 1024^{2} cells and freeflow boundaries everywhere. An overpressured (p_{s} = 1) and overdense (ρ_{s} = 10^{−2}) cylinder of radius r = 0.8 was placed at the center of the grid, filled with a homogeneous ambient medium of density ρ_{a} = 10^{−4} and pressure p_{a} = 5 × 10^{−3}. A transition layer between r = 0.8 and r = 1 was established to smooth the initial jump and avoid numerical problems when the flow starts to propagate outward. All velocities were initially set to zero, and the magnetic field was aligned with the xaxis in the whole grid. For this test, we considered two degrees of magnetization: moderate (B = 0.1) and strong (B = 1.0). The adiabatic coefficient was set to the relativistic value, Γ = 4/3. The solution of the problem at t = 4.0 is shown in Fig. A.5. The difference in pressure between the cylindrical region and the ambient medium produces the expansion of the central region, which is delimited by a strong forward shock propagating radially near the speed of light. When we increased the magnetization of the medium (B = 1.0), the strong sideways magnetic confinement produced an elongated flow structure. The test preserves the symmetry with good accuracy, and no numerical artifacts or instabilities appear in our simulation.
Fig. A.5. Logarithmic gas pressure (top left), logarithmic density (top right), logarithmic Lorentz factor (bottom left), and logarithmic magnetic pressure (bottom right) at t = 0.4 for the 2D cylindrical magnetized blast wave with moderate (two left columns) and strong (two right columns) magnetization. Magnetic field lines are superposed on the magnetic pressure. We consider the Cartesian grid [ − 6, 6]×[−6, 6] with 1024 cells per spatial dimension. 
Rotor We consider the relativistic rotor problem of Del Zanna et al. (2003) in a unitary Cartesian 2D grid with a resolution of 1024^{2} cells and freeflow boundary conditions in all directions. The initial setup consisted of a disk with radius r = 0.1 rotating with an angular relativistic velocity ω = 0.95. The disk was ten times denser than the ambient medium (ρ_{r} = 10, ρ_{a} = 1) and the entire system was in pressure equilibrium with p = 1. The background was initially at rest (v = 0), and the magnetic field was aligned with the x direction in the whole grid, B_{x} = 1.0, B_{y} = 0, B_{z} = 0. An adiabatic coefficient of Γ = 5/3 was used for this test. The solution of the problem at t = 0.4 is shown in Fig. A.6. The results show complex wave patterns and torsional Alfvén waves that are generated due to the rotation of the disk. At the end of the simulation, the magnetic field lines inside the rotor are almost perpendicular with respect to the background field. The entire system is embedded in a fast rarefaction. The initial overdensity was swept away by the torsional waves and distributed in a thin oblique shell.
Fig. A.6. Logarithmic gas pressure (top left), logarithmic density (top right), logarithmic Lorentz factor (bottom left), and logarithmic magnetic pressure (bottom right) at t = 0.4 for the 2D relativistic rotor problem. Magnetic field lines are superposed on the magnetic pressure. We consider the unit square [ − 0.5, 0.5]×[−0.5, 0.5] with 1024 cells per spatial dimension. 
Shockcloud interaction A multidimensional test with relevance in astrophysical applications is the interaction of a strong shock wave with a density clump. We considered the relativistic version of the problem proposed by Mignone & Bodo (2006), where the magnetic field is orthogonal to the plane and carries a rotational discontinuity. We set a Cartesian 2D grid with dimensions [0, 1]×[−0.5, 0.5] and a resolution of 1024^{2} cells, with outflow boundaries in all directions. The shock wave was initially located at x = 0.6 such that the preshocked values (i.e., x > 0.6) are (ρ, W_{x}, p_{g}, B_{z}) = (1.0, 10, 10^{−3}, 0.5) (with the flow propagating to the left) and the shocked gas at x < 0.6 is given by (ρ, W_{x}, p_{g}, B_{z}) = (42.5942, 1.0, 127.9483, −2.12917), where ρ is the density, W_{x} is the Lorentz factor, p_{g} is the gas pressure, and B_{z} is the zcomponent of the magnetic field. The transverse components of the velocity v_{y} and v_{z} and the components of the magnetic field B_{x} and B_{y} were initially equal to zero in the whole domain. The cloud was characterized as a cylinder with radius r = 0.15 and density ρ = 10.0, centered at x = 0.8 in pressure equilibrium with the preshocked material. An adiabatic coefficient of Γ = 4/3 was used for this test. The solution of the problem at t = 1.0 is shown in Fig. A.7. Immediately after the impact between the cloud and the shock wave, the sphere experiences a strong compression that increases the density of the clump significantly. As a byproduct of this collision, a bow shock propagates to the left in the shocked material and a reverse shock is transmitted to the right, penetrating the cloud and producing a mushroomshaped structure. The overall evolution of the cloud after the impact of the shocked wave and the development of a mushroomshaped shell agrees well with Mignone & Bodo (2006). However, the higher resolution employed in our work and the use of the HLLD Riemann solver produce a less diffusive result where several complex wave patterns and eddies are noticeable in the arms of the shell.
Fig. A.7. Gas pressure (top left), density (top right), yvelocity (bottom left), and magnetic pressure (bottom right) at t = 1.0 for the 2D shockcloud interaction test. We consider the unit square [0, 1]×[−0.5, 0.5] with 1024 cells per spatial dimension. 
Relativistic OrszagTang vortex The OrszagTang vortex problem was originally proposed in Orszag & Tang (1979), and it has become a classical test for Newtonian MHD applications (Ryu et al. 1998; Londrillo & Del Zanna 2000). We adapted the relativistic version of the test proposed by Castro et al. (2017) into a two dimensionalCartesian grid with dimensions [0, 6]×[0, 6] and a resolution of 1024^{2} cells, with periodic boundary conditions in all directions. Initially, density (ρ) and pressure (p) were set to 1.0 with an adiabatic index of Γ = 4/3. The velocity field in the laboratory frame was given by
while the proposed magnetic field configuration was
Figure A.8 shows the solution to the problem at t = 4.0. Our results agree well with those obtained in Castro et al. (2017). This test validates the performance of the code describing the transition to supersonic MHD turbulence and its ability to handle the formation of shocks and shockshock interactions in the relativistic domain.
Fig. A.8. Logarithmic pressure (top left), logarithmic density (top right), logarithmic Lorentz factor (bottom left), and magnetic pressure (bottom right) at t = 4.0 for the relativistic OrszagTang vortex test. We consider the 2DCartesian grid [0.0, 6.0]×[0.0, 6.0] with 1024 cells per spatial dimension. 
Relativistic KelvinHelmholtz instability The last 2D test that we propose is the linear growth phase of the relativistic KelvinHelmholtz instability (KHI). Like the shockcloud interaction problem, this test is also relevant for astrophysical application since KHIs are commonly present in nature. In particular, this test may be useful to probe the performance of the code describing turbulence zones and the description of the growth of smallscale perturbations. KHIs develops when there is a velocity shear in a continuous flow or when there is a velocity difference across two separated flow states. Thus, it is an useful setup to test the response of the code when unstable initial conditions are provided. We considered a 2D Cartesian grid with dimensions [ − 0.5, 0.5]×[−1.0, 1.0] and resolution 512 × 1024 cells, with periodic boundaries. We followed the initial configuration given by Beckwith & Stone (2011), which in turn was adapted from Mignone et al. (2009) and exploited by other authors (Castro et al. 2017). The shear velocity and density profiles are given by
where a = 0.01 is the characteristic thickness of the shear layer and V_{shear} = 0.5 determines the velocity profile. The characteristic densities are ρ_{0} = 0.505, and ρ_{1} = 0.495. The instability was triggered by a perturbation in the transverse velocity V_{y} of the form
where A_{0} = 0.1 is the amplitude of the perturbation and σ = 0.1 the characteristic length scale. The gas pressure was constant everywhere, p = 1, and we considered an ideal EoS with Γ = 4/3. The magnetic field was aligned with the xaxis, B = (10^{−3}, 0, 0). Figure A.9 shows the solution of the problem at t = 3.0, which is qualitatively similar to the one present in the aforementioned papers. In our simulation, we also resolved the secondary vortex near x = 0.1, although Beckwith & Stone (2011) showed that it can be annihilated by using more diffusive Riemann solvers such as HLL. As pointed out by Castro et al. (2017), these secondary instabilities may be nonphysical, and their appearance depends on the Riemann solver and on the resolution employed in the simulation.
Fig. A.9. Logarithmic pressure (top left), logarithmic density (top right), logarithmic Lorentz factor (bottom left), and magnetic pressure (bottom right) at t = 3.0 for the relativistic KelvinHelmholtz test. Magnetic field lines are superposed on the magnetic pressure. We consider the 2D Cartesian grid [ − 0.5, 0.5]×[−1.0, 1.0] with 512 × 1024 cells. 
A.4.3. Three dimensions
Spherical magnetized blast wave The spherical magnetized blast wave allows us to test the ability of the code to handle parallel and perpendicular strong 3D shocks in magnetized plasma. The initial setup is similar to the one described in the cylindrical blast wave problem in 2D (see Sec. A.4.2). We considered the domain [ − 6.0, 6.0]^{3} with 256^{3} cells and free boundaries everywhere. An overpressured and overdense sphere of radius r = 0.8, density ρ_{s} = 10^{−2}, and pressure p_{s} = 1 was placed in the center of the grid, where ρ_{a} = 10^{−4}, p_{a} = 5 × 10^{−3} are the density and pressure of the ambient medium, respectively. A transition layer between r = 0.8 and r = 1 was established to smooth the initial jump and to avoid numerical problems when the flow starts to propagate. All velocities were initially set to zero. The adiabatic coefficient was set to the relativistic value, Γ = 4/3. In this test, we chose a magnetic field oblique to the grid as in Castro et al. (2017),
where B_{0} = 1 (high magnetization). The solution of the problem at t = 4.0 is shown in Fig. A.10. Due to the high magnetization of the test, the blast wave propagates outward, following the oblique field lines. The test preserves the symmetry with good accuracy, and no numerical artifacts or undesired instabilities appear in our simulation.
Fig. A.10. Logarithmic density (left), logarithmic pressure (middle), and logarithmic magnetic pressure (right) at t = 4.0 for the spherical magnetized blast wave in the xy plane at z = 0 (top) and xz plane at y = 0 (bottom). We consider the 3D Cartesian grid [ − 6.0, 6.0]×[−6.0, 6.0]×[−6.0, 6.0] with 256 cells per spatial dimension. 
3D rotor problem A 3D version of the relativistic rotor problem was considered (Mignone et al. 2009). The computational box covered the domain [ − 0.5, 0.5]^{3}, with 256^{3} cells and free boundary conditions everywhere. The initial configuration consisted of a sphere of density ρ_{s} = 10 and radius r_{s} = 0.1 that rotated around the zaxis with relativistic velocity components (v_{x}, v_{y}, v_{z}) = ω (−y, x, 0), where ω = 9.95 is the angular velocity. The sphere was in pressure equilibrium with the ambient medium, whose pressure and density were p_{a} = 1 and ρ_{a} = 1, respectively. The magnetic field was aligned with the xdirection in the whole domain, B = (1, 0, 0), and we used the nonrelativistic adiabatic index, Γ = 5/3. The solution of the problem at t = 0.4 is shown in Fig. A.11. When the sphere starts to spin, complex torsional waves and shocks propagate outward, carrying angular momentum from the sphere to the medium, producing an octogonlike disk in the xy plane (z = 0) surrounded by a shell of higher density and lower magnetic pressure, all embedded in a spherical fast rarefaction. The overall shape and the internal distribution of the flow agree well with the results presented in Mignone et al. (2009) using the HLLD Riemann solver. The test also shows the same flow distortions in the xy plane along the xaxis, which are 3D numerical pathologies that do not manifest in the 2D version of the problem. This demonstrates that our version of HLLD works similar to the version in PLUTO for 3D applications.
Fig. A.11. Logarithmic density (left), logarithmic pressure (middle), and logarithmic magnetic pressure (right) at t = 0.4 for the 3D relativistic rotor problem in the xy plane at z = 0 (top) and xz plane at y = 0 (bottom). We consider the 3D Cartesian grid [ − 0.5, 0.5]×[−0.5, 0.5]×[−0.5, 0.5] with 256 cells per spatial dimension. 
3D shockcloud interaction Finally, we considered a 3D version of the shockcloud collision problem of Mignone & Bodo (2006), following the method described in Mignone et al. (2012). We performed the simulation in the computational box [0, 1]×[−0.5, 0.5]×[−0.5, 0.5], with 256^{3} cells and free boundary conditions in the whole domain. A shockwave was initially located at x = 0.6. Upstream (x > 0.6), the preshocked material was defined by ρ = 1, p_{g} = 10^{−3}, W_{x} = 10 (v_{x} = −0.995), B_{z} = 0.5, while downstream, the medium was given by ρ = 42.5942, p_{g} = 127.9483, W_{x} = 1, B_{z} = −2.12971. In this version, we reduced the magnitude of the zcomponent of the magnetic field with respect to the 2D test by and we set B_{y} = B_{z}. Transversal velocities v_{y}, v_{z} and the component of the magnetic field parallel to the xaxis, B_{x}, were initially set to zero everywhere. At x = 0.8, an overdense sphere of radius r = 0.15 and density ρ = 10 was set in pressure equilibrium with the preshocked material and was advected with the flow. For the ideal EoS, we considered the relativistic adiabatic factor, Γ = 4/3. The solution of the problem at t = 1.0 is shown in Fig. A.12. As in the 2D simulation, immediately after the impact between the cloud and the shock wave, the sphere experiments a strong compression that increases the density of the clump significantly. As a byproduct of this collision, a bow shock propagates to the left in the shocked material, and a reverse shock is transmitted to the right, penetrating the cloud and producing a mushroomshaped structure like in the 2D shockcloud interaction test. However, due to the lower resolution employed in this case, the pattern of waves that appear at the end of the test is more diffusive.
Fig. A.12. Density (left), pressure (middle), and logarithmic magnetic pressure (right) at t = 1.0 for the 3D relativistic shockcloud interaction in the xy plane at z = 0 (top) and xz plane at y = 0 (bottom). We consider the 3D Cartesian grid [0, 1.0]×[−0.5, 0.5]×[−0.5, 0.5] with 256 cells per spatial dimension. 
All Tables
All Figures
Fig. 1. Scheme of the simulation setup in the XY plane. The radial velocity field that represents the stellar wind fills the box from the beginning. We also show a top view of the injection plane at Y = 0, where we distinguish the cylindrical nozzle (blue), the shear layer (red), and the jet cocoon (orange), where boundaries are reflecting. 

In the text 
Fig. 2. Propagation of a fiducial jet (jet B) during the first seconds after injection. The dimensions of the box are , with an effective resolution of 6 cells/R_{j}. The density distribution of the ambient medium is represented in logarithmic scale, where we limited the maximum density to ρ_{max} = 10 ρ_{w} to highlight the wind clumps. Magnetic field lines are represented in the head of the jet, where the density is highest. The 3D jet render is constructed using the jet tracer function and is colored with the same colour scale as the density distribution (although we do not show the tracer legend to avoid redundancy), where jet particles (f = 1) are represented in light blue and white and the ambient medium (f = 0) in yellow and white. A gas pressure contour (faint yellow) is also included to show the position of the jet bow shock. 

In the text 
Fig. 3. Toroidal magnetic field vectors in the same time frame as in Fig. 2. Stellar wind clumps are represented with gold isosurfaces at ρ = 3 ρ_{w}. 

In the text 
Fig. 4. Axial cuts of jet A simulation. Vertical panels: Logarithmic restmass density, logarithmic gas pressure, and logarithmic magnetic pressure at t = 1020 (R_{j}/c). Ten tracer contours are overplotted together with the gas pressure, from f ≈ 0 (white to blue) to f = 1 (yellow to white). Horizontal panels: evolution of the jet velocity at t = 310 (top), t = 610 (middle), and t = 1020 (bottom) (R_{j}/c). 

In the text 
Fig. 5. 3D model of jet A simulation. Top panel: 3D render of tracer at t = 1010 (R_{j}/c) and stellar wind clumps. A gas pressure contour (faint yellow) is included to show the position of the jet bow shock. Bottom panel: density isocontours at ρ = 0.1 ρ_{w} (red), ρ = 2.5 ρ_{w} (gold), ρ = 5.0 ρ_{w} (dark blue), and ρ = 7.5 ρ_{w} (light blue) together with a collection of magnetic field lines. 

In the text 
Fig. 6. Axial cuts of jet B simulation. Vertical panels: Logarithmic restmass density, logarithmic gas pressure, and logarithmic magnetic pressure at t = 560 (R_{j}/c). Ten tracer contours are overplotted together with the gas pressure, from f ≈ 0 (white to blue) to f = 1 (yellow to white). Horizontal panels: Evolution of the jet velocity at t = 110 (top), t = 360 (middle), and t = 560 (bottom) (R_{j}/c). 

In the text 
Fig. 7. 3D model of jet B simulation. Top panel: 3D render of tracer at t = 560 (R_{j}/c) and stellar wind clumps. A gas pressure contour (faint yellow) is included to show the position of the jet bowshock. Bottom panel: density isosurfaces at ρ = 0.1 ρ_{w} (red), ρ = 2.5 ρ_{w} (gold), ρ = 5.0 ρ_{w} (dark blue), and ρ = 7.5 ρ_{w} (light blue) together with a collection of magnetic field lines. 

In the text 
Fig. 8. Axial cuts of jet C simulation. Vertical panels: Logarithmic restmass density, logarithmic gas pressure, and logarithmic magnetic pressure at t = 720 (R_{j}/c). Ten tracer contours are overplotted together with the gas pressure, from f ≈ 0 (white to blue) to f = 1 (yellow to white). Horizontal panels: Evolution of the jet velocity at t = 210 (top), t = 490 (middle), and t = 720 (bottom) (R_{j}/c). 

In the text 
Fig. 9. 3D model of jet C simulation. Top panel: 3D render of tracer at t = 720 (R_{j}/c) and stellar wind clumps. A gas pressure contour (faint yellow) is included to show the position of the jet bow shock. Bottom panel: density isosurfaces at ρ = 0.1 ρ_{w} (red), ρ = 2.5 ρ_{w} (gold), ρ = 5.0 ρ_{w} (dark blue), and ρ = 7.5 ρ_{w} (light blue) together with a collection of magnetic field lines. 

In the text 
Fig. 10. Time evolution of the logarithm of the energy in the simulation for jet B (left) and jet C (right) at the three time frames of Figs. 6 and 8, respectively. In the right panel (jet C), we have included two additional time frames at t = 10 (R_{j}/c) and t = 110 (R_{j}/c). Energy stored by the jet plasma (f ≠ 0) is represented with dots and dashed lines, ambient medium energies (f = 0) with dots and dashdotted lines, and total energies (jet+ambient medium) with squares and solid lines. Continuous solid lines (without marks) are the theoretical curves (i.e., the injected values) for each type of energy: kinetic (red), internal (black), magnetic (blue), and total (green). Energy units are . 

In the text 
Fig. 11. Fraction of magnetic energy density with respect to total energy density for jet C in the XY (top) and ZY (bottom) plane. 

In the text 
Fig. A.1. Riemann fan for the HLLfamily of Riemann solvers: HLL (blue), HLLC (blue+black), and HLLD (blue+black+red). 

In the text 
Fig. A.2. Gas density (top), pressure (middle), and xvelocity (bottom) at t = 0.35 for the hydrodynamical sinusoidal perturbation test. The shock tube contains 2000 grid points in the xdirection. 

In the text 
Fig. A.3. From left to right (top): Density, total pressure, xvelocity, and yvelocity. From left to right (bottom): zvelocity, B_{y}, B_{z}, and logarithmic Lorentz factor at t = 0.4 for the test BA3 (red dots). The analytic solution of the Riemann problem is overplotted (solid blue line). The shock tube contains 1600 grid points in the xdirection. 

In the text 
Fig. A.4. From left to right (top): Density, total pressure, xvelocity, and yvelocity. From left to right (bottom): zvelocity, B_{y}, B_{z}, and logarithmic Lorentz factor at t = 0.4 for the test BA4 (red dots). The analytic solution of the Riemann problem is overplotted (solid blue line). The shock tube contains 1600 grid points in the xdirection. 

In the text 
Fig. A.5. Logarithmic gas pressure (top left), logarithmic density (top right), logarithmic Lorentz factor (bottom left), and logarithmic magnetic pressure (bottom right) at t = 0.4 for the 2D cylindrical magnetized blast wave with moderate (two left columns) and strong (two right columns) magnetization. Magnetic field lines are superposed on the magnetic pressure. We consider the Cartesian grid [ − 6, 6]×[−6, 6] with 1024 cells per spatial dimension. 

In the text 
Fig. A.6. Logarithmic gas pressure (top left), logarithmic density (top right), logarithmic Lorentz factor (bottom left), and logarithmic magnetic pressure (bottom right) at t = 0.4 for the 2D relativistic rotor problem. Magnetic field lines are superposed on the magnetic pressure. We consider the unit square [ − 0.5, 0.5]×[−0.5, 0.5] with 1024 cells per spatial dimension. 

In the text 
Fig. A.7. Gas pressure (top left), density (top right), yvelocity (bottom left), and magnetic pressure (bottom right) at t = 1.0 for the 2D shockcloud interaction test. We consider the unit square [0, 1]×[−0.5, 0.5] with 1024 cells per spatial dimension. 

In the text 
Fig. A.8. Logarithmic pressure (top left), logarithmic density (top right), logarithmic Lorentz factor (bottom left), and magnetic pressure (bottom right) at t = 4.0 for the relativistic OrszagTang vortex test. We consider the 2DCartesian grid [0.0, 6.0]×[0.0, 6.0] with 1024 cells per spatial dimension. 

In the text 
Fig. A.9. Logarithmic pressure (top left), logarithmic density (top right), logarithmic Lorentz factor (bottom left), and magnetic pressure (bottom right) at t = 3.0 for the relativistic KelvinHelmholtz test. Magnetic field lines are superposed on the magnetic pressure. We consider the 2D Cartesian grid [ − 0.5, 0.5]×[−1.0, 1.0] with 512 × 1024 cells. 

In the text 
Fig. A.10. Logarithmic density (left), logarithmic pressure (middle), and logarithmic magnetic pressure (right) at t = 4.0 for the spherical magnetized blast wave in the xy plane at z = 0 (top) and xz plane at y = 0 (bottom). We consider the 3D Cartesian grid [ − 6.0, 6.0]×[−6.0, 6.0]×[−6.0, 6.0] with 256 cells per spatial dimension. 

In the text 
Fig. A.11. Logarithmic density (left), logarithmic pressure (middle), and logarithmic magnetic pressure (right) at t = 0.4 for the 3D relativistic rotor problem in the xy plane at z = 0 (top) and xz plane at y = 0 (bottom). We consider the 3D Cartesian grid [ − 0.5, 0.5]×[−0.5, 0.5]×[−0.5, 0.5] with 256 cells per spatial dimension. 

In the text 
Fig. A.12. Density (left), pressure (middle), and logarithmic magnetic pressure (right) at t = 1.0 for the 3D relativistic shockcloud interaction in the xy plane at z = 0 (top) and xz plane at y = 0 (bottom). We consider the 3D Cartesian grid [0, 1.0]×[−0.5, 0.5]×[−0.5, 0.5] with 256 cells per spatial dimension. 

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