AGN jet acceleration

Accelerating AGN jets to parsec scales using general relativistic MHD simulations

K. Chatterjee, M. Liska, A. Tchekhovskoy, S.B. Markoff
Anton Pannekoek Institute for Astronomy, University of Amsterdam, Science Park 904, 1098 XH Amsterdam, The Netherlands
Center for Interdisciplinary Exploration & Research in Astrophysics (CIERA), Physics & Astronomy, Northwestern University, Evanston, IL 60202, USA
Departments of Astronomy and Physics, Theoretical Astrophysics Center, University of California Berkeley, Berkeley, CA 94720-3411, USA
Lawrence Berkeley National Laboratory, 1 Cyclotron Rd, Berkeley, CA 94720, USA
Kavli Institute for Theoretical Physics, Kohn Hall, University of California at Santa Barbara, Santa Barbara, CA 93106, USA
Gravitation Astroparticle Physics Amsterdam (GRAPPA) Institute, University of Amsterdam, Science Park 904, 1098 XH Amsterdam, The Netherlands
Accepted XXX. Received YYY; in original form ZZZ

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

galaxies: black hole physics – accretion, accretion discs, jets – galaxies: individual (M87) – magnetohydrodynamics (MHD) – methods: numerical
pubyear: 2019pagerange: Accelerating AGN jets to parsec scales using general relativistic MHD simulationsC

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 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 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 . In its current version, in addition to CPUs, H-AMR also runs on graphical processing units (GPUs), achieving zone cycles per second on an NVIDIA Tesla V100 GPU and shows excellent parallel scaling to thousands of GPUs.

Full Model Short Floor Resolution factor Initial field loop LAT
name name [] [] [] [] configuration
R36A93B10 B10 36 73.97 10 5.0 62,14 multiple on
R36A93B10-S B10-S 36 73.97 10 1.2 96, 24 single on
R25A93B10 B10-R 25 50 10 1.2 36, 10 multiple on
R36A93B10-S-LR B10-SLR 36 73.97 10 0.8 20, 4 single on
R36A93B10-S-LR-L B10-SLRL 36 73.97 10 0.8 20, 4 single off
R36A93B10-HR B10-HR 36 73.97 10 0.5 -,- multiple on
R36A93B3 B3 36 73.97 3 1.5 32,10 multiple on
R36A93B50 B50 36 73.97 50 1.9 32,10 multiple on
R36A93B100 B100 36 73.97 100 2.0 40, 12 multiple on
Table 1: Simulations carried out for this work. Model names shorthand disc inner radius (R), black hole spin (A; , common for all models) and the floor value (B; in other words, the maximum magnetisation ) 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 , outer grid radius , grid resolution, simulation time , MRI quality factor at (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).

2.1 Numerical grid

We use units such that . This sets both the characteristic timescale, , and spatial scale, , to unity. In fact, our simulations are scale-free, i.e., if we provide the black hole mass and mass accretion rate , we can rescale our simulation to the corresponding black hole system. Our grid is axisymmetric, extending from to , where is the event horizon radius, , where we set the dimensionless black hole spin parameter to . We carry out the simulations on a uniform grid in internal coordinates (, see Appendix B) that are transformations of the spherical polar coordinates () in the Kerr-Schild foliation specifically optimised to follow the collimating jets. To resolve the jets, we typically use a numerical resolution of cells in the radial direction and cells in the polar direction (Table 1). We use the following boundary conditions (BCs): in , we use outflow BCs at the inner and outer grid radii; at the poles, we reflect the 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 and a minimum internal energy of , where , , and 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 . In all models we start with a Fishbone & Moncrief (1976) (FM76) torus in hydrostatic equilibrium around a rapidly spinning Kerr black hole. We place the torus inner edge at and density maximum at . For all but one of our simulations, we set and . We choose the density scale by setting . We adopt the ideal gas law equation of state and set the adiabatic index to that of a non-relativistic gas . We set the magnetic field vector potential as described in Secs. 3.1 and 3.2 and normalise it such that , where and are the gas and magnetic pressure, respectively.

3.1 Single field loop setup

