Accelerating AGN jets to parsec scales using general relativistic MHD simulations

TL;DR

This study uses high-resolution general relativistic MHD simulations to explore how black hole jets interact with their environment, revealing mechanisms for jet collimation, variability, and particle acceleration consistent with observations.

Contribution

It provides the first detailed GRMHD simulation of AGN jets spanning over five orders of magnitude, elucidating jet collimation, instabilities, and mass-loading processes.

Findings

Jets attain a parabolic shape similar to M87

Pinch instabilities cause jet variability and heating

Ambient medium squeezing enhances jet mass-loading

Abstract

Accreting black holes produce collimated outflows, or jets, that traverse many orders of magnitude in distance, accelerate to relativistic velocities, and collimate into tight opening angles. Of these, perhaps the least understood is jet collimation due to the interaction with the ambient medium. In order to investigate this interaction, we carried out axisymmetric general relativistic magnetohydrodynamic simulations of jets produced by a large accretion disc, spanning over 5 orders of magnitude in time and distance, at an unprecedented resolution. Supported by such a disc, the jet attains a parabolic shape, similar to the M87 galaxy jet, and the product of the Lorentz factor and the jet half-opening angle, $\gamma\theta\ll 1$, similar to values found from very long baseline interferometry (VLBI) observations of active galactic nuclei (AGN) jets; this suggests extended discs in AGN. We find that the interaction between the jet and the ambient medium leads to the development of pinch instabilities, which produce significant radial and lateral variability across the jet by converting magnetic and kinetic energy into heat. Thus pinched regions in the jet can be detectable as radiating hotspots and may provide an ideal site for particle acceleration. Pinching also causes gas from the ambient medium to become squeezed between magnetic field lines in the jet, leading to enhanced mass-loading of the jet and potentially contributing to the spine-sheath structure observed in AGN outflows.

Figures

Tables

Table 1. Table 1 : Simulations carried out for this work. Model names shorthand disc inner radius r in subscript 𝑟 in r_{\rm in} (R), black hole spin (A; 0.9375 0.9375 0.9375 , common for all models) and the b 2 / ρ superscript 𝑏 2 𝜌 b^{2}/\rho floor value (B; in other words, the maximum magnetisation σ 0 subscript 𝜎 0 \sigma_{0} ) with additional parameters including S for single loop magnetic field and LR and HR for the low and high resolution versions respectively. We also mention the pressure-maximum radius r max subscript 𝑟 max r_{\rm max} , outer grid radius r out subscript 𝑟 out r_{\rm out} , grid resolution, simulation time t sim subscript 𝑡 sim t_{\rm sim} , MRI quality factor Q x superscript 𝑄 𝑥 Q^{x} at t = 2 × 10 4 ​ t g 𝑡 2 superscript 10 4 subscript 𝑡 𝑔 t=2\times 10^{4}t_{g} (see text in Sec.  4.1 ) and magnetic field configuration (indicated by the number of field line loops) as well as the use of local adaptive time stepping (LAT).

Full Model	Short	$r_{\mathrm{in}}$	$r_{\max }$	$r_{\mathrm{out}}$	Floor	Resolution	$t_{\mathrm{sim}}$	$Q$ factor	Initial field loop	LAT
name	name	[$r_g$]	[$r_g$]	[$r_g$]	$b^2/\rho$	$N_r\times N_{\theta }$	[${10}^5{t}_g$]	$Q^r,Q^{\theta }$	configuration
R36A93B10	B10	36	73.97	${10}^6$	10	$6000\times 800$	5.0	62,14	multiple	on
R36A93B10-S	B10-S	36	73.97	${10}^5$	10	$4000\times 800$	1.2	96, 24	single	on
R25A93B10	B10-R	25	50	${10}^5$	10	$3000\times 800$	1.2	36, 10	multiple	on
R36A93B10-S-LR	B10-SLR	36	73.97	${10}^4$	10	$640\times 256$	0.8	20, 4	single	on
R36A93B10-S-LR-L	B10-SLRL	36	73.97	${10}^4$	10	$640\times 256$	0.8	20, 4	single	off
R36A93B10-HR	B10-HR	36	73.97	${10}^5$	10	$18000\times 1200$	0.5	-,-	multiple	on
R36A93B3	B3	36	73.97	${10}^5$	3	$4000\times 800$	1.5	32,10	multiple	on
R36A93B50	B50	36	73.97	${10}^5$	50	$4000\times 800$	1.9	32,10	multiple	on
R36A93B100	B100	36	73.97	${10}^5$	100	$4000\times 800$	2.0	40, 12	multiple	on

Equations

