# General Relativistic Magnetohydrodynamic Simulations of Accreting Tori:
Resolution Study

Angelos Karakonstantakis  1,a Debora Lančová  []2,b and Miljenko Čemeljić
[]2,1,3
1Nicolaus Copernicus Astronomical Centre  Polish Academy of Sciences  [1] ul.
Bartycka 18  00-716 Warsaw  Poland
2Research Centre for Computational Physics and Data Processing  [3] Institute
of Physics  Silesian University in Opava  [2] Bezručovo nám. 13  CZ-746 01
Opava  Czech Republic
3Academia Sinica  Institute of Astronomy and Astrophysics  [3] P.O. Box 23-141
Taipei 106  Taiwan
<EMAIL_ADDRESS>
<EMAIL_ADDRESS>
###### Abstract

We present two-dimensional general relativistic radiative
magnetohydrodynamical simulations of accretion disks around non-rotating
stellar-mass black hole. We study the evolution of an equilibrium accreting
torus in different grid resolutions to determine an adequate resolution to
produce a stable turbulent disk driven by magneto-rotational instability. We
evaluate the quality parameter, $Q_{\theta}$, from the ratio of MRI wavelength
to the grid zone size and examine the effect of resolution in various
quantitative values such as the accretion rate, magnetisation, fluxes of
physical quantities and disk scale-height. We also analyse how the resolution
affects the formation of plasmoids produced in the magnetic reconnection
events.

###### keywords:

accretion, accretion disks – black hole physics – radiation GRMHD

A. Karakonstantakis, D. Lančová, M. ČemeljićGeneral Relativistic
Magnetohydrodynamic Simulations of Accreting Tori: Resolution Study

##  1 Introduction

Global General Relativistic (GR) Radiative (R) Magnetohydrodynamical (MHD)
simulations are an excellent tool to model the violent and turbulent
environment in the vicinity of an accreting compact object. Even though this
method can capture the interplay of many different physical processes
governing the behaviour of magnetised plasma, it also has many limitations.

One significant limitation is that the GRRMHD solves the flux of quantities on
a fixed discrete grid, neglecting the structure inside the grid cells and
approximating them as homogeneous. The solution to the related Riemann problem
is always only approximate; thus, some information on the low-scale
fluctuations is lost. At the same time, accretion in magnetised disks is known
to be driven by low-scale turbulences caused by magneto-rotational instability
(MRI, Balbus and Hawley, 1991, 1998).

The smaller cell size leads to a smaller time step. The issue becomes even
more complex in the case of the GR due to singularities in the commonly used
coordinate systems and limitation by the speed of light in the case of
radiative MHD. A suitable choice of computational grid parameters is
fundamental in GRMHD simulations, as it must have a sufficient resolution to
resolve the low-scale turbulences while considering computational demands.

The challenge of the GRMHD method is the scalability of turbulences present in
the accretion flow. This makes it difficult to achieve a fully converged
solution. Convergence in an accretion disk simulation is often determined by
the stability of the flow parameters (e.g., mass accretion rate, density scale
height, radiative luminosity). Recently, White et al. (2019) studied numerical
convergence in simulations of magnetically arrested disks (MAD) and its
effects in the lunch Jets. They studied the effect of 4 different resolutions.
They found that turbulent structure does not appear in their lowest
resolution, and MRI is suppressed. Also, Ripperda et al. (2022) conducted
extreme-resolution $10^{4}$ grid shell resolution 3D simulations of MAD flows
and analyzed the appearance of blob structures. They showed that magnetic flux
bundles can escape from the event horizon through reconnection and complete a
full orbit as low-density hot spots consistent with interferometric
observations.

Quality parameters were derived to assess simulation convergence in studies by
Hawley et al. (2011, 2013). They are based on the assumption that the
characteristic wavelength of the MRI must be well resolved in order to
maintain the MRI-induced turbulences. Although satisfying these conditions is
important for obtaining reliable and physical results, other limitations and
approximations of the GRMHD must still be considered.

Here, we investigate the influence of increasing resolution on simulation
results. The structure of this Proceeding is as follows: in Section 2, we
describe the employed code KORAL and the numerical setup of our simulation. In
Section 3, we describe the diagnostic tools we considered; in Section 4, we
present the analysis of the simulation results, and in Section 5, we comment
on our finding and draw conclusions.

##  2 Numerical setup and initial condition

