11institutetext: Center for Theoretical Physics, Polish Academy of Sciences,
Al. Lotnikow 32/46, 02-668 Warsaw, Poland

# Magnetically arrested accretion disks launching structured jets in
application to GRB and AGN engines

Agnieszka Janiuk<EMAIL_ADDRESS>Bestin James 11

###### Abstract

Context. We explore the formation, energetics, and geometry of relativistic
jets along with the variability of their central engine. We study both fast
and slowly rotating black holes and address our simulations to active galaxy
centers as well as gamma ray burst engines.

Aims. The structured jets are postulated to account for emission properties of
high energy sources across the mass scale, launched from stellar mass black
holes in gamma ray bursts (GRBs) and from supermassive black holes in active
galactic nuclei (AGNs). Their active cores contain magnetized accretion disks
and the rotation of the Kerr black hole provides a mechanism for launching
relativistic jets. This process works most effectively if the mode of
accretion turns out to be magnetically arrested. In this mode, the modulation
of jets launched from the engine is related to internal instabilities in the
accretion flow that operate on smallest time and spatial scales. As these
scales are related to the light-crossing time and the black hole gravitational
radius, the universal model of jet-disk connection is expected to scale with
the mass of the black hole.

Methods. We investigated the jet-disk connection by means of 3D general
relativistic magneto-hydrodynamical simulations of the magnetically arrested
disk accretion in Kerr geometry. We also quantified the variability of the
disk by means of a Fourier analysis.

Results. We found that the system evolution is governed by the physical
parameters of the engine, such as the black hole spin and disk size, as well
as disk magnetization, and we applied our scenarios to typical types of
sources in AGN and GRB classes. We found that the magnetically arrested disk
(MAD) scenario is applicable to AGN engines and supports persistent jet
emissions. It can also be applied to GRBs, as it gives the variability pattern
roughly consistent with observations. However, in some cases, strong magnetic
fields may lead to jet quenching, and this effect is found to be important
mainly for GRB jets. We speculate that it may be related to the strength of
magnetically driven winds from the GRB engines.

###### Key Words.:

black hole physics – gamma-ray burst: general – galaxies: active – stars:
winds, outflows

##  1 Introduction

The central engine of an active galactic nucleus (AGN) is powered by accretion
flow mediated by magnetic fields and centrifugal forces, where the matter
accretes onto a rotating or non-rotating black hole (Krolik, 1999). Accretion
powers the relativistic jets of AGN, which manifest their complex structure on
pc, kpc, and Mpc scales (Harris & Krawczynski, 2006). The power of these jets,
in many radio loud galaxies, is comparable to their accretion rates,
$P_{jet}\sim\dot{M}c^{2}$ (Punsly, 2007; Sikora et al., 2013). It has been
argued that for this level of efficiency among jets produced by the Blandford-
Znajek mechanism, (Blandford & Znajek, 1977), the magnetic flux accumulated on
the black hole horizon must be extremely high (Sikora & Begelman, 2013), so
that the accretion proceeds in the magnetically dominated (MAD) mode
(Igumenshchev, 2008; Tchekhovskoy et al., 2011; McKinney et al., 2012; Narayan
et al., 2012; Tchekhovskoy & Giannios, 2015). A similar engine (on scales of
one order of magnitude smaller in size) operates in gamma ray bursts (GRBs),
where it launches baryon-free, relativistically expanding, and collimated jets
(Sari et al., 1999; Rhoads, 1999).

An open question abounds regarding the main physical parameter responsible for
the jets efficiency and the driver of the AGN radio-loudness (Rusinek et al.,
2020). Previous works (e.g. Moderski et al. (1998)) have suggested that the
black hole spin is an important quantity in this context. However, more recent
works (e.g., van Velzen & Falcke (2013)) have claimed that the spin is not the
main driver and contributes only to the scatter of radio-loudness distribution
of active galaxies. As found by Avara et al. (2016) in their gamma-ray
magnetohydrodynamic (GR MHD) simulations, the jets launched from moderately
spinning black holes, $a=0.5$, are twice as efficient for the geometrically
thin disks than for thick ones. However, the radiation-transport simulations
presented by Morales Teixeira et al. (2018) have shown that the radiative
efficiency of luminous disks is in line with the Novikov-Thorne model
prediction, while it is unclear how much of the radiation and thermal energy
trapped in the outflows could ultimately escape. This effect should be
confronted with the spectral state transitions in stellar black hole binaries
and their jets being quenched in their soft states (Fender et al., 2004).

The significant part of the jet acceleration after it leaves the black hole
region is inside the hot cocoon, whose presence may alter the internal
energetics (e.g., Lazzati et al. (2017)). Such a cocoon may plausibly be
formed from an uncollimated wind outflow that was launched from the accretion
disk and driven, for instance, by magnetic fields (Murguia-Berthier et al.,
2021). On the other hand, the imprint of the MRI within the accretion disk at
the base of the jet will manifest itself in the variability timescales
(Sapountzis & Janiuk, 2019). Hence, the interplay between this variability and
the ultimate jet Lorentz factor should reveal the properties of the magnetic
field configuration preserved in the accretion torus. From an observational
perspective, this variability results from the periodic structures seen in the
radio brightness of some large-scale jets (Godfrey et al., 2012) or in the
minute-timescales variability seen in mid-infrared in black hole X-ray
binaries (Baglio et al., 2018). Central engine properties can be also
imprinted on the observed variability of GRB jets, including correlations
between the peak energy and power density spectral slope (Dichiara et al.,
2016).

The maximum jet power in the MAD state was estimated to be on the order of
betwen three and four times the accretion power, that is,
$P_{MAD}=1.3h_{0.3}a^{2}{\dot{M}}c^{2}$ (Davis & Tchekhovskoy, 2020). In terms
of the magnetic pressure driven onto the black hole, which is proportional to
the square of the magnetic flux, $\Phi^{2}_{BH}$, the balance of the ram
pressure of matter, given by the accretion rate, should be satisfied. The
relative strength of the ram and magnetic pressure scales, thus, the disk-jet
connection, while the dimensionless quantity,
$\phi_{BH}=\Phi_{BH}/{\sqrt{\dot{M}r_{g}c^{2}}}$, defined at the black hole
horizon, gives the characteristic threshold for magnetically arrested state.
In this state, the dimensionless magnetic flux $\phi_{BH}$ can be in the range
between 10-20, which means that it can be stronger than the gravitational
force pulling the matter to the black hole and can actively push it back from
the horizon. While the material gets accreted, the magnetic flux cannot escape
due to the inward pressure. The process is regulated via interchange
instabilities and magnetic field reconnections. The detailed picture of this
MAD state is described by strong non-linear interactions between magnetic
field and the gas, and comprises the subject of complex numerical simulations.

Gamma ray bursts are the transient sources resulting from a collapse of
massive stars, or compact binary mergers, where the relativistic jets launched
from the central engine are responsible for the observed emission (see reviews
in e.g., Lee & Ramirez-Ruiz (2007); Corsi & Lazzati (2021)). The prompt
emission in gamma rays is modulated on very short timescales, which can
reflect the variable activity of the central engine (Kobayashi et al., 1997).
In the MAD state, this variability is governed by the free-fall timescale of
matter in the region where it is stalled by magnetic fields. The parametric
model proposed by Lloyd-Ronning et al. (2018) relates the degree of
arrestedness (the so-called ”MAD-ness parameter”), the size of arrested
region, and the mass of black hole with the observable variability timescales
– which are on the order of $\sim 1$ second or less, for typical parameters.
Observed correlations between the minimum timescale of variability and the jet
Lorentz factor (Wu et al., 2016), and between the computed power-density slope
and the black hole spin (Janiuk et al., 2021) seem to confirm the crucial role
of the central engine in governing this variability.

In this work, we perform a set of numerical simulations of the central engine
as composed of the rotating black hole surrounded by accretion torus. Our
models can be applicable to both AGN cores and GRB centers. We address here
the widely-adopted MAD paradigm, which has been aimed at explaining the power
released in the jets launched from the magnetically arrested disks and their
time variability. The novelty of our approach allows us to infer the
observationally testable results for a wide range of sources, which scale from
stellar mass black holes to the supermassive ones. In contrast to the previous
studies which were focused, for instance, on the specific aspects of the state
transitions in black-hole X-ray binaries (Dexter et al., 2014) or the
changing-look AGNs (Scepi et al., 2021), here we model the short-term
variability of the accreting black holes across their mass scale.

Our numerical setup is based on the full general relativistic MHD framework
and the simulation results are scalable with the black hole mass over many
orders of magnitude. We show fully non-axisymmetric 3D simulation results in
the fixed Kerr background metric that are sufficient for a negligible mass and
spin change of the black hole (for other treatment, see e.g., Janiuk et al.
(2018)). The torus matter is highly magnetized and its accurate evolution is
followed when the MHD turbulence (Balbus & Hawley, 1991) drives the mass
inflow. In addition, the magnetic field is responsible for launching the
uncollimated outflows of the plasma, in the form of magnetically driven winds,
while the rotation of the black hole powers the bipolar jets.

The numerical scheme is based on the GRMHD code HARM (Gammie et al., 2003;
Noble et al., 2006). The code has been extended to 3D and efficiently
parallelized with a hybrid MPI-OpenMP technique (Janiuk et al., 2018). Our
code is publicly available under the name HARM-COOL
111https://github.com/agnieszkajaniuk/HARM_COOL, as it contains the module for
neutrino cooling of dense and hot plasma (this specific module currently only
operates in 2D). The current paper is a follow-up of our previous studies
(Sapountzis & Janiuk, 2019), which presented the axisymmetric model. For
efficient 3D runs, we adopted a technique of minimizing the total entropy in
the cells where the magnetic pressure is extremely high. The initial setup of
our code is also changed due to the breaking of axial symmetry in the flow via
random perturbations of internal energy. In this way, we allow the plasma to
start leaking onto the black hole through the finger-like structures, within
the magnetically arrested regions. As a result, we are able to probe much more
highly magnetized configurations.