$$\text{centering}{A}_{\varphi }=\{\begin{array}{ll}(\rho -0.05{)}^2{r}^2,& \text{if r<300rg and ρ>0.05}.\\ 0,& \text{otherwise}.\end{array}\text{@add@centering}$$

$$\text{centering}{A}_{\varphi }=\{\begin{array}{ll}f,& \text{if r<300rg and ρ>0.05},\\ 0,& \text{otherwise},\end{array}\text{@add@centering}$$

$$\mu =\frac{F_{\mathrm{E}}}{F_{\mathrm{M}}}=\frac{(\rho +u_{\mathrm{g}}+p_{\mathrm{g}}+b^2)u^r{u}_t-b^r{b}_t}{-\rho u^r},$$

$$\text{centering}\sigma =\frac{F_{\mathrm{EM}}}{F_{\mathrm{K}}}=\frac{b^2{u}^r{u}_t-b^r{b}_t}{\rho u^r{u}_t},\text{@add@centering}$$

$$\frac{\gamma }{\mu }\approx 1-\frac{\pi B_p{R}^2}{\Phi }=1-a_{bp},$$

$$\text{centering}{\gamma }_f{v}_f={\left(\frac{b^2}{\rho }\right)}^{1/2},\text{@add@centering}$$

$$\text{centering}{\theta }_{\mathrm{M}}=\frac{1}{M_f}=\frac{{\gamma }_f{v}_f}{\gamma v}\stackrel{\text{eqref}eqn:fast}{=}\frac{\sqrt{b^2/\rho }}{\gamma v}\stackrel{\text{eqref}eqn:sigma}{\approx }\frac{\sqrt{\sigma }}{\gamma v}.\text{@add@centering}$$

$$\text{centering}\mu =\gamma (\sigma +h+1),\text{@add@centering}$$

$$\text{centering}1-\cos \theta ={\left(\frac{r+r_0}{r_{\mathrm{H}}+r_0}\right)}^{-\nu }.\text{@add@centering}$$

$$\text{centering}{A}_{\varphi }={\left(\frac{r+r_0}{r_{\mathrm{H}}+r_0}\right)}^{\nu }(1-\cos \theta ).\text{@add@centering}$$

$$\text{centering}s=\frac{1}{\Gamma -1}{\log }_{10}\left(\frac{p_{\mathrm{g}}}{{\rho }^{\Gamma }}\right),\text{@add@centering}$$

$$\text{centering}{\gamma }_f^2{v}_f^2\stackrel{\text{eqref}eqn:sigma}{\approx }\sigma \stackrel{\text{eqref}eqn:energy}{\approx }\frac{\mu }{\gamma }-1,\text{@add@centering}$$

$$\text{centering}{\gamma }_f\approx {\mu }^{1/3}.\text{@add@centering}$$

$$\text{centering}{r}_f\approx {\left(\frac{{\mu }^{1/3}}{C\Omega }\right)}^{1/(1-\zeta )}.\text{@add@centering}$$

$$t$$

$$r$$

$$\theta$$

$$\varphi$$

Full text

††pubyear: 2019††pagerange: Accelerating AGN jets to parsec scales using general relativistic MHD simulations–C

1 Introduction

Powered by magnetic fields brought inwards by infalling gas, relativistic jets form in a variety of astrophysical black hole systems such as active galactic nuclei (AGN), black hole X-ray binaries (XRBs), tidal disruption events (TDEs), and neutron star mergers. Through the exchange of energy with the ambient medium, jets heat up gas in the interstellar medium (ISM) creating a feedback loop between the AGN and its environment (e.g., Bower et al., 2006; Fabian, 2012). AGN feedback plays a key role in regulating the growth of galaxies and star formation (e.g., Silk & Rees 1998; Magorrian et al. 1998; for a recent review, see Harrison et al. 2018) and therefore, its implementation in cosmological simulations (e.g., Springel et al., 2005; Anglés-Alcázar et al., 2017; Weinberger et al., 2018) warrants an accurate understanding of jet energetics (e.g., Sijacki et al., 2007; Bourne & Sijacki, 2017). Addressing how jets build up the magnetic and kinetic energy required to explain the radiative emission seen from jet observations is still an open problem. Therefore, by studying how jets accelerate and interact with the ambient gas, we can better understand jet emission and its relation to jet-ISM interactions.

Very long baseline interferometry (VLBI) imaging of AGN jets (e.g., Lister et al., 2016) has made it possible to track radio emission features of the jet over long time periods and estimate the bulk kinematic properties. As an example, Asada & Nakamura (2012) show that the parsec scale jet of the galaxy M87 is roughly parabolic in structure. At $\simeq 10^{5}$ gravitational radii away from the central black hole, the jet becomes conical, with the transition appearing to coincide with the bright HST-1 feature (Biretta et al., 1999). The HST-1 knot may be a result of self-collimation (e.g., Polko et al., 2010) or a changing density profile of the ISM (Asada & Nakamura, 2012; Nakamura & Asada, 2013). The change in the confining pressure may cause the jets to over-collimate and activate magnetic instabilities or result in the formation of internal shocks, leading to particle acceleration and the appearance of knots, stationary or moving features (Tchekhovskoy & Bromberg, 2016; Barniol Duran et al., 2017). The same phenomena also seem to occur in the jets of stellar mass black holes in XRBs (e.g., Markoff et al., 2001, 2005; Russell et al., 2013; Romero et al., 2017). Accelerated particles produce high energy emission, an important observational probe of these outflows. Due to the intricate relationship between jet structure, kinematics and particle acceleration, a variety of theoretical outflow models have been constructed to meaningfully interpret jet observations (for a review, see Meier, 2012).

In spite of the large body of theoretical work on outflows, understanding their dynamics remains a challenge, particularly because the jet structure depends on the interaction of the outflow with the ambient medium in a complex and non-linear way. Understanding this coupling at any scale requires the knowledge of the magnetic field configuration as well as the physical conditions at the jet base and the ambient medium, since both constrain the jet’s final energy content and Lorentz factor. However, constructing analytic models of jet acceleration is difficult due to the highly non-linear nature of the governing equations. Idealised semi-analytic models, or SAMs (e.g., Vlahakis & Königl, 2003a; Beskin & Nokhrina, 2006; Broderick & Loeb, 2006; Lyubarsky, 2009; Pu et al., 2015), provide reasonable estimates of jet properties and have been used to constrain the behaviour of the jets both near and far away from the event horizon.

Both the black hole and the surrounding accretion disc can launch outflows via two popular mechanisms. The Blandford-Znajek (BZ77; Blandford & Znajek, 1977) mechanism taps into the rotational energy of the black hole, via the magnetic field lines connected to the black hole event horizon, and leads to relativistic Poynting flux dominated jets (Beskin & Kuznetsova, 2000; Komissarov, 2001; Gammie et al., 2003; Komissarov, 2005; McKinney, 2005; Tchekhovskoy et al., 2010b). Field lines anchored in the accretion disc can launch sub-relativistic mass dominated winds by means of the Blandford-Payne mechanism (BP82; Blandford & Payne 1982, see also Meier et al. 1997; Mizuno et al. 2004). While SAMs can explain the basic physics of energy conversion from magnetic to kinetic form in jets, they are time-independent solutions that often neglect accretion disc physics as well as general relativistic effects of the space-time geometry around a (spinning) black hole. Due to these simplifications, SAMs are not able to completely capture the complexities of real jets, which is why general relativistic magneto-hydrodynamic (GRMHD) models are required (e.g., De Villiers et al., 2003; McKinney, 2006). GRMHD simulations typically model accretion starting from an initial gas torus, which, when threaded with large scale vertical or toroidal magnetic flux, launches both powerful BZ77- and BP82-type outflows (e.g., McKinney, 2005; Hawley & Krolik, 2006; Tchekhovskoy et al., 2011; Liska et al., 2018a).

There is currently a disagreement in the literature regarding jet acceleration. While semi-analytic models (e.g., Beskin et al., 1998; Beskin & Nokhrina, 2006) and idealised jet simulations (i.e. the simulations that model the ambient medium by placing a conducting wall at the jet’s outer boundary; e.g., Komissarov et al., 2007, 2009; Tchekhovskoy et al., 2010b) found efficient acceleration of jets to nearly the maximum Lorentz factor by utilising the jet’s entire energy budget, similar high efficiencies have never been seen in GRMHD simulations of a disc-jet system (e.g., McKinney, 2006; Bromberg & Tchekhovskoy, 2016; Barniol Duran et al., 2017), where internal pinch/kink instabilities are seen to convert a considerable amount of the jet’s energy content to heat by dissipating magnetic energy (e.g., Eichler, 1993; Spruit et al., 1997; Begelman, 1998; Giannios & Spruit, 2006). Does this mean that realistic systems are incapable of producing efficiently accelerating jets? In order to answer this question, we need to evolve jets over large time and distance scales since jet acceleration is an inherently time-dependent process in the presence of an ambient medium.

In this work, we revisit the problem of jet acceleration and collimation using a new state-of-the-art, GPU-accelerated GRMHD code H-AMR (Liska et al., 2018c) to carry out high-resolution axisymmetric disc-jet simulations. This was not possible for a long time, since as jets collimate, their magnetic field lines bunch up towards the jet axis, necessitating an extremely high resolution in the polar region. However, due to several algorithmic improvements, we can produce simulations spanning more than 5 orders of magnitude in both time and distance with unparalleled resolutions, presenting us with a unique opportunity to understand the jet physics. In Sec. 2, we give an overview of our problem setup. In Sec. 3, we describe our initial conditions. In Sec. 4, we present our results for disc-jet models. In Sec. 5, we compare our disc-jet models with an idealised one. In Sec. 6, we discuss our results and we conclude in Sec. 7.

2 Numerical setup

We use the H-AMR code (Liska et al., 2018a, b, c) that builds upon harmpi111Freely available at https://github.com/atchekho/harmpi and HARM2D (Gammie et al., 2003; Noble et al., 2006) and evolves the GRMHD equations on a fixed spacetime. It uses a Harten-Lax-van Leer (HLL) Riemann solver (Harten, 1983) to calculate fluxes at cell faces and a staggered grid akin to Gardiner & Stone (2005) to evolve the magnetic fields. H-AMR performs third order accurate spatial reconstruction at cell faces from cell centres using a piece-wise parabolic method (PPM; Colella & Woodward, 1984) and is second order accurate in time. The novelty of H-AMR lies in the use of advanced features such as adaptive mesh refinement (AMR, not utilised in this work) and a local adaptive time-step (LAT). The LAT reduces the number of conserved to primitive variable inversions and thereby increases the accuracy of the simulation (Appendix A); additionally, it speeds up the code by a factor of $\sim 3{-}5$. In its current version, in addition to CPUs, H-AMR also runs on graphical processing units (GPUs), achieving $10^{8}$ zone cycles per second on an NVIDIA Tesla V100 GPU and shows excellent parallel scaling to thousands of GPUs.

2.1 Numerical grid

We use units such that $G=M=c=1$. This sets both the characteristic timescale, $t_{g}=GM/c^{3}$, and spatial scale, $r_{g}=GM/c^{2}$, to unity. In fact, our simulations are scale-free, i.e., if we provide the black hole mass $M$ and mass accretion rate $\dot{M}$, we can rescale our simulation to the corresponding black hole system. Our grid is axisymmetric, extending from $0.85r_{\rm H}$ to $10^{5}r_{g}$, where $r_{\rm H}$ is the event horizon radius, $r_{\rm H}=r_{g}(1+\sqrt{1-a^{2}})$, where we set the dimensionless black hole spin parameter to $a=0.9375$. We carry out the simulations on a uniform grid in internal coordinates ($x^{\mu}$, see Appendix B) that are transformations of the spherical polar coordinates ($t,r,\theta,\phi$) in the Kerr-Schild foliation specifically optimised to follow the collimating jets. To resolve the jets, we typically use a numerical resolution of $3,000-18,000$ cells in the radial direction and $800$ cells in the polar direction (Table 1). We use the following boundary conditions (BCs): in $r-$, we use outflow BCs at the inner and outer grid radii; at the poles, we reflect the $\theta-$ component of the velocity and magnetic field.

2.2 Density floors

In the jet funnel matter either falls towards the black hole due to gravity or gets flung out due to magnetic forces depending on its location with respect to the stagnation surface, at which the inward pull of gravity balances the outward centrifugal force. The vacuum region thus created at this surface is a common numerical issue for grid-based MHD code, as gas density drops too low to be handled accurately. To avoid this, we replenish gas density and internal energy if they fall too low. We may physically motivate our floors as approximating the poorly understood processes leading to particle creation around the stagnation surface. Namely, mass flow divergence around the stagnation surface may lead to charge separation followed by particle creation (e.g., Hirotani & Okamoto, 1998; Broderick & Tchekhovskoy, 2015; Hirotani & Pu, 2016; Ptitsyna & Neronov, 2016; Levinson & Segev, 2017; Chen et al., 2018; Parfrey et al., 2019). In our simulations, we adopt the approach of Ressler et al. (2017) and mass-load the jets in the drift frame of the magnetic field. In this method, the component of the fluid momentum along the magnetic field is conserved. Our floor model consists of setting a minimum rest mass density $\rho_{\rm fl}=\max[b^{2}/\sigma_{0},2\times 10^{-4}(r/r_{g})^{-2.5}]$ and a minimum internal energy of $u_{g,\rm fl}=\max[b^{2}/750,2\times 10^{-5}(r/r_{g})^{-2.5\Gamma}]$, where $b=\sqrt[]{b^{\mu}b_{\mu}}$, $b^{\mu}$, $\sigma_{0}$ and $\Gamma$ are the co-moving magnetic field strength, magnetic four-vector, maximum magnetisation and the ideal gas law adiabatic index, respectively.

3 Simulation models

We have carried out 9 simulations (including 2 models from the Appendix) to elucidate the physics of jet acceleration and interaction with the ambient medium (see Table 1). Most of our models ran for $t>10^{5}t_{g}$. In all models we start with a Fishbone & Moncrief (1976) (FM76) torus in hydrostatic equilibrium around a rapidly spinning $a=0.9375$ Kerr black hole. We place the torus inner edge at $r_{\rm in}$ and density maximum at $r_{\rm max}$. For all but one of our simulations, we set $r_{\rm in}=36r_{g}$ and $r_{\rm max}=73.97r_{g}$. We choose the density scale by setting $\max\rho=1$. We adopt the ideal gas law equation of state and set the adiabatic index to that of a non-relativistic gas $\Gamma=5/3$. We set the magnetic field vector potential as described in Secs. 3.1 and 3.2 and normalise it such that $\max p_{\rm g}/\max p_{\rm B}=100$, where $p_{\rm g}$ and $p_{\rm B}$ are the gas and magnetic pressure, respectively.

3.1 Single field loop setup

We first start with a single poloidal magnetic field loop, B10-S model, where ‘S’ stands for ‘single loop’, as seen in the left panel of Fig 1. We choose the magnetic vector potential of the form:

[TABLE]

As the black hole starts accreting the gas, it drags along the magnetic flux, which starts threading the event horizon. Over time, the accumulated magnetic flux becomes large enough to push the gas away and prevent it from accreting, as seen in Fig. 2. This is an axisymmetric variant (Proga & Zhang, 2006) of the magnetically arrested disc (MAD) state (Igumenshchev et al., 2003; Narayan et al., 2003; Tchekhovskoy et al., 2011), which leads to highly efficient outflows: their efficiency, defined as total outflow power normalised by the time-average accretion rate, can exceed unity: $P_{\rm outflow}/\langle\dot{M}c^{2}\rangle>1$. Here, $P_{\rm outflow}=\dot{M}c^{2}-\dot{E}$, where $\dot{E}=\iint F_{\rm E}dA_{\theta\phi}$ is the energy accretion rate (defined to be positive when the energy flows into the black hole, see Eq. 3 for $F_{\rm E}$ definition), $dA_{\theta\phi}=\sqrt{-g}d\theta d\phi$ is the surface area element, and $g=|g_{\mu\nu}|$ is the determinant of the metric. Similarly, we can define the magnetic flux on the black hole as $\Phi_{\rm BH}=0.5\iint|B^{r}|dA_{\theta\phi}$, where the integral is over the entire event horizon. We also define its dimensionless counterpart normalised by the mass accretion rate, $\phi_{\rm BH}=\Phi_{\rm BH}/\langle\dot{M}r_{g}^{2}c\rangle$.

In axisymmetry, the magnetic fields on the black hole and the surrounding gas have no way of exchanging places: the magnetic interchange instability requires a third dimension. Thus constrained by 2D symmetry, the fight between gravitational and magnetic forces degenerates into the gas bouncing in and out on top of the magnetic barrier, as seen in Fig. 2. Figure 3 shows that this results in large (up to an order of magnitude) fluctuations in the mass accretion rate $\dot{M}$ (Fig. 3a), dimensionless magnetic flux $\phi_{\rm BH}$ (Fig. 3b) and the dimensionless total outflow power (Fig. 3c). Such oscillations make it difficult to extract the physics from the simulations and thereby make this configuration undesirable.

3.2 Fiducial setup

We initialise our fiducial model B10 with a disc threaded with a large enough magnetic flux such that a powerful jet forms while taking care to avoid over-saturating the black hole and getting a MAD (Fig. 1, right). The magnetic field configuration in the torus consists of two poloidal field loops described by the following vector potential,

[TABLE]

where $f=0.1x_{1}(\rho-0.05)^{1/2}+0.9x_{1}^{2}(\rho-0.05)^{2}\sin^{2}[\pi(x_{1}-2)/2]\sin(\pi x_{2}/2)$. The first term in $f$ describes a large-scale field loop and the second term a pair of oppositely polarised smaller loops embedded within the large loop (in terms of the internal coordinates: $x_{1}$ and $x_{2}$; Appendix B). The magnetic fluxes within the pair cancel exactly such that the total magnetic flux is set by the large-scale loop. This cancellation is convenient, because it gives us fine-grained control over the amount of net magnetic flux in the initial conditions and allows us to choose the positive polarity to dominate only slightly over the negative one. Unlike the model B10-S, which showed violent variability in the black hole mass accretion rate $\dot{M}$, dimensionless magnetic flux $\phi_{\rm BH}$ and the total outflow power, model B10 shows a steady behaviour of all three quantities (Fig. 3), primarily because it has a smaller magnetic flux in the disc. This makes model B10 ideal for studying long duration steady outflows.

4 Fiducial model results

4.1 Global evolution

In this section, we focus on our fiducial model B10. With time, in model B10, the accretion disc develops turbulence through the magneto-rotational instability (MRI; Balbus & Hawley, 1991), which leads to accretion onto the black hole and launching of the jets on both sides of the disc (Figure 4). For accretion to take place, angular momentum needs to be redistributed to the outer parts of the disc via the MRI and therefore, it is important for simulations to properly resolve the MRI turbulence. To quantify this, we calculate the quality factors $Q^{r,\theta}$, where $Q^{i}=\langle 2\pi v_{A}^{i}/(\Delta^{i}\Omega)\rangle_{w}$ measures the number of cells per MRI wavelength in direction $i=$[$r,\theta$], where $v_{A}^{i}$ is the Alfvén velocity, $\Delta^{i}$ the cell size, $\Omega$ the angular velocity of the fluid. $Q^{i}$ is averaged over the inner disc ($r<50r_{g}$) and weighted by $w=\sqrt[]{b^{2}\rho}$. We achieve $Q^{\theta}\sim 14$ (Table 1) at $t=2\times 10^{4}t_{g}$, fulfilling the numerical convergence criteria (see e.g., Hawley et al., 2011). $Q$ values decrease over time as expected for axisymmetric systems (Cowling, 1934), but since we focus on the physics of the jet, it is not a significant concern that we do not resolve the MRI well in the disc at late times.

Initially, the jet expands as if there were no confinement (i.e., ballistically) until its ram pressure drops below the confining pressure of the disc-wind, which snaps back on the jet. Figure 4 shows that this unstable interaction between the disc-wind and the jet leads to oscillations of the jet-wind interface: that we refer to as pinches. The pinches also give rise to small scale eddies that mass-load the jet at late times222Movie showing magnetised jet formation of model B10: https://youtu.be/4MeLZZPYsfc (see Sec. 6.4 for a more detailed discussion). Interestingly, even at an early time, the jet on one side of the disc shows qualitative differences in behaviour compared to the other jet. Namely, the upper jet ($z>0$ in Fig. 4) is strongly affected by the interface instabilities, while the lower jet ($z<0$ in Fig. 4) remains much more stable.

To understand how the oscillating interface affects jet dynamics, we look at how the jet evolves along a field line, which we define as a surface of constant enclosed poloidal magnetic flux, $\Phi(r,\theta)=\int^{\theta}_{0}B_{p}dA_{\theta\phi}$, where $B_{p}$ is the poloidal field strength. We choose a field line whose foot-point makes an angle of $\theta_{j,\rm H}=0.44$ rad with the black hole spin axis at the event horizon. This makes up about $40$% of the jet’s half-opening angle $\theta_{\rm jet}$. Figure 5 shows the evolution of various quantities along our chosen field line at $t=2\times 10^{5}t_{g}$, a long enough time for the jet to establish a quasi-steady state solution out to $10^{5}r_{g}$. First, we consider the total specific energy $\mu$, which is the maximum Lorentz factor a jet could attain if it converted all forms of energy into kinetic energy, and equals the ratio of the total energy flux, $F_{\rm E}$, and the rest-mass flux, $F_{\rm M}$,

[TABLE]

where $\rho$ is gas density, $u_{\rm g}$ is internal energy density, $p_{\rm g}=(\Gamma-1)u_{\rm g}$ is gas pressure, $b^{\mu}$ and $u^{\mu}$ are the magnetic and velocity four-vectors respectively. Figure 5 shows that the radial profile of $\mu$ stays approximately constant for both jets, with small oscillations due to the lateral movement of gas in response to the jet pushing against the confining pressure of the disc-wind (Lyubarsky, 2009; Komissarov et al., 2015). The flip in the sign of $\mu$ around $8r_{g}$ is caused by the presence of a stagnation surface (see Sec. 2.2), where $u^{r}=0$ and thus $\mu\rightarrow{\infty}$. Downstream of the stagnation surface, the floor values (see Sec. 2.2) set the value of $\mu$. To quantify the conversion efficiency of magnetic to kinetic energy, we calculate the ratio of the Poynting flux, $F_{\rm EM}$, to the mass energy flux, $F_{\rm K}$, called magnetisation $\sigma$:

[TABLE]

where $u_{t}\sim-\gamma$ for $r\gg r_{g}$. Figure 5 shows that $\sigma$ decreases to below unity as magnetic energy is converted into kinetic energy. This means that even though the jets started out strongly magnetically-dominated near the black hole, the process of acceleration converted their magnetic energy into kinetic form to the point where the jets end up becoming kinetically-dominated. Whether $\sigma$ will keep decreasing even further, leading to unmagnetised jets, or will level off around unity, leading to somewhat magnetised jets, will require a simulation extending to even larger distances. In the following sections, we take a closer look at the jet acceleration profile as well as the dissipation due to interface instabilities.

4.1.1 Jet structure and acceleration

Figure 6(a) shows that both the upper and the lower jets collimate similarly from $\theta_{j}\sim 0.3$ rad (or $17.2\deg$) at $8r_{g}$ to $\theta_{j}\sim 8\times 10^{-4}$ rad (or $0.046\deg$) at $10^{5}r_{g}$ while displaying a power-law shape $\theta_{j}\propto(r/r_{g})^{-0.63}$. Remarkably, the outer jet also displays a continuous power-law collimation profile and resembles that of the M87 jet (see Sec. 6.1).

Beyond the stagnation surface, the Lorentz factor (Fig. 5) increases smoothly until $r\lesssim 200r_{g}$ for both the lower and upper jets. Such a $\gamma$ profile is typical for highly magnetised jets in the poloidal field dominated regime (e.g., Beskin et al., 1998; Tchekhovskoy et al., 2008) where the acceleration occurs on the scale of the light cylindrical radius, $R_{\rm L}=c/\Omega$, where $\Omega$ is the conserved field line angular velocity. The Lorentz factor in this regime behaves as $\gamma\approx(1+(R/R_{\rm L})^{2})^{1/2}$, where $R$ is the cylindrical radius of the field line.

Once the jet becomes super-fast magnetosonic (i.e., the jet velocity becomes larger than the fast magnetosonic wave speed), adjacent field lines must shift in the transverse direction in a non-uniform manner for the jet to efficiently accelerate (Begelman & Li, 1994; Chiueh et al., 1998; Vlahakis, 2004; Tchekhovskoy et al., 2009; Komissarov et al., 2009). This can be thought of as differential bunching of field lines towards the jet axis. To quantify this bunching, we can use Eq. (26) of Tchekhovskoy et al. (2009),

[TABLE]

where we define the field line bunching parameter $a_{bp}$ as the ratio of the local poloidal field strength $B_{p}$ and the mean poloidal field strength, $\Phi/\pi R^{2}$. For efficient acceleration, the poloidal field must decrease faster with distance than the mean poloidal field, hence creating a pressure gradient that exceeds the hoop stress, which slows down the jet, thereby accelerating the jet. Indeed, Fig. 6(b) shows that the bunching parameter $a_{bp}$ does decrease and therefore, $\gamma$ increases to values approaching, within a factor of few, the maximum Lorentz factor $\mu$. The Lorentz factor in this simulation reaches $\gamma\sim 6$, which is consistent with typical AGN jets (e.g., Pushkarev et al., 2017).

4.1.2 Transverse jet causality

As we discussed in Sec. 4.1.1, for efficient jet acceleration, magnetic field lines need to move across the jet and bunch up towards the axis. This requires the jet to be laterally causally connected (Tchekhovskoy et al., 2009; Komissarov et al., 2009). We approximate that a field line is in lateral causal connection with the jet axis as long as its Mach cone (i.e., the range of directions that a point on the field line can communicate with) crosses the axis of the jet. Thus, for efficient communication with the axis, the jet half-opening angle $\theta_{j}$ must be smaller than the Mach cone half-opening angle ($\theta_{\rm M}=1/M_{f}$, where $M_{f}$ is the fast Mach number). Using the definition of the fast magnetosonic wave velocity $v_{f}$ and Lorentz factor $\gamma_{f}$ (e.g., Gammie et al., 2003),

[TABLE]

and focusing on the asymptotic regime of the jet, i.e., $R\gg R_{\rm L}$ and $b_{t}\approx 0$, we can compute the Mach cone half-opening angle as,

[TABLE]

Approximating the flow velocity $v\approx 1$ (for a relativistic jet), we have $\theta_{\rm M}\approx\sqrt{\sigma}/\gamma$. Figure 6(c) shows the transverse causality parameter $\theta_{j}/\theta_{\rm M}$, the ratio of the jet and Mach cone half-opening angles. For $\theta_{j}/\theta_{\rm M}>1$, we expect acceleration to slow down as causal contact is lost (Tomimatsu, 1994; Beskin et al., 1998; Tchekhovskoy et al., 2009). However Fig. 6(c) demonstrates that the flow along the field line always remains in causal contact with the polar axis, which is consistent with VLBI observations of AGN jets (e.g., Jorstad et al., 2005; Clausen-Brown et al., 2013) that show $\gamma\theta_{j}\sim 0.1-0.3$ (Fig. 6d). The acceleration slow down occurs as the jet becomes matter-dominated, i.e., $\sigma<1$, with additional deceleration due to pinching which we discuss next.

4.1.3 Toroidal pinch instabilities

The upper jet has a distinctly different acceleration profile compared to the lower jet. As the jet propagates through the ambient medium, it expands and adiabatically cools down. Due to this, the specific enthalpy $h=(u_{\rm g}+p_{\rm g})/\rho$ (Fig. 5, black line) decreases initially, slightly at $r\lesssim 100r_{g}$. However, beyond $100r_{g}$, pinch instabilities cause magnetic dissipation that raises the jet specific enthalpy to order unity, effectively creating a thermal pressure gradient directed against the direction of the flow at $r\lesssim 10^{3}r_{g}$. This pressure gradient significantly slows down the outflow to $\gamma\lesssim 2$ by $r\approx 10^{3}r_{g}$. At larger radii, $h$ drops, and the jet re-accelerates under the action of both the magnetic and thermal pressure forces. We can interpret this behaviour also through energy conservation. Using the definition of magnetisation (Eq. 4), specific enthalpy and approximating $u_{t}\sim-\gamma$, from Eq. (3) we get

[TABLE]

a useful form of the energy equation. It clearly shows that for $\sigma\simeq$constant, an increase in enthalpy to $h\sim 1$ results in $\gamma$ decreasing, which is seen for the upper jet around $10^{3}r_{g}$. The upper jet shows stronger pinch activity and thus collimates slightly more than the lower jet: $\theta_{j}$ is smaller by a factor of $\lesssim 2$; Fig. 6a, solid line), and other quantities like the maximum Lorentz factor $\mu$ and bunching parameter $a_{bp}$ strongly oscillate.

4.1.4 Energetics across the jet

Figure 7 shows how the different components of jet energy flux, $\mu$, $\sigma$, $\gamma$ and $h$, vary across the jet at different distances along the upper jet. A fiducial field line with the foot-point at $\theta_{j,\rm H}=0.44$ rad (see also Fig. 5), shown with filled circles, collimates faster than the jet-edge, indicated by $\theta_{\rm jet}$. We take the jet-edge to be at $\mu=1.2$, which is reasonable since $\mu$ should drop rapidly at the jet-edge 333We also note that $\mu$ is roughly constant throughout the width of the jet as opposed to increasing as $\mu\propto\sin^{2}\theta$, seen in previous simulations (e.g., Komissarov et al., 2007; Tchekhovskoy et al., 2008), which might be a consequence of the density floor model coupled to the stagnation surface.. This means that there is internal reconfiguration of the flow within the body of the jet that leads to the formation of a fast magnetised inner core. At the jet-edge, mass-loading via pinching (Sec. 6.4) causes specific total energy $\mu$, the magnetisation $\sigma$ and the Lorentz factor $\gamma$ to drop gradually, forming a slower sheath that surrounds the core, resulting in a structure similar to the spine-sheath seen in AGN jets.

4.2 Dependence on mass-loading at the stagnation surface

In Sec 4.1.4, we showed that pinching instabilities lead to dissipation and therefore, reduces jet acceleration efficiency. Would the acceleration efficiency change if the jet had a different specific total energy? In order to answer this question, we compare results from four models different only by the density floors (Sec. 2.2; also see Table 1), namely the maximum magnetisation $\sigma_{0}$ values of 3 (model B3), 50 (model B50) and 100 (model B100), along with model B10. Figure 8 shows that in all models the Lorentz factor $\gamma$ increases steeply until $\sim 10^{3}r_{g}$, beyond which the acceleration slows down. Models with higher jet base magnetisation (and hence, larger $\mu$ values) accelerate slightly faster and reach higher Lorentz factors, but get affected by pinch instabilities at roughly the same distance as models with lower $\mu$. For model B3 the pinch instability is weaker and leads to smaller oscillations in, for example, the Lorentz factor (see discussion in Sec. 6.3). Similar to Komissarov et al. (2009), we find that models with smaller $\mu$ achieve smaller $\sigma$ at large distances. For instance, model B3 achieves $\sigma\simeq 0.2$ at $r=10^{5}r_{g}$. It is interesting to note that efficient heating via relativistic shocks requires $\sigma\lesssim 0.1$ (e.g., Kennel & Coroniti, 1984; Komissarov, 2012). Such low $\sigma$ values are very difficult to achieve for collimated flows as $\sigma$ drops very gradually in the toroidally dominated regime (known as the $\sigma$ problem; e.g., Tchekhovskoy et al., 2009; Komissarov et al., 2009). This suggests that magnetic reconnection might be a more efficient mechanism for particle acceleration in jets.

4.3 Acceleration of a jet collimated by a small disc

So far, we discussed models with large discs that collimate the jets out to large distances. Here, we consider model B10-R with a small disc. In this case, the disc wind collimates the jets out to smaller distances444Movie showing the difference in collimation and acceleration between B10 and B10-R: https://youtu.be/2C4re4aiuQM. This happens because larger, more radially extended discs launch disc-winds over an extended range of radii: the winds launched from small radii collimate off of those launched further out, and off of the disc itself, leading to a radially extended collimation profile of the jets. B10-R, with a smaller disc extending only up to $500r_{g}$, is embedded with the same magnetic field configuration as model B10 (Eq. 2). The left panel in Fig. 9 shows the vertical slice through the initial conditions and the right panel shows the system at $t=2\times 10^{5}t_{g}$. Model B10-R, compared to B10 (Fig. 4, right), has a much wider jet as the weaker confining pressure of the disc-wind enables the jet to expand laterally.

In Fig. 10 we show the radial profiles of quantities along two magnetic field lines in the upper jet, one located in the inner jet (foot-point half-opening angle of $\theta_{j,\rm H}=0.41$ rad; highlighted in pink in Fig. 9, right panel) and the other in the outer jet ($\theta_{j,\rm H}=1.17$ rad; magenta in Fig. 9, right panel), at $t=2\times 10^{5}t_{g}$. Initially, similar to our fiducial model B10, the disc and disc-wind collimate the jet into a parabolic shape. However, beyond $10^{3}r_{g}$, the confining pressure of the disc drops and the outer field lines in the jet become conical (Fig. 10b). The deconfinement leads to a drop in the confining pressure, causing field lines to diverge and experience a quicker acceleration due to the outwards pressure gradient. This boost in $\gamma$ due to smooth deconfinement of the jet has been shown to occur by previous idealised simulations (e.g., Tchekhovskoy et al., 2010a; Komissarov et al., 2010), though the increase in acceleration is not quite as significant as in Fig. 2 of Tchekhovskoy et al. (2010a). Beyond $10^{3}r_{g}$, the quick lateral expansion of the jet suppresses pinch instabilities to a large extent (consistent with e.g., Moll et al., 2008; Granot et al., 2011; Porth & Komissarov, 2015) and adiabatically cools the jet leading to an order of magnitude smaller enthalpy $h$ compared our fiducial model B10 (Fig. 5).

Even though the outer field line expands ballistically and maintains an approximately constant opening angle, the inner field lines continue to collimate off of the outer ones into a parabolic shape similar to our fiducial model (Fig. 5). The outer jet experiences a relatively larger change in the bunching parameter $a_{bp}$ compared to the inner field line (Fig. 10c), in accordance with Eq. (5) and reaches $\sigma\simeq 0.2-0.3$ (Fig. 10a). Upon loss of collimation, the outer jet also loses transverse causal connection ($\theta_{j}/\theta_{\rm M}>1$; Fig. 10d), and the acceleration ceases (see Sec. 4.1.2). From Fig. 10(e), $\gamma\theta_{j}\sim 0.7-1$ for the outer jet, while $\gamma\theta_{j}$ is between 0.1 and 0.4 within $1000r_{g}$ for the inner jet. See Sec. 6.1 for further discussion of $\gamma\theta_{j}$ values in our models.

5 Comparisons to idealised jet simulations

Here we aim to study jet dynamics in the absence of pinching instabilities by constructing smooth idealised outflows and maintaining fine control over the jet shape by confining the flow using a conducting collimating wall (Komissarov et al., 2007, 2009; Tchekhovskoy et al., 2010b). Such a setup also removes the shear-induced turbulence and dissipation at the jet edge-disc wind layer.

5.1 Model setup

We set up an outflow bound by a perfectly conducting wall, mimicking a jet collimated by an external medium. We refer to this setup as a wall-jet simulation, in contrast to disc-jet simulations in which the disc-wind collimates the jets. The field lines threading the event horizon initially follow the shape of the wall that collimates in a parabolic fashion:

[TABLE]

This gives us the initial poloidal field configuration:

[TABLE]

Here, the outermost field line touching the wall is given by $A_{\phi}=1$: it starts out at the intersection of the event horizon ($r=r_{\rm H}$) and the equatorial plane ($\theta=\pi/2$), is initially radial for $r\lesssim r_{0}$ and asymptotically collimates as $\theta\approx r^{-\nu/2}$. In this setup, $\nu=0$ results in a monopolar field shape ($\theta_{j}=$constant), while $\nu=1$ gives us the parabolic field shape ($\theta_{j}\propto r^{-1/2}$). We set a transitional radius $r_{0}=10r_{g}$ and employ $\nu=0.8$ as these values give a good match to the disc-jet shape, as we discuss below. We define the physical coordinates ($r$, $\theta$) as functions of the internal coordinates ($x_{1}$, $x_{2}$) as $r=\exp(x_{1})$ and $x_{2}=\mathrm{sign}(\theta)|A_{\phi}|^{1/2}$. We use a resolution of $12800\times 400$ cells. Our computational domain range extends radially from $0.85r_{\rm H}$ to $10^{6}r_{g}$. We employ the same polar reflective boundary conditions at $x_{2}=0$ as for the disc-jet simulations (see Sec. 2.1). At the wall, $x_{2}=1$, the boundary conditions are also reflective so that the gas and fields follow the wall (Tchekhovskoy et al., 2010b). The density floors are the same as in model B10. The left panel in Fig. 11 shows the resulting wall-jet solution. For comparisons with the disc-jet model B10, we choose a field line in the inner jet of B10 model as field lines near the jet-edge are strongly affected by mass-loading (Fig. 11, right). At time $t=0$, we start the wall-jet simulation with a purely poloidal magnetic field given by Eq. (10), which then develops a toroidal component due to the rotation of the black hole. Our simulation time extends to $10^{7}t_{g}$, which ensures that the outflow reaches a steady state up to at least a distance of $10^{6}r_{g}$. To speed up the simulation, from $t=1000t_{g}$ onwards, we freeze out cells that reached steady state, i.e., the cells located at $r<0.1ct$, where $t$ is the simulation time (similar to Tchekhovskoy et al. 2008; see also Komissarov et al. 2007).

5.2 Disc-jets vs. idealised wall-jets

Figure 11 shows the comparison between the ideal wall setup (left) and disc-jet setup (right). For the disc-jet setup, the presence of a pressure imbalance between the jet and the accretion disc-wind gives rise to oscillations in the jet shape. In contrast, for the idealised wall-jet, the boundary is rigid and hence, the pinches are absent. The energy flux components, $\mu$, $\sigma$, and $\gamma$ agree between the disc-jet and the wall-jet rather well (Fig. 12a), showing that wall-jet models capture most of the time-average steady state dynamical properties of disc-jet models with the same shape (Fig. 12b), especially in the absence of pinches. For the wall-jet, the specific enthalpy $h$ decreases with increasing $r$ as expected due to adiabatic jet expansion. For the disc-jet, $h$ increases substantially at $r\sim 200r_{g}$ due to the onset of the pinch instabilities that convert the poloidal field energy into enthalpy. Free of pinches, the wall-jet smoothly accelerates as $\gamma\propto R$ until a few times $10^{3}r_{g}$, followed by the slower acceleration as the field lines slowly become cylindrical when they enter the jet core (Fig. 12a). The acceleration is more rapid for field lines closest to the wall as field lines in this region diverge away from each other more (Fig. 11, left). For the disc-jet, the presence of the pinches causes $B_{p}$ to dissipate (Fig. 12c), along with a slight drop in $\theta_{j}/\theta_{\rm M}$ (Fig. 12d). The product of the Lorentz factor and jet opening angle $\gamma\theta_{j}\gtrsim 0.1$ for $r<3000r_{g}$(Fig. 12e), similar to the values found for the inner jet (see Sec. 4.1.2).

6 Discussion

6.1 Comparison to the M87 jet

We find that all of our simulated jets with large discs propagate with a parabolic shape over at least 5 orders of magnitude in distance. In this section, we consider whether our fiducial model jet behaves the same way as the jets seen in nature, looking at the shape and acceleration profiles inferred from multiple Very Large Array (VLA)/ Very Long Baseline Interferometry (VLBI) observations of the M87 jet. Figure 13(a) compares the jet geometry for our fiducial model B10 with the observed M87 jet. The shape of the field line near the jet-edge fits very well with the observed data, displaying a parabolic collimation profile up to $10^{5}r_{g}$, close to the location of HST-1 (Asada & Nakamura, 2012; Hada et al., 2013; Mertens et al., 2016; Kim et al., 2018b; Nakamura et al., 2018).

Mertens et al. (2016) showed using VLBI measurements that the acceleration profile of M87 follows $\gamma\propto R\propto z^{0.58}$ till $10^{3}r_{g}$, changing to $\gamma\propto z^{0.16}$ up to the HST-1 knot. These power-law profiles agree reasonably well with the Lorentz factor profile along a field line in the mid-jet as can be seen in Fig. 13(b). The discrepancy at small radii may result from systematic measurement errors in the VLBI observations, the jet opening angle becoming comparable to the viewing angle or our preference for a particular field line. Additionally, in Fig. 13(a), the D12 (Doeleman et al., 2012) and A15 (Akiyama et al., 2015) data points representing the Event Horizon Telescope Core at 230 GHz depend considerably on the assumed black hole mass and viewing angle along with additional uncertainties on the position (Nakamura et al., 2018). Hence, it is possible that instead of the 230GHz core being smaller than the jet interior as the figure suggests, the emission might come from the disc, i.e., outside the jet-edge, where the Lorentz factor $\sim 1$. The outer field line Lorentz factor is close to 1.1 (Fig. 7), similar to the velocities found for the M87 jet sheath (Mertens et al., 2016) and agrees with previous GRMHD simulations modified for M87 (Nakamura et al., 2018). Indeed, as the jet gets mass-loaded via the jet-wind interaction, we expect a gradual decrease in the Lorentz factor as we go from the inner jet to the jet-edge, which may explain the wide distribution of the Lorentz factors across the M87 jet in Fig. 15 of Mertens et al. (2016).

The HST-1 knot in M87 is a region where the jet is deemed to over-collimate and transitions from parabolic to conical structure (Asada & Nakamura, 2012). Unfortunately, we do not find such a dissipative feature in any of our simulations, nor do we see the jet turn conical around $10^{5}r_{g}$. One possible reason may be that HST-1 lies very close to the Bondi radius of M87 ($\sim 7.6\times 10^{5}r_{g}$; Nakamura & Asada, 2013) where the shallow density profile of the ISM prevails. If there is an increase of confining pressure from the ISM beyond $10^{5}r_{g}$, it is possible that the jet becomes over-pressured, perhaps forming a re-collimation feature. Further, as is the case for model B10-R (see Sec. 4.3), if the jet pressure subsequently becomes larger than the confining pressure, the jet would open up and turn conical. This suggests that it is important to consider a more realistic ISM pressure profile in future work (e.g., Barniol Duran et al., 2017).

The product of the Lorentz factor and the jet opening angle $\gamma\theta_{j}$ is an important quantity we use for comparison to AGN jets. It is clear from Fig. 6 that the inner jet exhibits very low values of $\gamma\theta_{j}$ ($<0.01$) for distances larger than $10^{3}r_{g}$, compared to those observed ($\sim 0.1-0.3$, e.g., Jorstad et al. 2005; Pushkarev et al. 2009; Clausen-Brown et al. 2013; Jorstad et al. 2017), whereas for the jet-edge, we find $\gamma\theta_{j}<0.1$. It is possible that the difference in measurements of $\gamma\theta_{j}$ between our models and observed jets might be a result of the latter assuming a conical jet as well as the uncertainty of attributing the Lorentz factor of the underlying jet flow to emission features (e.g., there might be standing shocks in the jet). In the case of model B10-R, the outer jet becomes conical and causally disconnected: $\gamma\theta_{j}\gtrsim 1$ (see Fig. 10), which may be more applicable for jets in gamma-ray bursts. Additionally, the peak Lorentz factor for many of the observed jets is over 10, a value only one of our models achieves (model B100), suggesting highly magnetised jets.

6.2 Causal structure of jets

When a jet becomes super-fast magnetosonic (i.e., downstream of the fast magnetosonic surface), perturbations in the jet cannot be communicated upstream and thus, the jet loses causal connection along its flow. However, beyond the fast surface, globally, causal contact can still be maintained as the jet can communicate upstream via the sub-fast jet axis. Thus, when the flow along the axis turns super-fast, the jet reaches a magnetosonic horizon and causal connection is fully lost. The location where this causal breakdown occurs is the fast magnetosonic separatrix surface (FMSS; for a review, see Meier, 2012). Self-similar models (such as e.g., Vlahakis, 2004; Polko et al., 2010; Ceccobello et al., 2018) predict that a jet collapses on its axis once the jet reaches the FMSS and may form a highly radiating hot-spot. Could then bright features in the jets, such as HST-1, be powered by such over-collimation seen in self-similar models?

To test if the FMSS can explain bright jet features, we have developed an algorithm that determines the FMSS location. This algorithm calculates the Mach cone angle (Eq. D5 in Tchekhovskoy et al., 2009) for each cell assuming approximate magnetosonic fast wave velocity (Gammie et al., 2003). We track the left and right edges of the Mach cone to check if a fast magnetosonic wave can travel to the sub-fast region near the the jet’s axis. Figure 14(a), salmon line) shows that the FMSS in model B10 travels inwards across the jet from the outer boundary, before joining with the fast surface at the jet’s axis. The FMSS does not coincide with dissipative features, which are shown in Fig. 14(a), via the entropy $s$,

