# Magnetorotational dynamo can generate large-scale vertical magnetic fields in 3D GRMHD simulations of accreting black holes

Jonatan Jacquemin-Ide, François Rincon, Alexander Tchekhovskoy, Matthew Liska

<figure><img src="image_1.png" /></figure>

## To cite this version:

Jonatan Jacquemin-Ide, François Rincon, Alexander Tchekhovskoy, Matthew Liska. Magnetorotational dynamo can generate large-scale vertical magnetic fields in 3D GRMHD simulations of accreting black holes. Monthly Notices of the Royal Astronomical Society, 2024, 532 (2), pp.1522-1545. 10.1093/mnras/stae1538. hal-04289057

HAL Id: hal-04289057

https://hal.science/hal-04289057v1

Submitted on 26 Sep 2024

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Magnetorotational dynamo can generate large-scale vertical magnetic fields in 3D GRMHD simulations of accreting black holes

Jonatan Jacquemin-Ide,<fnref n="1" />* François Rincon,<fnref n="2" /> Alexander Tchekhovskoy<fnref n="1" /> and Matthew Liska<fnref n="3" />,<fnref n="4" />

Accepted 2024 June 10. Received 2024 May 16; in original form 2024 January 4

## ABSTRACT

Jetted astrophysical phenomena with black hole engines, including binary mergers, jetted tidal disruption events, and X-ray binaries, require a large-scale vertical magnetic field for efficient jet formation. However, a dynamo mechanism that could generate these crucial large-scale magnetic fields has not been identified and characterized. We have employed three-dimensional global general relativistic magnetohydrodynamical simulations of accretion discs to quantify, for the first time, a dynamo mechanism that generates large-scale magnetic fields. This dynamo mechanism primarily arises from the non-linear evolution of the magnetorotational instability (MRI). In this mechanism, large non-axisymmetric MRI-amplified shearing wave modes, mediated by the axisymmetric azimuthal magnetic field, generate and sustain the large-scale vertical magnetic field through their non-linear interactions. We identify the advection of magnetic loops as a crucial feature, transporting the large-scale vertical magnetic field from the outer regions to the inner regions of the accretion disc. This leads to a larger characteristic size of the, now advected, magnetic field when compared to the local disc height. We characterize the complete dynamo mechanism with two time-scales: one for the local magnetic field generation, $t_{gen}$, and one for the large-scale scale advection, $t_{adv}$. Whereas the dynamo we describe is non-linear, we explore the potential of linear mean field models to replicate its core features. Our findings indicate that traditional $\alpha$-dynamo models, often computed in stratified shearing box simulations, are inadequate and that the effective large-scale dynamics is better described by the shear current effects or stochastic $\alpha$-dynamos.

**Key words:** accretion, accretion discs – dynamo – black hole – neutron star mergers – neutron star mergers – transients: tidal disruption events.

## 1 INTRODUCTION

### 1.1 Compact mergers, disruptions, and binaries

The first direct detection of gravitational waves from a binary neutron star (BNS) merger occurred in 2017 with GW170817. This remarkable event was accompanied by the detection of a gamma-ray burst (GRB) and a kilonova, marking a significant milestone in multimessenger observations (for reviews, see Nakar 2020; Margutti & Chornock 2021). Black hole neutron star (BHNS) mergers also hold promise for multimessenger astronomy, as they are expected to exhibit similar electromagnetic emission components as BNS mergers, including kilonovae and jet emissions (Paczynski 1991; Mochkovitch et al. 1993; Janka et al. 1999; Rosswog 2005; Metzger et al. 2010; Fernández et al. 2015; Foucart et al. 2015; Gompertz et al. 2023).

GRBs are often associated with jet launching from the central engine (Rezzolla et al. 2011; Ruiz et al. 2016), be it a neutron star or a black hole (BH), during BHNS or BNS mergers (see Metzger 2019, for a recent review). The formation of a relativistic jet

during a BHNS merger, contingent on the presence of an accretion disc, depends on various system characteristics such as mass ratio, neutron star radius, and orbital dynamics. However, accretion discs are anticipated to form post-merger for specific mass ratios and non-excessively compact neutron stars (Shibata & Uryū 2006, 2007; Foucart et al. 2011; Foucart 2012; Hayashi et al. 2021; Biscoveanu et al. 2023).

If the jet engine is a BH, it likely operates through the Blandford & Znajek (1977, hereafter BZ) mechanism. This mechanism, in turn, necessitates the presence of a large-scale vertical magnetic field threading the accretion disc. Magnetic fields also play a significant role in driving accretion and mass ejection within post-merger discs (Siegel & Metzger 2017; Christie et al. 2019; Fernández et al. 2019). Large-scale vertical fields are responsible for driving magnetohydrodynamic (MHD) outflows that effectively remove angular momentum and mass (Blandford & Payne 1982; Ferreira & Pelletier 1995), while both small and large-scale magnetic fields facilitate angular momentum transport by generating MHD turbulence through the magnetorotational instability (MRI; Balbus & Hawley 1991). Consequently, understanding the emergence, regeneration, and evolution of magnetic fields is essential for understanding and constraining the dynamics of mergers.

<footnote n="1">¹Center for Interdisciplinary Exploration & Research in Astrophysics (CIERA), Physics and Astronomy, Northwestern University, Evanston, IL 60202, USA</footnote>

<footnote n="2">²Institut de Recherche en Astrophysique et Planétologie (IRAP), Université de Toulouse, CNRS, UPS, F-31400, Toulouse, France</footnote>

<footnote n="3">³Institute for Theory and Computation, Harvard University, 60 Garden Street, Cambridge, MA 02138, USA</footnote>

<footnote n="4">⁴Center for Relativistic Astrophysics, Georgia Institute of Technology, Howey Physics Bldg, 837 State Street NW, Atlanta, GA 30332, USA</footnote>

<footnote>* E-mail: jonatan.jacqueminide@northwestern.edu</footnote>
In the context of BNS mergers, it is highly probable that magnetic
fields undergo significant amplification during the collision phase.
This violent event likely establishes the initial magnetic field config-
uration of the resulting accretion disc. The underlying mechanisms
for this magnetic amplification are thought to primarily involve the
Kelvin–Helmholtz instability (KHI) and magnetic winding processes
(Price & Rosswog 2006; Kiuchi et al. 2015; Aguilera-Miret et al.
2020; Aguilera-Miret, Viganò & Palenzuela 2022; Aguilera-Miret
et al. 2023). Increasingly, it is becoming clear that during the
initial stages of the merger evolution (*t* < 100 ms), a dynamically
important toroidal magnetic field perpendicular to the axis of rotation,
is generated. This phenomenon is observed not only in simulations
of BNS mergers but also in those of BHNS mergers after the initial
disruption phase (Ruiz, Shapiro & Tsokaros 2018; Most et al. 2021).
However, current simulations do not completely capture the magnetic
field amplification dynamics (Kiuchi et al. 2018) and often rely on
arbitrarily large magnetic fields, subgrid models, or prescription for
the magnetic field amplification (Aguilera-Miret et al. 2023; Hayashi
et al. 2023; Most & Quataert 2023; Gottlieb et al. 2023b).

The bulk of research efforts on magnetic field evolution in mergers
has been concentrated on this initial phase. While the generated
toroidal magnetic field can be exceptionally strong, it is incapable of
launching jets. In this study, we aim to investigate the transformation
of an initially toroidal magnetic field into a large-scale vertical
magnetic field that can launch jets.

The question of magnetic field generation is also relevant to
other astrophysical sources with less constrained initial magnetic
topologies. Jetted tidal disruption events (TDEs) share similarities
with BHNS mergers but occur on much larger spatial and temporal
scales. These phenomena have received considerable attention, with
Swift J1644 + 57 serving as a prominent example (Bloom et al. 2011;
Burrows et al. 2011). Tchekhovskoy et al. (2014) suggested that
the broad-spectrum emission, spanning radio, X-rays, and γ-rays,
observed in Swift J1644 + 57, could be explained by a dynamically
significant magnetic flux anchored at the BH event horizon, resulting
in remarkably efficient jet formation. However, the main challenge
arises from the limited magnetic flux available on a stellar object
to facilitate such highly efficient jet launching (Giannios & Metzger
2011; Kelley, Tchekhovskoy & Narayan 2014). Therefore, finding
the mechanism to amplify the magnetic field is imperative. During
the tidal disruption process, the magnetic field of the disrupted star
can undergo significant amplification due to the strong shear, similar
to BHNS mergers. However, without the generation of a large-scale
vertical magnetic field, the resulting magnetic field configuration
remains primarily radial or toroidal, and hence inadequate for jet
launching (Bonnerot et al. 2017). This underscores the critical
importance of generating large-scale vertical magnetic fields in TDEs
to facilitate powerful jet formation to account jetted TDEs like Swift
J1644 + 57.

X-ray binaries undergo impressive outbursts where their luminos-
ity changes by several orders of magnitude, and their spectrum is
significantly modified. At least some of this dramatic evolution is
associated with changes in the accretion disc (Done, Gierliński &
Kubota 2007). During these outbursts, the radio emission, a proxy
for jet activity, is highly correlated with the X-ray luminosity, a
proxy for the accretion rate, linking both processes (Corbel et al.
2003, 2013). It is believed that the magnetic field plays a natural role
in connecting both mechanisms. Furthermore, the radio luminosity
dramatically changes during the outbursts, suggesting that jets turn
on and off. Models explaining this behaviour assume that the vertical
magnetic field is advected and diffused during the outburst, mediating
the change in jet luminosity (Ferreira et al. 2006; Marcel et al. 2019).

A mechanism generating a large-scale vertical magnetic field would
complement those models, as it could provide a source for the large-
scale vertical magnetic field necessary to power the jets.

Finally, recent polarization observations of Sgr A* during a bright
near-infrared flare have emphasized the necessity of a substantial
poloidal magnetic field that has vertical component (Jiménez-Rosales
et al. 2020). This magnetic field component is crucial for matching
the periodicity of the flare. Furthermore, observations from the
Event Horizon Telescope focused on polarized emissions around
the supermassive BH in M87 find better agreement with a strong and
organized poloidal magnetic field (EHT Collaboration 2021).

1.2 The MRI dynamo

An efficient mechanism that generates a large-scale poloidal mag-
netic field from the dominant toroidal magnetic field could be the
MRI. The MRI turbulence was investigated as a dynamo mechanism
soon after its discovery (Brandenburg et al. 1995; Hawley, Gammie &
Balbus 1996). An initially zero net vertical flux magnetic field leads to
self-sustaining extended duration turbulence with an effective angu-
lar momentum transport coefficient, of at most α<sub>v</sub> ~ 0.01 (Shakura &
Sunyaev 1973; Fromang & Papaloizou 2007). The turbulence can be
sustained because the initial MRI-unstable magnetic field can self-
generate and self-sustain (Rincon, Ogilvie & Proctor 2007).

This self-sustaining process has been studied in detail in non-
stratified shearing box simulations (Lesur & Ogilvie 2008b; Herault
et al. 2011; Riols et al. 2013, 2015, 2017; Mamatsashvili et al. 2020;
Held & Mamatsashvili 2022). The self-sustaining process can be
broken down into three main ingredients: (1) First, a large-scale
axisymmetric toroidal field is generated through Keplerian shear
action on a weak large-scale poloidal field. (2) This large-scale
toroidal field is unstable to the non-axisymmetric MRI, driving
perturbations in the form of MHD shearing waves (Goldreich &
Lynden-Bell 1965; Balbus & Hawley 1992; Johnson 2007). (3) Non-
linear wave interactions generate a large-scale axisymmetric poloidal
magnetic field. This large-scale poloidal can then be sheared again
into a toroidal field to restart the process. It must be stressed that
the MRI dynamo is intrinsically nonlinear and very different from
a typical kinematic dynamo; a finite initial seed field is needed to
kick-start the whole process (Rincon 2019).

In stratified shearing boxes, the MRI dynamo self-organizes into
quasi-periodic reversals of the large-scale magnetic field (Bran-
denburg et al. 1995; Davis, Stone & Pessah 2010; Gressel 2010;
Gressel & Pessah 2015; Salvesen et al. 2016). These reversals
propagate from the disc mid-plane to the disc corona and look
qualitatively similar to the butterfly diagram of the solar dynamo.
Recently, the self-sustaining process described above was extensively
characterized in stratified shearing boxes (Held, Mamatsashvili &
Pessah 2024). To corroborate the relevance of the MRI dynamo
for mergers, the butterfly diagrams was also observed in dynamical
general relativistic MHD simulations of BHNS mergers (Hayashi
et al. 2022, 2023).

While there has been some discussion of the dynamo-generating
mechanism operating in global 3D simulations, most of the work has
focused on finding reduced mean-field models (Flock et al. 2012;
Hogg & Reynolds 2018; Dhang et al. 2020). While useful, mean field
models linearize the dynamics and cannot capture the full nonlinear
dynamo physics. Furthermore, it is crucial to understand the physical
mechanism before extrapolating mean field models to astrophysical
regimes.

We hope the analysis we present here can inform the construction
of mean-field models that could help computing merger simulations
with weaker initial magnetic fields. The relevance of this mechanism for mergers and disruptions means that a quantitative analysis of the self-sustaining MRI dynamo mechanism in global simulations is long overdue. Recently, Liska, Tchekhovskoy & Quataert (2020) simulated an accretion disc with an initially toroidal magnetic field and found that a large-scale poloidal field emerged naturally. Then, this large-scale poloidal magnetic field accumulates into the central BH until it reaches the magnetically arrested state (MAD; Igumen-shchev, Narayan & Abramowicz 2003; Narayan, Igumenshchev & Abramowicz 2003; Tchekhovskoy, Narayan & McKinney 2011), in which the vertical magnetic flux is dynamically important. Where MADs can be thought of as the natural end state of magnetized accretion that possesses a large amount of large-scale vertical magnetic flux (Jacquemin-Ide, Lesur & Ferreira 2021).

This manuscript elucidates the self-sustaining MRI dynamo in global 3D general relativistic MHD (GRMHD) simulations of radiatively inefficient accretion flows (RIAFs). This mechanism appears to be a generic feature of thick and highly turbulent MHD accretion discs, making it applicable to a wide range of astrophysical sources. In this study, we specifically concentrate on the mergers of compact objects, as they offer a compelling context for a dynamo mechanism capable of producing large-scale vertical magnetic fields. The manuscript is organized as follows. Section 2 introduces the numerical experiment and our Reynolds averaging method. Section 3 qualitatively describes the main features of the numerical experiment and their time evolution. Section 4 quantitatively describes the magnetic field generation mechanism. Section 5 finds the best mean-field model to describe the dynamo mechanism. Section 6 summarizes and discusses our results.

# 2 METHODS

## 2.1 Numerical set-up

We employ the H-AMR code (Liska et al. 2022) to solve the set of ideal GRMHD equations on a spherical polar grid (r, θ, φ) in the Kerr-Schild coordinates (Gammie, McKinney & Tóth 2003). We also define for convenience the cylindrical radius R and the vertical height z. We adopt dimensionless units such that G = M = c = 1, where M is the mass of the BH. This implies that in our units, the gravitational radius is unity, $r_g = GM/c^2 = 1$. For the magnetic field, we use Lorentz-Heaviside units, such that the magnetic pressure is given by $b^2/2$ in terms of the fluid-frame magnetic field strength, b.

In this work, we reanalyse the simulation first presented by Liska et al. (2020). An accretion disc orbiting a nearly maximally rotating BH (a = 0.9) is simulated. The disc is initialized with a sub-Keplerian Chakrabarti (1985) torus, where the specific angular momentum profile follows $l \propto R^{1/4}$. The torus has an inner edge at $r_{in} = 6r_g$, a pressure maximum at $r_{max} = 13.792r_g$, and an outer edge at $r_{out} \approx 4 \times 10^4 r_g$. The grid extends to $r = 10^5 r_g$ and the inner boundary is causally disconnected from the flow inside the event horizon. We use a polytropic equation of state with $\gamma = 5/3$, which gives a torus scale height of $H/r \sim 0.2$ at $r_{max}$ and $\sim 0.5$ at $r_{out}$. We use transmissive BCs at the poles, periodic BCs in the φ-direction, and outflowing BCs at the inner and outer radial boundaries are used. We use the piecewise parabolic method (Colella & Woodward 1984) for spatial reconstruction and second-order time-stepping. The simulation uses a base grid of resolution $N_r \times N_\theta \times N_\phi = 1872 \times 624 \times 128$ that is uniform in $\log r$, $\theta$, and $\phi$, respectively. On top of the base grid, we use three levels of static mesh refinement (SMR) in the φ-direction: this progressively increases the effective φ-resolution from the pole to equator and