We study models with several values of the black hole spin parameter, broadly
ranging from fast to non-rotating black holes. In addition, we probe two
limits for the accretion disk size, as determined by the initial radius of the
pressure maximum. The study is addressed to both stellar mass black hole-jet
systems in the GRB engines and to the AGN engines.

This article is organized as follows. In Section 2, we present the simulation
scheme and the initial configuration of our model. The results showing the
general properties of the flow and the results of the jet power and efficiency
analysis are presented in Section 3. Our discussion and conclusions make up
Section 4.

##  2 The simulation setup

We use the general relativistic magneto-hydrodynamic code, HARM (Gammie et
al., 2003; Noble et al., 2006; Sapountzis & Janiuk, 2019), with a fixed
background Kerr metric, that is, we neglected the effects of self gravity and
the BH spin changes (cf. Janiuk et al. (2018)). The HARM code is a finite-
volume, shock-capturing scheme to solve the hyperbolic system of the partial
differential equations of GR MHD. The numerical scheme is based on GR MHD
equations where the energy-momentum tensor, $T^{\mu\nu}$, is contributed by
the gas and electromagnetic fields:

$\displaystyle{T_{\left(m\right)}}^{\mu\nu}=\rho\xi
u^{\mu}u^{\nu}+pg^{\mu\nu},$
$\displaystyle{T_{\left(em\right)}}^{\mu\nu}=b^{\kappa}b_{\kappa}u^{\mu}u^{\nu}+\frac{1}{2}b^{\kappa}b_{\kappa}g^{\mu\nu}-b^{\mu}b^{\nu},$
$\displaystyle
T^{\mu\nu}={T_{\left(m\right)}}^{\mu\nu}+{T_{\left(em\right)}}^{\mu\nu}.$

Here $u^{\mu}$ is the four-velocity of gas, $u$ denotes internal energy
density, $b^{\mu}$ is the magnetic four-vector and $\xi$ is the fluid specific
enthalpy. The continuity and momentum conservation equation is expressed as:

$\ (\rho u^{\mu})_{;\mu}=0,\hskip 28.45274ptT^{\mu}_{\nu;\mu}=0.$ (1)

These are brought into a conservative form by implementing a Harten, Lax, and
van Leer (HLL) solver to numerically calculate the corresponding fluxes. In
terms of the Boyer-Lindquist coordinates, $\left(r,\theta,\phi\right)$, the
black hole is located at $0<r\leq r_{\rm h}$, where $r_{\rm
h}=\left(1+\sqrt{1-a^{2}}\right)r_{\rm g}$ is the horizon radius of a rotating
black hole with mass, $M,$ and angular momentum, $J,$ in geometrized units,
$r_{\rm g}=GM/c^{2}$, and $a$ is the dimensionless Kerr parameter,
$a=J/(Mc),0\leq a\leq 1$. In our simulations, we adopt the rotating black
hole, with $a=0.9$.

The HARM code doesn’t perform the integration in the Boyer-Lindquist
coordinates, but instead in the so called Modified Kerr-Schild ones:
$\left(t,x^{(1)},x^{(2)},\phi\right)$ (Noble et al., 2006). The transformation
between the coordinate systems is given by:

$\displaystyle r=R_{0}+\exp\left[{x^{(1)}}\right]$
$\displaystyle\theta=\frac{\pi}{2}\left(1+x^{(2)}\right)+\frac{1-h}{2}\sin\left[\pi\left(1+x^{(2)}\right)\right]$

where $R_{0}$ is the innermost radial distance of the grid, $0\leq x^{(2)}\leq
1$, and $h$ is a parameter that determines the concentration of points at the
mid-plane. In our models, we use $h=0.3$ (noting that for $h=1$ and a uniform
grid on $x^{(2)}$, we obtain an equally spaced grid on $\theta$, while for
$h=1$ the points concentrate on the mid plane). The exponential resolution in
the $r$-direction leads to higher resolution and it is adjusted to resolve the
initial propagation of the outflow.

In the current simulation, we use additional modifications of the original
grid, to optimally resolve the regions closest to the black hole horizon and
the polar axis. The cylindrified coordinates are used inside the radius
$R_{\rm cyl}=3.5$ and inside the angle $\theta_{\rm cyl}=-1+1/N2$. This
restricts grid reshaping to the one single grid cell closest to the pole and
minimizes the resulting change in resolution at the poles.

Our computational domain is spanning $10^{5}$ gravitational radii, while the
grid cells are scaled with distance from the center. The inner and outer
radius are set to $R_{\rm in}=0.87\times(1.+\sqrt{(}1.-a^{2}))$, and $R_{\rm
out}=10^{5}$. The progressively sparser grid starts at $R_{\rm br}=400$, with
$x^{(1)}_{i}=x^{(1)}_{i-1}+1/C\times\delta x$, where $\delta
x=x^{(1)}_{i}-(\log(R_{br}-R_{0}))^{4}$. We use $C=1$ to limit the hyper-
exponentiation. Our grid resolution is $N1\times N2\times N3=288\times
256\times 256$ in the $r,\theta$, and $\phi$ directions.

The boundary conditions implemented in the code are free outflow boundary in
the radial direction, transmitting and reflecting in the polar direction, and
periodic in the azimuthal direction. In the radial direction, primitive
variables (density, internal energy, and radial magnetic field component), are
copied from the boundary zones to the two ghost zones. The radial velocity and
the polar and azimuthal components of the magnetic field are also
extrapolated. In the polar direction, the primitive variables are copied to
the ghost zones and we have two such zones, located on the opposite side of
the polar axis. The $\theta$-components of velocity and magnetic field change
signs at the polar axis and, in addition, the transverse velocity at the polar
axis is linearly interpolated over two innermost cells in $\theta$.

###  2.1 The torus initial configuration

The accreting material is modeled following Fishbone & Moncrief (1976)
(hereafter, FM) based on the analytic solution of a constant specific angular
momentum, in a steady-state configuration of a pressure-supported ideal fluid
in the Kerr black hole potential. The initial configuration of our fiducial
model is shown in Figure 1.

Figure 1: Initial state of the model GRB-HS100, in the polar and equatorial
cuts, as shown in the top and bottom panels, respectively. Maps present the
density distribution in the FM torus configuration and poloidal magnetic
field.

The location of the material reservoir is determined by the radial distance of
the innermost cusp of the FM torus, $r_{\rm in}$, and the distance where the
maximum pressure occurs, $r_{\rm max}$. In the current simulations, we adopt
values 6 (12) and 13 (25), for $r_{\rm in}$ and $r_{\rm max}$, respectively.
The relative difference of the two radii determines also the dimension of the
torus. Subsequently, we obtain the angular momentum value and the distribution
of the angular velocity along the torus. For our parameters, and back hole
spin of $a=0.9$, it is equal to $l_{\rm spec}=4.397$, and $l_{spec}=5.556$, in
a smaller and larger torus, respectively. We also note that the torus size and
the value of specific angular momentum in the FM solution depend on the Kerr
parameter. For small black hole spin, the torus size decreases, while the
specific angular momentum is increased.

The initial torus is embedded in a poloidal magnetic field, prescribed with
the vector potential,

$A_{\varphi}=\max(({\bar{\rho}\over\rho_{\rm max}}-\rho_{0})\times r^{5},0),$
(2)

where $\bar{\rho}$ is the density in the torus (averaged over 2 neighboring
cells) and $\rho_{\rm max}$ is the density maximum, and we use $\rho_{0}=0.2$,
to restrict the initial field to the regions where density in the torus is
larger than $0.2\rho_{\max}$. The factor of $r^{5}$ ensures that higher
magnetic flux will be brought onto the black hole horizon from larger
distances, as the evolution proceeds. Other spatial components of the initial
vector potential are zero, $A_{r}=A_{\theta}=0$; hence, the magnetic field
vector will have only $B_{r}$ and $B_{\theta}$ components The vector potential
is related to the Faraday tensor,
$F_{\mu\nu}=\partial_{\mu}A_{\mu}-\partial_{\nu}A_{\mu}$, and we thus obtain
the magnetic field four-vectors with $b^{\mu}=-^{*}F^{\mu\nu}u_{\nu}$. It is
related to the normal observers magnetic field vector, as
$B^{\mu}=\alpha^{*}F^{\mu t}$, where $\alpha=1/g_{tt}$ is the lapse (noting
that there is no time-component for $B^{\mu}$).

The magnetized flows are, by definition, not in equilibrium and the angular
momentum is transported. The plasma $\beta-$parameter is defined as the ratio
of the fluid’s thermal to the magnetic pressure, $\beta\equiv p_{\rm g}/p_{\rm
mag}$. We normalize the magnetic field in the FM torus to have
$\beta=(\gamma-1.)*u_{max}/(0.5*b^{2}_{max})$, where $u_{max}$ is the internal
energy at the torus pressure maximum radius. We examined several models with
different minimal initial $\beta$ of the plasma (with tested values of 30, 50,
and 100).

For the thermal pressure, we adopt the adiabatic equation of state, of the
form $p=(\gamma-1)u$, and we use the adiabatic index of $\gamma=4/3$. Finally,
in order to introduce a non-axisymmetric perturbation and allow for the
generation of azimuthal modes, permitting the gas to overpass the magnetic
barrier, we impose an initial perturbation of internal energy. It is given by
$u=u_{0}(0.95+0.1C)$, where $C$ is a random number generated in the range
(0;1). This perturbation on the order of less than 5% is typically used in
other 3D simulations of weakly magnetized flows (see Mizuta et al. (2018)).

