© A. Charlet et al. 2022

1. Introduction

Jets are an ubiquitous manifestation of the activity of compact objects that are at the origin of the microquasar phenomenon (Romero et al. 2017). High-mass microquasars (HMMQ) are a subclass of high-mass X-ray binaries (HMXRB) and are composed of a black hole (BH) and a massive star companion. HMMQ launch powerful collimated jets (e.g., Mirabel & Rodriguez 1999; Gallo et al. 2005) at relativistic speeds either from the accretion disk of the compact object through the Blandford-Payne magneto-centrifugal ejection mechanism (Blandford & Payne 1982), or the BH magnetosphere through the Blandford-Znajek mechanism (Blandford & Znajek 1977). Jets are of interest as integral parts of the astrophysical objects harboring them, but also because their impressive stability due to collimation allows them to extend orders of magnitude farther than their injection scale (Migliori et al. 2017; Gourgouliatos & Komissarov 2018a), offering a very powerful way to study their environment or their contribution to observed thermal (radio), but also (nonthermal) emissions (Malzac 2014; Zdziarski et al. 2014; Rodriguez et al. 2015; Molina et al. 2019; Albert et al. 2021; Motta et al. 2021).

The HMMQ jets closely resemble to scaled-down jets from active galactic nuclei (AGN) in regard of the overall energy released by accretion. However, in HMMQs, the ambient medium is dominated by the powerful winds of the stellar companion, which are often the source of accretion for the compact object. This stellar wind dominates the environment in which the jet will be launched and evolve, which makes the jet propagation in HMMQ fundamentally different from jets in AGNs and low-mass microquasars.

The effects of stellar wind on jets were studied by Perucho et al. (2010a), who performed 3D simulations of relativistic hydrodynamical jets with a simulation box scale shorter than the orbital separation. They suggested that jets with a power of a few 10³⁶ erg s⁻¹ can be disrupted by the wind through the effect of the Kelvin-Helmholtz instability (KHI). A companion paper (Perucho et al. 2010b) highlighted the formation and evolution of the recollimation shock and its potential role in particle acceleration in HMMQs. A larger-scale nonrelativistic study was performed by Yoon & Heinz (2015) with a simulation box scale of ∼15 orbital distance, focusing on jet bending at larger scales and obtaining a simple analytical formula for the asymptotic bending angle. A follow-up study by Yoon et al. (2016) reconsidered the formation of a recollimation shock, emphasizing that such a shock is likely present in Cygnus X-1, while the situation in Cygnus X-3 is less clear. Several papers pointed to the fact that stellar winds are more clumpy than homogeneous and explored associated consequences for the jets (Perucho & Bosch-Ramon 2012; de la Cita et al. 2017).

This paper adds to this picture with a set of 3D hydrodynamical simulations that distinguish themselves from existing work in that they are relativistic, include radiative cooling, and follow the jet evolution over comparatively large spatial distances of about 20 to 75 orbital separations. The comparatively large spatial domain allows us to follow jet evolution from an initial smooth phase through the nonlinear growth of instabilities to a turbulent, autosimilar state, thereby creating a larger-scale perspective for some of the results cited above. The parameter study we perform with this general setup is anchored at system parameters inspired by Cygnus X-1 (Orosz et al. 2011), where cooling is moderate, and Cygnus X-3 (Zdziarski et al. 2013), where a strong cooling effect occurs due to the combination of a stronger stellar wind, magnetic field, and luminosity with an orbital separation that is one order of magnitude smaller than in Cygnus X-1.

We obtain typical timescales for the initial instability growth depending on the various parameters and the presence of cooling, highlighting the importance of the beam internal shocks in the growth of the KHI and therefore jet structure and dynamics. Cooling is found to play a role only in Cygnus X-3 on the scales covered by our simulations. Once cooling becomes dominant, the jet cocoon is immediately blown away by the stellar wind. Our simulations further suggest that the parameter sensitivities we explored somewhat diminish or are more difficult to clearly diagnose later on, when the jet has become fully turbulent. We confirm earlier findings that the jet is broken when the wind power is too strong. We find a strong instability developing at the jet beam – cocoon interface that destroys the beam. A turbulent expanding region develops subsequently that eventually expands away from the orbital plane, and the jet is recovered.

The layout of the article is as follows: in Sect. 2 we present the physics of hydrodynamical relativistic jets, our models for radiative losses, and the numerical setup and methods we used in our parameter study of jet outbreak. In Sect. 3 we discuss the jet propagation, destabilization, and structure with the cooling and parameters. In Sect. 4 we discuss the validity and limitations of our results before we summarize and conclude in Sect. 5.

2. Physical scenario and numerical methods

2.1. Relativistic equations

The equations for special relativistic hydrodynamics (SRHD) can be written in the form

(1)

where 𝒰 contains the “conservative” variables and ℱⁱ the corresponding fluxes in the i direction, which are given by

(2)

D = γρ, Sⁱ = γ²ρhvⁱ, and τ = γ²ρh − p are the conservatives variable, with γ the Lorentz factor and h the specific enthalphy. ρ, vⁱ, and p are the rest-mass density, velocity, and thermal pressure, respectively. They are called the “primitive” variables. The derivation of these equations can be found in textbooks such as Landau & Lifshitz (1959) and Mihalas & Mihalas (2013). Additionally, a passive tracer J distinguishing the jet material (J = 1) from the ambient medium (J = 0) is advected independently, following

(3)

The system is closed by an adiabatic equation of state (EoS) with constant adiabatic index Γ taken equal to 5/3 for both wind and injected material. A value of Γ = 4/3 is better suited to model flows with high Lorentz factor. We chose the classical value of 5/3 for our mildly relativistic jets. A 1D comparison between these two values is given in Appendix C.5. We find that changing the adiabatic index has an almost negligible impact on the jet head propagation speed, but jets with Γ = 5/3 value present a more advanced front shock than jets with Γ = 4/3. This may translate into a more extended cocoon in 2D and 3D jets. The adapted inversion method to recover the primitive from the conservative variables is taken from Del Zanna & Bucciantini (2002) and is detailed in Appendix C.2. In the rest of this paper, quantities with subscript b, ic, and oc refer to the beam, inner cocoon, and outer cocoon, respectively. They are defined properly in Sect. 3.1. Values with subscript j refer to injection parameters, and w refers to the stellar wind.

2.2. Jet propagation

We generalize the model for the propagation of a relativistic jet as derived by Martí et al. (1997) and Mizuta et al. (2004), for example, where multidimensional effects are neglected and 1D momentum balance between the beam with velocity v_b and the ambient gas is assumed in the rest frame of the contact discontinuity at the head of the jet, to the case of an ambient medium with its own (nonrelativistic) flow speed v_w. We obtain the following expression for the jet head velocity (details are given in Appendix A.1):

(4)

with η^* the injected-to-ambient ratio of inertial density γ²ρh = γ²(ρ + Γ₁p/c²), Γ₁ ≡ Γ/(Γ − 1). Parameters η and η^* are therefore linked by the relation

(5)

Equation (4) recovers the equation derived in Martí et al. (1997) and Mizuta et al. (2004), for instance, by taking v_w = 0. This model is useful to describe the early jet evolution, but several effects such as the growth of instabilities limit its applicability to longer-term dynamics. We can also cite jet propagation models that take deceleration into account, such as the extended Begelman-Cioffi model from Scheck et al. (2002) and the decelerated momentum balance from Mukherjee et al. (2020), but these models are not adapted to fit the dynamics of our HMMQ jets as they were developed in the context of AGN jets.

2.3. Linear growth of the Kelvin-Helmholtz instability

During the jet propagation, various hydrodynamical instabilities can be triggered and in time perturb the beam, reducing the effective beam speed at the front shock and decelerating the jet. An overview is given in Appendix A.2. Here we only mention the KHI at a relativistic flow interface. For the relativistic case we are interested in, we derive this dispersion equation from the resonance condition in Hanasz & Sol (1996) (for details, see again Appendix A.2.1),

(6)

with η_c ≡ ρ_b/ρ_ic the ratio of rest mass densities of the beam and inner cocoon, R ≡ r_ic/r_b the radius ratio of the sheet and the core, k_x the wave number in the jet propagation direction, and n an integer number. To derive growth time from our simulations, densities are measured as a volume-averaged value over each jet zone, while beam and inner cocoon radius are derived from the respective measured volume and length in the propagating direction of the zone by approximating both as coaxial cylinders. This equation is solved using a Newton-Raphson method and leads to the calculation of the wavenumber k_x, maximizing ω with the densities and radii as parameters for the first four modes (n ≤ 4). We then use the maximum corresponding ω to derive the linear growth time for the KHI. These timescales correspond to the linear growth times, whereas the observed growth rate in our simulations will be significantly higher due to nonlinear effects. They are still of interest when we compare them for different runs, as the relation between different linear growth time is the same as the observed runs.

2.4. Radiative processes

Following Bodo & Tavecchio (2018), radiative losses can be added in SRHD equations by introducing a source term in Eq. (1),

(7)

Radiative losses occur by means of four main phenomena: inverse Compton (IC) scattering, free-free (or Bremsstrahlung) emission, synchrotron emission, and line and recombination cooling.

When an external photon field acts as seeds for IC scattering, the emission pattern is anisotropic in the comoving frame of the emitting region, causing the jet to recoil. According to Ghisellini & Tavecchio (2010), however, this recoil can be neglected in the case of an ion-electron plasma because the majority of the jet momentum is transported by the ions. Both free-free emission and line cooling are dominated by the collisions between electrons and ions, which are isotropic in the fluid rest frame. This logic also applies to synchrotron losses because the pitch angle (between electron speed and magnetic field direction) distribution is isotropic, which results in no global momentum loss due to these radiative processes. We can then model the effect of the various radiative losses as a single energy-loss term,

where P_rad = P_IC + P_syn + P_ff + P_line is the volumic power losses due to inverse Compton scattering, synchrotron emission, Bremsstrahlung emission, and line and recombination cooling. Detailed expressions for each individual loss term including the derivation of the first two from a Maxwell-Jüttner distribution of the electrons and the necessary adaptations we made for them to be compatible with SRHD are given in Appendix B. The power-law exponents of each process are compiled in Table B.1. The evolution of the different loss terms with temperature for both types of runs is illustrated Fig. 1.

Fig. 1. Evolution of the loss processes with rest frame temperature in Cygnus X-1 (left) and Cygnus X-3 (right). The thin black lines shows the various temperature scalings detailed Table B.1. Line recombination losses (“line”) are not drawn for T > 107.7 as they are disabled above this temperature. Colored shading shows synchrotron losses (“syn”) when the stellar magnetic field B⋆ is either multiplied or divided by 5, and the shading around inverse Compton losses (“ic”) illustrates a ±10% uncertainty on T⋆. The values for the physical quantities have been chosen in the outer cocoon leeward side in fiducial runs. They allow for a clear showcasing of the loss scaling, but are not representative of the losses over the simulation. We refer to Fig. 7 to compare them over the whole jet. For Cygnus X-1: ρ = 10−15 g cm−3, B⋆ = 10 G, T⋆ = 3 × 104 K, (x, y, z) = (1.5 × 1012, 4.05 × 1013, 4 × 1013) cm, (vx, vy, vz) = (108, 108, 108) cm s−1. For Cygnus X-3: ρ = 10−14 g cm−3, B⋆ = 100 G, T⋆ = 8 × 104 K, (x, y, z) = (1.5 × 1012, 3.05 × 1013, 3 × 1013) cm, (vx, vy, vz) = (108, 5 × 108, −103) cm s−1.

The corresponding cooling time can be written as

(8)

with the isobaric (ṗ = 0) cooling time and the isochoric () cooling time, see Appendix B.6 for the related derivation. Taking γ ≈ 1 as a first approximation, t_c, p ≈ 10²¹ρ/P_rad and t_c, ρ ≈ 1.5p/P_rad. In the range of density and pressure observed in our simulations, the isochoric cooling time is consistently the shortest by about two orders of magnitude.

This work is limited to the case of thermal optically thin plasmas: photons produced in these processes can freely escape without interacting with the gas and thus carry all the energy away from the jet. This approximation is true for X-rays, but does not hold at every frequency. Verifying this hypothesis would require more specific investigations. We also applied our cooling model to the ambient medium, even though it is optically thick for optical and UV lines in Cygnus X-1 due to the properties of O-type star winds, and to optical and UV continuum in Cygnus X-3 because the companion star is a Wolf-Rayet (WR) star. We modeled this effect by setting up a temperature floor slightly lower than the surface temperature of the star to ensure that the medium did not cool to nonphysical values. This value is arbitrary and has no impact on the jet dynamics as it has close to no effect on the ambient inertial density, and the jet overpressure is such that modifying the ambient pressure has no effects.

The radiative processes have different scalings with rest mass density, temperature, and distance to the star in the orbital frame. Therefore the dominance of one or two processes over the others may vary with time and space. The power-law exponents of each process detailed Appendix B are compiled in Table B.1. The evolution of the different loss terms with temperature for both types of runs is shown Fig. 1.

2.5. A-MaZe simulation toolkit

We performed 3D simulations using the hydrodynamical module from the A-MaZe simulation toolkit (Walder & Folini 2000; Folini et al. 2003; Melzani et al. 2013). It uses the method of lines, a semidiscretized finite-volume method: after discretizing in space, the resulting system of ordinary differential equations is solved with a forward Euler scheme. Fluxes are computed by the (stabilized) Lax-Friedrichs approximation using a second-order reconstruction based on min-mod limiters. The equations, the solution method, and a benchmark for the accuracy of this method are all reported in Appendix C. A relativistic solver was implemented by adding the recovery of primitive variables using the inversion method detailed Appendix C.2.