maintains the cell aspect ratio close to unity. This leads to an effective resolution, $N_r \times N_\theta \times N_\phi = 1872 \times 624 \times 1024$, within 60 deg of the equator that corresponds to 70–90 cells per disc scale height. The MRI is seeded with a strong large-scale and uniform toroidal field, with an initial plasma beta $\beta_{ini} = 2P/b^2 = 5$; where P is the fluid-frame gas pressure. We also analyse a supplementary simulation identical to the first, except that the toroidal magnetic field polarity is inverted in the south hemisphere (see Appendix B).

An initially axisymmetric toroidal magnetic field might seem like an unrealistic initial condition. However, as described in Section 1(a) dominant toroidal magnetic field is naturally generated due to azimuthal shear in dynamical space-time GRMHD simulations after the early collision/disruption phase. Furthermore, Aguilera-Miret et al. (2023) find that the strong toroidal magnetic fields organize into axisymmetric structures.

Liska et al. (2020) evolved the simulation out to $t = 120,000 r_g/c$ to verify that the emergent poloidal field was stable on long lived. The analysis presented here focuses solely on the emergence of the poloidal magnetic field, which happens at $t < 10,000 r_g/c$. We record the data snapshots at a cadence of $\Delta t = 10 r_g/c$. At $r = 20 r_g$, this provides $\sim 9$ data snapshots per Keplerian time-scale, $\Omega_K^{-1} = R^{3/2}$, without the $2\pi$.

## 2.2 Averaging procedure

Throughout this manuscript we use a Reynolds decomposition to decompose all quantities as

$$X = \langle X \rangle + \delta X, \quad (1)$$

where $\langle X \rangle$ is the average or large-scale component and $\delta X$ is the turbulent component. For equation (1) to be consistent, the turbulent components must vanish under averaging, $\langle \delta X \rangle = 0$. Only non-linear terms, $\langle \delta X \delta Y \rangle$, do not vanish under averaging due to their correlations. Those correlation terms represent the backreaction of the turbulence on the mean fields.

The Reynolds decomposition is defined here with the azimuthal average,

$$\langle X \rangle = \frac{1}{2\pi} \int_{0}^{2\pi} X d\phi. \quad (2)$$

Consistent with this definition, we will also refer to average quantities as axisymmetrized. The effect of MHD turbulence on large-scale magnetic field generation is perceived by Reynolds averaging the induction equation.

We also define a disc average in the latitudinal coordinate,

$$X|_d = \frac{1}{\int_{\theta_1}^{\theta_2} \sqrt{-g} d\theta} \int_{\theta_1}^{\theta_2} \sqrt{-g} \langle X \rangle d\theta, \quad (3)$$

where g is the determinant of the metric. The disc average will only be applied to already azimuthally averaged quantities ($\langle X \rangle$ or $\langle XY \rangle$). We choose

$$\theta_{1,2} = \frac{\pi}{2} \pm \arctan \left( 3 \frac{h}{R} \right), \quad (4)$$

where h is the disc geometrical thickness, which is $h/R \simeq 0.35$ in the regions of interest. This choice of disc average avoids including the jet funnel and includes most of the turbulent signal. We verify that this choice captures most of the turbulent signal in Section 4. We experimented with other choices and found identical trends.
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Snapshots of the axisymmetrized toroidal magnetic field, $R \langle B_{\phi} \rangle$, in colour as a function of R and z. We also show the poloidal magnetic field lines through the poloidal magnetic flux (equation 9): dashed lines show negative polarity and solid lines show positive polarity. The magnetic field structures become larger and larger with time, reflecting the generation of large-scale poloidal magnetic flux.</figcaption></figure>

We define the poloidal magnetic field as

$$B_p = B_r e_r + B_{\theta} e_{\theta}. \qquad (5)$$

Within this manuscript all quantities will be in physical units unless stated otherwise, thus the polar field above is actually $B_{\theta} = \sqrt{g_{\theta\theta}} B^{\theta}$. Throughout the rest of the manuscript we drop the hats for clarity.

The Reynolds decomposition can be readily performed on the Newtonian induction equation,

$$\frac{\partial B}{\partial t} = \nabla \times (u \times B), \qquad (6)$$

to make the analysis more tractable. As we will see the magnetic field is generated at large distances far from the BH, where this approximation is accurate. We check that energy is conserved, a posteriori, by comparing the left- and right-hand sides of equation (6). If energy is well conserved our Newtonian approximation is accurate at the radii in question and farther away. The average of equation (6) gives

$$\frac{\partial \langle B \rangle}{\partial t} = \nabla \times (\langle u \rangle \times \langle B \rangle + \mathcal{E}), \qquad (7)$$

where the gradient now only depends on r and θ. We define the turbulent electromotive force (EMF)

$$\mathcal{E} = \langle \delta u \times \delta B \rangle. \qquad (8)$$

The turbulent EMF is the feedback of turbulence on the mean magnetic field, it dissipates or generates mean magnetic fields and is at the core of large-scale dynamo theory. For magnetic field generation and sustenance to be possible, there needs to be an energy feedback loop connecting the mean toroidal and poloidal magnetic fields. The turbulent EMF plays a critical role connecting the two magnetic field components.

## 3 DESCRIPTION OF THE SIMULATION

In this section, we set the stage with a description of the qualitative evolution of the magnetic field during the simulation. We then quantitatively characterize the saturated turbulent state of the simulation.

### 3.1 Qualitative evolution of the magnetic field

Fig. 1 shows $R \langle B_{\phi} \rangle$ in colour and the poloidal magnetic field lines through the magnetic potential

$$\Phi(r, \theta, t) = 2\pi \left( \int_{0}^{\theta_{\text{NH}}(\theta)} \langle B^r \rangle \sqrt{-g} \, d\theta - \int_{\pi}^{\theta_{\text{SH}}(\theta)} \langle B^r \rangle \sqrt{-g} \, d\theta \right), \qquad (9)$$

where $B^r$ is the contravariant component of the radial magnetic field, and $\theta_{\text{NH}}$ and $\theta_{\text{SH}}$ are the latitudinal variables in the north and south hemisphere, respectively ($0 > \theta_{\text{NH}} < \pi/2$ and $\pi > \theta_{\text{SH}} > \pi/2$). Fig.
<figure><img src="image_4.png" /></figure>

**Figure 2.** Panel (a): Axisymmetrized poloidal magnetic flux (equation 9), $\Phi$, evaluated at $\theta = \pi/2$ as a function of $t$ and $R$. The vertical dashed lines show the times of the snapshots of Fig. 1. Notice the large-scale advection of magnetic field structures. The dominant large-scale flux loop appears at $R \sim 60r_g$ and $t \sim 4.5 \times 10^3 r_g/c$. The dashed vertical lines show the time at which the snapshots in Fig. 1 are shown. Panel (b): We show the MADness parameter, $\Phi_{BH}/\sqrt{\dot{m}c^2}$ where $\Phi_{BH} = \Phi(r=r_H, \theta=\pi/2, t)$, as a function of time. Notice that the net magnetic field connects to the BH at $t = 10^4 r_g/c$ and the MAD state appears at $t = 2.8 \times 10^4 r_g/c$. Panel (c): We show the poloidal magnetic flux at late times, $\Phi(r, \theta = \pi/2, t = 3 \times 10^4)$, as a function of radii, we notice that there is not net polarity in the system taken as a whole, the net flux emerges only in the inner regions.

2(a) shows a space-time diagram of the magnetic flux through the disc mid-plane, $\Phi(r, \theta = \pi/2, t)$. At $t = 4000 \ r_g/c$, the system has already formed axisymmetric poloidal magnetic fields, $\langle B_p \rangle$, in loop structures (see Fig. 1a). The growth of the axisymmetric poloidal field occurs on a relatively short time-scale.

In Fig. 3, we show the local time-evolution of the different magnetic energy components averaged in the disc (see Section 2.2), evaluated at $r = 20r_g$, as functions of the local shear time. As we will see below, the first large-scale poloidal magnetic loop that connects to the BH is originally advected from $r = 20r_g$; subsequent dominant poloidal magnetic field polarities originate from even farther out.

We can distinguish two clear phases: (1) an initial exponential growth of the turbulent field, $\langle \delta B^2 \rangle$, and the axisymmetric poloidal field, $\langle B_p \rangle^2$. This exponential growth lasts up to $t \simeq 15 \ \Omega_K (r = 20r_g)^{-1}$ when the turbulent magnetic energy, $\langle \delta B^2 \rangle$, reaches saturation. (2) Afterwards, the magnetic fields at the local radius reach a nonlinear statistically steady-state characterized by no net growth and a well-defined average. The axisymmetric poloidal magnetic energy, $\langle B_p \rangle^2$, saturates simultaneously with the turbulent magnetic component, showing their dynamical link. The axisymmetric toroidal magnetic energy, $\langle B_\psi \rangle^2$, decreases once the turbulent magnetic energy becomes comparable to it. The backreaction of turbulence on the seed field then erases the initial condition.

On average, the poloidal loop structures have opposite polarities on opposite sides of the disc mid-plane (see Fig. 1a). However, this trend is only roughly true, and disappears at later times; the system is

<figure><img src="image_5.png" /></figure>

**Figure 3.** Time evolution of different axisymmetrized magnetic-field energies averaged in the disc and evaluated at $r = 20r_g$. We observe a clear initial exponential growth for $\langle B_p \rangle^2$ and $\langle \delta B^2 \rangle$. After the fluctuations become of the same order as for the average toroidal field, $\langle B_\psi \rangle^2$, the axisymmetric toroidal field energy rapidly decreases.

quite stochastic. Within the inner regions the toroidal magnetic field changes signs through shearing of the poloidal magnetic field loops. The topology of the magnetic field is analysed in greater detail in Section 4.6.

We notice that at high latitudes above the disc mid-plane the poloidal field forms elongated structures in the form of escaping flux tubes. These are reminiscent of the magnetic towers described by Lynden-Bell (2003, 2006) and show the vertical shedding of magnetic-field structures. Lynden-Bell (2006) show that a poloidal magnetic field line anchored at two distinct radii ($r_1$ and $r_2$, where $r_2 > r_1$) will undergo vertical expansion due to differential shear. This expansion is driven by the differential rotation introducing a twist in the field line (owing to $\Omega(r_1) > \Omega(r_2)$), which results in the formation of a strong toroidal field capable of buoyantly rising and piercing through the surrounding medium. Consequently, this process forms vertically elongated poloidal magnetic flux structures known as magnetic towers, illustrated in Fig. 1(a,b) (see also fig. 2 of Lynden-Bell 2006). We anticipate that similar ballooning effects will occur for field lines connecting the BH to the disc, analogous to phenomena observed in interactions between stellar dipoles and discs (Lynden-Bell 2003; Zanni & Ferreira 2013).

For $t < 4000 \ r_g/c$, there is no coherent poloidal field connected to the central BH but coherent poloidal loops are still present (Fig. 2). At $t \sim 4000 \ r_g/c$, a coherent poloidal structure attaches to the BH. This coherent poloidal loop originates far from the BH, at $r \sim 20r_g$. We show a zoom in on the axisymmetrized poloidal and toroidal field structures connected to the BH in Fig. 4. The topology of how this magnetic field connects to the BH is distinct from the one observed in magnetically arrested discs (MAD; Tchekhovskoy et al. 2011). Fig. 4 shows open and closed field lines, with the closed field lines threading the disc mid plane. This structure drives jets at low-energy efficiencies (Christie et al. 2019; Liska et al. 2020; Gottlieb
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Axisymmetrized toroidal magnetic field, $R\langle B_\psi \rangle$ in colour as a function of $R$ and $z$ at $t = 8 \times 10^3 r_g/c$. We also show the poloidal magnetic field lines through the poloidal magnetic flux (equation 9), dashed lines show negative polarity and solid lines show positive polarity. In the initial stages of evolution, $t < 10^4 r_g/c$, the inner magnetic field structure is distinct from what is typically seen in MADs.</figcaption></figure>

et al. 2023b). A similar structure was studied in 2D GR particle in cell simulations by El Mellah et al. (2022), who showed that it can drive efficient particle acceleration. The fundamental reason for this different field structure is that the magnetic field loop connected to the BH has a smaller size, leading to less magnetic flux opening up through the action of the light cylinder (see Section 4.5 and Chashkina, Bromberg & Levinson 2021).

The system continues to evolve, for $t > 4000 r_g/c$, and forms loops of larger and larger size (see the movie and Fig. 1b). A large-scale loop emerges at $R \sim 50 r_g$ and $t \sim 5000 r_g/c$ (see Fig. 2) and at $t \sim 8000 r_g/c$, this poloidal field structure is threading the disc mid-plane (Fig. 1b). The ballooning instability (Lynden-Bell 2003) can then take hold on both hemispheres. This instability expands the poloidal loop, which pushes away all other magnetic loops, leading to a dominant polarity that connects to the BH (see Fig. 1c). Fig. 2(b) depicts the Madness parameter, $\Phi_{BH}/\sqrt{\dot{m}c^2}$, where $\Phi_{BH} = \Phi(r = r_H, \theta = \pi/2, t)$, as a function of time. We note that although a large-scale magnetic field structure connects to the BH at $t = 10^4$, $r_g/c$, the system only reaches the MAD state at $t \sim 2.8 \times 10^4 r_g/c$. The polarity emerging at $t \sim 10^4 r_g/c$ then persists for the remainder of the simulation and is still present at $t \sim 1.5 \times 10^5 r_g/c$ (Liska et al. 2020). In this study, we will focus on the emergence of this large-scale polarity, $t < 10^4 r_g/c$.

<figure><img src="image_7.png" /><figcaption>**Figure 5.** Panel (a): $\theta-t$ diagram of the axisymmetrized toroidal magnetic field, $\langle B_\psi \rangle$, evaluated at $r = 20 r_g$, normalized to the initial toroidal field. We observed the ejection of inverting toroidal structures, reminiscent of the typical butterfly diagram of MRI turbulence. However, the structures in our simulation are more stochastic than those typically seen in MRI turbulence. Panel (b): $r-t$ diagram for the axisymmetrized toroidal magnetic field, $\langle B_\psi \rangle$, evaluated at the surface of the disc, $\theta = \pi/2 + \arctan(3h/R)$, normalized to the initial field. Large-scale (as large as $\sim 70 r_g$) toroidal field structures of differing polarity are ejected from the system. There are inversions in the polarity of the field, but no clear periodicity is observed.</figcaption></figure>

Even though the ballooning instability can naturally explain the runaway growth of a magnetic field loop whose vertical length scale, $l_{loop}$, is larger than the disc scale height, $h$, it cannot explain how a poloidal magnetic field larger than the disc scale height appears in the first place. Furthermore, it is surprising that the loop that connects to the BH at $t \sim 4000 r_g/c$ does not become the dominant polarity, while the one that connects at $t \sim 10^4 r_g/c$ does. In Section 4, we discuss the mechanism that generates the poloidal magnetic field and how the dominant polarity takes hold.

## 3.2 Butterfly diagram

Fig. 5(b) shows a $\theta-t$ diagram, also called the butterfly diagram, of the toroidal magnetic field at $r = 20 r_g$. The toroidal magnetic field is ejected from the system and the ejection happens in cycles with polarity inversions. Shearing box simulations show similar flux inversions with a vertical propagation (Gressel & Pessah 2015). However, periodic cycles in shearing box simulations are far more regular than those observed in Fig. 5(a). This is consistent with
Hogg & Reynolds (2018), who showed that thicker discs have a more incoherent dynamo cycle.

