# Papaloizou–Pringle instability suppression by the magnetorotational instability in relativistic accretion discs

Matteo Bugli, Jérôme Guilet, Ewald Müller, Luca del Zanna, Niccolò

Bucciantini, Pedro J. Montero

► To cite this version:

Matteo Bugli, Jérôme Guilet, Ewald Müller, Luca del Zanna, Niccolò Bucciantini, et al.. Papaloizou–Pringle instability suppression by the magnetorotational instability in relativistic accretion discs. Monthly Notices of the Royal Astronomical Society, 2018, 475 (1), pp.108-120. 10.1093/mn-ras/stx3158. hal-01704910

HAL Id: hal-01704910

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

Submitted on 11 Nov 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_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Papaloizou–Pringle instability suppression by the magnetorotational instability in relativistic accretion discs

M. Bugli,<sup>1,2*</sup> J. Guilet,<sup>1,3,4</sup> E. Müller,<sup>1,2</sup> L. Del Zanna,<sup>5,6,7</sup> N. Bucciantini<sup>5,6,7</sup>
and P. J. Montero<sup>1</sup>

Accepted 2017 December 4. Received 2017 December 3; in original form 2017 July 6

## ABSTRACT

Geometrically thick tori with constant specific angular momentum have been widely used in the last decades to construct numerical models of accretion flows on to black holes. Such discs are prone to a global non-axisymmetric hydrodynamic instability, known as Papaloizou–Pringle instability (PPI), which can redistribute angular momentum and also lead to an emission of gravitational waves. It is, however, not clear yet how the development of the PPI is affected by the presence of a magnetic field and by the concurrent development of the magnetorotational instability (MRI). We present a numerical analysis using three-dimensional GRMHD simulations of the interplay between the PPI and the MRI considering, for the first time, an analytical magnetized equilibrium solution as initial condition. In the purely hydrodynamic case, the PPI selects as expected the large-scale $m = 1$ azimuthal mode as the fastest growing and non-linearly dominant mode. However, when the torus is threaded by a weak toroidal magnetic field, the development of the MRI leads to the suppression of large-scale modes and redistributes power across smaller scales. If the system starts with a significantly excited $m = 1$ mode, the PPI can be dominant in a transient phase, before being ultimately quenched by the MRI. Such dynamics may well be important in compact star mergers and tidal disruption events.

**Key words:** accretion, accretion discs – instabilities – MHD – plasmas – turbulence – waves.

## 1 INTRODUCTION

Accretion on to black holes provides one of the most efficient mechanisms to power up high-energy astrophysical sources such as active galactic nuclei (Rees 1984; Marconi et al. 2004; Reynolds 2014), X-ray binaries (Narayan & Yi 1995; Fender, Belloni & Gallo 2004; Remillard & McClintock 2006), and gamma-ray bursts (Woosley 1993; Piran 1999; Kumar & Zhang 2015), just to cite a few. The conservation of angular momentum commonly leads to the formation of accretion discs, which in order to convert gravitational binding energy into thermal, kinetic or magnetic energy need to lose their angular momentum sufficiently fast.

Since black hole-disc models were proposed as central engine for quasars by Lynden-Bell (1969), there has been a continuous interest

in the physics underlying accretion on to black holes (for a general review, see Abramowicz & Fragile 2013). The seminal papers by Shakura & Sunyaev (1973) and Lynden-Bell & Pringle (1974) first described what now is known as the *standard disc model*: a geometrically thin, optically thick Keplerian disc where the accretion process is driven by a local turbulent viscosity that is parametrized by the quantity $\alpha$. Given the great success of the standard disc model in providing a self-consistent way to enable accretion along with accurate predictions for the observed emission, the actual nature and physical mechanism behind the parameter $\alpha$ (which essentially gives an estimate of the efficiency of the angular momentum transport in the disc) has been the object of numerous studies and is currently still under investigation. The capability of hydrodynamic Keplerian disc models to explain from first principles the onset of a turbulent accretion flow has been debated for decades. They are indeed stable to local linear perturbations, since their distribution of specific angular momentum increases with distance from the

<footnote n="1"><sup>1</sup>Max-Planck-Institut für Astrophysik, Karl-Schwarzschild-Str. 1, D-85741 Garching, Germany</footnote>

<footnote n="2"><sup>2</sup>Technische Universität München, Physik-Department, James-Franck-Str. 1, D-85748 Garching, Germany</footnote>

<footnote n="3"><sup>3</sup>Max Planck/Princeton Center for Plasma Physics, Max-Planck-Institut für Astrophysik, Karl-Schwarzschild-Str. 1, D-85741 Garching, Germany</footnote>

<footnote n="4"><sup>4</sup>Laboratoire AIM, CEA/DRF-CNRS-Université Paris Diderot, IRFU/Département d'Astrophysique, CEA-Saclay, F-91191, France</footnote>

<footnote n="5"><sup>5</sup>Dipartimento di Fisica e Astronomia, Università degli Studi di Firenze, Via G. Sansone 1, I-50019 Sesto Fiorentino (Firenze), Italy</footnote>

<footnote n="6"><sup>6</sup>INAF – Osservatorio Astrofisico di Arcetri, Largo E. Fermi 5, I-50125 Firenze, Italy</footnote>

<footnote n="7"><sup>7</sup>INFN – Sezione di Firenze, Via G. Sansone 1, I-50019 Sesto Fiorentino (Firenze), Italy</footnote>

<footnote>* E-mail: matteo@mpa-garching.mpg.de</footnote>
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Schematic representation of the mechanism responsible for the PPI unstable modes. The negative-energy waves from the inner part of the disc (where $\omega < \Omega$) interact with the positive-energy ones from the outer disc (where $\omega > \Omega$), redistributing angular momentum across the corotation radius with a positive feedback by the reflective boundaries. The grey area represents the *forbidden zone* surrounding the corotation radius (dashed circle), where the waves become evanescent.</figcaption></figure>

central object (this is the well-known Rayleigh stability criterion), and in general, a small displacement of a fluid element will lead to epicyclic oscillations (Kato 2016). The non-linear stability of hydrodynamic discs is to this date a matter of debate. The first series of shearing box simulations tackling this problem (Balbus, Hawley & Stone 1996; Hawley, Balbus & Winters 1999) suggested that hydrodynamic accretion flows do not present a transition to a turbulent state, but it has been argued by Longaretti (2002) that their results were due to a lack of resolution. More recently Lesur & Longaretti (2005) showed that at high enough Reynolds numbers non-linear perturbations could lead to self-sustained turbulence, which, however, would be too weak to explain observed accretion rates.

An important class of hydrodynamic discs often used in the context of accretion on to black holes is constituted by accretion tori (Abramowicz, Jaroszynski & Sikora 1978), also referred to as *Polish doughnuts*. These thick discs have a large internal energy, and they rely on pressure gradients to support the disc together with centrifugal forces, resulting in a significant vertical thickening of the disc and a departure from a Keplerian distribution of specific angular momentum. Despite their local stability, Papaloizou & Pringle (1984) discovered that they are prone to develop a global non-axisymmetric instability [known as Papaloizou-Pringle instability (PPI); see Fig. 1], which is able to transport angular momentum outwards. Although capable of triggering some accretion, the PPI cannot explain in a satisfactory way the ubiquity of accreting systems, since it mainly affects nearly constant angular momentum tori (Goldreich, Goodman & Narayan 1986; Blaes & Hawley 1988). The PPI does also not fit well in the standard disc model (which assumes locally generated turbulence to enable accretion), as the instability tends to saturate in strong spiral pressure waves in radially wide, nearly

constant angular momentum tori (Hawley 1991; De Villiers & Hawley 2002).

The breakthrough in accretion theory was the realization that magnetic fields are the key to explain how discs can get rid of their angular momentum. The discovery of the *magnetorotational instability* (MRI) in astrophysics by Balbus & Hawley (1991) provided a local mechanism, efficient for a wide range of magnetic field strength, which leads to a growth of linear perturbations on dynamical time-scales and naturally develops MHD turbulence. Since then, the properties of the MRI have been studied in great detail, from both a local and a global point of view (Balbus & Hawley 1998; Fromang 2013; Blaes 2014).

Despite the fundamental importance of magnetic fields in providing a general and universal mechanism to enable accretion in astrophysical discs, the PPI is still quite relevant as an agent of global non-axisymmetric instability, since thick discs with sub-Keplerian angular momentum distributions are expected to form in binary neutron stars (Rezzolla et al. 2010; Kiuchi et al. 2010) or black hole-neutron star (Shibata & Uryū 2006; Foucart et al. 2012) mergers, after the rotational gravitational collapse of massive stars (MacFadyen & Woosley 1999; Aloy et al. 2000) and in tidal disruption events (Loeb & Ulmer 1997; Coughlin & Begelman 2014). The stability of such wide tori has been studied from both an analytical (Goldreich et al. 1986; Glatzel 1987) and numerical (Blaes & Hawley 1988; Hawley 1991; De Villiers & Hawley 2002) point of view, and they have been proven to be quite generally unstable to some non-axisymmetric mode induced by the PPI. Moreover, a residual kick velocity of the central black hole after the merger can excite large-scale spiral shocks in the accretion torus, even for an almost Keplerian distribution of angular momentum (Zanotti et al. 2010). In recent years, there have also been several studies that included self-gravity of the disc (Kiuchi et al. 2011; Korobkin et al. 2011; Mewes et al. 2016). They have shown how the non-axisymmetric structures that arise from the instability can lead to a significant emission of gravitational waves.