Our simulations were set in a static grid made of five refinement levels centered on the jet injection nozzle, as shown in Fig. 2. Cells from the coarse grid had a 4 × 10¹¹ cm edge, and the edge of the highest-level cells was 64 times lower for a maximum resolution of 6.25 × 10⁹ cm. The number of coarse level grid cells was 250 × 200 × 200 and 250 × 150 × 150 for Cygnus X-1 and Cygnus X-3, respectively. The associated physical domain sizes are given in Table 1. The cfl number was set to 0.15. The time step was refined along with the spatial grid. On the coarse grid, it was about 2s for Cygnus X-1 and 5s for Cygnus X-3. The jet was injected perpendicular to the orbital plane (y − z plane) by fixing (ρ_j, v_j, T_j) on a few cells at x = 0, always imposing at least 20 cells of the finest grid to fix the diameter of the beam. The environment was set by fixing the wind velocity and density at the stellar surface, resulting in an isotropic wind with constant speed modulus and density in r⁻². The boundary condition the at x = 0 plane was reflective, while the other boundaries of the simulation grid had outflow conditions.

Fig. 2. Structure of the computational grid over the whole domain, illustrated for a rest mass density slice along the plane containing the star and jet center. The density scale extends from 10−19 (deep blue) to 10−12 g cm−3 (red), same as Fig. 3. The grid contains five refinement levels: the whole domain is refined twice (green and blue interfaces) by a factor of 4 and then twice more (cyan and magenta interfaces) by a factor of 2 to attain a factor of 64 in the finest levels, which are better shown in the zoom into the jet injection at the orbital scale in the bottom right corner of the image.

2.6. Covered parameter space

We defined runs CygX1 and CygX3 as our fiducial runs for Cygnus X-1 and Cygnus X-3, respectively. The main parameter values for these two runs are given in Table 1. The choice of physical values was inspired by Orosz et al. (2011) and Yoon & Heinz (2015) for Cygnus X-1 and by Orosz et al. (2011) and Dubus et al. (2010) for Cygnus X-3. The parameter choices for the various sensitivity studies are listed in Tables D.1 and D.3.

Table 1 shows the value of the environment parameters relevant for jet radiative losses as well as the parameters of the respective fiducial runs. The characteristics of the Cygnus X-3 system mean the radiative losses will be stronger overall. We chose a stronger magnetic field base value for Cygnus X-3 to compensate for the addition of the distance scaling detailed in Appendix B.5. The luminosity of the two companion stars are similar because the companion star in Cygnus X-3 is hotter but smaller. Thus the smaller orbital distance implies stronger synchrotron and inverse Compton losses by a factor 100. The beam density ρ_b was chosen to be ten times greater than in the Cygnus X-1 runs, implying stronger line and free-free losses by a factor 100 also. Both beams roughly have the same internal energy density, but will cool a ∼100 times faster in Cygnus X-3 case.

2.7. Post-processing

To perform a quantitative analysis of our simulations, we identified each computational cell according to the following rules for the various interfaces of the jet: the separation between ambient material and outer cocoon was made at p = .01 Ba and T = 10⁷ K, the working surface between inner and outer cocoon was defined where J = 0.05, following the definition for the mixing layer in Perucho et al. (2004a), and cells were considered part of the beam if ζ ≡ (v_x/v_j)J > 0.8. This criterion was defined in Yoon & Heinz (2015). We found that choosing 0.8 as the threshold value identified the beam up to the reverse shock with more success than a criterion that was purely based on J, especially in the later phases of the jet outbreak when beam and inner cocoon mix. The value J = 0.05 was found to correctly separate the low-density high-temperature inner cocoon from the outer cocoon. These criteria are deemed correct in the sense that their limits correspond to the jumps in the various physical quantities between jet zones. Redundant cells between the different refinement levels were ignored to avoid errors.

After identifying each computational cell as part of a zone, we measured various quantities relative to each zone in postprocessing, such as the length, the volume, the volume-averaged quantities, and the probability density functions over the jet. We also defined a proxy for the jet aspect ratio as (πl³/V)^1/2, where l and V are the total length and volume of the jet, respectively. This is equivalent to the ratio l/R_j, eff, where R_j, eff is the radius of a cylindrical jet of length l and volume V.

3. Results

The results we present below are, to the best of our knowledge, the first 3D simulations of jets in HMMQs that are relativistic and include radiative cooling in parameterized form. They cover the evolution of the jet from its launch over the onset of instabilities and radiative cooling to the turbulent phase at the end of our simulations.

More specifically, we discuss the propagation of the jet through the stellar wind from its outburst close to the BH up to scales of about 6 × 10¹³ cm for Cygnus X-1. This corresponds to about 20 times the separation between the two stellar components d_orb, and 2 × 10¹³ cm ≈75 d_orb for Cygnus X-3 (the values for d_orb are consigned in Table 1). A brief overview of previous studies is given in Sect. 3.1. Details about our fiducial simulations are given in Sect. 3.2, with particular focus on the development of the KHI and its role in different phases of jet propagation (Sect. 3.2.1), as well as cocoon evolution (Sect. 3.2.2). Then, the impact of radiative losses on jet structure and dynamics (Sect. 3.3) is investigated, before we perform a short parameter study (Sect. 3.4) of the jet temperature (Sect. 3.4.1), kinetic power (Sect. 3.4.2), and stellar wind (Sect. 3.4.3).

3.1. Previous studies

The interaction of a jet with its ambient medium typically results in a richly structured flow field. For any discussion of this complex flow field, it is useful to resort to its basic idealized morphology, which is characterized by three surfaces that separate the flow into four zones. The innermost zone, the jet beam or spine, consists of unprocessed jet material. A combination of discontinuities (including a reverse or terminal shock and potentially shear layers, or also reconfinement shocks) separates the jet beam from the inner cocoon, which is composed of shocked jet material (also called “shear layer”, “jet sheath”, or simply “cocoon”). A contact discontinuity or working surface marks the transition from inner to outer cocoon (or “cavity”), the latter containing shocked ambient medium. A third discontinuity, a bow or forward shock at the jet head, marks the transition to the ambient medium (good sketches of this structure can be found in other works, e.g., Matsumoto & Masada 2019, for a recent example).

Martí et al. (1997) identified five parameters to completely specify a relativistic jet propagating into a homogeneous medium: the density ratio η ≡ ρ_j/ρ_w, the pressure ratio K ≡ p_j/p_w, the beam flow velocity v_j (or its associated Lorentz number γ_b), the beam Mach number M_j, and the polytropic index Γ.

Quantities with subscript j are relative to the values at injection, while the subscript w denotes quantities of the (in our case, wind-dominated) ambient medium. As pointed out by Martí et al. (1997), the propagation efficiency of a relativistic jet is mostly determined by the inertial mass density ξ = γ²ρh introduced in Sect. 2.1, and especially by the ratio η^* between beam and ambient medium inertial mass densities because the latter determine the momentum balance at the contact discontinuity in the jet head.

Sufficiently light (η < 0.1) supersonic jets, as considered in this work, display extended cocoons as the high pressure of the shocked gas drives a backflow toward the source. These jets also display a series of internal oblique shocks in the beam, whose strength and spacing are determined by the Mach number and the gradient in the pressure external to the beam (see, e.g., Gómez et al. 1995, 1997, and references therein): the higher these numbers, the stronger and closer to each other these oblique shocks. Increasing the Mach number also intensifies the expansion of the cocoon.

The structure of this cocoon is determined by the adiabatic index Γ: for models with Γ = 5/3, the cocoon is stable at first, but eventually evolves into vortices, producing turbulent structures and generating perturbations at the beam boundary. This enriches the internal structure of the beam. For models with Γ = 4/3, the first internal conical shock is strong enough for the resulting beam collimation to accelerate the flow. During this acceleration phase, the beam gas is less efficiently redirected to the cocoon downstream, accumulating at jet head. Once this acceleration is over, the continuous flow reestablishes itself, forming small turbulent vortices in the cocoon. The cocoon structure reflects this history, the parts of it formed before and after this beam acceleration phase presenting different morphologies. See Martí et al. (1997) for a more in-depth discussion.

Bodo et al. (1994) observed that during its propagation, a jet will present different structures tightly linked to the evolution of the KHI modes, identifying three phases: in a first “linear phase” the various modes grow following the linear behaviour, and no shock is present in the beam. The apparition of biconical shocks centred on the beam axis marks the beginning of the “expansion phase”, during which the strength of these shocks grow and the jet radius expands in the post-shock region. Finally, the evolution of the shocks leads to mixing between the jet and the external material, marking the start of the “mixing phase”.

The stellar wind can influence the propagation of a relativistic jet greatly: jet disruption by a constant and perpendicular wind were observed in Perucho & Bosch-Ramon (2008) and Perucho et al. (2010a) even for moderate jet kinetic luminosities. In particular, the presence of a wind strengthens the initial recollimation shock in the beam which can also strengthen jet asymmetric KHI produced at the wind/jet contact discontinuity. These may in turn disrupt the jet.

3.2. Cygnus X-1 and Cygnus X-3 fiducial cases

We start with a description of the fiducial cases for Cygnus X-1 and Cygnus X-3, respectively, against which all other sensitivity studies will be compared later on. Converting the numerical values of the parameters to the dimensionless quantities introduced in Sect. 3.1 (Tables D.2 and D.4) places our jets in the supersonic case with extended, turbulent cocoons and a beam with rich internal structure. Our fiducial runs indeed follow these expectations: the evolution of the jet is shown Figs. 3 and 4 for the CygX1 run, Figs. 5 and 6 for CygX3. Several features catch the eye, which we further elaborate on below. First, there is qualitative change in the appearance of the jet, from an early ‘well ordered’ state to a turbulent state later on. This change is also reflected in the propagation of the jet head and three phases of the jet evolution can be identified. Second, the aspect ratio of the jet is different for Cygnus X-1 and Cygnus X-3. Third, jet bending due to the lateral wind impact is observed in all simulations. Fourth, the jet is asymmetric due to the wind of the companion star.

Fig. 3. Rest-mass density slices of run CygX1 at times (top to bottom) t = 2000, 6000, and 12 000 s, showing the three evolutionary phases detailed in the text in Sect. 3.2.

Fig. 4. Dynamical and structural evolution of fiducial run CygX1. Left: position and speed of the jet head for fiducial run CygX1. Right: jet volume (beam and inner and outer cocoon) and aspect ratio. The limits of each evolutionary phase are marked by the vertical dash-dotted lines in the various panels.

Fig. 5. Rest-mass density slices of run CygX3 at times (top to bottom) t = 400, 1200, and 2500 s, showing the three evolutionary phases detailed in the text in Sect. 3.2.

Fig. 6. Dynamical and structural evolution of fiducial run CygX3. Left: position and speed of the jet head for fiducial run CygX3. Right: jet volume (beam and inner and outer cocoon) and aspect ratio. The limits of each evolutionary phase are marked by the vertical dash-dotted lines in the various panels.

3.2.1. Instability growth and phases of jet propagation

The structure of our jets goes through three phases, common to both Cygnus X-1 and Cygnus X-3 runs. We refer to these three evolution stages as the “smooth”, “instability growth” and “turbulent” phases respectively according to their inner structure. An illustration via rest-mass density slices is given in Figs. 3 and 5. The three phases also leave an imprint on the time-series data shown in Figs. 4 and 6.

During the first phase, the beam flow is surrounded by a smooth cocoon that is symmetrical at its head. A few internal shocks are present in the beam, starting with a strong recollimation shock situated at a few ∼10¹² cm downstream of the injection. Its existence and position are coherent with the criterion and analytical prediction from Yoon et al. (2016) obtained by equating wind ram pressure and lateral ram pressure in the beam. The aspect ratio of the jet defined in Sect. 2.7 increases gradually. It roughly follows a power law in time in the early propagation phase with an exponent of ∼0.6. For Cygnus X-3, the aspect ratio in the same phase first decreases before it increases with what could be a power law with a similar exponent as for CygX1. This may be explained by the strong asymmetry of the cocoon during this phase due to the strong stellar winds. A small deviation of the CygX3 beam can already be observed at this point. The jet head position and velocity follow the theoretical 1D result from Appendix A.1. Deviations indicative of the transition from phase one (smooth) to phase two (instability growth) occur after roughly 2000 s in the case of CygX1 and much earlier, after a few hundred seconds, in the case of CygX3. In particular, the speed diagram for CygX3 breaks almost immediately from the theoretical profile. This may be a consequence of the already existing bending of the beam.

In the second phase, instabilities grow in the jet, which perturb the flow in both inner cocoon and beam head. While the jet volume tends to grow faster now than during the first phase, the growth of the aspect ratio slows down with a ∼0.4 exponent for CygX1 case, and the jet head velocity overall decreases while the position breaks from the theoretical values. In CygX1, the number of over- and underpressure regions in the beam (see Figs. E.1 and E.2) stays approximately constant before increasing after about 5000 s. Ultimately, the growing instabilities cause oscillations of the beam head perpendicular to its propagation direction. For CygX3, these oscillations induce speed fluctuations even though the beam still retains its structure.

In the last phase, after about 6000 s in CygX1 and 2000 s in CygX3, the perturbations have reached the beam core. They modify the beam structure at the jet head severely, while the inner cocoon has become turbulent. This also marks a change in shape of both the cocoon and the jet head. The modification of the jet head shape can be linked to oscillations of the beam region that end in the reverse shock, as a beam head that is misaligned with the general jet propagation direction causes beam material to flow at higher speed and in same direction as the cocoon expansion, which deforming the jet. The jet head position evolves with an approximately constant mean velocity. Fluctuations up to roughly 30% are visible in the speed plots, which is in line with the persisting motion of the jet head position perpendicular to the jet axis. The volume of the outer cocoon evolves roughly as a power law in time with an exponent of about three for CygX3 and half as much for CygX1. The volume of the beam features a similar time dependence in the case of CygX1, but a shallower dependence for CygX3, with a powe- law exponent of about two instead of three. The aspect ratio decreases somewhat before becoming constant, at least in the case of CygX1.