In Fig. 5(a), we show an r-t diagram for the toroidal magnetic field evaluated at the surface of the disc, $\theta = \pi/2 + \arctan(3h/R)$. We see magnetic flux inversions, as in the butterfly diagram. However, in this case the inversions show a certain degree of coherence for multiple radii. We believe that these eruptions of magnetic field are due to the ballooning instability discussed above. As the poloidal magnetic field loops grow, they can shed larger loop structures vertically. This leads to the radially coherent field structures seen in Fig. 5(b).

### 3.3 Turbulent steady state of the system

In the quasi-statistically steady state, we have an energetic ordering, $\langle \delta B^2 \rangle \gg \langle B_\phi \rangle^2 \gg \langle B_p \rangle^2$, with about an order of magnitude difference between each component (see Fig. 3). The turbulent magnetic energy saturates at roughly the energy of the initial condition. We define different plasma betas for every energetic component

$$ \beta_{\delta} = \frac{8\pi \langle P \rangle}{\langle \delta B^2 \rangle}, \quad (10) $$

$$ \beta_p = \frac{8\pi \langle P \rangle}{\langle B_p \rangle^2}, \quad (11) $$

$$ \beta_{\phi} = \frac{8\pi \langle P \rangle}{\langle B_{\phi} \rangle^2}, \quad (12) $$

and the average turbulent stress

$$ \alpha_{\nu} = - \frac{\langle \delta B_{\phi} \delta B_{p} \rangle}{4\pi \langle P \rangle}, \quad (13) $$

where we only compute the Maxwell stress because the Reynolds stress is negligible in our simulations. In Fig. 6, we show the inverse of the different plasma betas and the turbulent stress vertically averaged in the disc (equation 3) and temporally averaged between $t = 5000 \ r_g/c$ and $t = 8000 \ r_g/c$. The vertical average is computed after squaring the magnetic field or else opposite polarity fields would cancel out, leading to an artificially small measurement (see Fig. 1). The shaded region highlights the regions that have not converged and do not follow the same trend as the inner radii.

We measure $\beta_p \simeq 300$ and $\beta_\phi \simeq 40$. The toroidal field lost an order of magnitude in magnetic energy; most of that energy was converted into or dissipated by the turbulent fluctuations, which saturate at a $\beta_\delta \simeq 4$, of the order of $\beta_{ini} = 5$. At large radii, $r > 50r_g$, the different components trend towards their initial values ($\beta_\delta^{-1} \to 0$, $\beta_p^{-1} \to 0$, $\beta_\phi^{-1} \to \beta_{ini}^{-1}$) as the dynamo has not reached saturation at those radii.

We measure an $\alpha_\nu \simeq 0.1$ a factor of 10 larger than in typical zero-net flux shearing box simulations (Hawley et al. 1996). We attribute this large value to the initially strong toroidal field. Hawley et al. (1995) measured the following scalings for fully developed MRI turbulence in shearing boxes with an initially toroidal field (see also Section 8.1.1 in Lesur 2020)

$$ \alpha_{H95} \simeq 0.51 \beta_{\delta}^{-1}, \quad (14) $$

$$ \alpha_{H95} \simeq 0.35 \beta_{\text{ini}}^{-1/2}. \quad (15) $$

For our initial condition equation (15) gives $\alpha_{H95} \simeq 0.16$ which is remarkably close to the value we compute. Equation (14) implies a factor of 2 difference between $\beta_\delta^{-1}$ and $\alpha_\nu$ that can be read from Fig. 6. Overall, we find that the saturated turbulent configuration of our simulation is consistent with scalings derived from shearing box simulations of MRI. The large value of $\alpha_\nu$ could instead be attributed to the generated axisymmetric poloidal magnetic field driving its own

<figure><img src="image_8.png" /><figcaption>Figure 6. Inverse plasma betas for the different magnetic field components, $\beta_\delta^{-1}$ for the turbulent field, $\beta_p^{-1}$ for the axisymmetrized poloidal field, $\beta_\phi^{-1}$ for the axisymmetrized toroidal field, and the angular momentum transport $\alpha_\nu$ coefficient. They are vertically averaged in the disc and temporally averaged between $t = 5000 \ r_g/c$ and $t = 8000 \ r_g/c$. The shaded regions highlights the non-saturated outer radii, $r > 70r_g$. For $r < 70r_g$ the plasma betas are relatively constant and have reached a saturated quasi-steady state. The vertical dashed line highlights the point where $\langle u_r \rangle = 0$, as shown in Fig. 10(b).</figcaption></figure>

axisymmetric MRI, as proposed by Bai & Stone (2013) to explain the large $\alpha_\nu$ in global zero net flux simulations. Indeed, the turbulent transport driven by axisymmetric MRI and mediated by a vertical magnetic field is much more efficient than non-axisymmetric MRI mediated by a toroidal field. For a vertical magnetic field, the scaling for $\alpha_\nu$ yields (Salvesen et al. 2016; Lesur 2020)

$$ \alpha_{S16} \simeq 10 \beta_p^{-1/2} \sim 0.14, \quad (16) $$

which is also consistent with the magnitude of the turbulent transport displayed in Fig. 6. The fact that $\alpha_{H95}(\beta_{ini}) \simeq \alpha_{S16}(\beta_p)$ is remarkable and may be an interesting feature to investigate for lower initial fields.

## 4 ANALYSIS OF THE DYNAMO MECHANISM

We now proceed to analyse the dynamical mechanism underlying the magnetic evolution described above.

### 4.1 Toroidal magnetic field budget

To better track the dynamics and sources of magnetic energy we perform a Reynolds average decomposition of the magnetic energy induction equation. First, we dot equation (8) with the average toroidal field to find

$$ \frac{1}{2} \frac{\partial \langle B_{\phi} \rangle^2}{\partial t} = S_{\phi} + A_{\phi} + \delta A_{\phi}, \quad (17) $$

where

$$ S_{\phi} = \langle \mathbf{B}_{\phi} \rangle \cdot \nabla \times (\langle \mathbf{u}_{\phi} \rangle \times \langle \mathbf{B}_{p} \rangle), \quad (18) $$
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Different contributions to the toroidal magnetic energy equation (equation 17), averaged within the disc and evaluated at $r = 20r_g$, and are normalized to the initial magnetic energy and the Keplerian disc frequency at that radius. Toroidal magnetic energy is only generated through shear of the poloidal field, $S_\phi$. The turbulent component, $\delta A_\phi$, always dissipates the toroidal field through turbulent mixing, while the contribution of the, poloidal fluid motions, $A_\phi$, is to eject the toroidal field from the disc. That the black and grey dotted lines agree, validates the accuracy of our Newtonian approximation and Reynolds decomposition (equation 17).</figcaption></figure>

is the stretching of poloidal into toroidal field by the large-scale shear (Ω-effect),

$$A_{\phi} = \langle \mathbf{B}_{\phi} \rangle \cdot \nabla \times (\langle \mathbf{u}_p \rangle \times \langle \mathbf{B}_{\phi} \rangle), \quad (19)$$

describes the large-scale axisymmetric advection of toroidal field, and

$$\delta A_{\phi} = \langle \mathbf{B}_{\phi} \rangle \cdot \nabla \times (\mathcal{E}_p), \quad (20)$$

captures the turbulent generation or dissipation, depending on the sign, of axisymmetric magnetic field.

We show the time dependence of the different, vertically averaged, components of equation (17) evaluated at $r = 20r_g$ in Fig. 7 (the different terms in Fig. 7 are vertically averaged within the disc and normalized to the initial toroidal magnetic field and the Keplerian frequency). The right- and left-hand sides of equation (17) are observed to be equal, validating the use of the Newtonian approximation and showing that our Reynolds decomposition is accurate.

We notice that after $10\Omega_k^{-1}(r = 20r_g)$ the system reaches a quasi-statistically steady state, consistent with Fig. 3. The only component that produces toroidal magnetic energy is the shear, as in shearing box simulations of the MRI dynamo (Lesur & Ogilvie 2008b; Riols et al. 2017). The turbulent term, $\delta A_\phi$, always dissipates toroidal magnetic energy. Shearing box simulations also find that turbulence always dissipates toroidal magnetic energy via turbulent mixing. The advective term, absent in shearing boxes, transports toroidal magnetic energy outwards.

To better understand the vertical profile of the magnetic field source terms, we compute the time average of the different components of equation (17) and evaluate them at $r = 20$. We average them between $20\Omega_k^{-1}(r = 20r_g)$ and $110\Omega_k^{-1}(r = 20r_g)$ and show them as functions of $\theta$ in Fig. 8, normalized to the initial magnetic energy and Keplerian frequency. The turbulent magnetic component is only important on the disc mid-plane and tends to 0 outwards. The shaded regions in Fig. 8 show the bounds of the disc average defined in Section 2.2. We see that our choice of integral bounds includes most of the turbulent signal.

Overall, the Ω-effect compensates for the dissipation due to turbulence and outward advection present within the body of the disc. Above the disc, in the shaded region, shear and advection balance each other.

<figure><img src="image_10.png" /><figcaption>**Figure 8.** Different contribution to the toroidal magnetic energy equation (equation 17) averaged between $20\Omega_k^{-1}(r = 20r_g)$ and $110\Omega_k^{-1}(r = 20r_g)$ and evaluated at $r = 20r_g$, and normalized to the initial magnetic energy and the Keplerian disc frequency at that radius. The shaded regions show $|\theta - \pi/2| > \arctan(3h/R)$. The toroidal field loses energy to the advection term, $A_\phi$, vertically escaping toroidal fields, and the turbulent dissipation, acting within the bulk of the disc. Both terms are cumulatively balanced by the Ω effect, $S_\phi$.</figcaption></figure>

## 4.2 Poloidal magnetic field budget

We have found that the sustainment of toroidal field is only possible thanks to the shearing of a large-scale poloidal field. Thus, the dynamics and generation of the poloidal field is essential to the dynamo mechanism as a whole. In particular, since the only term that produces toroidal magnetic energy is the Ω-effect, poloidal field generation is necessary for toroidal field production above the disc.
We derive an equation for the energetics of the axisymmetrized poloidal field

$$ \frac{1}{2} \frac{\partial \langle B_p \rangle^2}{\partial t} = A_p + \delta A_p, \qquad (21) $$

where

$$ A_p = \langle \mathbf{B}_p \rangle \cdot \nabla \times (\langle \mathbf{u}_p \rangle \times \langle \mathbf{B}_p \rangle) \qquad (22) $$

describes the large-scale or axisymmetric advection of poloidal field, and

$$ \delta A_p = \langle \mathbf{B}_p \rangle \cdot \nabla \times (\mathcal{E}_\varphi \mathbf{e}_\varphi), \qquad (23) $$

represents the turbulent generation or dissipation of poloidal magnetic field.

In Fig. 9(a), we show the different components of equation (21). We average all components within the disc and evaluate them at $r = 20r_g$. We normalize them to the initial magnetic energy and the Keplerian disc frequency at that radius. We notice that our Newtonian Reynolds decomposition is also accurate for the poloidal magnetic energy equation. At $r = 20r_g$, the large-scale advection of the magnetic field is positive, i.e. it locally brings poloidal magnetic energy. We explore below how this quantity varies with radius. We show the turbulent EMF term in Fig. 9(b). This turbulent term is incoherent; it stochastically changes sign, producing and dissipating large-scale magnetic energy. To understand the mechanism generating the magnetic field, we must go a step further, and make a Fourier decomposition of the different modes contributing to the turbulent magnetic production/dissipation. Using Parseval's theorem, we decompose the turbulent term as follows:

$$ \delta A_p = \sum_{m=1}^{\infty} \delta A_p^{(m)}, \qquad (24) $$

where

$$ \delta A_p^{(m)} = \langle \mathbf{B}_p \rangle \cdot \nabla \times (\mathcal{E}_\varphi^{(m)} \mathbf{e}_\varphi), \qquad (25) $$

and

$$ \mathcal{E}_\varphi^{(m)} \mathbf{e}_\varphi = \mathcal{R} [\mathbf{u}_p^m \times \mathbf{B}_p^m]^*, \qquad (26) $$

$\mathcal{R}$ represents taking the real part, $X|^*$ is the complex conjugate of $X$, and $\mathbf{u}_p^m$ and $\mathbf{B}_p^m$ are the spectral coefficients of the $m$-mode in the fourier expansion of $\mathbf{u}_p$ and $\mathbf{B}_p$,

$$ \mathbf{u}_p = \sum_{m=0}^{\infty} \mathbf{u}_p^m e^{im\varphi} \qquad (27) $$

$$ \mathbf{B}_p = \sum_{m=0}^{\infty} \mathbf{B}_p^m e^{im\varphi}. \qquad (28) $$

Note that $\mathcal{E}_\varphi^{(m)}$ represent the axisymmetric effect of nonlinear interaction of non-axisymmetric $m$-modes.<fnref n="1" /> The $^{(m)}$superscript denotes nonlinear-correlations on the mean ($\varphi$-averaged) EMF due to non-zero $m$ fluctuations of $\mathbf{u}_p$ and $\mathbf{B}_p$. This notation should not be confused with the superscript $^m$ (no parenthesis) denoting the spectral coefficients of the fields $u$ and $B$ in the Fourier expansions (equations (27, 28)). We stress that the Fourier modes $\mathbf{B}_p^m$, $\mathbf{u}_p^m$ and their non-linear correlation $\mathcal{E}_\varphi^{(m)}$ are solely functions of $r$ and $\theta$. When referring to the complete $m$-mode non-axisymmetric motion, we will explicitly include the $\varphi$ dependence as $\mathbf{u}_p^m e^{im\varphi}$.

To avoid double counting of the $m = 0$ advection term, $A_p$, we only sum modes with $m \ge 1$ in equation (24). This approach yields a definition of $\delta A_p$ that aligns perfectly with equation (23).

This procedure is inspired by work done in shearing box simulation, where only the large-scale, low $m$ modes contribute to poloidal field production (Lesur & Ogilvie 2008b; Riols et al. 2015). Here too we measure that only the correlations of the largest scale modes, $m = 1, 2, 3$, generate $\varphi$-average poloidal magnetic energy at all radii; we refer to these modes as active (see also Appendix A). All correlations of modes with $m > 3$ dissipate $\varphi$-averaged poloidal magnetic energy on average, we refer to them as turbulent dissipative or passive. In Fig. 9(b), we show the sum of active, magnetic field-generating, $\sum_{m=1}^3 \delta A_p^{(m)}$ versus dissipative/passive modes, $\sum_{m=4}^{\infty} \delta A_p^{(m)}$.

The sum of active modes always generates axisymmetric poloidal magnetic fields. The passive modes have a more complicated behaviour, generating magnetic energy for $t < 20\Omega_k^{-1}$ but dissipating it for $t > 20\Omega_k^{-1}$. The time-scale when this term goes from generating to dissipating is equal to the saturation time scale shown in Fig. 3. Therefore, we infer that this early behaviour is related to the initial linear growth phase of such modes, while the later one is associated with their dissipative-action through a non-linear turbulent cascade. Due to the short duration of this phase, it does not affect the general trends; the non-linear effect of $m > 3$ modes is to dissipate the energy of the axisymmetric poloidal field. Furthermore, for $t < 20\Omega_k^{-1}$ the active component ($m < 3$) is comparable or larger than the turbulent passive component. Therefore, we conclude that the active modes generate basically all the axisymmetric poloidal magnetic energy.

We have also computed the mode decomposition for the toroidal field. In Appendix A, we show that the nonlinear effect of all modes with $m \ge 1$ is to dissipate $\varphi$-averaged toroidal magnetic energy. In particular, the $\varphi$-averaged toroidal field loses energy to the non-linear correlations of active modes, $m = 1, 2, 3$. We also show in Appendix A that the active modes transport angular momentum. We can therefore interpret the active modes as MRI-driven perturbations mediated by the axisymmetrized toroidal field.