Figure 1: We start out with a magnetised disc around a spinning black hole (at origin). We show two different initial conditions for our accretion disc with the fluid-frame density in colour (red shows high and blue low density; see the colour bar), representing an equilibrium hydrodynamic torus around the spinning black hole, with black lines showing the initial magnetic field configuration of a large poloidal field loop (left panel, model B10-S) and two poloidal loops of opposite polarity (right panel, model B10). Dashed black lines show field lines containing negative magnetic flux. The simulation grid extends out to or , depending on the model (see Table 1).
Figure 2: Large magnetic flux near the black hole can stop accretion from the disc. We show a time snapshot of the model B10-S within the innermost ten gravitational radii. The central black hole is shown in black along with the rest mass matter density (in colour) and the magnetic field lines (in black). Black hole gravity attempts to pull in matter from the accretion disc (orange-red region), while the accumulated strong magnetic flux pushes gas away, resulting in the magnetically arrested disc (MAD) state. In axisymmetry, the accretion rate is highly variable for MADs, which in turn affects the jet (Fig. 3).

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:


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: . Here, , where is the energy accretion rate (defined to be positive when the energy flows into the black hole, see Eq. 3 for definition), is the surface area element, and is the determinant of the metric. Similarly, we can define the magnetic flux on the black hole as , where the integral is over the entire event horizon. We also define its dimensionless counterpart normalised by the mass accretion rate, .

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 (Fig. 3a), dimensionless magnetic flux (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.

Figure 3: The MAD model B10-S (orange) shows an oscillating mass accretion rate (panel a), magnetic flux (panel b) and total outflow power (panel c). In contrast, model B10 (blue) shows steady behaviour in all. The normalised is computed over . Free of violent variability, model B10 is suitable for studying large scale jet dynamics.

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,


where . The first term in 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: and ; 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 , dimensionless magnetic flux 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.

Figure 4: The disc in model B10 launches opposing jets with field lines anchored in the black hole event horizon as well as the disc. Mixing between the jet and the disc-wind mass-loads the jet over time via eddies generated from the wind-jet interaction. We show vertical slices though the density (see the colour bar) at an early time (left panel) and at late time (right panel). Pinch instabilities, in the form of finger-like projections, significantly contribute to mass-loading and plays a vital role in determining jet dynamics. In order to gauge the influence of pinching on the jet, we extract useful information about energetics along a field line (indicated by the pink line) in both jets as shown in Fig. 5.
Figure 5: Magnetised jets accelerate by converting Poynting and thermal energy into kinetic energy: Lorentz factor increases at the expense of decreasing magnetisation and specific enthalpy , while maintaining a near constant specific total energy flux . We show the radial profile of these quantities along a field line with a foot-point half-opening angle  rad (Fig. 4, pink) for both the upper (Fig. 4, ) and lower (Fig. 4, ) jets (solid and dash-dotted lines, respectively) of model B10 at . The two jets do not have the same acceleration profile, as the upper jet is affected by strong pinches, which lead to gas moving across the field lines in a non-uniform way and contributes to mass-loading the jet. The added inertia in the jet causes a drop in the Lorentz factor and a rise in the specific enthalpy. We radially average the plotted quantities over .

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 , where measures the number of cells per MRI wavelength in direction [], where is the Alfvén velocity, the cell size, the angular velocity of the fluid. is averaged over the inner disc () and weighted by . We achieve (Table 1) at , fulfilling the numerical convergence criteria (see e.g., Hawley et al., 2011). 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: (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 ( in Fig. 4) is strongly affected by the interface instabilities, while the lower jet ( 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, , where is the poloidal field strength. We choose a field line whose foot-point makes an angle of  rad with the black hole spin axis at the event horizon. This makes up about % of the jet’s half-opening angle . Figure 5 shows the evolution of various quantities along our chosen field line at , a long enough time for the jet to establish a quasi-steady state solution out to . First, we consider the total specific energy , 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, , and the rest-mass flux, ,


where is gas density, is internal energy density, is gas pressure, and are the magnetic and velocity four-vectors respectively. Figure 5 shows that the radial profile of 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 around is caused by the presence of a stagnation surface (see Sec. 2.2), where and thus . Downstream of the stagnation surface, the floor values (see Sec. 2.2) set the value of . To quantify the conversion efficiency of magnetic to kinetic energy, we calculate the ratio of the Poynting flux, , to the mass energy flux, , called magnetisation :


where for . Figure 5 shows that 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 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.

Figure 6: The collimation of the jet has a profound effect on the acceleration profile. Continuing from Fig. 5, here we show more quantities along the field line for model B10. Panel(a) shows the jet half-opening angle in radians, (b) the bunching parameter (see text), (c) the transverse causality parameter and (d) a second causality parameter , useful for observed jets. Both jets show a continuous parabolic collimation profile. Strong pinching in the upper jet forces reconnection between two nearby field lines, and result in poloidal flux dissipation. The lower jet is also affected by pinches but they are weaker comparatively, as illustrated by the difference in the Lorentz factor between the two jets at approximately as shown in Fig. 5. Both jets are casually connected throughout their length, in agreement to VLBI images of AGN jets (e.g., Jorstad et al., 2005).
Figure 7: Transverse cross-sections of the jet at different distances (see the legend) show that as the jet accelerates, the poloidal flux surfaces differentially bunch up towards the axis and build up a fast inner jet core. The figure shows (panel a) the specific total energy , (panel b) magnetisation , (panel c) Lorentz factor and (panel d) the specific enthalpy at at different distances ). We take the jet-edge to be . The corresponding jet-edge half-opening angles () for the different distances are (, , , ) rad. We also indicate with circles the opening angle of the field line shown in Fig. 5. The peak in the profile shifts towards the jet axis with increasing distance as a result of differential field line bunching. The jet-edge experiences mass-loading from the wind and thus with increasing distance, the specific energies decrease at the edge. Pinching causes magnetic dissipation and hence the specific enthalpy tends to increase with distance (see also Fig. 5).


Figure 8: Jets with higher values of accelerate to higher Lorentz factors. We compare the jet acceleration profile of models B3 (dashed double dotted), B10 (solid), B50 (dashed) and B100 (dotted), with maximum magnetisation of 3, 10, 50 and 100 respectively, i.e., only varying the jet base magnetisation (see Sec. 2.2). We show the evolution of specific total energy , magnetisation and the Lorentz factor along field lines in the mid-jet ( rad) at . Due to a larger jet base magnetisation, B100 accelerates to while B3 only accelerates to . Except for B3, pinching significantly affects all models around . Evidently, the jet base magnetisation plays a role in determining pinch activity.

4.1.1 Jet structure and acceleration

Figure 6(a) shows that both the upper and the lower jets collimate similarly from  rad (or ) at to  rad (or ) at while displaying a power-law shape . 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 for both the lower and upper jets. Such a 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, , where is the conserved field line angular velocity. The Lorentz factor in this regime behaves as , where 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),


where we define the field line bunching parameter as the ratio of the local poloidal field strength and the mean poloidal field strength, . 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 does decrease and therefore, increases to values approaching, within a factor of few, the maximum Lorentz factor . The Lorentz factor in this simulation reaches , 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 must be smaller than the Mach cone half-opening angle (, where is the fast Mach number). Using the definition of the fast magnetosonic wave velocity and Lorentz factor (e.g., Gammie et al., 2003),


and focusing on the asymptotic regime of the jet, i.e., and , we can compute the Mach cone half-opening angle as,


Approximating the flow velocity (for a relativistic jet), we have . Figure 6(c) shows the transverse causality parameter , the ratio of the jet and Mach cone half-opening angles. For , 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 (Fig. 6d). The acceleration slow down occurs as the jet becomes matter-dominated, i.e., , 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 (Fig. 5, black line) decreases initially, slightly at . However, beyond , 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 . This pressure gradient significantly slows down the outflow to by . At larger radii, 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 , from Eq. (3) we get


a useful form of the energy equation. It clearly shows that for constant, an increase in enthalpy to results in decreasing, which is seen for the upper jet around . The upper jet shows stronger pinch activity and thus collimates slightly more than the lower jet: is smaller by a factor of ; Fig. 6a, solid line), and other quantities like the maximum Lorentz factor and bunching parameter strongly oscillate.