[TABLE]

a proxy for identifying shocks and magnetic dissipation (Barniol Duran et al., 2017). Instead, the fast surface (Fig. 14(a), blue) coincides with the steady rise in entropy at the jet-edge (Fig. 14b). Outside of the jet, in the disc-wind, Fig. 14(b) shows that entropy begins increasing in the sub-fast regime (due to shearing between the disc and the wind) and continues to rise smoothly till the jet’s fast surface. The above results suggest that the fast surface in the jet plays a role in triggering events which cause dissipation.

While we can compare the jet structure with radially self-similar models (e.g., Blandford & Payne, 1982; Vlahakis, 2004; Polko et al., 2010; Ceccobello et al., 2018), the lack of an over-collimation in our simulations suggests that there is a difference in the way the FMSS manifests in the jet within a self-similar approximation. Namely, in radially self-similar models, the FMSS is located where the flow achieves super-magnetosonic speeds towards the polar axis, which our models never reach. That radially self-similar models restrict the radial dependence of quantities to fixed power-laws, which is not the case in our simulations, might be the crucial difference that leads to different nature of FMSS in the self-similar models. Asymptotically in our jets, field lines join with the jet “core”, by which point they become almost cylindrical and stop accelerating efficiently. Perhaps this asymptotic behaviour can be explored in self-similar models by placing the FMSS at infinity (e.g., Li et al., 1992; Vlahakis & Königl, 2003b) or via $\theta-$self-similar models (e.g., Sauty et al., 2004).