We now look at the large-scale global behaviour of the poloidal energy equation, equation (21), by vertically and temporally averaging it. The temporal average is performed between $5000r_g/c$ and $8000r_g/c$. Fig. 10(a) shows the vertical average of the different terms of equation (21) divided by the vertical average of $\langle B_p \rangle^2$ and normalized to the local Keplerian frequency. Note that here, we show the opposite of the magnetic energy dissipation term (e.g. multiplied by $-1$). The three terms have similar magnitudes, $0.1\Omega_K^{-1}$, and are roughly constant for $r < 30r_g$ but follow different trends at large radii. First, let us note that in the inner regions of the disc, $r < 50$, their signs are consistent with Fig. 9. The magnetic generation term is roughly constant for the inner radii, only increasing for large radii, $r > 70$. We attribute this increase to a lack of saturation in the outer regions, leading to a lower average axisymmetric poloidal magnetic energy (see Fig. 6). By dividing by the axisymmetric poloidal magnetic energy, we can interpret the magnitude of the quantity, $2\sum_{m=1}^3 \delta A_p^{(m)}/\langle B_p \rangle^2$, as the typical frequency of magnetic field regeneration. We define a regeneration time scale,

$$ t_{\text{gen}} = \frac{\frac{1}{2} \langle B_p \rangle^2}{\frac{3}{\sum_{m=1}^{3} \delta A_p^{(m)}}} \approx 10\Omega_K^{-1}. \qquad (29) $$

This time-scale is consistent with Fig. 3 up to a factor of $\sim 2$, it is also roughly consistent with the cycles of Fig. 5. The time-

<footnote n="1">¹In practice, we never compute the Fourier modes of the EMF. We only extract the Fourier modes of the velocity and magnetic fields to calculate $\mathcal{E}_\varphi^{(m)}$ in equation (26).</footnote>
<figure><img src="image_11.png" /><figcaption>**Figure 9.** Panel (a): Different contributions to the poloidal magnetic energy equation (equation 21) averaged within the disc and evaluated at $r = 20r_g$, and normalized to the initial magnetic energy and the Keplerian disc frequency at that radius. Panel (b): The sum, $\delta A_p$, (in red) of magnetic field-producing (or active, shown in blue), $\sum_{m=1}^{3} \delta A_p^{(m)}$, and dissipative (or passive, shown in orange), $\sum_{m=4}^{\infty} \delta A_p^{(m)}$, modes. The sum of active modes is always positive, while the passive term is mostly negative. Large-scale non-axisymmetric structures, $m = 1, 2, 3$, therefore generate axisymmetric poloidal magnetic fields through their non-linear self-interactions.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Figure 10.** Panel (a): Vertically and temporally averaged poloidal energy equation (equation 21), normalized to the local Keplerian frequency. The temporal average is performed between $5000r_g/c$ and $8000r_g/c$. The different terms of equation (21) are divided by the vertical and temporal average of $\langle B_p \rangle^2$. The opposite of the magnetic energy dissipation term is shown. Advection, $A_p$, and dissipation, $\sum_{m=4}^{\infty} \delta A_p^{(m)}$, of axisymmetric poloidal field balance each other in the inner regions. The vertical dashed line highlights the point, where $\langle u_r \rangle = 0$. Panel (b): Vertically and temporally averaged radial velocity $\langle u_r \rangle / V_K$, notice that the trends of $A_p$ are well explained by the behaviour of the radial velocity (see the text). We also show an affine fit of the normalized radial velocity. Finally, we show the magnitude of the magnetic field generating velocity, $u_{sw}$.</figcaption></figure>

scale, $t_{gen}$, should not be interpreted as a growth rate, as it possesses more information than a growth rate, and encapsulates the saturation energy of the axisymmetric poloidal field.

The large-scale axisymmetric poloidal advection term is positive in the inner regions of the disc, consistent with the inward advection of the poloidal magnetic field (Fig. 10a). We have checked that the

advection term, $A_p$, is dominated by radial advection and radial velocity divergence. In Fig. 10(b), we show the disc and temporal average (between $5000r_g/c$ and $8000r_g/c$) of the radial velocity, $\langle u_r \rangle$. This velocity is negative in the inner regions, consistent with inward advection.
At the distance, $r_{\text{rev}} \sim 50$, the radial velocity changes sign. The inversion of the radial velocity is a fossil feature of our initial condition. In standard accretion theory, the radial velocity first moves outwards in the outer regions and then changes sign after a few viscous time-scales (Lynden-Bell & Pringle 1974). Hence, this feature is related to the outer regions not being in inflow equilibrium. The radii at which we have $\langle u_r \rangle = 0$, viscously moves out over the course of our simulation, starting at $r_{\text{rev}} = 13r_g$ and moving out to $r_{\text{rev}} \sim 50$ by $t \sim 6 \times 10^3 r_g/c$.

At $r \simeq r_{\text{rev}} \sim 50$, the advection term, $A_p$, begins to trend towards negative values. As expected, the advective flux depends on the radial velocity. The advection term is negative for radii $50 < r < 100$, even though the radial velocity is positive at those radii. This is because the radial dependency of the radial velocity avoids a pile-up of the magnetic field in those regions. Magnetic field bundles farther away from $r_{\text{rev}}$ travel faster than bundles closer to $r_{\text{rev}}$, leading to a chase. The chase of magnetic bundles leads to a net depletion of the magnetic field at those radii, $50 < r < 100$. In the outer regions ($r > 10^2$), the radial velocity can become quite large, with $u_r \sim 0.1V_K$, resulting in efficient ejection of axisymmetric poloidal magnetic fields. This feature will be revisited in Section 4.4.

In the inner regions, $r < 50$, the opposite of the turbulent dissipation term ($-\sum_{4}^{\infty} \delta A_p^{(m)}$) closely tracks the large-scale advection term. The two terms roughly balance each other, indicating a statistically steady state of advection and diffusion. This equilibrium does not imply that the axisymmetric poloidal magnetic field is static, but rather that any reconnecting or escaping field is constantly replenished by advection and the local regeneration of axisymmetric poloidal magnetic field.

Turbulent dissipation plays a different role for small-scale magnetic field loops than it does for a large-scale vertical magnetic field. For a large-scale vertical magnetic field, local dissipation in the disc can only diffuse the field towards larger radii (Lubow, Papaloizou & Pringle 1994; Guilet & Ogilvie 2012). In contrast, for magnetic field loops smaller than the disc, dissipation acts as a local reconnecting agent.

In the outer regions, $r > 50$, the turbulent term $\sum_{4}^{\infty} \delta A_p^{(m)}$ is positive; this sign change is attributed to the linear growth shown in Fig. 9(b). During the linear growth stage, before the local saturation time-scale and a turbulent cascade develops, the 'passive' modes generate magnetic energy. However, it is clear that the outer regions are not yet in a statistical steady state at this stage for $r > 50$ (see Fig. 6).

## 4.3 Shearing wave structures

As shown above, the dynamo mechanism relies on non-linear correlations of non-zero-m-modes, that are the result of shearing wave structures, and it is thus interesting to analyse these structures in more detail. Fig. 11 shows the poloidal velocity streamlines, $\mathcal{R}(\mathbf{u}_p^m=2e^{i2\varphi})$, of the $m = 2$ mode at $\phi = \frac{4}{3}\pi$ and $t = 6.5 \times 10^3 r_g/c$ superimposed on the axisymmetrized cylindrical radial field, $\langle B_R \rangle$. Notice that the $m = 2$ mode takes the form of rolling structures in the poloidal plane. These structures should not be interpreted as axisymmetric turbulent eddies; they are the rolling poloidal motions of traveling non-axisymmetric shearing waves. The centres of the rolling structures are well correlated with the places where the radial magnetic field changes sign. We interpret this rolling motion of the shearing modes as the field-generating action (Herault et al. 2011). Fig. 11 can also be compared to Fig. 7 of Riols et al. (2017), who found similar poloidal plane structures for the field-generating waves in local shearing box simulations (without the complexity of a global

<figure><img src="image_13.png" /><figcaption>**Figure 11.** Poloidal velocity streamlines, $\mathcal{R}(\mathbf{u}_p^m=2e^{i2\varphi})$, of the $m = 2$ mode at $\phi = \frac{4}{3}\pi$ and $t = 6.5 \times 10^3 r_g/c$ superimposed on the axisymmetrized cylindrical radial field, $\langle B_R \rangle$. They take the form of rolling structures, whose centres are well correlated with the regions where $\langle B_R \rangle = 0$. Poloidal magnetic-field generating motions thus takes the form of spiral rolling motions, which should not to be interpreted as axisymmetric eddies, as visualized here.</figcaption></figure>

structure). They also found a correlation between the centre of the wave structures and the radial magnetic field.

In the previous section, we estimated the regeneration time scale of the poloidal field as a function of radius, $t_{\text{gen}} = 10\Omega_k^{-1}(R)$. To show the self-consistency of our analysis, we compute $t_{\text{gen}}$ using the amplitude of the wave velocities. We use the following estimate,

$$u_{\text{sw}} = \sum_{m=1}^{m=3} \sqrt{\langle (\mathcal{R}(\mathbf{u}_p^m e^{im\varphi})|_d)^2 \rangle}, \quad (30)$$

where $X|_d$ represents the disc average defined in Section 2.2, $\mathcal{R}$ represents taking the real value. We square $\mathbf{u}_p^m e^{im\varphi}$ before taking its azimuthal average because, $\langle \mathbf{u}_p^m e^{im\varphi} \rangle = 0$, due to it being a non-axisymmetric quantity.

Fig. 10(b) shows the estimate for the average wave amplitude, $u_{\text{sw}}$, normalized to the local Keplerian velocity and time averaged between $t = 5 \times 10^3 r_g/c$ and $t = 8 \times 10^3 r_g/c$ as function of radius.
We see that in the inner regions, which are in quasi-steady state,<fnref n="2" />

$$u_{sw} \approx 2.3 \times 10^{-2} V_K(r). \quad (31)$$

If we assume that the wave structures generate magnetic fields on length scales the size of half the disc vertical thickness (see Fig. 1a), we get

$$t_{sw} = \frac{h/2}{u_{sw}} \sim 10\Omega_k^{-1}(r). \quad (32)$$

This time-scale is completely consistent with $t_{gen}$ computed in the previous section. This consistency is satisfying, as we can predict the field regeneration rate based on the amplitude of the waves. A computation using the vorticity,

$$\omega_{sw} = \frac{1}{2} \sum_{m=1}^{m=3} \sqrt{\langle \nabla \times (\mathcal{R}(u_p^m e^{im\varphi})) \rangle_d^2} \approx 0.2\Omega_k \quad (33)$$

of the wave structures yields a similar time-scale, $t_{vort} = 1/\omega_{sw} \approx 20\Omega_k^{-1}$. However, this calculation does not show the role of the length scale in the dynamics.

By choosing $h/2$ as the length scale, we have illuminated a subtle issue. In Section 3, we described that for one polarity to take hold, one magnetic field loop needs to be larger than the size of the disc scale height. Hence, if every magnetic field loop generated is of maximal size $\sim h(r)/2$ there is no way for one polarity to take hold. A secondary mechanism is required to enhance the size of magnetic field structures. In Section 4.5, we argue that the advection serves this purpose.

## 4.4 Field transport

As demonstrated earlier, advection plays a crucial role in shaping the large-scale dynamics of the poloidal magnetic field by transporting magnetic flux from outer to inner regions. In this section, we provide a quantitative analysis of the velocities and time-scales associated with this field transport, elucidating their interaction with local magnetic field regeneration.

A simple affine fit for the radial velocity gives,

$$\frac{\langle u_r \rangle}{V_K(r)} = a_1 r - a_2, \quad (34)$$

where $a_1 = 7.8 \times 10^{-4}$ and $a_2 = 4.1 \times 10^{-2}$ (see Fig. 10b). The value of $a_2$ is around a factor of 4 larger than the value predicted by standard accretion theory<fnref n="3" />, $-\frac{\langle u_r \rangle}{V_K} \approx \alpha_v (\frac{h}{r})^2 \sim 10^{-2}$. We extrapolate equation (34) into two simpler limits for the radial velocity

$$\langle u_r \rangle \approx a_1 r V_K \quad \text{for } r \gg r_{\text{rev}}, \quad (35)$$

$$\langle u_r \rangle \approx -a_2 V_K \quad \text{for } r \ll r_{\text{rev}}. \quad (36)$$

Following Jacquemin-Ide et al. (2021), we use the expressions above to compute the trajectory of magnetic field lines. To do this we solve for the characteristics of field transport, given by $dR_\Phi/dt = \langle u_r \rangle$.

This leads to the following trajectories

$$t_{\text{eje}} = \frac{2}{a_1} (R_{\Phi}^{1/2} - R_0^{1/2}) + t_0 \quad \text{for } r \gg r_{\text{rev}}, \quad (37)$$

$$t_{\text{adv}} = -\frac{2}{3a_2} (R_{\Phi}^{3/2} - R_0^{3/2}) + t_0 \quad \text{for } r \ll r_{\text{rev}}, \quad (38)$$

<figure><img src="image_14.png" /><figcaption>**Figure 12.** Axisymmetrized poloidal magnetic field flux (equation 9), $\Phi$, evaluated at $\theta = \pi/2$, as a function of $t$ and $R$. We also show the growth time-scale of the magnetic-field loop structures, $t_{gen} = 10\Omega_k^{-1}$ (equation 29) and the magnetic field advection and ejection trajectories (see the text). The largest poloidal magnetic field structures appear along the growth time-scale line, then clearly follow the advection and ejection trajectories.</figcaption></figure>

where $t_{adv}$ represents the poloidal field advection time-scale, $t_{eje}$ is the poloidal field ejection time-scale, $R_\Phi$ is the anchoring radii of the field line at time $t_{adv}$ (or $t_{eje}$) and $R_0$ is the initial anchoring radii at time $t_0$. Using the expressions above, we can compute the transport of the large-scale loop that emerges as the dominant polarity at $t \sim 10^4 r_g/c$. In Fig. 12, we show an $r-t$ diagram for the poloidal magnetic flux (equation 9), same as Fig. 2. We show the regeneration time-scale computed in equation (29) on Fig. 10, $t_{gen} \approx 10\Omega_k^{-1}$. The large-scale poloidal field structures appear along the line following the regeneration time-scale; $t_{gen}$ is a good estimate for the growth of large-scale structures. After axisymmetric poloidal field loop forms, it advects either inwards or outwards.

The ejection of poloidal magnetic flux explains the existence of magnetic field structures at large radii within time-scales $t < t_{gen}$. Specifically, at large radii, $t_{eje} < t_{gen}$, poloidal field structures originate at smaller radii before being ejected outwards.

We show the trajectories of inwards and outwards advecting loops in Fig. 12. We compute the trajectory of the inwards (outwards) advecting loop for $R_0 = 50r_g$ ($R_0 = 70r_g$) and $t_0 = 4.5 \times 10^3 r_g/c$. It is clear that our simple model, equations (37)–(38), reproduces the large-scale advection of the magnetic field. Furthermore, magnetic field loops at earlier times, for $t < 4 \times 10^4 r_g/c$, follow the same kind of trajectories. We find that the flux is constantly advected at a velocity that follows a Keplerian scaling, in the inner regions, and a velocity that follows $\langle u_r \rangle \propto R^{1/2}$, in the outer regions. We conclude that the evolution of the axisymmetric poloidal flux and the emergence of the large-scale field are well described by the two time-scales $t_{gen}$ and $t_{adv}$.

<footnote n="2">²Other estimates for $u_{sw}$ lead to similar values, for example $u_{sw} = |\mathcal{E}_\Phi|_d / \langle B_p \rangle^2 |_d \sim 10^{-2} V_K$.</footnote>

<footnote n="3">³We note that this approximation ignores the radial gradient of $R^2 \langle P \rangle$.</footnote>
## 4.5 Advection and size of poloidal field structures

