# Improved LARS algorithm for adaptive LASSO in the linear regression model

Manickavasagar Kayanan111Corresponding author: Department of Physical Science,
University of Vavuniya, Vavuniya, Sri Lanka. E-mail<EMAIL_ADDRESS>Department of Physical Science, University of Vavuniya, Vavuniya, Sri Lanka;
Pushpakanthie Wijekoon  National Institute of Fundamental Studies, Kandy, Sri
Lanka

###### Abstract

The adaptive LASSO has been used for consistent variable selection in place of
LASSO in the linear regression model. In this article, we propose a modified
LARS algorithm to combine adaptive LASSO with some biased estimators, namely
the Almost Unbiased Ridge Estimator (AURE), Liu Estimator (LE), Almost
Unbiased Liu Estimator (AULE), Principal Component Regression Estimator
(PCRE), r-k class estimator, and r-d class estimator. Furthermore, we examine
the performance of the proposed algorithm using a Monte Carlo simulation study
and real-world examples.

###### keywords:

Adaptive LASSO, LARS, Biased estimators, Monte Carlo simulation.

draft

##  1 Introduction

Let us consider a linear regression model

$\bm{y=X\beta+\epsilon}.$ (1)

Here, $\bm{y}$ represents the $n\times 1$ vector of observations on the
dependent variable, $\bm{X}$ is the $n\times p$ matrix of observations on the
non-stochastic predictor variables, $\bm{\beta}$ stands for a $p\times 1$
vector of unknown coefficients, and $\bm{\epsilon}$ denotes the $n\times 1$
vector of random error terms. These errors are assumed to be independent and
identically normally distributed with mean zero and common variance
$\sigma^{2}$.

It is widely acknowledged that the Ordinary Least Squares Estimator (OLSE)
serves as the Best Linear Unbiased Estimator (BLUE) for determining the
unknown parameter vector in model (1.1), expressed as:

$\begin{split}\hat{\bm{\beta}}_{OLSE}&=\arg\min_{\bm{\beta}}\left\\{(\bm{y-X\beta})^{\prime}(\bm{y-X\beta})\right\\}\\\
&=(\bm{X}^{\prime}\bm{X})^{-1}\bm{X}^{\prime}\bm{y}.\end{split}$ (2)

Nevertheless, the OLSE demonstrates instability and yields parameter estimates
with high variance in the presence of multicollinearity within $\bm{X}$. To
mitigate this multicollinearity issue, many researchers resort to biased
estimators.

As per Kayanan and Wijekoon (2017), the generalized representation of biased
estimators including Ridge Estimator (RE), Almost Unbiased RidgeEstimator
(AURE), Liu Estimator (LE), Almost Unbiased Liu Estimator (AULE), Principal
Component Regression Estimator (PCRE), r-k class estimator and r-d class
estimator can be expressed as:

$\hat{\bm{\beta}}_{G}=\bm{G}\hat{\bm{\beta}}_{OLSE}$ (3)

where

$\hat{\bm{\beta}}_{G}=\begin{cases}\hat{\bm{\beta}}_{RE}&\quad\text{if
}\bm{G}=(\bm{X}^{\prime}\bm{X}+k\bm{I})^{-1}\bm{X}^{\prime}\bm{X}\\\
\hat{\bm{\beta}}_{AURE}&\quad\text{if
}\bm{G}=\left(\bm{I}-k^{2}(\bm{X}^{\prime}\bm{X}+k\bm{I})^{-2}\right)\\\
\hat{\bm{\beta}}_{LE}&\quad\text{if
}\bm{G}=(\bm{X}^{\prime}\bm{X}+\bm{I})^{-1}(\bm{X}^{\prime}\bm{X}+d\bm{I})\\\
\hat{\bm{\beta}}_{AULE}&\quad\text{if
}\bm{G}=\left(\bm{I}-(1-d)^{2}(\bm{X}^{\prime}\bm{X}+\bm{I})^{-2}\right)\\\
\hat{\bm{\beta}}_{PCRE}&\quad\text{if }\bm{G}=\bm{T}_{h}\bm{T}_{h}^{\prime}\\\
\hat{\bm{\beta}}_{rk}&\quad\text{if
}\bm{G}=\bm{T}_{h}\bm{T}_{h}^{\prime}(\bm{X}^{\prime}\bm{X}+k\bm{I})^{-1}\bm{X}^{\prime}\bm{X}\\\
\hat{\bm{\beta}}_{rd}&\quad\text{if
}\bm{G}=\bm{T}_{h}\bm{T}_{h}^{\prime}(\bm{X}^{\prime}\bm{X}+\bm{I})^{-1}(\bm{X}^{\prime}\bm{X}+d\bm{I})\end{cases}$

Kayanan and Wijekoon (2017) demonstrated that the r-k class estimator and r-d
class estimator outperform other estimators within a specified range of
regularization parameter values when multicollinearity exists among the
predictor variables. However, biased estimators can introduce substantial bias
when the number of predictor variables is high, potentially leading to the
inclusion of irrelevant predictor variables in the final model. To address
this issue, Tibshirani (1996) proposed the Least Absolute Shrinkage and
Selection Operator (LASSO) as

$\hat{\bm{\beta}}_{LASSO}=\arg\min_{\bm{\beta}}\left\\{(\bm{y-X\beta})^{\prime}(\bm{y-X\beta})\right\\}\text{
subject to }\sum_{j=1}^{p}|\beta_{j}|\leq t,$ (4)

where $t\geq 0$ is a turning parameter. The LASSO solutions has been obtained
by the Least Angle Regression (LARS) algorithm.

According to Zou and Hastie (2005), LASSO failed to outperform Ridge Estimator
if high multicollinearity exists among predictors, and it is unsteady when the
number of predictors is higher than the number of observations. To overcome
this problem, Zou and Hastie (2005) proposed Elastic Net (ENet) estimator by
combining LASSO and RE as

$\hat{\bm{\beta}}_{Enet}=\arg\min_{\bm{\beta}}\left\\{(\bm{y-X\beta})^{\prime}(\bm{y-X\beta})+k\sum_{j=1}^{p}\beta^{2}_{j}\right\\}\text{
subject to }\sum_{j=1}^{p}|\beta_{j}|\leq t.$ (5)

The LARS-EN algorithm, which is a modified version of the LARS-LASSO
algorithm, has been used to obtain solutions for ENet.

Further, Zou and Hastie (2005) noted that LASSO does not care about variable
importance when a group of variables among which the pairwise correlations are
very high.

To handle this problem, Zou (2006) proposed adaptive LASSO by giving different
weights to regression coefficients in L1 penalty of LASSO. By taking weight
vector $\hat{\bm{w}}=|\hat{\bm{\beta}}_{OLSE}|^{-\alpha}$ for any $\alpha>0$,
the adaptive LASSO is defined as

$\hat{\bm{\beta}}_{adpLASSO}=\arg\min_{\bm{\beta}}\left\\{(\bm{y-X\beta})^{\prime}(\bm{y-X\beta})\right\\}\text{
subject to }\sum_{j=1}^{p}|w_{j}\beta_{j}|\leq t.$ (6)

In addition to that Zou and Zhang (2009) proposed adaptive Enet estimator by
combining adaptive LASSO and RE, and it is defined as

$\hat{\bm{\beta}}_{adpEnet}=\arg\min_{\bm{\beta}}\left\\{(\bm{y-X\beta})^{\prime}(\bm{y-X\beta})+k\sum_{j=1}^{p}\beta^{2}_{j}\right\\}\text{
subject to }\sum_{j=1}^{p}|w_{j}\beta_{j}|\leq t,$ (7)

where where $\hat{\bm{w}}=|\hat{\bm{\beta}}_{Enet}|^{-\alpha}$.

Kayanan and Wijekoon (2020) proposed the generalized version of LARS (GLARS)
algorithm that combines LASSO and biased estimators such as RE, AURE, LE,
AULE, PCRE,r-k class estimator and r-d class estimator. Further, they have
shown that the combination of LASSO and r-d class estimator performed well in
the high dimensional linear regression model when high multicollinearity exits
among the predictor variables.

In this article, we propose improved version of GLARS algorithm that can be
combined adaptive LASSO with other biased estimators such as AURE, LE, AULE,
PCRE, r-k class and r-d class estimators. Further, we compared the prediction
performance of the proposed algorithm with existing algorithms of adaptive
LASSO and adaptive ENet using a Monte-Carlo simulation study and a real-world
example.

The rest of the article is organized as follows: Section 2 presents the
proposed adaptive GLARS algorithm, Section 3 evaluates the performance of the
proposed algorithm, and Section 4 concludes the article.

##  2 Adaptive GLARS algorithm for LASSO

Based on the methodology outlined by Kayanan and Wijekoon (2020), we propose
the adaptive GLARS algorithm as follows:

Algorithm 1 Adaptive GLARS

1:Standardize the predictor variables $\bm{X}$ to have a mean of zero and a
standard deviation of one, and the response variable $\bm{y}$ to have a mean
of zero.

2:Define $\hat{\bm{w}}=\hat{|\bm{\beta}}_{G}|^{-\alpha}$ for $\alpha>0$, where
$\hat{\bm{\beta}}_{G}$ is the general form of the biased estimators defined in
equation (3), and $\bm{X}=\dfrac{\bm{X}}{\hat{\bm{w}}}$.

3:Initialize the estimated value of $\bm{\beta}$ as $\hat{\bm{\beta}}=0$, and
set the residual $\bm{r}_{0}=\bm{y}$.

4:Identify the predictor variable most correlated with $\bm{r}_{0}$ by:

  * •