The code works in dimensionless units of $c=G=M=1$. For the purpose of these
computations, we keep the scaling, and the black hole mass $M$ does not change
during the simulation. In physical scales, the GRB engines are modeled with
the following global parameters: mass of the black hole, $M_{\rm BH}$; the
torus mass; and black hole dimensionless spin, $a$. Typical short GRB
parameters are $M_{\rm BH}=3\,M_{\odot}$, and $a=0.6-0.9$, while the torus
mass is about 0.1 $M_{\odot}$. The latter depends on the configuration of our
torus (i.e., the inner and pressure maximum radii) and can be scaled to
physical units by assuming a density normalization. For the GRB case, we adopt
the density scaling to keep the mass of torus enclosed in the volume, so that
$D_{unit}=M_{unit}/L_{unit}^{3}=3.4\times 10^{11}$ g cm-3. The magnetic flux
is then normalized with $\Phi_{unit}=\sqrt{4\pi D_{unit}c^{2}}L_{unit}^{2}$
(see Table 1).

Source | $M_{BH}$ | $M_{disk}^{unit}$ | Timeunit | $\dot{M}^{unit}$ | $\dot{E}^{unit}$ | $D^{unit}$ | Bunit | $\Phi_{B}^{unit}$
---|---|---|---|---|---|---|---|---
| ($M_{\odot}$) | ($M_{\odot}$) |  |  | (erg s-1) | (g cm-3) | (G) |
Short GRB | 3 | 1.5$\cdot 10^{-5}$ | 1.5$\times 10^{-5}$ s | 1 $M_{\odot}$s-1 | 1.8 $\cdot 10^{54}$ | 3.4 $\cdot 10^{11}$ | 6.2 $\cdot 10^{16}$ |  $1.2\cdot 10^{28}$ G cm2
Long GRB | 10 | 1.5$\cdot 10^{-4}$ | 4.9$\times 10^{-5}$ s | 3 $M_{\odot}$s-1 | 5.5 $\cdot 10^{54}$ | 9.5 $\cdot 10^{10}$ | 3.2 $\cdot 10^{16}$ |  $7.0\cdot 10^{28}$ G cm2
Sgr A$\star$ | $4\cdot 10^{6}$ | 6.2$\cdot 10^{-10}$ | 6.2$\times 10^{-7}$ yr | 10-3 $M_{\odot}$yr-1 | 5.6 $\cdot 10^{43}$ | 6.0 $\cdot 10^{-12}$ | 2.6 $\cdot 10^{5}$ |  $9.0\cdot 10^{-8}$ G pc2
M 87 | $5\cdot 10^{9}$ | 7.8$\cdot 10^{-7}$ | 7.8$\times 10^{-4}$ yr | 10-3 $M_{\odot}$yr-1 | 5.6 $\cdot 10^{43}$ | 3.8 $\cdot 10^{-18}$ | 2.6 $\cdot 10^{2}$ |  $1.2\cdot 10^{-5}$ G pc2
Table 1: Conversion units for various astrophysical sources

In case of AGN, we can apply our models to low-luminosity AGN, such as Galaxy
center SgrA$\star$ or M87. In these sources, the accretion rates are very
small and our radiatively inefficient models are still a viable approximation.
We therefore adopted two similar scalings for accretion rates, on the order of
$10^{-5}$ $M_{\odot}$yr-1, for these AGN (Quataert et al., 1999; Cuadra et
al., 2005; Prieto et al., 2016). We adopted the black hole mass of $4\times
10^{6}M_{\odot}$ and $5\times 10^{9}M_{\odot}$, for SgrA$\star$ and M87,
respectively, and we derive our length units, according to their gravitational
radius. We apply the scaling of SgrA$\star$ to the simulation with a non-
spinning black hole (Melia et al., 2001; Broderick et al., 2011), while the
model with black hole spin of $a=0.9$ represents the M87 case (Tamburini et
al., 2020). The magnetic flux on the black hole horizon is expected to be on
the order of $10^{-10}-0.01$ G pc2 (McKinney et al., 2012).

Density scaling in the accretion disk of AGNs is more arbitrary than for GRBs,
as the spacial extension of the disk is much larger than in case of compact
binaries and the density drops by many orders of magnitude, from the black
hole up to the broad line region (Czerny et al., 2016). At the outskirts of
the disk, it can be as low as the insterstellar medium (ISM) density,
$\rho_{\rm ISM}\sim 10^{-27}$ g cm-3. Also, the disk thickness in AGN is much
smaller than its radius, while the self-gravity effects may further reduce the
surface density. In all our simulations, the total mass of the disk enclosed
in our computational domain, $R_{out}=10^{5}r_{g}$, will be a fraction of the
solar mass. Magnetic flux is again normalized to physical units with the
scaling of density, where now we adopt $D_{unit}=6.03\times 10^{-12}$ g cm-3
and $\Phi_{unit}=9.08\times 10^{28}$ G cm2 and then $D_{unit}=8.85\times
10^{-18}$ g cm-3 and $\Phi_{unit}=1.15\times 10^{32}$ G cm2, for SgrA$\star$
and M87, respectively.

##  3 Results

The bipolar Poynting-dominated jets emerge in our GRB and most of AGN models.
They are collimated close to the black hole by magnetic stress, created by
rotation of the disk and associated with a disk wind (Globus & Levinson,
2016), as well as supported by toroidal magnetic field (Levinson & Globus,
2017). Figure 2 shows an example of such a structure, where the central engine
is at the base of the jet.

Figure 2: 3D view of the inner parts of torus in the end of simulation with $\beta=100$ and black hole spin of $a=0.9$ (model GRB-HS100). The color scale presents density of the flow, while the thick blue line is a visualization of the magnetic field structure within the jet, with a chosen field line originating in the region close to black hole horizon (sphere of radius $r=2$). We note that the plot does not show the full computational domain, but it is zoomed to the innermost $\sim 100$ $r_{g}$. For a GRB model with a black hole of three solar masses in the engine, the physical size in this plot is about $4.5\times 10^{7}$ cm.  Model |  $M_{torus}$ (t=0) |  $M_{torus}$ (t=$t_{end}$) | BH spin | $r_{\rm in}$ | $r_{\rm max}$ | $\beta_{0}$ | $P_{BZ}$ | $\langle\dot{M}_{in}\rangle$ | $\langle\Phi_{B}\rangle$ | $\Delta M_{out}$
---|---|---|---|---|---|---|---|---|---|---
GRB-HS100 | 26443 | 12500 | 0.9 | 6.0 | 13.0 | 100 | 0.0102 | 0.147 | 2.121 | 6.464
GRB-HS50 | 26466 | 12189 | 0.9 | 6.0 | 13.0 | 50 | 0.0013 | 0.143 | 1.338 | 12.89
GRB-HS30 | 26498 | 15634 | 0.9 | 6.0 | 13.0 | 30 | - | 0.269 | 1.317 | 21.32
GRB-LS100 | 6562 | 2443 | 0.3 | 6.0 | 13.0 | 100 | 0.0319 | 0.082 | 2.191 | 0.0004
GRB-LS50 | 6562 | 2392 | 0.3 | 6.0 | 13.0 | 50 | 0.0228 | 0.075 | 2.359 | 0.0007
GRB-LS30 | 6562 | 2534 | 0.3 | 6.0 | 13.0 | 30 | 0.0105 | 0.065 | 2.284 | 0.0013
AGN-HS100 | 532647 | 309765 | 0.9 | 12.0 | 25.0 | 100 | 0.120 | 3.765 | 8.801 | $8.3\times 10^{-5}$
AGN-HS50 | 532647 | 274132 | 0.9 | 12.0 | 25.0 | 50 | 0.109 | 3.465 | 9.177 | $2.5\times 10^{-3}$
AGN-HS30 | 532647 | 293182 | 0.9 | 12.0 | 25.0 | 30 | 0.069 | 3.104 | 8.129 | $2.9\times 10^{-3}$
AGN-NS100 | 298344 | 195455 | 0.0 | 12.0 | 25.0 | 100 | – | 2.284 | 8.421 | $4.9\times 10^{-5}$
AGN-NS50 | 298344 | 210646 | 0.0 | 12.0 | 25.0 | 50 | – | 1.976 | 6.442 | $6.1\times 10^{-5}$
AGN-NS30 | 298344 | 210083 | 0.0 | 12.0 | 25.0 | 30 | – | 2.043 | 6.197 | $7.2\times 10^{-5}$
Table 2:  Summary of the models. The inner radius of the torus, $r_{\rm in}$,
the radius of pressure maximum, $r_{\rm max}$, and the plasma-$\beta$ are
given as the initial state parameters. The last two columns give the magnetic
flux averaged over the simulation time, and the total mass lost through the
outer boundary. All quantities are given in code units. The BZ power is given
at the mid-time of simulation, at $[25000t_{g}]$. The final time for all
simulations is 50000 $t_{g}$

###  3.1 General structure of the accretion flow

Our simulations are divided into two sets with the GRB class of models
referring to a system composed of a stellar mass black hole surrounded by a
remnant accretion disk formed following a compact binary merger, while the AGN
class of models refer to an engine of active galaxy, where a supermassive
black hole accretes matter of a diluted torus in the circum-nuclear region.
The second class of models adopts larger spatial extension of the tori, in
terms of their pressure maximum radius. The inner radius of the torus is also
located further out from the black hole horizon in AGN models (in terms of the
dimensionless units). We further consider fast and slow spinning black holes,
including non-spinning ones for an AGN case to describe radio-quiet sources.
Finally, we also vary the initial plasma $\beta$-parameter.

Table 2 presents a summary of the models and their initial parameters, the
total mass of the disk at initial state, the averaged mass accretion rate, the
average magnetic flux of the black hole horizon, and the total mass lost in
the outflows. We also calculate the jet power at the mid-state of the
simulations.