6.3 Origin of pinch instabilities

Pinch instabilities forming at the jet-wind interface are easily excited in 2D GRMHD simulations of black hole accretion (see e.g. McKinney 2006; Barkov & Baushev 2011; Nakamura et al. 2018). Whereas previous work found that pinch instabilities do not survive beyond $\sim 10^{3}r_{g}$ (e.g., McKinney, 2006), we observe them significantly affecting jet dynamics throughout the length of the jet. We hypothesise that this difference in pinch activity stems from the small disc size in McKinney (2006), as a smaller disc would lead to a conical jet in which pinching instabilities are suppressed (we see the same behaviour for model B10-R, see Sec 4.3). However, the amount of dissipation seen in McKinney (2006) is far larger ($>2$ orders of magnitude) than in any of our models. To address this discrepancy, we ran an additional low-resolution simulation B10-SLR. Because the B10-SLR model under-resolves the pinches at large radii, this leads to enhanced dissipation at $r\gtrsim 10^{3}r_{g}$ (Fig. 15). Additionally, the floor model used in McKinney (2006) could also contribute to larger dissipation.

To better quantify the location at which the pinches start to affect the jet, we show in Fig. 16(a) the standard deviation in the jet opening angle along a field line for several models. From Fig. 16(b), we see that the field lines tend to wobble significantly due to pinching, achieving $\sqrt{\langle(\Delta\theta_{j})^{2}\rangle}/\langle\theta_{j}\rangle>10\%$ close to the fast surface. The fast surface could play a role here since beyond the fast surface, the ram pressure of the wind may become the dominant pressure component and may produce shock-like events at the jet-disc wind interface (e.g., Komissarov, 1994; Bromberg & Levinson, 2007). However, in our simulations, we find that the wind remains subsonic in the $\theta-$ direction, suggesting ram pressure in the $\theta-$ direction is not prominent. Indeed, we see a smooth increase of entropy across the jet-disc interface, indicating no prominent shocks (Fig. 14). In Fig. 16(c), we see that the inner jet crosses the fast surface at a very large distance, while the entropy begins to rise much earlier. However, the fast surface coincides with the increase of entropy for the mid- and outer-jet. These results suggest that the oscillations in the jet-wind interface (that form very close to the black hole) might give rise to the pinch instabilities, which grow significantly once the jet-edge becomes super-fast.