This classification can be compared to the classification from Bodo et al. (1994) given in Sect. 3.1, but their “linear” phase escapes our data analysis because we dump data frames only every 100 s and 25 s for CygX1 and CygX3, respectively, while the estimated KHI linear growth timescale is typically on the order of a few tens of seconds for Cygnus X-1 and a few seconds for Cygnus X-3 (see Table 2). No value was found with this method for run CygX1_mP, where the beam is heavily disrupted by the stellar wind and the approximations made are no longer valid. Our first two phases (smooth and instability growth) appear to be subdivisions of their “expansion” phase, while our turbulent and their “mixing” phase match.

These phases are also visible in the speed diagram of the jet head that display the same trend for the Cygnus X-1 and Cygnus X-3 fiducial runs in Figs. 4 and 6: we can link the smooth phase with the initial acceleration, the deceleration and the concave part with the instability growth, followed by the turbulent phase. The first two of these three phases are of interest in the context of dedicated studies on instability onset and growth. Although a large body of associated literature exists, we are not aware of any such studies for relativistic jets in HMMQs including radiative cooling.

The link between the internal structure and the dynamic was discussed in Martí et al. (2016), suggesting that the growth of the KHI is related to the strength of these oblique internal shocks inside the beam: the KHI grows as the sound wave travels back and forth between the beam surface and the contact discontinuity, therefore more and stronger internal shocks produce a greater number of reflections within a given time or distance. This ultimately accelerates the growth of the KHI. The ripple-like structures observed in the cocoon, similar to pressure perturbations in Perucho et al. (2004b), could be viewed as markers of such sound waves.

Returning to the density in Figs. 3 and 5, although the beam and cocoon mix together at the jet head, the flow is not slowed down until the very end of the jet. The jets are bent away from the star almost as soon as the jet is established for the Cygnus X-3 runs, but also in Cygnus X-1 runs after enough lifetime of the jet. This bending angle ψ, defined in Yoon & Heinz (2015) as the angle between the local and the initial velocity vector, can be compared to the analytical value derived in the same paper. For run CygX1, we find ψ = 0.1 rad for a beam end at x = 6.8 × 10¹³ cm, which is close to ψ = 0.09 found analytically. For CygX3, we find ψ = 0.04 rad for a beam end at x = 1.52 × 10¹³ cm and 0.03 analytically, showing good agreement of our runs with the analytical estimate.

3.2.2. Cocoon evolution and radiative losses

Over the course of the initial jet outburst, the outer cocoon expands in the directions perpendicular to the jet propagation due to its overpressure compared to the ambient medium. On the windward side closest to injection, the interface between cocoon and wind is a bow shock and its dynamic are determined at first order by the balance between wind ram pressure and internal thermal pressure of the cocoon: depending on this balance, the interface will move either away from or closer to the beam. Farther away from the plane of orbit, the wind ram pressure becomes negligible and the interface dynamics is driven by the balance between internal and external thermal pressure. On the leeward side, the cocoon expands in the same direction as the wind speed. The resulting asymmetry of the cocoon is apparent at early times in both CygX1 and CygX3. At later times, the difference between windward and leeward side diminishes as the wind speed is increasingly aligned with the propagation direction of the jet (Figs. 3 and 5).

As the cocoon cools down with time either adiabatically due to expansion and/or from radiative cooling, the thermal pressure of the cocoon diminishes, which increases the influence of the wind on its dynamics. Figure 7 displays the volumic power losses per jet zone per process for the fiducial runs CygX1 and CygX3, measured over all the jet cells for two data points per evolutionary phase. In both cases, free-free losses dominate the cooling in the beam and inner cocoon, with a stronger cooling in the beam than in the inner cocoon. The colder outer cocoon is dominated by the very efficient line recombination cooling. This result holds true for all our simulated runs.

Fig. 7. Time evolution of radiative losses for the fiducial runs CygX1 (top) and CygX3 (bottom), derived at each cell and summed per zone, at time steps t = (500, 2000, 4000, 6000, 10 000, 15 000) s for CygX1 and t = (250, 500, 750, 1500, 3000, 6000) s for CygX3 (two data points per evolutionary phase). In both cases, free-free losses dominate cooling in the inner cocoon and beam, while line cooling dominates the outer cocoon.

Moreover, the gas in the cocoon has a velocity component in the positive x-direction. Thus the cocoon moves outward of the system with the beam, but at a slower pace. Ultimately, no trace of the original cocoon is left in the innermost parts of the jet. A thin interface of shocked stellar wind only a few r_b wide has instead formed between the wind and the beam (late times in Figs. 3 and 5). This “naked beam” is of interest because it represents a (quasi-) stationary state structure studied in the literature (e.g., Wilson 1987; Komissarov et al. 2015 for hydrodynamical jets, Martí et al. 2016; Bodo & Tavecchio 2018 for MHD jets) and can be related to direct observations.

3.3. Effects of losses on jet structure and dynamics

Cooling times in our two fiducial cases are such that over the time covered by our simulations, radiative losses have no significant impact on CygX1, and a clear influence on CygX3. The comparison of the position and speed diagrams of Figs. 8 and 9 shows that the first cooling effects are seen after 7600 s for Cygnus X-1, and the first cooling effects are visible after 100 s for Cygnus X-3. The jet head velocity evolves remarkably similarly with and without radiative losses during the first two phases of jet evolution (Fig. 8, bottom). In the case of CygX3, by contrast, radiative losses lead to a loss of much of the outer cocoon at distances of a few 10¹² cm (Fig. 9, top) and slow the jet head dwon (Fig. 9, bottom). Therefore, we restrict the following discussion mostly to CygX3.

Fig. 8. Effects of losses on run CygX1 structure and dynamics. Top: temperature slices of runs CygX1 and CygX1_noLoss at time t = 15 000 s. Both jets display similar structures, with the exception of a slightly larger outer cocoon at the head of the non-cooled jet. Bottom: jet head propagation and speed of the same runs, CygX1 in red and CygX1_noLoss in blue. The theoretical 1D propagation from Appendix A.1 is drawn as dotted lines following the same color-coding. The propagation is identical in both runs until the start of the turbulent phase, after which speed fluctuations differ, but the average propagation speed is identical in the two runs with almost no difference in the jet head position plot.

Fig. 9. Effects of losses on run CygX3 structure and dynamics. Top: temperature slices of runs CygX3 and CygX3_noLoss at time t = 9000 s. Two main differences appear: (1) The cooled beam is thinner. Its envelope closely follows the internal shocks structure in contrast to the non-cooled case. (2) The cocoon of the non-cooled jet expands farther at its basis, almost wrapping around the star, whereas in the cooled case, the cocoon has almost disappeared because ambient material has cooled enough to be blown back by the wind. Bottom: same as Fig. 8. The cooled jet (in red) is initially faster, but leaves the smooth phase earlier, after which point it is slower on average, as seen in the propagation plot.

3.3.1. Beam destabilizing effect through cocoon pressure

The addition of the loss terms has a destabilizing effect on the beam through its interaction with the inner cocoon: free-free cooling, shown in Fig. 7 to be the dominant process in both the beam and the inner cocoon, diminishes pressure in the jet with a different intensity depending on the jet zone: the beam cools faster than the inner cocoon, causing a stronger pressure gradient between inner cocoon and beam. This strengthens the oblique internal shocks (Fig. 10, left panels), which in turn accelerates the growth of KHI, as detailed Appendix A.2. Thus KHI grows faster in the cooled case, changing the dynamical behavior of the cooled jet, as shown in Fig. 9. Pressure in each zone is derived from the mean rest mass density and temperature measured over the corresponding marked cells defined in Sect. 2.7.

Fig. 10. Effect of radiative losses on the CygX3 beam structure. Top: ratio of volume-averaged pressure of the inner cocoon and beam for runs CygX3 (red) and CygX3_noLoss (blue). The overpressure is always higher in the cooled case. Bottom: tracer density at t = 750 s for runs CygX3 and CygX3_noLoss. The cooled jet features stronger oblique shocks.

The cocoon-to-beam pressure ratio for runs CygX3 and CygX3_noLoss is shown top panel of Fig. 10. It is to be noted that this ratio is greater than 1 at all time, ensuring a pressure collimation of the beam. In the non-cooled case, after the initial decrease that is caused by the settling in of the jet structure, the overpressure grows in two phases with a transition around ∼2000 s. In the cooled case, the growth rate seems constant from the start with the exception of a strong increase starting at t = 650 s, peaking at ∼1000 s, and joining the overall trend at 1350 s. This occurs as the jet transitions from the instability growth to the turbulent phase: t = 600 s indeed marks the apparition of strong oblique shocks in the cooled case. This stronger overpressure explains the difference in beam structure that is shown in the bottom panel of Fig. 10, displaying a longitudinal slice of the jet material tracer at t = 750 s including the plane containing the stellar center: a higher gradient of inner cocoon to beam pressure causes the stronger oblique shocks in CygX3 runs.

3.3.2. Effects on the outer cocoon expansion

Volumes of individual jet zones are affected by radiative cooling in different ways (Fig. 11). As explained in Sect. 3.2.2, the dynamics of the wind-cocoon interface near injection zone is mostly controlled by the inner thermal pressure of the outer cocoon. Therefore, the more efficient the cooling, the smaller the outer cocoon and the faster the evolution of the cocoon up to the “naked beam” situation. In Cygnus X-3 case, the effect on the cocoon is shown in the volume diagram in Fig. 11: the outer cocoon very quickly evolves to be consistently of greater volume in the non-cooled case. In contrast, the volumes of the two inner cocoons are similar in the smooth and turbulent phases: in the first phase, the cooling effects on the inner cocoon are negligible because it is both the hottest and least dense part of the jet, while in the second phase, the instability-induced turbulences dominate the cocoon flow. The period during which the volume of the inner cocoon differs is likely due to the different starting time of the mixing phase between these runs, as explained in Sect. 3.3.1.

Fig. 11. Volumes of the different jet zones as a function of time are affected differently by radiative cooling. The outer cocoon is larger without losses, while the volumes are comparable for the inner cocoon. The power-law dependence on time of the cocoon is robust for the cocoon, but not for the beam.

When the turbulent phase is reached, the inner and outer cocoon volumes show a similar power law dependence on time of roughly t³, regardless of whether radiative losses are included. By contrast, the beam volume displays a different power-law dependence in the cooling and non-cooling case. The relative volume of outer to inner cocoon is much larger in the no-loss case than in the loss case. This may be an issue if radiative losses are diagnosed only during postprocessing from adiabatic solutions.

The cocoon form is also affected: The comparison of runs CygX3 and CygX3_noLoss in the top panel of Fig. 9 shows that the expansion is strong enough in the non-cooled case to cause the cocoon to almost wrap around the star before it is blown back as the cocoon pressure diminishes, while in the cooled case, the cocoon is almost immediately blown back to a thin shell by the strong stellar winds.

3.4. Parameter sensitivities

We start with sensitivities to the assumed beam temperature (Sect. 3.4.1), which has been covered comparatively little in the literature and is thus somewhat more extensively dealt with here. Sensitivities to beam power and wind parameters follow (Sects. 3.4.2 and 3.4.3).

3.4.1. Effects of the jet temperature on instabilities growth

An increase in the jet injection temperature T_j lowers the beam Mach number. We expect the jet to display a smaller cocoon and to be less stable as the distance between internal shocks in the beam diminishes with it. This is true for the step from 10⁸ to 10⁹ K, but not for the step from 10⁷ to 10⁸ K. We ascribe this difference to the action of the first recollimation shock, which heats the beam of CygX1_T7 to similar values as are found in CygX1. Increasing T_j also results in a higher overpressure between the inner cocoon and the beam, which further strengthens the effects of oblique shocks we described above and leads to a faster destabilization of the beam, as explained in Sect. 3.2.1.

This destabilization is visible in the speed diagrams in Fig. 12, showing colder jets to be more stable than hot jets: CygX1_T7 shows similar dynamics as the fiducial run with a longer instability growth phase. In contrast, the run CygX1_T9 displays different dynamics in this phase from the other two: the jet propagation speed slows down to a plateau instead of exhibiting a progressive acceleration. This difference in dynamical regime can be linked to the mean beam temperature in the top panel of Fig. 13: beams associated with runs with T_j = 10⁷ and 10⁸ K display almost the same temperature as a result of the heating at the initial recollimation shock, which raises them to a few 10⁹ K with very little differences. They begin to deviate from each other around the same time as the jet propagation speed does. With T_j = 10⁹ K, the temperature upstream of the shock is about the same as the downstream temperature, resulting in an higher effective beam temperature. This in turn means that the internal shocks that are closer to each other. This observation is confirmed by drawing the probability density function (PDF) of the temperature in Fig. 14. At t = 4000 s (full lines, during the instability growth phase), the PDFs for CygX1 and CygX1_T7 are identical, while CygX1_T9 differs for temperatures higher than ∼4 × 10⁹ K. At t = 10 000 s, the PDFs of the three runs differ by roughly the same amount, especially at the peak around 2 × 10¹⁰ K. This can be interpreted as showing that the turbulence in the cocoon distributes the available thermal energy and therefore makes the difference in injected temperature visible.

Fig. 12. Jet head propagation and speed for runs CygX1_T7 (red), CygX1 (blue) and CygX1_T9. Run CygX1 slows down to the turbulent phase earlier than CygX1_T7, but display the same average speed during the turbulent phase, as shown by the almost parallel propagation curves. Run CygX1_T9 displays a different behaviour, decelerating to a plateau in the instability growth phase and with a lower average speed than the other two runs, which appears to decelerate after t = 12 000 s.