4.1.4 Energetics across the jet

Figure 7 shows how the different components of jet energy flux, , , and , vary across the jet at different distances along the upper jet. A fiducial field line with the foot-point at  rad (see also Fig. 5), shown with filled circles, collimates faster than the jet-edge, indicated by . We take the jet-edge to be at , which is reasonable since should drop rapidly at the jet-edge 333We also note that is roughly constant throughout the width of the jet as opposed to increasing as , 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 , the magnetisation and the Lorentz factor to drop gradually, forming a slower sheath that surrounds the core, resulting in a structure similar to the spine-sheath seen in AGN jets.

Figure 9: Vertical slices through the density profile model B10-R, which has a small disc, show that smaller discs have wider jets and broader wind regions as compared to larger discs (compare to Fig. 4). Density and field lines are labelled as in Fig.4. The initial disc (left panel) is much smaller than in Fig. 4, enabling the jet to freely expand laterally into the low density ambient medium as seen in the right panel at . As the jet undergoes rapid lateral expansion, pinches do not have enough time to develop and therefore, the jet mass-loading is much lower (see the main text). We highlight two representative field lines with pink and magenta colours and show their properties in Fig. 10.

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 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 increases steeply until , beyond which the acceleration slows down. Models with higher jet base magnetisation (and hence, larger values) accelerate slightly faster and reach higher Lorentz factors, but get affected by pinch instabilities at roughly the same distance as models with lower . 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 achieve smaller at large distances. For instance, model B3 achieves at . It is interesting to note that efficient heating via relativistic shocks requires (e.g., Kennel & Coroniti, 1984; Komissarov, 2012). Such low values are very difficult to achieve for collimated flows as drops very gradually in the toroidally dominated regime (known as the 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.

Figure 10: Radial profiles of quantities along two magnetic field lines in the upper jet at for model B10-R, one in the inner jet ( rad, solid lines; the field line is highlighted in pink in Fig. 9, right panel) and the other in the outer jet ( rad, dashed-dotted; magenta in Fig. 9, right panel). Refer to Fig. 5 and 6 for the notations used. Panel(a): The outer field line accelerates slightly faster than the inner one. The enthalpy remains on average beyond . Panel(b): The jet is initially parabolic and becomes conical beyond , with becoming constant for the outer field line. This expansion causes the outer jet to exhibit (panel c) a sudden drop in the bunching parameter and (panel d) loss of causal contact in the lateral direction as . Panel(e): The inner jet remains causally connected while the outer jet becomes conical with . Thus, deconfinement has a notable influence on jet dynamics, reinforcing the notion that jet acceleration is coupled to the jet collimation profile.

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: 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 , 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 . 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  rad; highlighted in pink in Fig. 9, right panel) and the other in the outer jet ( rad; magenta in Fig. 9, right panel), at . Initially, similar to our fiducial model B10, the disc and disc-wind collimate the jet into a parabolic shape. However, beyond , 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 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 , 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 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 compared to the inner field line (Fig. 10c), in accordance with Eq. (5) and reaches (Fig. 10a). Upon loss of collimation, the outer jet also loses transverse causal connection (; Fig. 10d), and the acceleration ceases (see Sec. 4.1.2). From Fig. 10(e), for the outer jet, while is between 0.1 and 0.4 within for the inner jet. See Sec. 6.1 for further discussion of 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.