The overall evolution pattern is as follows. The FM torus initial state
retains its structure for a relatively long integration time. Its innermost,
densest part is enclosed with $\sim 50r_{g}$ and forms a narrow cusp through
which the matter sinks under the black hole horizon. The surface layers of the
torus, which extend to about 200 $r_{g}$, form a kind of corona. This region
is the base of sub-relativistic outflows that are launched from the center as
soon as the initial condition of the simulation is relaxed (which is at about
2000 M). The models were evolved until the time of $t_{\rm f}=50,000$ $t_{g}$
(geometrical units, $t_{g}=GM_{\rm BH}/c^{3}$). This corresponds to 0.74
seconds in the GRB case and to about 9.4 months in our AGN case.

In Figure 3, we show the profiles of density, and magnetic field streamlines,
in the $r-\theta$ polar plane (XZ), and in the equatorial $r-\phi$ plane (XY),
as obtained at the mid-time of the simulation. We show three different times
during the simulation, $t=20,000M$ $t=25,000M$, and $t=30,000M$, in the left,
middle, and right panels, respectively. This sequence of snapshots shows that
the dense torus wobbles about the equatorial plane, while the jet funnel is
slightly asymmetric with respect to the vertical axis. The model shown here is
GRB-HS100, representing a weakly magnetized, small torus with $\beta=100$ and
black hole spin of $a=0.9$. The top panels are scaled to distance of 100
$r_{g}$, while the middle and bottom panels show a zoom-in on the accretion
region, within the inner 10 gravitational radii. As shown in these maps, the
magnetically arrested flow appears in this simulation, forcing the accretion
to proceed via the azimuthally distributed ”fingers” surrounding the black
hole and entering its horizon. The zoomed-in equatorial plots show that the
magnetic field lines are pulled into black hole within the low-density
regions, while the high-density fingers are pushing out the magnetic fields.
Thus, the magnetic tension can locally push the gas out and suppress the
accretion rate, while in denser regions, the matter prevents the flux from
entering the horizon. We checked that the total mass accretion rate integrated
over the horizon region in these three time intervals is equal to 0.0208,
0.0224, and 0.0110, respectively, while the magnetic flux integrated over the
horizon is 1.11, 1.53, and 1.35 (in code units). Hence, the smallest net
accretion rate and highest rate characterizing the magnetically arrested state
are found in the final time snapshot shown. It is only very high resolution
simulations (see e.g., (Ripperda et al., 2022)) that can verify whether the
MAD regions with the dense magnetic flux getting pushed out can serve as sites
of magnetic reconnection. Our simulations have moderate resolution in 3D, so
that we encompass this effect only qualitatively. We also note that when the
jet forms, hot luminous regions appear near the black hole rotation axis,
which have a very low baryon density and a vertically oriented magnetic field.
These polar regions represent bipolar jets, with a luminosity that depends on
the black hole spin (cf. Table 2).

|  |
---|---|---
|  |
|  |
Figure 3: Maps of the density isocontours in the equatorial plane XY (bottom
row), and polar plane XZ (top row). The model is $GRB-HS100$ and its
parameters are $\beta=100$ and $a=0.9$. The color maps are in logarithmic
scale and are taken at the three times of the simulation, at $t=20,000M$,
$t=25,000M$, and $t=30,000M$. Streamlines show the magnetic three-vector
configuration. We show the equatorial cuts of a small region of $10r_{g}$ from
the black hole (bottom panels) and the polar cuts are shown for small and
large regions, up to $100r_{g}$ (upper panels).

In the model of an AGN disk around a non-spinning black hole, the weakly
magnetized torus is of a larger size and its inner radius and radius of
pressure maximum are 12 and 25 $r_{g}$, respectively. We may notice that in
this model, the large mass of the torus and zero black hole spin does not
allow for the formation of magnetically dominated, low-density funnels at the
polar region. The non-spinning black holes in AGN engines do, however, also
accrete in the MAD mode. As shown recently by Mościbrodzka et al. (2021), the
specific features of this accretion mode lead to the polarized radiation
emission from the region of the black hole horizon. Here, in contrast to the
emission from remote jets, the magnitude of the circular polarization is high.
Because of high Faraday thickness, the circular polarization is particularly
sensitive to dynamics of toroidal magnetic field in the accretion flow.

The case of a fast-rotating black hole engine of a gamma ray burst, with
smaller values of initial $\beta$, equal to 50, and 30, was studied in models
GRB-HS30 and GRB-HS50. Here, also the polar funnel and equatorial disk
structure have been developed, with accretion proceeding through the equator,
while bi polar, magnetically dominated jets are launched. The qualitative
differences between these runs and $\beta=100$ case are seen in the density
range in the equatorial plane, which is (on average) larger than small $\beta$
models. The clear asymmetry in the azimuthal distribution of matter, which is
otherwise developing in the GRB-HS30 model, results from stronger magnetic
fields which are pushing the accretion flow out from the horizon. The barrier
shape is also not symmetric, so the gas must overpass it and reach the black
hole on the opposite side.

###  3.2 Mass inflow and the MAD state

Figure 4: Time dependence of the mass flux through the black hole horizon (top
panel), horizon-integrated magnetic flux (middle panel), and the so-called
MAD-ness parameter, $\phi_{BH}$, given by Eq. 5 (bottom panel). Blue lines
denote the models with initial gas-to-magnetic pressure ratio, $\beta=100$,
cyan lines present the results for $\beta=50$, and green lines denote
$\beta=30$. Black hole spin parameter in all models is $a=0.9$. The disk inner
radius is located at 6 and maximum pressure radius at 13 $r_{g}$ (GRB-type
models).

In Figure 4, we show the time-dependence of so called MAD-ness parameter,
namely, the ratio of the mass and magnetic flux through the black hole
horizon. These are given by the relations:

${\dot{m}}_{BH}=-\int{\rho u^{r}\sqrt{-g}d\theta d\phi},$ (3)
${\Phi}_{BH}=0.5\int{|{B^{r}}|\sqrt{-g}d\theta d\phi},$ (4)

and

${\phi}_{BH}={\Phi_{BH}\over\sqrt{\dot{m}_{BH}}}.$ (5)

As shown in Figure 4, accretion rate is large at the beginning of the
simulation and it is higher then for smaller $\beta$ value. At later times,
the magnetic flux that has been dragged to the black hole vicinity, growing
high enough to capture the accretion flow. Hence, the ratio of these two
fluxes remains constant and greater than $\phi_{BH}\sim 10$ for most of the
simulation after the initial condition is relaxed. The quantitative difference
in $\phi_{BH}(time)$ profile between the models with $\beta=100$ and
$\beta=50$ is very small in the case of the average value. However, the time
variability of the more magnetized flow appears to be increased; this relation
is quantified in more detail below (i.e., Sect. 3.3).

Figure 5: Time dependence of the mass flux through the black hole horizon (top
panel), the horizon integrated magnetic flux (middle panel), and the so-called
MAD-ness parameter, $\phi_{BH}$, given by Eq. 5 (bottom panel). Red lines
denote the models with initial gas-to-magnetic pressure ratio $\beta=100$,
orchid lines present results for $\beta=50$, and magenta lines denote
$\beta=30$. Black hole spin parameter in all models is $a=0.3$. The disk inner
radius is located at 6 and maximum pressure radius at 13 $r_{g}$ (GRB-type
models).

In Figure 5, we present the results of simulation with a black hole spin of
$a=0.3$ and the same small accretion disk, with $r_{in}=6$ and $r_{max}=13$
$r_{g}$, referring to the GRB-LS100, GRB-LS50, and GRB-LS30 models. The three
values for the $\beta$-parameter were chosen to differentiate between the
magnetic field strength in these models, which are plotted with red, orchid,
and magenta lines. We observe that for this small value of black hole spin,
the magnetic flux on the horizon, $\Phi_{B}$, is almost the same for all three
values of $\beta$. It is also more stable than it was in the high-spin models,
so that it has a smaller value at the beginning of the run, but does not drop
so much over time, as compared to the models with a spin of $a=0.9$. As the
time evolution for the accretion rates is rather similar in both these
simulation runs, the average profile of MAD-ness parameter remains almost
constant with time, namely, at the level of $\phi_{BH}=20$. The timescales of
variability seem to be rather similar to each other and not depend much on the
magnetization parameter (see, however, the quantitative analysis in Sect.
3.3.)

Figure 6: Time dependence of the mass flux through the black hole horizon (top
panel), the horizon integrated magnetic flux (middle panel), and the so-called
MAD-ness parameter, $\phi_{BH}$, given by Eq. 5 (bottom panel). Blue lines
denote the models with initial gas-to-magnetic pressure ratio $\beta=100$,
cyan lines present results for $\beta=50$, and green lines denote $\beta=30$.
Black hole spin parameter in all models is $a=0.9$. The disk inner radius is
located at 12 and maximum pressure radius at 25 $r_{g}$ (AGN-type models).

Finally, we present in Figure 6 a similar set of models, assuming, for the
black hole spin, $a=0.9,$ and three-disk magnetization, normalized to
$\beta=100$, 50, and 30, while embedding the larger torus. The radii of inner
disk cusp and its pressure maximum are taken as $r_{in}=12$ and
$r_{max}=25r_{g}$, respectively. Now, both the mass accretion rate and
magnetic flux on the horizon, remain large for the time of the simulation. The
MAD state is prolonged, while the MAD-ness parameter, $\phi_{BH}$, rises
smoothly after an initial relaxation period of about a 1000 tg. The
differences between the results obtained for three values of $\beta$-parameter
are now less pronounced. From the above three sets of simulations, we conclude
therefore that the effect of magnetization changing the variability timescale
depends on the black hole spin and the size and position of the torus. In our
GRB models with high values for the black hole spin, the variability pattern
is affected by the disk magnetization. In models with either small spin or a
large torus, referring to the AGN case, the magnetic fields do not seem to
alter the variability timescales as significantly.