Fig. 13. Comparison of temperature and overpressure of runs CygX1_T7, CygX1 and CygX1_T9, the color-coding is the same as in Fig. 12. Top: evolution of the zone-averaged temperature with simulation time. CygX1 and CygX1_T7 differ only slightly at first, but then start to evolve differently after the 5000 s mark. CygX1_T9 shows an overall higher temperature, with a small peak in inner cocoon temperature around t = 6000 s. Bottom: evolution of the ratio of the inner cocoon to beam pressure, as defined Sect. 3.3.1. CygX1 and CygX1_T7 present similar values up to t ∼ 5000 s, while CygX1_T9 presents a stronger overpressure.

Fig. 14. Comparison between runs CygX1 (red), CygX1_T7 (blue) and CygX1_T9 (magenta). Probability density function of the temperature at t = 4000 s (full lines) and t = 10 000 s (dash-dotted lines). In the early stages, CygX1 and CygX1_T7 are indistinguishable from each other, while CygX1_T9 is slightly hotter and more of its volume is hotter than ∼5 × 109 K. In the mixing phase, the temperature repartition is more in line with the injected temperature.

The bottom panel of Fig. 13 shows the pressure ratio between inner cocoon and beam for the same three runs, where the pressure is derived from the mean temperature and from the rest mass density obtained by averaging over the marked cells. The pressure ratio in run CygX1_T9 displays a different behavior from the other two runs, showing higher values as soon as t = 1000 s. This results in stronger internal shocks in the beam. The values for runs CygX1_T8 and T7 are similar up to time t ∼ 6000 s, however, meaning that the internal structure of these two jets are similar during this period. These two effects both accelerate the growth of KHI modes, which is confirmed by the derived values of t_KHI of 207, 71.4 and 26.8 s found in Table 2 for runs CygX1_T7, CygX1, and CygX1_T9, respectively.

For the turbulent phase, the effect of the different beam temperatures is weaker. The velocity of the jet head is comparable to within its fluctuation range, except possibly for the very late time that is still covered by our simulations when the jet head velocity for the hottest jet seems to slow down slightly as compared to the two simulations with cooler jets. This would suggest that as soon as a more generic turbulent behavior takes over the dynamics, the beam injection temperature becomes less important.

3.4.2. Effects of the injected power

A decrease in the jet kinetic power decreases its propagation efficiency as well as its stability. Starting from our fiducial test cases CygX1 and CygX3, we reduced the jet power by a factor of 10 (CygX1_mP) and 2 (CygX3_mP) via reducing the jet density at constant beam speed, as detailed in Tables D.2 and D.4. In these modified settings, the jets are expected to propagate more slowly and to be more prone to instabilities because the inertial mass density is lower as long as the jet is not disrupted by the stellar wind, as pointed out in Perucho et al. (2010a). This occurs for run CygX1_mP, as shown in Fig. 16. At constant beam speed v_b, the amplitude of the speed variations along the trend defined in the beginning of this section and the timescales at which they occur are controlled by the injected kinetic power, but the trend itself is not affected.

Figure 15 shows propagation plots for runs CygX1 (red) and CygX1_mP (blue), with the same parameters except for ρ_b, which is ten times smaller in the second case. The weak jet displays a different behavior as the beam is strongly bent away by the wind and is then broken down by instabilities after t = 8000 s, as shown in red in Fig. 16, drawn at t = 10 000 s. The beam mixes with the cocoon farther away from the contact discontinuity, and the momentum flowing from the reverse shock at beam head is partly dissipated in the cocoon. This smoothes the effects of beam head oscillations on the jet head dynamics and also slows the jet propagation down.

Fig. 15. Sensitivity of the jet head position and speed to kinetic power and wind ram pressure for Cygnus X-1. The jet kinetic power is divided by 10 from CygX1 (red) to CygX1_mP (blue). In the weak case, the jet is slower, but the position fits the theoretical evaluation for a longer time. The speed diagram shows no oscillations. The wind speed is increased by 50% for run CygX1_wind (green), which results in a higher starting speed, a shorter reacceleration phase, and a weaker bump in the second deceleration phase. The average speed in the turbulent phase is weaker. The speed and position plots differ faster from the theoretical 1D values in the strong wind case.

Fig. 16. vx/vj slice for run CygX1_mP at time t = 10 000 s. The jet is heavily bent from the wind effects, and its velocity breaks down before it arrives at the head.

The propagation of runs CygX3 and CygX3_mP is shown in the left panel of Fig. 17. ρ_b is here divided by 2 between the two runs. In this case, the remark about the shape of the speed plot holds also: weaker jets show a similar evolution with smaller amplitudes in the global variations of the jet velocity. In contrast to the Cygnus X-1 case, however, a decrease in density destabilizes the jet: the first speed peak that occurs when the jet propagation breaks from the momentum balance model occurs earlier in the weak case, and the speed fluctuations begin earlier.

Fig. 17. Sensitivity of the jet head position and speed to kinetic power (left) and wind momentum (right) for the Cygnus X-3 runs. Top left: a division of the jet power by 2 decreases all dynamical values, starting from the mean propagation speed. In particular, the acceleration is far less efficient: it shows a gain of ∼40% in speed between the starting point and the peak against a gain of ∼120% in the CygX3 run. Top right: a reduction of Ṁ⋆ by 25% and vw by 33% (keeping η constant) lengthens the smooth phase and strengthens the speed fluctuations in the following phases. The jets settle to the same median speed in the turbulent phase. Bottom right: same reduction as in the top right panel, with half the jet power. In the CygX3_mPmW case (right, blue), the initial acceleration phase is followed by a speed plateau.

3.4.3. Wind effects on the jet propagation

An increase in the wind ram pressure shortens the instability growth phase. This shows that this deceleration and the wind impact on the beam are linked. In the Cygnus X-1 runs, the impact of a stellar wind speed that is 50% higher and therefore causes an increase of 50% in the wind ram pressure at a constant mass-loss rate is shown in Fig. 15 by comparing the CygX1 (red) and CygX1_wind (green) runs: the stronger the wind speed, the shorter the instability growth phase. The beginning of this phase also starts slightly earlier: at 1800 s for CygX1_wind versus 2200 s in the fiducial case. The ambient density, moreover, drops slightly from CygX1 to CygX1_wind because Ṁ^⋆ is constant between the two runs. This means a higher η and therefore an easier jet propagation through the medium, and it explains the difference in starting propagation speed between the two runs. The theoretical 1D estimates for position and speed deviate from the measured values earlier in the strong wind case because the multidimensional effects are stronger. A higher wind speed induces a small plateau in jet speed at the very beginning of its propagation, which cannot be modeled by our 1D theoretical estimate.

For the Cygnus X-3 runs, the right panel of Fig. 17 compares the fiducial run CygX3 (red) with run CygX3_mW (blue), in which the wind is slower. All other parameters were kept constant, with the exception of mass-loss rate, to ensure same η value between the runs and halving the wind ram pressure on the jet. The right panel shows the same modification using the weaker jet setup (runs CygX3_mP and CygX3_mPmW). In the first case, the initial accelerating phase is twice as long for run CygX3_mW than run CygX3, and the second phase in run CygX3_mW consists only of a global deceleration because no reacceleration is observable. After this deceleration, run CygX3_mW shows a slightly higher median propagation speed.

When we compare runs CygX3_mP and CygX3_mPmW in the bottom right panel of Fig. 17, a new effect arises: instead of a deceleration, the initial acceleration phase of run CygX3_mPmW is followed by a plateau. The velocity fluctuations around the trend also appear later, but with a greater amplitude when the wind is weaker. After this speed plateau, the jet decelerates to a median value similar to the strong wind case, while the speed fluctuation timescale also decreases to a similar value as the strong wind case. This occurs after t ∼ 3000 s, when the jet has propagated to a distance of ∼5 × 10¹² cm (=20 d_orb). At this point, the wind is almost colinear with the jet propagation, and its lateral ram pressure on the jet is negligible.

4. Discussion

We have simulated jet outbreak and propagation in a large spatial and temporal frame based on the physical model developed in Sects. 2.1 and 2.4. The simulation particularly included radiative cooling by different processes. We now critically discuss some underlying points of this model and the limitations of the results, and we describe the directions in which future work should be oriented. This point is of peculiar importance if we wish to compare simulation results to observations.

4.1. Assumption of a hydrodynamical framework

High-energy nonthermal photons are a prominent feature of microquasar observations. They are present in the low and high state, but in different forms. They are thought to be produced by nonthermal processes that occur in the jet (Molina & Bosch-Ramon 2018; Poutanen & Veledina 2014; Malzac et al. 2018), in which high-energy particles are accelerated to relativistic speeds and adopt a power-law spectrum. The mechanisms invoked to accelerate particles are stochastic acceleration (Fermi mechanism) at shocks, magnetic reconnection, or wave turbulence; see the recent reviews by Marcowith et al. (2020), Matthews et al. (2020). All these acceleration mechanisms imperatively demand the plasma to be collisionless to a high degree. The Coulomb-logarithm of a typical jet beam in microquasars is approximately 15 and even higher in large regions of the cocoon. Consequently, kinetic timescales and kinetic inertial length scales are more than 10 orders of magnitude smaller than hydrodynamical scales. In short, the mean free path of thermal particles of a jet may easily be as long as the dynamical length scale of the jet and its cocoon. This poses the question to which degree jets can be understood on the basis of a hydrodynamical model.

Attempts to develop kinetical models of jets have been pursued; see Nishikawa et al. (2021) for a recent review. These models provide good insight into the various microphysical processes, and show how instabilities such as the Kelvin-Helmholtz, Rayleigh-Taylor, and kink instabilities translate into the collisionless regime. They also show that jet-like structures can develop on kinetic scales. However, it remains an open question whether the results found on kinetic scales can be scaled up to hydrodynamical scales and thus to scales on which the objects are seen in the skies. Some first steps toward the connection between kinetic and macroscopic scales have recently been made: Dieckmann et al. (2019) showed that for even weak magnetic guide fields, kinetic jets can develop a structure resembling the hydrodynamical structure of a jet over scales much larger than the kinetic scales. An electromagnetic piston-like structure is coherently formed, acting like the contact discontinuity between jet and ambient material in a hydrodynamical jet (see also Dieckmann et al. 2017). We note that in the same context, Dieckmann et al. (2019) also showed that a leptonic jet propagating into an ambient ion-electron plasma can accelerate positrons to speeds of several hundred MeV. In this way, microquasars candidate sources can explain the positron population in the cosmic ray spectrum.

4.2. Electron temperature and the single-fluid assumption

In single-fluid simulations, the dynamics is dominated by the ions, while the electrons are responsible for most of the radiative losses. Studies such as Vink et al. (2015), Zhdankin et al. (2021) have shown that processes such as shocks and radiative relativistic turbulences may create and/or maintain a difference between ion and electron temperature. This temperature difference has been observed in supernova remnants (Vink et al. 2003). In particular, Vink et al. (2015) suggested from thermodynamical consideration at shocks that the sonic Mach number of the upstream flow is the main parameter governing this temperature difference: at low Mach numbers (M ≤ 2), the shock will not induce a temperature difference, while at a high Mach number above M ∼ 60, the ratio of the electron and ion temperature will be equal to the mass ratio T_e/T_i = m_e/⟨m_i⟩. Between these extremal values, this ratio will roughly vary with M⁻².

In relativistic jets, at least two strong shocks present an upstream Mach number strong enough to induce a difference between ion and electron temperature downstream of it: the recollimation shock after injection, and the reverse shock at the beam end. To ensure that the calculated energy losses are valid, we need to determine the time taken by the electrons to thermalize to the flow temperature. Using formulas from Trubnikov (1965) found in Huba (2016), we obtained an estimate for the thermalization timescale considering the extremal case for a high Mach number. Table 3 presents the rest mass density and temperature values after the recollimation shock and reverse shock in CygX1 at t = 5000 s and CygX3 at t = 1000 s, the derived thermalization timescale, and the associated characteristic length using the local flow speed value. In all cases except downstream of the CygX3 reverse shock, the associated length scale is shorter than the grid resolution, meaning that the electrons would have thermalized to the flow temperature after a single simulation grid downstream of the shock. In the last case, we can estimate the volume within which an error is made using the flow temperature to derive losses as , where c_s is the local sound speed. With v_h = 3 × 10⁹ cm s⁻¹ at t = 1000 s, we obtain V_therm = 3.3 × 10³³ cm³, which accounts for 0.37% of the inner cocoon volume at that time. We then consider our T_e = T_i = T approximation to be valid in our hydrodynamical framework. In a more realistic framework in which the plasma is collisionless (see the discussion in Sect. 4.1), Coulomb interactions cannot equilibrate the electron and ion temperatures, but other processes such as plasma instabilities may fill that role.

4.3. Thermal diffusion

Another important process that in not included in our numerical model is thermal diffusion by relativistic thermal and nonthermal electrons or X-ray photons. This process is likely important as it creates hot shock-precursors, decreases peak temperatures, smoothes out contact interfaces, and enhances cooling. These features are suggested when looking at work which includes the process in the context of supernova remnants (Chevalier 1975; Tilley et al. 2006) and colliding winds in binary systems (Myasnikov & Zhekov 1998; Motamen et al. 1999). Including heat transfer in a simulation is numerically demanding as it demands (semi-) implicite solvers due to the very stiffness of the process. Nevertheless, a few attempts have been made to develop performing solvers (Balsara et al. 2008; Palenzuela et al. 2009; Viallet et al. 2011; Kupka et al. 2012; Commerçon et al. 2014; Bucciantini & Del Zanna 2013; Dubois & Commerçon 2016; Viallet et al. 2016). Future work may expand these attempts to jet simulations.