However at this point, it is not clear whether the oscillating interface at small radii and pinches at large radii are due to the same underlying physical phenomenon. Both the oscillating interface and the pinches appear to be the response of the jet to the pressure of the surrounding disc-wind (Sobacchi & Lyubarsky, 2018a). Indeed, both the oscillations and pinches disappear in the case of the idealised wall-jet, where the rigid wall prevents a dynamic jet-wind boundary. The jet becomes susceptible to pinch instabilities when the toroidal field dominates over the poloidal field (given by the Tayler criterion: Tayler 1957; also see Sobacchi & Lyubarsky 2018b). We note that the jet-edge field line shown in Fig. 16(c) becomes strongly toroidal at a very small distance. The growth rate of the pinch/kink instability scales with the $\phi-$ component of the Alfvén velocity, which is proportional to the toroidal field strength (e.g., Moll et al., 2008). As Moll et al. (2008) also notes, jet expansion restricts the growth of the pinch, seen in the case for model B10-R, where the small disc allows rapid jet de-collimation, and hence, the jet exhibits weak pinching.

Interestingly, pinches begin to noticeably heat up the jet within $200-800r_{g}$ from the central black hole (agreeing with e.g., Giannios & Spruit, 2006). This is similar to the distances at which the synchrotron break is estimated to occur for both AGN and X-ray binary jets (e.g., Markoff et al., 2001, 2005; Russell et al., 2013; Lucchini et al., 2018). The break arises when the synchrotron emission of a compact jet shifts from its characteristic power-law profile to a flat/inverted spectrum due to self-absorption, transitioning from an optically thin regime at higher frequencies to optically thick (Blandford & Königl 1979; for a review, see Markoff 2010; Romero et al. 2017) and is generally attributed to non-thermal emission from particle acceleration caused by e.g., shocks (e.g., Sironi & Spitkovsky, 2009) or magnetic reconnection (e.g., Spruit et al., 2001; Drenkhahn & Spruit, 2002; Sironi & Spitkovsky, 2014; Sironi et al., 2015) in the jet. Given that pinching sets off magnetic reconnection in our simulations, we suggest that the start of the pinch region may potentially be an ideal site for particle acceleration to occur for the first time and hence, can manifest itself as a break in the synchrotron spectrum. Additionally, pinching may cause variation in the optical depth for the synchrotron self-absorption, leading to variability in the observed jet depth at a given frequency over time, which might have consequences for the radio core shift in AGN jets (Blandford & Königl 1979; Plavin et al. 2019; for M87: Hada et al. 2011). Using GRMHD simulations that extended out to $100r_{g}$, Nakamura et al. (2018) found higher values of Lorentz factor in pinched regions and suggested that if such compressions lead to dissipation, pinches can be associated with superluminal blobs observed in the jets.