As the poloidal magnetic field is generated within the disc, the size of the magnetic field structures will be limited to the local disc scale height, ~h(r)/2. A simple way to build field structures larger than the local scale height without invoking a reverse cascade is provided by the large-scale transport, $A_p$. Within the disc the poloidal magnetic field is advected towards the inner regions by the large-scale accretion flow, $\langle u_r \rangle$. Hence, if a loop is born at $r \simeq 60 r_g$ it needs to get advected to $r \simeq 30 r_g$ to be larger than the disc thickness. We can see this happen in Fig. 12, a magnetic field loop generated at $r \simeq 60 r_g$ slowly gets advected towards $r \simeq 20 r_g$, and then experiences runaway growth by the ballooning instability (see Section 3 and Lynden-Bell 2003, 2006). A additional complication is that, if there are two or more loops of different polarities, they will compete over which can be advected. This competition will add some stochasticity to the evolution of the system as imbalances determine the dominant magnetic field polarity. The reason why the dominant polarity emerges at $r = 60r_g$, instead of $r < 60r_g$, is unclear, but we propose three complementary arguments:

(1) Magnetic field loops born before $t \sim t_{gen}(60r_g) \sim 4.5 \times 10^4 r_g/c$ do not have enough runway to be advected before being attached to the BH (see Fig. 4). If a loop appears at $r = 30r_g$ it is larger than the disc scale height at $r \sim 15r_g$. However, at this distance the loop might already be too close to the ISCO and is thus doomed to connect to the BH before becoming larger. Thus there is a critical radius where the large-scale field is large enough to experience the ballooning instability without being trapped by the BH.

(2) We define the poloidal magnetic flux that the MRI dynamo can generate, $\Phi_{loop} \simeq \langle B_p \rangle l_{loop}^2$, where $l_{loop}$ is the typical size of the loop. We have

$$ \langle B_p \rangle = \sqrt{8\pi\beta_p^{-1}} \sqrt{\langle P \rangle} \propto r^{-1}, \quad (39) $$

since $\langle P \rangle \propto r^{-2}$, and with $l_{loop} \sim h/2$ we find that

$$ \Phi_{loop} \simeq \sqrt{8\pi\beta_p^{-1}} \sqrt{\langle P \rangle} (h/2)^2 \propto r, \quad (40) $$

an increasing function of radius. Thus, the farther from the BH the generated loop is, the larger the magnetic flux it contains. While equation (40) is an approximation, we have verified its effectiveness in describing the time-averaged flux.<fnref n="4" /> From this it is clear that there will be a critical radius that depends on the initial magnetic field profile, where the magnetic flux generated will be sufficient to saturate to a MAD. The previous argument implies that while a large-scale magnetic field (larger than the disc scale height) might be inevitable, the maximal magnetic flux that can be generated might depend on the subtle details like disc scale-height, radial extent of the disc and initial magnetic field strength. The consequences of the maximal magnetic flux depending on disc geometry is discussed in Section 6.2.

(3) Using GRMHD simulations initialized with multiple poloidal field loops, Chashkina et al. (2021) argue that for the jet to have significant power, the advected poloidal loops need to be larger than ~10$r_g$ so that a sufficient amount of magnetic flux opens due to the action of the light cylinder. We compute that for poloidal field loops of size $l_{loop} \sim h/2$, a loop of size ~10$r_g$ must be generated at $r \sim 60r_g$, which is consistent with the description above and Fig. 12. Thus, the analysis of Chashkina et al. (2021) corroborates the existence of a critical poloidal field loop size (or equivalent poloidal flux strength) at which the system transitions into the MAD state.

For our simulation, the large-scale field then advects into the BH until it reaches the MAD state (Tchekhovskoy et al. 2011). This advection phase can be seen in Fig. 2(b) for $t > 10^4 r_g/c$. The dynamo mechanism then probably shuts off once the large-scale poloidal magnetic field is large enough (Salvesen et al. 2016).

## 4.6 Advection of large-scale magnetic field symmetries

To illustrate the impact of advection on the large-scale magnetic field structure more clearly, we compute the polarity of the poloidal and toroidal magnetic fields. Following the approach outlined by Flock et al. (2012), we perform a latitudinal averaging of the magnetic field components in the Southern hemisphere, denoted as $B_i^{NH}$, and the Northern hemisphere, denoted as $B_i^{SH}$. Subsequently, we define the symmetric component, $B_i^S = \frac{1}{2}(B_i^{NH} + B_i^{SH})$, and the antisymmetric component, $B_i^{AS} = \frac{1}{2}(B_i^{NH} - B_i^{SH})$, contributing to the magnetic field structure. Dipolar and quadrupolar magnetic energies can then be constructed for both the toroidal and poloidal magnetic fields:

$$ E_p^D = (B_r^{AS})^2 + (B_\theta^{S})^2 \quad (41) $$

$$ E_\psi^D = (B_\psi^{AS})^2, \quad (42) $$

$$ E_p^Q = (B_r^S)^2 + (B_\theta^{AS})^2, \quad (43) $$

$$ E_\psi^Q = (B_\psi^S)^2. \quad (44) $$

The dipolar configuration, with its energy denoted as $E^D$, exhibits the following topology: a vertical magnetic field that has the same sign above and below the disc mid-plane, in contrast to the radial and toroidal fields which display anti-symmetry. This configuration resembles a MAD (Tchekhovskoy et al. 2011) and also self-similar solutions (Ferreira & Pelletier 1995).

On the other hand, the quadrupolar configuration, with its energy labelled as $E^Q$, represents the converse arrangement. It features an antisymmetric latitudinal field, while $B_r$ and $B_\psi$ exhibit symmetric configurations. This configuration can be visualized as comprising two poloidal field loops with opposite polarities located, respectively, on the Northern and Southern hemispheres of the disc.

Finally, we define the parity of the poloidal and toroidal magnetic field

$$ P_i = \frac{E_i^D - E_i^Q}{E_i^D + E_i^Q}. \quad (45) $$

If $P_i = 1$ the configuration is dipolar, while if $P_i = -1$ the configuration is quadrupolar.

Fig. 13(a) shows the parity, $P_p$, of the poloidal magnetic field as a function of radius and time. Similar to Fig. 12, we display the field generating time-scale, $t_{gen}$, derived previously, along with the advection ($t_{adv}$) and ejection ($t_{eje}$) trajectories. At large radii and short time-scales, we observe that the magnetic field structures emerge in a quadrupolar configuration, as depicted in Fig. 1(a), which illustrates pairs of poloidal field loops with opposite polarities above and below the disc mid-plane. Fig. 13(a) indicates that the inner regions swiftly, $t \simeq 2.5 \times 10^3 r_g/c$, transition to a dipolar topology (also observed in Figs A3 and 1b), resulting from the advection of larger-scale structures into the inner regions. Consistent with the findings above, we conclude that the initial parity, and thus the initial vertical structure of the poloidal magnetic field, is broken by large-scale transport.

At $t = 10^4 r_g/c$ within the inner regions, a larger poloidal field structure is advected from the outer regions. This briefly results in a

<footnote n="4">⁴Where the time average allows to reduce the contribution from advection.</footnote>
<figure><img src="image_15.png" /><figcaption>**Figure 13.** Panel (a): Polarity of the poloidal magnetic field (equation 45) as a function of $t$ and $R$. Quadrupolar topology is associated with $\mathcal{P}_p = -1$, while dipolar topology is associated with $\mathcal{P}_p = 1$. We also show the growth time-scale of the magnetic-field loop structures, $t_{\text{gen}} = 10\Omega_k^{-1}$ (equation 29) and the magnetic field advection and ejection trajectories (see Section 4.4 and Fig. 12). The large-scale symmetries of the poloidal field are determined by its large-scale advection. Panel (b): Polarity of the toroidal magnetic field (equation 45) as a function of $t$ and $R$. While the toroidal polarity shows no clear signs of advection, the polarity of the toroidal field changes once the large-scale poloidal fields are advected inwards.</figcaption></figure>

quadrupolar configuration as the newly advected larger loop and the old smaller loop co-exist. However, eventually, the smaller loop is displaced, returning the configuration to a dipolar state (see Fig. 1c). At large radii, quadrupolar field structures emerge faster than the generation time-scale $t_{\text{gen}}$. This phenomenon is elaborated in Section 4.4, stemming from the ejection of these structures from the inner regions to the outer scales. They travel outward at a pace exceeding the generation time-scale, as indicated by $t_{\text{eje}} < t_{\text{gen}}$ for large $R$.

Fig. 13(b) shows an analogous figure for the toroidal field. The parity of the toroidal field, $\mathcal{P}_{\varphi}$, exhibits a different evolution compared to the poloidal field. Our initial condition for the toroidal magnetic field is quadrupolar, which is evident at large radii and short times. The evolution of, $\mathcal{P}_{\varphi}$, contrary to $\mathcal{P}_p$ shows no clear signs of advection. It is clear that different dynamics command the larger scales of the toroidal field, one of the being the shear and the buoyancy shown in Fig. 5. The lack of inward advection for the toroidal magnetic field was also shown in Fig. 7 (see also Fig. Ab), where the advection term, $\mathcal{A}_{\varphi}$, is shown to be negative.

For $t < 10^4$, we observe a combination of dipolar and quadrupolar configurations in the toroidal field, with rapid oscillations between the two. We attribute these oscillations to the large-scale oscillatory patterns depicted in Fig. 5(b), which are caused by the escape of large-scale toroidal flux from the accretion disc.

For $t > 10^4$, the toroidal field in the inner regions abruptly shifts to a radially coherent dipolar configuration. We attribute this change

to the advection of the large-scale poloidal magnetic loop occurring at that moment. Given that the large-scale toroidal field is generated by the poloidal field, it follows its evolution. Therefore, even though the toroidal field exhibits a distinct evolution without clear inward advection, its final state is ultimately determined by the large-scale inward advection of the poloidal field.

## 5 POSSIBLE DYNAMO CLOSURES

Mean-field models are useful tools to simplify the dynamics of complex systems by linearizing them. However, the MRI dynamo described in Section 4 is inherently non-linear, so mean field models may be inadequate for accurately describing its dynamics. None the less, we attempt here to distinguish which of two common mean field models may be the most accurate for our non-linear dynamics. In this section, we focus on specific radii, but the trends described have been verified to be independent of this choice.

The toroidal magnetic field is only generated by the $\Omega$-effect, while the turbulence only acts as a dissipating agent of it and could be modelled as a turbulent resistivity. In contrast, the poloidal magnetic field can be generated by turbulence and requires a more careful analysis. As we showed in Section 4, in statistical steady state only the non-linear correlations of the large-scale MRI-unstable waves, $m = 1, 2, 3$, generate magnetic fields, while non-linear correlations of smaller scale modes dissipate the magnetic field. Thus, we only
attempt to find a mean-field model for

$$ \mathcal{E}_{\psi}^{\text{sw}} = \sum_{m=1}^{3} \mathcal{E}_{\psi}^{(m)}. \qquad (46) $$

Focusing on large-scale field structures helps filter out turbulent mixing and small-scale dynamos, aiding on the construction of the meanfield model. Similarly, recent shearing box computations use frequency filtering to remove small-scale dynamos, reducing contamination (Dhang, Bendre & Subramanian 2024).

## 5.1 Alpha versus shear-current effects

There are two main competing mechanisms for mean-field models of the MRI dynamo. First, the classical *alpha effect* that has the following form,

$$ \mathcal{E}_{\psi}^{\text{sw}} = \alpha_d \langle B_{\psi} \rangle, \qquad (47) $$

where $\alpha_d$ is the alpha effect. This effect has been measured in various stratified shearing box simulations (Brandenburg et al. 1995; Gressel 2010; Gressel & Pessah 2015) and global simulations (Flock et al. 2012; Hogg & Reynolds 2018; Dhang et al. 2020). Although all authors agree that the sign changes between the upper and lower hemispheres, there is no clear consensus on which sign should be measured in which hemisphere. The rough maximum magnitude measured in simulations is roughly consistent in the literature, $\alpha_d \sim \pm 5 \times 10^{-3} h \Omega_K$.

An alternative version of the $\alpha_d$ effect is the *stochastica* effect (Moffatt 1978). In this model, the mean $\bar{\alpha}_{st} = 0$ but large fluctuations in $\sqrt{\bar{\alpha}_{st}^2}$ can lead to magnetic-field growth (Brandenburg et al. 2008), where $\bar{X}$ is a temporal average of $X$. This effect involves fluctuations in the mean field itself, which are thought to naturally occur when there is not a large-scale separation between the mean-field generating structures and the mean field (Brandenburg 2018). Heinemann, McWilliams & Schekochihin (2011) showed that a collection of shearing waves could excite a stochastic dynamo and generate large-scale magnetic field. Although the dynamo they computed is completely kinematic, and there is no backreaction of the magnetic field on the flow, the main features are reminiscent of the dynamics of our simulation. The fact that the dynamo is dominated by large-scale shearing waves might also explain its stochasticity.

The shear current effect is another possible model for the MRI dynamo, it is expressed as

$$ \mathcal{E}_{\psi}^{\text{sw}} = -\eta_{\text{sc}} \langle J_r \rangle, \qquad (48) $$

where $\langle J_r \rangle$ is the electric current.<fnref n="5" /> This effect was identified analytically (Lesur & Ogilvie 2008a) and also measured in non-stratified shearing boxes (Lesur & Ogilvie 2008b; Squire & Bhattacharjee 2015; Shi, Stone & Huang 2016; Squire & Bhattacharjee 2016) and statistical models of MRI turbulence (Mondal & Bhat 2023). Crucially, $\eta_{sc} < 0$ is a necessary condition for this dynamo to operate (Squire & Bhattacharjee 2015, 2016; Rincon 2019).

Lesur & Ogilvie (2008b) found that the shear current effect requires a time delay, $\tau$, to model the dynamo cycle. This time delay is understood as the time it takes for the mean field to excite waves that will amplify the magnetic field: they assumed $\tau \sim \Omega_k^{-1}$.

To determine the model that fits our non-linear dynamo, we introduce the following time-delayed moving-averaged Pearson

correlation function:

$$ C_p(X, Y(t-\tau)) = \frac{\int_{t-\delta t}^{t} X(t')Y(t'-\tau) dt'}{\sqrt{\int_{t-\delta t}^{t} X^2(t') dt' \int_{t-\delta t}^{t} Y^2(t'-\tau) dt'}} \qquad (49) $$

where $\tau$ is the time delay and $\delta t = n\Omega_K^{-1} + \tau$ is the averaging window.<fnref n="6" /> We choose $n = 4$ to ensure that we average over a few Keplerian shear time-scales, we find no perceptible differences for different values of $n$ as long as $1 < n < 10$.

We show the time-averaged Pearson correlations of equation (49) in Fig. 14(a), between $t = 20\Omega_K^{-1}(r = 20r_g)$ and $t = 100\Omega_K^{-1}(r = 20r_g)$, for $r = 20r_g$. In practice, the correlation function $C_p$ is computed for all $\theta$ and then averaged within the Northern hemisphere of the disc, between $\theta_1 = \frac{\pi}{2} - \arctan(\frac{h}{R})$ and $\theta_2 = \frac{\pi}{2} - \arctan(0.1\frac{h}{R})$. The order of operations is crucial because computing the correlation on the $\theta$-averaged data could lead to error. We first compute the correlation functions for all angles $\theta$ and then average it. We show the correlations of both models equations (47) and (48) and the correlation of $\langle B_r \rangle$ and $\langle B_{\psi} \rangle$, which we use as a baseline of comparison.

We see that the correlation of the radial and toroidal fields peaks at $\tau = 0$ with a maximum value of $-0.5$. The position of the peak at $\tau = 0$ is consistent with zero-lag shear, which is what is expected from the induction equation, as the $\Omega$-effect does not include any lag. However, the maximum amplitude of only $-0.5$ shows how difficult it is to find correlations in global GRMHD simulations, a perfect correlation would lead to value of $C_p = \pm 1$. Other global GRMHD simulations also find correlations of the order of 0.3, although for neutron star dynamos (Kiuchi et al. 2024).

The correlation for a shear-current-effect-like prescription (equation 48) shows a peak at $\tau = 0.5\Omega_K^{-1}(r = 20r_g)$, consistent with Lesur & Ogilvie (2008b). In other words, the $\mathcal{E}_{\psi}^{\text{sw}}$ reacts to the radial current with a time delay $\tau$ of the order of the MRI growth time-scale. The maximal value of $C_p(\mathcal{E}_{\psi}^{\text{sw}}(t), \langle J_r \rangle(t-\tau))$ is 0.25, a marginally significant correlation. However, the correlation for an $\alpha_d$ effect prescription (equation 47) is far smaller, showing a maximum amplitude of only $-0.05$, and it also shows a maximum at a time lag of $\tau = \Omega_K^{-1}(r = 20r_g)$.