4.4. Wind structure and jet bending

In high-mass microquasars, jets are launched into winds originating from orbiting and rotating stars. This causes a circumbinary environment structured in Archimedean spirals in colliding-wind binaries (Walder 1995) and accreting binaries (Walder et al. 2008, 2014). These spirals, also called corotating interaction regions (CIRs), are bound by shocks that confine high-density, high-temperature regions. These structures will likely have an impact on the jet propagation and its stability, but will also likely lead to flashing outbursts in emission.

The bending of the jet away from the mass-shedding donor star, investigated, for example, by Perucho & Bosch-Ramon (2008) and subsequent papers, Yoon & Heinz (2015), Bosch-Ramon & Barkov (2016), and again confirmed in this paper, will lead to helical jet trajectories. In the most extreme cases, an observer will see the jet under different angles during an orbit of the system (see, e.g., Horton et al. 2020). Future works will lead to a more quantitative statement of this prediction. As a side note, we showed in Appendix C.4 that jets performed with a relativistic scheme display a higher propagation speed, even in the mildly relativistic case. Because the jet bending depends on the transverse momentum accumulated by the jet during its propagation, we would expect a weaker bending in relativistic simulations than in Newtonian simulations that are performed with the same parameters.

Another unknown is wind clumping. We know that winds from hot massive stars are clumped (e.g., Oskinova et al. 2012), with density contrasts ⟨ρ²⟩/⟨ρ⟩² ≈ 2. Poutanen et al. (2008) suggested that this is likely the case for Cygnus X-1 based on dips in the X-ray light curve (Bałucińska-Church et al. 2000; Miskovicová et al. 2016; Hirsch et al. 2019). It is unknown, however, where these clumps are formed: in the stellar atmosphere, or farther away from the star, for instance, by supersonic turbulence in the wind. Whether these clumps are small and dense, an intermittent density fluctuation of compressible turbulence, or a mix of both is unknown as well; see the discussion in Walder & Folini (2003). Perucho & Bosch-Ramon (2012) suggested that even moderate wind clumping has strong effects on jet disruption, mass loading, bending, and likely energy dissipation. de la Cita et al. (2017) suggested that the standing shocks introduced in the jet flow by its interaction with a clumpy wind would generate a higher apparent gamma-ray luminosity through inverse Compton scattering of the stellar photons, as well as efficient synchrotron cooling. In light of these results, we may expect that introducing clumpiness in the wind might significantly enhance the radiative cooling of the jet and modify its dynamics even more strongly.

4.5. Wind driving and dynamics

Velocity and density profiles in winds shed from hot massive stars in binaries are prone to a large uncertainty. This uncertainty pushed us to choose a grid of fixed wind parameters for the study presented in this paper and to make no attempts to model the acceleration of the wind. Winds from massive stars are driven by the radiative pressure on free electrons (to about one-third of the driving force) as well as the scattering of stellar photons in millions of UV and optical lines. This wind driving is expected to be vastly different in supergiants (as in Cygnus X-1) and WR stars (as in Cygnus X-3), in that the winds from WR stars are optically thick in the subsonic acceleration phase, whereas they are relatively optically thin in supergiant winds.

These winds are well understood for single stars, with the exception of the subsonic phase of WR winds; see the review by Kudritzki & Puls (2000). In binaries, the situation is more complex. Hainich et al. (2020) presented an observation-based study of winds in HMXRBs and basically confirmed that the wind driving from stars in those binaries is the same as for single stars of the same type. However, the situation is more complicated by the secondary radiation source, either another star or a compact object and its environment. This second source leads to an ionization structure within the wind that in contrast to the single-star situation does not fit the frequencies of the photons emitted by the star that produces the wind.

This in particular affects the region between the wind-shedding star and its companion. There, the driving by photon scattering in spectral lines can be suppressed and the wind acceleration may be weaker or even inhibited compared to the winds of a single star (Stevens & Kallman 1990; Stevens 1991; Blondin 1994). Models that take this inhibition into account can be fit to observations (Gies et al. 2008). The case of Cygnus X-1 was investigated through complex radiation hydrodynamics simulations by Hadrava & Čechura (2012) and Krtička et al. (2015), who suggested that strong inhibition of the wind speed might occur at the BH location. In this region, wind density and velocity profiles in binaries are very difficult to access observationally and are largely unknown. For Cygnus X-3, the situation is even more unknown because the wind is still optically thick at the BH location. Studies such as those by Vilhu et al. (2021) find a strong wind-suppressing effect in particular in the extreme-UV region. However, we can firmly state that a wind is present in both systems, based on the modulation of the X-ray light curve over the binary orbit, for instance. This modulation is due to the different optical paths that X-ray photons have to travel through the wind and thus to the different attenuation of the X-ray light through absorption in the wind (Bonnet-Bidaud & van der Klis 1981; Poutanen et al. 2008; Grinberg et al. 2015). Finally, the effect of the BH gravitational field must be considered, which may focus the stellar wind and modify the density and velocity profile that are encountered by the propagating jet, but might also modify the evolution of the jet cocoon at its base.

Another important open problem is the influence of magnetic fields on the wind because strong fields could substantially influence the mass-shedding process, for example, by a strong enhancement of the mass loss in the equatorial plane (ud-Doula et al. 2006; Ud-Doula et al. 2008; Bard & Townsend 2016). A general discussion of magnetic fields of massive stars and how they influence their environment can be found in Walder et al. (2012).

4.6. Nonthermal components and radiative processes emission

Our losses are modeled as optically thin at every frequencies. This may not be the case, especially in UV, optical, and radio ranges at which our thermal electrons radiates. Reheating due to photon absorption (Belmont et al. 2008) or synchrotron self-absorption as well as effects such as synchrotron self-Compton were not investigated and may have a strong impact on the jet temperature and therefore structure and dynamics. This may not hold for Cygnus X-1, for which cooling and its effects are moderate but it might play a strong role in Cygnus X-3 case, especially because winds from WR stars are optically thick over higher frequency bands than winds from O-type stars. Moreover, the temperatures attained in our simulations exceed the electron pair-production threshold and reach a regime in which proton distribution becomes relativistic, which will impact the radiative losses.

Another point we did not consider is the particle acceleration and the nonthermal processes: particles can be accelerated at magnetic reconnections (Giannios 2010; Giannios et al. 2009; McKinney & Uzdensky 2012; Melzani et al. 2014) and shocks (Araudo et al. 2009; Bordas et al. 2009; de la Cita et al. 2017) or by stochastic interaction with magnetized turbulence. Shearing flow acceleration as also been invoked as potential injection process of relativistic particles (Rieger & Duffy 2019). An extensive review of these processes is given in Marcowith et al. (2020). As all processes will significantly contribute to the nonthermal emission, they can inject a substantial fraction of the kinetic jet energy into nonthermal components and probably change the dynamics because part of the gas energy will be dissipated by the accelerated particles. Alternatively, the pressure imparted to these particles can directly modify the flow dynamics.

4.7. Jet outburst dynamics and emerging structures

The dynamics of the second (instability) phase is not well understood, and neither is the physical phenomenon that sets the fluctuation speeds. The physics of instability growth is deep and rich, and we may have overlooked some effects, in particular, concerning the nonlinear phase. Another limitation of this study is our grid resolution, which may not accurately capture all the shocks that form in such a system or the growth of small-amplitude modes, for which the impact was highlighted here. Furthermore, the adiabatic index Γ was chosen to be constant and equal to 5/3, where a variable index such as the method suggested in Mignone et al. (2005) and Mignone & McKinney (2007) would more accurately reflect the relativistic gas EoS (Synge 1957) and lead to more realistic jet structures.

We included thermal radiative losses in our study of jets structure and dynamics. These effects were not included in similar works. No focus was placed on jet bending, but it is observable in the jets of Cygnus X-1 and Cygnus X-3, in agreement with similar studies such as Yoon & Heinz (2015). Moreover, the jets are injected with no inclination relative to the orbital plane, in contrast to works such as Dubus et al. (2010), who suggested that the inclination angle for Cygnus X-3 lies between 20 and 80°. Another strong point of this study is the simulation box size of 10¹⁴ cm, which is 33 times the orbital separation of Cygnus X-1 and close to 400 times that of Cygnus X-3. This allowed us to observe the outbreak and early dynamics details as well as the emergence of a steady “naked beam” structure for later timesteps that may be compared to similar works and translated into synthetic observational datas.

Last, our study is purely hydrodynamical and does not cover the potential stabilizing effects of the magnetic field as well as the development of MHD instabilities. The study of MHD jets is recent (e.g. Mizuno et al. 2015; Martí et al. 2016; Mukherjee et al. 2020) and the field will continue to grow as time passes. The influence of the stellar magnetic field may also lead to further modifications on the structure and dynamics of the jet.

5. Summary and conclusions

We performed large-scale 3D relativistic hydrodynamical simulations of jet outbreak and early propagation in HMMQs. We added radiative cooling effects in the energy equation using a Maxwell-Jüttner distribution for the electrons. Two fiducial cases inspired by the HMMQs Cygnus X-1 and Cygnus X-3 were considered, along with parameter sensitivity studies. The jets of the two systems are in the same domain of dimensionless parameters, but with different instability growth and cooling timescales (about ∼10² shorter in the latter case), enabling us to better highlight the importance of these processes in jet structure and dynamics. In particular, on the timescales covered by our simulations, radiative losses, mainly free-free mechanism, play a relevant role for Cygnus X-3, but not for Cygnus X-1.

We investigated the impact of radiative losses on jet structures and dynamics and showed their strengthening effect on KHI growth, which was a focal point of our analysis. We identified three main dynamical phases: (1) an initial self-similar propagation in line with 1D momentum balance arguments, followed by (2) a modification of the inner jet structure and a phase of instability growth, and finally, (3) a turbulent cocoon and destabilized beam. We explored the effects of losses on the dynamics in the second phase through their modification of the KHI growth rate as well as the impact on the jet structure, and we examined the impact of the various parameters.

In Cygnus X-3, where radiative losses are observed on the timescales covered by our simulations, these losses are found to affect the volume ratio of the outer to inner cocoon. This means that simple postprocessing of simulations that do not take radiative cooling into account to study emissions from HMMQ jets must be met with caution. Likewise, the beam volume is found to obey a different power law in time depending on whether or not cooling is active.

We find the sensitivities to jet power and wind parameters to be in line with the literature. We add to this picture by demonstrating that an increasing beam temperature results in a faster instability growth and a slower jet head velocity. In view of the wide variety of known HMMQs, but also given the difficult estimation of their system parameters, further and more extended parameter studies are clearly desirable in the future. Based on the results presented here, we would advocate that such studies not be limited to adiabatic simulations.

Acknowledgments

We acknowledge support from the French National Program for High Energies (PNHE). Simulations have been performed at the Pôle Scientifique de Modélisation Numérique (PSMN) hosted by the ENS de Lyon, and at the centers of the Grand Equipement National de Calcul Intensif (GENCI) under grant number A0090406960. We thank the referee for the helpful review and insights, which helped improve the paper significantly.

References

Appendix A: Propagation and instabilities

A.1. Model for jet propagation

A model for the propagation of a relativistic jet was derived, for instance, by Martí et al. (1997) and Mizuta et al. (2004), who neglected multidimensional effects, but assumed 1D momentum balance between the beam and the ambient gas in the rest frame of the contact discontinuity at the head of the jet. In the case of an ambient medium with its own (nonrelativistic) flow speed v_w, this gives

(A.1)

where γ are Lorentz factors, S are cross-section areas, and v are velocities. Subscripts b refer to beam values, w refers to the ambient medium, and h to the interface at jet head. Solving Eq. A.1 for v_h, we obtain

(A.2)

where , A is the cross-section ratio (A ≡ S_b/S_w), and K is the pressure ratio defined in Sect. 3.1. This equation is consistent with the one derived by Mizuta et al. (2004) when taking v_w = 0. The term including K can be neglected as the sound speed in the stellar wind is much slower than the beam velocity. Assuming A = 1, Eq. A.2 therefore becomes

(A.3)

A.2. Instabilities in relativistic jet flows

During the jet propagation, various hydrodynamical instabilities can be triggered and in the end perturb the beam, reducing the effective beam speed at the front shock and decelerating the jet. One of them is the KHI at the relativistic flow interface, which was extensively studied: Turland & Scheuer (1976), Blandford & Pringle (1976), Ferrari et al. (1978), Hardee (1979), Bodo et al. (1994), Hanasz & Sol (1998), Hardee et al. (1998), Hardee et al. (2001), Perucho et al. (2004b), Perucho et al. (2004a), Perucho et al. (2005), Mizuno et al. (2007), Rossi et al. (2008), Perucho et al. (2010c).

Then, the Rayleigh-Taylor instability (RTI) can be triggered when a lighter fluid supports a heavier one against gravity, or equivalently, if the lighter fluid accelerates the heavier one (see Norman et al. 1982; Allen & Hughes 1984; Duffell & MacFadyen 2011 for nonrelativistic flows and Matsumoto et al. 2017 for relativistic flows). This impacts the jet structures stability when the jet expands radially due to the centrifugal force, as studied by Meliani & Keppens (2007), Meliani & Keppens (2009), Millas et al. (2017), for example, or oscillates radially because of a pressure gradient, see Matsumoto & Masada (2013), Toma et al. (2017).