###  3.3 Variability analysis

We quantitatively analyzed the variability of the accretion rate onto the
black hole by means of a Fourier transformation.

We calculated the power-density spectra of the time-series of GRB models,
depicted in Figures 4 and 5; and for AGN models, shown in Fig. 6, as well as
for non-spinning black hole models. In the analysis, we first skipped $\delta
t=10,000$ tg and we focused on the variability of the flow when the initial
transient phase is relaxed. Figure 7 shows the power density spectra of four
chosen models, with the smallest initial magnetization in the disk.

Table 3 gives the values of power-law fitting to the power density spectra
(PDS) of individual models. We notice that the PDS slope depends on the black
hole spin, while only weakly on the disk magnetization. For both classes of
models, GRB and AGN, the PDS slopes tend to be steeper for lower spin of the
black hole. For highly spinning black hole models, the GRB models present
flatter PDS spectra than the AGN models. We also calculated the minimum
variability timescale, estimated as the minimum peak width at its half
maximum, over the time series (here, we also skipped the first $10,000$ tg).
There is an anticorrelation between this timescale and the disk magnetization
for GRB and AGN models. Also, models with highly spinning black holes tend to
present shorter variability timescales than those with slowly or non-spinning
black holes.

Figure 7:  PDS and their power-law fits presented for four exemplary simulations. The models are: GRB-LS100, GRB-HS100, AGN-NS100, and AGN-HS100, as marked in the figures. The light curves were binned logarithmically to obtain better fits.  Model | a | Initial $\beta$ | PDS slope | MTS
---|---|---|---|---
GRB-HS100 | 0.9 | 100 | -1.78 $\pm 0.02$ | 29.15
GRB-HS50 | 0.9 | 50 | -1.67 $\pm 0.02$ | 30.95
GRB-HS30 | 0.9 | 30 | -1.97 $\pm 0.01$ | 39.18
GRB-LS100 | 0.3 | 100 | -1.64 $\pm 0.03$ | 70.06
GRB-LS50 | 0.3 | 50 | -1.80 $\pm 0.01$ | 66.64
GRB-LS30 | 0.3 | 30 | -1.66 $\pm 0.03$ | 62.79
AGN-HS100 | 0.9 | 100 | -1.68 $\pm 0.01$ | 29.95
AGN-HS50 | 0.9 | 50 | -1.68 $\pm 0.01$ | 30.78
AGN-HS30 | 0.9 | 50 | -1.70 $\pm 0.01$ | 26.64
AGN-NS100 | 0.0 | 100 | -1.80 $\pm 0.02$ | 72.84
AGN-NS50 | 0.0 | 50 | -1.66 $\pm 0.03$ | 90.63
AGN-NS30 | 0.0 | 30 | -1.87 $\pm 0.02$ | 77.39
Table 3: Results of power density spectral fitting. Model parameters are black
hole spin $a$ and $\beta=p_{gas}/p_{mag}$. The PDS slope is fitted power law
index. MTS is the minimum variability timescale, in $t_{g}$ units.

###  3.4 Jet power and energetic structure

The jets in GRB models are Poynting energy-dominated and confined by the
toroidal magnetic fields. Figure 2 shows the 3D structure of the magnetic
field, as depicted by an arbitrarily chosen field line, tangential to the
magnetic field vectors and footed at the black hole horizon. The model shown
here is GRB-HS100 with black hole spin $a=0.9$.

In Table 2, we give the values of jet power, as calculated in code units in
the mid-time of the simulation, at $25,000\leavevmode\nobreak\ t_{g}$. The
power is computed using an estimate of the Blandford-Znajek process, where the
rotating black hole is surrounded by a force-free, magnetized plasma. We use
the formula given by McKinney & Gammie (2004):

$P_{\rm BZ}=-\int_{r=r_{h}}(b^{2}u^{r}u_{t}-b^{r}b_{t})dA,$ (6)

so that we integrate the radial part of electromagnetic flux over the surface
of black hole horizon. In general, the integral over the radial energy flux,
$T^{r}_{t}$, may be subdivided into matter and electromagnetic parts. In the
above formula, we neglect the matter component, assuming a force-free
approximation. For a steady, uniform outflow, this expression gives the radial
Poynting energy flux measured by a stationary observer at large distance from
the black hole. We notice that the BZ power is positive only if the energy can
be extracted to the remote jets from the black hole, so only the flux with a
positive sign will contribute to it. Otherwise, there is no energy extraction
from the black hole via the BZ process. These jets, if produced, will be
dominated by thermal energy.

As listed in Table 2, at the mid-time of the simulation, a strong, Poynting-
dominated bipolar jets have already developed in most of the models. We note
that for the models with zero black hole spin, such jets have not emerged;
therefore, the value computed from Eq. 6 is set to zero. We also notice that
in the model GRB-HS30, with a high black hole spin, the Poynting-dominated jet
is quenched. The polar funnels are having a low density, however, the MAD
state is not operating, according to a standard definition, with $\phi_{BH}$
value being smaller than 10.

For all other simulations, the jet power obtained in the Blandford-Znajek
process is non-zero for the black hole spin, $a>0,$ and at the mid-term of the
simulation, it ranges from $\sim 10^{42}$ to $\sim 10^{52}$ erg s-1, for AGN
and GRB, respectively. We calculate this power and we normalize it with the
adopted physical scalings, related to the black hole mass. The characteristic
quantities in physical units and unit conversions, are summarized in Table 1.
We note that the jet power in GRB scalings is very large and it could be
relevant only for some extremely bright sources. It is thus possible that the
MAD scenario is not capable of explaining faint GRBs.

In Figure 8, we show the jet power as a function of time for GRB and AGN
models. We chose the simulations with black hole spin of $a=0.9$ in both plots
and we show the power (in code units) for various degrees of disk
magnetization. We find that strong powerful jets are emitted in all AGN models
with high black hole spin and that the jet power does not depend very much on
the initial disk magnetization. For smaller $\beta$, the jet power evolves
slightly faster at the beginning of the simulation to reach its equilibrium
value. The power does not drop much until the end of the runs, that is, until
time 50,000 $t_{g}$. In our GRB models, the initial power of the jet is
smaller (in code units) than in the AGN case and it also drops by at least
three to four orders of magnitude by the end of the runs. For smallest initial
$\beta$-parameter, the jet is quenched at the mid-time of simulation and its
power drops to zero. The same occurs with the jet emitted in $\beta=50$ case,
albeit at the end of the run. This behavior is correlated with the decrease of
the dimensionless magnetic flux, $\phi_{BH}$ (cf. Fig. 4).

Figure 8: Power of the Blandford-Znajek process, calculated from Eq. 6, as a
function of time. We present models GRB-HS and AGN-HS, shown in the top and
bottom panels, respectively. Three colors denote different initial disk
$\beta$-parameters, shown with blue, red, and green points for $\beta=100$,
50, and 30, respectively.

The jet is driven by toroidal magnetic fields that operate in the accretion
disks. As shown by Mizuta et al. (2018), the Alvfenic pulses in the jets are
driven by the fluctuating toroidal field in the disk equatorial plane. Their
toroidal component dominates over the poloidal field as a result of the disk’s
differential rotation. Here, we plot the ratio of the toroidal to poloidal
magnetic field components for two simulations, with both high spin and zero
spin.

Figure 9: Ratio of the toroidal to poloidal magnetic fields,
${B_{\phi}/B_{z}}$, for models AGN-NS50 and AGN-HS50, shown in the top and
bottom panels, respectively. Color maps are in logarithmic scale

The maps in Fig. 9 are taken at the mid-time of the simulations, while both
components of the magnetic field fluctuate over time. We checked that (on
average), the toroidal-to-poloidal magnetic field ratio in the jet is larger
by two orders of magnitude in the model with high spin of the black hole than
for a non-spinning one. This toroidal field is enhanced by twisting due to the
rotation of the central black hole. In addition, the $B_{\phi}/B_{z}$ ratio in
this model is high at the disk equator up to the distance of $\sim 10r_{g}$,
because of the energy transfer from the rotating black hole to the torus via
the magnetic coupling (Janiuk & Yuan, 2010).

In Figure 10 we plot the jets energetics structure, according to the
definition:

$\mu=-{T^{r}_{t}\over\rho u^{r}},$ (7)

which is equivalent to the total (Poynting and thermal) energy of the jet,
which can be transformed to the bulk kinetic content. As can be seen from the
plot, the jets are not axisymmetric, due to the central engine structure, and
not uniform. The blobs of energy expanding with time lead to the jets
variability (cf. Janiuk et al. (2021)). The plots compare shape of the jet
base for two GRB models, with $\beta=100$ and $\beta=50$. We checked the
trends also for $\beta=30$ and models with smaller black hole spins. The main
relation between $\beta$ and the jet opening angle is that for larger gas-to-
magnetic pressure ratio, the jets seem to be more symmetric, while for smaller
$\beta,$ the asymmetric shape makes the jet narrower at one viewing angle,
while becoming broader at the opposite. The maximum energetics parameter
depends mostly on black hole spin. The azimuthally averaged jet profiles might
appear more uniform (Kathirgamaraju et al., 2019).

Figure 10: Profiles of jet energetics, defined by $\mu$-parameter (see Eq. 7)
for two GRB models with high spin, $a=0.9$, and $\beta=100$ or $\beta=50$, on
the two top, and the two bottom panels, respectively. The maps show profiles
taken at two opposite azimuthal angles, $\Phi=0$ and $\Phi=\pi$. All snapshots
are taken at a time of $t=2.5\times 10^{4}t_{g}$