Fig. 14(b) shows the time variability of the Pearson correlation functions, vertically averaged in the north hemisphere of the disc for $r = 20r_g$ for their respective $\tau$ maximum measured in Fig. 14(a):

(i) shear: $\tau = 0$ for $C_p(\langle B_r \rangle(t), \langle B_{\psi} \rangle(t-\tau))$ (ii) shear current: $\tau = 0.5\Omega_K^{-1}(r = 20r_g)$ for $C_p(\mathcal{E}_{\psi}^{\text{sw}}(t), \langle J_r \rangle(t-\tau))$

(iii) $\alpha$-effect: $\tau = 1.1\Omega_K^{-1}(r = 20r_g)$ for $C_p(\mathcal{E}_{\psi}^{\text{sw}}(t), \langle B_{\psi} \rangle(t-\tau))$

All correlations are large, with $C_p > 0.4$, and show strong variability. We notice that even though the average correlation between the mean toroidal field $\langle B_{\psi} \rangle$ and the electromotive force $\mathcal{E}_{\psi}^{\text{sw}}$, $C_p(\mathcal{E}_{\psi}^{\text{sw}}(t), \langle B_{\psi} \rangle(t-\tau))$, is small, it can reach high correlation values on shorter time-scales. This large deviation from the average could be a sign of a stochastic $\alpha_{st}$ dynamo effect. We find identical trends for the Northern hemisphere of the disc. For the sake of conscience, we show the $t-\tau$ diagram in Appendix C (see Fig. C1).

<footnote n="6">⁶The inclusion of $\tau$ in $\delta t$ is to ensure that our average window is always larger than the time shift.</footnote>

<footnote n="5">⁵We compute $\langle J \rangle$ in its Newtonian approximation, neglecting the displacement current, a valid approximation at $r = 20r_g$.</footnote>
<figure><img src="image_16.png" /><figcaption>**Figure 14.** Panel (a): Time and vertically averaged Pearson correlations at $r = 20r_g$. The averages are computed between $t = 20\Omega_K^{-1}(r = 20)$ and $t = 100\Omega_K^{-1}(r = 20)$, and between $\theta_1 = \pi/2 - \arctan(\frac{\pi}{R})$ and $\theta_2 = \pi/2 - \arctan(0.1\frac{\pi}{R})$. The correlations are shown for both models, equations (47) and (48), and for the $\Omega$-effect which is used as a baseline of comparison. Only the $\alpha_d$-effect does not show a significant time average correlation. [Panel (b)]: Time variability of the Pearson correlation function, vertically averaged in the north hemisphere of the disc for $r = 20r_g$. For every correlation, we show the $t$ dependence at the $\tau$ maximum measured in panel (a). Notice that on short time-scales the $\alpha_d$ effect is comparable to the shear-current effect.</figcaption></figure>

## 5.2 Mean-field EMFs

For large enough correlations, it makes sense to approximate the $\alpha_{st}$ and $\eta_{sc}$ coefficient with

$$ \alpha_{st} \approx \frac{\int_{t-\delta t}^{t} \mathcal{E}_{\varphi}^{\text{sw}}(t') \langle B_{\varphi} \rangle (t' - \tau) dt'}{\int_{t-\delta t}^{t} \langle B_{\varphi} \rangle^2(t') dt'} \quad (50) $$

$$ \eta_{sc} \approx \frac{\int_{t-\delta t}^{t} \mathcal{E}_{\varphi}^{\text{sw}}(t') \langle J_r \rangle (t' - \tau) dt'}{\int_{t-\delta t}^{t} \langle J_r \rangle^2(t') dt'} \quad (51) $$

Due to variability of $C_p$, this equation only makes sense for small $\delta t > \tau$. We choose $\delta t = 5\Omega_K^{-1}(r = 20r_g)$. We use the same values of $\tau$ as the ones described above (see also Fig. 14b).

Fig. 15(a) shows the time evolution of the full $\mathcal{E}_{\varphi}^{\text{sw}}$, compared to the two different models. With $\alpha_{st}$ and $\eta_{sc}$ coefficients computed from equations (50) and (51), respectively. We compute $\alpha_{st}$ and $\eta_{sc}$ as functions of $\theta$ and $t$ before computing latitudinal or temporal averages. Both models roughly reproduce $\mathcal{E}_{\varphi}^{\text{sw}}$; the $\alpha_{st}$ model appears to do a better job at reproducing the solid curve. We tried different values of $\delta t$: the fit deteriorates for larger values of $\delta t$, but the $\alpha_{st}$ model is always a marginally better fit to the full $\mathcal{E}_{\varphi}^{\text{sw}}$.

Fig. 15(b) shows the $\theta$-time diagram of $\alpha_{st}$ evaluated at $r = 20r_g$. As expected, the alpha coefficient is highly variable and does not show a clear top down symmetry, in contrast with other works (Flock et al. 2012; Hogg & Reynolds 2018). We do not show the temporal evolution $\eta_{sc}$, as it shows similar variability to the stochastic $\alpha_{st}$.

Instead, we show in Fig. 16 the temporal average ($\bar{\alpha}_{st}$, $\bar{\eta}_{sc}$) and rms value ($\sqrt{\bar{\alpha}_{st}^2 + \bar{\eta}_{sc}^2}$) of both coefficients, computed between $t_1 = 20$ and $t_2 = 110\Omega_K^{-1}$. We find no clear trends for the polar profile of $\bar{\alpha}_{st}$, although there might be a slight change of sign around the disc mid-plane. This signal is likely overwhelmed by the large fluctuations, as the rms value is around a factor of 5 larger. We compute a rms

<figure><img src="image_17.png" /><figcaption>**Figure 15.** Panel (a): Time evolution of the full $\mathcal{E}_{\varphi}^{\text{sw}}$, compared to mean-field EMFs using $\alpha_{st}$ and $\eta_{sc}$; coefficients computed from equations (50) and (51), respectively. We have latitudinally averaged $\mathcal{E}_{\varphi}^{\text{sw}}$ and its models over the Northern hemisphere of the disc at $r = 20r_g$. Both model roughly reproduce the EMF with a slight edge for the stochastic alpha model. [Panel (b)]: $\theta - t$ diagram of the stochastic alpha, $\alpha_{st}$, computed using equation (50) normalized to $h\Omega_k$. The coefficient $\alpha_{st}$ is temporally incoherent and lacks the antisymmetry observed in most stratified shearing box simulations.</figcaption></figure>
<figure><img src="image_18.png" /><figcaption>**Figure 16.** Time-averaged $\bar{\alpha}_{st}$ and $\bar{\eta}_{sc}$ coefficients as functions of the polar coordinate, normalized to $h\Omega_k$ and $h^2\Omega_k$, respectively, and evaluated at $r = 20r_g$. We also show their respective rms values $\sqrt{\bar{\alpha}_{st}^2}$ and $\sqrt{\bar{\eta}_{sc}^2}$. For the sake of visibility, the shear current coefficients are multiplied by 5. The coefficients are computed using equations (50) and (51), respectively. No clear trend is observed in the polar profile of $\bar{\alpha}_{st}$ and its rms value is about a factor of 2 larger. The shear-current prescription, $\bar{\eta}_{sc}$, shows a clearer trend hitting a maximum just above the disc mid-plane.</figcaption></figure>

value for the stochastic $\alpha$ prescription of $\sqrt{\bar{\alpha}_{st}^2} \simeq 10^{-3}h\Omega_k$, which is roughly constant in $\theta$.

Our computation of $\bar{\alpha}_{st}$ is around an order of magnitude smaller than fiducial values found in the literature in stratified shearing boxes (Brandenburg et al. 1995; Gressel 2010; Gressel & Pessah 2015; Dhang et al. 2024) or global simulations (Flock et al. 2012; Hogg & Reynolds 2018; Dhang et al. 2020).

As far as the shear-current effect is concerned, we compute a clearly negative shear-current coefficient $\bar{\eta}_{sc} \simeq -2 \times 10^{-4}h\Omega_k^2$ at its maximal amplitude; the sign of the $\bar{\eta}_{sc}$ coefficient thus allows for magnetic field amplification (Rincon 2019). In contrast with $\bar{\alpha}_{st}$, the shear current effect shows a slight trend in $\theta$, reaching its maximum slightly above the disc mid-plane. Furthermore, in this case the maximal value is of the same order as the rms fluctuations. Our computed value of $\bar{\eta}_{sc}$ aligns with findings from non-stratified shearing box simulations (Lesur & Ogilvie 2008b; Shi et al. 2016). However, it is inconsistent with results from stratified shearing box simulations, which find positive shear current effects (Gressel 2010; Gressel & Pessah 2015; Dhang et al. 2024).

From the previous analysis, it is difficult to conclude which model is the better fit for the MRI dynamo. While the stochastic alpha model better fits the behaviour of the EMF when compared the shear current effect, it has a worse average correlation. However, the stochastic alpha model is more intuitively appealing as it is reminiscent of the wave phenomenology detailed in Section 4 (Heinemann et al. 2011).

The other main conclusion of this mean-field analysis is that the antisymmetric (with respect to the disc mid-plane) alpha coefficient computed in stratified shearing box simulations is incompatible with our 3D global simulations (Brandenburg et al. 1995; Gressel 2010). There are two possible explanations for this discrepancy: (1) Stratified shearing boxes are known to have too many symmetries (Lesur, Ferreira & Ogilvie 2013; Lesur 2020). The additional symmetries of shearing boxes could lead to the organization of the large-scale wave structures, which could lead to a more coherent $\alpha$ effect. Furthermore, both global simulations where a coherent alpha effect was measured had a limited latitudinal extent (Flock et al. 2012; Hogg & Reynolds

2018). This limited latitudinal extent could lead to a similar effect of organizing the wave structures. (2) A thinner disc could also lead to a more coherent $\alpha$ effect by vertically confining the wave structures. Indeed, Hogg & Reynolds (2018) found that thicker discs had more irregular dynamo cycles. Future work will focus on how this dynamo changes in thinner discs.

Finally, Dhang et al. (2020) used singular value decomposition to measure mean field dynamo coefficients in RIAFs. Although they use a different initial field, the overall structure looks similar to ours before $t = 8 \times 10^3 r_g/c$. They were only able to fit $\alpha_d$ coefficient at the surface of the disc, and not within the interior of the disc, and attributed this lack of convergence to the disc being too turbulent. This observation might be related to the intrinsic variability that we show in Fig. 15(b).

# 6 CONCLUSIONS AND DISCUSSION

## 6.1 Main conclusions

Binary mergers are theorized to drive jets as a component of the GRB and afterglow emission. If the central engine is a BH, the jet will be driven by the Blandford & Znajek (1977) process. However, the BZ mechanism requires a large-scale poloidal magnetic field. In this paper, we analysed simulations of GRMHD accretion disc to investigate the dynamical mechanism underlying the generation, transport and self-organization of the magnetic field. We have here isolated a mechanism that can generate the large-scale poloidal magnetic fields required for jet launching.

We elucidated that, excluding large-scale transport, the production and dissipation of magnetic energy relies on the following ingredients: (1) The toroidal magnetic field is generated through the shear of the poloidal field (Fig. 7). (2) Large-scale non-axisymmetric MRI wave structures of wavenumber, $m = 1, 2, 3$, non-linearly generate the poloidal field (Fig. 9). (3) The toroidal field loses energy to those non-axisymmetric field fluctuations (Fig. 7 and Appendix A). (4) Finally, both field components lose energy to the small-scale, $m > 3$, non-axisymmetric structures through turbulent dissipation. This behaviour is highly reminiscent of previous results obtained in simpler shearing box simulations (Lesur & Ogilvie 2008b; Herault et al. 2011; Riols et al. 2015, 2017).

As described in the introduction, this dynamo mechanism is interpreted as an interplay between the non-axisymmetric MRI and large-scale axisymmetric dynamics of the magnetorotational instability. The low-$m$ modes are (non-)linear non-axisymmetric MRI modes that feed on the shear and develop on the axisymmetric toroidal magnetic field. The low $m$-modes take the form of MRI-amplified shearing waves (see Fig. 11), which are sheared spiral waves (Goldreich & Lynden-Bell 1965; Johnson 2007). In the statistical steady state, the large-$m$ modes ($m > 4$) are non-linearly generated by the turbulent cascade and mix the large-scale field, dissipating it all small scales. We computed a local growth time scale for the magnetic-field-generation mechanism and found $t_{gen} \simeq 10\Omega_k^{-1}$. The velocity of the spiral modes, $u_{sw}$, can be used to predict this saturation time-scale, $t_{gen} \sim t_{sw} = h/u_{sw}$, with good precision (see Fig. 10b). Hence, magnetic field loops are generated locally at a specific radii with a frequency 10 times slower than the Keplerian frequency.

We also found that large-scale advection of the magnetic field, poloidal or toroidal, is an essential feature of the full magnetic-field dynamics. We showed that the toroidal magnetic field vertically escapes the accretion disc (Figs 5 and 8). Most importantly, the poloidal magnetic field is advected inwards by the accretion flow. The
advection serves as a scale amplification mechanism by transporting the larger-scale loops, generated in the outer regions of the accretion disc, towards the inner regions, with a smaller disc scale height. Thus, advection can transport a magnetic field larger than the local disc scale height. We computed the advection time-scale, $t_{\text{adv}}$, over which this happens using a simple fit for the accretion flow velocity (see equation 38). We have determined that this straightforward calculation of the advection time-scale is remarkably consistent with magnetic field transport and, as a result, possesses predictive capabilities. When putting everything together we can formulate a description for the emergence of a dominant large-scale vertical magnetic field using only two time-scales, $t_{\text{gen}}$ and $t_{\text{adv}}$. We have also approximated the magnitude and profile of the magnetic flux and found that it depends mostly on the geometry of the disc.

We have corroborated the role of advection in transporting the large-scale fields by examining the transport of polarity in the axisymmetrized poloidal and toroidal magnetic fields. We demonstrate that advection, acting on the time-scale $t_{\text{adv}}$, plays a crucial role in breaking the initial quadrupolar symmetry of the system and generating a dipolar symmetry by advecting the large-scale poloidal field structures from the outer regions (see Fig. 13). Furthermore, although the toroidal field is not advected by the accretion flow, its evolution is ultimately dictated by the advection of the poloidal field. Indeed, the polarity of the toroidal field is controlled by the shear of the poloidal field, thus once the poloidal field reaches a dipolar symmetry, the toroidal field swiftly follows.

Finally, even though the dynamo described here is completely non-linear, for the purpose of global modelling and better closures, we attempted to determine which linear mean-field model is capable of reproducing the main features of our simulations. We found that a simple antisymmetric $\alpha$ dynamo model, often computed in stratified shearing box simulations, cannot reproduce the dynamo mechanism observed in our simulation. However, a shear current effect or a stochastic $\alpha$ dynamo could reproduce some of the main features. We believe this discrepancy with shearing box simulations is due to the thickness of our accretion disc. Thinner disc, which are well modelled by shearing boxes, might confine and organize the field-generating structures, and the organization of such structures might lead to a more coherent $\alpha$ effect.

## 6.2 Discussion

### 6.2.1 Dependence on initial field geometry and strength

In this work, we have only considered an initially axisymmetric toroidal magnetic field. It is not entirely clear how sensitive our results are to these initial conditions. However, as noted earlier in Section 2.1, an approximately axisymmetric toroidal field seems to emerge in the initial phase of evolution of BHNS mergers, $t < 100$ ms (Aguilera-Miret et al. 2023). This is encouraging, as this axisymmetric toroidal field plays a key mediating role in MRI mechanism described in this work. Furthermore, Gottlieb et al. (2023b) found that their non-axisymmetric initial condition, consisting of two large-scale toroidal field features with opposite polarities, still led to the generation of large-scale poloidal fields. The authors have verified that the main features of the mechanism described in this manuscript are also present in the simulation of Gottlieb et al. (2023b, private communication). This suggests that the MRI dynamo described in this work is robust to non-axisymmetric initial conditions.