Other instabilities include the Richtmyer-Meshkov instability (RMI); see Nishihara et al. (2010) for a review and Matsumoto & Masada (2013) for a numerical study. The centrifugal instability (CFI) was studied in Gourgouliatos & Komissarov (2018b). Finally, we cite vortex formation and shedding at the contact discontinuity of jet head, as in Norman et al. (1982), Scheck et al. (2002), or Mizuta et al. (2004), for instance, which increases the cross-section of the jet head. This therefore diminishes the jet propagation speed. The vortices may also perturb the beam flow when they encounter it while propagating in the inner cocoon. Norman et al. (1982) noted the importance of the beam Mach number in the vortex-shedding phenomenon, observing that lower Mach numbers produce higher vortices. This work will focus specifically on the KHI as other instabilities develop at the beam radius scale, meaning that most of their modes will not appear due to grid resolution limitations.

A.2.1. Kelvin-Helmholtz instability

To numerically derive the linear growth time of the Kelvin-Helmholtz instability t_KHI, we chose to use the approach introduced in Hanasz & Sol (1996), describing the jet with a 2D slab geometry. This choice was made because the solutions for the slab and cylindrical geometries behave in a similar way in a wide range of physical parameters, with only slight numerical differences (Ferrari et al. 1982). In this analogy, symmetric and antisymmetric modes in the slab correspond to pinching and helical modes in the cylinder. Only high-order fluting modes do not have counterparts in slab jets.

Hanasz & Sol (1996) considered a core-sheet jet made of three layers: a beam with a relativistic flow, a cocoon with a nonrelativistic flow speed, and an ambient medium at rest. The transition layers at all interfaces are described in the vortex-sheet approximation. In the following, quantities with subscript b refer to the beam, subscript c refers to the cocoon, and w refers to the ambient medium. We introduce the quantities η_c ≡ ρ_b/ρ_c and η_w ≡ ρ_c/ρ_w as the beam-cocoon and cocoon-ambient medium density contrast, respectively, and R ≡ r_c/r_b as the radius ratio of the sheet and the core.

The deformations of the internal interface generate a sound wave that travels in the cocoon layer with the same frequency and wavenumber. As the maxima of the growth rate coincide with acoustic waves that travel an integer n times their wavelength on their path back and forth in the cocoon, we chose to focus on these resonances to derive an estimate for the KHI growth time.

For a sound wave traveling between the internal and external interfaces with a wavelength and a propagation angle α defined by tanα = k_x/k_z, the distance traveled by the wave between the interfaces is L = (R − 1)/cosα. From the full derivation of the dispersion relation, we have in the cocoon

(A.4)

(A.5)

Substituting these expressions in the resonance condition 2L = nλ_c, we obtain the following expression:

(A.6)

This equation can then be solved using a Newton-Raphson method and leads to the calculation of the wavenumber k_x maximizing ω with the densities and radii as parameters. The linear growth time for the KHI is obtained from the corresponding ω by taking t_KHI = ω⁻¹. To derive the growth time from our simulations, densities were measured as a volume-averaged value over each jet zone (more on this method in Sect. 2.7), while the beam and inner cocoon radius were derived from the respective measured volume and length in the propagating direction of the zones by approximating them as coaxial cylinders. These timescales correspond to the “linear” growth times, whereas the observed growth rate in our simulations will be significantly higher due to nonlinear effects. They are still of interest when we compare them for different runs because the relation between different linear growth time is the same as for the observed runs.

A.2.2. Rayleigh-Taylor instability

Neglecting the effects of interface curvature and assuming the perturbation of physical variables have the Wentzel-Kramers-Brillouin spatial and temporal dependence in exp[i(ks−ωt)], where s is the local coordinate tangent to the interface and perpendicular to jet propagation, we can derive the linear growth time using the dispersion relation from Matsumoto et al. (2017),

(A.7)

We take γ_ic = 1 (which is verified in all our runs), and g is the acceleration of the surface assuming that the amplitude of the oscillation is roughly equal to the beam radius and p₀ is the pressure at the interface, we take p₀ = (p_b + p_ic)/2 as a first approximation. We deduce from Eqs. A.7 that the RTI timescale t_RTI = Imω⁻¹ is proportional to the square root of the wavelength λ = k⁻¹, meaning that the temporal growth of the shorter-wavelength modes is faster. However, when the instability to wavelengths smaller than the jet radius and in a numerical grid is neglected, only the modes with wavelengths of several cell length will matter. This restricts the growth of the RTI to only a few possible wavelengths. This explains the large difference to the growth time of the KHI: we find t_KHI to be mainly in the range 10^1 − 2s, while t_RTI is found to have values of a few 10⁴ to 10⁵s using this derivation.

Appendix B: Radiative processes

A thermal plasma distribution of the electrons can be written in a general way using a Maxwell-Jüttner distribution of the electrons in the plasma rest frame from Jüttner (1911), as in Wardzinski & Zdziarski (2000),

(B.1)

where K₂ is the modified Bessel function of the second kind of order 2 and Θ = k_BT/m_ec² is the normalized temperature. This formulation was adapted to our model because temperatures can exceed 10¹⁰ K behind shocks. The proton distribution can at first be assumed to be nonrelativistic at the temperatures reached by our simulations as Θ_p = k_BT/m_pc² ≪ 1 for T < 10¹² K, where the Maxwell-Jüttner and Boltzmann distributions coincide. An improved model should include a relativistic treatment of protons and e^± pair production because some shock zones may approach temperatures at whic either protons start to become relativistic or where pair creation starts to be effective.

The following expressions were derived in the plasma rest frame, but the calculated volumic power loss is invariant with frame because the volume dilatation is compensated for by the time dilatation of the boost. The numerical approximation of the Bessel K function used in A-MaZe is detailed in Appendix C.6.

B.1. Free-free radiation

Free-free radiation is emitted as an electron is accelerated by the Coulomb field of an ion. The expression of the volumic power loss due to the free-free mechanism of electrons in a hydrogen - helium plasma, including relativistic corrections, can be found in Rybicki & Lightman (2008),

(B.2)

This expression shows a dependence on ρ² and T^1/2 at classical electron temperatures, which becomes T^3/2 at high temperatures due to relativistic corrections. Rybicki & Lightman (2008) suggested as a good numerical approximation for the frequency-averaged Gaunt factor in all temperature ranges, which is the value we used for the Cygnus X-1 runs. This approximation is good up to T ∼ 10⁹ K, but a better approximation is needed because the temperature in the jet can become relativistic. A better evaluation of this term was given in van Hoof et al. (2015), who provided the values for and as functions of a parameter ∝T⁻¹. An illustration of the coefficients is given in Fig. B.1. These functions were implemented in the Cygnus X-3 runs.

Fig. B.1. Gaunt factor for hydrogen and helium across the equivalent tabulated range of temperature covered in van Hoof et al. (2015) and used in our Cygnus X-3 runs. Above 109 K, . The dotted line shows the mean value used for the Cygnus X/-1 runs.

B.2. Synchrotron

Ghisellini (2013) derived the synchrotron power emitted by a single electron with a Lorentz factor and pitch angle θ (the angle between its velocity and magnetic field lines) in the flow frame,

(B.3)

with σ_T the Thomson scattering cross-section and U_B = B²/8π the magnetic field energy density. This power can be averaged over the pitch angle θ because the electron distribution is assumed here to be isotropic,

(B.4)

The volumic power loss is derived by integrating over the electron distribution,

(B.5)

with K₃ the modified Bessel function of the second kind of order 3. In our simulations, the magnetic energy density is the sum of the contributions from the jet and the star, weighted by the tracer J.

These two contributions were modeled differently in the Cygnus X-1 and Cygnus X-3 runs: the magnetic field was chosen to be constant in Cygnus X-1 as a first attempt to model losses, at the risk of overestimating synchrotron losses in the long run. Seeing no impact of radiative losses until late simulation times in the Cygnus X-1 runs, we chose to launch runs based on Cygnus X-3 and updated our assumptions on the magnetic field. Because many authors such as Perucho (2019) suggested that the jet magnetic field structure is presumably toroidal, we chose a linear decrease with distance for the jet inner field to reflect this assumption. The stellar magnetic field was assumed to be a dipole, decreasing as r⁻³, with r the distance to the stellar center. This assumption does not take effects such as increased magnetic field downstream of shocks in consideration and may cause us to underestimate the synchrotron losses in the beam and inner cocoon. A better treatment of synchrotron losses would require magnetohydrodynamical simulations, which are beyond the scope of this study.

B.3. Inverse Compton scattering

Following Ghisellini (2013), the radiated power per electron in the flow frame is

(B.6)

with ψ the incident photon angle and U_ph the radiative energy density in the flow frame. In this frame, the electron distribution is assumed to be isotropic, therefore (1 − β_ecosψ)² can be averaged over the solid angle, from which we obtain . The power loss of a single electron is then

(B.7)

The volumic power loss in the fluid frame is obtained by integrating over the electron distribution, which yields

(B.8)

Considering the star as the sole source of seed photons for inverse Compton scattering, we derive the radiative energy density in the rest frame of the flow moving with a speed v_j = β_jc. When we define θ as the angle between the photon direction and the flow direction in the star rest frame, the radiative energy density is

(B.9)

where r is the distance to the stellar center in the stellar rest frame. Synchrotron and inverse Compton cooling follow the same law, and their ratio is equal to the ratio of the magnetic and stellar photon energy density, U_B and U_ph, respectively.

B.4. Line and recombination cooling

This term accounts for the collisional excitation of resonance lines and dielectronic recombination, where an ion captures an electron into a high-energy level and then decays to the ground state. We assumed solar photospheric abundances and the thermodynamical equilibrium of the plasma (Saha equilibrium), although in reality, the recombination may be delayed and may not correspond to the actual temperature of the plasma. This term follows the law

(B.10)

with i the different ion species. To facilitate the calculations, the various ions were taken into account in a single parameter Λ(T) from Cook et al. (1989), such that we can take . This parameterization was then extended in temperature range and implemented numerically in Walder & Folini (1996) and subsequent works. We chose an upper temperature of 10^7.7 K for this process, which corresponds to the recombination of fully ionized iron and the Fe-α line. This very efficient process is only effective in the coolest and most external parts of the cocoon.

B.5. Scalings of radiative processes

The scalings of the various radiative losses with density, temperature, and distance to the star are listed in Table B.1. For free-free losses, the dependence of the Gaunt factors on temperature modifies the high-temperature scaling from T^3/2 in the Cygnus X-1 runs to T² in the Cygnus X-3 runs. For the synchrotron and inverse Compton losses, the term K₃/K₂(Θ⁻¹) is constant at low temperatures and proportional to T at relativistic temperatures, which explains the evolution of the scaling with temperature from a linear to a square power-law for these two processes. Last, the line and recombination losses have a non-power-law dependence on T. At low temperatures, which are found only in the outer cocoon, the main cooling process is line recombination. At higher temperatures, free-free losses take over. At the highest temperatures, the cooling is dominated by synchrotron or inverse Compton processes, depending on our choice for the magnetic field intensity.

B.6. Cooling time

The cooling time in the observer’s frame of a fluid particle with rest frame temperature T and Lorentz factor γ is defined as t_cool = γT/Ṫ, where the dot marks the derivation with respect to the proper time of the fluid. For a perfect gas, T = p/ℛρ, with ℛ the gas constant divided by the molar mass of the fluid. Therefore

(B.11)

with all thermodynamic quantities measured in the comoving frame of the flow. Two extreme cases can be considered: the isobaric case (ṗ = 0), where , and the isochoric case (), where t_c,ρ = γp/ṗ. From the definitions given in Sect. 2.1,

(B.12)

with Γ₁ = Γ/(Γ − 1). Then, using and considering as an approximation in the weakly relativistic case, we can approximate these timescales as

(B.13)

(B.14)

Assuming γ² = 1 to approximate these cooling times, t_c, p ∝ 10²¹ρ/P_rad and t_c, ρ ∝ 1.5p/P_rad. In all the cases considered in the article, the isochoric cooling time is the shortest by about two orders of magnitude.

Appendix C: Detailed numerical methods

We detail here the numerical methods we used to perform the simulations described in this paper. We use the hydrodynamical framework of the A-MaZe simulation toolkit as described in Popov et al. (2019), on a Cartesian static mesh and without the well-balanced option.

C.1. Integration scheme

Semidiscretization of Eq. 7 in space results in

(C.1)

Here, dx, dy, and dz represent the spatial discretization in the x-, y-, and z-direction, and U_i, j, k is the vector of the discrete conserved variables at cell centers (i, j, k)∈(1, …, N_x, 1, …, N_y, 1, …, N_z), with N_x, N_y, and N_z the number of cells in the x−, y−, and z−direction of the computational space. Half indices denote cell faces. F_{i ± 1/2, j, k}, G_{i, j ± 1/2, k}, and H_{i, j, k ± 1/2} denote the fluxes through the cell faces in the x-, y-, and z-direction. Ψ_i, j, k represents the source terms that are also evaluated at the cell centers.

We performed a time integration of the N_x × N_y × N_z-D system of ordinary differential equations, Eq. C.1 with a first-order Runge-Kutta method (forward Euler method), although A-MaZe also offers strong stability-preserving (SSP) higher-order integration schemes (Shu & Osher 1988; Gottlieb et al. 2001).

We used a simple central scheme to evaluate the fluxes, as detailed here for the flux through the right x-interface of cell (i, j, k),

(C.2)

are the limited reconstructed variable values to the left and right of the cell interface i + 1/2, j, k. We used linear reconstruction and minmod limiters. λ_max is the highest characteristic speed. This integrator is relatively diffuse, but easy to implement for any hyperbolic system of equations and for an arbitrary EoS.

C.2. Inversion scheme

To solve these equations for the conservative variables (D, S^j, τ), primitive variables (ρ, v^j, p) are also necessary to compute the fluxes ℱⁱ. The following system is obtained from the definition of the conservative variables:

(C.3)

(C.4)

(C.5)