As the jet base magnetisation increases (see Sec. 4.2), Fig. 16(a) shows that that the fast surface moves out, away from the launchpoint of the jet (consistent with results from radially self-similar models: Fig. 11, left panel of Ceccobello et al. 2018). Starting from Eq. (6), we can derive the approximate distance at which the fast surface resides. We have,

[TABLE]

assuming a cold jet, i.e., specific enthalpy $h\ll 1$. At the fast surface ($\gamma=\gamma_{f}$), we then have for the jet,

[TABLE]

Assuming that $\gamma\approx\Omega R$ and $\theta_{j}=C(r/r_{g})^{-\zeta}$, we arrive at the location of the fast surface,

[TABLE]

From Eq. (14), we can indeed say that $r_{f}$ increases with increase in $\mu$ for a given field line. For $\mu=10$, Eq. (13) gives $\gamma_{f}=2.154$ which is about $25\%$ off from the simulation value ($\approx 2.7$). However, the assumption of $\gamma\approx\Omega R$ is not valid in the outer jet, due to stronger effects of mass-loading. The time variability of $\theta_{j}$ decreases for model B10 through to model B100, which suggests that larger magnetisation stabilises against pinching activity (consistent with results of e.g., Mizuno et al., 2015; Fromm et al., 2017; Kim et al., 2018a). With higher magnetisation, the Alfvén speed and subsequently the magnetosonic speed increases (and hence, the fast surface moves to a larger radii: Fig. 16a), which means that the wave takes less time to travel across the jet. Therefore the oscillations have smaller wavelengths and the jet exhibits a small standard deviation in the shape over time. On the other hand, if the jet base magnetisation is low enough, current driven instabilities are not fully triggered, which is the case for model B3, where the jet is mildly magnetised and pinches are absent.