####  3.4.1 Mass outflow

We also checked the mass outflow through the outer boundary. It is very small
at the beginning of the simulation, when the initial condition of the pressure
equilibrium torus is being relaxed. However, after a time of about 20000
$t_{g}$, the outflow rate value is larger and some mass is left in the system.
At the end of the simulation, the outflow rate is again small, below
$10^{-3}-10^{-7}$ (in code units). We noticed a general trend that both higher
black hole spin and higher magnetic pressure result in more powerful outflows
and the maximal mass loss rate during the simulation scales with these
parameters.

We found that the mass loss rate is somewhat larger for the GRB models, while
it is always negligibly small for our AGN models. The latter result seems in
contradiction to the large mass outflow rates reported by McKinney et al.
(2014). In their simulations, the radiatively driven outflows were more
powerful, however, the difference might also be attributed to a different
value of the outer radius in the simulations (our outer radius is five times
larger). The total mass lost during the simulation via the unbound wind
outflow in our models is given in Table 2.

####  3.4.2 Resolution of the simulations

The resolution of our simulations is moderate and would not be sufficient for
the detailed studies of magnetic field reconnections. However, we checked that
the number of cells is adequate to resolve the magnetorotational instability
(MRI). In Figure 11, we plot the quantity $Q_{MRI}$, defined as the number of
grid cells in polar direction, per wavelength of the fastest growing mode (cf.
Siegel & Metzger (2018)):

$\lambda_{\rm MRI}={2\pi\over\Omega}{b\over\sqrt{4\pi\rho h+b^{2}}}.$ (8)

The MRI is well resolved, when $Q_{\rm MRI}={\lambda_{\rm
MRI}\over{\Delta\theta}}$ is larger than 10. This requirement is satisfied in
our initial conditions, even in models with a weak magnetic field, as shown in
the figure.

|  |
---|---|---
|  |
Figure 11: Resolution of MRI instability, in the initial conditions and
evolved state. The models are $GRB-HS100$ (top row) and $AGN-HS100$ (bottom
row). The color maps are in logarithmic scale and are taken at $t=0$ (left)
and at the intermediate time of the simulation, at $t\sim 25,000M$ (right).

We also note that the results of the simulation might be somewhat sensitive to
the density and energy floors adopted to set their minimum values handled by
the numerical scheme. The floor values in our calculations were set to
$10^{-8}$ and $10^{-10}$ for the minimum density and internal energy (in code
units), respectively. These values, however, mainly influence the simulation
results at the initial time, when the density contrasts are high. As shown by
Sapountzis & Janiuk (2019), the time-averaged distributions, for instance,
those of the jet energetics, are saturated at the same level for different
floor values.

##  4 Discussion

We studied fully 3D general relativistic numerical simulations of accretion
flows that lead to a magnetically arrested (MAD) state and that can describe
the active centers of quasars or engines of gamma ray bursts. Our simulations
adopt an adiabatic equation of state that neglects the microphysics effects
and that are performed in a dimensionless setup. Hence, with the scaling of
the black hole mass and the torus density, this universal setup can be used
for both supermassive and stellar mass accreting black holes. In particular,
the timescales of flares and the time variability of energy injection to the
jets, are scaled with geometric time, $t_{g}=GM_{BH}c^{-3}$. Therefore, they
can be directly compared to both supermassive and stellar mass accreting black
holes, in AGN and GRBs. The timescales corresponding to the flaring activity
of blazars, such as that observed in Mrk 421 (Chatterjee et al., 2021) or in
Mrk 501 (Giommi et al., 2021), are less than 1 day. Similarly, the variability
timescales in Sgr A$\star$ can reveal about 1 flare per day, as reported by
(Neilsen et al., 2013) for the Chandra observations, or up to 2.5 flares per
day on average, found in Chandra and XMM-Newton data, as reported by Ponti et
al. (2015). The observed variability of gamma ray emission in GRBs down to
millisecond timescales is correlated with the burst duration (MacLachlan et
al., 2013). It has also been found that a correlation exists between the burst
variability and its intrinsic luminosity, suggesting a Cepheid-like relation
for these sources, for a sample observed mainly by BATSE (Reichart et al.,
2001). The absolute variability timescales are related to the black hole mass
and this an approach was already taken in our previous study (Janiuk et al.,
2021) to explain and observed correlation between the jet Lorentz factor and
minimum variability timescale across the mass scale (Wu et al., 2016).

The variability of prompt emission of GRBs is related to both variability of
the central engine itself and to the interaction of the jet with ambient
medium, such as disk wind, the post-merger ejecta in a short GRB (Murguia-
Berthier et al., 2021), or with the progenitor star envelope in a long GRB
case (Morsony et al., 2010). The latter may introduce an observed variability
even if the central engine fluctuations are very slight.

In the current simulations, we focus on the engine, and jet base variability,
as driven by magneto-rotational instability, or interchange instabilities in
MAD disks. The jet propagation inside the star is a complex non-linear
process, where the turbulence develops and energy may be dissipated in the
form of shocks. As this aspect is ignored in our simulations, it basically
represents the situation where the engine activity timescale is much longer
than the jet propagation timescale.

Our analysis shows that the accretion rate variations at the inner edge of the
disk may be responsible for the observed variability of the sources. The PDS
have slopes between 1.67-1.97 for the highly spinning black hole engine in the
GRB models, while the PDS slopes for moderate black hole spins are in the
range 1.64-1.8. Our simulated PDS spectra are comparable the observed
variability of GRBs. The calculated PDS slopes given in Table 3 are roughly
consistent with the canonical value of power-law PDS with $\alpha=5/3$
reported by Beloborodov et al. (2000) for the bolometric light curves of long-
duration GRBs. More recent observations show that the mean value of PDS slope
found by Ukwatta et al. (2011) was equal to $\alpha=1.4\pm 0.6$, while the
variability studies based on Swift observations of individual GRBs presented
by Guidorzi et al. (2016) revealed a much broader range of PDS slopes, namely,
$\alpha=1.3-3.2$. These authors also found that the number of pulses in the
GRB light curve is strongly anticorrelated with the PDS slope. We note,
however, that most of the variability studies have been based on the long GRB
data, which are easier to analyze thanks to their long time series. The
outlier short GRB source reported in Dichiara et al. (2016) was found to have
a very steep PDS slope of about $\alpha=3.5$, while Dichiara et al. (2013)
found the PDS slope for the sample of 44 short GRBs to be between 1.4 and 2.5.

The variability of supermassive accreting black holes measured in the X-ray,
radio, and NIR bands shows the variations of synchrotron radiation. As found
by Witzel et al. (2018), the NIR light curve variations during the low flux
density phase of Sgr A$\star$ are best fitted by PDS with a slope of
$\alpha=1.6$. We notice that this value is in agreement with our model results
for the AGN-HS simulations, namely, the required spin of the black hole would
be $a=0.9$. If the black hole in Sgr A$\star$ is a non-spinning one, then the
only simulation consistent with the observed variability of this source is the
AGN-NS30 model, namely, the one with largest magnetization.

The minimum variability timescale inferred from our models is between 25-30
and 70-90 gravitational timescales, depending on the black hole spin. For the
supermassive black hole in M87, the timescale of 30 $t_{g}$ is equal to about
8.5 days (cf. Tab. 1). Hence, our AGN models with highly spinning black hole
may account for the flaring activity of this source, as observed in the very
high energies (Aharonian et al., 2003; Acciari et al., 2009).

The jet may interact with the disk winds, which are expected to accompany the
GRB prompt phase. Such interaction can occur only when the wind mass loss rate
is large. The total mass loss rate in our models is given in Table 2. We note
that this is only a lower limit, as the simulation duration is not long enough
to account for the mass loss calculation through the outer boundary.
Nevertheless, we show that there should be a correlation between mass loss
rate and disk magnetization as well as with the black hole spin. Therefore, we
expect that the jet can be partially or completely choked, in the most
magnetized scenario. The results displayed in Figure 8 confirm this effect. In
the models with $\beta=30$, the wind density is high, while the jet power is
smaller than for the higher $\beta$ parameter. The jets launched from a
moderately spinning black hole, however, seem to have more chance to
successfully break through the wind, than those launched from highly spinning
black holes. This is because the wind density is small. The interaction with
disk winds should impose an additional variability in the observed emission,
as the outflow propagation within a denser medium will produce irregularities
of the Rayleigh-Taylor instability. As shown by Lazzati et al. (2021), the
jet’s centroid oscillates around the axis of the system, due to
inhomogeneities encountered in the propagation, and the jet slows down on a
very short timescale.

In the current 3D simulation, we use the simplified microphysics treatment in
the case of GRB disks. This is done in part to allow for a unified modeling of
both GRB and AGN types of sources, but also partially reflects our
computational limitations. The equation of state of dense matter and neutrino
cooling effects have been implemented into our numerical scheme, (Janiuk,
2019); however, because of the computational resource limitations, they
effectively work in 2D geometry. The full 3D numerical simulation of the
magnetically arrested disk in the GRB engine has been postponed to a future
work. At this stage, we speculate that the quantitative differences in the
results will be important mostly for the physics of the unbound outflows,
namely, winds that are affected by both neutrinos and magnetically driven
acceleration. In this way, the winds may be more dense and powerful if the
neutrino-driven mechanism supports or substitutes the magnetically driven wind
(Murguia-Berthier et al., 2021).

In the initial setup, we adopted the poloidal configuration of magnetic
fields, with the field lines that follow the iso-contours of constant density.
Our choice is similar to that of other authors. Initially, McKinney et al.
(2012) noticed that initially toroidally-dominated models lead to only patches
of magnetic flux and transient outflows or weak jets (see also Beckwith et al.
(2008); Paschalidis et al. (2015)). More recent works have shown that jets
form even from purely toroidal magnetic field initial conditions (Christie et
al., 2019; Liska et al., 2020). Due to the dynamo mechanism, the toroidal
component of magnetic field develops over time during the simulation and
exceeds the vertical field component. The specific choice of the magnetic
field configuration also assumes that the magnetic field is additionally
scaled with the distance to allow for the magnetic pressure to build up as the
simulation continues.

