© A. Janiuk et al. 2022

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 ∼ Ṁc² (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 between three and four times the accretion power, that is, P_MAD = 1.3 h_0.3a² Ṁc² (Davis & Tchekhovskoy 2020). In terms of the magnetic pressure driven onto the black hole, which is proportional to the square of the magnetic flux, ${\mathrm{\Phi }}_{\mathrm{BH}}^2$, 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, ${\varphi }_{\mathrm{BH}}={\mathrm{\Phi }}_{\mathrm{BH}}/\sqrt{\dot{M}{r}_{\mathrm{g}}{c}^2}$, defined at the black hole horizon, gives the characteristic threshold for magnetically arrested state. In this state, the dimensionless magnetic flux ϕ_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 ∼1 s 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. 2018a). 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. 2018b). Our code is publicly available under the name 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 Sect. 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 Sect. 3. Our discussion and conclusions make up Sect. 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. 2018a). 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^μν, is contributed by the gas and electromagnetic fields:

$$\begin{array}{cc}& {T_{\left(m\right)}}^{\mu \nu }=\rho \xi u^{\mu }{u}^{\nu }+p{g}^{\mu \nu },\hfill \\ \hfill & {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 },\hfill \\ \hfill & T^{\mu \nu }={T_{\left(m\right)}}^{\mu \nu }+{T_{\left(em\right)}}^{\mu \nu }.\hfill \end{array}$$

Here u^μ is the four-velocity of gas, u denotes internal energy density, b^μ is the magnetic four-vector and ξ is the fluid specific enthalpy. The continuity and momentum conservation equation is expressed as:

$$\begin{array}{c}\hfill {(\rho u^{\mu })}_{;\mu }=0,T_{\nu ;\mu }^{\mu }=0.\end{array}$$(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, (r,θ,ϕ), the black hole is located at 0 <  r ≤ r_h, where $r_{\mathrm{h}}=(1+\sqrt{1-a^2})r_{\mathrm{g}}$ is the horizon radius of a rotating black hole with mass, M, and angular momentum, J, in geometrized units, r_g = GM/c², and a is the dimensionless Kerr parameter, a = J/(Mc), 0 ≤ a ≤ 1. In our simulations, we adopt the rotating black hole, with a = 0.9.

The HARM code does not perform the integration in the Boyer-Lindquist coordinates, but instead in the so called Modified Kerr-Schild ones: (t,x⁽¹⁾,x⁽²⁾,ϕ), (Noble et al. 2006). The transformation between the coordinate systems is given by:

$$\begin{array}{cc}& r=R_0+exp\left[x^{(1)}\right]\hfill \\ \hfill & \theta =\frac{\pi }{2}(1+x^{(2)})+\frac{1-h}{2}sin\left[\pi (1+x^{(2)})\right],\hfill \end{array}$$

where R₀ is the innermost radial distance of the grid, 0 ≤ x⁽²⁾ ≤ 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⁽²⁾, we obtain an equally spaced grid on θ, 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_cyl = 3.5 and inside the angle θ_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⁵ gravitational radii, while the grid cells are scaled with distance from the center. The inner and outer radius are set to $R_{\mathrm{in}}=0.87\times (1.+\sqrt{(}1.-a^2))$, and R_out = 10⁵. The progressively sparser grid starts at R_br = 400, with $x_i^{(1)}=x_{i-1}^{(1)}+1/C\times \delta x$, where $\delta x=x_i^{(1)}-{(log(R_{\mathrm{br}}-R_0))}^4$. We use C = 1 to limit the hyper-exponentiation. Our grid resolution is N1 × N2 × N3 = 288 × 256 × 256 in the r, θ, and ϕ 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 θ-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 θ.

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 Fig. 1.

Fig. 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_in, and the distance where the maximum pressure occurs, r_max. In the current simulations, we adopt values 6 (12) and 13 (25), for r_in and r_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_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,

$$\begin{array}{c}\hfill A_{\phi }=max((\frac{\overline{\rho }}{{\rho }_{\max }}-{\rho }_0)\times r^5,0),\end{array}$$(2)

where $\overline{\rho }$ is the density in the torus (averaged over 2 neighboring cells) and ρ_max is the density maximum, and we use ρ₀ = 0.2, to restrict the initial field to the regions where density in the torus is larger than 0.2 ρ_max. The factor of r⁵ 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_θ = 0; hence, the magnetic field vector will have only B_r and B_θ components The vector potential is related to the Faraday tensor, F_μν = ∂_μA_μ − ∂_νA_μ, and we thus obtain the magnetic field four-vectors with b^μ = −^*F^μνu_ν. It is related to the normal observers magnetic field vector, as B^μ = α^*F^μt, where α = 1/g_tt is the lapse (noting that there is no time-component for B^μ).

The magnetized flows are, by definition, not in equilibrium and the angular momentum is transported. The plasma β-parameter is defined as the ratio of the fluid’s thermal to the magnetic pressure, β ≡ p_g/p_mag. We normalize the magnetic field in the FM torus to have $\beta =(\gamma -1.)\ast u_{\max }/(0.5\ast b_{\max }^2)$, where u_max is the internal energy at the torus pressure maximum radius. We examined several models with different minimal initial β 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 = (γ − 1) u, and we use the adiabatic index of γ = 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.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_BH; the torus mass; and black hole dimensionless spin, a. Typical short GRB parameters are M_BH = 3 M_⊙, and a = 0.6 − 0.9, while the torus mass is about 0.1 M_⊙. 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_{\mathrm{unit}}=M_{\mathrm{unit}}/L_{\mathrm{unit}}^3=3.4\times {10}^{11}$ g cm⁻³. The magnetic flux is then normalized with ${\mathrm{\Phi }}_{\mathrm{unit}}=\sqrt{4\pi D_{\mathrm{unit}}{c}^2}{L}_{\mathrm{unit}}^2$ (see Table 1).

In case of AGN, we can apply our models to low-luminosity AGN, such as Galaxy center Sgr A^⋆ or M 87. 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⁻⁵ M_⊙ yr⁻¹, for these AGN (Quataert et al. 1999; Cuadra et al. 2005; Prieto et al. 2016). We adopted the black hole mass of 4 × 10⁶ M_⊙ and 5 × 10⁹ M_⊙, for Sgr A^⋆ and M 87, respectively, and we derive our length units, according to their gravitational radius. We apply the scaling of Sgr A^⋆ 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 M 87 case (Tamburini et al. 2020). The magnetic flux on the black hole horizon is expected to be on the order of 10⁻¹⁰ − 0.01 G pc² (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, ρ_ISM ∼ 10⁻²⁷ g cm⁻³. 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⁵ 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 × 10⁻¹² g cm⁻³ and Φ_unit = 9.08 × 10²⁸ G cm² and then D_unit = 8.85 × 10⁻¹⁸ g cm⁻³ and Φ_unit = 1.15 × 10³² G cm², for Sgr A^⋆ and M 87, 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.

Fig. 2. 3D view of the inner parts of torus in the end of simulation with β = 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 ∼100 rg. 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 × 107 cm.

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 β-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 ∼50 r_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_f = 50 000 t_g (geometrical units, t_g = GM_BH/c³). This corresponds to 0.74 s in the GRB case and to about 9.4 months in our AGN case.

In Fig. 3, we show the profiles of density, and magnetic field streamlines, in the r − θ polar plane (XZ), and in the equatorial r − ϕ plane (XY), as obtained at the mid-time of the simulation. We show three different times during the simulation, t = 20 000 M t = 25 000 M, and t = 30 000 M, 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 β = 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).

Fig. 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 β = 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 000 M, t = 25 000 M, and t = 30 000 M. Streamlines show the magnetic three-vector configuration. We show the equatorial cuts of a small region of 10 rg from the black hole (bottom panels) and the polar cuts are shown for small and large regions, up to 100 rg (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 β, 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 β = 100 case are seen in the density range in the equatorial plane, which is (on average) larger than small β 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

In Fig. 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:

$$\begin{array}{cc}& {\dot{m}}_{\mathrm{BH}}=-{\displaystyle \int \rho u^{\mathrm{r}}\sqrt{-g}\mathrm{d}\theta \mathrm{d}\varphi },\hfill \end{array}$$(3)

$$\begin{array}{cc}& {\mathrm{\Phi }}_{\mathrm{BH}}=0.5{\displaystyle \int |B^{\mathrm{r}}|\sqrt{-g}}\mathrm{d}\theta \mathrm{d}\varphi ,\hfill \end{array}$$(4)

and

$$\begin{array}{c}\hfill {\varphi }_{\mathrm{BH}}=\frac{{\mathrm{\Phi }}_{\mathrm{BH}}}{\sqrt{{\dot{m}}_{\mathrm{BH}}}}·\end{array}$$(5)

Fig. 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, ϕBH, given by Eq. (5), (bottom panel). Blue lines denote the models with initial gas-to-magnetic pressure ratio, β = 100, cyan lines present the results for β = 50, and green lines denote β = 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 rg (GRB-type models).

As shown in Fig. 4, accretion rate is large at the beginning of the simulation and it is higher then for smaller β 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 ϕ_BH ∼ 10 for most of the simulation after the initial condition is relaxed. The quantitative difference in ϕ_BH(time) profile between the models with β = 100 and β = 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).

In Fig. 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 β-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, Φ_B, is almost the same for all three values of β. 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 ϕ_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).