6.4 Gas entrainment and jet mass-loading

In most of our simulations, the specific total energy flux $\mu$ oscillates and drops by a small amount in the pinched jet region (see e.g., Fig. 12). Pinching forces the gas to move across field lines in a non-uniform fashion, which disrupts the jet’s outward movement as well as causes mass-loading. Figure 17 shows the effect of jet mass-loading over time, as the B10 upper jet changes substantially over $t=(2-5)\times 10^{5}t_{g}$. Namely, eddies trap matter in the disc-wind and travel inwards through the jet boundary during pinching, forming finger-like structures. These fingers bend as they interact with the fast moving jet interior, dissipating poloidal field lines through reconnection (Figure 17, middle panels), and finally depositing surrounding gas into the jet body555Movie of mass-loading via gas entrainment in the B10 jet model: https://youtu.be/1aBoNormcS0. Figure 18 shows that the mass-flux through the jet indeed increases over time. It will be interesting to test whether explicit resistivity (e.g., Ripperda et al. 2019; we rely on numerical dissipation in H-AMR) brings any changes to the mass-loading in jets.

The entrainment mechanism we see here might be a manifestation of the Kruskal Schwarzschild instability (KSI, Kruskal & Schwarzschild 1954), a magnetised analogue of the Rayleigh-Taylor instability. Possibly, as the jet gets pinched, acceleration towards the axis (which can be seen as an effective gravity term) causes the heavy wind to push against the jet funnel leading to the onset of the KSI and producing the finger-like structures protruding into the jet. The growth rate of KSI is proportional to the square-root of the effective gravity term (Lyubarsky, 2010; Gill et al., 2017), which increases only when the pinches are fully developed. The small gravity term is why KSI may appear at such a late stage ($t>10^{5}r_{g}$) in the simulation. The fingers, when extended long enough, may be susceptible to secondary Kelvin-Helmholtz instabilities and bend, much like the distortion of a gas blob immersed in a magnetised fluid under gravity, shown in Fig. 12 of Gill et al. (2017). These bent fingers eventually collapse on themselves and lead to magnetic field reconnection. In fact, laboratory experiments of a plasma jet by Moser & Bellan (2012) also exhibit similar hybrid kink-KSI behaviour, which sets off magnetic reconnection events.