where ξ = γ²ρh needs to be determined to derive the primitive variables. Del Zanna & Bucciantini (2002) suggested a method adapted to the polytropic EoS. With Eqn. C.5 and Γ₁ ≡ Γ/(Γ − 1), we obtain

(C.6)

Using the definitions of the Lorentz factor and S,

(C.7)

The two expressions for ξ are combined to obtain the final equation for γ,

(C.8)

which is solved numerically using the Brent method (Brent 1973). Primitive variables were then computed using the formulas

This method is quite efficient, but is only valid for a constant Γ-law EoS. Mignone & McKinney (2007) suggested a discussion of the validity of a constant Γ-law EoS and an inversion method suitable for all EoS, but we find that the method described above is suitable in our case.

C.3. Benchmark for the scheme

Central schemes similar the one used in this paper have been widely used to perform (magneto-)hydrodynamical simulations (e.g., Del Zanna & Bucciantini 2002; van der Holst et al. 2008; Del Zanna et al. 2007). These schemes are easy to implement and very robust, but are relatively diffuse (see, e.g., Tóth & Odstrčil 1996 for a discussion). As these schemes are not based on (even partial) characteristic decomposition, contact interfaces in particular are smeared out relatively strongly. This has consequences for the growth of instabilities along these interfaces.

Central schemes have become more popular again because more sophisticated Riemann solvers – in particular exact solvers – are very CPU costly. Moreover, they are also not really adapted to the situation when more complex physics is involved in addition to (M)HD. Flows that include radiation, gravity, and particles show a different wave pattern, and waves have different velocities than pure (M)HD waves.

There are two ways to overcome the large diffusivity of central schemes: 1) higher-order spatial reconstruction schemes as proposed in Del Zanna et al. (2007), for example, or 2) meshes with a finer spatial discretization. This second approach was chosen for this work, where we used concentrated fine meshes along the beam of the jet where the instabilities develop. Ideally, the two approaches and the adaptive mesh algorithm implemented in A-MaZe may be combined for a more relevant mesh refinement.

C.3.1. Basic tests of the adiabatic scheme

The central scheme has been used by the authors for other work (Folini et al. 2004). We tested the implementation of SR by performing about 20 tests as proposed in the literature and found that we can reproduce these results well. Here, we only show two examples. The first is the relativistic blast wave problem as originally proposed by Donat et al. (1998). This Riemann problem is defined by setting the state to the left or right of the original interface located at 0.5 to (ρ, v, p)_L = (1, 0, 1000) and (ρ, v, p)_R = (1, 0, 0.01), resulting in a γ = 6 blast shock propagating to the right and a strong rarefaction fan propagating to the left. The solution at t = 0.35 on a very fine mesh of 25600 cells is shown in the left panel of Fig. C.1. The problem is tough and demonstrates why relativistic hydrodynamics is a numerical challenge.

Fig. C.1. Solution of the relativistic blast problem with a Lorentz factor γ = 6 shock detailed in the text at time t = 0.35. Left: Pressure (P), density (D), and velocity (V) as computed on a very fine mesh (25600 cells) showing the shock, the contact discontinuity, and the rarefaction fan. Right: Zoom into the thin high-density layer between the shock and the contact discontinuity. This feature is the most difficult to resolve. We show the solutions based on different resolutions: 400 cells (black), 1600 cells (blue), 6400 cells (pink), 12800 cells (green), and 25600 cells (orange).

Donat et al. (1998) presented a solution based on a third-order scheme combined with the Marquina solver, which usses the full spectral decomposition. Del Zanna & Bucciantini (2002) presented two solutions of the same problem based on a mesh of 400 cells. The first solution was computed with a third-order scheme based on the Harten, Lax, van Leer (HLL) solver (which uses only a part of the spectral information) and a monotized central (MC) limiter (their third-order convex essentially non-oscillatory - or CENO3 - scheme). The second solution is based on the same method as used in this paper, the second-order Lax-Friedrichs-scheme and minmod limiters. The hard part to compute is the thin high-density shell between the shock wave and the contact interface. These shells are typical for relativistic flows. In Fig. C.1, at t = 0.35, it is located between x = 0.84 and x = 0.85, where the density jumps by about two orders of magnitude in the shock wave and by three orders of magnitude in the contact interface. Most of the mass is concentrated within a region covering only about 1% of the domain.

Based on a discretization of 400 cells, none of the described schemes resolved the shell: The third-order schemes of Donat et al. (1998) and Del Zanna & Bucciantini (2002) reached a density of about 7.3 in 1-2 cells, and the second-order LF-scheme in Del Zanna & Bucciantini (2002) reached a density of about 6.5 over 1-2 cells. The correct value is about 10.5, however. The convergence of our scheme to the correct density value is demonstrated in the right panel of Fig. C.1. With 400 cells, our scheme is in line with that of Del Zanna & Bucciantini (2002). When 1600 cells are used, the density peaks at about 9.0 in 1-2 cells and at almost the correct value when 6400 cells are used. With 12800 cells, the density peak is well resolved, but the contact interface is still somewhat smeared out. The simulation using 25600 cells fully resolves the thin shell with some tens of cells, and the transition to the contact is quite sharp. We note that the computational costs for our scheme for 1600 cells is probably not (much) more than using a third-order scheme and a Riemann solver based on spectral decomposition of 400 cells. This demonstrates that a strategy based on fine meshes and a simple solver can be efficient. Admittedly, data files are more heavy than those produced on a 400 cell mesh, however.

The second test is the jet test-case proposed in Del Zanna & Bucciantini (2002): In cylindrical geometry, a γ = 7.1 jet is launched into a uniform environment with a low pressure, corresponding to a relativistic Mach number of about 17.9. This test is harder to simulate than the jets presented in this paper. Twenty cells covered the beam width, and the mesh in the domain was 160x400 in radial- and z-direction, respectively. Comparing the result obtained with our code (see Fig. C.2) with Fig. 9 of Del Zanna & Bucciantini (2002), we observe an excellent agreement in the position of the front bow-shock, the position of the Mach stem at the end of the beam, the position of the cross-shocks in the beam, and the general shape of the cocoon. In our case, the interface between inner and outer cocoon is more smeared out. The smaller modes in the instability developing along this interface are less resolved than in Del Zanna & Bucciantini (2002). This discrepancy is natural because Del Zanna & Bucciantini (2002) used the more accurate CENO3 scheme, while our result is based on the second order in the space LF method. However, this drawback can be overcome by using a finer mesh (not shown).

Fig. C.2. γ = 7.1 (Mach ≈17.9) jet simulated in Del Zanna & Bucciantini (2002) and reproduced by the scheme used for this work. The image shows time = 40 and can be directly compared with the last panel in Fig. 9 of Del Zanna & Bucciantini (2002).

C.3.2. Uncertainty for simulations of turbulent and cooling flows

The exactness of the simulation presented in this paper is harder to estimate. The flows are turbulent, and cooling introduces more instabilities, waves, and interfaces. The turbulent region of the cocoon has no fixed boundary, but is connected by shocks and material interfaces to the environment. Interior turbulent fluctuations will impact the shape of the interfaces and, inversely, the dynamics of the interfaces will act on the interior turbulence. Based on these arguments, we cannot expect to find a converged solution in the sense demonstrated in Fig. C.1. We have to trust the general correctness of the scheme and have to give some reasons why the presented solutions are close to correct. A rigorous error analysis based on statistical analysis of many simulations that differ slightly in their initial conditions would be desirable, but is sophisticated, complex, and computationally expensive and thus beyond reach for this study. A step in this direction, nevertheless, is presented in Fig. C.3. We list in the following some points that shed some light on the uncertainties.

Fig. C.3. Simulation of the fiducial case of CygX1 on a mesh that is twice as fine as the simulation presented in Fig. 3. We show the simulation at 6000 s, corresponding to the middle panel of Fig. 3.

1. Turbulence: Reynolds numbers are too high to resolve the turbulent cascade with any numerical scheme. Moreover, ideal hydrodynamics does not treat diffusion explicitly. A numerical scheme implicitly introduces a certain diffusion (Hirsch 2006; LeVeque et al. 2002), however, which is much higher in astrophysical rarefied flows than the physical diffusion. However, as pointed out by Boris et al. (1992) and further explored by Porter et al. (1992) and Porter & Woodward (1994), finite-volume methods such as ours cut the turbulent cascade in a way that does not lead to an energy pile-up or -sink at the numerical diffusion scale, thus cutting the cascade correctly, at least to first order. This approach is called monotone integrated large-eddy simulation (MILES). A more rigorous study of the MILES approach is given in Garnier et al. (1999), for example, and a summary of the idea and more references can be found in Folini & Walder (2006). These studies show that turbulent flows are relatively well captured by finite-volume methods and do not introduce large errors.

2. Cooling: Radiative shocks are prone to an overstability whenever the slope of the cooling law is sufficiently shallow or negative. For a radiative cooling parameterized as a function of density and temperature with Λ(T) = Λ₀T^β, which applies for free-free and line cooling, Chevalier & Imamura (1982) and Bertschinger (1986) have shown that the overstability is present whenever β ≲ 0.4 (fundamental mode), and β ≲ 0.8 (first-overtone mode). We have shown (Walder & Folini 1996) that the presence and amplitude of the overstable modes in a numerical study critically depend on the numerical resolution because smeared-out interfaces radiate more than better resolved interfaces. The resolution we chose for the simulations is sufficient to resolve the overstability (not shown). We add two remarks, however. First, the numerical model we used does not include mass diffusion and, in particular, heat diffusion, which physically determines the smearing of the interface. It is thus not clear whether we under- or over-estimate this particular effect. Second, radiative multidimensional shocks can generate and drive turbulence (Walder & Folini 1998) and turbulent thin shells (Folini & Walder 2006; Folini et al. 2014).

3. Resolution comparison: The fiducial case of CygX1 was simulated on a mesh twice finer than the generic mesh up to about 10’000 s. The snapshot at 6000 s is shown in Fig. C.3. This can be compared to the snapshot of the generic case shown in the middle panel of Fig. 3. The comparison shows that the instability sets in at about the same time in both simulations. However, in the simulation on the finer mesh, the jets propagate about 10% faster than in the simulation on the generic mesh. We also observe similar effects, of the same order, in 1D test simulations. This is expected on the basis of the arguments given in the point above. Better resolving the contact interface at the head of the jet will reduce cooling there, leaving slightly more energy to push the bow shock to a larger distance. An error of 10% is quite a good result for most large-scale fluid-dynamical simulations. We stress again that the correct jet speed depends on the physical diffusion.

To even improve the confidence in the solutions presented in this paper, we ran a 2D resolution study, using the generic and the finer mesh in the 3D case. This study also covered the turbulent phase. Again we find that the essential features of the jet (number and location of the cross shocks, cocoon shape, and time at which the instability and the turbulence set in) are independent of the resolution. However, the 2D simulations cannot directly be compared to the 3D simulations because the character of the turbulence is different in 2D and 3D.

C.4. Necessity of relativistic simulations

Because computational costs are considerably higher for a relativistic simulation, it might be questioned whether it is necessary to perform relativistic simulations to obtain correct solutions for the mildly relativistic problems presented in this paper, with γ_b ≈ 1.06 for CygX1 and γ_b ≈ 1.51 for CygX3. However, even these small Lorentz factors lead to a significant difference in the jet propagation between a relativistic and a Newtonian simulation. This is illustrated in Fig. C.4, which shows 1D simulations at (observer) time t = 6500 s of the jet propagation of CygX3, including all cooling terms. The jet head is located in the thin high-density shell. In the Newtonian case, this shell in the observers frame is located at about x = 45 ⋅ 10¹² cm. The shell in the relativistic case is located at x = 71 ⋅ 10¹² cm. Both simulations used 12800 cells. Thus, the relativistic jet head propagates about 1/3 faster than the Newtonian jet head. This can be explained on the basis of Eq. 5: the ratio of η^* and η for CygX3 (see Table D.4) is about 2.3, resulting in a difference in the jet-propagation speed of about 45%. The difference for CygX1 is smaller, but still about 5%. In 1D, the shocked beam will cool down, in contrast to the multidimensional simulations, in which the beam is regularly reheated by the cross shocks.

Fig. C.4. Comparison in density profiles of 1D simulations for parameters similar to the fiducial cases. We show in the observer frame the comparison between a Newtonian and a relativistic simulation for CygX1 at 30 000 s (blue) and for CygX3 for 6 500 s (red). The simulations include all cooling terms. The pattern with forward and reverse shock and a thin layer of cooled gas to the left of the contact interface between beam and environment material is similar for all simulations. The mesh consists of 12800 cells, about as many as the mesh covering the beam in the 3D simulations presented in the paper.

C.5. Impact of the adiabatic index choice on jet propagation

Because the simulated flows are relativistic, the choice of using a constant adiabatic index of value 5/3 may be discussed. We show in Fig. C.5 a 1D test similar to the one presented in Sect. C.4, comparing jets with Γ = 5/3 and Γ = 4/3 with the relativistic and the Newtonian solver. Changing the adiabatic index has an almost negligible effect on the jet head (the thin high-density shell) propagation, but jets with Γ = 5/3 present a more advanced front shock by 25% in the Newtonian case, and it is more advanced by 16% with the relativistic solver. This is caused by the higher post-shock densities for Γ = 4/3, which result in a stronger cooling of the shocked gas. The situation is different in 2 and 3D, however, because Mignone & McKinney (2007) reported that jets with a smaller adiabatic index propagate faster. This is due, as mentioned Sect. 3.1 (following Martí et al. 1997), to the first recollimation shock in the beam, which is strong enough to reaccelerate the beam flow at Γ = 4/3.