Fig. 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, ϕBH, given by Eq. (5), (bottom panel). Red lines denote the models with initial gas-to-magnetic pressure ratio β = 100, orchid lines present results for β = 50, and magenta lines denote β = 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 rg (GRB-type models).

Finally, we present in Fig. 6 a similar set of models, assuming, for the black hole spin, a = 0.9, and three-disk magnetization, normalized to β = 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 = 25 r_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, ϕ_BH, rises smoothly after an initial relaxation period of about a 1000 t_g. The differences between the results obtained for three values of β-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.

Fig. 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, ϕBH, given by Eq. (5), (bottom panel). Blue lines denote the models with initial gas-to-magnetic pressure ratio β = 100, cyan lines present results for β = 50, and green lines denote β = 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 rg (AGN-type models).

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 Figs. 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 δt = 10 000 t_g 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.

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

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 t_g). 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.

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 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):

$$\begin{array}{c}\hfill P_{\mathrm{BZ}}=-{\displaystyle {\int }_{r=r_{\mathrm{h}}}}(b^2{u}^{\mathrm{r}}{u}_{\mathrm{t}}-b^{\mathrm{r}}{b}_{\mathrm{t}})\mathrm{d}A,\end{array}$$(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_{\mathrm{t}}^{\mathrm{r}}$, 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 ϕ_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 ∼10⁴² to ∼10⁵² erg s⁻¹, 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 Fig. 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 β, 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 β-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 β = 50 case, albeit at the end of the run. This behavior is correlated with the decrease of the dimensionless magnetic flux, ϕ_BH (cf. Fig. 4).

Fig. 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 β-parameters, shown with blue, red, and green points for β = 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.

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_ϕ/B_z ratio in this model is high at the disk equator up to the distance of ∼10 r_g, because of the energy transfer from the rotating black hole to the torus via the magnetic coupling (Janiuk & Yuan 2010).

Fig. 9. Ratio of the toroidal to poloidal magnetic fields, Bϕ/Bz, for models AGN-NS50 and AGN-HS50, shown in the top and bottom panels, respectively. Color maps are in logarithmic scale.

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

$$\begin{array}{c}\hfill \mu =-\frac{T_{\mathrm{t}}^{\mathrm{r}}}{\rho u^{\mathrm{r}}},\end{array}$$(7)

Fig. 10. Profiles of jet energetics, defined by μ-parameter (see Eq. (7)) for two GRB models with high spin, a = 0.9, and β = 100 or β = 50, on the two top, and the two bottom panels, respectively. The maps show profiles taken at two opposite azimuthal angles, Φ = 0 and Φ = π. All snapshots are taken at a time of t = 2.5 × 104 tg.

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 β = 100 and β = 50. We checked the trends also for β = 30 and models with smaller black hole spins. The main relation between β and the jet opening angle is that for larger gas-to-magnetic pressure ratio, the jets seem to be more symmetric, while for smaller β, 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).

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 20 000 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⁻³ − 10⁻⁷ (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 Fig. 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):