Hawley (2000) systematically investigated for the first time the evolution of three-dimensional (3D) magnetized tori, presenting a series of models with different gravitation potential (Newtonian and pseudo-Newtonian), distribution of angular momentum, magnetic field strength and topology (toroidal and poloidal), and azimuthal range (from one quadrant to the full $2\pi$ angle). Despite considering models unstable to the PPI, almost each model failed to display a significant growth of the PPI, as the relatively faster development of the MRI led to a suppression of the hydrodynamic instability. The only model that displayed a non-negligible (although still not dominant) growth of the PPI was the one threaded by a subthermal constant toroidal field (named CT2 in the paper). This set-up triggers, in fact, a non-axisymmetric MRI mode whose growth is slower than the one selected by a vertical magnetic field, hence enabling an early linear development of the PPI.

Fu & Lai (2011) tried to analytically establish the effect of magnetic fields on the development of the PPI. Their analysis (which assumes an incompressible fluid) suggests instead that sufficiently strong magnetic fields can actually further destabilize the torus and reinvigorate the hydrodynamic instability. Despite this, in the last two decades several studies have investigated the dynamics of thick magnetized tori accreting on to black holes with 3D global simulations (Arlt & Rüdiger 2001; Hawley 2001; Hawley & Krolik 2001; De Villiers & Hawley 2003; De Villiers, Hawley & Krolik 2003; Machida & Matsumoto 2003; Machida, Nakamura & Matsumoto 2004; De Villiers et al. 2005; Kigure & Shibata 2005; Fragile et al. 2007; McKinney & Blandford 2009;
McKinney, Tchekhovskoy & Blandford 2012; Wielgus et al. 2015; Fragile & Sądowski 2017) and none of them reported any significant growth of the PPI once the presence of weak magnetic fields of different topologies and strengths were taken into account.

These results suggest that the MRI can efficiently inhibit the onset of the PPI in accretion tori. However, it is yet not completely clear what is the extent of their validity and whether there are serious limitations to them in astrophysically relevant scenarios. Indeed, the vast majority of the previously cited numerical studies considered equilibrium hydrodynamic solutions threaded by a poloidal magnetic field, hence setting an environment favourable for the development of axisymmetric MRI modes (also referred to as channel modes), which generally have a faster growth rate than the most unstable PPI mode (Hawley 2000). Only few of previous works considered models with purely toroidal magnetic fields and the full azimuthal range (Hawley 2000; Hawley & Krolik 2001; Machida & Matsumoto 2003; Machida et al. 2004; McKinney et al. 2012). However, such models are necessary conditions to study the interaction between the PPI in wide tori (where the fastest growing mode has an azimuthal number $m = 1$) and the non-axisymmetric MRI, whose modes have much slower growth rates than their axisymmetric counterparts and could in principle allow for a significant (although transient) growth of the PPI (as suggested by model CT2 in Hawley 2000).

A few recent numerical works employed as initial condition the analytical equilibrium solution provided by Komissarov (2006), where the extra-support of a purely toroidal magnetic field is consistently taken into account in the initial state of the disc. For instance, Wielgus et al. (2015) present a detailed study of the onset of non-axisymmetric MRI modes in relativistic accretion tori, but the models considered only covered a limited azimuthal range ($\phi \in \{0, \pi/2\}$). Therefore, these authors only consider the dynamics of modes with an azimuthal number $m$ being a multiple of 4, filtering out a possibly dominant $m = 1$ mode that would be selected for a hydrodynamical wide torus. The same limited azimuthal range was employed in Fragile & Sądowski (2017), who also considered stronger magnetic fields. Moreover, in most studies, the initial solution is perturbed with a random fluctuation in order to trigger the growth of instabilities. This choice naturally favours the early dominance of the fastest growing instability, hence preventing a more general study of the dynamics that characterizes the simultaneous growth of such two fast instabilities as the PPI and MRI.

In this work, we present 3D GRMHD simulations of magnetized tori with the goal to establish the effect of the MRI on the PPI development starting from the equilibrium solution by Komissarov (2006), which consistently includes the support of a purely toroidal magnetic field in the disc’s initial condition. The goal of this paper is to assess whether the PPI can significantly grow (and if so, to what extent) in a more favourable environment, where a purely toroidal field triggers the growth of slow non-axisymmetric MRI modes and where both random and monochromatic initial perturbations are considered. We first analyse the behaviour of an unmagnetized torus to characterize the standard (and well known) development of the hydrodynamic instability. We then consider equilibria with different strengths of the toroidal field, perturbation seeds, and grid resolutions to assess how the growth of non-axisymmetric global modes depends on these parameters.

The plan of the paper is as follows. The disc models and the numerical set-up are provided in Section 2, while in Section 3 we introduce the diagnostics used to analyse the data produced by the simulations. We present then our results and discuss them in Section 4, and we finally give our conclusions in Section 5.

<table><caption>**Table 1.** List of the parameters defining the initial unperturbed state of the hydrodynamic torus. The same parameters were used to initialize the magnetized models (see Table 2 for the additional ones related to the magnetic field).</caption><thead><tr><th><i>M</i><sub>BH</sub></th><th><i>a</i></th><th><i>r</i><sub>in</sub></th><th><i>r</i><sub>c</sub></th><th><i>l</i></th><th><i>ρ</i><sub>c</sub></th></tr></thead><tbody><tr><td>1</td><td>0</td><td>6.16</td><td>10.17</td><td>3.97</td><td>1</td></tr></tbody></table>

## 2 DISC MODEL AND NUMERICAL SET-UP

For all our simulations we used the code ECHO, which integrates the full set of GRMHD equations and has been used in the past to study non-Keplerian discs around black holes (Zanotti et al. 2010; Bugli, Del Zanna & Bucciantini 2014). A complete description of the code’s algorithms and structure is provided in Del Zanna et al. (2007) and Bucciantini & Del Zanna (2011).

Although ECHO can use a covariant closure for Ohm’s law that includes magnetic resistivity and mean-field dynamo action [see Bucciantini & Del Zanna (2013) and Del Zanna et al. (2014) for more details on the modelling, Bugli et al. (2014) for a study of $\alpha - \Omega$ dynamos in accretion discs, and Del Zanna et al. (2016) for an application to relativistic reconnection], here we assume the plasma to be a perfect conductor by using the condition of vanishing electric field in the fluid comoving frame:

$$E = -v \times B, \quad (1)$$

where $E, B$, and $v$ are, respectively, the electric field, magnetic field, and fluid velocity measured in the frame of the Eulerian observer (also referred to as the Zero Angular Momentum Observer, ZAMO). The effects of magnetic dissipation will be included in a forthcoming work.

Since we are interested in the dynamics taking place within the torus and not in proximity of the black hole event horizon, we consider for the sake of simplicity a non-rotating black hole of mass $M = 1$ in a spherical coordinate system $(r, \theta, \phi)$. The black hole is surrounded by a thick torus (Abramowicz et al. 1978) whose inner edge and centre are located at $r_{in} = 6.16$ (close to the last marginally stable orbit located at $r_{ms} = 6$) and $r_c = 10.17$, respectively. This choice sets the specific angular momentum of the disc to $l = 3.97$. The orbital period at the disc centre is $P_c \sim 207$. We adopt the Cowling approximation, i.e. assuming a time-independent metric, and neglect the self-gravity of the disc, disregarding also any change in the central black hole mass and spin due to accretion. The rest mass density $\rho$ can therefore be rescaled to the central peak value at $r_c$, that is $\rho_c$ (in the following density is always to be intended in this way). We assume an ideal EoS for a perfect gas with adiabatic index set everywhere to $\Upsilon = 4/3$.

For the magnetized tori, we initialize our simulations with the stationary solution provided by Komissarov (2006), keeping the same parameters as in the hydrodynamic models for black hole spin, disc’s inner edge and centre location, and density normalization (see Table 1). Indicating the gas thermal pressure with $p$, we vary the value of the magnetization $\sigma = B^2/2p$ at the disc’s centre, i.e. $\sigma_c$, to investigate the role of the magnetic field strength on the system’s stability. Further details on the initialization of the magnetized thick disc can be found in Del Zanna et al. (2007).