Fig. C.5. Comparison in density profiles of 1D simulations similar to the fiducial case CygX3. Everything is the same as in Fig. C.4, with the exception of the adiabatic index Γ. Little effects are observable on the jet head (the thin high-density shell) propagation when Γ is varied, but choosing Γ = 5/3 results in a more advanced front shock by 25% in the Newtonian case and 16% in the relativistic one. This is caused by a stronger cooling of the shocked gas due to higher densities.

C.6. Numerical approximation of Bessel K functions

Equations B.1 and therefore Eqns. B.5 and B.8 use the modified Bessel function of the second kind (also called Bessel K function or Macdonald function), especially the ratio K₃/K₂. A Fortran 90 implementation of this function by Moreau (2005) was ported to A-MaZe, but as both functions tend to zero at low temperature, a simple division of K₃(Θ⁻¹) by K₂(Θ⁻¹) caused underflows during calculations. We therefore modified the method to derive the ratio directly. Figure C.6 compares our Fortran method with the built-in Bessel K functions from the SciPy package and shows the stability of our method over the whole temperature range compared to a simple division.

Fig. C.6. Comparison of our Fortran method for the ratio K3/K2 with the simple division of both terms using the functions defined in Python package SciPy. Our Fortran method avoids underflows for high values of Θ−1 with a relative error that never exceeds 10−8.

Appendix D: Simulations setups

Tables D.1 and D.3 show the input parameters for each run in CGS units. Tables D.2 and D.4 show the same inputs in terms of the a-dimensional parameters introduced Sect. 3.1.

Appendix E: Instability growth phase

Figures E.1 and E.2 show pressure slices of the runs CygX1 and CygX3 runs, respectively. They highlight the internal structure of the beam with alternating over- and under-pressured regions compared to the fluid in the surrounding cocoon.

Fig. E.1. Pressure slices during the instability growth phase of the fiducial Cygnus X-1 run CygX1. The color scale is fixed from 1 (blue) to 1000 Ba (red) to better highlight the beam structure. The beam shows alternating over- and underpressured zones whose number has risen at the 5000 s mark. The inner cocoon shows ripple-like structures alternating on either side of the beam with increasing intensity as the jet evolves.

Fig. E.2. Pressure slices during the instability growth phase of the fiducial Cygnus X-3 run CygX3. The color scale is fixed from 103 (blue) to 105 Ba (red) to better highlight the beam structure.

All Tables

All Figures

Fig. 1. Evolution of the loss processes with rest frame temperature in Cygnus X-1 (left) and Cygnus X-3 (right). The thin black lines shows the various temperature scalings detailed Table B.1. Line recombination losses (“line”) are not drawn for T > 107.7 as they are disabled above this temperature. Colored shading shows synchrotron losses (“syn”) when the stellar magnetic field B⋆ is either multiplied or divided by 5, and the shading around inverse Compton losses (“ic”) illustrates a ±10% uncertainty on T⋆. The values for the physical quantities have been chosen in the outer cocoon leeward side in fiducial runs. They allow for a clear showcasing of the loss scaling, but are not representative of the losses over the simulation. We refer to Fig. 7 to compare them over the whole jet. For Cygnus X-1: ρ = 10−15 g cm−3, B⋆ = 10 G, T⋆ = 3 × 104 K, (x, y, z) = (1.5 × 1012, 4.05 × 1013, 4 × 1013) cm, (vx, vy, vz) = (108, 108, 108) cm s−1. For Cygnus X-3: ρ = 10−14 g cm−3, B⋆ = 100 G, T⋆ = 8 × 104 K, (x, y, z) = (1.5 × 1012, 3.05 × 1013, 3 × 1013) cm, (vx, vy, vz) = (108, 5 × 108, −103) cm s−1.

In the text

Fig. 2. Structure of the computational grid over the whole domain, illustrated for a rest mass density slice along the plane containing the star and jet center. The density scale extends from 10−19 (deep blue) to 10−12 g cm−3 (red), same as Fig. 3. The grid contains five refinement levels: the whole domain is refined twice (green and blue interfaces) by a factor of 4 and then twice more (cyan and magenta interfaces) by a factor of 2 to attain a factor of 64 in the finest levels, which are better shown in the zoom into the jet injection at the orbital scale in the bottom right corner of the image.

In the text

	Fig. 3. Rest-mass density slices of run CygX1 at times (top to bottom) t = 2000, 6000, and 12 000 s, showing the three evolutionary phases detailed in the text in Sect. 3.2.
In the text

	Fig. 4. Dynamical and structural evolution of fiducial run CygX1. Left: position and speed of the jet head for fiducial run CygX1. Right: jet volume (beam and inner and outer cocoon) and aspect ratio. The limits of each evolutionary phase are marked by the vertical dash-dotted lines in the various panels.
In the text

	Fig. 5. Rest-mass density slices of run CygX3 at times (top to bottom) t = 400, 1200, and 2500 s, showing the three evolutionary phases detailed in the text in Sect. 3.2.
In the text

	Fig. 6. Dynamical and structural evolution of fiducial run CygX3. Left: position and speed of the jet head for fiducial run CygX3. Right: jet volume (beam and inner and outer cocoon) and aspect ratio. The limits of each evolutionary phase are marked by the vertical dash-dotted lines in the various panels.
In the text

Fig. 7. Time evolution of radiative losses for the fiducial runs CygX1 (top) and CygX3 (bottom), derived at each cell and summed per zone, at time steps t = (500, 2000, 4000, 6000, 10 000, 15 000) s for CygX1 and t = (250, 500, 750, 1500, 3000, 6000) s for CygX3 (two data points per evolutionary phase). In both cases, free-free losses dominate cooling in the inner cocoon and beam, while line cooling dominates the outer cocoon.

In the text

Fig. 8. Effects of losses on run CygX1 structure and dynamics. Top: temperature slices of runs CygX1 and CygX1_noLoss at time t = 15 000 s. Both jets display similar structures, with the exception of a slightly larger outer cocoon at the head of the non-cooled jet. Bottom: jet head propagation and speed of the same runs, CygX1 in red and CygX1_noLoss in blue. The theoretical 1D propagation from Appendix A.1 is drawn as dotted lines following the same color-coding. The propagation is identical in both runs until the start of the turbulent phase, after which speed fluctuations differ, but the average propagation speed is identical in the two runs with almost no difference in the jet head position plot.

In the text

Fig. 9. Effects of losses on run CygX3 structure and dynamics. Top: temperature slices of runs CygX3 and CygX3_noLoss at time t = 9000 s. Two main differences appear: (1) The cooled beam is thinner. Its envelope closely follows the internal shocks structure in contrast to the non-cooled case. (2) The cocoon of the non-cooled jet expands farther at its basis, almost wrapping around the star, whereas in the cooled case, the cocoon has almost disappeared because ambient material has cooled enough to be blown back by the wind. Bottom: same as Fig. 8. The cooled jet (in red) is initially faster, but leaves the smooth phase earlier, after which point it is slower on average, as seen in the propagation plot.

In the text

	Fig. 10. Effect of radiative losses on the CygX3 beam structure. Top: ratio of volume-averaged pressure of the inner cocoon and beam for runs CygX3 (red) and CygX3_noLoss (blue). The overpressure is always higher in the cooled case. Bottom: tracer density at t = 750 s for runs CygX3 and CygX3_noLoss. The cooled jet features stronger oblique shocks.
In the text

	Fig. 11. Volumes of the different jet zones as a function of time are affected differently by radiative cooling. The outer cocoon is larger without losses, while the volumes are comparable for the inner cocoon. The power-law dependence on time of the cocoon is robust for the cocoon, but not for the beam.
In the text

Fig. 12. Jet head propagation and speed for runs CygX1_T7 (red), CygX1 (blue) and CygX1_T9. Run CygX1 slows down to the turbulent phase earlier than CygX1_T7, but display the same average speed during the turbulent phase, as shown by the almost parallel propagation curves. Run CygX1_T9 displays a different behaviour, decelerating to a plateau in the instability growth phase and with a lower average speed than the other two runs, which appears to decelerate after t = 12 000 s.

In the text

Fig. 13. Comparison of temperature and overpressure of runs CygX1_T7, CygX1 and CygX1_T9, the color-coding is the same as in Fig. 12. Top: evolution of the zone-averaged temperature with simulation time. CygX1 and CygX1_T7 differ only slightly at first, but then start to evolve differently after the 5000 s mark. CygX1_T9 shows an overall higher temperature, with a small peak in inner cocoon temperature around t = 6000 s. Bottom: evolution of the ratio of the inner cocoon to beam pressure, as defined Sect. 3.3.1. CygX1 and CygX1_T7 present similar values up to t ∼ 5000 s, while CygX1_T9 presents a stronger overpressure.

In the text

Fig. 14. Comparison between runs CygX1 (red), CygX1_T7 (blue) and CygX1_T9 (magenta). Probability density function of the temperature at t = 4000 s (full lines) and t = 10 000 s (dash-dotted lines). In the early stages, CygX1 and CygX1_T7 are indistinguishable from each other, while CygX1_T9 is slightly hotter and more of its volume is hotter than ∼5 × 109 K. In the mixing phase, the temperature repartition is more in line with the injected temperature.

In the text

Fig. 15. Sensitivity of the jet head position and speed to kinetic power and wind ram pressure for Cygnus X-1. The jet kinetic power is divided by 10 from CygX1 (red) to CygX1_mP (blue). In the weak case, the jet is slower, but the position fits the theoretical evaluation for a longer time. The speed diagram shows no oscillations. The wind speed is increased by 50% for run CygX1_wind (green), which results in a higher starting speed, a shorter reacceleration phase, and a weaker bump in the second deceleration phase. The average speed in the turbulent phase is weaker. The speed and position plots differ faster from the theoretical 1D values in the strong wind case.

In the text

	Fig. 16. vx/vj slice for run CygX1_mP at time t = 10 000 s. The jet is heavily bent from the wind effects, and its velocity breaks down before it arrives at the head.
In the text

Fig. 17. Sensitivity of the jet head position and speed to kinetic power (left) and wind momentum (right) for the Cygnus X-3 runs. Top left: a division of the jet power by 2 decreases all dynamical values, starting from the mean propagation speed. In particular, the acceleration is far less efficient: it shows a gain of ∼40% in speed between the starting point and the peak against a gain of ∼120% in the CygX3 run. Top right: a reduction of Ṁ⋆ by 25% and vw by 33% (keeping η constant) lengthens the smooth phase and strengthens the speed fluctuations in the following phases. The jets settle to the same median speed in the turbulent phase. Bottom right: same reduction as in the top right panel, with half the jet power. In the CygX3_mPmW case (right, blue), the initial acceleration phase is followed by a speed plateau.

In the text

	Fig. B.1. Gaunt factor for hydrogen and helium across the equivalent tabulated range of temperature covered in van Hoof et al. (2015) and used in our Cygnus X-3 runs. Above 109 K, . The dotted line shows the mean value used for the Cygnus X/-1 runs.
In the text

Fig. C.1. Solution of the relativistic blast problem with a Lorentz factor γ = 6 shock detailed in the text at time t = 0.35. Left: Pressure (P), density (D), and velocity (V) as computed on a very fine mesh (25600 cells) showing the shock, the contact discontinuity, and the rarefaction fan. Right: Zoom into the thin high-density layer between the shock and the contact discontinuity. This feature is the most difficult to resolve. We show the solutions based on different resolutions: 400 cells (black), 1600 cells (blue), 6400 cells (pink), 12800 cells (green), and 25600 cells (orange).

In the text

	Fig. C.2. γ = 7.1 (Mach ≈17.9) jet simulated in Del Zanna & Bucciantini (2002) and reproduced by the scheme used for this work. The image shows time = 40 and can be directly compared with the last panel in Fig. 9 of Del Zanna & Bucciantini (2002).
In the text

	Fig. C.3. Simulation of the fiducial case of CygX1 on a mesh that is twice as fine as the simulation presented in Fig. 3. We show the simulation at 6000 s, corresponding to the middle panel of Fig. 3.
In the text

Fig. C.4. Comparison in density profiles of 1D simulations for parameters similar to the fiducial cases. We show in the observer frame the comparison between a Newtonian and a relativistic simulation for CygX1 at 30 000 s (blue) and for CygX3 for 6 500 s (red). The simulations include all cooling terms. The pattern with forward and reverse shock and a thin layer of cooled gas to the left of the contact interface between beam and environment material is similar for all simulations. The mesh consists of 12800 cells, about as many as the mesh covering the beam in the 3D simulations presented in the paper.

In the text

Fig. C.5. Comparison in density profiles of 1D simulations similar to the fiducial case CygX3. Everything is the same as in Fig. C.4, with the exception of the adiabatic index Γ. Little effects are observable on the jet head (the thin high-density shell) propagation when Γ is varied, but choosing Γ = 5/3 results in a more advanced front shock by 25% in the Newtonian case and 16% in the relativistic one. This is caused by a stronger cooling of the shocked gas due to higher densities.

In the text

	Fig. C.6. Comparison of our Fortran method for the ratio K3/K2 with the simple division of both terms using the functions defined in Python package SciPy. Our Fortran method avoids underflows for high values of Θ−1 with a relative error that never exceeds 10−8.
In the text

Fig. E.1. Pressure slices during the instability growth phase of the fiducial Cygnus X-1 run CygX1. The color scale is fixed from 1 (blue) to 1000 Ba (red) to better highlight the beam structure. The beam shows alternating over- and underpressured zones whose number has risen at the 5000 s mark. The inner cocoon shows ripple-like structures alternating on either side of the beam with increasing intensity as the jet evolves.

In the text

	Fig. E.2. Pressure slices during the instability growth phase of the fiducial Cygnus X-3 run CygX3. The color scale is fixed from 103 (blue) to 105 Ba (red) to better highlight the beam structure.
In the text