We also found no difference in the magnetic field generation mechanism between a net large-scale toroidal field and an antisymmetric,

zero net flux toroidal field (Appendix B). This independence on the net toroidal flux is encouraging, as it could be that the mechanism can be excited with smaller scale toroidal field structures.

Dhang, Bai & White (2023) found that for an initial weak, $\beta_{\text{ini}} = 800$, zero net flux poloidal field, the accretion disc was unable to reach the MAD state after $120\,000r_g/c$. Although no MAD was present, a coherent large-scale magnetic field appeared, as shown in their Fig. 7 and the last panels of Fig. 2. However, this field did not have enough flux to reach the MAD state. Similarly, Rodman & Reynolds (2024) recently found that initially weak toroidal magnetic field, $\beta_{\text{ini}} = 200$, would not enter the MAD state after $120000r_g/c$. Their results show that weaker magnetic fields might substantially weaken the efficiency of the non-linear magnetic field generation mechanism shown here.

However, in both studies the authors either imposed a null field for $r > 200r_g$ or considered a limited radial domain size $r_{\text{max}} < 300r_g$. As was argued on Section 4.5 the maximal magnetic flux that can be generated is set by the size of the accretion disc, including its radial extent. It would be interesting to consider more weakly magnetized discs like the ones simulated by Rodman & Reynolds (2024) and Dhang et al. (2023), but with a larger radial extent. Following the arguments of Section 4.5, this could lead to the formation of an MAD. Rodman & Reynolds (2024) also ran a simulation with a stronger initial magnetic field, $\beta_{\text{ini}} = 5$, that reached a state more reminiscent of an MAD and might be qualitatively similar to the simulation presented here.

### 6.2.2 Reaching the MAD state

An intriguing implication of our study is that the maximum magnetic flux generated is determined by the radial extent and geometrical thickness of the accretion disc (see e.g. Rodman & Reynolds 2024; Dhang et al. 2023, where different final steady states are computed). This constraint on the maximal magnetic flux that can be generated is specially relevant for mergers, where the radial extent of the disc is limited. At first glance this might be bad news for the MAD state. However, it should be noted that the transition to this MAD state also depends on the accretion rate. For a finite mass reservoir, $\dot{m}$ will drop, and the MAD state will be eventually reached regardless of the maximal magnetic flux that can be generated (Tchekhovskoy 2015; Tchekhovskoy & Giannios 2015; Gottlieb et al. 2023b). However, it is crucial to note that if the system attains the MAD state towards the end of its existence, then the MAD state becomes irrelevant to the system's dynamics.

The study by Gottlieb et al. (2023a) suggests that the X-ray power-law decay observed in the extended emission of compact merger GRBs may be attributed to the system transitioning to the MAD state. If this hypothesis holds, the maximal poloidal magnetic flux that can be generated becomes the primary controlling parameter, along with the initial mass of the disc, for the duration of the GRB event. An exploration of the potential maximal magnetic fluxes based on various parameters of the binary merger, such as the central engine, disc mass and size, could be approximated using our work as a foundation (see Section 4.5). This analysis is deferred to future investigations but has the potential to provide complementary insights into understanding the unified evolution of merger GRBs (Gottlieb et al. 2023a). This analysis could also be useful for models of TDEs (Teboul & Metzger 2023). In such models, the jet efficiency is free parameter, one could use our prescriptions for the generation and transport of magnetic flux to self-consistent compute the jet efficiency.
Finally, the dynamo mechanism outlined in this study, relying heavily on advection, may encounter challenges when the accretion disc becomes geometrically thin. Thinner discs are notably less proficient in the advection of magnetic fields (Lubow et al. 1994; Jacquemin-Ide et al. 2021). It is a well-established fact that thin discs, such as those found during the soft state in X-ray binaries (Done et al. 2007), exhibit reduced tendencies to launch relativistic jets. Consequently, it becomes imperative to investigate the mechanism's dependence on the disc's thickness, as this insight may shed light on why thinner discs tend to be less predisposed to jet launching.

### 6.2.3 Stochasticity and flares

Due to the stochastic nature of the dynamo mechanism, simulations initially set with toroidal magnetic fields exhibit an initial bursty behaviour in jet power and efficiency, preceding the emergence of a dominant polarity (see Fig. 12 and Christie et al. 2019; Gottlieb et al. 2023b). This early activity mirrors the quasi-periodic precursor flares observed in short Kilonova gamma-ray bursts (GRBs; Xiao et al. 2022). Our findings indicate that this behaviour is attributed to the early attachment of a magnetic loop to the BH (see Fig. 4). Although the topology of this attached loop differs from that observed in magnetically arrested discs (MADs), it still facilitates the launching of a Blandford–Znajek (BZ) jet. It is worth noting that the jet produced by this small-scale loop is less efficient compared to those driven by MADs (Christie et al. 2019). The cancellations of magnetic flux on the BH event horizon occur as additional magnetic loops of opposite polarities are advected, like the simulations of Parfrey, Giannios & Beloborodov (2015). These cancellations give rise to oscillations in magnetic flux, which may explain the quasi-periodic precursor flares observed in short GRBs. Such variability could also explain the initially bursty X-ray luminosity curve of the TDE Swift J1644 + 57, that is normally modelled with precession of the accretion disc (Tchekhovskoy et al. 2014; Teboul & Metzger 2023).

As time progresses, a truly large-scale magnetic field is generated and the MAD state is reached. This large-scale magnetic flux may be responsible for the extended emission observed for $t > 1$s in short GRBs (Norris & Bonnell 2006; Gottlieb et al. 2023a). Further investigation of this scenario is warranted, particularly with simulations initialized with weaker magnetic fields.

El Mellah et al. (2022) examined a hybrid magnetic topology, with closed and open field lines connected to the BH, as shown in Fig. 4. They used 2D GR particle-in-cell simulations to find that this hybrid configuration could drive efficient jets and particle acceleration through reconnection sheets. In this work, we showed that such a configuration naturally emerges in the early stages of the evolution of the MRI dynamo. This provides a natural justification for their *previously ad-hoc* configuration.

Overall, the work presented here, anchored in a detailed analysis of 3D GRMHD dynamics, provides new insights into the physical mechanism underlying the global magnetized dynamics and magnetic field generation in black-hole accretion. In particular, simulations of mergers are usually initialized with artificially large magnetic fields due to constraints on resolving the magnetized turbulence. While subgrid models might alleviate some of the issues related to resolution, they rely on assumptions and fine-tuned coefficients. Understanding the complicated nonlinear dynamics therefore paves the way towards better mean-field models that are complementary to the subgrid approaches. Overall, the work presented here provides a contribution to better global large-scale models of BH accretion.

## ACKNOWLEDGEMENTS

JJ thanks Geoffroy Lesur and Jonathan Ferreira for their helpful feedback on the results and Antoine Riols, Nick Kaaz, Omer Bromberg and Beverly Lowell for valuable discussions. We thank the referee for their insightful suggestions.

JJ and AT acknowledge support by the NSFAST-2009884, NASA80NSSC21K1746, and NASAXMM-Newton80NSSC22K0799 grants. AT also acknowledges support by NSF grants AST-2107839, AST-1815304, AST-1911080, OAC-2031997, and AST-2206471. Support for this work was also provided by the National Aeronautics and Space Administration through Chandra Award Number TM1-22005X issued by the Chandra X-ray Center, which is operated by the Smithsonian Astrophysical Observatory for and on behalf of the National Aeronautics Space Administration under contract NAS8-03060. This research was also made possible by NSF PRAC award no. 1615281 at the Blue Waters sustained-petascale computing project and supported in part under grant no. NSFPHY-1125915. This research used resources of the Oak Ridge Leadership Computing Facility, which is a DOE Office of Science User Facility supported under Contract DE-AC05-00OR22725 via ALCC, INCITE, and Director Discretionary allocations PHY129. ML was supported by the John Harvard, ITC, and NASA Hubble Fellowship Program fellowships.

## DATA AVAILABILITY

The data underlying this article will be shared on reasonable request to the corresponding author.

## REFERENCES

Aguilera-Miret R., Viganò D., Carrasco F., Miñano B., Palenzuela C., 2020, Phys. Rev. D, 102, 103006

Aguilera-Miret R., Viganò D., Palenzuela C., 2022, ApJ, 926, L31

Aguilera-Miret R., Palenzuela C., Carrasco F., Viganò D., 2023, Phys. Rev. D, 108, 103001

Bai X.-N., Stone J. M., 2013, ApJ, 767, 30

Balbus S. A., Hawley J. F., 1991, ApJ, 376, 214

Balbus S. A., Hawley J. F., 1992, ApJ, 400, 610

Biscoveanu S., Burns E., Landry P., Vitale S., 2023, RNAAS; Res. Notes AAS, 7, 136

Blandford R. D., Payne D. G., 1982, MNRAS, 199, 883

Blandford R. D., Znajek R. L., 1977, MNRAS, 179, 433

Bloom J. S. et al., 2011, Science, 333, 203

Bonnerot C., Price D. J., Lodato G., Rossi E. M., 2017, MNRAS, 469, 4879

Brandenburg A., 2018, J. Plasma Phys., 84, 735840404

Brandenburg A., Nordlund A., Stein R. F., Torkelsson U., 1995, ApJ, 446, 741

Brandenburg A., Rädler K. H., Rheinhardt M., Käpylä P. J., 2008, ApJ, 676, 740

Burrows D. N. et al., 2011, Nature, 476, 421

Chakrabarti S. K., 1985, ApJ, 288, 1

Chashkina A., Bromberg O., Levinson A., 2021, MNRAS, 508, 1241

Christie I. M., Lalakos A., Tchekhovskoy A., Fernández R., Foucart F., Quataert E., Kasen D., 2019, MNRAS, 490, 4811

Colella P., Woodward P. R., 1984, J. Comput. Phys., 54, 174

Corbel S., Nowak M. A., Fender R. P., Tzioumis A. K., Markoff S., 2003, A&A, 400, 1007

Corbel S., Coriat M., Brocksopp C., Tzioumis A. K., Fender R. P., Tomsick J. A., Buxton M. M., Bailyn C. D., 2013, MNRAS, 428, 2500

Davis S. W., Stone J. M., Pessah M. E., 2010, ApJ, 713, 52
Dhang P., Bendre A., Sharma P., Subramanian K., 2020, MNRAS, 494, 4854

Dhang P., Bai X.-N., White C. J., 2023, ApJ, 944, 182

Dhang P., Bendre A., Subramanian K., 2024, MNRAS, 530, 2778

Done C., Gierliński M., Kubota A., 2007, A&AR, 15, 1

EHT Collaboration T. E. H. T., 2021, ApJ, 910, L13

El Mellah I., Cerutti B., Crinquand B., Parfrey K., 2022, A&A, 663, A169

Fernández R., Quataert E., Schwab J., Kasen D., Rosswog S., 2015, MNRAS, 449, 390

Fernández R., Tchekhovskoy A., Quataert E., Foucart F., Kasen D., 2019, MNRAS, 482, 3373

Ferreira J., Pelletier G., 1995, A&A, 295, 807

Ferreira J., Petrucci P.-O., Henri G., Saugé L., Pelletier G., 2006, A&A, 447, 813

Flock M., Dzyurkevich N., Klahr H., Turner N., Henning T., 2012, ApJ, 744, 144

Foucart F., 2012, Phys. Rev. D, 86, 124007

Foucart F., Duez M. D., Kidder L. E., Teukolsky S. A., 2011, Phys. Rev. D, 83, 024005

Foucart F. et al., 2015, Phys. Rev. D, 91, 124021

Fromang S., Papaloizou J., 2007, A&A, 476, 1113

Gammie C. F., McKinney J. C., Tóth G., 2003, ApJ, 589, 444

Giannios D., Metzger B. D., 2011, MNRAS, 416, 2102

Goldreich P., Lynden-Bell D., 1965, MNRAS, 130, 125

Gompertz B. P., Nicholl M., Smith J. C., Harisankar S., Pratten G., Schmidt P., Smith G. P., 2023, MNRAS, 526, 4585

Gottlieb O. et al., 2023a, ApJ, 958, L33

Gottlieb O. et al., 2023b, ApJ, 954, L21

Gressel O., 2010, MNRAS, 405, 41

Gressel O., Pessah M. E., 2015, ApJ, 810, 59

Guilet J., Ogilvie G. I., 2012, MNRAS, 424, 2097

Hawley J. F., Gammie C. F., Balbus S. A., 1995, ApJ, 440, 742

Hawley J. F., Gammie C. F., Balbus S. A., 1996, ApJ, 464, 690

Hayashi K., Kawaguchi K., Kiuchi K., Kyutoku K., Shibata M., 2021, Phys. Rev. D, 103, 043007

Hayashi K., Fujibayashi S., Kiuchi K., Kyutoku K., Sekiguchi Y., Shibata M., 2022, Phys. Rev. D, 106, 023008

Hayashi K., Kiuchi K., Kyutoku K., Sekiguchi Y., Shibata M., 2023, Phys. Rev. D, 107, 123001

Heinemann T., McWilliams J. C., Schekochihin A. A., 2011, Phys. Rev. Lett., 107, 255004

Held L. E., Mamatsashvili G., 2022, MNRAS, 517, 2309

Held L. E., Mamatsashvili G., Pessah M. E., 2024, MNRAS, 530, 2232

Herault J., Rincon F., Cossu C., Lesur G., Ogilvie G. I., Longaretti P. Y., 2011, Phys. Rev. E, 84, 036321

Hogg J. D., Reynolds C. S., 2018, ApJ, 861, 24

Igumenshchev I. V., Narayan R., Abramowicz M. A., 2003, ApJ, 592, 1042

Jacquemin-Ide J., Lesur G., Ferreira J., 2021, A&A, 647, A192

Janka H. T., Eberl T., Ruffert M., Fryer C. L., 1999, ApJ, 527, L39

Jiménez-Rosales A. et al., 2020, A&A, 643, A56

Johnson B. M., 2007, ApJ, 660, 1375

Kelley L. Z., Tchekhovskoy A., Narayan R., 2014, MNRAS, 445, 3919

Kiuchi K., Cerdá-Durán P., Kyutoku K., Sekiguchi Y., Shibata M., 2015, Phys. Rev. D, 92, 124034

Kiuchi K., Kyutoku K., Sekiguchi Y., Shibata M., 2018, Phys. Rev. D, 97, 124039

Kiuchi K., Reboul-Salze A., Shibata M., Sekiguchi Y., 2024, Nat. Astro., 8, 298

Lesur G., 2021, J. Plasma Phys., 87, 205870101

Lesur G., Ogilvie G. I., 2008a, MNRAS, 391, 1437

Lesur G., Ogilvie G. I., 2008b, A&A, 488, 451

Lesur G., Ferreira J., Ogilvie G. I., 2013, A&A, 550, A61

Liska M., Tchekhovskoy A., Quataert E., 2020, MNRAS, 494, 3656

Liska M. T. P. et al., 2022, ApJS, 263, 26

Lubow S. H., Papaloizou J. C. B., Pringle J. E., 1994, MNRAS, 267, 235

Lynden-Bell D., 2003, MNRAS, 341, 1360

Lynden-Bell D., 2006, MNRAS, 369, 1167

Lynden-Bell D., Pringle J. E., 1974, MNRAS, 168, 603

Mamatsashvili G., Chagelishvili G., Pessah M. E., Stefani F., Bodo G., 2020, ApJ, 904, 47

Manikantan V., Kaaz N., Jacquemin-Ide J., Musoke G., Chatterjee K., Liska M., Tchekhovskoy A., 2024, ApJ, 965, 175

Marcel G. et al., 2019, A&A, 626, A115

Margutti R., Chornock R., 2021, ARA&A, 59, 155

Metzger B. D., 2019, Living Rev. Rel., 23, 1

Metzger B. D., Arcones A., Quataert E., Martínez-Pinedo G., 2010, MNRAS, 402, 2771