The atmosphere is initialized as a Michel’s radial inflow (Michel 1972), a stationary solution in the Schwarzschild metric determined by the adiabatic index $\Upsilon = 4/3$ and the value of the atmospheric density at distance $r_c$, i.e. $\rho_{atm} = 10^{-6}$. To provide stability for the integration, we set a numerical floor value for the density equal to $\rho_{fl} = 10^{-9}$.
<table><caption>**Table 2.** List of the models considered in our study. $N_r$, $N_\theta$, and $N_\phi$ are the number of grid points in radial, polar, and azimuthal direction, and $A(v^\phi)$ is the amplitude of the initial perturbations. $\sigma_c$ and $|B|_c$ are, respectively, the initial value of the magnetization and magnetic field magnitude at the centre of the disc. $\langle Q_\phi \rangle_V$ is the initial volume average of the quality metric $Q_\phi$ (weighted by the rest mass density $\rho$).</caption><thead><tr><th></th><th>N<sub>r</sub></th><th>N<sub>θ</sub></th><th>N<sub>φ</sub></th><th>A(v<sup>φ</sup>)</th><th>Excitation</th><th>σ<sub>c</sub></th><th>|B|<sub>c</sub></th><th>⟨Q<sub>φ</sub>⟩<sub>V</sub></th></tr></thead><tbody><tr><td>H64r3</td><td>256</td><td>256</td><td>64</td><td>10<sup>-3</sup></td><td>Random</td><td>0</td><td>0</td><td>0</td></tr><tr><td>H256r-4</td><td>256</td><td>256</td><td>256</td><td>10<sup>-4</sup></td><td>Random</td><td>0</td><td>0</td><td>0</td></tr><tr><td>H{32...256} m1-4</td><td>256</td><td>256</td><td>32, ..., 256</td><td>10<sup>-4</sup></td><td>m = 1</td><td>0</td><td>0</td><td>0</td></tr><tr><td>H64m{2..5} -2</td><td>256</td><td>256</td><td>64</td><td>10<sup>-2</sup></td><td>m = 2, ..., 5</td><td>0</td><td>0</td><td>0</td></tr><tr><td>H256r</td><td>256</td><td>256</td><td>256</td><td>10<sup>-3</sup></td><td>Random</td><td>0</td><td>0</td><td>0</td></tr><tr><td>M1256r</td><td>256</td><td>256</td><td>256</td><td>10<sup>-3</sup></td><td>Random</td><td>10<sup>-2</sup></td><td>8.06 × 10<sup>-3</sup></td><td>6.94</td></tr><tr><td>M1512r</td><td>256</td><td>256</td><td>512</td><td>10<sup>-3</sup></td><td>Random</td><td>10<sup>-2</sup></td><td>8.06 × 10<sup>-3</sup></td><td>13.89</td></tr><tr><td>Mm256r</td><td>256</td><td>256</td><td>256</td><td>10<sup>-3</sup></td><td>Random</td><td>3 × 10<sup>-2</sup></td><td>1.38 × 10<sup>-2</sup></td><td>11.88</td></tr><tr><td>Mm512r</td><td>256</td><td>256</td><td>512</td><td>10<sup>-3</sup></td><td>Random</td><td>3 × 10<sup>-2</sup></td><td>1.38 × 10<sup>-2</sup></td><td>23.75</td></tr><tr><td>Mh256r</td><td>256</td><td>256</td><td>256</td><td>10<sup>-3</sup></td><td>Random</td><td>10<sup>-1</sup></td><td>2.43 × 10<sup>-2</sup></td><td>20.79</td></tr><tr><td>Mh512r</td><td>256</td><td>256</td><td>512</td><td>10<sup>-3</sup></td><td>Random</td><td>10<sup>-1</sup></td><td>2.43 × 10<sup>-2</sup></td><td>41.57</td></tr><tr><td>H512m1</td><td>256</td><td>256</td><td>512</td><td>10<sup>-3</sup></td><td>m = 1</td><td>0</td><td>0</td><td>0</td></tr><tr><td>M1512m1</td><td>256</td><td>256</td><td>512</td><td>10<sup>-3</sup></td><td>m = 1</td><td>10<sup>-2</sup></td><td>8.06 × 10<sup>-3</sup></td><td>13.89</td></tr></tbody></table>

We adopt Kerr–Schild coordinates (see the Appendix) to allow for an inner radial boundary inside the black hole event horizon located at a radius $r_h = 2$, thus preventing numerical artefacts due to boundary effects that could otherwise propagate through the domain and affect the simulation at $r > r_h$. In radial direction, the numerical domain ranges from $r_{min} = 0.97 r_h$ to $r_{max} = 100$ covered by $N_r = 256$ grid points, with outflow boundary conditions (zeroth-order extrapolation) applied at both radial extrema. The radial mesh is non-uniform to increase the resolution towards the black hole event horizon by defining our radial grid points $r_i$ as

$$r_i = r_{\min} + \frac{r_{\max} - r_{\min}}{\varepsilon} \tan[\arctan(\varepsilon)x_i], \quad (2)$$

where $x_i = (i - 0.5)/N_r$ and the stretching parameter $\varepsilon$ is set to 10. The polar domain extends from $\theta_{\min} = 0$ to $\theta_{\max} = \pi$ with a resolution of $N_\theta = 256$ points. For simplicity, we impose axisymmetric reflection at the boundaries, since we do not expect our 3D models to be affected by the dynamics in the low-density regions close to the rotation axis. To better resolve the disc, the polar mesh is refined towards the equatorial mid-plane by setting the polar grid points $\theta_i$ to

$$\theta_i = \frac{\pi}{2}[1 + (1 - \zeta)(2y_i - 1) + \zeta(2y_i - 1)^n], \quad (3)$$

where $y_i = (i - 0.5)/N_\theta$, $\zeta = 0.6$, and $n = 29$. This gives a roughly constant and fine grid spacing across the disc and a rapidly decreasing resolution towards the rotational axis (Noble, Krolik & Hawley 2010). Finally, we consider the full azimuthal range $\phi \in [0, 2\pi]$ with uniformly distributed cells and periodic boundaries to be able to resolve global azimuthal modes with mode number $m = k_\phi r = 1$, which are expected to develop and to be also the fastest growing modes for the PPI in our disc model.

An important aspect to consider in any numerical experiment is its convergence, i.e. whether or not the results depend on the grid resolution. For a simulation involving magnetized accretion flows, the key aspect that needs to be properly resolved is the MRI turbulence that appears whenever a differentially rotating fluid is threaded by a magnetic field of any topology. Following Hawley, Guan & Krolik (2011), we define a quality metric as the ratio of the characteristic wavelength of the MRI mode $\lambda_{\text{MRI}} = 2\pi|u_A|/\Omega$ (which corresponds to the distance travelled by an Alfvén wave during an orbital period) and the grid zone size, where $u_A = B/\sqrt{\rho h + B^2}$ is the relativistic Alfvén velocity, $h$ is the specific enthalpy, $\Omega = u^\phi/u^\prime$

is the disc’s angular velocity, and $u^\mu$ is the fluid 4-velocity. Since we start with a purely toroidal magnetic field, the relevant quality metric should consider the wavelength along the $\phi$ direction, that is

$$Q_{\phi} = \frac{\lambda_{\text{MRI}}}{\Delta\phi\sqrt{\gamma_{\phi\phi}}} = \frac{2\pi|u_{\Delta\phi}|}{\Omega\Delta\phi\sqrt{\gamma_{\phi\phi}}}. \quad (4)$$

The volume average of $Q_\phi$ for each magnetized model at the beginning of the simulation is reported in Table 2. Hawley et al. (2011) suggest that $Q_\phi \gtrsim 20$ should provide a sufficiently good description of the non-linear phase of MHD turbulence. Note, however, that the recent stratified shearing box simulations of Ryan et al. (2017) suggest that none of the current simulations may actually be converged, even at much higher resolution than achievable in a global model. For all our simulations, we use the Harten–Lax–van Leer Riemann solver instead of the more dissipative Lax–Friedrichs scheme, together with a PPM reconstruction scheme.

To trigger the growth of non-axisymmetric modes, we introduce inside the torus a small perturbation $\delta v^\phi$ of the equilibrium azimuthal velocity $v_0^\phi$, with either random noise or cosine waves of the form:

$$\delta v^{\phi} = A v_{0}^{\phi} \cos(m\phi), \quad (5)$$

with $m = 1, ..., 5$ and amplitudes $A$ ranging from $10^{-6}$ to $10^{-2}$ depending on the simulation.

### 3 DIAGNOSTICS

We now introduce the quantities that we calculate from each simulation to probe the dynamical evolution of the models and determine the relative importance of the PPI and MRI in these models.

#### 3.1 Power of azimuthal modes

In analogy with De Villiers & Hawley (2002) and Wielgus et al. (2015), for any given azimuthal number $m$ we calculate the power contained in an azimuthal mode for a generic quantity $Q$ as

$$P_{m,Q}(t) = \frac{\int_{\text{rin}}^{\text{rout}} \int_{0}^{\pi} \left| \frac{1}{2\pi} \int_{0}^{2\pi} Q e^{im\phi} d\phi \right|^2 w(r, \theta) \sqrt{\gamma} d\theta dr}{\int_{\text{rin}}^{\text{rout}} \int_{0}^{\pi} w(r, \theta) \sqrt{\gamma} d\theta dr}, \quad (6)$$
where $w(r, \theta)$ is an axisymmetric weight function. For $Q = u_A^\phi$, the integrals in radial and polar direction in equation (6) are weighted by the rest mass density, i.e.

$$w = \langle \rho \rangle_{\phi} = \frac{1}{2\pi} \int_{0}^{2\pi} \rho \, d\phi, \qquad (7)$$