###  2.1 KORAL implementation

We use the GRRMHD code KORAL to perform global two-dimensional (2D)
axisymmetric simulations of accretion onto a $M=10\,{\mathrm{M}}_{\odot}$ non-
rotating black hole (BH)111${\mathrm{M}}_{\odot}$ denoting the solar mass. The
grid resolution is logarithmic in radius and stretches up to
$r_{\mathrm{out}}=1000\,r_{\mathrm{g}}$, where $r_{\mathrm{g}}\equiv GM/c^{2}$
is the gravitational radius. We considered four increasing-resolution
simulation grids described in Table 1.

KORAL solves the conservation equations for a fluid with rest-mass density
$\rho$, 4-velocity $u^{\mu}$, and stress-energy tensor $T^{\mu}{}_{\nu}$,
coupled with the radiation tensor ($R^{\mu}{}_{\nu}$) via the radiation
4-force density $G_{\nu}$. The equations are

$\nabla_{\mu}\left(\rho
u^{\mu}\right)=0,\quad\nabla_{\mu}T^{\mu}{}_{\nu}=G_{\nu},\quad\nabla_{\mu}R^{\mu}{}_{\nu}=-G_{\nu}.$
(1)

The MHD stress-energy tensor $T_{\mu}{}_{\nu}$ is given by

$T^{\mu}{}_{\nu}=\left(\rho+{u_{\mathrm{int}}}+{p_{\mathrm{gas}}}+b^{2}\right)u^{\mu}u_{\nu}+\left({p_{\mathrm{gas}}}+\frac{1}{2}b^{2}\right)\delta^{\mu}_{\nu}-b^{\mu}b_{\nu},$
(2)

where $b^{\mu}$ is the magnetic field four-vector, ${u_{\mathrm{int}}}$ and
${p_{\mathrm{gas}}}=(\gamma-1){u_{\mathrm{int}}}$ are the internal energy and
pressure of the gas in the comoving frame, with $\gamma$ being the adiabatic
index.

The implementation of the various processes in the KORAL code was described in
an extensive collection of published papers, including Sadowski et al. (2013,
2014); Sadowski and Narayan (2015); Sadowski et al. (2015, 2017); Chael et al.
(2017), and recently also summarised in Lančová (2023).

####  2.1.1 Mean-field dynamo

MRI-driven turbulences are responsible for angular momentum transport in BH
accretion disks. A physical mechanism is required to sustain the magnetic
fields in time against dissipation. Local and global simulations of accretion
disks have shown that shearing due to differential rotation induces a
turbulent dynamo capable of amplifying and sustaining magnetic fields against
dissipation (see Del Zanna et al., 2022, and references therein).

MHD dynamo cannot operate in axisymmetric plasma configurations, as follows
from Cowling’s anti-dynamo theorem (Brandenburg et al., 1995). Therefore,
capturing the growth of MRI self-consistently requires a three-dimensional
(3D) simulation, which is computationally expensive given the long evolution
times and requirements for high resolution. However, the KORAL code
implementation allows for replenishment of the magnetic field with an
artificial mean-field dynamo term in the case of a 2D axisymmetric
configuration, as introduced in Sadowski et al. (2015). The dynamo efficiently
counterbalances the magnetic field dissipation.

###  2.2 Initial state

All simulations presented in this work are initialised by a torus in
hydrodynamical equilibrium following Penna et al. (2013). Fig. 1 shows the
initial distribution of matter and magnetic field topology in the meridional
plane. The inner edge of the initial torus is located at $R_{\rm
in}=22\,r_{\mathrm{g}}$. It is initially threaded with a poloidal magnetic
field, such that the magnetization in the torus centre is
$\beta={p_{\mathrm{mag}}}/{p_{\mathrm{gas}}}=20$, where $p_{\mathrm{mag}}$ is
the magnetic pressure.

The initial angular momentum on the equatorial plane is set up as a fraction
of the Keplerian $\Omega$ multiplied by a factor $\xi=0.975$ for
$R>R_{1}=30\,r_{\mathrm{g}}$. Inside $R_{1}$, the angular momentum is
constant. The initial state corresponds to simulation r300a0 (but with single
loop magnetic field loop as shown in Fig. 1) in Sadowski et al. (2015), where
the results were compared to a 3D simulation with the same parameters but
without the dynamo term, leading to comparable time-averaged results. The
effects of different magnetic field configurations, including a similar looped
configuration, have been studied by KoloĹA and Janiuk (2020).

