1. Introduction
Jets and outflows represent a common and crucial process in astrophysical systems; from forming stars (Herbig-Haro objects) to supermassive black holes in active galactic nuclei (AGN). In each of these systems, the properties of the jets (velocities, densities, composition…) are different, but they share a common association with accretion and magnetohydrodynamical processes that explain their formation, collimation and acceleration (see, e.g. Komissarov Reference Komissarov2012).
We know that extragalactic, AGN jets, can be modelled as fluids or plasmas because the typical Larmor radius of the jet particles is orders of magnitude smaller than the spatial scales of the system (Blandford & Rees Reference Blandford and Rees1974). The same is true for the case of jets generated in collapsars and X-ray binary stars. Despite technical advances, it is still impossible to reproduce the most extreme of these systems, i.e. to accelerate plasmas to relativistic velocities, in laboratories to study their properties (e.g. Bellan Reference Bellan2018).Footnote 1 This leaves numerical simulations as the only option to approach the nature of these systems, which can be central in the evolution of the interstellar medium in host galaxies and clusters (e.g. McNamara & Nulsen Reference McNamara and Nulsen2007; Fabian Reference Fabian2012), or in the generation of high and very high energy radiation in the Universe (e.g. Rieger & Levinson Reference Rieger and Levinson2018).
Relativistic jets are formed in the environment of compact objects by the extraction of rotational energy from the accreting body and/or the accretion disk (Blandford & Znajek Reference Blandford and Znajek1977; Blandford & Payne Reference Blandford and Payne1982). The idea is, in both cases, related to the twisting of magnetic field lines, which, by resisting this deformation, are able to extract rotational energy from the system. The magnetic field is supposed to be dragged onto the central object by the accreting gas itself. In the case of Kerr black holes, it anchors to the hole's ergosphere and is thus forced to rotate faster inside than outside of this region, which creates the aforementioned twisting (Blandford & Znajek Reference Blandford and Znajek1977; Komissarov Reference Komissarov2012).
In all cases, jet evolution involves many orders of magnitude in distance. In the case of AGN jets, jet formation takes place at spatial scales of the order of the gravitational radius of the central supermassive black hole, i.e. ${\sim }10^{10}\,{\rm m}$, or, if the inner regions of the accretion disk are involved, it may extend to a factor 10–100 larger. They reach, in contrast, distances up to ${\sim }1\,{\rm Mpc}$, i.e. ${\sim }10^{22}\,{\rm m}$, as revealed by low frequency surveys (see, e.g. Icke Reference Icke1991; Molina et al. Reference Molina, Bassani, Malizia, Bird, Bazzano, Ubertini and Venturi2014; Oei et al. Reference Oei, van Weeren, Hardcastle, Botteon, Shimwell, Dabhade, Gast, Röttgering, Brüggen and Tasse2022; Simonte et al. Reference Simonte, Andernach, Brüggen, Best and Osinga2023). As a consequence, the computational resources required to study the whole evolution are prohibitive, even for modern architectures, and studies focus on different regions. Depending on this, the properties of the jets vary and the numerical codes have to be adapted to them.
When ejected, jets consist mainly of a Poynting flux, which may be loaded by a dilute, hot gas with an electron–positron pair composition the pairs would be generated by photon–photon collisions at the black hole corona and/or an electron–proton gas if loaded from the accretion disk. These jets are launched with a relatively small radial velocity. Taking into account the large values achieved by magnetic fields in these regions, their velocities must initially be sub-Alfvènic. Nevertheless, relativistic speeds are reported at sub-parsec to parsec scales, which means that jets need to be accelerated to super-Alfvènic (super-fast-magnetosonic) speeds. The mechanisms that are typically invoked to explain this acceleration are both related to jet expansion with magnetic and internal energy contributions (e.g. Vlahakis & Königl (Reference Vlahakis and Königl2004), Komissarov et al. (Reference Komissarov, Barkov, Vlahakis and Königl2007) and Ricci et al., in preparation).
At kiloparsec scales, jets show different morphologies and clear signs of interaction with the warm–hot intergalactic medium (WHIM, e.g. Croston et al. Reference Croston, Hardcastle, Mingo, Best, Sabater, Shimwell, Williams, Duncan, Röttgering and Brienza2019). A coarse, initial dichotomy was established by Fanaroff & Riley (Reference Fanaroff and Riley1974) between Fanaroff–Riley type I (FRI) sources and Fanaroff–Riley type II (FRII) sources. The former are brighter within their host galaxies, jets become dimmer with distance and show brightness symmetry at large distances (e.g. Laing & Bridle Reference Laing and Bridle2014). The FRIIs are brighter at the impact sites with the WHIM, the so-called hotspots, and show brightness asymmetry along their length between the jets and their counter-jets. The common interpretation is related to jet velocity: on the one hand, FRII jets seem to keep their collimation and mildly relativistic velocities up to the interaction site with the WHIM, which means that the flow is supersonic at impact, triggering strong shocks that show up as hotspots, where Doppler boosting explains the brightness asymmetry; on the other hand, in the case of FRI jets, the loss of collimation and brightness symmetry at large scales point towards deceleration and weakening of the Doppler asymmetry (Bridle Reference Bridle1984; Bridle & Perley Reference Bridle and Perley1984).
In the case of microquasars, although jets evolve on time scales much shorter than AGN jets (see, e.g. Mirabel (Reference Mirabel2010) and Fender (Reference Fender, Lewin and van der Klis2006), for a review), a large number of unknowns remain about fundamental aspects such as their composition (Migliari, Fender & Méndez Reference Migliari, Fender and Méndez2002), internal dynamics and its relation to the formation mechanism and changes in the accretion flow (Fender, Belloni & Gallo Reference Fender, Belloni and Gallo2004), or the mechanisms that generate the production of high- to very-high-energy radiative output within the binary scales (see, e.g. Perucho & Bosch-Ramon Reference Perucho and Bosch-Ramon2008; Bosch-Ramon & Khangulyan Reference Bosch-Ramon and Khangulyan2009a; Perucho, Bosch-Ramon & Khangulyan Reference Perucho, Bosch-Ramon and Khangulyan2010a; Bosch-Ramon et al. Reference Bosch-Ramon, Barkov, Khangulyan and Perucho2012a; Perucho & Bosch-Ramon Reference Perucho and Bosch-Ramon2012; Bosch-Ramon, Barkov & Perucho Reference Bosch-Ramon, Barkov and Perucho2015; Paredes-Fortuny et al. Reference Paredes-Fortuny, Bosch-Ramon, Perucho and Ribó2015; de la Cita et al. Reference de la Cita, Bosch-Ramon, Paredes-Fortuny, Khangulyan and Perucho2017; López-Miralles et al. Reference López-Miralles, Perucho, Martí, Migliari and Bosch-Ramon2022), even far from the sub-milliarscsecond resolution achieved by very large baseline interferometry at cm–mm wavelengths in the radio band. Furthermore, the transition of the jet from these compact regions (${\sim }10^{10}\,{\rm m}$) to parsec scales (${\sim }10^{16}\,{\rm m}$) has been barely studied by numerical simulations, and only few works have been devoted to microquasar jet evolution at the largest scales (e.g. Bordas et al. Reference Bordas, Bosch-Ramon, Paredes and Perucho2009; Monceau-Baroux et al. Reference Monceau-Baroux, Porth, Meliani and Keppens2015; Charlet et al. Reference Charlet, Walder, Marcowith, Folini, Favre and Dieckmann2022).
In all scenarios in which relativistic outflows are involved, and once accelerated, the jets are subject to dissipative processes such as shocks or mass entrainment via interaction with stars and the development of instabilities (Perucho Reference Perucho2019). The remarkable stability of jets on the different scales through which many of them are observed to propagate has been a matter of study during the last decades. The mechanisms by which extragalactic FRI jets are decelerated have also been deeply investigated (e.g. Perucho et al. Reference Perucho, Martí, Laing and Hardee2014; Massaglia et al. Reference Massaglia, Bodo, Rossi, Capetti and Mignone2016, Reference Massaglia, Bodo, Rossi, Capetti and Mignone2019, Reference Massaglia, Bodo, Rossi, Capetti and Mignone2022; Gourgouliatos & Komissarov Reference Gourgouliatos and Komissarov2018; Perucho Reference Perucho2020; Anglés-Castillo et al. Reference Anglés-Castillo, Perucho, Martí and Laing2021). Attention has been paid to the scenarios of generation of high-energy particles (e.g. neutrinos, see Murase & Stecker Reference Murase and Stecker2022; Buson et al. Reference Buson, Tramacere, Oswald, Barbano, Fichet de Clairfontaine, Pfeiffer, Azzollini, Baghmanyan and Ajello2023) and cosmic rays in these systems, both in galactic and extragalactic sources (e.g. Rieger & Levinson Reference Rieger and Levinson2018; Matthews, Bell & Blundell Reference Matthews, Bell and Blundell2020; Seo, Ryu & Kang Reference Seo, Ryu and Kang2022 and references therein). The particles could be accelerated via shock or shear acceleration, or also magnetic reconnection processes, i.e. as a byproduct of jet development. Recent works based on particle-in-cell simulations suggest a relation between synchrotron light polarization and jet composition (e.g. Meli & Nishikawa Reference Meli and Nishikawa2021). It is relevant to stress that a baryonic component is necessary to explain the generation of neutrinos, for instance. This component could be either provided by a disk-generated wind surrounding the inner leptonic jet spine, or by entrainment. Finally, the role of jets in heating galactic environments, avoiding the fast cooling and collapse of the WHIM onto the jet's host galaxy, also known as a feedback process, could be very relevant in galaxy evolution and is not only part of research on jet physics, but also on cosmological simulations (e.g. Tremmel et al. Reference Tremmel, Karcher, Governato, Volonteri, Quinn, Pontzen, Anderson and Bellovary2017; Weinberger et al. Reference Weinberger, Ehlert, Pfrommer, Pakmor and Springel2017; Vazza et al. Reference Vazza, Wittor, Di Federico, Brüggen, Brienza, Brunetti, Brighenti and Pasini2023).
In this paper, we summarize our recent work on jet physics, which is mainly focused on jet propagation, stability and interaction with the ambient medium at different scales in both extragalactic sources and X-ray binaries, via numerical simulations, and also present new results either based on revisiting published work or on new simulations. In § 2, we present the basic equations that need to be solved to simulate the evolution and physics of these systems, and the numerical techniques that we use to solve them. Section 3 is devoted to the different applications that we have performed in the last years, both in the case of relativistic hydrodynamics (RHD) and magnetohydrodynamics (RMHD), and the main results that we have obtained. Among these, we include new results such as the pseudo-synchrotron images derived from AGN jet simulations (figures 1 and 2 in § 3.1.1), and the new simulations on jet–interstellar medium (ISM) interactions (§ 3.1.3) as well as a new analysis of recent RMHD simulations (§ 3.2). Finally, our conclusions are given in § 4.
2. Numerical codes and methods
2.1. The RMHD equations
The stress–energy tensor for an ideal plasma is $T^{\mu \nu } = T^{\mu \nu }_G + T^{\mu \nu }_{EM}$. On the one hand, the gas contribution is given by
where the Greek letters take values $\mu,\ \nu =0,1,2,3$ and stand for the dimensions of the Minkowski space–time, $\rho$ is the flow rest-mass density, $h$ the gas specific enthalpy, $u$ is the four-velocity, $p$ is the gas pressure and $g$ the Minkowski metric, and we have used $c=1$. On the other hand, the field contribution is given by
where $F^{\alpha \beta }=\epsilon ^{\alpha \beta \mu \nu } b_\mu u_\nu$ are the elements of the electromagnetic field tensor satisfying the Maxwell equations
In this expression $\epsilon ^{\alpha \beta \mu \nu }$ is the Levi-Civita tensor, $b_\mu$ are the components of the magnetic field in the fluid rest frame and $J^\nu$ is the charge four-current. Collecting the expressions above, we obtain the stress–energy tensor of a magnetized, perfect fluid
where $h^*= h + |b|^2 /\rho$ is the total enthalpy, $p^*$ the total (gas plus magnetic) pressure ($p^*=p_{\rm gas} + |b|^2/2$), $b^\mu$ is the relativistic magnetic field four-vector, with $b^0=\gamma (\boldsymbol {v}\boldsymbol {\cdot }\boldsymbol {B})$, and $b^i={B^i}/{W}+v^ib^0$, so $|b|^2 = b_\alpha b^\alpha = B^2/\gamma ^2 + (\boldsymbol {v}\boldsymbol {\cdot }\boldsymbol {B})^2$, with $\boldsymbol {B}$, $\boldsymbol {v}$ the magnetic field and velocity vectors in the laboratory frame, and $\gamma$ is the corresponding Lorentz factor. We note that a $\sqrt {4{\rm \pi} }$ factor is included in the definition of the magnetic fields.
The conservation laws that govern the dynamics of the plasma are
These equations can be written in matrix form as
where $\boldsymbol {U}=\{D,S^j,\tau _e,B^j\}$ is a vector of conserved variables, $D$ is the relativistic rest-mass density, $S^j$ is the momentum density of the magnetized fluid, $\tau _e$ is the energy density and $\boldsymbol {F}^i$ are the fluxes in each spatial direction. The relation of those conserved quantities and fluxes to their primitive physical variables is given by
Finally, from (2.3a,b), we obtain the divergence-free condition for the magnetic field, $\boldsymbol {\nabla } \boldsymbol {\cdot } \boldsymbol {B} = 0$.
This system of equations needs to be closed by means of an equation of state. It is common to use the ideal gas equation of state for this purpose. However, in the case of relativistic outflows propagating through their host galaxies, the different compositions of the jet flow (electron/positron pairs and/or electron/proton) and the ambient medium (ionized hydrogen in the intergalactic medium, but also atomic hydrogen in two phase media such as the ISM) recommend the use of the relativistic gas equation of state (Synge Reference Synge1957).
Source terms can be added to the equations, depending on the problems to be solved. Typical source terms used are gravity (or pseudo-gravity to keep a density/pressure profile in equilibrium, e.g. Perucho & Martí Reference Perucho and Martí2007; Perucho, Martí & Quilis Reference Perucho, Martí and Quilis2019), thermal cooling terms (e.g. Perucho, Bosch-Ramon & Barkov Reference Perucho, Bosch-Ramon and Barkov2017a), hydrogen ionization and recombination effects (e.g. Vaidya et al. (Reference Vaidya, Mignone, Bodo and Massaglia2015), Perucho et al. (Reference Perucho, López-Miralles, Reynaldi and Labiano2021), Perucho et al., in preparation) or radiation feedback (López-Miralles, Martí & Perucho Reference López-Miralles, Martí and Perucho2023).
It is worth noting that, on the one hand, relativistic effects arise when $v\rightarrow c$ and/or when $\varepsilon \rightarrow c^2$, and, on the other hand, magnetic fields can be dynamically relevant for even small values of the magnetic parameter $\sigma =|b^2|/(\rho h)$, which gives the ratio of magnetic to gas kinetic energy.
2.2. The numerical codes
Numerical codes are aimed at solving the RMHD equations (see Martí & Müller (Reference Martí and Müller2015) and Martí (Reference Martí2019) for reviews). In our case, we use multi-dimensional high resolution shock-capturing methods, with a finite volume scheme. Time advance of the equations is done via dimensional splitting to compute the numerical fluxes across cell boundaries, using the integral form of the conservation equations (see e.g. Komissarov Reference Komissarov1999) and Runge–Kutta methods (Shu & Osher Reference Shu and Osher1989).
In the case of RHD simulations, we use our code ratpenat, a hybrid parallel code – MPI $+$ OpenMP – (see Perucho et al. (Reference Perucho, Martí, Cela, Hanasz, de La Cruz and Rubio2010b) and references therein) in which: (i) primitive variables within numerical cells are reconstructed using piecewise-parabloic method (PPM) routines, (ii) numerical fluxes across cell interfaces are computed with the Marquina flux formula, (iii) advance in time is performed with third-order total variation diminishing (TVD)-preserving Runge–Kutta methods.
The RMHD simulations are run with the code lóstrego (López-Miralles et al. Reference López-Miralles, Perucho, Martí, Migliari and Bosch-Ramon2022), which uses the same basic strategy, although with different schemes: (i) primitive variables are reconstructed using piecewise linear methods: MinMod (Roe Reference Roe1986), monotonized central (MC) (van Leer Reference van Leer1977) or VanLeer (van Leer Reference van Leer1974), with slope limiters that preserve monotonicity (Suresh & Huynh Reference Suresh and Huynh1997), (ii) the approximate solvers used to compute the fluxes are Harten-Lax-Van Leer (HLL) (Harten, Lax & van Leer Reference Harten, Lax and van Leer1983), HLLC (Mignone & Bodo Reference Mignone and Bodo2006) and HLLD (Mignone, Ugliano & Bodo Reference Mignone, Ugliano and Bodo2009). In this case, the divergence-free condition of the field is maintained by using the constrained transport technique (Evans & Hawley Reference Evans and Hawley1988; Ryu et al. Reference Ryu, Miniati, Jones and Frank1998; Balsara & Spicer Reference Balsara and Spicer1999; Gardiner & Stone Reference Gardiner and Stone2005).
These codes are run on supercomputing resources at the Universitat de València (Tirant, Vives) or within the Spanish Supercomputing Network (Red Española de Supercomputación, RES). The parallelization used in both codes permits a three-dimensional splitting of the grid for distributed memory architectures. In this contribution, we present simulations which typically use 512–1024 cores and require ${>}100\,{\rm khr}$ computing hours.
3. Applications
3.1. The RHD simulations
3.1.1. Long-term evolution
One line of our work has been aimed at understanding the evolution of powerful jets through the density/pressure profile of their host galaxies, and what role they play in heating the intergalactic medium. Our simulations were set up as an initial grid covering ${\sim }100\unicode{x2013}200\,{\rm kpc}$ in the three-dimensional case (Perucho et al. Reference Perucho, Martí and Quilis2019; Perucho, Martí & Quilis Reference Perucho, Martí and Quilis2022) to ${\sim }1\,{\rm Mpc}$ in two-dimensional simulations (Perucho, Quilis & Martí Reference Perucho, Quilis and Martí2011; Perucho et al. Reference Perucho, Martí, Laing and Hardee2014). The jets are injected as a boundary condition, with properties expected for jets at the injection point, typically 1 kpc from the central formation engine. The ambient density at this distance is ${\sim }10^{5}\,{\rm m}^{-3}$, and it drops with distance following a profile derived from modelling X-ray observations of radio galaxy 3C 31 (Hardcastle et al. Reference Hardcastle, Worrall, Birkinshaw, Laing and Bridle2002).
At these scales, the jet is expected to be dominated by kinetic energy flux, and the magnetic field to have a dominant disordered component, thus only contributing to pressure, but triggering minor tension forces. This is an assumption that we apply to simulate these scenarios with RHD codes.
The jets are injected into the grid with a mildly relativistic velocity $0.9\,c$ and specific enthalpies around $c^2$. All these simulations show that collimated relativistic jets generate high-pressure regions of shocked jet gas. The high sound speed in this region facilitates equilibrium within the whole shocked region, feeding pressure-driven shock expansion. Altogether, this is a very efficient process in terms of energy transfer from the jet to the ambient medium, as explained in Perucho et al. (Reference Perucho, Martí, Quilis and Borja-Lloret2017b). The main point that the simulations revealed is that, because the post-shock pressure depends on the momentum flux density, collimated relativistic jets transfer most of their energy flux (which is not dominated by the rest-mass energy of the particles, as happens in classical jets) to the ambient medium via shock heating.
Jet collimation is maintained due to its propagation through a decreasing density/pressure environment, which favours a faster decrease in the pressure of the shocked region, allowing (i) the jet to expand with small opening angles, and (ii) the increase of the relative value of jet-to-ambient inertia, which limits the growth rates of unstable modes (Hardee Reference Hardee1982). Furthermore, the jet advance velocity through such an ambient medium also increases with distance as compared with a homogeneous constant density jet. All these factors also contribute to a decrease of the growth rates of unstable modes (Perucho, Martí & Hanasz Reference Perucho, Martí and Hanasz2005; Perucho Reference Perucho2019), and to faster expansion, thus allowing the jet to develop to large scales while keeping its collimation.
In summary, our simulations have shown that powerful relativistic outflows ($L_j \geq 10^{38}\,{\rm W}$) propagating through galactic atmospheres are extremely efficient in heating the interstellar and intergalactic gas because of (i) their energy flux is not dominated by rest-mass energy of the particles, and (ii) they can keep their collimation.
Figures 1 and 2 show the results from two of the described simulations, with different ambient medium core densities: $10^{-7}\,{\rm m}^{-3}$ in the former and $4\times 10^{-7}\,{\rm m}^{-3}$ in the latter (Perucho et al. Reference Perucho, Martí and Quilis2019, Reference Perucho, Martí and Quilis2022). The images show cuts of rest-mass density and pressure weighted with tracer and a rendering of jet mass fraction (tracer) that show the jet structure for the last snapshot of each simulation at times ${\simeq }5.4\,{\rm Myr}$ and ${\simeq }8\,{\rm Myr}$, respectively. In both cases, the jets generate thin backflow/cocoon regions, as caused by the fast propagation through the density decreasing galactic atmosphere. The increased ambient density triggers slightly inflated backflow regions around the jet's terminal shock.
Figure 3 shows projected pseudosynchrotron emissivity images (at viewing angles of $15^\circ$) of the three-dimensional RHD simulations displayed in figures 1 and 2, from a dynamical point of view. The expression used to compute the emissivity at a given frequency is the same as that used by Hardee (Reference Hardee2003) (see also Clarke, Norman & Burns Reference Clarke, Norman and Burns1989)
where $n$ is the particle number density, $p$ is pressure, $B$ is the assumed magnetic field strength, taken to be proportional to $\sqrt {p}$, $\theta _B$ is the angle between the field lines and the viewing angle assumed to be a constant, $D$ is the Doppler factor and $\alpha$ is the spectral index (defined as $S_\nu \propto \nu ^{-\alpha }$), taken as constant and equal to 0.5 (see, e.g. Heavens & Drury Reference Heavens and Drury1988; Kirk & Duffy Reference Kirk and Duffy1999). This expression represents a raw approximation, with strong simplifications as that relative to the magnetic field structure or spectral indeces. This assumption cancels differences in emissivity caused by changes in $\theta _B$, but this is expected to be a second-order contribution.Footnote 2 Therefore, the approach is enough to show relevant features. Furthermore, a detailed distribution of the magnetic field across the grid would require a fully RMHD simulation, which is left for future simulations of AGN jets (see § 3.2). The images show the effect of the Doppler factor in generating the jet-to-counter-jet brightness asymmetry and reproduce the canonical FRII jet morphology (see, e.g. Bridle Reference Bridle1984; Bridle & Perley Reference Bridle and Perley1984; Harwood et al. Reference Harwood, Croston, Intema, Stewart, Ineson, Hardcastle, Godfrey, Best, Brienza and Heesen2016), with emissivity enhancements at recollimation shocks, and bright hotspots surrounded by radio lobes.
3.1.2. Mass load
Less powerful jets ($L_j \leq 10^{36}\,{\rm W}$) can be decelerated by mass entrainment (Bicknell Reference Bicknell1984). This process may take place in different ways, but mainly through the development of (Kelvin–Helmholtz, current driven, centrifugal…) instabilities at the jet boundary (e.g. Meliani & Keppens Reference Meliani and Keppens2007; Perucho & Martí Reference Perucho and Martí2007; Matsumoto & Masada Reference Matsumoto and Masada2013; Massaglia et al. Reference Massaglia, Bodo, Rossi, Capetti and Mignone2016, Reference Massaglia, Bodo, Rossi, Capetti and Mignone2019, Reference Massaglia, Bodo, Rossi, Capetti and Mignone2022; Matsumoto, Aloy & Perucho Reference Matsumoto, Aloy and Perucho2017; Gourgouliatos & Komissarov Reference Gourgouliatos and Komissarov2018) or by interactions between the jet and stars or clouds that cross it while orbiting the galactic centre (e.g. Komissarov (Reference Komissarov1994), Bowman, Leahy & Komissarov (Reference Bowman, Leahy and Komissarov1996), Laing & Bridle (Reference Laing and Bridle2002), Hubbard & Blackman (Reference Hubbard and Blackman2006), Perucho et al. (Reference Perucho, Martí, Laing and Hardee2014), Perucho (Reference Perucho2019), Anglés-Castillo et al. (Reference Anglés-Castillo, Perucho, Martí and Laing2021) and references therein).
Modelling of FRI jets has revealed that the region in which jets seem to become symmetric in brightness with respect to their counter-jets is also a region in which they become bright and that deceleration takes place progressively, from the jet boundaries towards their axes (Laing & Bridle Reference Laing and Bridle2014). This indicates that deceleration is driven by small-scale processes that take place at the boundaries and points towards either small wavelength instability modes arising at the shear transition with the jet environment or towards other kinds of perturbative interactions (Perucho et al. Reference Perucho, Bosch-Ramon and Barkov2017a; Perucho Reference Perucho2020).
In the case of clouds and stars interacting with the jet, a bow shock is formed against the jet flow, and the mixing and loading process takes place at the cometary tail formed downstream of the interaction site (Komissarov Reference Komissarov1994). It has been suggested that this process alone could decelerate low power jets ($L_j \leq 10^{35}\,{\rm W}$, Bowman et al. Reference Bowman, Leahy and Komissarov1996; Perucho et al. Reference Perucho, Martí, Laing and Hardee2014) and also change the properties (composition, magnetization…) in intermediate power jets ($L_j \sim 10^{36}\,{\rm W}$, Anglés-Castillo et al. Reference Anglés-Castillo, Perucho, Martí and Laing2021). These changes would be forced, according to Anglés-Castillo et al. (Reference Anglés-Castillo, Perucho, Martí and Laing2021), by the relevant relative contribution of the rest-mass density of protons with respect to the jet's injected electron–positron pairs, on the one hand, and the corresponding increase in kinetic energy flux and gas pressure, which reduces the dynamical relevance of the magnetic field, on the other hand.
Furthermore, these interactions can host particle acceleration and generate the production of high-energy radiation and, altogether, may produce the brightness flaring in decelerating jets (e.g. Bednarek & Protheroe Reference Bednarek and Protheroe1997; Bosch-Ramon, Perucho & Barkov Reference Bosch-Ramon, Perucho and Barkov2012b; Wykes et al. Reference Wykes, Croston, Hardcastle, Eilek, Biermann, Achterberg, Bray, Lazarian, Haverkorn and Protheroe2013, Reference Wykes, Hardcastle, Karakas and Vink2015; Vieyro, Torres-Albà & Bosch-Ramon Reference Vieyro, Torres-Albà and Bosch-Ramon2017; Torres-Albà & Bosch-Ramon Reference Torres-Albà and Bosch-Ramon2019; Perucho Reference Perucho2020).
Numerical simulations have allowed us to characterize these scenarios and estimate the radiative output from them at different jet scales, e.g. parsecs and hundreds of parsecs (Bosch-Ramon et al. Reference Bosch-Ramon, Perucho and Barkov2012b; Perucho et al. Reference Perucho, Bosch-Ramon and Barkov2017a). Future work should be focused on understanding the role of magnetic fields in these collisions and the study of particle acceleration and cosmic ray production (e.g. Wykes et al. Reference Wykes, Croston, Hardcastle, Eilek, Biermann, Achterberg, Bray, Lazarian, Haverkorn and Protheroe2013; Murase, Oikonomou & Petropoulou Reference Murase, Oikonomou and Petropoulou2018).
3.1.3. Jet–ISM interactions
Another aspect that has improved during the last decade in the field of AGN jet simulations is the characterization of galactic ambient media, both inside and outside the host galaxies (see, e.g. Wagner, Bicknell & Umemura Reference Wagner, Bicknell and Umemura2012; Mukherjee et al. Reference Mukherjee, Bicknell, Sutherland and Wagner2016; Bicknell et al. Reference Bicknell, Mukherjee, Wagner, Sutherland and Nesvadba2018; Mandal et al. Reference Mandal, Mukherjee, Federrath, Nesvadba, Bicknell, Wagner and Meenakshi2021; Perucho et al. Reference Perucho, López-Miralles, Reynaldi and Labiano2021). Galactic environments have started to be set up by using distributions obtained from cosmological simulations (e.g. Heinz et al. Reference Heinz, Brüggen, Young and Levesque2006; Yates-Jones, Shabala & Krause Reference Yates-Jones, Shabala and Krause2021). Within host galaxies, the two-phase nature of the ISM, with cold denser gas hosted in clouds and hot more dilute gas between them represents a non-negligible property to consider. This is especially relevant within the inner kiloparsec of host galaxies, where the broad and narrow line regions are located.
In addition, observational results at radio-to-UV ranges show that the jets trigger line emission and gas motions in their host galaxies (Morganti & Oosterloo Reference Morganti and Oosterloo2018). Based on previous work by Vaidya et al. (Reference Vaidya, Mignone, Bodo and Massaglia2015), we have included hydrogen ionization and recombination feedback effects in our RHD code with the aim of tracing the effect of jets on the cold gas in the inner kiloparsec (Perucho et al. Reference Perucho, López-Miralles, Reynaldi and Labiano2021). We have now improved the set up, using realistic ambient/cloud densities and temperatures. In these simulations, a purely leptonic ($e^-/e^+$) jet is injected into an inhomogeneous medium composed of a mixture of clouds that host atomic hydrogen, at temperatures ${\sim }100\,{\rm K}$ and maximum density $n\sim 10^9\,{\rm m}^{-3}$, and an ionized medium with $T\sim 10^6\,{\rm K}$ and $n\sim 10^{5}\,{\rm m}^{-3}$. The numerical box reproduces the inner 500 pc in the host galaxy, with the jet injected as a boundary condition at a certain distance from the forming region with relativistic velocity $0.98\,c$, density $\rho _j=1.67\times 10^{-26}\,{\rm kg}\,{\rm m}^{-3}$ and a very high specific internal energy $\varepsilon \sim 10^3\,c^2$ to simulate the high pressure expected in this region (the simulation does not include a magnetic field so this approach could only account for a magnetic pressure generated by a disordered field configuration). In the simulations, the ambient medium is also given a transverse velocity of 100 km s$^{-1}$ to emulate the rotation around the galactic nucleus.
Altogether, the parameter range in the simulations spans more than eight orders of magnitude in density and pressure, which makes the simulation extremely challenging. Moreover, time-step limiters need to be included to allow for the short hydrogen recombination rates in cold regions or ionization in shocked cells. The simulation uses two equations of state: a relativistic (Synge Reference Synge1957) and a non-relativistic (see Vaidya et al. Reference Vaidya, Mignone, Bodo and Massaglia2015) ideal gas. The criterion to choose between them in each cell is based on the composition, namely, on the fact that neutral hydrogen is present or not in the cell. The top panels of figure 4 show cuts of rest-mass density (left, in code units $\rho _a=10^7\,{\rm m}_{\rm p}\,{\rm m}^{-3}$, with an upper limit to the colour scale set at $\rho _a=10^8\,{\rm m}_{\rm p}\,{\rm m}^{-3}$), pressure (centre, in code units $\rho _a c^2=1.5\times 10^{-3}\,{\rm Pa}$, with a lower limit set at $\rho _ac^2=1.5\times 10^{-8}\,{\rm Pa}$) and velocity field (right, in code units $c$, limited at $10^{-5}\,c$) at the last snapshot obtained for the simulation, which is still being run, at $t=460\,{\rm yr}$ after injection. These plots reveal a complex shock structure around the jet as forced by the strong inhomogeneity in density found in the ambient medium. The shocked region is also highly inhomogeneous in density, revealing the richness in the interaction between the jet and the ISM shocked gas. The high sound speed contributes to homogenizing the pressure within the shocked region with the exception of some knots of denser (ambient) gas, and the jet's hotspot. Finally, the velocity field shows a wide range of values, from $10^{-4}$ to $0.1\,c$. Obviously, the smaller velocities correspond to the denser regions (see the comparison with the density panel). Although there are regions with values that can fall to several tens of km/s (blue colour within the shocked area, ${\sim }10^{-4}\,c$), the dominating blue–white transition implies typical velocities of the order of hundreds to a thousand km/s. These values are in agreement with the typical ones measured from line emission in jetted active galaxies (see, e.g. Morganti & Oosterloo Reference Morganti and Oosterloo2018; Schulz et al. Reference Schulz, Morganti, Nyland, Paragi, Mahony and Oosterloo2021).
The bottom left panel in figure 4 shows a limited box centred on the region occupied by the jet structure with the rendering of a tracer that indicates the regions where originally dense cloud atomic gas can be found. From the image, it is evident that this gas is fixed in clouds outside the shocked volume, but it is completely disrupted and mixed by the shock, and appears concentrated towards the shock region. Finally, the bottom right panel displays, for the same box, a collection of isosurfaces of pressure (red) to highlight the shock wave (red, $7\times 10^{-7}\rho _a c^2 \simeq 10^{-10}\,{\rm Pa}$), leptonic number (dark blue, $\rho _{e^-/e^+}/\rho$ at $10^{-3}$, as compared with ${\simeq }5.4\times 10^{-4}$ for $e^-/p$ gas) to show the jet and still lepton dominated unmixed regions, atomic hydrogen density outside the shock (orange, for $10^8\,{\rm m}_{\rm p}\,{\rm m}^{-3}$) to show the cloud distribution and inside the shock (bluish, for $10^5\,{\rm m}_{\rm p}\,{\rm m}^{-3}$) to show hints of atomic hydrogen within the shocked cavity. The low density (bluish) hydrogen density contours follow clearly the shock surface and thus reveal a very fast and complete ionization of the atomic gas. Therefore, we observe that the shocked gas is completely ionized, so despite our simulation recovering the observed velocities, as stated above, it also shows that the shocks produced by powerful jets ionize atomic (and therefore molecular) gas very rapidly and this forbids line detection, unlike what the observations clearly reveal. The origin of this discrepancy probably lies in the fact that observed lines correspond to radio sources more evolved than our simulation at the current position (several kpc vs 200 pc) and that post-shock temperatures (typically $10^6\unicode{x2013}10^7\,{\rm K}$) are too high to allow for recombination of shock-ionized hydrogen within the simulated time. Nevertheless, we observe white spotted areas in the cuts (top left panel) which imply densities between $10^7$ and $10^8\,{\rm m}_{\rm p}\,{\rm m}^{-3}$ and this could translate in cooling times of ${\sim }10^{3\unicode{x2013}4}\unicode{x2013}10^{4\unicode{x2013}5}\,{\rm yr}$ for gas between solar and 0 metallicity at the aforementioned temperatures (Dopita & Sutherland Reference Dopita and Sutherland2003). Although the simulations are limited to pairs, and atomic+ionized hydrogen in the composition of the gas, they show that the conditions for line emission with the observed velocities could be recovered once the jet has evolved for at least twice the simulated time so far. Therefore our simulations represent a promising first step in the understanding of the ISM dynamics driven by jet injection in the host galaxies.
3.2. The RMHD simulations
3.2.1. Jet–stellar wind interactions
Jets in X-ray binaries propagate through the binary environment before carving their way towards the ISM. In this region, especially in the case of massive binaries, jets travel through a dense, slow stellar wind that can produce a strong lateral impact on the jet flow, triggering internal shocks and, in the case of low power jets, even disruption. This was proposed as a plausible scenario for triggering gamma-ray radiation and thus explaining the emission physics of the very few gamma-ray emitting binaries observed in the Galaxy (see e.g. Rieger, Bosch-Ramon & Duffy Reference Rieger, Bosch-Ramon and Duffy2007; Bosch-Ramon & Khangulyan Reference Bosch-Ramon and Khangulyan2009b and references therein).
Following a number of papers that studied RHD simulations of jet–massive wind interactions (Perucho & Bosch-Ramon Reference Perucho and Bosch-Ramon2008; Perucho et al. Reference Perucho, Bosch-Ramon and Khangulyan2010a; Perucho & Bosch-Ramon Reference Perucho and Bosch-Ramon2012), we run a set of dedicated three-dimensional numerical simulations using our new code lóstrego, including the dynamical effects of magnetic fields (López-Miralles et al. Reference López-Miralles, Perucho, Martí, Migliari and Bosch-Ramon2022). The results of these simulations showed that the latter could play a stabilizing role – and thus, to provide extra collimation within the binary scale – in the jet evolution when the flux of magnetic energy is low compared with the flux of kinetic energy, while a non-negligible magnetic energy flux makes the jet prone to the development of current-driven instabilities within the scale of the binary.
The dynamics of the jet model with a dynamically relevant field in López-Miralles et al. (Reference López-Miralles, Perucho, Martí, Migliari and Bosch-Ramon2022) is particularly interesting from several points of view. The jet evolution during the interaction with the stellar wind is significantly different with respect to a hydrodynamical jet with the same total power: whereas a powerful hydrodynamical jet evolves without losing collimation, the magnetized jet is prone to the development of current-driven instability pinching and kink modes triggered by overpressure and the toroidal field. Therefore, the jet shows a chain of recollimation shocks within the binary scale, followed by the development of a strong kink that spreads the jet momentum throughout a large area and decelerates its advance. The magnetic field structure is toroidal from injection to the development of the kink, while it becomes highly entangled near the head, filling the cavity.
In figure 5 we present the resulting configuration of the magnetic field in this simulation by focusing on the field lines. This plot shows relevant information about the field structure in the jet spine, which was hidden by the mixed plasma in the cocoon in our representation in López-Miralles et al. (Reference López-Miralles, Perucho, Martí, Migliari and Bosch-Ramon2022). Since the magnetic field within the cocoon is low compared with that in the jet, we have applied a linear transparency to the field lines, weighted by the field vector module. In the mid–low part of the jet, the magnetic field preserves the toroidal structure of injection, which is reinforced at the recollimation shocks, as expected. The jet tracer three-dimensional render (top panel), which is limited to $f>0.9$ to highlight the jet core, shows a pinched morphology and a set of annuli of jet plasma surrounding the core, which are deposited in the cavity due to the jet dynamics during the first stages of propagation. This morphology is very particular and could give rise to different emission patterns that deserve to be studied in detail, but this is out of the scope of this contribution. In the mid–upper part of the simulation, both the new representation of the field and the tracer function allow us to distinguish a well-resolved precessing morphology, triggered by the development of a kink instability. The inset plot of figure 5 zooms in on the jet head, showing that the field is highly reinforced at the elbow of the twisted trace, which is also directly impacting the shocked environment formed by the stellar wind.
Further analysis will be required to investigate the consequences of this precessing morphology on the jet dynamics (within and beyond the binary), on the one hand, and its implications in terms of non-thermal emission, on the other. This kink could contribute to the observed signatures of precessing jets that have been resolved at radio wavelengths for different X-ray binaries with jet-like structures. The most relevant example of jet precession is the well-known microquasar SS 433, but there are many other examples (see e.g. Mioduszewski et al. Reference Mioduszewski, Rupen, Hjellming, Pooley and Waltman2001; Massi, Ros & Zimmermann Reference Massi, Ros and Zimmermann2012; Miller-Jones et al. Reference Miller-Jones, Tetarenko, Sivakoff, Middleton, Altamirano, Anderson, Belloni, Fender, Jonker and Körding2019; Luque-Escamilla, Martí & Martínez-Aroza Reference Luque-Escamilla, Martí and Martínez-Aroza2020). The origin of precession in microquasars is still debated for most of these sources. Although the classical theoretical models tend to explain precession by invoking relativistic effects of the inner disc (i.e. Lense–Thirring precession, Liska et al. Reference Liska, Hesp, Tchekhovskoy, Ingram, van der Klis and Markoff2018; Motta et al. Reference Motta, Franchini, Lodato and Mastroserio2018) or the coupled effect of stellar winds and orbital motion (Barkov & Bosch-Ramon Reference Barkov and Bosch-Ramon2022), our simulation shows that current-driven instabilities can also trigger helical patterns during the jet–wind interaction, although the periods associated with this kind of precession should be accurately related. For example, current-driven instabilities have been claimed to explain the large-scale morphology of the jet structures observed in GRS 1758-258 (Luque-Escamilla et al. Reference Luque-Escamilla, Martí and Martínez-Aroza2020).
This simulated jet thus represents an extraordinary candidate to produce high-energy radiation in microquasars, since energy dissipation may not only occur at the several reconfinement shocks where particles can be accelerated, but also a kink instability that can lead to kink-driven magnetic reconnection processes further out from the injection base (see e.g. Bodo, Tavecchio & Sironi Reference Bodo, Tavecchio and Sironi2021; Bodo et al. Reference Bodo, Mamatsashvili, Rossi and Mignone2022). The characterization of the periodic behaviour of these processes (i.e. shocks and kink-driven precession) can be also relevant to interpreting the rapid X-ray variability of some X-ray binaries/microquasars. Therefore, future work will include consideration of those dissipative processes in the acceleration of non-thermal particles and the triggering of high-energy radiative output. Although magnetic reconnection is not reproduced in ideal RMHD codes like ours, recipes such as those presented in Bodo et al. (Reference Bodo, Tavecchio and Sironi2021) will allow us to characterize the formation of current sheets and magnetic dissipation regions into our simulation. We expect that the combination of such analysis with numerical simulations such as those presented in this contribution, will provide new insights into the physical processes that lead to non-thermal emission in microquasars jets within the scale of the binary.
4. Conclusions
In this paper, we have presented a number of scenarios that require the use of numerical codes to study due to their extreme nature, which makes this research prohibitive for current laboratory experimental capabilities. All these scenarios have in common the need for relativistic (magneto-)hydrodynamics, and they use similar methods to solve the system of equations that is used to model them. Different improvements have been made in the codes used in our group through the past years, including optimization and parallelization, modelling of ambient media, etc. These improvements have allowed us to tackle more realistic simulations, which, although still far from reproducing the studied scenarios accurately, have brought about relevant advances in our knowledge of the physics of relativistic outflows in astrophysics and their impact on their environment. Future work should address in detail the role of magnetic fields in the evolution and impact of jets in the scenarios previously studied with RHD simulations while, at the same time, improving the physical processes that the codes take into account and that can be relevant for studying those scenarios. One such improvement consists of including the radiation dynamics, which could allow us to compare the role of shocks and radiative output from the central engines (accretion disks) in shaping the heating mechanisms of the ISM and the intergalactic medium (IGM/WHIM). Future improvements of codes will be related to linking multi-scale processes and adding relevant physics to improve the simulations of astrophysical scenarios that are still extremely challenging for current supercomputing capabilities and far from realistic for laboratory experiments.
Acknowledgements
M.P. acknowledges M. Perucho-Franco, his father, for standing up as an example through his whole life. The authors thank the referees for their constructive comments.
Editor Luís O. Silva thanks the referees for their advice in evaluating this article.
Funding
Computer simulations have been carried out in the Servei d'Informàtica de la Universitat de València (Tirant). This work has been supported by the Spanish Ministry of Science through Grants PID2019-105510GB-C31/AEI/10.13039/501100011033, PID2019-107427GB-C33, from the Generalitat Valenciana through grant PROMETEU/2019/071, it was also funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) – project number 443220636, and it forms part of the Astrophysics and High Energy Physics programme and was supported by MCIN with funding from European Union NextGenerationEU (PRTR-C17.I1) and by Generalitat.
Declaration of interests
The authors report no conflict of interest.
Author contributions
Both authors contributed to analysing data and reaching conclusions, and in writing the paper.