where the operator $\langle\rangle_\phi$ represents an azimuthal average. This choice avoids overestimating the contribution of the rarefied atmosphere enveloping the disc and at the same time leads to a quantity that relates to the azimuthal components of magnetic energy. For $Q = \rho$, we set $w = 1$. Since we are interested mostly in the relative importance of the $m = 1$ mode with respect to other higher order non-axisymmetric modes, we decided to normalize the mode power in equation (6) at any given time by $P_{0,Q}(t)$, i.e. by the instantaneous power in the axisymmetric mode. By doing so, we factor out any decrease of the azimuthal power introduced by the significant mass-loss occurring in all the models we considered. The evident drawback of this time-dependent normalization is a loss of information in absolute terms on the evolution of the power within a specific azimuthal mode, once a non-negligible fraction of the disc's mass has been accreted on to the black hole. However, we verified that this normalization did not affect the estimate of growth rates and the value at which the modes power saturates, since most mass is accreted during later stages.

With these diagnostics, we estimate growth rates and saturation levels. We also quantify the time evolution of the relative importance of the various modes resolved by the numerical simulations by constructing spectrograms. At any given time $t$, we compute the power in equation (6) of modes with azimuthal number up to $m = 50$, and we plot the power in a $m$ versus time diagram. By time-averaging over the full duration of the simulations, we also compute spectra to characterize the power distribution across different modes.

We display information on the frequency components present in the fastest growing PPI mode with frequency-radius diagrams. We consider the complex amplitude of the $m = 1$ mode of the density in the equatorial plane (since most of the dynamics takes place in this region):

$$M(r, t) = \frac{1}{2\pi} \int_{0}^{2\pi} \rho(r, \pi/2, \phi) e^{i\phi} d\phi. \qquad (8)$$

For each radius $r$, we compute $\tilde{M}(r, \omega) = FFT(M(r, t))$, i.e. the Fourier Transform in the frequency domain.

## 3.2 Turbulence and accretion

We keep track of the development of turbulence in the system by considering the evolution of the $r-\phi$ components of the Reynolds and Maxwell stress tensors, defined, respectively, as

$$W_{\text{Re}} = \rho \delta u^r \delta u^\phi \sqrt{\gamma_{rr}} \sqrt{\gamma_{\phi\phi}}, \qquad (9)$$

$$W_{\text{Ma}} = B^r B^\phi \sqrt{\gamma_{rr}} \sqrt{\gamma_{\phi\phi}}, \qquad (10)$$

where $\delta u^r = u^r - \langle u^r \rangle_\phi$ and $\delta u^\phi = u^\phi - \langle u^\phi \rangle_\phi$. We compute their volume averages by considering only those regions of the computational domain where the rest mass density $\rho$ exceeds a threshold value set to $\rho_{\text{th}} = \sqrt{\rho_c \rho_{\text{atm}}}$ to track the dynamics of the disc and exclude that of the atmosphere.

Still related to the stresses, we compute the disc alpha parameter (not to be confused with the lapse function) as the ratio of the volume

average of the total stress $W_{\text{tot}} = W_{\text{Re}} + W_{\text{Ma}}$ and the volume average of the thermal pressure:

$$\alpha_{\text{turb}} = \frac{\langle W_{\text{tot}} \rangle_V}{\langle p \rangle_V}. \qquad (11)$$

As a further diagnostics of the efficiency of angular momentum transport in the disc, and thus of the overall accretion process, we also monitor the evolution of the radial distribution of the disc's orbital angular velocity $\Omega$. The radial dependence of $\Omega$ is usually described with a power law:

$$\Omega \propto r^{-q}, \qquad (12)$$

where the parameter $q$ can range from 3/2 (for a Keplerian disc) to 2 (constant specific angular momentum). However, in the relativistic case and for a non-rotating black hole $\Omega = -lg_{\text{tt}}/g_{\phi\phi} = l(r-2)/r^3$, and the value $q = 2$ can be assumed only if the disc extends sufficiently far away from the black hole. Since our disc model extends from $r_{\text{in}} = 6.16$ to $r_{\text{out}} = 21.6$, initially $q$ ranges from $q(r_{\text{in}}) = 1.52$ to $q(r_{\text{out}}) = 1.90$. Since we want to relate the redistribution of specific angular momentum to the classical case where the slope in equation (12) starts as a constant equal to 2 across the disc, instead of $q$ we monitor the evolution of the quantity

$$\tilde{q} = 2 - \left| \frac{d \log l}{d \log r} \right|, \qquad (13)$$

which is evaluated performing a least-squares fit of the power law describing the radial dependence of the specific angular momentum $l$.

# 4 RESULTS AND DISCUSSION

## 4.1 Hydrodynamic disc

First, we focus on the development and saturation of the PPI in the absence of magnetic fields to have an initial benchmark for a later comparison with the results for magnetized models. Fig. 1 illustrates the core mechanism responsible for the PPI growth, but a more detailed and complete discussion on the topic can be found in Narayan & Goodman (1989). As shown in Table 2, we performed a set of simulations that differ by resolution in the azimuthal direction, and amplitude and spectrum of the initial perturbation. We followed the evolution of these models up to 20 orbital periods at the disc centre, which is sufficient for the hydrodynamic instability to reach saturation in terms of azimuthal mode power.

Our results confirm that the $m = 1$ azimuthal mode is the fastest growing one. Therefore, it was selected by the system independently of the initial perturbation. We ran a series of simulations with the same monochromatic $m = 1$ perturbation but with different resolutions in the azimuthal direction (models H32m1-4 to H256m1-4). Fig. 2 shows the time evolution of the azimuthal power of the density fluctuations for the $m = 1$ and $m = 2$ modes. Even with a modest resolution of 32 zones, ECHO is capable of capturing the dynamical evolution of the PPI, since the most unstable mode has a quite large wavelength. Fig. 2 also shows the linear phase of the instability during the first 10 orbital periods, with growth rates for the $m = 1$ and $m = 2$ modes of 0.080 $\Omega_c$ and 0.042 $\Omega_c$, respectively. The power of the two modes differs also in the value of saturation levels, which is more than an order of magnitude larger for the $m = 1$ mode.

The minimum resolution required to properly resolve the fastest growing mode of the PPI can be much higher in radial direction than in azimuthal direction (i.e. higher than 32 points), because the $m = 1$ mode developing from wide tori is not (as in the case
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Time evolution of the power in density for the *m* = 1 (solid curves) and *m* = 2 (dashed curves) modes, as computed from equation (6) for models H32m1-4 to H256m1-4. All models are initialized with an *m* = 1 perturbation.</figcaption></figure>

of slender tori) the so-called *principal mode*. This mode is the result of the interaction of two nodeless (across the radius) surface gravity waves (also referred to as *edge-waves*) that are generated at the disc's inner and outer boundaries and are advected by the shear flow with respect to the corotation radius. While this mode is essentially incompressible, the one that we see in our simulations is the outcome of the interaction between a pressure wave in the outer part of the disc and an edge wave from the interior. The former has multiple nodes in radial direction (two, in the case of our models), which require an adequate radial resolution. If one has insufficient resolution, one systematically underestimates the instability growth rates [see Blaes & Hawley (1988) for a detailed discussion].

Fig. 3 (left-hand panel) shows an equatorial slice of the rest mass density for model H256r after about 15 orbital periods. The dominant *m* = 1 mode is clearly visible as an overdensity that corotates with the disc, while the flow still maintains overall a smooth profile. The region between the black hole horizon and

<figure><img src="image_5.png" /><figcaption>**Figure 4.** Time evolution of the power in density for the *m* = 1 (solid curves) and *m* = 2 (dashed curves) modes for models H256r (black), M1512r (blue), Mm512r (green), and Mh512r (red). All models are initialized with a random perturbation.</figcaption></figure>

the disc's inner edge is relatively depleted of mass, apart from an inspiraling flow that detached from the main body of the disc.

From the time evolution of the mode power in Fig. 4 and the spectrogram in Fig. 5 (top panel), it is clear that the *m* = 1 mode dominates since very early times, and no other small-scale perturbation grows as much during the linear phase of the instability. After 10 orbital periods, a further deposition of energy occurs on smaller scales presumably because of non-linear interactions, but the *m* = 1 mode remains the strongest one. This interpretation is confirmed by the density power spectrum (the top panel of Fig. 6), which shows also most power in low-order modes, peaking at *m* = 1, and decaying as a power law $m^{-4}$ for *m* ≲ 3.

All these results are in full agreement with the predictions of the linear theory (Goldreich et al. 1986; Glatzel 1987) as well as previous numerical simulations (Blaes & Hawley 1988; De Villiers & Hawley 2002).

<figure><img src="image_6.png" /><figcaption>**Figure 3.** Equatorial cut of the rest mass density $\rho$ for models H256r (left) and M1512r (right) at $t = 3000 \simeq 15 P_c$. The maximum value of $\rho$ is normalized to 1 in each plot. The solid black curve represents the black hole event horizon, while the dotted curve indicates the radius of the last marginally stable orbit $r_{ms}$.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Rest mass density spectrograms of the azimuthal mode power, colour-coded in logarithmic scale. The panels refer, from top to bottom, respectively to models H256r, M1512r, Mh512r, H512m1, and M1512m1.</figcaption></figure>

## 4.2 Magnetized discs