Figure 1: Snapshot of the initial state in our simulations, showing logarithm
of gas density. Contours show the initial poloidal magnetic field lines.

The simulations were run for $20\,000$ $t_{\mathrm{g}}$ in
total222$t_{\mathrm{g}}=r_{\mathrm{g}}/c$. After $t\sim
10\,000\,t_{\mathrm{g}}$ from the beginning, the simulations become more
stable as the accretion rate stabilises to a value close
to$\dot{M}=\int_{0}^{\pi}\rho u^{r}\sqrt{g}\mathrm{d}\theta\simeq
10\,\dot{M}_{\mathrm{Edd}}$333$\dot{M}_{\mathrm{Edd}}=L_{\mathrm{Edd}}/\eta
c^{2}$, where, $\eta=-L_{\rm tot}/\int_{0}^{\pi}\rho u^{r}\mathrm{d}\theta
c^{2}$ is the radiative efficiency with $L_{\rm tot}$ defined in eq. (6)..
From this point to the end of the simulation, we averaged the data to study
the mean properties of the accretion disk.

##  3 Diagnostic tools for testing resolution

Inadequate resolution can lead to a number of numerical artefacts, crashing of
the simulations, or even wrong results. For example, it may suppress the
growth of MRI when the wavelength (of at least) the fastest-growing
instability mode is longer than the grid shell size. The MRI produces
turbulences, but when the resolution is low, only a small range of wave
numbers can be captured, distorting the properties of those turbulences. If
the simulation fails to capture the turbulent structure initiated and
sustained by MRI fluctuations, angular momentum is not transported outwards,
which prevents mass accretion. However, high-resolution simulations are
computationally very expensive. Thus, it is useful to determine whether a
given simulation has adequate resolution to produce accurate quantitative
results.

One would presume that numerical results converge to an exact solution as the
grid size decreases. An approach to determine the accuracy of a given
simulation is to run again for a higher (e.g. doubled) resolution and check
for convergence in the results. As part of this investigation, we ran a series
of simulations with increasing resolution to test the convergence. However,
different quantities have different convergence rates and are subject to the
limitations of the numerical algorithm used. Still, additional effects beyond
resolution, such as resistivity or viscosity terms, significantly impact
magnetic stress for the same resolution. Moreover, the chaotic nature of the
turbulent processes introduces an additional layer of complexity in
interpreting the convergence. These intricacies make it challenging to reach a
definitive conclusion regarding the effect of resolution alone.

###  3.1 MRI quality parameter

The first diagnostic tool considered for resolution is derived from the MRI
characteristic wavelength, $\lambda^{\rm MRI}\propto|u_{\mathrm{A}}|/\Omega$,
where $u_{\mathrm{A}}$ is the AlfvĂŠn velocity and $\Omega$ is the angular
velocity. Sano et al. (2004) studied the linear evolution of MRI in local
shearing box simulation and found that at least 6 resolution grid zones should
resolve MRI wavelength ($\lambda^{\rm MRI}_{i}/\Delta x_{i}\geq 6$).

Common indicators for the proper resolution to resolve the MRI turbulences are
the quality parameters $Q_{i}\equiv\lambda^{\rm MRI}_{i}/\Delta x_{i}$ (first
introduced in Noble et al., 2010; Hawley et al., 2011), where $\Delta x_{i}$
denotes the grid size in the $i$-th direction. Noble et al. (2010) studied the
properties of the turbulence through global GRMHD simulations and suggested
that $Q_{\theta}\gtrsim 20$ is adequate for resolving MRI.

Furthermore, Hawley et al. (2011), who performed global MHD simulations using
a pseudo-Newtonian potential and examined previous results of shearing box
simulations, emphasise the importance of considering a resolution parameter in
the toroidal direction (i.e. $Q_{\phi}\gtrsim 20$ and $Q_{\theta}\gtrsim 10$).
The authors noted that the two parameters ($Q_{\phi},\ Q_{\theta}$) are not
independent of each other: higher values in one direction can compensate for
lower values in the other. For example $Q_{\theta}\gtrsim 10-15$ can
compensate for a lower value of $Q_{\phi}\approx 20$ (Sorathia et al., 2012).
Hence, the product of $Q_{\theta}Q_{\phi}\gtrsim 200$ has also been suggested
as a criterion for convergence (Narayan et al., 2012).