Mochkovitch R., Hernanz M., Isern J., Martin X., 1993, Nature, 361, 236

Moffatt H. K., 1978, Magnetic Field Generation in Electrically Conducting Fluids. Cambridge Univ. Press, Cambridge

Mondal T., Bhat P., 2023, Phys. Rev. E, 108, 065201

Most E. R., Quataert E., 2023, ApJ, 947, L15

Most E. R., Papenfort L. J., Tootle S. D., Rezzolla L., 2021, MNRAS, 506, 3511

Nakar E., 2020, Phys. Rep., 886, 1

Narayan R., Igumenshchev I. V., Abramowicz M. A., 2003, PASJ, 55, L69

Norris J. P., Bonnell J. T., 2006, ApJ, 643, 266

Paczynski B., 1991, AcA, 41, 257

Parfrey K., Giannios D., Beloborodov A. M., 2015, MNRAS, 446, L61

Price D. J., Rosswog S., 2006, Science, 312, 719

Rezzolla L., Giacomazzo B., Baiotti L., Granot J., Kouveliotou C., Aloy M. A., 2011, ApJ, 732, L6

Rincon F., 2019, J. Plasma Phys., 85, 167

Rincon F., Ogilvie G. I., Proctor M. R. E., 2007, Phys. Rev. Lett., 98, 254502

Riols A., Rincon F., Cossu C., Lesur G., Longaretti P. Y., Ogilvie G. I., Herault J., 2013, J. Fluid Mech., 731, 1

Riols A., Rincon F., Cossu C., Lesur G., Ogilvie G. I., Longaretti P.-Y., 2015, A&A, 575, A14

Riols A., Rincon F., Cossu C., Lesur G., Ogilvie G. I., Longaretti P.-Y., 2017, A&A, 598, A87

Rodman P. E., Reynolds C. S., 2024, ApJ, 960, 97

Rosswog S., 2005, ApJ, 634, 1202

Ruiz M., Lang R. N., Paschalidis V., Shapiro S. L., 2016, ApJ, 824, L6

Ruiz M., Shapiro S. L., Tsokaros A., 2018, Phys. Rev. D, 98, 123017

Salvesen G., Simon J. B., Armitage P. J., Begelman M. C., 2016, MNRAS, 457, 857

Shakura N. I., Sunyaev R. A., 1973, A&A, 24, 337

Shi J.-M., Stone J. M., Huang C. X., 2016, MNRAS, 456, 2273

Shibata M., Uryū K., 2006, Phys. Rev. D, 74, 121503

Shibata M., Uryū K., 2007, Class. Quantum Gravity, 24, S125

Siegel D. M., Metzger B. D., 2017, Phys. Rev. Lett., 119, 231102

Squire J., Bhattacharjee A., 2015, Phys. Rev. Lett., 115, 175003

Squire J., Bhattacharjee A., 2016, J. Plasma Phys., 82, 535820201

Tchekhovskoy A., 2015, in Contopoulos I., Gabuzda D., Kylafis N., eds, Astrophysics and Space Science Library Vol. 414, The Formation and Disruption of Black Hole Jets. Springer, Cham., Switzerland, p. 45

Tchekhovskoy A., Giannios D., 2015, MNRAS, 447, 327

Tchekhovskoy A., Narayan R., McKinney J. C., 2011, MNRAS, 418, L79

Tchekhovskoy A., Metzger B. D., Giannios D., Kelley L. Z., 2014, MNRAS, 437, 2744

Teboul O., Metzger B. D., 2023, ApJ, 957, L9

Xiao S. et al., 2022, preprint (arXiv:2205.02186)

Zanni C., Ferreira J., 2013, A&A, 550, A99

## APPENDIX A: ENERGY GENERATION OR DISSIPATION OF THE FIRST EIGHT NON-AXISYMMETRIC MODES

In this appendix, we confirm that the first three non-zero *m*-modes, referred to as active modes, are the sole contributors to the generation of poloidal magnetic energy. Furthermore, we show that these modes dissipate toroidal magnetic energy and facilitate angular momentum
<figure><img src="image_19.png" /><figcaption>**Figure A1.** Vertically and temporally averaged poloidal (a) (equation 21), and toroidal (equation 17), (b) energy equation normalized to the local Keplerian frequency. The temporal average is performed between 5000$r_g/c$ and 8000$r_g/c$. The different terms are divided by the vertical and temporal average of $\langle B_p \rangle^2$ (a) or $\langle B_\varphi \rangle$ (b). (a) The coloured solid lines show the first eight non-axisymmetric $m$-modes, $\delta A_p^{(m)}$ (a) or $\delta A_\varphi^{(m)}$ (b). The dotted line shows the large-scale axisymmetric transport of poloidal (a) or toroidal (b) field. The dashed line shows the small-scale turbulent transport without the first eight non-axisymmetric modes, $\delta A_p^{small}$ (a) or $\delta A_\varphi^{small}$. The dotted-dashed line shows the shear of poloidal into toroidal field, $S_\varphi$.</figcaption></figure>

transport. In essence, the active modes are identified as MRI-driven shearing waves, mediated by the axisymmetric toroidal field.

## A1 Poloidal magnetic energy

In Fig. A1(a), we present the various source terms for axisymmetric poloidal magnetic energy (equation 21). These terms are vertically averaged within the disc, temporally averaged between 5000 $r_g/c$ and 8000 $r_g/c$, and normalized by the averaged axisymmetric poloidal magnetic energy, $\frac{1}{2}\langle B_p \rangle^2$. The sources of axisymmetric poloidal magnetic energy are categorized into three components: (1) The contributions of the first eight non-linear correlations of $m$-modes, denoted as $\delta A_p^{(m)}$, as defined in equations (25) and (26). (2) The large-scale axisymmetric transport of $\varphi$-averaged poloidal field, represented by $A_p$. (3) The small-scale turbulent transport, excluding

<figure><img src="image_20.png" /><figcaption>**Figure A2.** Vertically and temporally averaged of the different component of the magnetic stresses computed from equations (14, A7, A8). The temporal average is performed between 5000$r_g/c$ and 8000$r_g/c$ and the division between pressure and stress tensor is performed after both averages. The coloured solid lines show the first eight non-axisymmetric $m$-modes, $\alpha_n u^{(m)}$. The dotted line shows the large-scale axisymmetric transport of angular momentum by $-\langle B_r \rangle \langle B_\varphi \rangle$. The dashed line shows the small-scale turbulent transport of all non-axisymmetric $m$-modes, $\alpha_\varphi = \sum_{m=1}^\infty \alpha_\varphi^{(m)}$.</figcaption></figure>

the first eight non-linear correlations of $m$-modes,

$$ \delta A_p^{\text{small}} = \delta A_p - \sum_{m=1}^{8} \delta A_p^{(m)}. \qquad (A1) $$

Fig. A1(a) confirms that the only $m$-mode correlations that generate axisymmetric poloidal magnetic energy on average are $m = 1, 2, 3$, where $m = 1, 2$ are the main contributors. We notice that $\delta A_p^{small} \simeq \sum_{m=3}^\infty \delta A_p^{(m)}$, as this term is dominated by the more numerous large $m$-mode correlations (see Fig. 10). Consequently, we refrain from providing a detailed analysis of the magnitude and behaviour of $A_p$ and $\delta A_p^{small}$, as this information is discussed in Section 4.

## A2 Toroidal magnetic energy

We analyse the sources of axisymmetric toroidal magnetic energy by decomposing them into their non-linear correlations of non-zero-$m$-modes. Extending the same decomposition applied to the poloidal magnetic field, we obtain

$$ \delta A_\varphi = \sum_{m=1}^{\infty} \delta A_\varphi^{(m)}, \qquad (A2) $$

where

$$ \delta A_\varphi^{(m)} = \langle \mathbf{B}_\varphi \rangle \mathbf{e}_\varphi \cdot \nabla \times (\mathcal{E}_p^{(m)}), \qquad (A3) $$

and

$$ \mathcal{E}_p^{(m)} = \mathcal{E}^{(m)} - \mathcal{E}_\varphi^{(m)} \qquad (A4) $$

with $\mathcal{E}_\varphi^{(m)}$ defined in equation (26), or more explicitly

$$ \mathcal{E}_p^{(m)} = \mathcal{R}[u^m \times B^m]^* - \mathcal{R}[u_p^m \times B_p^m]^*. \qquad (A5) $$
<figure><img src="image_21.png" /><figcaption>**Figure A3.** Snapshots of the axisymmetrized toroidal magnetic field, $R \langle B_{\varphi} \rangle$ in colour as a function of R and z. We also show the poloidal magnetic field lines through the poloidal magnetic flux (equation 9, dashed lines show negative polarity and solid lines show positive polarity). Magnetic field structures become larger and larger with time. Notice that in this case the initial condition, visible at large radii in panel (a) is antisymmetric with respect to the disc mid-plane.</figcaption></figure>

We also define the small scale, large m, non-linear correlations as

$$ \delta A_{\varphi}^{\text{small}} = \delta A_{\varphi} - \sum_{m=1}^{8} \delta A_{\varphi}^{(m)}. \qquad (A6) $$

In Fig. A1(b), we present the same decomposition as shown in Fig. A1(a), but now for the source terms of toroidal magnetic energy. These terms are calculated using equations (17–A2-A4), and all of them are subject to vertical averaging within the disc. Temporal averaging is performed between 5000, $r_g/c$ and 8000, $r_g/c$, and the values are normalized by the averaged axisymmetric toroidal magnetic energy, $\frac{1}{2} \langle B_{\varphi} \rangle^2$.

In line with Section 4, we demonstrate that the sole term responsible for generating axisymmetric toroidal magnetic energy is $S\varphi$, associated with the large-scale shear of the axisymmetric poloidal field. The large-scale advection term, $A\varphi$, results in a loss of magnetic energy throughout the entire disc, corresponding to the vertical transport described in Sections 4 and 3. All large-scale non-linear m-mode correlations dissipate axisymmetric toroidal magnetic energy, in contrast to the decomposition of axisymmetric poloidal magnetic energy, here even the large-scale active modes ($m = 1, 2, 3$) consume toroidal magnetic energy. This observation aligns with findings from shearing box simulations (Lesur & Ogilvie 2008b; Riols et al. 2015, 2017). Once again, we observe that dissipation of the large scales is

primarily driven by the more numerous small-scale m-modes, leading to the approximation that $\delta A\varphi^{\text{small}} \approx \delta A\varphi$.

### A3 Angular momentum transport

To show that the active modes transport angular momentum, we decompose the Maxwell stress using Parseval’s theorem

$$ \alpha_{\nu} = \sum_{m=1}^{\infty} \alpha_{\nu}^{(m)}, \qquad (A7) $$

where

$$ \alpha_{\nu}^{(m)} = -\frac{1}{4\pi \langle P \rangle} \mathcal{R} [B_{r}^{m} \times B_{\varphi}^{m}]^{*}. \qquad (A8) $$

Then, we define the large-scale laminar stress (Jacquemin-Ide et al. 2021; Manikantan et al. 2024)

$$ \alpha_{\nu}^{\text{la}} = -\frac{1}{4\pi \langle P \rangle} \langle B_r \rangle \langle B_{\varphi} \rangle. \qquad (A9) $$

In Fig. A2, we present the m-mode decomposition of the magnetic angular momentum stresses using equations (17, A2, A4). All terms are subjected to vertical averaging within the disc and temporal averaging between 5000 $r_g/c$ and 8000 $r_g/c$; the division between pressure and stress tensor is performed after both averages. The active
<figure><img src="image_22.png" /><figcaption>**Figure B1.** Vertically and temporally averaged poloidal (a) (equation 21), and toroidal (equation 17) (b) energy equation normalized to the local Keplerian frequency. The temporal average is performed between 4000$r_g/c$ and $7000r_g/c$. The different terms are divided by the vertical and temporal average of $\langle B_p \rangle^2$ (a) or $\langle B_\psi \rangle$ (b). (a) The coloured solid lines show the first eight non-axisymmetric $m$-modes, $\delta A_p^{(m)}$ (a) or $\delta A_\psi^{(m)}$ (b). The dotted line shows the large-scale axisymmetric transport of poloidal (a) or toroidal (b) field. The dashed line shows the small-scale turbulent transport without the first eight non-axisymmetric modes, $\delta A_p^{\text{small}}$ (a) or $\delta A_\psi^{\text{small}}$. The dotted-dashed line shows the shear of poloidal into toroidal field, $S_\psi$.</figcaption></figure>

modes ($m = 1, 2, 3$) contribute to angular momentum transport, but they are not the dominant stress term. The turbulent stress is primarily driven by the more numerous small-scale modes, this observation aligns with the contribution of the latter to turbulent dissipation.

## APPENDIX B: SIMULATION WITH AN ANTISYMMETRIC TOROIDAL FIELD

In this appendix, we describe the simulation initialized with an antisymmetric, with respect to the disc mid-plane, toroidal magnetic field. This simulation shows identical features to the one described in the main body of the manuscript.

Fig. A3 shows the axisymmetrized poloidal magnetic field lines and the axisymmetric toroidal magnetic field. We observe very similar evolution to that presented in Fig. 1, even though the initial condition is very different: (1) large-scale magnetic fields are generated and vertically shed (see Fig. A3a). (2) With time the size of the largest magnetic field structure increases, we notice advection of the larger field from the outer radii into the inner regions (Fig. A3b,c).

Fig. A3 illustrates the axisymmetrized poloidal magnetic field lines alongside the axisymmetric toroidal magnetic field. Despite the contrasting initial condition, we observe a very similar evolution to that depicted in Fig. 1: (1) The generation of large-scale magnetic fields, which are vertically expelled through ballooning instability (Lynden-Bell 2003, see also Fig. A3a). (2) Over time, an increase in the size of the magnetic field structures, with noticeable advection of the larger field from outer radii into the inner regions (see Fig. A3 b-c).

Fig. B1(a) shows the source terms for axisymmetric poloidal magnetic energy for the simulation initiated with an antisymmetric toroidal field, mirroring the content of Fig. A1(a) or the symmetric case. The trends observed in the generation of poloidal magnetic field source terms align with the descriptions in Section 4 and Appendix A.

Specifically, the findings reveal that in this case too only the largest non-linear correlations of $m$-modes ($m = 1, 2, 3$) generate axisymmetric poloidal magnetic fields. Large-scale advection, $A_p$, facilitates the transport of flux from outer regions, while the smaller non-linear correlations of $m$ modes ($m > 3$) contribute to the dissipation of the poloidal magnetic field. Moreover, the amplitudes of these terms exhibit striking similarity with the symmetric case, yielding an identical regeneration timescale of approximately $t_{\text{gen}} =$

$$ \frac{1}{2} \langle B_p \rangle^2 / \sum_{m=1}^{3} A_p^{(m)} \approx 10 \Omega_K^{-1}. $$

The patterns and amplitudes are also consistent for axisymmetric toroidal magnetic energy generation, as shown in Fig. B1(b). Similarly to the poloidal magnetic energy case, the observations remain unchanged: all non-linear correlations of $m$ modes contribute to magnetic energy dissipation, magnetic energy generation primarily results from the shear of the poloidal field, and the toroidal magnetic field is shed vertically through the advective term.

We conclude that the dynamo mechanism is independent of the global symmetry of initial conditions.

## APPENDIX C: TIME AND LAG DIAGRAM FOR CORRELATIONS
<figure><img src="image_23.png" /></figure>

**Figure C1.** Vertically averaged lag-time diagram of the Pearson correlations at $r = 20r_g$. The averages are computed between $t = 20\Omega_K^{-1}(r = 20)$ and $t = 100\Omega_K^{-1}(r = 20)$, and between $\theta_1 = \pi/2 - \arctan(h/r)$ and $\theta_2 = \pi/2 - \arctan(0.1 h/r)$. The correlations are shown for (a) the $\Omega$-effect which is used as a baseline and both models, (b) equations (47), and (c) equation (48). Notice that on short time-scales both models (b,c) show strong correlations, but the shear current effect (c) shows a more consistent correlation for $\tau = 0.5\Omega_K^{-1}$.

This paper has been typeset from a TEX/LATEX file prepared by the author.