We now consider the effect of a weak toroidal magnetic field. We use the analytical solution provided by Komissarov (2006), which allows one to avoid initial spurious transients by taking into account the extra pressure provided by a toroidal magnetic field. Recent models by Fragile & Sądowski (2017) suggest that strongly magnetized thick discs with a purely toroidal magnetic field may experience a fast drop of magnetization, mostly due to redistribution of gas and migration of magnetic fields into the funnel region. As we are more interested in the dynamics taking place within the torus, we limit ourselves to values of central magnetization up to $\sigma_c = 0.1$.

We focus for the moment on the magnetized models with the highest resolution along the $\phi$ direction (i.e. 512 grid points), leaving the discussion of the effects of a lower resolution to the next paragraph. As illustrated by the right-hand panel of Fig. 3, the distribution of rest mass density in the magnetized models show small-scale fluctuations not present in the hydrodynamic models, with the clear presence of MHD turbulence triggered by the MRI and the absence of a large-scale overdense non-axisymmetric structure.

The time evolution of the azimuthal mode power in the magnetized models (Fig. 4, blue, green, and red curves with increasing magnetic field strength) reveals an earlier growth of low order modes, without a clear distinction between the $m = 1$ and $m = 2$

<figure><img src="image_8.png" /><figcaption>**Figure 6.** Rest mass density (top) and azimuthal Alfvén velocity (bottom) spectra in azimuthal number $m$, averaged over the time $t \in \{0, 13.6\} P_c$. The black curve represents the hydrodynamic model, the other solid curves refer to magnetized models with high resolution along the $\phi$ direction (i.e. M1512r, Mm512r, and Mh512r), while the dashed ones stand for the low-resolution models M1256r, Mm256r, and Mh256r.</figcaption></figure>

modes as in the hydrodynamic case. The mode power reaches a maximum, which occurs at earlier times and at a slightly higher value than the magnetic field strength, although both the $m = 1$ and $m = 2$ modes still saturate at roughly the same level. This behaviour is confirmed also by the spectrograms in Fig. 5, which show a much broader range of excited modes in the magnetized cases (second and third panel from the top), and provides a picture which is fully consistent with the onset of the linear phase of the non-axisymmetric MRI (Balbus & Hawley 1992; Wielgus et al. 2015).

The time-averaged density spectra shown in the top panel of Fig. 6 provide a more quantitative confirmation of this trend. The hydrodynamic model H256r (black curve) develops predominantly large-scale modes, with a steep power-law decline from the $m = 1$ mode to $m \approx 5$, followed by a shallower decline up to $m = 10$ and again a steep drop. The magnetized models (coloured curves) have no strong excess of power at large scales, with a much shallower slope down to $m = 10$ followed by a steeper decrease due to the numerical dissipation introduced by the reconstruction scheme. The spectra of the high-resolution models (coloured solid curves) behave quite similarly to each other. For $\sigma_c = 0.01, 0.03$, the spectra computed from the orbital Alfvén velocity show no particular difference in their shape (the bottom panel of Fig. 6), as they overlap
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Reynolds and Maxwell stresses (solid and dotted curves, top panel) and $\alpha_{\text{turb}}$ (bottom panel) for models H256r, M1512r, Mm512r, and Mh512r (black, blue, green, and red as in the previous figure).</figcaption></figure>

quite well for both small and large azimuthal numbers. The high-magnetization model ($\sigma_c = 0.1$) exhibits a slightly steeper profile and also less power, which is probably due to the earlier transport of magnetic flux in the disc's atmosphere. In fact, by choosing as weighting function in equation (6) the rest mass density, a higher concentration of magnetic field in the rarefied atmosphere leads to an underestimate of the mode's power.