In our case of 2D axisymmetric simulations, only $Q_{\theta}$ can be defined.
Therefore, we consider $Q_{\theta}>20$ as adequate for resolving MRI. We
define $Q_{\theta}$ as:

$Q_{\theta}=\frac{2\pi}{\Omega\mathrm{d}\hat{x}^{\theta}}\frac{\left|\hat{b}^{\theta}\right|}{\sqrt{\rho}},$
(3)

where the magnetic field component $\hat{b}^{\theta}$ and the grid cell size
$\mathrm{d}\hat{x}^{\theta}$ are evaluated in the fluid frame (denoted with
hats). In addition, we consider the theta average of the quality parameter as
a function of radius. For this reason, we compute the density-weighted quality
parameter as a function of the cylindrical radius $R$:

$\left<Q_{\theta}\right>\left(R\right)=\frac{\int
Q_{\theta}\rho\mathrm{d}\theta}{\int\rho\mathrm{d}\theta}.$ (4)

##  4 Results

In this section, we describe selected properties of the accretion flow in the
simulations and their dependency on resolution. The quantities are shown
either as a function of time or their spatial or radial profile from time-
averaged data.

Fig. 2 shows the color-maps of the gas density and the radiation energy
density in the fluid frame ($\hat{E}$) from averaged data. The white solid
line shows the surface of the last scattering (photosphere) and the dashed
lines show the surface corresponds to density scale-height ($h_{\tau}$).
$h_{\tau}$, calculated in the $z$ direction. It is apparent that for a low
resolution, the funnel region is not formed properly, and the disk photosphere
is much higher than in the case of the high resolution. In the $(N_{r}\times
N_{\theta})=64^{2}$ grid shells simulation, the photosphere position cannot be
properly established at all, leading to the jump in the bottom half of the
domain. We also examined the scale-height $h_{\rho}$ of the simulated
accretion disk, calculated as

$h_{\rho}=\sqrt{\frac{\int\rho z^{2}\mathrm{d}z}{\int\rho\mathrm{d}z}},$ (5)

as shown with the dashed lines in Fig. 2 and as a function of radius in the
right panel of Fig. 5 with the various solid lines corresponding to different
resolutions. The time average scale-height is roughly stable for
$R>20\,r_{g}$, and a typical value is between 0.25–0.30.

Figure 2: The logarithm of radiation energy density in the fluid frame (left
panels) and time average density (right panels). The solid line surface shows
the position of the photosphere and the dashed line correspond to the density
scale-height.

Before discussing other results, we examine if the resolution is adequate for
capturing MRI following the discussion in Section 3. Fig. 3 shows the values
of $Q_{\theta}$ for each resolution. For higher resolutions, i.e.
$(N_{r}\times N_{\theta})=256^{2},\ 512^{2}$, the values are higher than 20
for the vast majority of the cells, which indicate adequate resolution (see
Section 3.1).

For the highest resolution, $Q_{\theta}>1000$, while for the lower resolution,
the values are $<10$. We can conclude that the resolution of $64^{2}$ is
inadequate to capture MRI. In the case of $(N_{r}\times N_{\theta})=128^{2}$,
$Q_{\theta}\sim 10$ for most cells. However, in the equatorial plane, the
values reach $Q_{\theta}\geq 20$. MRI may be properly resolved in that case,
especially if the resolution in the toroidal direction is high enough (i.e.
$Q_{\phi}\gtrsim 20$), which can compensate for lower values in $Q_{\theta}$.
We conclude that the resolution requirement is $(N_{r}\times
N_{\theta})=256^{2}$ or greater.

Figure 3: The values of resolution quality parameter $Q_{\theta}$ measured at
$t=15\times 10^{3}\,t_{\mathrm{g}}$.

In Fig. 4 (left panel), we show the time evolution of the $Q$ parameter
measured at$R=20\,r_{g}$ and theta as defined from eq. 4. We observe that the
values are initially very low and fluctuate until the disk is formed and
reaches inflow equilibrium for the measured radius, coinciding with the time
when the accretion rate becomes stable.