Calculate $X_{j1}=\max_{j}|Cor(X_{j},\bm{r}_{0}|$ for $j=1,2,...,p$.

  * •

Increase the estimate of $\hat{\beta}_{j1}$ from 0 until another predictor
$X_{j2}$ has a high correlation with the current residual as $X_{j1}$ does.

  * •

Proceed in the equiangular direction between $X_{j1}$ and $X_{j2}$.

  * •

Similarly, each subsequent variable $X_{ji}$ earns its way into the active
set, and proceed in the equiangular direction between all selected predictors

  * •

Update coefficient estimates using the formula:

$\hat{\bm{\beta}}_{ji}=\hat{\bm{\beta}}_{j(i-1)}+\rho_{i}\bm{u}_{i},$ (8)

where $\alpha_{i}$ is a value between 0 and 1 representing the distance the
estimate moves before another variable enters the model, and $\bm{u}_{i}$ is
the equiangular vector.

  * •

Calculate the direction $\bm{u}_{i}$ using:

$\bm{u}_{i}=\bm{G}_{E}(\bm{E}^{\prime}_{i}\bm{X}^{\prime}\bm{X}\bm{E}_{i})^{-1}\bm{E}^{\prime}_{i}\bm{X}^{\prime}\bm{r}_{i-1},$
(9)

where $\bf Ei$ is the matrix with columns $(e_{j1},e_{j2},...,e_{ji})$,
$e_{j}$ is the $j$-th standard unit vector in $\mathbb{R}^{p}$ with the
indices of selected variables, and $\bm{G}_{E}$ depends on the specific
estimator which can be substituted by respective expressions for any of
estimators of our interest as listed in Table 1.

  * •

Update $\rho_{i}$ as:

$\rho_{i}=\min\left\\{\rho^{+}_{ji},\;\rho^{-}_{ji},\;\rho^{*}_{ji}\right\\}\in[0,1]$
(10)

where

$\rho^{\pm}_{ji}=\frac{Cor(\bm{r}_{i-1},X_{ji})\pm
Cor(\bm{r}_{i-1},X_{j})}{Cor(\bm{r}_{i-1},X_{ji})\pm
Cor(\bm{X}\bm{u}_{i},X_{j})}\text{for any }j\text{ such that
}\hat{\bm{\beta}}_{j(i-1)}=0,$ (11)

and

$\rho^{*}_{ji}=-\frac{\hat{\bm{\beta}}_{j(i-1)}}{\bm{u}_{i}}\text{ for any
}j\text{ such that }\hat{\bm{\beta}}_{j(i-1)}\neq 0.$ (12)

  * •

If $\rho_{i}=\rho^{*}_{ji}$, update $\bm{E}_{i}$ by removing the column
$e_{j}$ from $\bm{E}_{i-1}$. Calculate the new residual $\bm{r}_{i}$ as:

$\bm{r}_{i}=\bm{r}_{i-1}-\rho_{i}\bm{X}\bm{u}_{i},$ (13)

and move to the next step where $j_{i+1}$ is the value of $j$ such that
$\rho_{i}=\rho^{+}_{ji}$ or $\rho_{i}=\rho^{-}_{ji}$ or $\rho^{*}_{ji}$.

  * •

End this step when $\rho_{i}=1$.

5:Output $\hat{\bm{\beta}}_{adp}=\dfrac{\hat{\bm{\beta}}}{\hat{\bm{w}}}$.

Table 1: $\bm{G}_{E}$ of the estimators for GLARS Estimators | $\bm{G}_{E}$
---|---
OLSE | $\bm{E}_{i}$
RE | $\bm{E}_{i}(\bm{E}^{\prime}_{i}(\bm{X}^{\prime}\bm{X}+k\bm{I})\bm{E}_{i})^{-1}(\bm{E}^{\prime}_{i}\bm{X}^{\prime}\bm{X}\bm{E}_{i})$
AURE | $\bm{E}_{i}(\bm{I}_{p_{E}}-k^{2}(\bm{E}^{\prime}_{i}(\bm{X}^{\prime}\bm{X}+k\bm{I})\bm{E}_{i})^{-2})$
LE | $\bm{E}_{i}(\bm{E}^{\prime}_{i}(\bm{X}^{\prime}\bm{X}+\bm{I})\bm{E}_{i})^{-1}(\bm{E}^{\prime}_{i}(\bm{X}^{\prime}\bm{X}+d\bm{I})\bm{E}_{i})$
AULE | $\bm{E}_{i}(\bm{I}_{p_{E}}-(1-d)^{2}(\bm{E}^{\prime}_{i}(\bm{X}^{\prime}\bm{X}+\bm{I})\bm{E}_{i})^{-2})$
PCRE | $\bm{T}_{h_{E}}\bm{T}^{\prime}_{h_{E}}\bm{E}_{i}$
r-k class | $\bm{T}_{h_{E}}\bm{T}^{\prime}_{h_{E}}\bm{E}_{i}(\bm{E}^{\prime}_{i}(\bm{X}^{\prime}\bm{X}+k\bm{I})\bm{E}_{i})^{-1}(\bm{E}^{\prime}_{i}\bm{X}^{\prime}\bm{X}\bm{E}_{i})$
r-d class | $\bm{T}_{h_{E}}\bm{T}^{\prime}_{h_{E}}\bm{E}_{i}(\bm{E}^{\prime}_{i}(\bm{X}^{\prime}\bm{X}+\bm{I})\bm{E}_{i})^{-1}(\bm{E}^{\prime}_{i}(\bm{X}^{\prime}\bm{X}+d\bm{I})\bm{E}_{i})$

In Table 1, $\bm{I}_{p_{E}}$ is the $p_{E}\times p_{E}$ identity matrix,
$p_{E}$ is the number of selected variables in each subsequent step, and
$\bm{T}_{h_{E}}=(t_{1},t_{2},...,t_{h_{E}})$ is the first $h_{E}$ columns of
the standardized eigenvectors of
$\bm{E}^{\prime}_{i}\bm{X}^{\prime}\bm{X}\bm{E}_{i}$.

The adaptive GLARS algorithm iteratively updates combined estimates of
adaptive LASSO and other estimators. Evaluation of prediction performance
relies on the Root Mean Square Error (RMSE) criterion, elaborated in Section
3. Adaptive GLARS facilitates the integration of adaptive LASSO with any
estimator listed in Table 1.

It is worth noting that when $\bm{G}_{E}$ corresponds to the expressions of
OLSE and RE, adaptive GLARS provides solutions akin to adaptive LASSO and
adaptive ENet, respectively. For ease of reference, we denote adaptive GLARS
as adpLARS-LASSO, adpLARS-EN, adpLARS-AURE, adpLARS-LE, adpLARS-AULE, adpLARS-
PCRE, adpLARS-rk, and adpLARS-rd when $\bm{G}_{E}$ corresponds to the
expressions of OLSE, RE, AURE, LE, AULE, PCRE, r-k class, and r-d class
estimators, respectively.

We can use two-dimensional cross-validation to find the suitable value of
$\alpha$ and shrinkage parameter $k$ or $d$ for adaptive GLARS.

##  3 Performance of the adaptive GLARS algorithms

Proposed algorithms are compared using the RMSE criterion, which is the
expected prediction error of the algorithms, and is defined as

$RMSE(\hat{\bm{\beta}})=\sqrt{\frac{1}{n}({\bm{y}_{new}}-\bm{X}_{new}\hat{\bm{\beta}})^{\prime}({\bm{y}_{new}}-\bm{X}_{new}\hat{\bm{\beta}})}$
(14)

where $(\bm{y}_{new},\bm{X}_{new})$ denotes the new data which are not used to
obtain the parameter estimates, and $\hat{\bm{\beta}}$ is the estimated value
of $\bm{\beta}$ using the respective algorithm. A Monte Carlo simulation study
and the real-world examples are used for the comparison.

###  3.1 Simulation study

According to McDonald and Galarneau (1975), first we generate the predictor
variables by using the following formula:

$x_{i,j}=\sqrt{(1-\rho^{2})}z_{i,j}+\rho
z_{i,m+1}\qquad;i=1,2,...,n.\;\;j=1,2,...,m.$ (15)

where $z_{i,j}$ is an independent standard normal pseudo random number, and
$\rho$ is the theoretical correlation between any two explanatory variables.

In this study, we have used a linear regression model of 100 observations and
20 predictors. A dependent variable is generated by using the following
equation

$y_{i}=\beta_{1}x_{i,1}+\beta_{2}x_{i,2}+...+\beta_{5}x_{i,20}+\epsilon_{i}\qquad;i=1,2,...,100.$
(16)

where $\epsilon_{i}$ is a normal pseudo random number with mean zero and
common variance $\sigma^{2}$. We choose
$\bm{\beta}=(\beta_{1},\beta_{2},...,\beta_{20})$ as the normalized
eigenvector corresponding to the largest eigenvalue of $\bm{X}^{\prime}\bm{X}$
for which $\bm{\beta}^{\prime}\bm{\beta}=1$. To investigate the effects of
different degrees of multicollinearity on the estimators, we choose
$\rho=(0.5,\;0.7,\;0.9)$, which represents weak, moderated and high
multicollinearity. For the analysis, we have simulated 50 data sets consisting
of 50 observations to fit the model and 50 observations to calculate the RMSE.

The Cross-validated RMSE of the adaptive GLARS algorithms are displayed in
Fig. 1 \- Fig. 3, and the median cross-validated RMSE of the algorithms are
displayed in Table 2 \- Table 4.

Figure 1:  Cross-validated RMSE values of the adaptive GLARS algorithms when $\rho=0.5$.  Figure 2:  Cross-validated RMSE values of the adaptive GLARS algorithms when $\rho=0.7$.  Figure 3:  Cross-validated RMSE values of the adaptive GLARS algorithms when $\rho=0.9$.  Table 2:  Median Cross-validated RMSE values of the adaptive GLARS algorithms when $\rho=0.5$. Algorithms | RMSE | ($k,\;d$) | $\alpha$ | $t$ | Selected variables
---|---|---|---|---|---
adpLARS-LASSO | 3.45489 | – | 1 | 6.6635 | 16
adpLARS-EN | 3.41614 | 0.2 | 1 | 7.5795 | 17
adpLARS-AURE | 3.44668 | 1.0 | 1 | 7.1685 | 17
adpLARS-LE | 3.34648 | 0.3 | 1 | 7.1018 | 15
adpLARS-AULE | 3.48312 | 0.2 | 1 | 8.0718 | 16
adpLARS-PCRE | 3.31719 | – | 1 | 6.5019 | 16
adpLARS-rk | 3.35712 | 0.2 | 1 | 6.0726 | 17
adpLARS-rd | 3.47994 | 0.99 | 1 | 6.5019 | 16
Table 3:  Median Cross-validated RMSE values of the adaptive GLARS algorithms when $\rho=0.7$. Algorithms | RMSE | ($k,\;d$) | $\alpha$ | $t$ | Selected variables
---|---|---|---|---|---
adpLARS-LASSO | 3.53553 | – | 1 | 8.7067 | 16
adpLARS-EN | 3.42320 | 0.3 | 0.5 | 8.9330 | 17
adpLARS-AURE | 3.53440 | 0.7 | 1 | 8.0610 | 17
adpLARS-LE | 3.45469 | 0.1 | 0.5 | 9.1520 | 15
adpLARS-AULE | 3.56472 | 0.1 | 1 | 8.0821 | 16
adpLARS-PCRE | 3.41530 | – | 1 | 9.5873 | 16
adpLARS-rk | 3.35452 | 0.1 | 1 | 8.9412 | 16
adpLARS-rd | 3.37755 | 0.2 | 1 | 8.8207 | 16
Table 4:  Median Cross-validated RMSE values of the adaptive GLARS algorithms when $\rho=0.9$. Algorithms | RMSE | ($k,\;d$) | $\alpha$ | $t$ | Selected variables
---|---|---|---|---|---
adpLARS-LASSO | 3.44950 | – | 0.5 | 4.0460 | 15
adpLARS-EN | 3.39404 | 1.0 | 1 | 8.2710 | 17
adpLARS-AURE | 3.50448 | 0.9 | 1 | 10.045 | 17
adpLARS-LE | 3.49651 | 0.1 | 0.5 | 10.005 | 15
adpLARS-AULE | 3.48735 | 0.1 | 0.5 | 8.0684 | 16
adpLARS-PCRE | 3.49078 | – | 0.5 | 10.682 | 17
adpLARS-rk | 3.42176 | 0.3 | 1 | 10.433 | 16
adpLARS-rd | 3.37842 | 0.99 | 0.5 | 7.0576 | 15

Based on the insights gathered from Fig. 1 to Fig. 3 and Table 2 to Table 4,
it is evident that the adpLARS-PCRE, adpLARS-rk, and adpLARS-rd algorithms
consistently demonstrate superior performance in terms of RMSE criterion
compared to other adaptive GLARS algorithms across varying degrees of
multicollinearity, from weak to moderate and high levels, respectively.

###  3.2 Real-world example

In our analysis, we utilized the Prostate Cancer Data (Stamey et al., 1989), a
well-established dataset explored by Tibshirani (1996), Efron et al. (2004),
and Zou and Hastie (2005) to evaluate the efficacy of LASSO, LARS algorithm,
and Enet.

The Prostate Cancer Data comprises eight clinical metrics: log cancer volume
(lcavol), log prostate weight (lweight), age, log of benign prostatic
hyperplasia volume (lbph), seminal vesicle invasion (svi), log capsular
penetration (lcp), Gleason score (gleason), and percentage Gleason score 4 or
5 (pgg45). The response variable is the log of prostate-specific antigen
(lpsa), with a dataset size of 97 observations. Notably, the predictor
variables exhibit Variance Inflation Factor (VIF) values of 3.09, 2.97, 2.47,
2.05, 1.95, 1.37, 1.36, and 1.32, indicating considerable multicollinearity,
as evidenced by a high condition number of 243. This dataset is readily
available within the ”lasso2” R package. Our analysis involved fitting the
model with 67 observations and computing the Root Mean Square Error (RMSE)
using 30 observations.

Table 5:  Cross-validated RMSE values of Prostate Cancer Data using adaptive GLARS. Algorithms | RMSE | ($k,\;d$) | $\alpha$ | $t$ | Selected variables
---|---|---|---|---|---
adpLARS-LASSO | 0.77653 | – | 0.2 | 1.57112 | 7
adpLARS-EN | 0.78716 | 0.3 | 1 | 0.80638 | 7
adpLARS-AURE | 0.80638 | 1.0 | 0.9 | 0.80638 | 7
adpLARS-LE | 0.80014 | 0.1 | 0.5 | 1.45884 | 7
adpLARS-AULE | 0.79046 | 0.2 | 1 | 1.31322 | 6
adpLARS-PCRE | 0.76890 | – | 0.9 | 1.44929 | 7
adpLARS-rk | 0.77698 | 0.2 | 0.9 | 1.36273 | 7
adpLARS-rd | 0.76854 | 0.7 | 0.9 | 1.44764 | 7

The cross-validated RMSE values obtained through the adaptive GLARS algorithms
are summarized in Table 5. Upon examining Table 5, it becomes evident that the
adpLARS-rd algorithm outperforms other algorithms when applied to the Prostate
Cancer Data.

##  4 Conclusions

In conclusion, this study clearly shows that the adpLARS-rk and adpLARS-rd
algorithms work well for dealing with high dimensional linear regression
problems, especially when there are many closely related independent
variables. These algorithms emerge as reliable tools for tackling high
dimensional regression models and offer promising avenues for future research
and practical application in data-driven environments.

## Acknowledgement

The authors are grateful to the anonymous referee for a careful checking of
the details and for helpful comments that improved this paper.

## References

  * [2] [] Efron B, Hastie T, Johnstone I and Tibshirani R (2004). Least angle regression, The Annals of statistics, 32(2), 407–499.
  * [3]
  * [4] [] Kayanan M and Wijekoon P (2017). Performance of existing biased estimators and the respective predictors in a misspecified linear regression model, Open Journal of Statistics, 7(5), 876–900.
  * [5]
  * [6] [] Kayanan M and Wijekoon P (2020). Variable selection via biased estimators in the linear regression model, Open Journal of Statistics, 10, 113–126.
  * [7]
  * [8] [] McDonald GC and Galarneau DI (1975). A monte carlo evaluation of some ridge-type estimators, Journal of the American Statistical Association, 70(350), 407–416.
  * [9]
  * [10] [] Stamey TA, Kabalin JN, et al. (1989). Prostate specific antigen in the diagnosis and treatment of adenocarcinoma of the prostate: Ii. radical prostatectomy treated patients, Journal of Urology, 141(5), 1076–1083.
  * [11]
  * [12] [] Tibshirani R (1996). Regression shrinkage and selection via the lasso, Journal of the Royal Statistical Society: Series B (Methodological), 58(1), 267–288.
  * [13]
  * [14] [] Zou H (2006). The adaptive lasso and its oracle properties, Journal of the American Statistical Association, 101(476), 1418–1429.
  * [15]
  * [16] [] Zou H and Hastie T (2005). Regularization and variable selection via the elastic net, Journal of the Royal Statistical Society: Series B, 67(2), 301–320.
  * [17]
  * [18] [] Zou H and Zhang HH (2009). On the adaptive elastic-net with a diverging number of parameters, The Annals of Statistics, 37(4), 1733–1751


# Magnetic Flux Transport in Radiatively Inefficient Accretion Flows and the
Pathway towards a Magnetically Arrested Disk

Prasun Dhang Institute for Advanced Study, Tsinghua University, Beijing
100084, China  IUCAA, Post Bag 4, Ganeshkhind, Pune, Maharashtra 411007, India
Xue-Ning Bai Institute for Advanced Study, Tsinghua University, Beijing
100084, China  Department of Astronomy, Tsinghua University, Beijing 100084,
China  Christopher J. White  Department of Astrophysical Sciences, Princeton
University, Princeton, NJ 08544

###### Abstract

Large-scale magnetic fields play a vital role in determining the angular
momentum transport and in generating jets/outflows in the accreting systems,
yet their origin remains poorly understood. We focus on radiatively
inefficient accretion flows (RIAF) around the black holes, and conduct three-
dimensional general-relativistic magnetohydrodynamic (GRMHD) simulations using
the Athena++ code. We first re-confirm that the dynamo action alone cannot
provide sufficient magnetic flux required to produce a strong jet. We next
investigate the other possibility, where the large-scale magnetic fields are
advected inward from external sources (e.g. the companion star in X-ray
binaries, magnetized ambient medium in AGNs). Although the actual
configuration of the external fields could be complex and uncertain, they are
likely to be closed. As a first study, we treat them as closed field loops of
different sizes, shapes and field strengths. Unlike earlier studies of flux
transport, where magnetic flux is injected in the initial laminar flow, we
injected the magnetic field loops in the quasi-stationary turbulent RIAF in
inflow equilibrium and followed their evolution. We found that a substantial
fraction ($\sim 15\%-40\%$) of the flux injected at the large radii reaches
the black hole with a weak dependence on the loop parameters except when the
loops are injected at high latitudes, away from the mid-plane. Relatively high
efficiency of flux transport observed in our study hints that a magnetically
dominated RIAF, potentially a magnetically-arrested disk, might be formed
relatively easily close to the black hole, provided that a source of the
large-scale field exists at the larger radii.

Accretion — Accretion disk — GRMHD

##  1 Introduction

Astrophysical accretion disks influence the systems over a large range of
scales spanning from the planet formation to the galaxy evolution. They also
energise the most powerful sources in the Universe. For example, disks
orbiting around the stellar-mass black holes (BHs) and neutron stars are
considered to be among the most luminous X-ray sources in the sky (Remillard &
McClintock 2006). Active galactic nuclei (AGNs), powered by the accretion of
matter onto a supermassive black hole at the centre of galaxies, are not only
the most powerful sources, but also energy released by the AGNs provide
feedback to the entire galaxy and determines its evolution (Silk & Rees 1998;
Harrison 2017; Morganti 2017).

Broadly speaking, accretion occurs via three different modes: (i)
geometrically thin and optically thick Keplerian disk (standard disk; Shakura
& Sunyaev 1973; Novikov & Thorne 1973), (ii) geometrically thick and optically
thin radiatively inefficient accretion flows (RIAF; Chakrabarti 1989; Narayan
& Yi 1994; Blandford & Begelman 1999) and (iii) geometrically and optically
thick slim disks (Abramowicz et al. 1988). Slim disks accrete matter at super-
Eddington rate, while mass accretion rate $\dot{m}$ is sub-Eddington in both
the standard disk ($10^{-4}\lesssim\dot{m}/\dot{M}_{\rm Edd}\lesssim 1$) and
in RIAFs ($\dot{m}/\dot{M}_{\rm Edd}\lesssim 10^{-4}$). In this paper, we
focus on the RIAFs, where disks likely span most of their time (e.g., Yuan &
Narayan 2014, such as disks around Sgr A∗ and the SMBH in M87), and their
dynamics is relatively simple compared to the other two states.

The structure and evolution of the rotationally supported accretion disks are
primarily determined by the process of angular momentum transport. The current
consensus is that the magnetorotational instability (MRI; Balbus & Hawley
1991) gives rise to angular momentum transport and vigorous turbulence in a
fully ionised accretion flow (e.g in X-ray binaries, inner part of the AGN
disk, sufficiently ionised part of proto-planetary disks). The MRI becomes
more efficient in angular momentum transport if the accretion disk is threaded
by a net vertical magnetic flux. This has been observed in the the local
shearing-box simulations that the presence of a net vertical magnetic flux
enhances the MRI turbulence and hence angular momentum transport (Bai & Stone
2013). Additionally, net flux threading the disk helps to launch the
winds/outflows (Bai & Stone 2013; Suzuki & Inutsuka 2014).

Large-scale magnetic field close to the central accretor (a BH or a neutron
star) is a necessary ingredient for jet production in accreting systems
(Blandford & Znajek 1977; Blandford & Payne 1982). It has been proposed that a
RIAF saturated with strong ploidal magnetic flux close to the BH provides an
ideal condition for the jet production (Bisnovatyi-Kogan & Ruzmaikin 1974;
Esin et al. 1997; Fender et al. 1999; Narayan et al. 2003; Meier 2005). The
idea has been verified in numerical simulations (Igumenshchev et al. 2003;
Narayan et al. 2012). The studies found that a strongly magnetized RIAF,
namely magnetically arrested disk (MAD) around spinning BH produces strong
jets, extracting net energy from the BH spin via the Penrose-Blandford-Znajek
process (Tchekhovskoy et al. 2011; McKinney et al. 2012).

MAD model predicts a correlation among the mass accretion rate, magnetic flux
threading the BH and jet power which is found to be in agreement with the
observations of radio-loud AGNs (Zamaninasab et al. 2014; Ghisellini et al.
2014). Recent polarization studies of M87 at 230 GHz from Even Horizon
Telescope EHT observations (Event Horizon Telescope Collaboration et al.
2021a, b; Yuan et al. 2022) also infer the presence of a dynamically important
near-horizon organized, poloidal magnetic flux consistent with GRMHD models of
MAD.

What could be the possible source of the magnetic flux close the BH? Most of
the numerical simulations of MAD start with a strong enough large-scale
poloidal flux which is eventually brought close to the BH and get accumulated
by flux-freezing (Tchekhovskoy et al. 2011; McKinney et al. 2012). However,
the source of the large-scale field is not entirely obvious. It can
potentially be generated in the disk itself by a dynamo action, or be advected
in from some external sources.

Efficiency of the dynamo action in generating coherent and strong large-scale
poloidal field required to produce strong jets is found to be different for
different numerical simulations. The $\alpha-$ effect (responsible for
generation of poloidal field in dynamo) is found to be weak in simulations
that start with small poloidal magnetic loops (Hogg & Reynolds 2018a; Dhang &
Sharma 2019; Dhang et al. 2020). The quasi-stationary states of those
simulations are in a weakly magnetized regime, popularly know as “standard and
normal evolution” (SANE; Narayan et al. 2012). However, recent simulations
with a very strong (with gas to magnetic pressure ratio $\beta\approx 5$) and
coherent initial toroidal field showed the production of large-scale poloidal
field loops of the size of scale-height $H\propto R$ and led to the MAD
eventually (Liska et al. 2020). Therefore, it is worth noting that the
simulations need to start with a strong and coherent large-scale field (either
poloidal or toroidal) to achieve MAD.

In addition to the in-situ generations of the magnetic field by a dynamo
process, it might be possible that an initially weak field supplied to the
disk (from the outer part of the disk or companion star in case of XRBs,
ambient medium of the AGNs) can in principle be amplified by flux freezing.
Flux accumulation near the BH depends on the relative efficiency between the
inward advection by the accretion flow and the outward diffusion due to a
turbulent resistivity (Lubow et al. 1994). Additionally, turbulent pumping can
also cause outward transport of the large scale magnetic field in a dynamo-
active accretion flow (Dhang et al. 2020). However, a few studies proposed
that vertical magnetic field accretion can be efficient in the hot tenuous
surface layer (coronal region, where radial velocity is comparatively higher
compared to that in the mid- plane) in a hot accretion flow (Beckwith et al.
2009). It is also interesting to note that the simulations of large-scale
accretion flow around the galactic centre fed by the magnetized winds of Wolf-
Rayet stars also show efficient inward transport of magnetic field towards the
centre (Ressler et al. 2020a, b).

This paper studies the magnetic flux transport in a fully turbulent RIAF,
unlike the previous studies where magnetic flux is injected in the initial
laminar condition. Therefore, first, we run a simulation to attain a quasi-
stationary RIAF in the SANE regime (weakly magnetized). Then we inject the
external magnetic flux on top of the existing magnetic field in this turbulent
SANE RIAF. In the later part of the paper, we refer it as the Initial RIAF
run. It is customary to use net vertical magnetic flux threading the disk to
investigate the flux transport (Beckwith et al. 2009; Zhu & Stone 2018; Mishra
et al. 2019. However, we argue that the geometry of the external magnetic
field is likely to be closed. In this paper, as a first step, we huse magnetic
field loops as the simplest possible form of the external magnetic field,
studying its transport in the turbulent RIAF and its possibility of saturating
the BH with magnetic flux towards the MAD regime.

The paper is organized as follows. In Section 2, we discuss the solution
method and physical set-up of the RIAF simulations. In Section 3, we discuss
the evolution of the flow, convergence and magnetic state of the Initial RIAF
run. We describe the method of flux injection and its results in section 4.
Finally, the key points of results are discussed and summarized in Sections 5
and 6.

##  2 Method

###  2.1 Equations solved

We solve the ideal general relativistic magnetohydrodyamic (GRMHD) equations

$\displaystyle\partial_{0}\left(\sqrt{-g}\rho
u^{0}\right)+\partial_{j}\left(\sqrt{-g}\rho u^{j}\right)=0$ (1)
$\displaystyle\partial_{0}\left(\sqrt{-g}T^{0}_{\mu}\right)+\partial_{j}\left(\sqrt{-g}T^{j}_{\mu}\right)=\frac{1}{2}\sqrt{-g}T^{\nu\sigma}\partial_{\mu}g_{\nu\sigma}$
(2)
$\displaystyle\partial_{0}\left(\sqrt{-g}B^{i}\right)+\partial_{j}\left(\sqrt{-g}\
F^{*ij}\right)=0$ (3)
$\displaystyle\frac{1}{\sqrt{-g}}\partial_{i}\left(\sqrt{-g}B^{i}\right)=0$
(4)

in a spherical-like Kerr-Schild coordinates ($t,\ r,\ \theta,\ \phi$) with
$G=c=M_{BH}=1$. All the length scales, time scales in this work are expressed
in units of the gravitational radius $r_{g}=GM_{BH}/c^{2}$ and $t_{g}=r_{g}/c$
respectively, unless stated otherwise. Here, $g_{\mu\nu}$ and $g$ are metric
coefficients and metric determinant respectively. Following convention, Greek
indices run through [0,1,2,3], while $i$ denotes a spatial index. Equations 1,
2, 3 & 4 describe conservation of particle number, conservation of energy-
momentum, source-free Maxwell equations and no-magnetic monopole constraint
respectively. Here

$T^{\mu\nu}=\left(\rho h+b^{2}\right)u^{\mu}u^{\nu}+\left(p_{\rm
gas}+\frac{b^{2}}{2}\right)g^{\mu\nu}-b^{\mu}b^{\nu}$ (5)

is the stress-energy tensor and

$F^{*\mu\nu}=b^{\mu}u^{\nu}-b^{\nu}u^{\mu}$ (6)

is the dual of the electromagnetic field tensor given that $\rho$ is the
comoving rest mass density, $p_{\rm gas}$ is the comoving gas pressure,
$u^{\mu}$ is the coordinate frame 4-velocity, $\Gamma=5/3$ is the adiabatic
index of the gas, $h=1+\Gamma/(\Gamma-1)p_{\rm gas}/\rho$ is comoving enthalpy
per unit mass, $B^{i}=F^{*i0}$ is the magnetic field in the coordinate frame.
Four magnetic field $b^{\mu}$ are related to 3-magnetic field $B^{i}$ as

$\displaystyle b^{0}=g_{i\mu}B^{i}u^{\mu},$ (7) $\displaystyle
b^{i}=\frac{B^{i}+b^{0}u^{i}}{u^{0}}.$ (8)

For diagnostics, we also use magnetic field components ($B_{r},\ B_{\theta},\
B_{\phi}$) defined in a spherical-polar like quasi-orthonormal frame as

$B_{r}=B^{1},\ B_{\theta}=rB^{2},\ B_{\phi}=r\sin\theta B^{3}.$ (9)

We use the GRMHD code Athena++ (White et al., 2016) to perform the
simulations. We employ the HLLE solver (Einfeldt, 1988) with a third-order
piecewise parabolic method (PPM; Colella & Woodward (1984)) for spatial
reconstruction. For time integration, a second-order accurate van Leer
integrator is used with the CFL number 0.3. We use a CT (Gardiner & Stone,
2005; White et al., 2016) update of the face-centered magnetic fields to
maintain the no magnetic monopole condition.

###  2.2 Initial Condition

Figure 1: Pictorial description of the grids and initial hydrostatic and
magnetic conditions of the simulation. Each block represents a meshblock of
size $16\times 14$ in poloidal plane in the $r$, $\theta$ directions
respectively. Each meshblock has 16 grid points in the azimuthal direction.
Disk aspect ratio of the initial disk is $H/R=0.23$ and shown by the white
dashed lines. Streamlines describe the initial poloidal magnetic field lines
with average plasma $\beta\approx 800$.

We initialise a geometrically semi-thick disk of aspect ratio
$\epsilon^{\text{*}}_{\rm in}=H_{G}/R=0.23$ embedded in a hot corona. Here
$H_{G}$ is the Gaussian scale-height. The rest mass density distribution of
the initial disk is given by

$\rho_{d}(r,\theta)=e^{-z^{2}/2H_{G}^{2}}\left(\frac{r_{\rm
in}}{R}\right)^{q_{d}}\delta(r_{\rm in});$ (10)

and the gas pressure is given by

$p_{\rm gas,d}=\rho_{d}c^{2}_{sd}=\rho_{d}\ \epsilon^{\text{*}2}_{\rm
in}\left(\frac{M_{\rm BH}}{R}\right).$ (11)

Here, $z=r{\rm cos}\theta$, $R=r{\rm sin}\theta$ and $\delta(r_{\rm
in})=1/[1+e^{-(R-r_{\rm in})/5\epsilon^{\text{*}}_{\rm in}H_{G}}]$ is a
tapping function with an inner disk radius $r_{\rm in}=15$. We consider
$q_{d}=-1.5$ and mass of the BH to be $M_{\rm BH}=1$. It is to be noted that
the Gaussian scale height $H_{G}$ is related to the density-weighted scale
height

$H=\frac{\int\sqrt{-g}\ \rho|\frac{\pi}{2}-\theta|\ d\theta\
d\phi}{\int\sqrt{-g}\ \rho\ d\theta\ d\phi}$ (12)

as $H=\sqrt{2/\pi}H_{G}$ and hence, disk aspect ratio
$\epsilon=\sqrt{2/\pi}\epsilon^{\text{*}}_{\rm in}$.

The disk is surrounded by an atmosphere defied by

$\displaystyle\rho_{c}=\rho_{\rm in}\left(\frac{r_{\rm
in}}{r}\right)^{q_{c}};~{}~{}~{}~{}p_{\rm gas,c}=\rho_{c}\frac{M_{BH}}{r}$
(13)

with $q_{c}=-1.5$ and $\rho_{\rm in}=10^{-5}$. The tenuous atmosphere is
static while the gas within the disk is rotating with a Keplerian speed given
by

$u^{3}=\frac{r}{r-2}R^{-3/2}$ (14)

in the Boyer-Lindquist coordinates.

In order to attain a quasi-stationary weakly magnetized RIAF (SANE; Narayan et
al. 2012), we initialise the multiple magnetic field loops using the vector
potential (Penna et al. 2013)

$A_{\phi}=\begin{cases}Q~{}\sin\left[f(r)-f(r_{\rm in})\right],~{}~{}Q>0\\\
0,{\rm otherwise}.\end{cases}$ (15)

Here,

$\displaystyle Q=C_{B}~{}{\rm sin}^{3}\theta\left(\frac{p_{1}}{p_{2}}-p_{\rm
cut}\right),$ (16) $\displaystyle
f(r)=\left(r^{2/3}+\frac{15}{8r^{2/5}}\right)\frac{1}{\lambda_{B}};$ (17)

with $p_{1}(r,\theta)=p_{\rm gas}(r,\theta)-p_{\rm gas}(r_{B0},\pi/2)$,
$p_{2}(r)=p_{\rm gas}(r,\pi/2)-p_{\rm gas}(r_{B0},\pi/2)$. The vector
potential $A_{\phi}$ vanishes for $r>r_{B0}=200$. We choose $C_{B}=0.5$,
$p_{\rm cut}=0.4$ and $\lambda_{B}=0.75$ giving rise to an average plasma
$\beta=800$ for the initial disk (averaging is done over the region within one
scale-height of the disk).

###  2.3 Numerical setup

We carry out all the simulations of RIAFs around a non-spinning BH (spin
parameter $a=0$). The computational domain spans over $r\in[1.94,300]$,
$\theta\in[0,\pi]$, $\phi\in[0,2\pi/3]$. It is to be noted that one grid point
is inside the event horizon in the radial direction at the root level. This
allows to obtain a causally disconnected inner boundary. Radial grids are
spaced logarithmically, while meridional grids are compressed towards the mid-
plane using

$\theta=\theta_{u}+\frac{1-s}{2}\ \sin 2\theta_{u}$ (18)

with $s=0.49$ which gives rise to $\Delta\theta_{\rm pole}/\Delta\theta_{\rm
eq}\approx 3.0$. Uniform grids are employed in the azimuthal direction. To
improve the effective resolution, we use two levels of static refinements with
a root grid resolution $160\times 56\times 32$ giving rise to $\Delta
r:r\Delta\theta:r\Delta\phi=1.2:1:2.3$ at the equator in the Newtonian limit.
While first level of refinement covers
$r_{L1}\in[1.94,180],\theta_{L1}\in[\pi/3,2\pi/3]$, a second level of
refinement is applied to the region $8<r<140$, $4\pi/9<\theta<5\pi/9$ such
that the number of $\theta$-cells per scale-height is $H/r\Delta\theta\approx
40$ in the quasi-stationary state.

We use a pure inflow boundary condition ($u^{1}\leq 0$) at the radial inner
boundary, while at the radial outer boundary, primitive variables are set
according to their initial radial gradients. Magnetic fields in the inner
ghost zones are copied from the nearest computation zone. On the other hand,
magnetic fields at the outer ghost zones are set according to
$B_{r},B_{\phi}\propto r^{-2}$ while keeping $B_{\theta}$ unchanged from the
last computation zone. Polar and periodic boundary conditions are used at the
meridional and azimuthal boundaries respectively.

We will introduce various diagnostics as we discuss simulation results, where
the data will be averaged in different ways. Here, for future reference, we
mention that the symbol ‘$\ \bar{}\ $’ is reserved for the azimuthally
averaged mean quantities, while any additional averaging (e.g. vertical or
time averaging) of the quantities will be indicated by $\langle.\rangle$ in
this paper.

##  3 Evolution of the Initial RIAF Run

To begin with, we would like to investigate the plausibility of the conversion
of a weakly magnetized RIAF (SANE) into a highly magnetized one (MAD) due to a
dynamo action. Therefore, we perform an Initial RIAF simulation as previously
mentioned. We run the simulation for the time $t=8\times 10^{4}$ to probe
whether a SANE to MAD conversion occurs. In this section, we describe the
evolution of the Initial RIAF towards the stationarity, its convergence and
magnetic state.

###  3.1 Flow evolution of the Initial RIAF

Figure 2: Evolution of the mean ($\phi$-averaged) poloidal
($\bar{\mathbf{B}}_{p}=\bar{\mathbf{B}}_{r}+\bar{\mathbf{B}}_{\theta}$) and
toroidal ($\bar{\mathbf{B}}_{\phi}$) magnetic fields with time in the Initial
RIAF run. The colour shows the mean toroidal field, while streamlines describe
the mean poloidal field lines. Poloidal field loops of alternate polarity
thread the initial disk. Shortly afterwards, toroidal magnetic field becomes
the dominant component due to background shear. The accretion flow largely
forgets the initial field configuration around $t=40200$, entering a quasi-
stationary phase with an MRI-driven dynamo in action.

Fig. 2 shows the evolution of the flow with time for the Initial RIAF run;
colors describe mean toroidal field ($\phi$-averaged, for definition see
equation 25), and streamlines depict mean poloidal fields
$\bar{\mathbf{B}}_{p}=\bar{\mathbf{B}}_{r}+\bar{\mathbf{B}}_{\theta}$. The
first panel of Fig. 2 shows the magnetic initial condition- poloidal field
loops of alternate signs, aiming to achieve a weakly magnetized RIAF (SANE;
Narayan et al. 2012) in the quasi-stationary phase. Shear in the accretion
flow converts poloidal field into the toroidal field, while MRI amplifies the
poloidal field. Therefore, both poloidal and toroidal fields grow
exponentially in a dynamical time ($t_{\rm dyn}\approx 1/\Omega\propto
R^{3/2}$) and after few dynamical time the system likely enters the non-linear
regime under the influence of parasitic instabilities (Goodman & Xu 1994) or
due to different super-Alfvénic rotational instabilities (SARIs; Goedbloed &
Keppens 2022).

Second panel of Fig. 2 shows the time $t=10050$ when MHD turbulence is fully
developed throughout the region of interest ($r\leq 120$). However, it is
worth noting that the system still remembers the initial field geometry as
indicated by the alternate signs of mean toroidal fields at different radii.
As time evolves, alternate polarity fields reconnect and the accretion flow
gradually removes the signature of initial field geometry as can be seen in
the third panel of Fig. 2. Around the time $t=40200$ (fourth panel of Fig. 2),
the accretion flow largely forgets its magnetic initial condition and the
magnetic fields generated due to an in-situ dynamo start to dominate. Finally,
the subsequent disk evolution is self-regulated with a combination of the MRI
turbulence, dynamo and angular momentum transport.

Magnetic field generated by dynamo in the quasi-stationary phase of RIAF are
not only of large-scales in the azimuthal direction (Dhang & Sharma, 2019),
but also in the radial direction, as can easily be seen by the radially
extended structures of mean poloidal and toroidal fields in the last panels of
of Fig. 2. However, the strength of the dynamo generated large-scale field is
insufficient to form a MAD, as we will discuss in the section 3.3.

###  3.2 Convergence

Figure 3: Time averaged radial profiles of the quality factors ($Q_{\theta}$
and $Q_{\phi}$) and the magnetic tilt angle $\theta_{B}$ in the mid-plane for
the Initial RIAF run. Time average is done over $t=(5-10)\times 10^{4}$. The
simulation is well resolved for $r\leq 120$ and resolvability starts declining
afterward.

Before analyzing simulations results, we first verify that our simulations
have achieved proper numerical convergence. In doing so, we calculate
different numerical metrics which were found to be useful in defining the
convergence of the MRI turbulence. In this work, we focus on the quality
factors

$\displaystyle
Q_{\theta}=\frac{2\pi}{\Omega}\frac{|\bar{b}^{\hat{\theta}}|}{\sqrt{{\overline{w}_{\rm
tot}}}}\frac{1}{dx^{\hat{\theta}}},$ (19) $\displaystyle
Q_{\phi}=\frac{2\pi}{\Omega}\frac{|\bar{b}^{\hat{\phi}}|}{\sqrt{{\overline{w}_{\rm
tot}}}}\frac{1}{dx^{\hat{\phi}}}$ (20)

and the magnetic tilt angle

$\theta_{B}=-\frac{\overline{b^{\hat{r}}b^{\hat{\phi}}}}{\overline{p}_{\rm
mag}}$ (21)

measured in an orthonormal fluid frame (White et al. 2019). Here the angular
velocity is defined as $\Omega(r,\theta)=\bar{u}^{3}/\bar{u}^{0}$, and total
entropy is given by $\overline{w}_{\rm
tot}(r,\theta)=\overline{\rho+\Gamma/(\Gamma-1)p_{\rm gas}+p_{\rm mag}}$. Line
elements are given by
$dx^{\hat{\theta}}=g_{\mu\nu}e^{\mu}_{\hat{\theta}}dx^{\nu}_{BL}$,
$dx^{\hat{\phi}}=g_{\mu\nu}e^{\mu}_{\hat{\phi}}dx^{\nu}_{BL}$, where
$dx^{\mu}_{BL}=\left[0,~{}\Delta r,~{}\Delta\theta,~{}\Delta\phi-
ar/(r^{2}-2Mr+a^{2})\Delta\right]$. Quality factors $Q_{\theta}$ and
$Q_{\phi}$ provide the information of number of cells across an wavelength of
the fastest growing mode in the $\theta$ and $\phi$\- directions respectively;
while $\theta_{B}$ measures the magnetic field anisotropy, a key factor behind
angular momentum transport. Earlier studies suggested that the toroidal and
poloidal resolutions are coupled and the product of the quality factors
$Q_{\theta}Q_{\phi}\geq 200-250$ is a good indicator for convergence in the
MRI simulations (Sorathia et al. 2012; Narayan et al. 2012; Dhang & Sharma
2019; Porth et al. 2019). In the meantime, we note that $\theta_{B}$ shows a
narrow range of value $10^{\circ}-14^{\circ}$ for the converged runs (e.g
Sorathia et al. 2012; Hogg & Reynolds 2018b; Dhang & Sharma 2019) and turned
out to be a better indicator of convergence.

The top and bottom panels of Fig. 3 show the radial profiles of the average
(averaged over $\phi$, $\theta$ and time) quality factors ($\langle
Q_{\theta}\rangle$, $\langle Q_{\phi}\rangle$), and magnetic tilt angle
($\langle\theta_{B}\rangle$) close to the mid-plane of the disk for the
Initial RIAF run. Meridional average is done over one scale-height above and
below the mid-plane, while azimuthal average is done over all cells. Time
average is done over the time interval $t=(5-10)\times 10^{4}$. While the
quality factors indicate that our simulation is marginally resolved with
$\langle Q_{\theta}\rangle\langle Q_{\phi}\rangle\gtrsim 180$ up to $r=100$,
the radial profile of $\langle\theta_{B}\rangle$ clearly shows that our
Initial RIAF simulation is well resolved till $r=120$ and resolvability starts
to decline afterwards because of the poor resolution at larger radii.

###  3.3 Characterizing the magnetic state of the Initial RIAF?

Figure 4: The time history of accreted specific angular momentum $j_{\rm net}$
and MAD parameter $\phi_{\rm BH}$ calculated at the event horizon and at the
ISCO for the RIAF run. For correspondence, time variation of the magnetic flux
$\Phi_{\rm NH}$ (equation 27) threading the BH is also shown.

In this section, we characterize our Initial RIAF simulation, particularly
examining the indicators that distinguish the SANE from the MAD state.
Following, Narayan et al. (2012), we study the time evolution of the specific
angular momentum of the accreting material

$j_{\rm net}(r,t)=-\frac{1}{\dot{m}}\int T^{1}_{3}\ dS_{r}$ (22)

and the MAD parameter

$\phi_{BH}(r,t)=\frac{\sqrt{4\pi}}{2\sqrt{\dot{m}}}\int|B^{1}|\ dS_{r}$ (23)

to investigate the magnetic state of the Initial RIAF run. Here, the mass
accretion rate at any radius $r$ is defined as

$\dot{m}(r,t)=-\int\rho u^{1}\ dS_{r},$ (24)

where area element is given by $dS_{r}=\sqrt{-g}\ d\theta\ d\phi$.

The MAD parameter $\phi_{BH}$ is a dimensionless number which is found to be
useful in characterizing the magnetic state of the simulations. Earlier
studies suggest that an accretion flow attains a MAD state once $\phi_{BH}$
reaches a critical value $\phi_{BH,c}\approx 40$ at the event horizon
(Tchekhovskoy et al. 2011). Additionally, $j_{\rm net}$ also shows a highly
sub-Keplerian nature at the event horizon for MAD simulations, where angular
momentum transport is highly efficient due to the large-scale Maxwell stress.
On the contrary, $j_{\rm net}$ maintains a slightly sub-Keplerian value at the
event horizon in the SANE simulations (Narayan et al. 2012). Moreover, $j_{\rm
net}$ is shown to be a good indicator of convergence in the MRI-active
turbulent accretion flow. For a converged simulation, $j_{\rm net}$ maintains
a sub-Keplerian value inside the ISCO with a non-decreasing trend in time
throughout the simulation (Hawley et al. 2013; Dhang & Sharma 2019).

Fig. 4 show the time evolution of $j_{\rm net}$ and $\phi_{BH}$ at two
different radii, at the ISCO and at the event horizon. We also plot the time
variation of signed flux threading the event horizon ($r_{H}=2$) of the BH in
northern hemisphere, $\Phi_{NH}(r_{H})$ (equation 27) for a future reference
in section 4. The value of $\phi_{\rm BH}$ always remains around one which is
well below the value ($\geq 40$) required for the MAD state. Such a low value
of $\phi_{\rm BH}$ implies that the magnetic state of the Initial RIAF run is
in the SANE regime. Slightly sub-Keplerian value of specific angular momentum
$j_{\rm net}$ is also an indicator of the SANE magnetic state of our initial
RIAF simulation.

###  3.4 Inflow equilibrium

Figure 5: Top panel: Time and spatial ($\theta$ and $\phi$) averaged accretion
rates as a function of radius for the Initial RIAF run. Time averages are done
over five different intervals $\Delta t_{0}=(3.75-7.5)\times 10^{3}$, $\Delta
t_{1}=(7.5-15)\times 10^{3}$, $\Delta t_{2}=(1.5-3)\times 10^{4}$, $\Delta
t_{3}=(3-6)\times 10^{4}$ and $\Delta t_{4}=(6-12)\times 10^{4}$. It can
roughly be inferred that the inflow equilibrium radius is $r_{\rm eq}=60$ for
the Initial RIAF run.
Bottom panel: disk aspect ratios $\epsilon=H/r$ for the RIAF run in the quasi-
stationary state. Simultaneously, we also plot its initial radial profile to
study the change over time. Time averages are done over $\Delta
t_{4}=(6-12)\times 10^{4}$.

We will inject external magnetic flux in the quasi-stationary turbulent RIAF
to study magnetic flux transport (ref. section 4). Therefore, it is important
to find out the inflow equilibrium radius - radius within which flow attains a
quasi-stationary state, for the Initial RIAF run. Following Narayan et al.
(2012), we investigate the variation of average mass accretion rate
$\langle\dot{m}(r)\rangle$ with time to find out the inflow equilibrium
radius. Spatial averages are done over all $\theta$ and $\phi$. We use five
different intervals $\Delta t_{0}=(3.75-7.5)\times 10^{3}$, $\Delta
t_{1}=(7.5-15)\times 10^{3}$, $\Delta t_{2}=(1.5-3)\times 10^{4}$, $\Delta
t_{3}=(3-6)\times 10^{4}$ and $\Delta t_{4}=(6-12)\times 10^{4}$ to do the
time average. The top panel of Fig. 5 shows $\langle\dot{m}(r)\rangle$ at
different time intervals for the Initial RIAF run. It can be inferred from the
radial profiles of $\dot{m}$ that the inflow equilibrium radius for the
Initial RIAF run reaches $r_{\rm eq}\approx 60$ at late times. This is the
radius that guides us to determine the injection radius for the external
magnetic field loops.

The bottom panel of Fig. 5 shows the radial variation of disk aspect ratio
$\epsilon=H/r$ in the quasi-stationary state. The disk aspect ratio $\epsilon$
slowly increases with increasing radius until the inflow equilibrium radius
and its value lies around $\epsilon=0.25$ for $r\gtrsim 20$, where general
relativistic effects are negligible. Such a variation of scale height in our
simulation is also in agreement with that observed in previous GRMHD
simulations of the SANE RIAF (e.g Narayan et al. 2012).

###  3.5 Large-scale magnetic field and dynamo

Figure 6: The butterfly diagram: Space-time plots of mean radial
$B_{r}(r=60,,\theta,t)$ (top panel) and toroidal
$\bar{B}_{\phi}(r=60,\theta,t)$ (bottom panel) for the Initial RIAF run. Time
is expressed in units of local orbit at $r=60$. Both radial and toroidal
fields show irregularity, which is typical feature of dynamo in a
geometrically thick RIAF.

We find that our Initial RIAF simulation is in the SANE state and an MRI
dynamo generates the large-scale magnetic fields and governs the magnetic
field evolution at late times as discussed in section 3.1. To characterize the
dynamo action, it is customary to visualise the spatio-temporal variation of
the mean magnetic field to investigate dynamo. We define the mean magnetic
field as the azimuthally averaged field

$\bar{B}_{i}(r,\theta)=\frac{1}{\phi_{\rm ext}}\int_{0}^{\phi_{\rm
ext}}B_{i}(r,\theta,\phi)d\phi,$ (25)

where $\phi_{\rm ext}$ is the extension in the $\phi$ direction, and
$i\in(r,\theta,\phi)$. Fig. 6 shows the variation of mean radial
$B_{r}(R_{0},,\theta,t)$ (top panel) and mean toroidal field
$\bar{B}_{\phi}(R_{0},\theta,t)$ (bottom panel) with latitude
($90^{\circ}-\theta$) and time at a radius $R_{0}=60$. This is also known as
the butterfly diagram. Both radial and toroidal fields show irregular
behaviour in their butterfly diagrams. Additionally, radial field is less
coherent compared to the toroidal field as observed in earlier studies of
dynamo in the SANE RIAF (Hogg & Reynolds, 2018b; Dhang et al., 2020). This
intermittent dynamo cycle in the RIAF is in contrast to the very regular
dynamo cycles observed in a thin Keplerian disk (e.g. see Flock et al.
(2012)). Irregularity in dynamo cycle arises because of slightly sub-Keplerian
angular velocity of the geometrical thick RIAF (Dhang & Sharma 2019).

Earlier studies found that while a large-scale dynamo generates large-scale
magnetic fields in the high latitudes, a fluctuation dynamo dominates close to
the disk mid-plane suppressing the production of large-scale magnetic field
there (Dhang & Sharma 2019). This can be qualitatively understood by the
looking at the large and coherent magnetic structures (specially for the
toroidal fields) in the high latitudes, while more patchy distribution near
the disk mid-plane ($90^{\circ}-\theta=0^{\circ}$) in the butterfly diagram in
Fig. 6 and also in the last two panels of Fig. 2. However it should be
emaphasized that although the MRI dynamo does produce large-scale magnetic
field, it is not strong enough to create a MAD, which is conducive for strong
jets (section 3.3). This inefficiency is likely to be due to the weak
$\alpha$-effect (Dhang et al. 2020) which is responsible for poloidal field
generation. Additionally, a strong turbulent pumping present in MRI-active
RIAF tends to prevent accumulation of large-scale magnetic field near the BH
as suggested in Dhang et al. (2020).

##  4 Transport of external magnetic field loops

Figure 7: Different types of external field loops injected at $t=8.06\times
10^{4}$. White dashed lines mark $z=1.5H$ above and below the mid-plane. For
details, see Table 1.

In this section, we study the accretion of external magnetic flux injected on
top of the fully turbulent SANE state obtained in Section 3. Our aim is to
investigate whether or not the system can bring in external magnetic flux
available at the outer radii, all the way to the central BH that may
eventually lead to a MAD state. While the actual configuration of the external
field is unknown and could be complex, we anticipate it likely to be closed.
Therefore, instead of the commonly used net vertical field, we inject the
poloidal magnetic field loops of different strengths, radial and vertical
sizes as shown in Fig. 7 and study their transport. As the controlled
experiments, these field loops confined between the radii $r_{l_{1}}$ and
$r_{l_{2}}$ are prescribed by

$A_{\phi,l}=\sqrt{\frac{2p_{\rm
gas}(r_{lc},\pi/2)}{C_{l}}}\left[\frac{\rho(r,\theta^{\prime})}{\rho(r,\pi/2)}-\delta_{l}\right]^{2}\sin\left[\kappa(R-r_{l1})\right]$
(26)

where $p_{\rm gas}$, $\rho$ are the initial pressure and density profiles
respectively, $\theta^{\prime}=\theta+\theta_{\rm shift}$,
$\kappa=\pi/(r_{l_{2}}-r_{l_{1}})$ and $r_{lc}=(r_{l_{1}}+r_{l_{2}})/2$.
Vanishing $\theta_{\rm shift}$ implies that the loop centre is at the mid-
plane, while a positive value of $\theta_{\rm shift}$ indicates that the loop
is off-centred. The vertical size of the loop is set by $\delta_{l}$. The
magnetization of the loop is controlled by the parameter $C_{l}$ and
characterized by $\beta_{l}=\langle p_{\rm gas}\rangle/\langle p_{\rm
mag,l}\rangle$, where $\langle p_{\rm gas}\rangle$ and $\langle p_{\rm
mag,l}\rangle$ are the gas pressure of the Initial RIAF at the time of loop
injection and magnetic pressure of the injected loop respectively.
Additionally note that average is performed within the loop. We choose
$r_{l_{1}}$ to be the inflow equilibrium radius $r_{eq}=60$, and different
values of $r_{l_{2}}$as tabulated in Table 1.

We restart the Initial RIAF run at $t=8.06\times 10^{4}$, inject the external
field loops (equation 26) and run till $t=t_{\rm end}$ as tabulated in Table
1. It is to be noted that we also run the Initial RIAF simulation longer to
compare it with the simulations with injected magnetic field loops. We
injected loops of different strengths with a wide range of plasma $\beta$
ranging from $\beta_{l}=7000$ (weak but stronger than the pre-existing mean
fields produced by MRI dynamo) to $\beta_{l}=70$ (very strong field typically
used in MAD simulations, but of much larger size than that used in our
simulations). We also explored the effects of other parameters such as radial,
vertical sizes and injection latitude of the loops on the flux transport
process while considering the loops of fiducial plasma $\beta$ values
$\beta_{l}=3500$ and $\beta_{l}=1500$ respectively. Additionally, we studied
the transport of big magnetic loops of strength (of $\beta_{l}=12200$) similar
to that of the mean fields produced by MRI dynamo in the quasi-stationary
phase of the Initial RIAF run. The configurations of injected field loops from
all these restart are illustrated in Figure 7.

Parameters
---
Name |  $C_{l}$ | $\delta_{l}$ |  $r_{l1}$ |  $r_{l2}$ |  $z_{l}$ |  $\beta_{l}$ |  $\theta_{\rm shift}$ |  $\Phi_{l,{\rm max}}$ |  $t_{\rm end}/10^{5}$
Initial RIAF |  - | - |  - |  - |  - |  - |  - |  - |  1.2
$\beta\\_7000$ |  $10^{-3}$ | 0.2 |  60 |  90 |  1.5 H |  7000 |  $0^{\circ}$ |  1.48 |  1
$\beta\\_3500$ |  $5\times 10^{-4}$ | 0.2 |  60 |  90 |  1.5 H |  3500 |  $0^{\circ}$ |  2.01 |  1.2
$\beta\\_1500$ |  $2.21\times 10^{-4}$ | 0.2 |  60 |  90 |  1.5 H |  1500 |  $0^{\circ}$ |  2.93 |  1.2
$\beta\\_700$ |  $10^{-4}$ | 0.2 |  60 |  90 |  1.5 H |  700 |  $0^{\circ}$ |  4.28 |  1.2
$\beta\\_200$ |  $2.78\times 10^{-5}$ | 0.2 |  60 |  90 |  1.5H |  200 |  $0^{\circ}$ |  7.95 |  1
$\beta\\_70$ |  $10^{-5}$ | 0.2 |  60 |  90 |  1.5 H |  70 |  $0^{\circ}$ |  13.14 |  1.2
$\beta\\_3500\\_{\rm tall}$ |  $1.168\times 10^{-3}$ | 0.0016 |  60 |  90 |  2.5H |  3500 |  $0^{\circ}$ |  2.04 |  1
$\beta\\_1500\\_{\rm tall}$ |  $5\times 10^{-4}$ | 0.0016 |  60 |  90 |  2.5 H |  1500 |  $0^{\circ}$ |  2.97 |  1
$\beta\\_3500\\_{\rm big}$ |  $1.429\times 10^{-5}$ | 0.2 |  60 |  120 |  1.5H |  3500 |  $0^{\circ}$ |  3.04 |  1
$\beta\\_1500\\_{\rm big}$ |  $6.12\times 10^{-5}$ | 0.2 |  60 |  120 |  1.5 H |  1500 |  $0^{\circ}$ |  4.42 |  1
$\beta\\_12200\\_{\rm big}$ |  $5\times 10^{-4}$ | 0.2 |  60 |  120 |  1.5 H |  12200 |  $0^{\circ}$ |  1.70 |  1.1
$\beta\\_3500\\_{\rm offc}$ |  $5.88\times 10^{-4}$ | 0.2 |  60 |  90 |  1.5H |  3500 |  $15^{\circ}$ |  1.56 |  1
$\beta\\_1500\\_{\rm offc}$ |  $2.55\times 10^{-4}$ | 0.2 |  60 |  90 |  1.5H |  1500 |  $15^{\circ}$ |  2.32 |  1
Table 1: Details of the injected external field loops characterized by plasma
$\beta_{l}$, vertical size $z_{l}$ and confinement radii between $r_{l1}$ and
$r_{l2}\geq r_{l1}$. $\theta_{\rm shift}$ and $\Phi_{l,{\rm max}}$ are tilt of
the loop with respect to the mid-plane and the total magnetic flux in the
injected loop in code units, respectively. Magnetic field loops are injected
at $t=8.06\times 10^{4}$ and run till $t=t_{\rm end}$.

###  4.1 Diagnostics

Before discussing the results in details, we define the following quantities
used to discuss the transport of external magnetic flux, including:

The radial magnetic flux threading the $r=const$ surface in the northern
hemisphere

$\Phi_{NH}(r)=\int_{\theta=0}^{\pi/2}\int_{0}^{\phi_{ext}}\sqrt{4\pi}B_{r}(r,\theta,\phi)\sqrt{-g}\
d\theta d\phi\ ,$ (27)

the vertical flux threading the mid-plane region

$\Phi_{\rm
mid}(r)=\int_{r=r_{H}}^{r}\int_{0}^{\phi_{ext}}-\frac{\sqrt{4\pi}}{r}B_{\theta}(r,\theta=\pi/2,\phi)\sqrt{-g}\
dr\ d\phi\ ,$ (28)

and the total flux available for accretion at different $r$ in the northern
hemisphere

$\Phi_{\rm tot}(r)=\Phi_{\rm NH}(r_{H})+\Phi_{\rm mid}(r);$ (29)

where $r_{H}=2r_{g}$ is the event horizon radius of the BH. We also define a
normalised flux representing the efficiency of the flux transport and defined
by

$f_{B}=\frac{\Phi_{NH}(r_{H})}{\Phi_{l,{\rm max}}}.$ (30)

Here, $\Phi_{l,{\rm max}}$ is the total flux at the loop centre at the time of
injection or in the beginning of the simulation (only for the Initial RIAF
run; also see Table 1).

###  4.2 Results for the fiducial parameter- plasma $\beta$

First we will discuss the dependence of the flux transport and the emergent
accretion properties on the strength of the loops, defined by plasma
$\beta_{l}$ (also see Table 1 and Fig. 7). We start by discussing the
qualitative picture on the evolution of the magnetic flux injected between the
radii $r=60$ and $r=90$ on top of the existing magnetic field in the quasi-
stationary Initial RIAF. It is worth noting that the plasma $\beta$ of the
total (mean + fluctuation) magnetic field is $\beta_{\rm tot}=70$, while that
of the mean field alone is $\beta_{\rm mean}=12250$, for the Initial RIAF run.
Fig. 8 shows the evolution of external magnetic field loops of three different
$\beta_{l}$-s. Colour shows the intensity of mean radial field $\bar{B}_{r}$,
while streamlines describe the mean poloidal fields
$\mathbf{\bar{B}_{p}}=\mathbf{\bar{B}_{r}}+\mathbf{\bar{B}_{\theta}}$.

Figure 8: Evolution of external magnetic field loops injected in a turbulent
quasi-stationary RIAF. Panels at the top, middle, and bottom show the
evolution of the magnetic field loops of plasma $\beta=7000$ (weak field
strength), $\beta=1500$ (moderate field strength ) and $\beta=70$ (strong
field), respectively. Colour represents the mean radial field $\bar{B}_{r}$;
while poloidal field distribution
$\mathbf{\bar{B}_{p}}=\mathbf{\bar{B}_{r}}+\mathbf{\bar{B}_{\theta}}$ is
described by streamlines. Note that a fraction of the injected magnetic flux
reaches the BH.

Top panels show the time evolution of the weakly magnetized loop of strength
$\beta_{l}=7000$. Injection of the weak external magnetic field loops re-
excite the MRI in the accretion flow, enhances accretion stresses (see Fig.
12) and hence lead to higher mass accretion rates (see Fig. 11). Poloidal flux
slowly drifts towards the BH and a fraction of the injected flux accumulates
near the BH (which is quantitatively shown in Fig. 9). We see an increase in
the radial magnetic flux threading the BH when compared to that in the Initial
RIAF run. This can be comprehended by comparing the snapshots at $t=80601$
(the flux level close to the BH does not change significantly in the quasi-
steady state of Initial RIAF as shown in Fig. 2) and last panel at $t=10^{5}$.

Next we discuss the transport of moderately strong magnetic field loops of
$\beta_{l}=1500$ as shown in the middle panels of Fig. 8. Magnetic flux
reaches the BH in a shorter time compared to that in the weak field case of
$\beta_{l}=7000$. This is due to the stronger accretion stresses produced (see
Fig. 12) in the accretion flow due to the injection of stronger magnetic field
loops. Radial magnetic field strength in the polar region is found to be
stronger than the weak field case at late times. This is because of the larger
amount of flux associated with the loop of $\beta_{l}=1500$ than that with the
loop of $\beta_{l}=7000$.

Finally, we examine the transport of the strong field loops with
$\beta_{l}=70$, similar to the strength of the total (mean + fluctuation)
magnetic field in the qausi-stationary phase of Initial RIAF. Unlike the
previous two weak-field cases, here the injected loops are so strong with most
unstable wavelengths comparable to the disk scale height, and that drives
strong channel flows over the entire vertical extent of the disk. This further
generates a spike in the large-scale Maxwell stress (see Fig. 12) producing a
strong inflow of mass and magnetic flux. Magnetic flux reaches the BH very
quickly and fills the polar region. The system remains to be in a strongly
turbulent state till the end of the simulation.

The qualitative pictures discussed above are representative among all our
simulations. The animations of the other runs with external magnetic field
loops along with the Initial RIAF run can be viewed in this playlist. In the
upcoming sub-sections, we will quantify with different metrics of magnetic
flux transport in greater detail for all the runs we performed listed in Table
1.

####  4.2.1 Evolution of magnetic flux

Figure 9: Left panel: Comparison of amount of radial magnetic flux threading
the event horizon in the northern hemisphere $\Phi_{\rm NH}(r_{H})$ for the
runs with injected loops of different plasma $\beta$. Right panel: Fraction of
the injected magnetic flux reaching the BH for runs with different plasma
$\beta$.

Fig. 9(a) shows the time evolution of the magnetic flux through the event
horizon in the northern hemisphere, $\Phi_{NH}(r_{H})$ for runs with different
strengths of injected magnetic field loops, which are also compared with that
of the Initial RIAF run. It is clearly visible that the injection of an
external loop enhances the amount of flux at the event horizon. In the runs
with low $\beta_{l}\lesssim 1000$, there is a transient rise of $\Phi_{\rm
NH}$ due to the fast transport from the strong MRI channel flow. In the more
extreme case of $\beta_{l}=70$, the transient phase is so extreme that it
leads to a strong initial spike in $\Phi_{\rm NH}$, followed by a gradual
decline towards a more steady flux level. For other runs with
$\beta_{l}\gtrsim 1000$, the build up of magnetic flux in the BH horizon is
more gradual, and the build up is slower for runs with with higher
$\beta_{l}$. It is worth noting that none of our simulations with injected
loops reach the MAD state, with MAD parameter ranging from $\phi_{BH}=2$ to
$\phi_{BH}=10$.

Figure 10: Space-time plot for $\Phi_{\rm tot}$ defining the amount of
poloidal magnetic flux at the disk midplane. Horizontal green dashed line
shows the location of ISCO. Vertical black dashed line in the top left panel
(the Initial RIAF) denotes the time of injection of external flux.

We show the spatial-temporal variation of total flux $\Phi_{\rm tot}(r,t)$
available for accretion in the northern hemisphere in Fig. 10 to obtain a more
complete picture of the flux transport at different radii. Each panel of Fig.
10 describes the evolution of the radial profile of $\Phi_{\rm tot}$ over time
for runs with different $\beta_{l}$. The first panel in the top row
corresponds to the Initial RIAF run. It again demonstrates that the system
forgets its initial magnetic field configuration after the time $3-4\times
10^{4}$. The rest of the panels show the spatial-temporal evolution of
$\Phi_{\rm tot}(r,t)$ for other runs after we inject external magnetic field
loops. In accordance with Figure 9, we see that there are two regimes of flux
transport depending on the strength of the injected loop. With very strong
external flux ($\beta_{l}=70,\ 200$), the external flux is quickly transported
both inwards and outwards, characteristic of the channel flows with flow
directions alternating over height, as seen in the last tow bottom panels of
Fig. 10. The channel flows lead to an initial transient transport of a large
fraction of initial flux into the BH, followed by subsequent relaxation and
diffusion towards a more steady flux level.

With weak external flux ($\beta_{l}\gtrsim 3500$), the initial external
magnetic flux gradually diffuses while being advected inwards. In the end, a
fraction the flux overcomes diffusion to reach the BH, which will be discussed
more quantitatively in the next subsection. In between between these two
regimes, there lies the case of moderately strong external flux
($\beta_{l}=700,\ 1500$), for which the flux transport by the channel flows
diffuses before reaching the BH, and subsequent transport is likely mediated
by a combination of advection and diffusion.

####  4.2.2 The Efficiency of transport

Till now, we considered the total flux as the primary diagnostics irrespective
of the amount of flux associated with the injected loops. However, it is worth
noting that different loops have different amount of fluxes. Therefore, a
normalised flux $f_{B}$ as defined in equation 30 would be the better
indicator on the efficiency of flux transport.

Fig. 9(b) shows the time variation of the fraction $f_{B}$ for different runs.
In the regime of very strong injected flux ($\beta_{l}=70,\ 200$), the
efficiency is quite high (up to $\sim 50\%$) during the initial phase when
channel flows dominate. Later, the efficiency goes down to around 15-20$\%$.
In the weak field regime ($\beta_{l}\geq 3500$), despite that flux in the BH
is accumulated gradually, the efficiency of flux transport is more or less
similar, which is around 15-20 percent. For comparison, we also show the
result for the Initial RIAF run, where we define $\Phi_{l,{\rm max}}$ by
calculating flux at a radius $r=75$ at $t=0$. Finally, it is interesting to
note that flux transport appears to be more efficient, reaching about
20-40$\%$, when the field strength is in between the two regimes, i.e for
$\beta_{l}=700$ and $\beta_{l}=1500$.

####  4.2.3 Effects on mass accretion rate and accretion stresses

Figure 11: Time evolution in mass accretion rate $\dot{m}_{\rm in}$ at the
event horizon after the injection of magnetic field loops of different
strengths defined by plasma $\beta_{l}$. Figure 12: Space-time plot of the
radial profile of the average midplane total accretions tress $\langle W_{\rm
Tot}\rangle=\langle W_{\rm Max}\rangle+\langle W_{\rm Rey}\rangle$ for
different runs, where average is done in the azimuthal and vertical directions
within one-scale height from the midplane. Vertical black dashed line in the
top left panel indicates the injection time of field loops. Horizontal black
dashed lines in the rest of the panels mark the radial range where field loops
are injected (between $r=60$ and $r=90$).

In this subsection, we study how the injection of external magnetic flux
influences the accretion properties such as the mass accretion rate and
accretion stresses. Fig. 11 shows time history of mass accretion rate at the
event horizon for runs with a range of $\beta_{l}$. We further show in Fig. 12
the space-time plot of the azimuthally and vertically (over one scale-height)
averaged total accretion stress $\langle W_{\rm Tot}\rangle$, which is a
combination of Maxwell and Reynolds stresses defined in the orthonormal fluid
frame (see section 3.2) as,

$\displaystyle W_{\rm Max}=2p_{\rm mag}\ u^{\hat{r}}\
u^{\hat{\phi}}-b^{\hat{r}}\ b^{\hat{\phi}},$ (31) $\displaystyle W_{\rm
Rey}=\left(\rho+\frac{\gamma}{\gamma-1}p_{\rm gas}\right)u^{\hat{r}}\
u^{\hat{\phi}},$ (32) $\displaystyle\langle W_{\rm Tot}\rangle=\langle W_{\rm
Max}\rangle+\langle W_{\rm Rey}\rangle\ ,$ (33)

where the Maxwell stress is the dominant component.

Fresh injection of external field loops reignites the linear MRI and lead to
higher accretion stresses and hence increase in mass accretion rate. We find
that with high field strength in the loop (i.e $\beta_{l}\lesssim 200$), there
is substantial enhanced accretion stress, leading to a rapid, strong and
transient increase of accretion rate. The stresses are reduced after the
transient phase, but are still much stronger than those in the Initial RIAF
run within simulation time, reflecting the prolonged influence of the initial
flux loop. Simulations with $\beta_{l}\gtrsim 3500$ show only modest increase
of accretion stress and the accretion rates compared to the Initial RIAF run,
indicating that the injection of external flux has only minor impact on disk
turbulence. For simulations with intermediate $\beta_{l}$, there is modest
enhancement of the accretion stress, resulting in modest enhancement of
accretion rate. We also note that after $t=1.1\times 10^{5}$, despite having a
higher level of magnetic flux (compared to Initial RIAF run; Fig. 10), mass
accretion rate in the runs with injected loops are very similar to that in the
Initial RIAF run. This happens due to the quick depletion of mass supply in
the disk at earlier times due to the enhanced stresses in the runs with
injected loops.

####  4.2.4 Disk and flow structures

Figure 13: Radial profiles of surface density ($\langle\Sigma\rangle$, left)
and radial velocity ($\langle u^{1}\rangle$) for the Initial RIAF run and runs
with injected loops with fiducial size and location with different strengths.
Time average is done over $t=9\times 10^{4}-10^{5}$. For $\langle
u^{1}\rangle$, spatial average is done over the azimuth and within one scale-
height about the midplane.

In this subsection, we further examine how external magnetic flux changes the
disk structure and flow properties. We start by considering the radial
profiles of surface density, defined as111We note that the standard definition
(34) asymptotes to $\rho r^{2}d\theta$ at large radii, with an extra $r$
factor compared to the Newtonian definition (assuming surface density is
defined by integrating along spherical shells).

$\Sigma(t,r)=\frac{1}{\phi_{\rm ext}}\int_{\phi=0}^{\phi_{\rm
ext}}\int_{\theta=0}^{\pi}\rho\ \sqrt{-g}\ d\theta\ d\phi\ ,$ (34)

and the radial velocity $\langle u^{1}(r)\rangle$, averaged within one scale
height about the midplane. The results for the Initial RIAF run and runs with
external field loops, time average is done over $t=9\times 10^{4}-10^{5}$, are
shown in Fig. 13.

In the Initial RIAF run, we see that the accretion velocity approaches the
free-fall velocity ($v_{ff}=\sqrt{2/r}$) within the ISCO, while accretion
velocity ranges between $0.01-0.5$ of the Keplerian velocity further out till
the radius of inflow equilibrium. Upon imposing external field, the higher
accretion stresses leads to higher accretion velocities. The enhancement can
be up to a factor of $\sim 10$ in the strong field case with $\beta_{l}=70$ at
the representative radii of $r\sim 30-60r_{g}$, while for weak field runs
(e.g., $\beta_{l}\gtrsim 3500$), the accretion velocity is only enhanced by a
modest factor of $\sim 2$. We also note that the profile of $\langle
u^{1}(r)\rangle$ also evolves over time accompanying magnetic flux transport,
but qualitatively, the profiles shown in Figure 13 are representative over the
duration of our simulations.

The altered accretion velocity profile $\langle u^{1}(r)\rangle$ further
modifies the surface density profile. Generally, after imposing an external
field loop, the surface density becomes steeper compared to the surface
density profile in the Initial RIAF run, though the deviation is only modest.
The surface density profile also evolves over time. We note that earlier RIAF
simulations of the SANE state already indicated that there might not be any
universal power law for the surface density profile and other flow properties
(White et al., 2020). When supplied with external magnetic flux in the outer
disk, our results suggest additional surface density variations during the
process of magnetic flux transport. In other words, the dynamics of RIAFs are
dependent on the magnetised mass reservoir at larger radii.

###  4.3 Results for the other parameters

In this section, we assess the robustness of our fiducial simulations results
by considering different geometries for the injected field loops. In
particular, we change loop sizes (both vertical and radial) and the injection
latitudes. We focus on loops of fiducial strength $\beta_{l}=3500$ and
$\beta_{l}=1500$ respectively. Additionally, we study the transport of a big
loop of similar strength ($\beta_{l}=12200$) to mean poloidal fields produced
by MRI dynamo in the Initial RIAF run in the quasi-stationary phase. Fig. 14
shows the time evolution of the radial magnetic flux threading the event
horizon in the northern hemisphere $\Phi_{\rm NH}(r_{H})$ (top panels) and
flux transport efficiency (bottom panels) for these additional simulations.

Figure 14: The time evolution of the amount of radial magnetic flux threading
the event horizon in the northern hemisphere $\Phi_{\rm NH}(r_{H})$ (top
panels) and its fraction over the injected magnetic flux (bottom panels) for
the runs with injected loops of larger vertical size (left panel), larger
radial size (middle panel), and off-centered injection latitudes (right). The
results are compared with runs with fiducial loop geometries and the Initial
RIAF run with the same level of magnetization ($\beta_{l}=1500$ and $3500$).

####  4.3.1 Vertical and Radial sizes

The left panels of Fig. 14 compare the flux transport for taller loops of
vertical size $z_{l}=2.5H$ with the loops of similar strength but of fiducial
size $z_{l}=1.5H$. We find that the vertical size of the loops does not affect
the amount of flux reaching the BH, and the efficiency of flux transport
remains almost unaltered with the change of the vertical size of the loops.

Similarly, the middle panels of Fig. 14 compares flux transport between loops
of different radial sizes, where we consider bigger loops of radial size
$\Delta r_{l}=60$ as opposed to the fiducial radial size of $\Delta r_{l}=30$.
We observed that a larger amount of flux reaches the BH for the bigger loops,
which is reasonable because more magnetic flux is available in these loops
compared to its smaller counterparts. However, the fraction of the flux
reaching the BH remains similar for both the smaller and bigger loop cases
with the same plasma $\beta_{l}$. This result also holds for our additional
run with $\beta_{l}=12200$. This indicates that the efficiency of flux
transport remains unaffected by the radial extent of the injected loops.

####  4.3.2  Injection latitude

In addition to study the effects of strength and size of the loops on the
transport process, we also consider injecting off-centred loops with plasma
$\beta_{l}=3500$ and $\beta_{l}=1500$ to examine whether loop injection away
from the mid-plane facilitates flux transport or not. The comparison with our
fiducial injection prescription are shown in the right panel of Fig. 14.
Surprisingly, injection of off-centred loops lead to a distinctly lower flux
level at the event horizon. While it is not entirely clear why this is the
case, we speculate that it is related to stronger magnetic reconnection in the
off-centred case that leads to more considerable destruction of magnetic flux,
which occurs during the interplay between the injected field and the dynamo-
generated background field.

The stark difference between the magnetic field evolution in off-centred and
the fiducial cases can be also seen by comparing the movies describing the
magnetic field loop evolution for the runs beta_1500 (movie-beta-1500) and
beta_1500_offc (movie-beta-1500-offc) respectively.

##  5 Discussion

###  5.1 Inefficiency of Dynamo in SANE/RIAF

We threaded the initial geometrically semi-thick disk ($H/R\approx 0.2$) with
small magnetic field loops of alternating polarity and attained a quasi-
stationary weakly magnetized RIAF (SANE; see Fig. 4), that does not remember
the initial field geometry (see sections 3.1 and 3.3). An MRI dynamo is
responsible for generating and sustaining magnetic fields (both small scale
and large-scale) in the quasi-stationary RIAF (Hawley et al. 2013; Hogg &
Reynolds 2018a). A large-scale dynamo does operate (Dhang & Sharma 2019) and
generate large-scale magnetic fields in the weakly magnetized RIAF (see last
two panels of Fig. 2), but not efficient enough to produce strong magnetic
fields that can convert a SANE to MAD. This result aligns with earlier works
which found dynamo action in a SANE RIAF does not lead to the jet formation
(Beckwith et al. 2008; Narayan et al. 2012). This inefficiency is likely to be
attributed to insufficient poloidal field generation (weak $\alpha$-effect)
and strong turbulent pumping which transports large scale magnetic field
radially outward in a RIAF(Dhang et al. 2020).

Recently, Liska et al. (2020) reported that when starting the simulation with
an unusually strong ($\beta\approx 5$) and coherent toroidal magnetic field,
the MAD state can be achieved at late times. They argued that an MRI dynamo
can produce strong poloidal field loops of size $H\propto R$ from the very
strong and coherent initial toroidal field. The further the creation location
is, the bigger the loops are. Most of the loops move outward, while a few
‘lucky’ loops created at large radii are somehow arrested and stretched inward
and lead to the MAD state. However, how the accretion disk can possess such a
coherent initial toroidal field of the same polarity spanning several decades
in radii at first place remains questionable.

Overall, we reaffirm that the MRI dynamo in standard SANE state does not
spontaneously generate strong coherent large-scale poloidal field to turn the
disk into the MAD state. In the absence of initial poloidal field, achieving
the MAD state may require unusually strong and coherent toroidal field that
may be unpractical in reality.

###  5.2 Plausible sources of external magnetic fields

In this work, we considered the possibility that the disk acquires external
poloidal field in the form of field loops of different sizes and shapes. What
can be the source of such external field loops? While definitive evidence is
lacking, we speculate that accreting such external field loops is plausible in
a variety of systems.

In the hard state of the XRBs, a RIAF close to the BH is proposed to be
connected to an outer thin disk Esin et al. 1997; Done et al. 2007, which can
supply large-scale magnetic flux to the inner RIAF. The outer thin disk can in
principle harbour large-scale magnetic field due to an efficient dynamo action
(Flock et al. 2012; Gressel & Pessah 2015) or due to the coronal accretion of
magnetic flux (Guilet & Ogilvie 2012) from the companion/donor star, or a
combination of both. The donor star in the low mass XRBs are likely to be
either K or M type dwarf stars (Fragos & McClintock 2015) or evolved stars
(e.g., as in GRS 1915+1105). The donor stars in the XRBs are supposed to be
tidally locked to the rotation period of the binaries with an orbital period
of hours to days (Coriat et al. 2012). These fast rotating dwarf stars show
vigorous magnetism with surface magnetic field of strength ($\sim 10^{3}$ G)
similar to sunspots (West et al. 2008; Davenport 2016). Additionally, in the
active region, the magnetic field is one order of magnitude stronger than the
average stellar magnetic field. Magnetized matter from the donor star passes
through the first Lagrange point ($L_{1}$) and enters the Roche lobe of
primary (accretor) almost ballistically and circularizes at the
circularization radius (Frank et al. 2002). We speculate that the mass-loss
from the L1-nozzle may proceed through a chain of mass blobs encircled by
field loops (e.g., as also considered in Ju et al. 2017), which may gets
amplified and become quasi-axisymmetric during the circularization process.
Thus, if this external flux can be brought in through the outer thin Keplerian
disk, then it may further feed the inner RIAF, where flux transport is
efficient and saturate the BH.

The accretion flow in a low-luminosity AGN is also thought be a RIAF. In this
case, the gas supplied by the ambient medium to the accretion flow is
magnetized. It can harbour a large-scale magnetic field as inferred from the
observation of the large scale poloidal flux in the Galactic centre (Nishiyama
et al. 2010). Recent numerical simulations by Ressler et al. (2020a, b) found
that the large-scale accretion flow around the galactic centre fed by the
winds of Wolf-Rayet stars can achieve the MAD state, with efficient inward
transport of magnetic field embedded in the accreting material. Their injected
magnetic fields have pure toroidal component with random orientation, thus we
may consider that such fields effectively enter the accretion disk in the form
of closed field loops from some random directions. Our results are in line
with their findings, while our controlled experiments further provide a
physical basis for better understanding the efficient flux transport around
SMBHs.

###  5.3 Transport efficiency in the SANE and its possibility of
transformation to MAD

In Section 4, we have seen that except for off-centered loops, transport of
externally injected magnetic flux loops is relatively efficient, with
typically $\sim 20\%$ of the available flux end up being accreted to the
central BH regardless of initial field strength and size.

Note that none of our simulations reach the MAD state, but if this relatively
high efficiency of flux transport obtained from our controlled experiments is
universal, we can estimate the requirement on the external flux to potentially
transform a SANE disk into MAD.

The MAD parameter $\phi_{\rm BH}$ (the normalised unsigned flux threading the
BH) is related to the magnetic flux (signed) $\Phi_{NH}$ threading the
northern hemisphere of the BH as follows

$\phi_{\rm BH}\approx\frac{\Phi_{NH}(r_{H})}{\sqrt{\dot{m}}r_{g}c^{1/2}}.$
(35)

We have found that a certain fraction $f_{B}=\Phi_{NH}(r_{H})/\Phi_{in}$
(equation 30) of the injected flux $\Phi_{\rm in}$ reaches the BH. Earlier
numerical experiments suggest that the MAD could be achieved if the MAD
parameter exceeds a critical value $\phi_{BH,c}$ at the event horizon (e.g.,
see Tchekhovskoy et al. 2011). This indicates to a critical value of the
injected flux $\Phi_{in,c}$ which is a plausible minimum flux required for the
MAD state and it is given by

$\Phi_{\rm in,c}=\left(\frac{r_{g}c^{1/2}}{f_{B}}\right)\phi_{{\rm
BH},c}\sqrt{\dot{m}}.$ (36)

Therefore, the minimum poloidal magnetic field required at the injection
location is given by

$B_{\rm in,c}=\left(\frac{r_{g}c^{1/2}}{2\pi
f_{B}}\right)\left(\frac{r}{\Delta r}\right)\left(\frac{\phi_{{\rm
BH},c}}{r^{2}}\right)\sqrt{\dot{m}},$ (37)

where we have estimated that for a loop centered on radius $r$ with half width
$\Delta r$, $\Phi_{\rm in,c}\approx 2\pi B_{\rm in,c}r\Delta r$. If we take
$\Delta r/r=0.2$, then the value of $B_{\rm in,c}$, in terms of Eddington
accretion rate $\dot{M}_{\rm Edd}=1.5\times 10^{19}\ M_{10}\ {\rm gm}\ {\rm
s}^{-1}$ is given by

$B_{\rm in,c}\approx\frac{10^{4}}{f_{B}}\ \phi_{40}\ r_{100}^{-2}\
\dot{m}_{-4}^{1/2}\ M_{10}^{-1/2}\ G,$ (38)

where $\phi_{40}=\phi_{{\rm BH},c}/40$, $r_{100}=r/(100r_{g})$,
$\dot{m}_{-4}=\dot{m}/10^{-4}\dot{M}_{\rm Edd}$, and
$M_{10}=M_{BH}/10M_{\odot}$.

Would this amount of magnetic field be available for accretion at the outer
radii of the RIAF? We will estimate the poloidal field strength available for
accretion in case of an XRB. In the low hard state, the RIAF close to the BH
is proposed to be connected to an outer thin disk. The plausible source of the
large-scale magnetic field in the thin disk could be the dynamo action.
Another scenario would be the advection of large-scale field loops from the
companion star as discussed in section 5.2. Independent of the mechanism, we
can estimate the characteristic poloidal field strength in the thin accretion
disk given that the accretion is driven by radial transport of angular
momentum in the disk as (e.g., Bai & Goodman 2009)

$-\overline{B_{r}B_{\phi}}\approx\frac{\dot{m}\Omega}{h_{a}}.$ (39)

Here, $h_{a}=\xi H_{\rm thin}$ is the thickness of the disk over which
accretion proceeds. Further, if we assume that $|B_{r}|\approx 1/5|B_{\phi}|$
(which is found to be consistent in the MRI simulations of accretion disks)
and $\xi=6$, then the total radial magnetic field in the thin disk of aspect
ratio $\epsilon_{\rm thin}=H_{\rm thin}/R$ is given by,

$B_{r,d}\approx 10^{4}\ \epsilon^{-1/2}_{0.05}\ m_{-4}^{1/2}\ r_{100}^{-5/4}\
M_{10}^{-1/2}\ G,$ (40)

where $\epsilon_{0.05}=\epsilon_{\rm thin}/0.05$. It is to be noted that in
the strongly magnetized coronal region of the thin disk, a large share of this
estimated total field $B_{r,d}$ will likely be in the mean coherent part of
the magnetic field. In reality, the mass accretion rate in the outer thin disk
is expected to be higher compared to that in a RIAF (Yuan & Narayan 2014), and
hence the total poloidal field will also be higher. Therefore, comparison of
$B_{{\rm in,c}}$ and $B_{r,d}$ leads to the inference that it is quite
possible that in an XRB, the outer thin disk reservoir can potentially supply
adequate amount of magnetic flux to the inner RIAF that eventually may form a
MAD close to the BH.

##  6 Summary

In this paper, we studied the magnetic field generation and transport in a
geometrically thick RIAF. We initialize the disk with magnetic field loops of
alternate polarity so that the quasi-stationary RIAF is in the weakly
magnetized, i.e in the SANE regime. In this quasi-sationary turbulent SANE
RIAF, we study the transport of external magnetic flux (in the form of loops)
of different strengths, sizes and shapes. Here we outline the key findings of
our work.

  * •

We reconfirm that the MRI dynamo in a standard SANE RIAF does not generate
strong coherent large-scale poloidal field to turn the SANE state into the MAD
state.

  * •

Magnetic flux transport is relatively efficient in the SANE RIAF: fifteen to
forty percent of the external magnetic flux injected at the outer radii is
able to reach the BH.

  * •

Flux transport efficiency is independent of the loop parameters such as
strength and size. However if the loops are injected at high latitudes rather
than at the mid-plane, the efficiency becomes poor.

We also find that accretion flow profiles (e.g surface density, accretion
velocity) are altered as external magnetic flux is injected in the disk. We
propose that the dynamics of the RIAF depends on the magnetized mass reservoir
at the outer radii.

Based on our results, we argue that it might be easier to transform a SANE
disk to a MAD by supplying external poloidal field loops at the outer disk
provided that the relatively high efficiency of flux transport obtained from
our controlled experiments is universal.

It should also be noted that we have studied transport of external magnetic
flux in the quasi-stationary turbulent RIAF in a limited parameter space. For
example, we have considered only one injection location with inner edge of the
loop being at $r=60$, whereas, in reality, the loops are supposed to be
available for accretion as far as in the disk truncation region in XRBs, or
even at a larger radii in Low-luminosity AGNs. We plan to explore magnetic
flux transport with different configurations and with larger dynamical ranges.
Furthermore, future work should extend this study to the thin disk regime,
which is applicable to regions beyond the truncation radius in the low/hard
state of the XRBs, as well as in luminous AGNs.

We thank Ramesh Narayan for initial inputs into this project. This research
was supported by NSFC grant 11873033. Numerical simulations are conducted on
TianHe-1 (A) at National Supercomputer Center in Tianjin, China, and on the
Orion cluster at Department of Astronomy, Tsinghua University.

All the movies of the simulations mentioned in Table 1 are available in this
link.

## References

  * Abramowicz et al. (1988) Abramowicz, M. A., Czerny, B., Lasota, J. P., & Szuszkiewicz, E. 1988, ApJ, 332, 646, doi: 10.1086/166683
  * Bai & Goodman (2009) Bai, X.-N., & Goodman, J. 2009, ApJ, 701, 737, doi: 10.1088/0004-637X/701/1/737
  * Bai & Stone (2013) Bai, X.-N., & Stone, J. M. 2013, ApJ, 769, 76, doi: 10.1088/0004-637X/769/1/76
  * Balbus & Hawley (1991) Balbus, S. A., & Hawley, J. F. 1991, ApJ, 376, 214, doi: 10.1086/170270
  * Beckwith et al. (2008) Beckwith, K., Hawley, J. F., & Krolik, J. H. 2008, ApJ, 678, 1180, doi: 10.1086/533492
  * Beckwith et al. (2009) —. 2009, ApJ, 707, 428, doi: 10.1088/0004-637X/707/1/428
  * Bisnovatyi-Kogan & Ruzmaikin (1974) Bisnovatyi-Kogan, G. S., & Ruzmaikin, A. A. 1974, Ap&SS, 28, 45, doi: 10.1007/BF00642237
  * Blandford & Begelman (1999) Blandford, R. D., & Begelman, M. C. 1999, MNRAS, 303, L1, doi: 10.1046/j.1365-8711.1999.02358.x
  * Blandford & Payne (1982) Blandford, R. D., & Payne, D. G. 1982, MNRAS, 199, 883, doi: 10.1093/mnras/199.4.883
  * Blandford & Znajek (1977) Blandford, R. D., & Znajek, R. L. 1977, MNRAS, 179, 433, doi: 10.1093/mnras/179.3.433
  * Chakrabarti (1989) Chakrabarti, S. K. 1989, ApJ, 347, 365, doi: 10.1086/168125
  * Colella & Woodward (1984) Colella, P., & Woodward, P. R. 1984, Journal of Computational Physics, 54, 174, doi: 10.1016/0021-9991(84)90143-8
  * Coriat et al. (2012) Coriat, M., Fender, R. P., & Dubus, G. 2012, MNRAS, 424, 1991, doi: 10.1111/j.1365-2966.2012.21339.x
  * Davenport (2016) Davenport, J. R. A. 2016, ApJ, 829, 23, doi: 10.3847/0004-637X/829/1/23
  * Dhang et al. (2020) Dhang, P., Bendre, A., Sharma, P., & Subramanian, K. 2020, MNRAS, 494, 4854, doi: 10.1093/mnras/staa996
  * Dhang & Sharma (2019) Dhang, P., & Sharma, P. 2019, MNRAS, 482, 848, doi: 10.1093/mnras/sty2692
  * Done et al. (2007) Done, C., Gierliński, M., & Kubota, A. 2007, A&A Rev., 15, 1, doi: 10.1007/s00159-007-0006-1
  * Einfeldt (1988) Einfeldt, B. 1988, 671.  https://ui.adsabs.harvard.edu/abs/1988stw..proc..671E
  * Esin et al. (1997) Esin, A. A., McClintock, J. E., & Narayan, R. 1997, ApJ, 489, 865, doi: 10.1086/304829
  * Event Horizon Telescope Collaboration et al. (2021a) Event Horizon Telescope Collaboration, Akiyama, K., Algaba, J. C., et al. 2021a, ApJ, 910, L12, doi: 10.3847/2041-8213/abe71d
  * Event Horizon Telescope Collaboration et al. (2021b) —. 2021b, ApJ, 910, L13, doi: 10.3847/2041-8213/abe4de
  * Fender et al. (1999) Fender, R., Corbel, S., Tzioumis, T., et al. 1999, ApJ, 519, L165, doi: 10.1086/312128
  * Flock et al. (2012) Flock, M., Dzyurkevich, N., Klahr, H., Turner, N., & Henning, T. 2012, ApJ, 744, 144, doi: 10.1088/0004-637X/744/2/144
  * Fragos & McClintock (2015) Fragos, T., & McClintock, J. E. 2015, ApJ, 800, 17, doi: 10.1088/0004-637X/800/1/17
  * Frank et al. (2002) Frank, J., King, A., & Raine, D. J. 2002, Accretion Power in Astrophysics: Third Edition
  * Gardiner & Stone (2005) Gardiner, T. A., & Stone, J. M. 2005, Journal of Computational Physics, 205, 509, doi: 10.1016/j.jcp.2004.11.016
  * Ghisellini et al. (2014) Ghisellini, G., Tavecchio, F., Maraschi, L., Celotti, A., & Sbarrato, T. 2014, Nature, 515, 376, doi: 10.1038/nature13856
  * Goedbloed & Keppens (2022) Goedbloed, H., & Keppens, R. 2022, ApJS, 259, 65, doi: 10.3847/1538-4365/ac573c
  * Goodman & Xu (1994) Goodman, J., & Xu, G. 1994, ApJ, 432, 213, doi: 10.1086/174562
  * Gressel & Pessah (2015) Gressel, O., & Pessah, M. E. 2015, ApJ, 810, 59, doi: 10.1088/0004-637X/810/1/59
  * Guilet & Ogilvie (2012) Guilet, J., & Ogilvie, G. I. 2012, MNRAS, 424, 2097, doi: 10.1111/j.1365-2966.2012.21361.x
  * Harrison (2017) Harrison, C. M. 2017, Nature Astronomy, 1, 0165, doi: 10.1038/s41550-017-0165
  * Hawley et al. (2013) Hawley, J. F., Richers, S. A., Guan, X., & Krolik, J. H. 2013, ApJ, 772, 102, doi: 10.1088/0004-637X/772/2/102
  * Hogg & Reynolds (2018a) Hogg, J. D., & Reynolds, C. S. 2018a, ApJ, 854, 6, doi: 10.3847/1538-4357/aaa6c6
  * Hogg & Reynolds (2018b) —. 2018b, ApJ, 861, 24, doi: 10.3847/1538-4357/aac439
  * Igumenshchev et al. (2003) Igumenshchev, I. V., Narayan, R., & Abramowicz, M. A. 2003, ApJ, 592, 1042, doi: 10.1086/375769
  * Ju et al. (2017) Ju, W., Stone, J. M., & Zhu, Z. 2017, ApJ, 841, 29, doi: 10.3847/1538-4357/aa705d
  * Liska et al. (2020) Liska, M., Tchekhovskoy, A., & Quataert, E. 2020, MNRAS, 494, 3656, doi: 10.1093/mnras/staa955
  * Lubow et al. (1994) Lubow, S. H., Papaloizou, J. C. B., & Pringle, J. E. 1994, MNRAS, 267, 235, doi: 10.1093/mnras/267.2.235
  * McKinney et al. (2012) McKinney, J. C., Tchekhovskoy, A., & Blandford, R. D. 2012, MNRAS, 423, 3083, doi: 10.1111/j.1365-2966.2012.21074.x
  * Meier (2005) Meier, D. L. 2005, Ap&SS, 300, 55, doi: 10.1007/s10509-005-1184-9
  * Mishra et al. (2019) Mishra, B., Begelman, M. C., Armitage, P. J., & Simon, J. B. 2019, arXiv e-prints.  https://arxiv.org/abs/1907.08995
  * Morganti (2017) Morganti, R. 2017, Frontiers in Astronomy and Space Sciences, 4, 42, doi: 10.3389/fspas.2017.00042
  * Narayan et al. (2003) Narayan, R., Igumenshchev, I. V., & Abramowicz, M. A. 2003, PASJ, 55, L69, doi: 10.1093/pasj/55.6.L69
  * Narayan et al. (2012) Narayan, R., SÄ dowski, A., Penna, R. F., & Kulkarni, A. K. 2012, MNRAS, 426, 3241, doi: 10.1111/j.1365-2966.2012.22002.x
  * Narayan & Yi (1994) Narayan, R., & Yi, I. 1994, ApJ, 428, L13, doi: 10.1086/187381
  * Nishiyama et al. (2010) Nishiyama, S., Hatano, H., Tamura, M., et al. 2010, ApJ, 722, L23, doi: 10.1088/2041-8205/722/1/L23
  * Novikov & Thorne (1973) Novikov, I. D., & Thorne, K. S. 1973, in Black Holes (Les Astres Occlus), 343–450
  * Penna et al. (2013) Penna, R. F., Kulkarni, A., & Narayan, R. 2013, A&A, 559, A116, doi: 10.1051/0004-6361/201219666
  * Porth et al. (2019) Porth, O., Chatterjee, K., Narayan, R., et al. 2019, ApJS, 243, 26, doi: 10.3847/1538-4365/ab29fd
  * Remillard & McClintock (2006) Remillard, R. A., & McClintock, J. E. 2006, ARA&A, 44, 49, doi: 10.1146/annurev.astro.44.051905.092532
  * Ressler et al. (2020a) Ressler, S. M., Quataert, E., & Stone, J. M. 2020a, MNRAS, 492, 3272, doi: 10.1093/mnras/stz3605
  * Ressler et al. (2020b) Ressler, S. M., White, C. J., Quataert, E., & Stone, J. M. 2020b, ApJ, 896, L6, doi: 10.3847/2041-8213/ab9532
  * Shakura & Sunyaev (1973) Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 24, 337.  https://ui.adsabs.harvard.edu/abs/1973A%26A....24..337S
  * Silk & Rees (1998) Silk, J., & Rees, M. J. 1998, A&A, 331, L1.  https://arxiv.org/abs/astro-ph/9801013
  * Sorathia et al. (2012) Sorathia, K. A., Reynolds, C. S., Stone, J. M., & Beckwith, K. 2012, ApJ, 749, 189, doi: 10.1088/0004-637X/749/2/189
  * Suzuki & Inutsuka (2014) Suzuki, T. K., & Inutsuka, S.-i. 2014, ApJ, 784, 121, doi: 10.1088/0004-637X/784/2/121
  * Tchekhovskoy et al. (2011) Tchekhovskoy, A., Narayan, R., & McKinney, J. C. 2011, MNRAS, 418, L79, doi: 10.1111/j.1745-3933.2011.01147.x
  * West et al. (2008) West, A. A., Hawley, S. L., Bochanski, J. J., et al. 2008, AJ, 135, 785, doi: 10.1088/0004-6256/135/3/785
  * White et al. (2019) White, C. J., Quataert, E., & Blaes, O. 2019, ApJ, 878, 51, doi: 10.3847/1538-4357/ab089e
  * White et al. (2020) White, C. J., Quataert, E., & Gammie, C. F. 2020, ApJ, 891, 63, doi: 10.3847/1538-4357/ab718e
  * White et al. (2016) White, C. J., Stone, J. M., & Gammie, C. F. 2016, ApJS, 225, 22, doi: 10.3847/0067-0049/225/2/22
  * Yuan & Narayan (2014) Yuan, F., & Narayan, R. 2014, ARA&A, 52, 529, doi: 10.1146/annurev-astro-082812-141003
  * Yuan et al. (2022) Yuan, F., Wang, H., & Yang, H. 2022, ApJ, 924, 124, doi: 10.3847/1538-4357/ac4714
  * Zamaninasab et al. (2014) Zamaninasab, M., Clausen-Brown, E., Savolainen, T., & Tchekhovskoy, A. 2014, Nature, 510, 126, doi: 10.1038/nature13399
  * Zhu & Stone (2018) Zhu, Z., & Stone, J. M. 2018, ApJ, 857, 34, doi: 10.3847/1538-4357/aaafc9