All the models we considered produce accretion on to the central black hole, but at different times and in different ways. In the hydrodynamic disc the development of the PPI is the sole responsible for angular momentum transport, hence accretion. After 10 orbital periods, when the $m=1$ mode approaches its maximum amplitude, the kinetic energy and stresses are large (the black curve in Fig. 7), which consequently lead to a significant redistribution of angular momentum (see the evolution of the parameter $\tilde{q}$ in Fig. 8) and mass-loss (almost 30 per cent of the disc's initial mass).

The situation is quite different in the magnetized models. The accretion is triggered at much earlier times (after only 2 orbital periods for the model with highest magnetization) with a steeper increase of $W_{\text{tot}}$. In model M1256r, the Reynolds and Maxwell stresses increase initially at the same growth rate, but then the magnetic component takes over. For model Mh256r, Maxwell stresses dominate from the very beginning and during the whole simulation. Consequently, accretion is enhanced and a higher mass-loss results, leading in the most highly magnetized model to a dramatic decrease of the disc mass down to 30 per cent of its initial value. Concerning the angular momentum distribution, Fig. 8 confirms these findings, showing a

<figure><img src="image_10.png" /><figcaption>**Figure 8.** Slope parameter $\tilde{q}$ for models H256r, M1512r, Mm512r, and Mh512r.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Figure 9.** Space-time diagrams of the azimuthally averaged radial distribution of specific angular momentum $l$ in the equatorial plane for models H256r (top), M1512r (middle), and Mh512r (bottom).</figcaption></figure>

faster and earlier decrease of $\tilde{q}$ (see equation 13) and also a lower saturation value with respect to the hydrodynamic case.

A space-time diagram of the radial profile of the specific angular momentum $l$ at the equator (the top panel in Fig. 9) clearly shows that in the hydrodynamic model the waves that constitute the unstable mode are increasing $l$ outside the corotation radius (starting from
the outmost parts of the disc) and at the same time decreasing it in the inner region of the disc. During this phase (up to $t \sim 13 P_c$), the angular momentum in the disc's central region does not change, since its transport is carried out by two coupled waves propagating in different parts of the torus. This reveals once again the inherently non-local nature of the PPI. The same diagrams for models M1512r and Mh512r (the middle and bottom panel of Fig. 9) show how in the magnetized case angular momentum is transported outwards from the inner regions of the disc much faster than in the absence of a magnetic field. The transport proceeds in this case from the highly magnetized regions of the disc towards the outer edge, with the angular momentum being deposited at increasingly larger radii, and hence producing a rather different pattern in the space-time diagram.

## 4.3 Mode frequency

The mode selected by the PPI as the fastest growing one is characterized not just by the azimuthal number $m=1$, but also by a specific angular frequency $\omega$, which we can measure and directly compare with the disc rotation rate $\Omega$. In fact, the mechanism behind the PPI growth (see once again Fig. 1) relies on the interaction between a negative energy wave (which has $\omega < \Omega$ and therefore is slower than the surroundings) and a positive energy wave (with $\omega > \Omega$). This interaction happens at the corotation radius $R_{\text{cor}}$ (defined by $\Omega(R_{\text{cor}}) = \omega$), so that negative energy waves from the inner region of the disc ($r < R_{\text{cor}}$) exchange energy and angular momentum with the positive energy ones from the outer parts ($r > R_{\text{cor}}$). A fundamental condition for this process to be efficient is the presence of a reflective boundary (Goldreich et al. 1986) that allows the waves to be reflected and approach the corotation radius where the interaction can take place, leading to a positive feedback and hence the development of the PPI (see Fig. 1).

The top panel of Fig. 10 clearly shows the presence of a specific spectral component at $\omega_0 \simeq 0.86\Omega_c$ (where $\Omega_c = 2\pi/P_c$ is the orbital angular frequency at the disc's centre) in the hydrodynamic model, which represents the $m=1$ mode selected by the PPI. Both positive and negative energy waves are present, and they interact through the corotation radius $R_{\text{cor}} \simeq r_c$, defined as the location where the mode frequency and the disc's orbital frequency are equal (the white curve in Fig. 10). The two waves cannot propagate through a narrow forbidden region surrounding the corotation radius, but they can still be transmitted by tunnelling.

In the magnetized case with high resolution in azimuthal direction (bottom panel), there is no clear selection of a single mode at a well-defined frequency. In the region beyond the corotation frequency the positive-energy waves are much less excited, and there seems to be a lack of modes with frequency below $\omega \sim 0.5 \Omega_c$. We also observed no significant dependence on the model's magnetization: models M1512r, Mm512r, and Mh512r all present a similar pattern in the $\omega$-radius diagram.

## 4.4 Dependence on resolution

We focus now on how the results we considered so far change when a lower resolution in azimuthal direction is employed in magnetized models (256 instead of 512 grid points).

As we can see from the last column of Table 2, some magnetized models have an average value of the quality metric $Q_{\phi} \lesssim 20$, suggesting that in these cases the MHD turbulence might not be sufficiently well resolved. Moreover, considering that the resolution constraint gets stricter when weaker magnetic field are employed,

<figure><img src="image_12.png" /><figcaption>**Figure 10.** Amplitude of the $m=1$ mode frequency components in a $\omega$ versus radius diagram calculated from the rest mass density $\rho$ and computed from the Fourier transform of (8) within the time interval $t \in \{0, 20\} P_c$. The white curve represents the azimuthally and time-averaged equatorial orbital frequency profile. The panels refer, respectively, to model H256r (top), M1256r (middle), and M1512r (bottom).</figcaption></figure>

one could also expect the low magnetized discs to exhibit greater differences between the two different grids (with respect to the models with higher magnetization).

This expectation is confirmed by Fig. 6. The rest mass density of model M1256r (the dashed blue curve in the top panel) shows a small excess of power in the $m=1$ mode, as the ratio $P_{1,\rho}/P_{2,\rho} \sim 3$ is roughly twice as large as in the corresponding high-resolution run. Moreover, at all $m$, the power is systematically lower than in the model with the finer grid, pointing to a less efficient MHD-driven mode excitation. This is corroborated by the bottom panel of Fig. 6, where the power in the spectrum of the azimuthal Alfvén velocity decreases by almost a factor of 2 on the coarser grid for all values of $m < 10$. The middle panel of Fig. 10 shows that at lower resolution a still significant excitation of positive-energy modes is present at frequency $\omega \sim 0.7 \Omega_c$, while there is almost none when a finer grid is employed (bottom panel). The model with low magnetization and $N_{\phi} = 256$ displays also a later onset of accretion and growth of stresses, consistent with a less efficient MHD turbulence.

Both models with higher magnetization ($\sigma_c = 0.03, 0.1$) show no such a significant dependence on the azimuthal resolution. The main difference visible from the averaged spectra in Fig. 6 is in the slopes at $m \gtrsim 10$, which can be ascribed to the increase in numerical dissipation on the coarser grid for a given accuracy of the reconstruction scheme.

Both the initial value of $Q_{\phi}$ and the spectral behaviour of the magnetized models suggest that the fastest growing mode of the MRI is underresolved in the low-resolution, low-magnetization model
<figure><img src="image_13.png" /><figcaption>Figure 11. Time evolution of the power in density for the $m = 1$ (solid curves) and $m = 2$ (dashed curves) modes for models H256r (black), M1512r (blue), H512m1 (green), and M1512m1 (red). The former two models are initialized with a random perturbation, while the others are initialized with an $m = 1$ perturbation.</figcaption></figure>

M1526r. In these models, the MHD turbulence is not well developed, and the PPI appears to be still capable of exciting to some extent the large-scale $m = 1$ mode.

Since the non-axisymmetric fastest growing modes selected by the MRI are typically associated with short vertical and radial wavelengths (Balbus & Hawley 1998), we also expect a dependence of our results on the resolution along the $\theta$-direction. Indeed, running the three magnetized models with a polar resolution decreased by a factor of 2 leads to results qualitatively similar to those obtained at lower azimuthal resolution. Underresolving the MHD turbulence produces a small excess of power in the low order modes and a general drop in power for modes with $m \gtrsim 5$.

## 4.5 Dependence on the initial perturbation

So far we considered the evolution of models that were started with a random perturbation of the azimuthal velocity, therefore exciting all modes with a low enough azimuthal number that could be resolved by the number of points in our numerical grid. Our high-resolution models are able to resolve the MRI fastest growing mode, which has azimuthal number typically much larger than 1, i.e.

$$m_{\text{MRI}} = \frac{2\pi R}{\lambda_{\text{MRI}}} \approx \frac{R\Omega}{u_{\text{A}}^{\phi}} \quad (14)$$

has an average value of $\langle m \rangle_V \sim 36$ for the models with magnetization $\sigma_c = 0.01$. Hence, a random perturbation will equally excite the PPI and the MRI fastest growing modes. However, the growth rate of the fastest MRI mode exceeds the typical growth rate of the hydrodynamic PPI, which plays a key role in the results from the previous sections (Hawley 2000).

It is not clear, though, whether the PPI could leave a clear signature on the disc structure if the initial perturbation were to favour the PPI fastest growing mode over the higher order ones excited by the MRI. To assess this, we initialized some of our models with a purely $m = 1$ perturbation, hence exciting from the beginning of the simulation the most unstable PPI mode.

When the hydrodynamic disc is initially excited with a $m = 1$ mode (model H512m1), its evolution resembles quite closely that of the randomly perturbed disc (see Fig. 11), with the growth of

the fastest growing mode's power peaking at much earlier times. Hence, there is also an earlier non-linear interaction between the $m = 1$ and higher order modes, leading to episodic redistributions of power across scales (compare the first and the fourth panel from the top in Fig. 5).

The other model we excited with a $m = 1$ perturbation is the low magnetized one, i.e. model M1512m1. From Fig. 11, we can see how in this case (red curve) the $m = 1$ mode develops during the first four orbital periods almost exactly as in the hydrodynamic counterpart (green curve). However, the mode's power starts then to drop and approaches the values reached in the randomly perturbed magnetized model (blue curve), i.e. it is damped by almost two orders of magnitude. The power in the $m = 2$ mode (dashed curves) is affected in a similar way: its growth significantly slows down with respect to the hydrodynamic case after about five orbital periods, and it reaches values that are roughly a factor of 2 below the unmagnetized model's ones.

The spectrograms in Fig. 5 show that high-order modes (i.e. with $m > 2$) are excited by the MRI in both models M1512r and M1512m1 (second from the top and bottom panels) roughly at the same time $t \sim 4P_c$, which marks also the beginning of the departure from the hydrodynamic model's behaviour. The onset of the MRI appears therefore to be rather insensitive of the initial perturbation employed. Model M1512m1 still shows a distinctive excitation of the $m = 1$ mode (contrary to the randomly perturbed case), which however disappears after $t \sim 10P_c$. There is no trace of the episodic excitation of the modes with $2 < m < 7$ that takes place in the unmagnetized disc between $t = 6P_c$ and $t = 10P_c$, but instead the power is more uniformly distributed across high-order modes, qualitatively resembling the behaviour of model M1512r.

## 4.6 Interaction between the two instabilities

From the previous analysis, it seems clear that the inclusion of a toroidal magnetic field deeply affects the development of the PPI, even when the field is highly sub-thermal. The $m = 1$ mode selected by the PPI is suppressed, leaving no clear signature on the azimuthal or frequency spectra.

One possible explanation could reside in the early accretion triggered by MHD turbulence, which leads to a faster redistribution of angular momentum across the disc. The free energy stored in the shearing flow decreases and the edge waves that should reach the corotation radius and transport energy are not efficient enough to sustain the unstable $m = 1$ mode. Hawley (2000) found that a torus threaded by a weak constant toroidal magnetic field (i.e. with magnetization $\sigma \sim 0.01$) can display an initial appearance of the PPI, though the later development of the MRI interrupts the growth of the hydrodynamic instability. The author suggests the redistribution of angular momentum produced by the MRI to be the main cause of the inability of the PPI to significantly affect the disc dynamics.

A comparison between model H512m1 and M1512m1 provides new insights in the dynamical quenching of the PPI by the MRI. For both models, the power in the $m = 1$ mode (Fig. 12) reaches its peak value at a similar time (roughly $t \sim 8P_c$), while the parameter $\tilde{q}$ decreases in both models at the same rate. This is consistent with the usual interpretation that the $m = 1$ mode ceases to be excited when $\tilde{q}$ reaches a sufficiently low value, which sets an unfavourable environment for the development of the PPI.

While it is true that in model M1512m1, the transport of angular momentum after $t \sim 10P_c$ further proceeds due to the MRI (hence leading $q$ to the same saturation value as the randomly perturbed model), the strong decrease in the $m = 1$ mode power between
<figure><img src="image_14.png" /><figcaption>**Figure 12.** Slope parameter $\tilde{q}$ for models H256r, M1512r, H512m1, and M1512m1.</figcaption></figure>

$t=8$ and $11$ cannot be attributed to this effect because $\tilde{q}$ has not yet dropped below the value of the corresponding hydrodynamical model. Furthermore, the first deviations of the $m=1$ mode from the hydrodynamic model appear even earlier at $t \sim 4 P_c$, at the same time as the rise of high-order modes excited by the MRI. We therefore argue that it is an early coupling between low- and high-order modes, and not the transport of angular momentum, that leads to the increasing difference in the $m=1$ mode's power between models H512m1 and M1512m1: energy is redistributed towards small scales, and the PPI fastest growing mode is damped until its power reaches the same value as in the randomly perturbed models. Hence, the disc evolves to a state which is ultimately independent of the initial perturbations: even when the PPI most unstable mode is initially enhanced, its later interaction with the higher order modes excited by the MHD instability ends the initial transient growth. The same coupling is also present in the hydrodynamic case, but since the high-order modes are less excited by the PPI than by the MRI, the resulting non-linear interaction is weaker and only leads to oscillations in the $m=1$ mode's power (see the first and fourth panels of Fig. 5).

To confirm that in presence of stronger magnetic fields, the disc dynamics is completely dominated by the MRI, we performed a run with the exact set-up used in model Mh256r but reduced the azimuthal range to $[0, \pi]$ to filter out the $m=1$ mode from the system. We observed no appreciable change in the evolution with respect to the results obtained using the full azimuthal range $[0, 2\pi]$. Hence, we conclude that in this case the dynamics of the system is independent of the large-scale $m=1$ mode, which does not play any significant role in the evolution of our most strongly magnetized model. Powerful accretion and small-scale turbulence set by the MRI prevent any significant development of the PPI fastest growing mode.

These results may seem to be in disagreement with what was reported by Fu & Lai (2011), who performed an analytical study of the influence of magnetic fields (in both toroidal and poloidal configurations) on the stability on accretion tori. They found that a sufficiently strong toroidal magnetic field can further destabilize the disc and enhance PPI development. However, they assumed incompressibility in their analysis, and therefore they excluded all those modes that, as in the case of the wide tori we considered, result from the interaction between a pressure wave outside the corotation radius coupled to an internal edge wave [see Glatzel

(1987) for a more complete analysis of the role of compressibility in the development of PPI].

## 5 CONCLUSIONS

We presented 3D GRMHD simulations of magnetized thick accretion discs around black holes with the goal of estimating the interaction between the PPI, which deeply affects hydrodynamic wide tori with constant specific angular momentum, and the MRI, which is omnipresent in ionized differentially rotating astrophysical discs threaded by magnetic fields. By starting from the magnetized equilibrium solution provided by Komissarov (2006) we could avoid initial transients in the disc's dynamics and assess how the linear growth and subsequent dynamics of the PPI is affected by the development of the MRI.

Consistently with previous works, we find that even in the presence of a sub-thermal toroidal field the growth of the non-axisymmetric $m=1$ mode (usually the fastest growing mode for wide hydrodynamic tori) is significantly quenched, i.e. no large-scale overdensity structure (planet) forms and the flow is more turbulent.

The inclusion of toroidal magnetic fields excites higher order modes, thus smaller length-scales, in contrast of the strong dominance of the $m=1$ mode in the purely hydrodynamic case. For all our high-resolution magnetized models, there is no clear sign of an excess of power in the $m=1$ mode, showing therefore that its growth has been suppressed by the action of MRI. The redistribution of angular momentum proceeds much faster than in the unmagnetized model, resulting in a lower saturation value for the parameter $\tilde{q}$ measuring the slope of the orbital angular velocity with radius. This is consistent with the fact that, as expected, higher magnetizations lead to stronger stresses and accretion rates on to the black hole.

Launching the model with the weakest magnetic field ($\sigma_c = 10^{-2}$) with a monochromatic $m=1$ azimuthal perturbation leads to a transient phase, where the PPI fastest growing mode behaves in the same way as in the hydrodynamic counterpart. For more than five orbital periods, even in presence of magnetic fields, the disc's dynamics is regulated by the PPI, with the power in the $m=1$ mode that greatly exceeds that of any other azimuthal mode. Then, the power drops and reaches in a very short time the same saturation level as in the randomly perturbed magnetized models. While the ultimate fate of the large-scale mode appears to be the same, independently of the particular spectrum of the perturbation, we observe nevertheless a transient PPI-dominated phase that lasts for a few orbital periods and that could be interesting in those astrophysical scenarios where the gravitational interaction with the central black hole can indeed excite the $m=1$ mode.

The comparison between model H512m1 and M1512m1 offers a very interesting insight into the physical mechanism behind the suppression of the PPI. The redistribution of specific angular momentum from a constant radial profile towards the Keplerian one occurs for both the PPI and the MRI, and it has been invoked in previous studies (Hawley 2000) as the sole responsible of the missing development of the PPI in 3D MHD global simulations, since the hydrodynamic instability is intrinsically susceptible to being stabilized if enough angular momentum is transported outwards in the disc. Our results are consistent with this mechanism being responsible for the halt in the growth of the PPI most unstable mode, since in both models H512m1 and M1512m1 the $m=1$ mode's power reaches a maximum at similar times, while displaying the same decrease in the parameter $\tilde{q}$. However, our simulations also show that
the general lack of a dominant $m = 1$ mode in magnetized discs is mainly due to the coupling with higher order modes that drains power from large scales and redistributes it towards smaller scales. Such a coupling is rather strong in the magnetized case (since the MRI excites modes with large $m$ more than the PPI does), and we identify it as the main mechanism that quenches the growth of the PPI. This process leads to a turbulent state completely dominated by the MRI and ultimately independent of the initial perturbation.

The initial PPI dominated evolution of the magnetized model M1512m1 is consistent with the findings of Nealon et al. (2018), who argued that a thick disc can undergo a significant accretion phase dominated by the PPI, before the onset of the MRI. Their analysis is based on numerical measurements of accretion rates in the evolution of the remnant of a hydrodynamic tidal disruption event around a supermassive black hole and a comparison between the expected duration of the linear phases of the PPI and MRI. They argue that for weak enough magnetic fields the saturation time of the MRI will significantly exceed that of the PPI, suggesting that in that case the disc could undergo an accretion phase regulated mostly by the PPI. However, they also notice that using this argument to discern between a PPI and MRI dominated regime can lead to a significant overestimate of the magnetic field strength required for the latter scenario. Indeed, we found that the quenching of the PPI fastest growing mode starts well before the saturation of either of the instabilities.

Our results suggest that with the inclusion of magnetic fields the dominance of the $m = 1$ mode (with respect to higher order ones) due to PPI should not hold in thick tori, but there are some caveats that need to be addressed. First, we neglected the self-gravity of the disc. This is expected to be a good approximation, for example, in X-ray binaries, but not in the case of the remnant of an NS–NS merger that produces a black hole–torus system. It has been shown by Korobkin et al. (2011) that the gravitational interaction between the disc and the central black hole can indeed further excite the $m = 1$ mode, leading hence possibly to a different outcome than ours once a magnetic field is included into the simulations.

In this work, we considered tori with a constant specific angular momentum radial profile, which are most susceptible to PPI and therefore represent the most favourable environment for its development. Despite the general suppression by the MRI, we find an initial transient phase where the PPI dominates the evolution of a magnetized disc excited with a weak $m = 1$ perturbation. It would be interesting to assess whether this would still be the case when using magnetized equilibrium solutions with a more realistic distribution of angular momentum, such as those recently presented in Gimeno-Soler & Font (2017).

Another important aspect is the role of turbulent resistivity. We found that if the dynamical evolution of the MRI is not resolved properly, the system shows a residual excess of power in the $m = 1$ mode and a clear signature of negative- and positive-energy waves coupled through the corotation radius. This suggests that for a strong enough magnetic dissipation, which could be present in a turbulent environment, the MRI could be quenched, and hence the PPI could still grow significantly and produce a non-negligible $m = 1$ like overdensity. As mentioned in Section 2, in a forthcoming work we will investigate how the system is affected by an explicit magnetic diffusivity.

## ACKNOWLEDGEMENTS

The authors would like to thank C. Gammie, O. Blaes, O. Sadowski, M. Wielgus, R. Nealon, and T. Font for providing helpful comments

and suggestions. A special thank to the anonymous referee, whose comments helped us improving this manuscript.

The authors gratefully acknowledge the Leibniz Supercomputing Centre (LRZ; http://www.lrz.de) and the Max Planck Computing and Data Facility (MPCDF; http://www.mpcdf.mpg.de) for funding this project by providing computing time, respectively, on the SuperMUC (project PR62LU) and Hydra clusters. MB thanks, in particular, F. Baruffa and M. Rampp for their invaluable assistance and support during the development of the parallelization scheme used in ECHO.

JG acknowledges support from the Max-Planck-Princeton Center for Plasma Physics and from the European Research Council (grant MagBURST–715368).

## REFERENCES

* Abramowicz M. A., Fragile P. C., 2013, Living Rev. Relativ., 16, 1

* Abramowicz M., Jaroszynski M., Sikora M., 1978, A&A, 63, 221

* Aloy M. A., Müller E., Ibáñez J. M., Martí J. M., MacFadyen A., 2000, ApJ, 531, L119

* Arlt R., Rüdiger G., 2001, A&A, 374, 1035

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

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

* Balbus S. A., Hawley J. F., 1998, Rev. Mod. Phys., 70, 1

* Balbus S. A., Hawley J. F., Stone J. M., 1996, ApJ, 467, 76

* Blaes O., 2014, Space Sci. Rev., 183, 21

* Blaes O. M., Hawley J. F., 1988, ApJ, 326, 277

* Bucciantini N., Del Zanna L., 2011, A&A, 528, A101

* Bucciantini N., Del Zanna L., 2013, MNRAS, 428, 71

* Bugli M., Del Zanna L., Bucciantini N., 2014, MNRAS, 440, L41

* Coughlin E. R., Begelman M. C., 2014, ApJ, 781, 82

* De Villiers J.-P., Hawley J. F., 2002, ApJ, 577, 866

* De Villiers J.-P., Hawley J. F., 2003, ApJ, 592, 1060

* De Villiers J.-P., Hawley J. F., Krolik J. H., 2003, ApJ, 599, 1238

* De Villiers J.-P., Hawley J. F., Krolik J. H., Hirose S., 2005, ApJ, 620, 878

* Del Zanna L., Zanotti O., Bucciantini N., Londrillo P., 2007, A&A, 473, 11

* Del Zanna L., Bugli M., Bucciantini N., 2014, in Pogorelov N. V., Audit E., Zank G. P., eds, ASP Conf. Ser. Vol. 488, 8th International Conference of Numerical Modeling of Space Plasma Flows (ASTRONUM 2013). Astron. Soc. Pac., San Francisco, p. 217

* Del Zanna L., Papini E., Landi S., Bugli M., Bucciantini N., 2016, MNRAS, 460, 3753

* Fender R. P., Belloni T. M., Gallo E., 2004, MNRAS, 355, 1105

* Foucart F., Duez M. D., Kidder L. E., Scheel M. A., Szilagyi B., Teukolsky S. A., 2012, Phys. Rev. D, 85, 044015

* Fragile P. C., Sądowski A., 2017, MNRAS, 467, 1838

* Fragile P. C., Blaes O. M., Anninos P., Salmonson J. D., 2007, ApJ, 668, 417

* Fromang S., 2013, in Hennebelle P., Charbonnel C., eds, EAS Publ. Ser. Vol. 62, Role and Mechanisms of Angular Momentum Transport During the Formation and Early Evolution of Stars. EDP Science, Les Ulis, France, p. 95

* Fu W., Lai D., 2011, MNRAS, 410, 1617

* Gimeno-Soler S., Font J. A., 2017, A&A, 607, A68

* Glatzel W., 1987, MNRAS, 225, 227

* Goldreich P., Goodman J., Narayan R., 1986, MNRAS, 221, 339

* Hawley J. F., 1991, ApJ, 381, 496

* Hawley J. F., 2000, ApJ, 528, 462

* Hawley J. F., 2001, ApJ, 554, 534

* Hawley J. F., Krolik J. H., 2001, ApJ, 548, 348

* Hawley J. F., Balbus S. A., Winters W. F., 1999, ApJ, 518, 394

* Hawley J. F., Guan X., Krolik J. H., 2011, ApJ, 738, 84

* Kato S. ed., 2016, Oscillations of Disks: Astrophysics and Space Science Library, Vol. 437. Springer, Japan, p. 30

* Kigure H., Shibata K., 2005, ApJ, 634, 879
Kiuchi K., Sekiguchi Y., Shibata M., Taniguchi K., 2010, Phys. Rev. Lett., 104, 141101

Kiuchi K., Shibata M., Montero P. J., Font J. A., 2011, Phys. Rev. Lett., 106, 251102

Komissarov S. S., 2006, MNRAS, 368, 993

Korobkin O., Abdikamalov E. B., Schnetter E., Stergioulas N., Zink B., 2011, Phys. Rev. D, 83, 043007

Kumar P., Zhang B., 2015, Phys. Rep., 561, 1

Lesur G., Longaretti P.-Y., 2005, A&A, 444, 25

Loeb A., Ulmer A., 1997, ApJ, 489, 573

Longaretti P.-Y., 2002, ApJ, 576, 587

Lynden-Bell D., 1969, Nature, 223, 690

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

MacFadyen A. I., Woosley S. E., 1999, ApJ, 524, 262

Machida M., Matsumoto R., 2003, ApJ, 585, 429

Machida M., Nakamura K., Matsumoto R., 2004, PASJ, 56, 671

Marconi A., Risaliti G., Gilli R., Hunt L. K., Maiolino R., Salvati M., 2004, MNRAS, 351, 169

McKinney J. C., Blandford R. D., 2009, MNRAS, 394, L126

McKinney J. C., Tchekhovskoy A., Blandford R. D., 2012, MNRAS, 423, 3083

Mewes V., Font J. A., Galeazzi F., Montero P. J., Stergioulas N., 2016, Phys. Rev. D, 93, 064055

Michel F. C., 1972, Ap&SS, 15, 153

Narayan R., Goodman J., 1989, in Meyer F., ed., Proc. of NATO Advanced Science Institutes (ASI) Series C, Vol. 290, Theory of Accretion Disks. Kluwer, Dordrecht, p. 231

Narayan R., Yi I., 1995, ApJ, 452, 710

Nealon R., Price D. J., Bonnerot C., Lodato G., 2018, MNRAS, 474, 1737

Noble S. C., Krolik J. H., Hawley J. F., 2010, ApJ, 711, 959

Papaloizou J. C. B., Pringle J. E., 1984, MNRAS, 208, 721

Piran T., 1999, Phys. Rep., 314, 575

Rees M. J., 1984, ARA&A, 22, 471

Remillard R. A., McClintock J. E., 2006, ARA&A, 44, 49

Reynolds C. S., 2014, Space Sci. Rev., 183, 277

Rezzolla L., Baiotti L., Giacomazzo B., Link D., Font J. A., 2010, Class. Quantum Gravity, 27, 114105

Ryan B. R., Gammie C. F., Fromang S., Kestener P., 2017, ApJ, 840, 6

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

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

Wielgus M., Fragile P. C., Wang Z., Wilson J., 2015, MNRAS, 447, 3593

Woosley S. E., 1993, ApJ, 405, 273

Zanotti O., Rezzolla L., Del Zanna L., Palenzuela C., 2010, A&A, 523, A8

## APPENDIX: COORDINATE TRANSFORMATION

To exploit the regularity property of the Kerr–Schild (KS) coordinates at the event horizon, we perform a transformation of all the vectorial variables computed for a stationary torus in Boyer–Lindquist (BL) coordinates. In the following, the labels KS and BL indicate quantities measured by a KS and BL *Eulerian observer*, respectively, while unprimed and primed indices refer to quantities expressed in KS and BL *coordinates*.

One has to consider the linear transformation:

$$A_{\mu'}^{\mu} = \begin{pmatrix} 1 & G & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & H & 0 & 1 \end{pmatrix} \quad \text{with} \quad \begin{cases} G = -\frac{2r}{\Delta} \\ H = -\frac{a}{\Delta} \end{cases} \qquad (A1)$$

where $\Delta = r^2 - 2Mr + a^2$, r is the radial coordinate, and a is the black hole spin (which in general will be different from 0). This

transformation relates vectors and tensors, respectively:

$$x^{\mu} = A_{\mu'}^{\mu} x^{\mu'}, \qquad (A2)$$

$$T^{\mu\nu} = A_{\mu'}^{\mu} A_{\nu'}^{\nu'} T^{\mu'\nu'}, \qquad (A3)$$

We also have to consider that *Eulerian observers* in the two coordinate systems are not identical in general; that is, we cannot simply apply the transformation in equation (A1) to the vectorial primitive variables (*v*, *B*, *E*), which represent quantities measured in the *Eulerian frame of reference*. For quantities like the fluid 4-velocity $u^\mu$ and the Faraday tensor $F^{\mu\nu}$, on the other hand, one must only apply equation (A1) to fully take into account the change of frame of reference (e.g. $A_{\mu'}^{\mu} u_{\mu'}^{\mu'} = u_{\mu}^{\mu}$).

To obtain for *v*, *E*, and *B* the correct relation between KS components measured by the KS *Eulerian observer* and BL components measured by the BL *Eulerian observer*, we first write the spatial velocity, magnetic field, and electric field in terms of 4-velocity and Faraday tensor in KS coordinates as measured by the KS *Eulerian observer*:

$$v_{\text{KS}}^{i} = \frac{1}{\alpha_{\text{KS}}} \left( \frac{u_{\text{KS}}^{i}}{u_{\text{KS}}^{\prime}} + \beta_{\text{KS}}^{i} \right), \qquad (A4)$$

$$B_{\text{KS}}^{i} = \alpha_{\text{KS}} F_{\text{KS}}^{*ti}, \qquad (A5)$$

$$E_{\text{KS}}^{i} = \alpha_{\text{KS}} F_{\text{KS}}^{ti}, \qquad (A6)$$

where $\alpha$ is the lapse function and $\beta^i$ is the shift vector. Next, we use equation (A1) to transform the components of $u_{\text{KS}}^i$, $F_{\text{KS}}^{*ti}$, and $F_{\text{KS}}^{ti}$ in BL coordinates, then we write these components in terms of $v_{\text{BL}}^{i}$, $E_{\text{BL}}^{i}$, and $B_{\text{BL}}^{i}$. The result gives the correct transformation rules:

$$v_{\text{KS}}^{i} = \frac{1}{\alpha_{\text{KS}}} \left[ \frac{A_{j'}^{i}(\alpha_{\text{BL}} v_{\text{BL}}^{j'} - \beta_{\text{BL}}^{j'})}{1 - A_{r'}^{0} \alpha_{\text{BL}} v_{\text{BL}}^{r'}} + \beta_{\text{KS}}^{i} \right], \qquad (A7)$$

$$B_{\text{KS}}^{i} = \frac{\alpha_{\text{KS}}}{\alpha_{\text{BL}}} \left\{ B_{\text{BL}}^{i} - A_{r'}^{0} \times \left[ B_{\text{BL}}^{r'} \beta_{\text{BL}}^{i} + \frac{\alpha_{\text{BL}}}{\gamma_{\text{BL}}^{1/2}} (\delta_{\phi'}^{i} E_{\theta'}^{\text{BL}} - \delta_{\theta'}^{i} E_{\phi'}^{\text{BL}}) \right] \right\}, \qquad (A8)$$

$$E_{\text{KS}}^{i} = \frac{\alpha_{\text{KS}}}{\alpha_{\text{BL}}} \left\{ E_{\text{BL}}^{i} - A_{r'}^{0} \times \left[ E_{\text{BL}}^{r'} \beta_{\text{BL}}^{i} - \frac{\alpha_{\text{BL}}}{\gamma_{\text{BL}}^{1/2}} (\delta_{\phi'}^{i} B_{\theta'}^{\text{BL}} - \delta_{\theta'}^{i} B_{\phi'}^{\text{BL}}) \right] \right\}. \qquad (A9)$$

This paper has been typeset from a TeX/LA TeX file prepared by the author.