Figure 4: The time evolution of the quality parameter $\langle Q_{\theta}\rangle(R=20\,r_{\mathrm{g}})$ (left panel) and of the accretion rate measured close to the black hole event horizon (right panel). The vertical line indicates$t\geq 10\,000\,t_{\mathrm{g}}$ for which we average the results. Dashed lines were chosen for better clarity between overlapping lines. Table 1: The $10-20\times 10^{3}\,t_{\mathrm{g}}$ average and standard variance of the accretion rate at event horizon and resolution parameter at $R=20\,r_{\mathrm{g}}$. ($N_{r}\times N_{\theta}$) | $\langle\dot{M}\rangle$ | $\langle Q_{\theta}\rangle$
---|---|---
($\dot{M}_{\mathrm{Edd}}$) | $(R=20\,r_{\mathrm{g}})$
$64^{2}$ | $11.3\pm 7.2$ | $1.5\pm 1.1$
$128^{2}$ | $12.6\pm 5.2$ | $06\pm 38$
$256^{2}$ | $11.2\pm 4.4$ | $25\pm 32$
$512^{2}$ | $09.4\pm 3.8$ | $100\pm 100$

###  4.1 Accretion Rate

The time evolution of the accretion rate measured close to the BH event
horizon is shown in Fig. 4 (right panel). The mean value and variance of the
accretion rate after $10\,000\,t_{\mathrm{g}}$ (indicated by the vertical
line) until the end of each simulation is given in the second column of Table
1 while the third column shows the time average value of the resolution
parameter. The variance of the accretion rate fluctuations is decreasing with
increasing resolution. The very high variance of the resolution parameter in
the highest resolutions is because a vary high value at $t\sim
18\,000\,t_{\mathrm{g}}$ . However, the variability of 2D simulations using
the artificial dynamo term is shown to be higher than in the 3D simulations
(Sadowski et al., 2015), so these results should be interpreted with caution.

###  4.2 Magnetization

We calculated the magnetization, $\sigma=b^{2}/\rho$, across the resolutions.
The results are plotted in Fig. 5 (left panel). The magnetization as a
function of radial coordinate decreases with increasing radius for all
resolutions and the values converge as the resolution is increased with
$256^{2}$ and $512^{2}$ grid shells resolution being almost identical. Thus,
we can conclude that magnetization values are affected by resolution and tend
to converge at the same value as the resolution increases.

Figure 5: The magnetization (left panel) and the disk scale-height (right
panel) as a function of radius from the time-averaged data.

###  4.3 Luminosities

An interesting application of numerical simulations is to compare the
luminosity or spectra with observations. Therefore, in this section we examine
if these results converge. The total luminosity in all forms of energy may be
defined by integrating the flux carried by gas and radiation:

$L_{\rm tot}\,\left(R\right)=-\int_{0}^{\pi}\left(T^{r}{}_{t}+R^{r}{}_{t}+\rho
u^{r}\right)\sqrt{-g}\mathrm{d}\theta.$ (6)

We define the radiative luminosity:

$L_{\rm
rad}\,\left(R\right)=-\int_{0}^{\pi}R^{r}{}_{t}\sqrt{-g}\mathrm{d}\theta.$ (7)

Fig. 6 shows the values of $L_{\rm tot}$ (left panel) and $L_{\rm rad}$ (right
panel) for each resolution. While these values are affected by resolution,
they remain within the same order of magnitude. In the left panel of Fig. 7,
the radiation energy density measured in the fluid frame, ($\hat{E}$), is
identical for all simulations. Interestingly, the resolution does not affect
the radiative quantities as much as those related to MHD fluid.

Figure 6: Left panel: The value of total (solid lines) and radiative
luminosity (dashed lines) from eqs. 6, 7 respectively. Right panel: The
radiation energy density in the fluid frame.

###  4.4 Formation of Plasmoids

Figure 7: Snapshot at $t\sim 10\,000\,t_{\mathrm{g}}$ showing the logarithm of
density and the formation of plasmoids (inserted panel) in the simulation with
resolution $N_{r}\times N_{\theta}=512^{2}$ grid cells. Contours show the
magnetic field lines.

Plasmoids, which in 2D show as magnetic islands, are bubbles of magnetised gas
formed in elongated magnetic sheets, and their formation is connected with
magnetic reconnection. Current sheets subject to turbulence induced in
accretion disks may become thin enough to break apart into plasmoids. Because
the thickness of the sheet is limited by numerical resolution, in low
resolutions, plasmoids are not formed. In our results, only the highest
resolution of $512^{2}$ grid cells allows the current sheet to become thin
enough for plasmoids to form, as shown in the right panel of Fig. 7. This
observation underscores the critical role of resolution in capturing intricate
magnetic reconnection processes.