The assumption of the same magnetic fields configuration for both AGN and GRB
models is motivated by some subtle arguments. In principle, there is no
obvious reason to believe that magnetic fields geometries should differ
between AGN and stellar mass-accreting black holes, while both types of
systems are able to launch relativistic jets on the cost of black hole
rotation, mediated by magnetic fields (Davis & Tchekhovskoy, 2020). On the
other hand, the GRBs jets are systematically stronger, while the magnetic
fields amplification may originate not only from the effective accretion of a
large-scale magnetic field to the inner regions, but may also be subject to
the dynamo process acting during the stellar collapse or neutron star merger.
Recent numerical simulations performed in this context have not been
conclusive, however, about the type of configuration that should be favored
(Ponce et al., 2014; Ruiz et al., 2020; Shibata et al., 2021).

Our simulations are addressed to both AGN and GRB engines, while the
calculations are limited to the innermost regions of the accretion flow and
spatial scales close to the base of the jets. Models were initiated from a
hydrostatic torus solution and evolution was triggered with a random
perturbation of the thermal pressure, which allows for the formation of non-
axisymmetric modes of Rayleigh-Taylor instabilities. These are important for
allowing the material to overpass the magnetic barrier and continue accretion
when large magnetic flux is accumulated on the black hole horizon. The
magnetic fields are mainly responsible for the transport of the angular
momentum, which enables accretion, but also drives the disk outflows. The
time-averaged radial magnetic flux on the black hole horizon, measured in our
simulations, is on the order of $(1.5-2.5)\cdot 10^{28}$ up to $1.6\cdot
10^{29}$ G cm2 for the GRB models, depending on the adopted scaling (cf. Table
1). This range is consistent with the amount of oriented magnetic flux in
supernova progenitors (Heger et al., 2005; Obergaulinger & Aloy, 2022) and
with the field strengths in coalescing binaries of neutron stars (Kuan et al.,
2021). Magnetic flux scaling for the supermassive black hole in our Galaxy
center gives the range of $5-7.5\cdot 10^{-7}$G pc2, which can be supported by
both the interstellar medium and massive magnetized O-type stars feeding the
galactic nucleus (cf. Witzel et al. (2021); Hubrig et al. (2020)).

The initial condition adopted for the density structure, namely, the
equilibrium torus, affects the simulations at early times, while it needs to
be relaxed via the action of magnetic fields. We overcome the drawback of such
an artificial initial condition by running a simulation for a sufficiently
long time, then we follow the time evolution of the system for long time, up
to $50,000$ $t_{g}$. We note that an alternative approach would be to use an
initial setup resulting from an evolved postmerger state provided by previous
numerical relativity simulations (Hossein Nouri et al., 2018). After the
formation of the large-scale open field lines, the Blandford-Znajek process
extracts the rotational energy of the BH. This process is common to both AGN
and GRB engines, while the power of jets obtained for these two types of
sources differs. The magnetic flux available for accretion in GRBs is on the
order of $10^{-9}$ G pc2, and for active galaxies, it is a few orders of
magnitude larger. Nevertheless, the density of accreting plasma is much larger
in GRB engines, therefore, due to the compactness of the system the total
power available for a jet is larger in these objects. In addition, the
neutrino annihilation can be a complementary process to power the GRB jets
(Mochkovitch et al., 1993; Aloy et al., 2005; Janiuk et al., 2013; Liu et al.,
2015; Janiuk, 2017), but this feature should not change our main results. It
is possible, however, that the GRB jets are choked in dense wind if it is
driven both magnetically and by neutrino heating. Such a simulation is
currently beyond the scope of this paper, but it is planned in the future,
based on the formalism derived in Janiuk (2019). Recent works on this topic,
for instance, by Siegel & Metzger (2018); Fernández et al. (2019), have shown
that powerful outflows emerge from BNS merger remnant disks and can be sites
of r-process nucleosynthesis, while the bipolar jets in these simulations have
only mildly relativistic velocities, with Lorentz factors of 1-10.

Finally, we note that in our models, the adiabatic index of $\gamma=4/3$ is
used for both AGN and GRB scenarios. In fact, the ideal gas law with
$\gamma=5/3$ could be used to describe low-luminosity AGN, while the
relativistically hot gas with $\gamma=4/3$ is more likely to describe GRB
engines which are cooled by neutrino emission. In super-Eddington accretion
disks, the radiative cooling is not efficient either, as lot of the dissipated
energy is advected toward the black hole (Kubota & Done, 2019). These disks
are geometrically thick and can readily be described by our simulations. The
effect of changing the adiabatic index may influence the frequency and
amplitude of accretion rate oscillations, especially if the flow undergoes a
shock compression (Palit et al., 2019).

##  5 Conclusions

We performed a series of full 3D general relativistic numerical simulations of
magnetically arrested disk and jets launched from these central engines. We
found that the system evolution is governed by the physical parameters of the
engine, such as the black hole spin, and disk size, as well as disk
magnetization, and we applied our scenarios to typical types of sources for
both AGN and GRB classes. We found that the MAD scenario is applicable to AGN
engines and supports persistent jet emissions. The variability pattern is,
however, consistent with currently available observations only if we assume
high black hole spins. We also find that while the MAD scenario can be applied
to GRBs, our jet power estimates are quite large and cannot explain faint
sources. Still, these simulations give a variability pattern that is roughly
consistent with GRB power density spectra. Finally, we found that in some
cases, strong magnetic fields may lead to jet quenching and this effect is
found to be important mainly for gamma ray burst jets. Finally we speculate
that this may be related to the strength of magnetically driven winds from the
GRB engines.

###### Acknowledgements.

We thank Bozena Czerny, Marek Sikora and Krzysztof Nalewajko for helpful
discussions. We also thank Kostas Sapountzis, for help in running simulations,
and preparing some plotting scripts in Jupyter Notebook. This work was
supported in part by the grant 2019/35/B/ST9/04000 from the Polish National
Science Center. We also acknowledge support from the Interdisciplinary Center
for Mathematical Modeling of the Warsaw University, and the PL-Grid
infrastructure, through the computational grant plggrb5.