$$\begin{array}{c}\hfill {\lambda }_{\mathrm{MRI}}=\frac{2\pi }{\mathrm{\Omega }}\frac{b}{\sqrt{4\pi \rho h+b^2}}·\end{array}$$(8)

Fig. 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 ∼ 25 000 M (right).

The MRI is well resolved, when $Q_{\mathrm{MRI}}=\frac{{\lambda }_{\mathrm{MRI}}}{\mathrm{\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.

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⁻⁸ and 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 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_BHc⁻³. 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^⋆ 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 α = 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 α = 1.4 ± 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, α = 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 α = 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^⋆ are best fitted by PDS with a slope of α = 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^⋆ 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 M 87, the timescale of 30 t_g is equal to about 8.5 days (cf. Table 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 Fig. 8 confirm this effect. In the models with β = 30, the wind density is high, while the jet power is smaller than for the higher β 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) × 10²⁸ up to 1.6 × 10²⁹ G cm² 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 × 10⁻⁷ G pc², 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⁻⁹ G pc², 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 γ = 4/3 is used for both AGN and GRB scenarios. In fact, the ideal gas law with γ = 5/3 could be used to describe low-luminosity AGN, while the relativistically hot gas with γ = 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.

Acknowledgments

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

All Tables

All Figures

	Fig. 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.
In the text

Fig. 2. 3D view of the inner parts of torus in the end of simulation with β = 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 ∼100 rg. 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 × 107 cm.

In the text

Fig. 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 β = 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 000 M, t = 25 000 M, and t = 30 000 M. Streamlines show the magnetic three-vector configuration. We show the equatorial cuts of a small region of 10 rg from the black hole (bottom panels) and the polar cuts are shown for small and large regions, up to 100 rg (upper panels).

In the text

Fig. 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, ϕBH, given by Eq. (5), (bottom panel). Blue lines denote the models with initial gas-to-magnetic pressure ratio, β = 100, cyan lines present the results for β = 50, and green lines denote β = 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 rg (GRB-type models).

In the text

Fig. 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, ϕBH, given by Eq. (5), (bottom panel). Red lines denote the models with initial gas-to-magnetic pressure ratio β = 100, orchid lines present results for β = 50, and magenta lines denote β = 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 rg (GRB-type models).

In the text

Fig. 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, ϕBH, given by Eq. (5), (bottom panel). Blue lines denote the models with initial gas-to-magnetic pressure ratio β = 100, cyan lines present results for β = 50, and green lines denote β = 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 rg (AGN-type models).

In the text

	Fig. 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.
In the text

	Fig. 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 β-parameters, shown with blue, red, and green points for β = 100, 50, and 30, respectively.
In the text

	Fig. 9. Ratio of the toroidal to poloidal magnetic fields, Bϕ/Bz, for models AGN-NS50 and AGN-HS50, shown in the top and bottom panels, respectively. Color maps are in logarithmic scale.
In the text

	Fig. 10. Profiles of jet energetics, defined by μ-parameter (see Eq. (7)) for two GRB models with high spin, a = 0.9, and β = 100 or β = 50, on the two top, and the two bottom panels, respectively. The maps show profiles taken at two opposite azimuthal angles, Φ = 0 and Φ = π. All snapshots are taken at a time of t = 2.5 × 104 tg.
In the text

	Fig. 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 ∼ 25 000 M (right).
In the text