##  5 Discussion and conclusions

In this study, we conducted a series of radiative General Relativistic
Magnetohydrodynamics (GRMHD) simulations, systematically varying the
resolution to assess convergence. The selected resolutions were tested using
quality parameters as resolution diagnostics, analyzed to examine the key
properties of the accretion flow and determine the optimal resolution for
magnetohydrodynamic simulations of the accretion disk and assessing result
convergence. We summarize our main results below:

  1. 1.

We observed that the accretion rate stabilizes when the Magnetorotational
Instability (MRI) is properly resolved at $t>3\times 10^{3}\,t_{\mathrm{g}}$
as indicated by the value of quality parameters ($Q_{\theta}>20$). The
variance of the average value of accretion rate is decreasing with increasing
resolution.

  2. 2.

Resolutions lower than $(N_{R}\times N_{\theta})<128^{2}$ are inadequate to
capture MRI suggesting this resolution being a potential threshold. We
observed higher resolutions (i.e. $256^{2},\,512^{2}$) not to have a
significant quantitative difference.

  3. 3.

The calculation of magnetization showed convergence as resolution increased
while radiation energy density remain unaffected by resolution changes,
luminosity and scale-height remained relatively stable across resolutions.

  4. 4.

Low resolution cannot resolve the fine structure of blobs, plasmoids are only
appear in the highest resolution.

This comprehensive analysis demonstrates the critical importance of resolution
in accurately simulating and interpreting the dynamics of accretion flows.
Adequate resolution is essential for capturing phenomena such as MRI or fine
structural details such as the formation of plasmoids, ensuring the
reliability of simulation results in understanding the intricate processes
governing accretion disks.

This work was supported in part by the Polish NCN grants 2019/33/B/ST9/01564
and 2019/35/O/ST9/03965. MČ acknowledges the Czech Science Foundation (GAČR)
grant No. 21-06825X. DL acknowledge the internal grant of Silesian University,
$\mathrm{SGS/31/2023}$. A.K. thanks Włodek Kluźniak and JiĹĂ HorĂAk for the
suggestions, and useful conversations. We gratefully acknowledge Polandâs
high-performance computing infrastructure PLGrid (HPC Centers: ACK Cyfronet
AGH), for providing computer facilities and support within computational grant
no. $\mathrm{PLG}/2023/016168$. Part of the computations for this article have
been performed using computer cluster at CAMK PAN.