## References

  * Acciari et al. (2009) Acciari, V. A., Aliu, E., Arlen, T., et al. 2009, Science, 325, 444
  * Aharonian et al. (2003) Aharonian, F., Akhperjanian, A., Beilicke, M., et al. 2003, A&A, 403, L1
  * Aloy et al. (2005) Aloy, M. A., Janka, H.-T., & Müller, E. 2005, A&A, 436, 273
  * Avara et al. (2016) Avara, M. J., McKinney, J. C., & Reynolds, C. S. 2016, MNRAS, 462, 636
  * Baglio et al. (2018) Baglio, M. C., Russell, D. M., Casella, P., et al. 2018, ApJ, 867, 114
  * Balbus & Hawley (1991) Balbus, S. A. & Hawley, J. F. 1991, ApJ, 376, 214
  * Beckwith et al. (2008) Beckwith, K., Hawley, J. F., & Krolik, J. H. 2008, ApJ, 678, 1180
  * Beloborodov et al. (2000) Beloborodov, A. M., Stern, B. E., & Svensson, R. 2000, ApJ, 535, 158
  * Blandford & Znajek (1977) Blandford, R. D. & Znajek, R. L. 1977, MNRAS, 179, 433
  * Broderick et al. (2011) Broderick, A. E., Fish, V. L., Doeleman, S. S., & Loeb, A. 2011, ApJ, 735, 110
  * Chatterjee et al. (2021) Chatterjee, R., Das, S., Khasnovis, A., et al. 2021, Journal of Astrophysics and Astronomy, 42, 80
  * Christie et al. (2019) Christie, I. M., Lalakos, A., Tchekhovskoy, A., et al. 2019, MNRAS, 490, 4811
  * Corsi & Lazzati (2021) Corsi, A. & Lazzati, D. 2021, New A Rev., 92, 101614
  * Cuadra et al. (2005) Cuadra, J., Nayakshin, S., Springel, V., & Di Matteo, T. 2005, MNRAS, 360, L55
  * Czerny et al. (2016) Czerny, B., Du, P., Wang, J.-M., & Karas, V. 2016, ApJ, 832, 15
  * Davis & Tchekhovskoy (2020) Davis, S. W. & Tchekhovskoy, A. 2020, ARA&A, 58, 407
  * Dexter et al. (2014) Dexter, J., McKinney, J. C., Markoff, S., & Tchekhovskoy, A. 2014, MNRAS, 440, 2185
  * Dichiara et al. (2016) Dichiara, S., Guidorzi, C., Amati, L., Frontera, F., & Margutti, R. 2016, A&A, 589, A97
  * Dichiara et al. (2013) Dichiara, S., Guidorzi, C., Frontera, F., & Amati, L. 2013, ApJ, 777, 132
  * Fender et al. (2004) Fender, R. P., Belloni, T. M., & Gallo, E. 2004, MNRAS, 355, 1105
  * Fernández et al. (2019) Fernández, R., Tchekhovskoy, A., Quataert, E., Foucart, F., & Kasen, D. 2019, MNRAS, 482, 3373
  * Fishbone & Moncrief (1976) Fishbone, L. G. & Moncrief, V. 1976, ApJ, 207, 962
  * Gammie et al. (2003) Gammie, C. F., McKinney, J. C., & Tóth, G. 2003, ApJ, 589, 444
  * Giommi et al. (2021) Giommi, P., Perri, M., Capalbi, M., et al. 2021, MNRAS, 507, 5690
  * Globus & Levinson (2016) Globus, N. & Levinson, A. 2016, MNRAS, 461, 2605
  * Godfrey et al. (2012) Godfrey, L. E. H., Lovell, J. E. J., Burke-Spolaor, S., et al. 2012, ApJ, 758, L27
  * Guidorzi et al. (2016) Guidorzi, C., Dichiara, S., & Amati, L. 2016, A&A, 589, A98
  * Harris & Krawczynski (2006) Harris, D. E. & Krawczynski, H. 2006, ARA&A, 44, 463
  * Heger et al. (2005) Heger, A., Woosley, S. E., & Spruit, H. C. 2005, ApJ, 626, 350
  * Hossein Nouri et al. (2018) Hossein Nouri, F., Duez, M. D., Foucart, F., et al. 2018, Phys. Rev. D, 97, 083014
  * Hubrig et al. (2020) Hubrig, S., Järvinen, S. P., Schöller, M., & Hummel, C. A. 2020, MNRAS, 491, 281
  * Igumenshchev (2008) Igumenshchev, I. V. 2008, ApJ, 677, 317
  * Janiuk (2017) Janiuk, A. 2017, ApJ, 837, 39
  * Janiuk (2019) Janiuk, A. 2019, ApJ, 882, 163
  * Janiuk et al. (2021) Janiuk, A., James, B., & Palit, I. 2021, ApJ, 917, 102
  * Janiuk et al. (2013) Janiuk, A., Mioduszewski, P., & Moscibrodzka, M. 2013, ApJ, 776, 105
  * Janiuk et al. (2018) Janiuk, A., Sapountzis, K., Mortier, J., & Janiuk, I. 2018, Supercomputing Frontiers and Innovations, 5
  * Janiuk et al. (2018) Janiuk, A., Sukova, P., & Palit, I. 2018, ApJ, 868, 68
  * Janiuk & Yuan (2010) Janiuk, A. & Yuan, Y. F. 2010, A&A, 509, A55
  * Kathirgamaraju et al. (2019) Kathirgamaraju, A., Tchekhovskoy, A., Giannios, D., & Barniol Duran, R. 2019, MNRAS, 484, L98
  * Kobayashi et al. (1997) Kobayashi, S., Piran, T., & Sari, R. 1997, ApJ, 490, 92
  * Krolik (1999) Krolik, J. H. 1999, Active galactic nuclei : from the central black hole to the galactic environment
  * Kuan et al. (2021) Kuan, H.-J., Suvorov, A. G., & Kokkotas, K. D. 2021, MNRAS, 508, 1732
  * Kubota & Done (2019) Kubota, A. & Done, C. 2019, MNRAS, 489, 524
  * Lazzati et al. (2017) Lazzati, D., Deich, A., Morsony, B. J., & Workman, J. C. 2017, MNRAS, 471, 1652
  * Lazzati et al. (2021) Lazzati, D., Perna, R., Ciolfi, R., et al. 2021, ApJ, 918, L6
  * Lee & Ramirez-Ruiz (2007) Lee, W. H. & Ramirez-Ruiz, E. 2007, New Journal of Physics, 9, 17
  * Levinson & Globus (2017) Levinson, A. & Globus, N. 2017, MNRAS, 465, 1608
  * Liska et al. (2020) Liska, M., Tchekhovskoy, A., & Quataert, E. 2020, MNRAS, 494, 3656
  * Liu et al. (2015) Liu, T., Hou, S.-J., Xue, L., & Gu, W.-M. 2015, ApJS, 218, 12
  * Lloyd-Ronning et al. (2018) Lloyd-Ronning, N., Lei, W.-h., & Xie, W. 2018, MNRAS, 478, 3525
  * MacLachlan et al. (2013) MacLachlan, G. A., Shenoy, A., Sonbas, E., et al. 2013, MNRAS, 432, 857
  * McKinney & Gammie (2004) McKinney, J. C. & Gammie, C. F. 2004, ApJ, 611, 977
  * McKinney et al. (2012) McKinney, J. C., Tchekhovskoy, A., & Blandford, R. D. 2012, MNRAS, 423, 3083
  * McKinney et al. (2014) McKinney, J. C., Tchekhovskoy, A., Sadowski, A., & Narayan, R. 2014, MNRAS, 441, 3177
  * Melia et al. (2001) Melia, F., Bromley, B. C., Liu, S., & Walker, C. K. 2001, ApJ, 554, L37
  * Mizuta et al. (2018) Mizuta, A., Ebisuzaki, T., Tajima, T., & Nagataki, S. 2018, MNRAS, 479, 2534
  * Mochkovitch et al. (1993) Mochkovitch, R., Hernanz, M., Isern, J., & Martin, X. 1993, Nature, 361, 236
  * Moderski et al. (1998) Moderski, R., Sikora, M., & Lasota, J. P. 1998, MNRAS, 301, 142
  * Morales Teixeira et al. (2018) Morales Teixeira, D., Avara, M. J., & McKinney, J. C. 2018, MNRAS, 480, 3547
  * Morsony et al. (2010) Morsony, B. J., Lazzati, D., & Begelman, M. C. 2010, ApJ, 723, 267
  * Mościbrodzka et al. (2021) Mościbrodzka, M., Janiuk, A., & De Laurentis, M. 2021, MNRAS, 508, 4282
  * Murguia-Berthier et al. (2021) Murguia-Berthier, A., Ramirez-Ruiz, E., De Colle, F., et al. 2021, ApJ, 908, 152
  * Narayan et al. (2012) Narayan, R., SÄ dowski, A., Penna, R. F., & Kulkarni, A. K. 2012, MNRAS, 426, 3241
  * Neilsen et al. (2013) Neilsen, J., Nowak, M. A., Gammie, C., et al. 2013, ApJ, 774, 42
  * Noble et al. (2006) Noble, S. C., Gammie, C. F., McKinney, J. C., & Del Zanna, L. 2006, ApJ, 641, 626
  * Obergaulinger & Aloy (2022) Obergaulinger, M. & Aloy, M. Á. 2022, MNRAS, 512, 2489
  * Palit et al. (2019) Palit, I., Janiuk, A., & Sukova, P. 2019, MNRAS, 487, 755
  * Paschalidis et al. (2015) Paschalidis, V., Ruiz, M., & Shapiro, S. L. 2015, ApJ, 806, L14
  * Ponce et al. (2014) Ponce, M., Palenzuela, C., Lehner, L., & Liebling, S. L. 2014, Phys. Rev. D, 90, 044007
  * Ponti et al. (2015) Ponti, G., De Marco, B., Morris, M. R., et al. 2015, MNRAS, 454, 1525
  * Prieto et al. (2016) Prieto, M. A., Fernández-Ontiveros, J. A., Markoff, S., Espada, D., & González-Martín, O. 2016, MNRAS, 457, 3801
  * Punsly (2007) Punsly, B. 2007, MNRAS, 374, L10
  * Quataert et al. (1999) Quataert, E., Narayan, R., & Reid, M. J. 1999, ApJ, 517, L101
  * Reichart et al. (2001) Reichart, D. E., Lamb, D. Q., Fenimore, E. E., et al. 2001, ApJ, 552, 57
  * Rhoads (1999) Rhoads, J. E. 1999, ApJ, 525, 737
  * Ripperda et al. (2022) Ripperda, B., Liska, M., Chatterjee, K., et al. 2022, ApJ, 924, L32
  * Ruiz et al. (2020) Ruiz, M., Tsokaros, A., & Shapiro, S. L. 2020, Phys. Rev. D, 101, 064042
  * Rusinek et al. (2020) Rusinek, K., Sikora, M., Kozieł-Wierzbowska, D., & Gupta, M. 2020, ApJ, 900, 125
  * Sapountzis & Janiuk (2019) Sapountzis, K. & Janiuk, A. 2019, ApJ, 873, 12
  * Sari et al. (1999) Sari, R., Piran, T., & Halpern, J. P. 1999, ApJ, 519, L17
  * Scepi et al. (2021) Scepi, N., Begelman, M. C., & Dexter, J. 2021, MNRAS, 502, L50
  * Shibata et al. (2021) Shibata, M., Fujibayashi, S., & Sekiguchi, Y. 2021, Phys. Rev. D, 104, 063026
  * Siegel & Metzger (2018) Siegel, D. M. & Metzger, B. D. 2018, ApJ, 858, 52
  * Sikora & Begelman (2013) Sikora, M. & Begelman, M. C. 2013, ApJ, 764, L24
  * Sikora et al. (2013) Sikora, M., Stasińska, G., Kozieł-Wierzbowska, D., Madejski, G. M., & Asari, N. V. 2013, ApJ, 765, 62
  * Tamburini et al. (2020) Tamburini, F., Thidé, B., & Della Valle, M. 2020, MNRAS, 492, L22
  * Tchekhovskoy & Giannios (2015) Tchekhovskoy, A. & Giannios, D. 2015, MNRAS, 447, 327
  * Tchekhovskoy et al. (2011) Tchekhovskoy, A., Narayan, R., & McKinney, J. C. 2011, MNRAS, 418, L79
  * Ukwatta et al. (2011) Ukwatta, T. N., Dhuga, K. S., Morris, D. C., et al. 2011, MNRAS, 412, 875
  * van Velzen & Falcke (2013) van Velzen, S. & Falcke, H. 2013, A&A, 557, L7
  * Witzel et al. (2018) Witzel, G., Martinez, G., Hora, J., et al. 2018, ApJ, 863, 15
  * Witzel et al. (2021) Witzel, G., Martinez, G., Willner, S. P., et al. 2021, ApJ, 917, 73
  * Wu et al. (2016) Wu, Q., Zhang, B., Lei, W.-H., et al. 2016, MNRAS, 455, L1