Figure 19 shows how gas entrainment and jet mass-loading affects model B10 quantities shown in Fig. 12. Compared to the wall-jet, specific energies $\mu$, $\sigma$ and $h$ drop and the jet slows down (panel a) as more and more gas enters the jet over time. From the marginal change in the field line shape (panel b) and the drop in the bunching parameter (panel c) as compared to Fig. 12(b,c), we conclude that the entrainment leads to significant poloidal field reconnection. Ultimately, due to the increasing jet density, the transverse profile of the Lorentz factor undergoes a dramatic change, as the slow sheath region decelerates from $\gamma\sim 3$ at $t=2\times 10^{5}t_{g}$ (Fig. 7, see also Sec. 4.1.4) to $\gamma\sim 1-2$ at $t>5\times 10^{5}t_{g}$, suggesting a dynamically changing sheath layer.

In order to understand how the jet behaves in various regimes, we have looked at jet quantities along field lines in different parts of the jet in model B10. In Sec. 4, we presented jet properties close to the axis (the foot-point jet opening angle $\theta_{j,\rm H}=0.44$ rad), where the mass-loading is low. We find that even though the inner jet, experiences a slowdown due to pinching, it achieves the peak Lorentz factor $\sim 5$. If we look at the mid-jet (Sec. 5.2, $\theta_{j,\rm H}=0.9$ rad), mass-loading is more prominent between $10^{3}-10^{4}r_{g}$ and eventually brings down $\gamma$ from 3 to near non-relativistic values at $t\sim(2-5)\times 10^{5}t_{g}$. When we look at the jet as a whole, we find that the distribution of gas from mass-loading and jet velocities is similar to the two-component (spine and sheath) structure of jets (e.g., Ghisellini et al., 2005; Ghisellini & Tavecchio, 2008) deduced from limb brightening in AGN jet observations (e.g., Giroletti et al., 2004; Nagai et al., 2014; Hada et al., 2016; Kim et al., 2018b).

7 Conclusions

In this work, we use our new GPU-accelerated GRMHD H-AMR code to investigate the largest extent disc-jet simulation performed till now, reaching over over 5 orders of magnitude in both distance and time in ultra high resolution. We start with a magnetised accretion disc around a spinning black hole that launches and accelerates a jet. This jet is self-consistently collimated by the disc-wind with the support of a large disc and qualitatively resembles the shape and acceleration profile of the M87 jet (Sec. 6.1). We find that the highly collimated jet maintains lateral causal connectivity with $\gamma\theta\lesssim 0.1$, consistent with VLBI observations of AGN jets ( 6.1).

Instead of the smooth outflow produced by jets in idealised models (e.g., Komissarov et al., 2007; Tchekhovskoy et al., 2010b), the interacting jet-wind interface exhibits oscillations from very near the jet origin. These oscillations appear to drive pinch instabilities at the jet’s outer boundary when the jet becomes super-fast (Sec. 6.3). Pinch instabilities significantly affect jet dynamics as they not only heat up the jet via magnetic reconnection, creating a thermal pressure gradient, but also lead to mass loading of the jet, both of which decelerate the jet. The dissipation due to pinching may lead to particle acceleration, potentially explaining non-thermal synchrotron and inverse Compton hot-spots in jets (e.g., Narayan et al., 2011; Sironi et al., 2015; Christie et al., 2018). The mass-loading, over time, helps to form a distinct slow moving layer with Lorentz factor $\gamma\sim 1$ that surrounds an inner fast moving jet core $\gamma\gtrsim 4$ (Sec. 6.4), resembling the spine-sheath structure seen in AGN jets (e.g., Kim et al., 2018b).

The instability causing the pinch modes is important to understand, since it may excite kink instabilities in 3D. As the jet kinks, it will interact with the ambient medium, leading to enhanced dissipation (e.g., Begelman, 1998; Giannios & Spruit, 2006), which disrupts the jet and potentially explains the FR-I/FR-II dichotomy (e.g., Bromberg & Tchekhovskoy, 2016; Tchekhovskoy & Bromberg, 2016; Barniol Duran et al., 2017; Liska et al., 2019). Our future work will thus focus on extending these results to full 3D, utilising the AMR capability of H-AMR to focus the resolution on the dissipative regions.

Acknowledgements

We thank O. Bromberg, A. Philippov, O. Porth and E. Sobacchi for helpful discussions. KC thanks M. Lucchini, T. Beuchert, F. Krauß and T.D. Russell for providing useful insights for the background study. This research was made possible by NSF PRAC award no. 1615281 and OAC-1811605 at the Blue Waters sustained-petascale computing project and supported in part under grant no. NSF PHY-1125915. KC and SM are supported by the Netherlands Organization for Scientific Research (NWO) VICI grant (no. 639.043.513), ML is supported by the NWO Spinoza Prize and AT by Northwestern University, the TAC and NASA Einstein (grant no. PF3-140131) postdoctoral fellowships and the National Science Foundation grant 1815304 and NASA grant 80NSSC18K0565.

The simulation data presented in this work is available on request to AT at atchekho@northwestern.edu.

Appendix A Local adaptive time-stepping

We have implemented in our block based AMR code H-AMR (Liska et al., 2018c) a so-called local adaptive time-stepping (LAT) routine. In addition to evolving higher spatial refinement levels with a smaller time-step, similar to AMR codes with a hierarchical time-stepping routine, LAT can also use different timesteps for blocks with the same spatial refinement level. Since most GRMHD simulations utilise a logarithmic spaced spherical grid, where cell sizes are small close to the black hole and large further away from the black hole, this can speed up the simulation by an additional factor $3-5$. In a future publication we will describe the detailed implementation of the LAT algorithm and show excellent scaling on pre-exascale GPU clusters.

LAT can also increase numerical accuracy by reducing the number of conserved to primitive variable inversions (Noble et al., 2006). Namely, as one moves away from the black hole on a logarithmic spaced spherical grid, the timescale of the problem increases. Evolving the outer grid with the same timestep as the inner grid leads to many unnecessary variable inversions. To illustrate that this produces noise in the outer grid, we produced two simulations with the same initial conditions at a resolution of $640\times 256\times 1$ (as lower resolutions naturally produce more noise): one with LAT enabled (Model B10-SLR) and one with LAT switched off (named B10-SLRL). Figure 20 shows that there are unphysical spikes in the Lorentz factor for model B10-SLRL caused by variable inversion failures in the outer jet region, which are absent in model B10-SLR. This confirms that LAT has the potential to increase speed and numerical accuracy.

Appendix B Grid shape

We have designed a grid that can track the shape of the jet over $5$ orders of magnitude in distance (Fig. 21, top). Furthermore, the grid keeps the cell aspect ratio in the outer jet below $10$ (Fig. 21, bottom), such that turbulent eddies at the disc-jet boundary remain resolved. The internal grid coordinates ($x_{0}$, $x_{1}$, $x_{2}$, $x_{3}$) are related to the real physical coordinates ($t$, $r$, $\theta$, $\phi$) as follows,

[TABLE]

where $A_{1}=[1+g_{1j}(\log_{10}r)^{g_{2j}}]^{-1}$, $A_{2}=g_{3j}\log_{10}r+g_{4j}$ and $A_{3}=0.5^{1-A_{2}}$. The parameters $g_{1j}=0.8$, $g_{2j}=3.0$, $g_{3j}=0.5$ and $g_{4j}=1.0$ are used to focus resolution on the jet, while $n_{r}=0.95$ is used to focus extra resolution on outer parts of the grid.

Appendix C Convergence of jet properties

Here we compare the time and spatial evolution of jets produced by the fiducial model B10 and its exceedingly high resolution version B10-HR. B10-HR has a resolution of $18000\times 1200\times 1$ and extends till $10^{5}r_{g}$, with all other parameters the same as B10. Figure 22 shows that the accretion properties of the B10 disc-jet system converges very well with respect to B10-HR. However, Fig. 23 shows that the jet in B10-HR gets mass-loaded ($\mu$ drops) earlier than the jet in B10, presumably because mass-loading is more efficient at higher resolutions, which capture the small scale eddies better.