Figure 11: Comparison between idealised wall-jet and disc-jet simulations shows that instabilities at jet-disc interface slow down the outflow. Left: Lorentz factor plot of a jet bound by a rigid parabolic wall. Right: combined Lorentz factor and density plot for disc-jet model B10. For the case of the wall-jet, the field line shape is smooth and the outflow quickly accelerates. However, for the disc-jet, the pinch instabilities distort the shape of the field line and slow down acceleration. This is better seen when we compare quantities along the indicated field lines (cyan) in Fig. 12.
Figure 12: Comparing the radial profile of quantities along field lines for the upper jet in disc-jet model B10 (solid lines,  rad) and the idealised wall-jet model (double dot-dashed lines,  rad) shows that the two models agree apart from the pinch instabilities that slow down the disc-jets and dissipate magnetic fields into heat. The disc-jet simulation is shown at , while the wall-jet is in steady state. Panel(a): The specific energy profiles (, , and ) for B10 match the idealised model reasonably well, with deviations arising in the pinched region of model B10, especially in the enthalpy . Panel(b): Both field lines have similar collimation profiles. Panel(c): Pinching causes the poloidal field in model B10 to fluctuate and dissipate into heat, thereby decreasing the bunching parameter.The jet mass-loading also plays a part as decreases with respect to the wall-jet . Panels(d and e): The values of and remain below 1, indicating lateral causal connection for both jets. Overall, the radial profiles between the two models match well except in the disc-jet pinched region.

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:


This gives us the initial poloidal field configuration:


Here, the outermost field line touching the wall is given by : it starts out at the intersection of the event horizon () and the equatorial plane (), is initially radial for and asymptotically collimates as . In this setup, results in a monopolar field shape (constant), while gives us the parabolic field shape (). We set a transitional radius and employ as these values give a good match to the disc-jet shape, as we discuss below. We define the physical coordinates (, ) as functions of the internal coordinates (, ) as and . We use a resolution of cells. Our computational domain range extends radially from to . We employ the same polar reflective boundary conditions at as for the disc-jet simulations (see Sec. 2.1). At the wall, , 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 , 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 , which ensures that the outflow reaches a steady state up to at least a distance of . To speed up the simulation, from onwards, we freeze out cells that reached steady state, i.e., the cells located at , where 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, , , and 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 decreases with increasing as expected due to adiabatic jet expansion. For the disc-jet, increases substantially at 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 until a few times , 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 to dissipate (Fig. 12c), along with a slight drop in (Fig. 12d). The product of the Lorentz factor and jet opening angle for (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

Figure 13: Comparison of the shape and Lorentz factor of a simulated jet with M87 observations shows remarkable resemblance. (Panel a) Jet radius along a field line near the jet-edge ( rad) for model B10 at . The data points are taken from Hada et al. (2011); Hada et al. (2013); Hada et al. (2016); Doeleman et al. (2012); Asada & Nakamura (2012); Nakamura & Asada (2013). The jet from model B10 fits very well with the M87 parabolic jet shape up to . (Panel b) Lorentz factor along panel(a) field line as well as a field line in the mid-jet ( rad) compared to the broken power-law profile for the M87 jet Lorentz factor as measured by VLBI (Mertens et al., 2016). There is a large distribution in across the jet, similar to Fig. 15 of Mertens et al. (2016).

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 , close to the location of HST-1 (Asada & Nakamura, 2012; Hada et al., 2013; Mertens et al., 2016; Kim et al., 2018b).

Mertens et al. (2016) showed using VLBI measurements that the acceleration profile of M87 follows till , changing to 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, the large observed jet radius close to the black hole might indicate that the emission comes from the disc, i.e., outside the jet-edge, where the Lorentz factor . The outer field line Lorentz factor is close to 1.5 (Fig. 7), similar to the velocities found for the M87 jet sheath (Mertens et al., 2016). 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 . One possible reason may be that HST-1 lies very close to the Bondi radius of M87 (; Nakamura & Asada, 2013) where the shallow density profile of the ISM prevails. If there is an increase of confining pressure from the ISM beyond , 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 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 () for distances larger than , compared to those observed (, 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 . It is possible that the difference in measurements of 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: (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.

Figure 14: (Panel a) Log-log plot of characteristic surfaces along with entropy in colour (arbitrary units) and magnetic field lines (black) for a highly magnetised jet time-averaged over in model B10. The lines shown are, starting from small radii, the event horizon (thick black), stagnation (magenta), Alfvén (green), classical fast magnetosonic surfaces (FMS; blue) and the fast magnetosonic separatrix surface (FMSS; salmon). We highlight a field line in silver and show the variability in its shape due to pinching in Fig. 16. The characteristic surfaces are only shown for the jet and the wind. The FMSS does not appear to coincide with dissipative features. (Panel b) We look at the entropy at , indicating the points where the fast surfaces for the jet (diamond) and the wind (circle) crosses the horizontal line, shown in panel(a). The entropy rises smoothly beginning from the sub-fast wind region, right up to the fast surface in the jet, suggesting that while the FMS might be relevant for dissipation in the jet, for the wind, the FMS is not so useful.
Figure 15: Sufficient resolution is required to properly capture the small scales of the pinches. We show the specific enthalpy along a field line in the inner jet for our fiducial model B10 and a low resolution model B10-SLR at . There is larger dissipation, i.e., higher , for B10-SLR as it is unable to resolve the micro-structures in the jet caused by toroidal pinch instabilities.
Figure 16: The jet shape changes over time due to the pinching between the jet and the disc-wind. (Panel a) We show the time averaged half-opening angle for field lines from different simulation models, with the same foot-point half-opening angle of  rad, along with their standard deviation (shaded area) over . To minimise crowding, the curves are shifted vertically. We indicate the position where the field line crosses the fast surface (FMS) with a circle. The field line for model B10 is shown in silver in Fig. 14. There is significant time variation in shape of the B10 and B50 field lines beyond the FMS due to the presence of toroidal pinch instabilities, which continue throughout the entire jet. The fast surface moves outward as the jet magnetisation increases (similar to Ceccobello et al., 2018). (Panel b) The relative deviation with respect to the time-averaged jet opening angle along the B10 field line shows that there is deviation in the pinched region. (Panel c) We show the entropy along all 3 field lines from Fig. 14, with  rad representing the inner, mid and outer jet respectively. We also indicate where the FMS crosses the field lines with circles. The rise in entropy appears to coincide very well with the FMS, except for the inner jet field line where the entropy rise may be due to round-off errors from machine-precision calculations end up affecting the smallest energy term, i.e., the internal energy, and consequently the gas pressure.

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 ,


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 self-similar models (e.g., Sauty et al., 2004).

Figure 17: Mixing due to pinch instabilities between the jet and the disc-wind leads to mass-loading of the jet. We show a vertical slice of the B10 upper jet-wind system with density in colour and magnetic field lines in black, at (left) and (right). We indicate the field line shown in Fig. 12 with pink. Inset panels: Zoom-in snapshots of the jet over at different times. Over time, pinch instabilities at the jet-wind interface set off reconnection events and mass-load the jet. The increase in density changes the energy distribution along the jet, reducing its specific energy content and greatly affecting the jet’s acceleration profile (Fig. 19).
Figure 18: Gas captured from the disc-wind via entrainment increases the mass flux of the jet, while the total energy flux is conserved on average, thus decreasing the specific total energy flux and reducing the jet energy budget. We show the jet total energy flux and the rest-mass energy flux over the length of the upper jet at and . The jet averaged fluxes are calculated as using the criterion for the magnetised jet.

6.3 Origin of pinch instabilities

Whereas previous GRMHD simulations (e.g., McKinney, 2006) found that pinch instabilities, forming at the jet-wind interface, do not survive beyond , 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 ( 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 (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 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 direction, suggesting ram pressure in the 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. 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). 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 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 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).

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,


assuming a cold jet, i.e., specific enthalpy . At the fast surface (), we then have for the jet,


Assuming that and , we arrive at the location of the fast surface,


From (14), we can indeed say that increases with increase in for a given field line. For , Eq. (13) gives which is about off from the simulation value (). However, the assumption of is not valid in the outer jet, due to stronger effects of mass-loading. The time variability of 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.

Figure 19: Over time, the disc-wind mixes with the jet and mass-loads it, significantly slowing the jet. We show the same as Fig. 12, except at a later time, . Panel(a): There is a large drop in the specific energy profiles for the B10 upper jet as compared to earlier (Fig. 12), directly affecting the acceleration profile. drops by almost an order of magnitude due to the increase in mass flux. Panel(b) shows that the jet collimation profile does not change by much. Panel(c): The poloidal field, on the other hand, drops due to reconnection of pinched field lines, reducing the bunching parameter. Since the jet slows down, causality is still maintained (panel e and f).

6.4 Gas entrainment and jet mass-loading

In most of our simulations, the specific total energy flux 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 . 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: 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. 2018; 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 () 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 , and 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 at (Fig. 7, see also Sec. 4.1.4) to at , 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  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 . If we look at the mid-jet (Sec. 5.2,  rad), mass-loading is more prominent between and eventually brings down from 3 to near non-relativistic values at . 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 , 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 that surrounds an inner fast moving jet core (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-AMRto focus the resolution on the dissipative regions.


We thank M. Nakamura, K. Asada and K. Hada for providing M87 data for our figures. We also thank O. Bromberg, A. Philippov and O. Porth for helpful discussion. 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 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


  • Anglés-Alcázar et al. (2017) Anglés-Alcázar D., Faucher-Gigure C.-A., Quataert E., Hopkins P. F., Feldmann R., Torrey P., Wetzel A., Kere D., 2017, MNRAS, 472, L109
  • Asada & Nakamura (2012) Asada K., Nakamura M., 2012, ApJL, 745, 5 pp
  • Balbus & Hawley (1991) Balbus S. A., Hawley J. F., 1991, ApJ, 376, 214
  • Barniol Duran et al. (2017) Barniol Duran R., Tchekhovskoy A., Giannios D., 2017, MNRAS, 469, 4957
  • Begelman (1998) Begelman M. C., 1998, ApJ, 493, 291
  • Begelman & Li (1994) Begelman M. C., Li Z.-Y., 1994, ApJ, 426, 269
  • Beskin & Kuznetsova (2000) Beskin V. S., Kuznetsova I. V., 2000, Nuovo Cimento B Serie, 115, 795
  • Beskin & Nokhrina (2006) Beskin V. S., Nokhrina E. E., 2006, MNRAS, 367, 375
  • Beskin et al. (1998) Beskin V. S., Kuznetsova I. V., Rafikov R. R., 1998, MNRAS, 299, 341
  • Biretta et al. (1999) Biretta J. A., Sparks W. B., Macchetto F., 1999, ApJ, 520, 621
  • Blandford & Königl (1979) Blandford R. D., Königl A., 1979, ApJ, 232, 34
  • Blandford & Payne (1982) Blandford R. D., Payne D. G., 1982, MNRAS, 199, 883
  • Blandford & Znajek (1977) Blandford R. D., Znajek R. L., 1977, MNRAS, 179, 433
  • Bourne & Sijacki (2017) Bourne M. A., Sijacki D., 2017, MNRAS, 472, 4707
  • Bower et al. (2006) Bower R. G., Benson A. J., Malbon R., Helly J. C., Frenk C. S., Baugh C. M., Cole S., Lacey C. G., 2006, MNRAS, 370, 645
  • Broderick & Loeb (2006) Broderick A. E., Loeb A., 2006, MNRAS, 367, 905
  • Broderick & Tchekhovskoy (2015) Broderick A. E., Tchekhovskoy A., 2015, ApJ, 809, 97
  • Bromberg & Levinson (2007) Bromberg O., Levinson A., 2007, ApJ, 671, 678
  • Bromberg & Tchekhovskoy (2016) Bromberg O., Tchekhovskoy A., 2016, MNRAS, 456, 1739
  • Ceccobello et al. (2018) Ceccobello C., Cavecchi Y., Heemskerk M. H. M., Markoff S., Polko P., Meier D. L., 2018, MNRAS, 473, 4417
  • Chen et al. (2018) Chen A. Y., Yuan Y., Yang H., 2018, ApJ, 863, L31
  • Chiueh et al. (1998) Chiueh T., Li Z.-Y., Begelman M. C., 1998, ApJ, 505, 835
  • Christie et al. (2018) Christie I. M., Petropoulou M., Sironi L., Giannios D., 2018, MNRAS, 482, 65
  • Clausen-Brown et al. (2013) Clausen-Brown E., Savolainen T., Pushkarev A. B., Kovalev Y. Y., Zensus J. A., 2013, A&A, 558, A144
  • Colella & Woodward (1984) Colella P., Woodward P. R., 1984, Journal of Computational Physics, 54, 174
  • Cowling (1934) Cowling T. G., 1934, MNRAS, 94, 768
  • De Villiers et al. (2003) De Villiers J.-P., Hawley J. F., Krolik J. H., 2003, ApJ, 599, 1238
  • Doeleman et al. (2012) Doeleman S. S., et al., 2012, Science, 338, 355
  • Drenkhahn & Spruit (2002) Drenkhahn G., Spruit H. C., 2002, A&A, 391, 1141
  • Eichler (1993) Eichler D., 1993, ApJ, 419, 111
  • Fabian (2012) Fabian A. C., 2012, Annual Review of Astronomy and Astrophysics, 50, 455
  • Fishbone & Moncrief (1976) Fishbone L. G., Moncrief V., 1976, ApJ, 207, 962
  • Fromm et al. (2017) Fromm C. M., Porth O., Younsi Z., Mizuno Y., De Laurentis M., Olivares H., Rezzolla L., 2017, Galaxies, 5
  • Gammie et al. (2003) Gammie C. F., McKinney J. C., Tóth G., 2003, ApJ, 589, 444
  • Gardiner & Stone (2005) Gardiner T. A., Stone J. M., 2005, Journal of Computational Physics, 205, 509
  • Ghisellini & Tavecchio (2008) Ghisellini G., Tavecchio F., 2008, MNRAS, 385, L98
  • Ghisellini et al. (2005) Ghisellini G., Tavecchio F., Chiaberge M., 2005, A&A, 432, 401
  • Giannios & Spruit (2006) Giannios D., Spruit H. C., 2006, A&A, 450, 887
  • Gill et al. (2017) Gill R., Granot J., Lyubarsky Y., 2017, MNRAS, 474, 3535
  • Giroletti et al. (2004) Giroletti M., et al., 2004, ApJ, 600, 127
  • Granot et al. (2011) Granot J., Komissarov S. S., Spitkovsky A., 2011, MNRAS, 411, 1323
  • Hada et al. (2011) Hada K., Doi A., Kino M., Nagai H., Hagiwara Y., Kawaguchi N., 2011, Nature, 477, 185
  • Hada et al. (2013) Hada K., et al., 2013, ApJ, 775, 70
  • Hada et al. (2016) Hada K., et al., 2016, ApJ, 817, 131
  • Harrison et al. (2018) Harrison C. M., Costa T., Tadhunter C. N., Flütsch A., Kakkad D., Perna M., Vietri G., 2018, Nature Astronomy, 2, 198
  • Harten (1983) Harten A., 1983, Journal of Computational Physics, 49, 357
  • Hawley & Krolik (2006) Hawley J. F., Krolik J. H., 2006, ApJ, 641, 103
  • Hawley et al. (2011) Hawley J. F., Guan X., Krolik J. H., 2011, ArXiv:1103.5987,
  • Hirotani & Okamoto (1998) Hirotani K., Okamoto I., 1998, ApJ, 497, 563
  • Hirotani & Pu (2016) Hirotani K., Pu H.-Y., 2016, ApJ, 818, 50
  • Igumenshchev et al. (2003) Igumenshchev I. V., Narayan R., Abramowicz M. A., 2003, ApJ, 592, 1042
  • Jorstad et al. (2005) Jorstad S. G., et al., 2005, AJ, 130, 1418
  • Jorstad et al. (2017) Jorstad S. G., et al., 2017, ApJ, 846, 98
  • Kennel & Coroniti (1984) Kennel C. F., Coroniti F. V., 1984, ApJ, 283, 694
  • Kim et al. (2018a) Kim J., Balsara D. S., Lyutikov M., Komissarov S. S., 2018a, MNRAS, 474, 3954
  • Kim et al. (2018b) Kim J.-Y., et al., 2018b, A&A, 616, A188
  • Komissarov (1994) Komissarov S. S., 1994, MNRAS, 266, 649
  • Komissarov (2001) Komissarov S. S., 2001, MNRAS, 326, L41
  • Komissarov (2005) Komissarov S. S., 2005, MNRAS, 359, 801
  • Komissarov (2012) Komissarov S. S., 2012, MNRAS, 422, 326
  • Komissarov et al. (2007) Komissarov S. S., Barkov M. V., Vlahakis N., Königl A., 2007, MNRAS, 380, 51
  • Komissarov et al. (2009) Komissarov S. S., Vlahakis N., Königl A., Barkov M. V., 2009, MNRAS, 394, 1182
  • Komissarov et al. (2010) Komissarov S. S., Vlahakis N., Königl A., 2010, MNRAS, 407, 17
  • Komissarov et al. (2015) Komissarov S. S., Porth O., Lyutikov M., 2015, Computational Astrophysics and Cosmology, 2, 9
  • Kruskal & Schwarzschild (1954) Kruskal M. D., Schwarzschild M., 1954, Proc. Roy. Soc. (London), A223, 348
  • Levinson & Segev (2017) Levinson A., Segev N., 2017, Phys. Rev. D, 96, 123006
  • Li et al. (1992) Li Z.-Y., Chiueh T., Begelman M. C., 1992, ApJ, 394, 459
  • Liska et al. (2018b) Liska M., Tchekhovskoy A., Ingram A., van der Klis M., 2018b, preprint, (arXiv:1810.00883)
  • Liska et al. (2018a) Liska M., Tchekhovskoy A., Quataert E., 2018a, preprint, (arXiv:1809.04608)
  • Liska et al. (2018c) Liska M., Hesp C., Tchekhovskoy A., Ingram A., van der Klis M., Markoff S., 2018c, MNRAS, 474, L81
  • Liska et al. (2019) Liska M., Hesp C., Tchekhovskoy A., Ingram A., van der Klis M., Markoff S., 2019, preprint, (arXiv:1901.05970)
  • Lister et al. (2016) Lister M. L., et al., 2016, ApJ, 152, 12
  • Lucchini et al. (2018) Lucchini M., Krauß F., Markoff S., Crumley P., Connors R. M. T., 2018, MNRAS, 482, 4798
  • Lyubarsky (2009) Lyubarsky Y., 2009, ApJ, 698, 1570
  • Lyubarsky (2010) Lyubarsky Y., 2010, ApJ, 725, L234
  • Magorrian et al. (1998) Magorrian J., et al., 1998, ApJ, 115, 2285
  • Markoff (2010) Markoff S., 2010, From Multiwavelength to Mass Scaling: Accretion and Ejection in Microquasars and AGN. Springer Berlin Heidelberg, Berlin, Heidelberg, pp 143–172, doi:10.1007/978-3-540-76937-8_6,
  • Markoff et al. (2001) Markoff S., Falcke H., Fender R. P., 2001, Astron. Astrophys., 372, L25
  • Markoff et al. (2005) Markoff S., Nowak M. A., Wilms J., 2005, ApJ, 635, 1203
  • McKinney (2005) McKinney J. C., 2005, ApJ, 630, L5
  • McKinney (2006) McKinney J. C., 2006, MNRAS, 368, 1561
  • Meier (2012) Meier D. L., 2012, Black Hole Astrophysics: The Engine Paradigm. Springer, Verlag Berlin Heidelberg
  • Meier et al. (1997) Meier D. L., Edgington S., Godon P., Payne D. G., Lind K. R., 1997, Nature, 388, 350
  • Mertens et al. (2016) Mertens F., Lobanov A. P., Walker R. C., Hardee P. E., 2016, A&A, 595, A54
  • Mizuno et al. (2004) Mizuno Y., Yamada S., Koide S., Shibata K., 2004, ApJ, 615, 389
  • Mizuno et al. (2015) Mizuno Y., Gómez J. L., Nishikawa K.-I., Meli A., Hardee P. E., Rezzolla L., 2015, The Astrophysical Journal, 809, 38
  • Moll et al. (2008) Moll R., Spruit H. C., Obergaulinger M., 2008, A&A, 492, 621
  • Moser & Bellan (2012) Moser A. L., Bellan P. M., 2012, Nature, 482, 379
  • Nagai et al. (2014) Nagai H., et al., 2014, ApJ, 785, 53
  • Nakamura & Asada (2013) Nakamura M., Asada K., 2013, ApJ, 775, 118
  • Narayan et al. (2003) Narayan R., Igumenshchev I. V., Abramowicz M. A., 2003, PASJ, 55, L69
  • Narayan et al. (2011) Narayan R., Kumar P., Tchekhovskoy A., 2011, MNRAS, 416, 2193
  • Noble et al. (2006) Noble S. C., Gammie C. F., McKinney J. C., Del Zanna L., 2006, ApJ, 641, 626
  • Parfrey et al. (2019) Parfrey K., Philippov A., Cerutti B., 2019, Phys. Rev. Lett., 122, 035101
  • Plavin et al. (2019) Plavin A. V., Kovalev Y. Y., Pushkarev A. B., Lobanov A. P., 2019, preprint, (arXiv:1811.02544)
  • Polko et al. (2010) Polko P., Meier D. L., Markoff S., 2010, ApJ, 723, 1343
  • Porth & Komissarov (2015) Porth O., Komissarov S. S., 2015, MNRAS, 452, 1089
  • Proga & Zhang (2006) Proga D., Zhang B., 2006, MNRAS, 370, L61
  • Ptitsyna & Neronov (2016) Ptitsyna K., Neronov A., 2016, A&A, 593, A8
  • Pu et al. (2015) Pu H.-Y., Nakamura M., Hirotani K., Mizuno Y., Wu K., Asada K., 2015, ApJ, 801, 56
  • Pushkarev et al. (2009) Pushkarev A. B., Kovalev Y. Y., Lister M. L., Savolainen T., 2009, A&A, 507, L33
  • Pushkarev et al. (2017) Pushkarev A. B., Kovalev Y. Y., Lister M. L., Savolainen T., 2017, MNRAS, 468, 4992
  • Ressler et al. (2017) Ressler S. M., Tchekhovskoy A., Quataert E., Gammie C. F., 2017, MNRAS, 467, 3604
  • Ripperda et al. (2018) Ripperda B., Porth O., Sironi L., Keppens R., 2018, preprint, (arXiv:1810.10116)
  • Romero et al. (2017) Romero G. E., Boettcher M., Markoff S., Tavecchio F., 2017, Space Science Reviews, 207, 5
  • Russell et al. (2013) Russell D. M., et al., 2013, MNRAS, 429, 815
  • Sauty et al. (2004) Sauty C., Trussoni E., Tsinganos K., 2004, A&A, 421, 797
  • Sijacki et al. (2007) Sijacki D., Springel V., Di Matteo T., Hernquist L., 2007, MNRAS, 380, 877
  • Silk & Rees (1998) Silk J., Rees M. J., 1998, A&A, 331, L1
  • Sironi & Spitkovsky (2009) Sironi L., Spitkovsky A., 2009, ApJ, 698, 1523
  • Sironi & Spitkovsky (2014) Sironi L., Spitkovsky A., 2014, ApJ, 783, L21
  • Sironi et al. (2015) Sironi L., Petropoulou M., Giannios D., 2015, MNRAS, 450, 183
  • Springel et al. (2005) Springel V., di Matteo T., Hernquist L., 2005, MNRAS, 361, 776
  • Spruit et al. (1997) Spruit H. C., Foglizzo T., Stehle R., 1997, MNRAS, 288, 333
  • Spruit et al. (2001) Spruit H. C., Daigne F., Drenkhahn G., 2001, A&A, 369, 694
  • Tayler (1957) Tayler R. J., 1957, Proc. Phys. Soc. B, 70, 1049
  • Tchekhovskoy & Bromberg (2016) Tchekhovskoy A., Bromberg O., 2016, MNRAS, 461, L46
  • Tchekhovskoy et al. (2008) Tchekhovskoy A., McKinney J. C., Narayan R., 2008, MNRAS, 388, 551
  • Tchekhovskoy et al. (2009) Tchekhovskoy A., McKinney J. C., Narayan R., 2009, ApJ, 699, 1789
  • Tchekhovskoy et al. (2010a) Tchekhovskoy A., Narayan R., McKinney J. C., 2010a, New Astronomy, 15, 749
  • Tchekhovskoy et al. (2010b) Tchekhovskoy A., Narayan R., McKinney J. C., 2010b, ApJ, 711, 50
  • Tchekhovskoy et al. (2011) Tchekhovskoy A., Narayan R., McKinney J. C., 2011, MNRAS, 418, L79
  • Tomimatsu (1994) Tomimatsu A., 1994, PASJ, 46, 123
  • Vlahakis (2004) Vlahakis N., 2004, ApJ, 600, 324
  • Vlahakis & Königl (2003a) Vlahakis N., Königl A., 2003a, ApJ, 596, 1080
  • Vlahakis & Königl (2003b) Vlahakis N., Königl A., 2003b, ApJ, 596, 1080
  • Weinberger et al. (2018) Weinberger R., et al., 2018, MNRAS, 479, 4056

Appendix A Local adaptive time-stepping

Figure 20: The Lorentz factor along a field line for our lowest resolution models, B10-SLR (with LAT implemented; in red solid) and B10-SLRL (without LAT, see text; in black dashed dotted) at . The spikes in correspond to inversion failures and are, therefore, unphysical. Using LAT, we reduce these failures and achieve a simulation speedup by a factor .

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 . 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 (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

Figure 21: Top: A colour map of density at for fiducial model B10 overplotted by grid lines (black). As designed, the grid follows the shape of the jet. Bottom: the cell aspect ratio for two field lines, one in the inner jet and the other near the jet-edge. To resolve the jet’s microstructure in both dimensions this ratio is ideally kept below .

We have designed a grid that can track the shape of the jet over orders of magnitude in distance (Fig. 21, top). Furthermore, the grid keeps the cell aspect ratio in the outer jet below (Fig. 21, bottom), such that turbulent eddies at the disc-jet boundary remain resolved. The internal grid coordinates (, , , ) are related to the real physical coordinates (, , , ) as follows,


where , and . The parameters , , and are used to focus resolution on the jet, while is used to focus extra resolution on outer parts of the grid.

Appendix C Convergence of jet properties

Figure 22: Accretion rate , normalised magnetic flux and normalised jet power from model B10 is compared with the high resolution B10-HR model. The parameters evolve similarly which gives us confidence that the B10 model is sufficiently converged. averaged over is taken as the normalisation.
Figure 23: The specific total energy , magnetisation and the Lorentz factor along a field line close to the upper jet axis from model B10 is compared with the high resolution B10-HR model at . Additionally, the parameters for B10 at is also shown. Resolving the pinch instabilities is thus important as the jet is mass-loaded over time, vastly changing dynamics.

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 and extends till , 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 ( drops) earlier than the jet in B10, presumably because mass-loading is more efficient at higher resolutions, which capture the small scale eddies better.

Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description