© 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 X-ray 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 strong-gravity and curved space-time; 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 shock-acceleration 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 Roche-lobe overflow (low-mass XRBs) or through the capture of stellar winds (high-mass 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 & Bosch-Ramon 2008, 2012; Perucho et al. 2010a; Bosch-Ramon & Barkov 2016; de la Cita et al. 2017; Barkov & Bosch-Ramon 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; Bosch-Ramon et al. 2011; Yoon et al. 2011; Monceau-Baroux et al. 2014).

For HMXBs, jet-wind 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 & Bosch-Ramon 2008; Perucho et al. 2010a). Moreover, the one-sided 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 proton-proton collisions (see, e.g., Romero et al. 2003; Molina et al. 2019). These processes may be at the origin of the gamma-ray 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 X-1 (e.g., Albert et al. 2007; Zanin et al. 2016), and Cygnus X-3 (e.g., Albert et al. 2007; Zanin et al. 2016).

Several attempts have been made in the past to describe jet-wind interactions from the point of view of theory and numerical simulations. Perucho & Bosch-Ramon (2008) performed numerical two-dimensional (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 three-dimensional (3D) hydrodynamical simulations, were found by Perucho et al. (2010a). This work also confirmed that jets with total power L_j ≲ 10³⁶Ṁ_w,−6 erg s⁻¹ (Ṁ_w,−6 ≡ Ṁ_w/10⁻⁶ M_⊙ yr⁻¹: the wind mass-loss rate) may be disrupted by the impact of the stellar wind, conditioning the radio-band detection of this population of XRBs. Perucho & Bosch-Ramon (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³⁷ Ṁ_w,−6 erg s⁻¹. Numerical calculations of the high-energy emission from one clump-jet 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 well-known XRBs, Cygnus X-1 and Cygnus X-3. This work was then extended by Yoon et al. (2016) and Bosch-Ramon & 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 large-scale 3D RHD simulations of jet outbreak and early propagation, finding that radiative cooling effects are more relevant for Cygnus X-3 than for Cygnus X-1. Furthermore, Barkov & Bosch-Ramon (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³⁷ Ṁ_w,−6 erg s⁻¹ 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 hoop-stress or by acting as a destabilising agent via the growth of current-driven 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 mass-loaded through mixing produced by small-scale 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; Cruz-Osorio et al. 2022) and the morphological characterization of kilo-parsec scale magnetized jets (see, e.g., Martí 2019, and references therein). Different authors focused on the properties of stationary one-dimensional 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; Moya-Torregrosa et al. 2021). Mizuno et al. (2007) studied the long-term stability of 3D magnetized spine-sheath relativistic jets. The first high-resolution 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 Poynting-flux-driven 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 & Bosch-Ramon (2012). These simulations can be used as a numerical benchmark upon which to build new high-resolution 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 jet-wind 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 high-resolution shock-capturing (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 X-1 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¹² 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⁹ cm at a distance to the compact object y₀ = 6 × 10¹⁰ cm, propagating along the y-axis 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 space-time) can be disregarded from these simulations. For the ambient medium, we considered an inhomogeneous stellar wind with a radial velocity field v₀ = 2 × 10⁸ cm s⁻¹, a mass-loss rate ∼10⁻⁶ M_⊙ yr⁻¹, a temperature T_w ∼ 10⁴ K, thermal pressure p_w ∼ 1, 6 × 10⁻³ erg cm⁻³, and a mean density ρ_w ∼ 3 × 10⁻¹⁵ g cm⁻³ (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 OB-type star (see, e.g., Perucho & Bosch-Ramon 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³⁵ and ∼10³⁷ erg s⁻¹, 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).

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 & Bosch-Ramon (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.

Time-resolved 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¹ (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 & Bosch-Ramon (2012), where clouds were modeled as Gaussian-shape 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 log-normal single-point 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 log-normal parent distribution was μ = 1.0 and the variance σ² = 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 wave-number 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₀ 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

(1)

(2)

(3)

(4)

where θ, ϕ are the characteristic angles of spherical coordinates with respect to the location of the star,

(5)

and the subscript j = 1 indicates that the values of density and velocity in the y-boundary 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 counter-jet, 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⁻⁸ (in units of ρ₀c²). 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)

(6)

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),

(7)

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,

(8)

where W is the Lorentz factor. When we integrate Eq. (8) by separation of variables, the gas pressure profile p(r) yields

(9)

where

(10)

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 divergence-free 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,

(11)

where h(r) is the specific enthalpy. The stellar wind enters the grid from the X-MAX boundary, while in all the remaining boundaries of the computational box (i.e., X-MIN, Y-MAX, Z-MIN, Z-MAX), we considered free-flow conditions.

The simulations presented in Sect. 3 were performed with Lóstrego v1.0, a new conservative, finite-volume 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 second-order Godunov-type 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 shock-flattening 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 third-order TVD-preserving Runge-Kutta (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 divergence-free 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 long-term survival of the whole simulation. This might become especially relevant when jets are injected as axisymmetric top-hat 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 top-hat 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)

(12)

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 rest-mass 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 rest-mass 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/Rj. 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¹² 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 one-dimensional relativistic approximation of Martí et al. (1997) The relativistic velocity of the jet head can be estimated based on a one-dimensional 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 rest-mass 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 rest-mass 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 low-density 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 large-scale 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 rest-mass density, logarithmic gas pressure, and logarithmic magnetic pressure at t = 1020 (Rj/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) (Rj/c).

Fig. 5. 3D model of jet A simulation. Top panel: 3D render of tracer at t = 1010 (Rj/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¹² 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 one-dimensional relativistic estimate, v_1D ≈ 0.43 c.

Fig. 6. Axial cuts of jet B simulation. Vertical panels: Logarithmic rest-mass density, logarithmic gas pressure, and logarithmic magnetic pressure at t = 560 (Rj/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) (Rj/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 long-term 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 Kelvin-Helmholtz 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 (Rj/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.

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¹² 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 one-dimensional 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 one-dimensional 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 rest-mass density, logarithmic gas pressure, and logarithmic magnetic pressure at t = 720 (Rj/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) (Rj/c).

The effect of the toroidal field on the jet dynamics is visible from the early jet evolution by the development of current-driven instabilities that destabilize the jet beyond the first reconfinement shock. The consequent increase in cross-section 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 Kelvin-Helmholtz 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 (Rj/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¹² vs. y ≃ 1.5 × 10¹² 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¹¹ 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 jet-to-cocoon 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¹² cm in t ≃ 110 s, and the second takes almost twice as long to reach y ≃ 2 × 10¹² 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 long-term 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 current-driven 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 long-term 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 large-angle 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 Kelvin-Helmholtz 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; Moya-Torregrosa 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 & Bosch-Ramon (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 wind-jet 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 long-term 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,

(13)

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 top-hat 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 (Rj/c) and t = 110 (Rj/c). Energy stored by the jet plasma (f ≠ 0) is represented with dots and dashed lines, ambient medium energies (f = 0) with dots and dash-dotted 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 dash-dotted 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 (dash-dotted 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%², 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 Bosch-Ramon & 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., Bosch-Ramon & Khangulyan 2009). In these systems, strong magnetic field dissipation through magnetic reconnection could also be relevant for particle acceleration (Bosch-Ramon 2012; Bosch-Ramon & 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 jet-wind 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 clump-jet interaction was simulated using relativistic hydrodynamics), the jet dynamics was semi-analytical. 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 frozen-in 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 space-wise 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 X-rays 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 gamma-ray reprocessing in the complex jet-wind 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 Bosch-Ramon & Barkov (2016) and shown by Barkov & Bosch-Ramon (2022), the impact of orbital motion n longer time-scales 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 & Bosch-Ramon 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 X-1 (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 (Bosch-Ramon & Barkov 2016; Barkov & Bosch-Ramon 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 quasi-steady 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 accretion-ejection 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 large-scale 3D hydrodynamical simulations (Charlet et al. 2022) show that for two fiducial HMXBs (namely, Cygnus X-1 and Cygnus X-3), thermal cooling (i.e., free-free 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 X-1, although they can be relevant for the evolution of Cygnus X-3 due to the denser stellar winds of the Wolf-Rayet 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 jet-wind simulations, including a toroidal magnetic field in transversal equilibrium with the gas to analyze its role for jet dynamics and long-term 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 non-negligible 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 & Bosch-Ramon (2008); Perucho et al. (2010a) it was concluded that this depended basically on the jet power, whereas in this work, we see that low-power jets threaded by relatively weak fields could propagate longer distances without being destabilized. However, more powerful jets, in contrast, might be disrupted by current-driven 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 & Bosch-Ramon 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 PID2019-105510GB-C31. M.P. and J.M-M. acknowledge support by the Spanish Ministerio de Ciencia through grants PID2019-107427GB-C33 and PGC2018-095984-B-I00, and from the Generalitat Valenciana through grant PROMETEU/2019/071. M.P. acknowledges additional support from the Spanish Ministerio de Ciencia through grant PID2019-105510GB-C31. V.B.-R. acknowledges financial support from the State Agency for Research of the Spanish Ministry of Science and Innovation under grant PID2019-105510GB-C31 and through the ”Unit of Excellence María de Maeztu 2020-2023” award to the Institute of Cosmos Sciences (CEX2019-000918-M), and by the Catalan DEC grant 2017 SGR 643. V.B-R. 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

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 shock-capturing (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 reconstruct-solve-average (RSA) strategy: first, cell-average 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) Runge-Kutta algorithms. Finally, an inversion scheme is applied to recover the primitive variables from the time-evolved solution. Several one-dimensional 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³ and Cartesian coordinates can be written as a system of conservation laws (Komissarov 1999b),

(A.1)

where U = {D, S^j, τ, B^j} is a vector of conserved variables, and Fⁱ are the vector of fluxes for each spatial direction. Both the vector of conserved variables U and the vector of fluxes Fⁱ can be expressed in terms of a set of primitive variables V = {ρ, v^j, p, B^j} using the following relations:

(A.2)

(A.3)

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 rest-mass density ρ, fluid three-velocity 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^*,

(A.4)

where ϵ is the specific internal energy. The total pressure p^* of the magnetic fluid is given by

(A.5)

In terms of the vectors v and B in the laboratory frame, the four-velocity u^α and magnetic field b^α covariant vectors are

(A.6)

(A.7)

(A.8)

with the normalizations

(A.9)

(A.10)

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 divergence-free constraint,

(A.11)

Finally, the system of RMHD equations is closed with an equation of state (EoS) of the form

(A.12)

which relates the thermodynamic primitive variables with each other. For the particular case of ideal gases, this relation becomes

(A.13)

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 cell-averaged 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 semi-discretization in a rectangular grid of size L_x × L_y × L_z and resolution N_x × N_y × N_z,

(A.14)

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,

(A.15)

while the numerical fluxes represent an approximation to the average in the surface of the cell,

(A.16)

(A.17)

(A.18)

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 cell-centered values by means of interpolation routines that are designed to reconstruct a piece-wise polynomial approximation inside each cell preserving the monotonicity and TV-stability 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 piece-wise linear method (PLM) in the x-direction follows

(A.19)

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² > 1, we reduce the method to first order (i.e., Godunov reconstruction). High-order 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 cell-centered 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

(A.20)

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 Harten-Lax-van 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 root-finding 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 Rankine-Hugoniot 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 five-wave 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 Rankine-Hugoniot 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 HLL-family 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 predictor-corrector methods. Numerical time integration is based on second- (RK2) and third-order (RK3) total variation dimishing (TVD) Runge-Kutta time integrators described by Shu & Osher (1989). The second-order method (RK2) follows

(A.21)

while the third-order method (RK3) is given by

(A.22)

where L(U⁽ⁿ⁾) is the upwind differencing operator,

(A.23)

A.2.5. Divergence-free constraint

In general, HRSC methods do not preserve the divergence-free 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 eight-wave method (Powell 1994), the hyperbolic-parabolic 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 divergence-free constraint is fulfilled up to machine round-off 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 Godunov-type 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,

(A.24)

(A.25)

(A.26)

where are the staggered field components, which represent an average of the magnetic field at cell surface,

(A.27)

(A.28)

(A.29)

The electromotive forces (emf) Ωⁱ are discrete representations of the components of Ω = v × B defined at cell corners. To compute them, we follow the flux-CT formalism (Balsara & Spicer 1999; Giacomazzo & Rezzolla 2007), taking the simple average of the neighboring upwind numerical fluxes,

(A.30)

(A.31)

(A.32)

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 plane-parallel algorithm for grid-align flows (Gardiner & Stone 2005; Mignone & Del Zanna 2021). Thus, we modified the definitions of Eqs. A.30 to A.32 following the CT-contact 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 grid-aligned planar flows. In the CT-contact approximation, the average electromotive forces are modified according to

(A.33)

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,

(A.34)

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,

(A.35)

The cell-center magnetic field can be recovered from the staggered representation using a simple linear interpolation at the end of the CT step,

(A.36)

(A.37)

(A.38)

The interpolations performed to obtain the required fluxes at cell edges in the CT-contact formalism, and to recover the cell-centered magnetic field from the staggered solution, limit the accuracy of the algorithm to second order. Other strategies such as the UCT-HLL(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 speed-up 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

(A.39)

(A.40)

where E′ = E − D, Z′ = Z − D, Z = DhW, S_B = S ⋅ B and |u|² = W²|v|². Equation A.39 can be solved for Z′ by means of a one-dimensional Newton-Raphson iterative method using Eq. A.40 to express |u|² as a function of Z′. The method involves the computation of the derivative

(A.41)

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

(A.42)

For ideal gases, ∂p/∂χ and ∂p/∂ρ become

(A.43)

Moreover, the variable χ can be written as a function of Z′ as

(A.44)

and its derivative with respect to Z′ as

(A.45)

where

(A.46)

On the other hand,

(A.47)

When Z′ and p have been found with some desirable degree of accuracy, the inversion scheme is completed by recovering the three-velocity and the density of the gas through

(A.48)

(A.49)

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|² using a bisections root finder with the following function:

(A.50)

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 cell-centered magnetic field makes the algorithm inconsistent because the conserved variables have been evolved according to the fluxes computed from the cell-centered 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)

(A.51)

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⁽¹⁾ to complete the relativistic correction of the momentum and energy,

(A.52)

(A.53)

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 shared-memory 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 CT-contact 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 direction-wise boundaries in the three spatial directions, but also the corners of each sub-box. Nevertheless, this information can be also transmitted direction-wise 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 one-dimensional 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 second-order 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 third-order TVD-preserving Runge-Kutta algorithm for time integration with CFL = 0.3. The magnetic field divergence-free constraint is preserved with the CT method, where electromotive forces were averaged according to the CT-contact 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 one-dimensional 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 x-velocity (bottom) at t = 0.35 for the hydrodynamical sinusoidal perturbation test. The shock tube contains 2000 grid points in the x-direction.

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, BA1-5), 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 one-dimensional 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 left-going fast rarefaction wave, a left-going compound wave (that only appears in the numerical solution), a contact discontinuity, a right-going slow shock, and a right-going 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 well-resolved with our code, and we clearly identify a left-going fast rarefaction wave, a left-going slow rarefaction, a contact discontinuity, a right-going slow shock, and a right-going 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 right-going shocks are under-resolved, 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 left-going and right-going, 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 left-going fast shock, a left-going Alfvén wave, a left-going slow rarefaction, a contact discontinuity, a right-going slow shock, a right-going Alfvén wave, and a right-going 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, x-velocity, and y-velocity. From left to right (bottom): z-velocity, By, Bz, 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 x-direction.

Fig. A.4. From left to right (top): Density, total pressure, x-velocity, and y-velocity. From left to right (bottom): z-velocity, By, Bz, 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 x-direction.

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]² with 1024² cells and free-flow boundaries everywhere. An overpressured (p_s = 1) and overdense (ρ_s = 10⁻²) cylinder of radius r = 0.8 was placed at the center of the grid, filled with a homogeneous ambient medium of density ρ_a = 10⁻⁴ and pressure p_a = 5 × 10⁻³. 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 x-axis 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² cells and free-flow 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.

Shock-cloud 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² 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⁻³, 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 z-component 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 mushroom-shaped structure. The overall evolution of the cloud after the impact of the shocked wave and the development of a mushroom-shaped 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), y-velocity (bottom left), and magnetic pressure (bottom right) at t = 1.0 for the 2D shock-cloud interaction test. We consider the unit square [0, 1]×[−0.5, 0.5] with 1024 cells per spatial dimension.

Relativistic Orszag-Tang vortex The Orszag-Tang 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² 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

(A.54)

while the proposed magnetic field configuration was

(A.55)

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 shock-shock 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 Orszag-Tang vortex test. We consider the 2DCartesian grid [0.0, 6.0]×[0.0, 6.0] with 1024 cells per spatial dimension.

Relativistic Kelvin-Helmholtz instability The last 2D test that we propose is the linear growth phase of the relativistic Kelvin-Helmholtz instability (KHI). Like the shock-cloud 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 small-scale 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

(A.56)

(A.57)

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.505, and ρ₁ = 0.495. The instability was triggered by a perturbation in the transverse velocity V_y of the form

(A.58)

where A₀ = 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 x-axis, B = (10⁻³, 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 Kelvin-Helmholtz 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]³ with 256³ cells and free boundaries everywhere. An overpressured and overdense sphere of radius r = 0.8, density ρ_s = 10⁻², and pressure p_s = 1 was placed in the center of the grid, where ρ_a = 10⁻⁴, p_a = 5 × 10⁻³ 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),

(A.59)

where B₀ = 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]³, with 256³ 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 z-axis 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 x-direction 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 octogon-like 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 x-axis, 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 shock-cloud interaction Finally, we considered a 3D version of the shock-cloud 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³ cells and free boundary conditions in the whole domain. A shock-wave was initially located at x = 0.6. Upstream (x > 0.6), the pre-shocked material was defined by ρ = 1, p_g = 10⁻³, 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 z-component 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 x-axis, 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 mushroom-shaped structure like in the 2D shock-cloud 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 shock-cloud 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/Rj. 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 rest-mass density, logarithmic gas pressure, and logarithmic magnetic pressure at t = 1020 (Rj/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) (Rj/c).

In the text

	Fig. 5. 3D model of jet A simulation. Top panel: 3D render of tracer at t = 1010 (Rj/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 rest-mass density, logarithmic gas pressure, and logarithmic magnetic pressure at t = 560 (Rj/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) (Rj/c).

In the text

	Fig. 7. 3D model of jet B simulation. Top panel: 3D render of tracer at t = 560 (Rj/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. 8. Axial cuts of jet C simulation. Vertical panels: Logarithmic rest-mass density, logarithmic gas pressure, and logarithmic magnetic pressure at t = 720 (Rj/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) (Rj/c).

In the text

	Fig. 9. 3D model of jet C simulation. Top panel: 3D render of tracer at t = 720 (Rj/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 (Rj/c) and t = 110 (Rj/c). Energy stored by the jet plasma (f ≠ 0) is represented with dots and dashed lines, ambient medium energies (f = 0) with dots and dash-dotted 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 HLL-family 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 x-velocity (bottom) at t = 0.35 for the hydrodynamical sinusoidal perturbation test. The shock tube contains 2000 grid points in the x-direction.
In the text

	Fig. A.3. From left to right (top): Density, total pressure, x-velocity, and y-velocity. From left to right (bottom): z-velocity, By, Bz, 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 x-direction.
In the text

	Fig. A.4. From left to right (top): Density, total pressure, x-velocity, and y-velocity. From left to right (bottom): z-velocity, By, Bz, 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 x-direction.
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), y-velocity (bottom left), and magnetic pressure (bottom right) at t = 1.0 for the 2D shock-cloud 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 Orszag-Tang 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 Kelvin-Helmholtz 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 shock-cloud 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