## References

  * Balbus and Hawley (1991) Balbus, S. A. and Hawley, J. F. (1991), A Powerful Local Shear Instability in Weakly Magnetized Disks. I. Linear Analysis, _The Astrophysical Journal_ , 376, p. 214.
  * Balbus and Hawley (1998) Balbus, S. A. and Hawley, J. F. (1998), Instability, turbulence, and enhanced transport in accretion disks, _Reviews of Modern Physics_ , 70(1), pp. 1–53.
  * Brandenburg et al. (1995) Brandenburg, A., Nordlund, A., Stein, R. F. and Torkelsson, U. (1995), Dynamo-generated Turbulence and Large-Scale Magnetic Fields in a Keplerian Shear Flow, _The Astrophysical Journal_ , 446, p. 741.
  * Chael et al. (2017) Chael, A. A., Narayan, R. and Sadowski, A. (2017), Evolving non-thermal electrons in simulations of black hole accretion, _Monthly Notices of the Royal Astronomical Society_ , 470(2), pp. 2367–2386, 1704.05092.
  * Del Zanna et al. (2022) Del Zanna, L., Tomei, N., Franceschetti, K., Bugli, M. and Bucciantini, N. (2022), General relativistic magnetohydrodynamics mean-field dynamos, _Fluids_ , 7(2), p. 87, ISSN 2311-5521, URL http://dx.doi.org/10.3390/fluids7020087.
  * Hawley et al. (2011) Hawley, J. F., Guan, X. and Krolik, J. H. (2011), Assessing Quantitative Results in Accretion Simulations: From Local to Global, _The Astrophysical Journal_ , 738(1), 84, 1103.5987.
  * Hawley et al. (2013) Hawley, J. F., Richers, S. A., Guan, X. and Krolik, J. H. (2013), Testing Convergence for Global Accretion Disks, _The Astrophysical Journal_ , 772(2), 102, 1306.0243.
  * KoloĹA and Janiuk (2020) KoloĹA, M. and Janiuk, A. (2020), Simulations of black hole accretion torus in various magnetic field configurations, 2004.07535.
  * Lančová (2023) Lančová, D. (2023), Computer modelling of accretion processes in binary systems with black holes and neutron stars, _arXiv e-prints_ , 2310.13152.
  * Narayan et al. (2012) Narayan, R., Sadowski, A., Penna, R. F. and Kulkarni, A. K. (2012), GRMHD simulations of magnetized advection-dominated accretion on a non-spinning black hole: role of outflows, _Monthly Notices of the Royal Astronomical Society_ , 426(4), pp. 3241–3259, 1206.1213.
  * Noble et al. (2010) Noble, S. C., Krolik, J. H. and Hawley, J. F. (2010), Dependence of Inner Accretion Disk Stress on Parameters: The Schwarzschild Case, _The Astrophysical Journal_ , 711(2), pp. 959–973, 1001.4809.
  * Penna et al. (2013) Penna, R. F., Kulkarni, A. and Narayan, R. (2013), A new equilibrium torus solution and GRMHD initial conditions, _A &A_, 559, A116, 1309.3680.
  * Ripperda et al. (2022) Ripperda, B., Liska, M., Chatterjee, K., Musoke, G., Philippov, A. A., Markoff, S. B., Tchekhovskoy, A. and Younsi, Z. (2022), Black hole flares: Ejection of accreted magnetic flux through 3d plasmoid-mediated reconnection, _The Astrophysical Journal Letters_ , 924(2), p. L32, URL https://dx.doi.org/10.3847/2041-8213/ac46a1.
  * Sano et al. (2004) Sano, T., Inutsuka, S.-i., Turner, N. J. and Stone, J. M. (2004), Angular Momentum Transport by Magnetohydrodynamic Turbulence in Accretion Disks: Gas Pressure Dependence of the Saturation Level of the Magnetorotational Instability, _The Astrophysical Journal_ , 605(1), pp. 321–339, astro-ph/0312480.
  * Sadowski and Narayan (2015) Sadowski, A. and Narayan, R. (2015), Photon-conserving Comptonization in simulations of accretion discs around black holes, _Monthly Notices of the Royal Astronomical Society_ , 454(3), pp. 2372–2380, 1508.04980.
  * Sadowski et al. (2014) Sadowski, A., Narayan, R., McKinney, J. C. and Tchekhovskoy, A. (2014), Numerical simulations of super-critical black hole accretion flows in general relativity, _Monthly Notices of the Royal Astronomical Society_ , 439(1), pp. 503–520, %****␣18_kar.bbl␣Line␣100␣****1311.5900.
  * Sadowski et al. (2015) Sadowski, A., Narayan, R., Tchekhovskoy, A., Abarca, D., Zhu, Y. and McKinney, J. C. (2015), Global simulations of axisymmetric radiative black hole accretion discs in general relativity with a mean-field magnetic dynamo, _Monthly Notices of the Royal Astronomical Society_ , 447(1), pp. 49–71, 1407.4421.
  * Sadowski et al. (2013) Sadowski, A., Narayan, R., Tchekhovskoy, A. and Zhu, Y. (2013), Semi-implicit scheme for treating radiation under M1 closure in general relativistic conservative fluid dynamics codes, _Monthly Notices of the Royal Astronomical Society_ , 429(4), pp. 3533–3550, 1212.5050.
  * Sadowski et al. (2017) Sadowski, A., Wielgus, M., Narayan, R., Abarca, D., McKinney, J. C. and Chael, A. (2017), Radiative, two-temperature simulations of low-luminosity black hole accretion flows in general relativity, _Monthly Notices of the Royal Astronomical Society_ , 466(1), pp. 705–725, 1605.03184.
  * Sorathia et al. (2012) Sorathia, K. A., Reynolds, C. S., Stone, J. M. and Beckwith, K. (2012), Global Simulations of Accretion Disks. I. Convergence and Comparisons with Local Models, _The Astrophysical Journal_ , 749(2), 189, 1106.4019.
  * White et al. (2019) White, C. J., Stone, J. M. and Quataert, E. (2019), A resolution study of magnetically arrested disks, _The Astrophysical Journal_ , 874(2), p. 168, URL https://dx.doi.org/10.3847/1538-4357/ab0c0c.