© ESO 2022

1. Introduction

The process of accreting matter towards the central object in an accretion disc requires angular momentum to be transported outwards (Lynden-Bell 1969; Shakura & Sunyaev 1973; Lynden-Bell & Pringle 1974; Balbus & Hawley 1998). Molecular viscosity is insufficient to generate the required angular momentum transport and, hence, turbulence is considered the source of the required stresses (see e.g., Shakura & Sunyaev 1973; Balbus & Hawley 1998). Various candidates for instabilities providing the necessary turbulence have been considered. Two prominent candidates are magneto-rotational instability (MRI; see e.g., Balbus & Hawley 1991, 1998; Hawley et al. 1995; Brandenburg et al. 1995) and gravitational instability (GI; Toomre 1964; Lynden-Bell & Kalnajs 1972; Gammie 2001; Kratter & Lodato 2016).

Magneto-rotational instability is proposed as a mechanism for generating turbulence in sufficiently ionised discs (Balbus & Hawley 1991; Blaes & Balbus 1994; Hawley et al. 1995) including accretion discs of binary objects (black holes, neutron stars, white dwarfs), active galactic nuclei (Balbus & Hawley 1998), and sufficiently ionised regions of protoplanetary discs (Armitage 2011). MRI relies on the coupling of adjacent fluid elements by magnetic field lines producing the necessary mechanism for outward angular momentum transport (see e.g., Balbus & Hawley 1998). The linear and nonlinear properties of MRI are studied in various numerical simulations and for a variety of different configurations (Hawley et al. 1995; Brandenburg et al. 1995; Stone et al. 1996; Suzuki & Inutsuka 2009; Guan & Gammie 2011; Shi et al. 2009; Simon et al. 2011; Bai & Stone 2013; Fromang et al. 2013; Bodo et al. 2014).

The GI plays an important role in galactic discs, contributing to the spiral structure (see e.g., Toomre 1964; Lin & Shu 1964; Lynden-Bell & Kalnajs 1972; Kratter & Lodato 2016), and it is also relevant for active galactic nuclei (Menou & Quataert 2001; Goodman 2003) and for sufficiently massive protoplanetary discs (see e.g., Armitage 2011; Kratter & Lodato 2016). A measure of stability is provided by the Toomre parameter (Toomre 1964):

$$\begin{array}{c}\hfill Q=\frac{c_{\mathrm{s}}{\mathrm{\Omega }}_0}{\pi G\mathrm{\Sigma }},\end{array}$$(1)

with the sound speed c_s, Kepler orbital frequency Ω₀, and mass surface density Σ. The system is gravitationally unstable for Q <  1 (see e.g., Toomre 1964). It is noted that this holds for the thin-disc limit, and three-dimensional effects can alter the stability criterion, typically stabilising the disc (Kratter & Lodato 2016; Binney & Tremaine 1987). Furthermore, turbulent viscosities and radiation physics can have an influence (Lin & Kratter 2016). Simulation studies address gravito-turbulence in both two-dimensional setups (Gammie 2001; Rice et al. 2011; Paardekooper 2012; Young & Clarke 2015; Löhnert et al. 2020) and three-dimensional configurations (see e.g., Rice et al. 2003; Lodato & Rice 2004; Boley et al. 2006; Stamatellos & Whitworth 2008; Cossins et al. 2009; Shi & Chiang 2014; Riols et al. 2017; Riols & Latter 2018b; Booth & Clarke 2019; Hirose & Shi 2019. The instability saturates in a nonlinear state that is a statistical balance between radiative cooling on the one hand and thermal energy production due to shocks generated by gravitational turbulence on the other hand (Gammie 2001; Kratter & Lodato 2016).

In addition to the ability to transport angular momentum, the disc’s ability to sustain large-scale magnetic fields via a dynamo is also attributed to turbulent properties (see e.g., Moffatt 1978; Brandenburg & Donner 1997; Vishniac & Brandenburg 1997; Balbus & Hawley 1998; Rüdiger & Pipin 2000). Magneto-rotational turbulence is often considered as a dynamo mechanism, sustaining large-scale magnetic fields (see e.g., Brandenburg et al. 1995; Balbus & Hawley 1998; Ziegler & Rüdiger 2000; Lesur & Ogilvie 2008; Gressel 2010; Guan & Gammie 2011; Käpylä & Korpi 2011). It was recently suggested that GI can also provide the means to sustain a dynamo (see e.g., Riols & Latter 2019).

The interplay of magneto-rotational and gravitational turbulence has previously been considered in the context of gravito-magneto limit cycles (see Armitage et al. 2001; Zhu et al. 2010; Martin & Lubow 2011; Martin et al. 2012), suggesting that the heating due to GI may enhance the ionisation level and in turn trigger MRI. A few studies also exist on the direct interplay of GI with magnetic fields or MRI (Fromang et al. 2004; Fromang 2005), and more recently Riols & Latter (2018a, 2019) investigated the interaction of GI and MRI more directly in the context of local shearing box simulations with a β-cooling prescription. A possible dynamo-mechanism is provided in Riols & Latter (2019). A global simulation for the interaction of both is provided by Deng et al. (2020).

The present study is intended to address the question of MRI and GI coexistence for ideal magnetohydrodynamic (MHD) cases with zero-net-flux initial conditions and for irradiated discs in the local shearing-box approximation.

This paper is structured as follows. In Sect. 2, the model equations are discussed and important quantities and averages are defined. In Sect. 3, the numerical scheme is described, as well as the boundary conditions, initial conditions and the method used to solve the Poisson equation for self-gravity. Section 4 briefly discusses benchmark simulations for magneto-rotational turbulence as well as gravitational-turbulence. Simulations including both self-gravity and MHD are studied in Sect. 5. The time evolution of field amplitudes, energy densities, and the saturated turbulent stresses are analysed in detail. In Sect. 6, we explain why a superposition of gravitational and magneto-rotational turbulence is likely to be present. The influence of irradiation is addressed in Sect. 7, thereby making a comparison to Riols & Latter (2018a, 2019). The importance of self-gravity in the dynamo mechanism is highlighted in Sect. 8. In Sect. 9, it is then shown that the findings are consistent with an α − Ω model for the dynamo. The main results are summarised in Sect. 11.

2. The model

The model consists of the MHD equations of motion in the local shearing box approximation, retaining the effect of self-gravity as well as a cooling (heating) term:

$$\begin{array}{cc}& {\partial }_t\rho +\mathrm{\nabla }·\left(\rho \mathit{v}\right)=0,\hfill \end{array}$$(2a)

$$\begin{array}{cc}\hfill & {\partial }_t(\rho \mathit{v})+\mathrm{\nabla }·(\rho \mathit{v}\mathit{v}-\frac{1}{{\mathrm{\mu }}_0}\mathbf{B}\mathbf{B}+(P+\frac{{\mathbf{B}}^2}{2{\mathrm{\mu }}_0})\underline{\mathbf{I}}+\underline{\mathbf{G}}),\hfill \\ \hfill & =-2\rho {\mathrm{\Omega }}_0{\mathbf{e}}_z\times \mathit{v}+3\rho {\mathrm{\Omega }}_0^{2}x{\mathbf{e}}_x-\rho {\mathrm{\Omega }}_0^{2}z{\mathbf{e}}_z\hfill \end{array}$$(2b)

$$\begin{array}{cc}\hfill & {\partial }_t\mathbf{B}-\mathrm{\nabla }\times (\mathit{v}\times \mathbf{B})=\eta {\mathrm{\nabla }}^2\mathbf{B},\hfill \end{array}$$(2c)

$$\begin{array}{cc}\hfill & {\partial }_{t}E+\mathrm{\nabla }·((E+P+\frac{{\mathbf{B}}^2}{2{\mathrm{\mu }}_0})\mathit{v}-\mathbf{B}(\mathbf{B}·\mathit{v}))=-\rho \mathit{v}·\mathrm{\nabla }\mathrm{\Phi },\hfill \\ \hfill & +3{\mathrm{\Omega }}_0^{2}x{v}_x-{\mathrm{\Omega }}_0^{2}z{v}_z+\rho \dot{q}.\hfill \end{array}$$(2d)

$$\begin{array}{cc}\hfill & {\mathrm{\nabla }}^2\mathrm{\Phi }=4\pi G\rho ,\hfill \end{array}$$(2e)

$$\begin{array}{cc}\hfill & \mathrm{\nabla }\times \mathbf{B}={\mathrm{\mu }}_0\mathbf{J}.\hfill \end{array}$$(2f)

In the equations above, ρ is the mass-density, v the fluid velocity, P the thermal pressure, B the magnetic field, E the total energy density, J the current density, and $\underline{\mathbf{G}}$ the gravitational stress tensor. The $\underline{\mathbf{G}}$-tensor is given by

$$\begin{array}{c}\hfill \underline{\mathbf{G}}=\frac{1}{4\pi G}(\mathrm{\nabla }\mathrm{\Phi }\mathrm{\nabla }\mathrm{\Phi }-\frac{1}{2}(\mathrm{\nabla }\mathrm{\Phi }·\mathrm{\nabla }\mathrm{\Phi })\underline{\mathbf{I}})\end{array}$$(3)

(see e.g., Lynden-Bell & Kalnajs 1972), with the potential of self-gravity Φ. The first term is a dyadic product and the second term contains the identity-matrix $\underline{\mathbf{I}}$. The total energy density is defined by

$$\begin{array}{c}\hfill E=\frac{\rho {\mathbf{v}}^2}{2}+\frac{P}{\gamma -1}+\frac{{\mathbf{B}}^2}{2{\mathrm{\mu }}_0}.\end{array}$$(4)

The cooling (heating) model $\rho \dot{q}$ includes both pure radiative cooling and a heating source:

$$\begin{array}{c}\hfill \rho \dot{q}=-\frac{P}{(\gamma -1){\tau }_{\mathrm{c}}}+\frac{\rho c_{\mathrm{s},0}^2}{\gamma (\gamma -1){\tau }_{\mathrm{c}}},\end{array}$$(5)

where c_s, 0 is the background sound speed. The first term on the right hand side corresponds to a simple β-cooling prescription with cooling timescale τ_c. The second term is the heating-source due to irradiation (e.g., from the central star). In general, the sound speed is defined here as $c_{\mathrm{s}}=\sqrt{\gamma P/\rho }$, with adiabatic index γ. Unless stated otherwise, γ = 1.64 is used. In the absence of turbulence, the interplay of cooling and heating generates a stationary, stratified state with ρ = ρ(z), P = P(z), and $c_{\mathrm{s}}=c_{\mathrm{s},0}=\sqrt{\gamma P(z)/\rho (z)}=\mathrm{const}$. Hence, irradiation heating would balance radiative cooling such that the disc settles into a state with constant temperature in the vertical direction. It is noted that in turbulent states, one may find ⟨ρ⟩_xy ≠ ρ(z) and ⟨P⟩_xy ≠ P(z), with a horizontal average ⟨…⟩_xy and $\sqrt{\gamma {⟨P⟩}_{x\mathit{y}}/{⟨\rho ⟩}_{x\mathit{y}}}\ne c_{\mathrm{s},0}$ in general.

To investigate turbulent states, averages of different types are used. Spatial averages are denoted by ⟨⟩_i, whereby the identifier i defines the average region; for example, ⟨f⟩_xy averages f over horizontal planes z = const, and the average is thus a function of z only. Without further designation, ⟨f⟩ is an average over the entire box-volume. Time averages are indicated by ⟨f⟩_t. A statement of the form ⟨⟨f⟩⟩_t then means that f is first averaged over the volume and the resulting values are then averaged over time.

To quantify the strength of gravitational instability, the Toomre-parameter is used (see e.g., Toomre 1964; Kratter & Lodato 2016; Gammie 2001):

$$\begin{array}{c}\hfill Q=\frac{c_{\mathrm{s}}{\mathrm{\Omega }}_0}{\pi G\mathrm{\Sigma }},\end{array}$$(6)

where Σ is the averaged mass surface density, Σ = ⟨ρ⟩ L_z, ⟨ρ⟩ is the volume-averaged mass density, and L_z is the vertical extend of the box volume. In the following discussions, two different definitions for the Toomre parameter are used, Q₀ and ⟨Q⟩. The difference arises due to different choices for the sound-speed. For Q₀, the irradiation value is used (c_s = c_s, 0), and for ⟨Q⟩, the volume average is used (c_s = ⟨c_s⟩). Hence,

$$\begin{array}{c}\hfill Q_0=\frac{c_{\mathrm{s},0}{\mathrm{\Omega }}_0}{\pi G\mathrm{\Sigma }},⟨Q⟩=\frac{⟨c_{\mathrm{s}}⟩{\mathrm{\Omega }}_0}{\pi G⟨\mathrm{\Sigma }⟩},⟨c_{\mathrm{s}}⟩=\sqrt{\gamma \frac{⟨P⟩}{⟨\rho ⟩}}·\end{array}$$(7)

Another frequently used quantity is the Alfven speed: $v_{\mathrm{A}}=\sqrt{B^2/({\mathrm{\mu }}_0\rho )}$ (see e.g., Jackson 2014; Balbus & Hawley 1998). When applying a volume average, $⟨v_{\mathrm{A}}⟩=\sqrt{⟨B^2⟩/({\mathrm{\mu }}_0⟨\rho ⟩)}$ is used.

The dimensionless turbulent stress parameter α (see Shakura & Sunyaev 1973) is defined by

$$\begin{array}{c}\hfill \alpha =\frac{2⟨S_{\mathit{xy}}⟩}{3\gamma ⟨P⟩}·\end{array}$$(8)

The additional factor 2/(3γ) is frequently used in the context of gravitational turbulence (see e.g., Gammie 2001) and is also used here. The turbulent stress is comprised of three contributions, the Reynolds stress, the Maxwell stress, and the gravitational stress:

$$\begin{array}{c}\hfill S_{\mathit{xy}}=\rho v_x\delta v_y-\frac{1}{{\mathrm{\mu }}_0}{B}_x{B}_y+\frac{1}{4\pi G}{\partial }_x\mathrm{\Phi }{\partial }_y\mathrm{\Phi }.\end{array}$$(9)

In dimensionless form (Eq. (8)), the corresponding stresses are referred to as α_r, α_m, and α_g, respectively.

In the simulations and plots, all quantities f are made dimensionless by using characteristic scales f_ch, yielding f = f̂ f_ch, with f̂ being dimensionless. Times are made dimensionless using the fiducial orbital angular frequency, that is, $t_{\mathrm{c}\mathrm{h}}={\mathrm{\Omega }}_0^{-1}$. Velocities scale with the background speed of sound v_ch = c_s, 0. Lengths are made dimensionless by using the scale height H = c_s, 0/Ω₀ and l_ch = H. It is noted that the exact definition of H varies in the literature as the additional factors $\sqrt{\gamma }$ or $\sqrt{2}$ may arise and H should rather be understood as a typical scale. Mass densities are made dimensionless by using ${\rho }_{\mathrm{c}\mathrm{h}}={\mathrm{\Omega }}_0^2/G$, and from there a typical pressure of $P_{\mathrm{c}\mathrm{h}}=c_{\mathrm{s},0}^2{\rho }_{\mathrm{c}\mathrm{h}}$ immediately follows. A typical magnetic field is conveniently defined as $B_{\mathrm{ch}}=\sqrt{{\mathrm{\mu }}_0{P}_{\mathrm{ch}}}$. Finally, the characteristic gravitational potential is ${\mathrm{\Phi }}_{\mathrm{c}\mathrm{h}}=c_{\mathrm{s},0}^2$.

3. Method and code

3.1. Code and algorithm

For all simulations shown here, the MHD code Athena ¹ is used (see Stone et al. 2008). The Riemann solver applied here is the Roe solver (Roe 1981), and the spatial reconstruction is third order (Colella & Woodward 1984). The integrator used is the corner transport upwind (CTU) integrator (see Colella 1990). The equations solved are Eqs. (2a)–(2f) (Stone & Gardiner 2010). The fast advection in rotating gaseous objects (FARGO) algorithm is also enabled (see Masset 2000; Stone & Gardiner 2010). The latter separates the velocity into background shear and perturbation v = −(3Ω₀/2)xe_y + δv. The equations are then solved for the perturbation and the shear advection separately, and the shear advection is added to the perturbations afterwards. In boxes spanning many scale heights in radial direction, the shear velocities can become large near the boundaries, and FARGO can therefore reduce the required time step. As velocities are expected to reach supersonic values, especially in self-gravitating cases, shock-capturing H correction is enabled (see Sanders et al. 1998).

The problem generator used builds upon the stratified-shearing-sheet generator that is based on the setup in Hawley & Balbus (1992), Stone et al. (1996). The time step is mediated via the Courant–Friedrichs–Lewy (CFL) number (see e.g., Courant et al. 1928), with a CFL value of 0.06.

3.2. Boundary conditions

All variables (velocities without background shear flow) are periodic in the y direction and shearing-periodic in the x direction. In the vertical direction, outflow boundary conditions are used for hydrodynamical variables, that is, the ghost-zone velocities are set to zero for inward motion and extrapolated constantly into the ghost zones otherwise. Mass density and pressure are extrapolated constantly into the ghost zones, whereby the pressure in the ghost zones is isothermally linked to the density. For the magnetic field, we utilised vertical field boundary conditions at the vertical boundaries. Hence, B_x = B_y = 0 and ∂_zB_z = 0 (see e.g., Brandenburg et al. 1995; Ziegler & Rüdiger 2001; Käpylä & Korpi 2011; Oishi & Low 2011) at the vertical boundaries (B_x, y is set to zero in the ghost zones and B_z is extrapolated constantly into the ghost zones). The mass density can be very small near the vertical boundaries, and in order to limit the time step, a lower boundary for the density is introduced (similarly to Shi & Chiang 2014). The limit is chosen as ρ₀(t = 0) × 10⁻⁴, with ρ₀ being the mid-plane density at initialisation, or ρ₀ = 1.0 for restarted simulations.

3.3. Initial conditions

The background equilibrium consists of a density and pressure distribution (ρ(z),P(z)), whereby the z dependence accounts for the vertical stratification. The background velocity is the Kepler shear flow, v₀ = −(3Ω₀/2)xe_y. Density and pressure are linked by $P=\rho c_{si,0}^2$, assuming a constant isothermal speed of sound, $c_{si,0}=c_{\mathrm{s},0}/\sqrt{\gamma }=\mathrm{const}$. This does not imply that an isothermal equation of state is used, and the condition is only imposed at the start of the simulation. It is noted that c_s, 0 corresponds to the irradiation value in Eq. (5), as this must be the case if the initial state is an equilibrium. With that, one finds ⟨Q⟩=Q₀ for the initial state. The equilibrium is then perturbed by random density perturbations that are introduced at the start of the simulation. Depending on the exact purpose, a seed magnetic field can also be introduced, which we describe in more detail in the respective sections.

The exact form of the stratification profile (ρ(z),P(z)) depends on whether self-gravity is used or not. Without self-gravity, the initialised equilibrium profile for density and pressure is Gaussian ($\rho ={\rho }_0\exp (-{\mathrm{\Omega }}_0^2{z}^2/(2c_{\mathrm{s},i}^2))$). Therein, ρ₀ (P₀ for the pressure) is the mid-plane value at z = 0.

With self-gravity, the stratification profile is given by the vertical force balance (see e.g., Shi & Chiang 2014):

$$\begin{array}{c}\hfill -\frac{1}{\rho }{\partial }_{z}P-{\partial }_z\mathrm{\Phi }-{\partial }_z{\mathrm{\Phi }}^{\star }=0,\end{array}$$(10)

whereby Φ is the gravitational potential due to self-gravity and ${\mathrm{\Phi }}^{\star }=0.5{\mathrm{\Omega }}_0^2{z}^2$ is the vertical contribution from the central star’s potential. Considering that the stratification is initially isothermal, one can rewrite the force balance, yielding

$$\begin{array}{c}\hfill {\partial }_z^2\rho =\frac{1}{\rho }{({\partial }_z\rho )}^2-\frac{4\pi G}{c_{\mathrm{s},i}^2}{\rho }^2-\frac{{\mathrm{\Omega }}_0^2}{c_{\mathrm{s},i}^2}\rho .\end{array}$$(11)

The latter equation is solved once at the start of the simulation using a Runge–Kutta fourth-order method.

3.4. Poisson solver for self-gravity

Athena provides a self-gravity solver based on fast Fourier transforms (FFTs) that can be used in three-dimensional shearing-box calculations. Though, the latter only allows for periodic (shearing periodic in x) boundary conditions in all directions. In order to simulate stratified cases, major changes to the method of solving the Poisson equation are thus required. A possible solution is provided in Koyama & Ostriker (2009), Shi & Chiang (2014). We adopted the code structure provided in Athena and implemented the method outlined in Koyama & Ostriker (2009). Here, a brief analytical explanation of the ansatz is provided, while more details about the numerical implementation are given in Appendix A. One starts by applying a Fourier transform only over (x, y), leaving z unchanged:

$$\begin{array}{cc}& \rho (x,y,z),\mathrm{\Phi }(x,y,z)\to \widehat{\rho }(\mathbf{k},z),\widehat{\mathrm{\Phi }}(\mathbf{k},z)\hfill \\ \hfill & \mathbf{k}={(k_x,k_y)}^T\mathrm{and}k=\sqrt{k_x^2+k_y^2}.\hfill \end{array}$$(12)

The Poisson equation then takes on the form

$$\begin{array}{c}\hfill -k^2\widehat{\mathrm{\Phi }}+{\partial }_z^2\widehat{\mathrm{\Phi }}=4\pi G\widehat{\rho }.\end{array}$$(13)

Solving Eq. (13) for $\widehat{\mathrm{\Phi }}$, assuming a mass density $\widehat{\rho }=\delta (z-z^{\prime })$ yields Green’s function for the vertical direction (see also Koyama & Ostriker 2009):

$$\begin{array}{c}\hfill G(k,(z-z^{\prime }))=2\pi G\{\begin{array}{c}-\frac{1}{k}exp(-k|z-z^{\prime }|)\mathrm{for}k\ne 0\hfill \\ |z-z^{\prime }|\mathrm{for}k=0.\hfill \end{array}\end{array}$$(14)

The solution for arbitrary $\widehat{\rho }(\mathbf{k},z)$ can then be obtained by applying a convolution $\widehat{\mathrm{\Phi }}(\mathbf{k},z)={\displaystyle \int \widehat{\rho }}(\mathbf{k},z^{\prime })G(k,(z-z^{\prime }))\mathrm{d}{z}^{\prime }$. The potential in real space then is

$$\begin{array}{cc}\hfill \mathrm{\Phi }(x,y,z)& =\frac{1}{4{\pi }^2}{\displaystyle \int \left({\int }_{-\infty }^{\infty }\widehat{\rho }(\mathbf{k},z^{\prime })G(k,(z-z^{\prime }))\mathrm{d}{z}^{\prime }\right)}\hfill \\ \hfill & \dots exp(i{k}_{x}x+i{k}_{y}y){\mathrm{d}}^2\mathbf{k}.\hfill \end{array}$$(15)

As is described in more detail in Appendix A, the convolution integral can be evaluated such that only three-dimensional FFTs of sizes N_x × N_y × N_z and N_x × N_y × (2N_z) occur. In order to evaluate the convolution for a finite domain size in z, vacuum boundary conditions (see Koyama & Ostriker 2009; Shi & Chiang 2014) are assumed. This implies that the mass density vanishes outside the vertical boundaries. The values in the vertical boundary cells are calculated by explicitly integrating ∇²Φ = 0 in the z direction. Finally, it is noted that this method is different from the method used in Riols & Latter (2018a). Both agree up to Eq. (13), but differences arise in the treatment of the z dependence. Instead of using Green’s function and convolutions, Riols & Latter (2018a) applied a finite difference scheme in order to solve Eq. (13) for all z values.

4. MRI and SG benchmarks

4.1. MRI turbulence

Pure MHD simulations were prepared as a benchmark for magneto-rotational instability. We solved Eqs. (2a)–(2f) without the self-gravity term $\underline{\mathbf{G}}$ in the Euler-equation and with the energy evolution equation (Eq. (2e)) replaced by an isothermal equation of state $P=\rho c_{\mathrm{s},i}^2$ (c_s, i = 0.780869 is set at simulation start). The initial condition for the magnetic field is a zero-net-flux type (B(t = 0)=B₀sin(2πx/L_x) e_z), with the field pointing in the z direction (similarly to Simon et al. 2011). The initial field amplitude is defined by the ratio of thermal-to-magnetic pressure β₀ at the start of the simulation. Hence, $B_0=\sqrt{2{\mathrm{\mu }}_0{P}_0/{\beta }_0}$, whereby P₀ is the mid-plane value of the initial Gaussian pressure stratification. Instability is initiated by adding small random perturbations to the mass density and the thermal pressure.

Three different configurations varying in resolution and β₀ are chosen (20pHbeta100, 30pHbeta100, and 20pHbeta200). The simulation specifics are listed in the first four lines of Table 1. All configurations share the same box size, whereby 20pHbeta100 and 20pHbeta200 are resolved by ∼20 points per scale height in all three directions and 30pHbeta100 by ∼30 points. The physical box size in the (x, y)-plane is chosen as equal to that used in Simon et al. (2011), Riols & Latter (2018a). The vertical domain size is 8H, which has turned out to reasonably capture the vertical disc structure. The vertical resolution in 20pHbeta100 is chosen similar to that used for GI (see Sect. 4.2) and is also similar to the (x, y)-resolution. It is noted that an additional constraint arises as the total number of processors used has to be divisible by 20 (see Appendix C).

All simulations immediately develop a butterfly diagram, as can be seen for 20pHbeta100 in Fig. 1, showing the horizontally averaged toroidal magnetic-field component ⟨B_y⟩_xy. This is a common feature in zero-net-flux MRI simulations (see e.g., Brandenburg et al. 1995; Brandenburg & Donner 1997; Miller & Stone 2000; Davis et al. 2010; Gressel 2010; Simon et al. 2011).

Fig. 1. (z, t)-diagram of the horizontally averaged toroidal magnetic field component ⟨By⟩xy for simulation 20pHbeta100. A butterfly pattern emerges shortly after the start of the simulation.

We also tested whether the magneto-rotational instability is resolved in the chosen configurations. An approximate criterion for this is given by the ratio of the most unstable MRI wavelength (∼2πv_A, i/Ω₀, see e.g., Balbus & Hawley 1998) to the discretisation length of the grid (δx_i = L_i/N_i). The index i refers to the ith Cartesian component. The resulting quality factor is often referred to as Q_MRI (see e.g., Noble et al. 2010; Simon et al. 2011):

$$\begin{array}{c}\hfill Q_{\mathrm{MRI},i}=\frac{2\pi v_{\mathrm{A},i}}{{\mathrm{\Omega }}_0\delta x_i}=\frac{2\pi N_i}{{\mathrm{\Omega }}_0{L}_i}\sqrt{\frac{⟨B_i^2⟩}{⟨\rho ⟩}}.\end{array}$$(16)

Hence, Q_MRI, i >  1 is required for the MRI to be resolved. Sano et al. (2004) found Q_MRI, z ≳ 6 as a criterion for MRI to be fully resolved. Time-averaged values ⟨Q_MRI, i⟩_t are shown for all three configurations in Table 1. In all cases, one finds ⟨Q_MRI, y⟩_t >  ⟨Q_MRI, x⟩_t >  ⟨Q_MRI, z⟩_t. The lowest value is ⟨Q_MRI, z⟩_t = 7.4, observed for 20pHbeta100 and above the threshold found by Sano et al. (2004). Furthermore, as a clearly visible butterfly diagram emerges, the simulations are considered sufficiently resolved.

In MRI calculations, the turbulent stress has two contributions, the Maxwell part and the Reynolds part, see Eqs. (8) and (9). Time-averaged values of the stress contributions α_r (Reynolds), α_m (Maxwell), and α = α_r + α_m are listed in Table 1. For incompressible and isothermal MRI simulations, the factor (2/(3γ)) in Eq. (8) is often not considered. Hence, for better comparison with other MRI-related literature, that additional factor is not used for the values listed in Table 1. Also, the stresses are sometimes normalised by the mid-plane pressure instead of the volume-averaged pressure. Hence, values normalised by the mid-plane pressure at the start of the simulation P₀(t = 0) are additionally provided in Table 1, designated ${\alpha }_{m,r}^{(0)}$. The obtained values agree well with previous studies. Stone et al. (1996) found (α_m = 0.0044, α_r = 0.00125) for run ‘IZ1’ (compared to values ${\alpha }_{m,r}^{(0)}$ in Table 1) and Simon et al. (2011) found (α_m = 0.022, α_r = 0.0058) for run ‘32Num’ (compared to values α_m, r in Table 1).

In all three configurations, one finds Maxwell-to-Reynolds stress ratios in the α_m/α_r ≈ 3.6−3.8 (Table 1) range. This is typical for magneto-rotational turbulence and is often observed in simulations (see e.g., Hawley et al. 1995; Blackman et al. 2008; Simon et al. 2011; Riols & Latter 2018a).

Another quantity that is found to yield typical values for MRI turbulence is the ratio of Maxwell stress to magnetic pressure: r_sp = ⟨⟨−B_xB_y⟩/⟨B²/2⟩⟩_t (Hawley et al. 1995; Blackman et al. 2008; Simon et al. 2011; Hawley et al. 2011). Observed values often lie in the 0.3−0.4 range. For all three cases studied here, one finds a value of r_sp ∼ 0.35 (see Table 1), which agrees well with the expected range for MRI. As we show below, this also holds for the GI-MRI combined state.

4.2. GI turbulence

In order to study purely gravito-turbulent states in stratified shearing box systems, Eqs. (2a)–(2f) are solved with B = J = 0. The cooling (heating) model used is that in Eq. (5), with the first term representing radiative cooling with cooling time τ_c and the second term accounting for irradiation heating with irradiation temperature ∝c_s, 0. It is noted that in code units, c_s, 0 = 1 (see Sect. 2).

We discuss four different configurations, which are shown as columns in Table 2. The four configurations differ in the total mass contained in the box volume. A convenient dimensionless measure for the mass content is the background Toomre value Q₀, Table 2. Background here means that the constant value c_s, 0 (irradiation) is used in Q₀ instead of c_s, and, hence, Q₀ only depends on the volume-averaged mass density ⟨ρ⟩ (see also Sect. 2). The vertical boundary conditions allow mass to leave the box, causing Q₀ to slightly evolve over time. As the latter effect is small, a time average ⟨Q₀⟩_t over an interval of $120{\mathrm{\Omega }}_0^{-1}$ is provided in Table 2. Listed separately from the background value is the saturated Toomre value ⟨⟨Q⟩⟩_t. The inner average ⟨Q⟩ is over the volume and calculated according to Eq. (7). The obtained values are then averaged over the respective time interval (outer brackets). It is noted that, in contrast to Q₀, ⟨Q⟩ also depends on the saturated sound speed ⟨c_s⟩ and is therefore also sensitive to the saturation level of the thermal energy density.

The box dimensions are the same for all four simulations and are listed in the first three rows of Table 2. Box dimensions of 20H in both horizontal directions turned out to be reasonable for GI, avoiding bursting behaviour due to unstable modes fitting the box size as well as being reasonable in terms of computational effort. The horizontal box size is also similar to that used in Riols & Latter (2018a). The corresponding number of grid points is a compromise between computational accessibility and resolution, whereby the z resolution is chosen to be similar to the horizontal. Hence, the resolution is also comparable to (slightly higher than) that used for 20pHbeta100. Similar to the MRI cases, one has to take into account that the total number of processors used to parallelise over (x, y) is ideally divisible by 20. Hence, 440 points in both horizontal directions is a convenient choice (see also Appendix C). The chosen CFL number in all cases is 0.06. All simulations listed in Table 2 are restarted from an initial simulation that has been brought into the nonlinear state. At the restart time, the mass densities have been rescaled to yield different total box masses (⇒ different Q₀; see Appendix B).

Mass-density snapshots for simulation GI072 are visualised in Fig. 2, with the first image showing a horizontal slice (x, y) and the second image depicting a vertical slice (x, z). In the bottom image, the base-ten logarithm of the mass density is used in order to provide a more convenient visual representation. In the (x, y) slice, one can see density waves and also shocks, separating the waves from lower density regions. Similarly to the findings in Riols & Latter (2018b), vortical motions above and below density waves are also observed. This is captured well in Fig. 3, which depicts an (x, z) slice of the mass density and the vertical stream lines for a fixed time point in the turbulent state of GI058. One can clearly observe four eddies symmetrically oriented around the density maximum near the centre of the image. The stream lines are traced only in the (x, z) plane, that is, neglecting the velocity component v_y. The line colouring is such that motion away from the mid-plane is positive and motion towards the mid-plane is negative (sign(z) v_z).

Fig. 2. Mass-density slices for GI072 in fully developed turbulence, $40{\mathrm{\Omega }}_0^{-1}$ after restart. First image: horizontal slice ρ(x, y, z = 0). Second image: vertical slice log10(ρ)(x, y = 0, z). The logarithm is used in the bottom image as the latter provides a more convenient visual guide for the vertical density profile.

Fig. 3. Poloidal mass-density slice ρ(x, y = 0, z) with stream lines traced in the (x, z) plane. Stream-line colours are mapping the vertical velocity component away from the mid-plane (sign(z) vz). The stream lines are traced for one fixed point in time (t = const.) and by considering only the velocity components (vx, vz).

In hydrodynamic, self-gravitating cases, the turbulent stress has two contributions: a Reynolds-part and a gravitational-part (see Eqs. (8)–(9)). Time averages of the dimensionless stress contributions α_r (Reynolds), α_g (gravitational), and α = α_r + α_g are provided in Table 2. The obtained values for α and α_g are shown as a function of the saturated Toomre parameter ⟨⟨Q⟩⟩_t in Fig. 4, whereby red dots correspond to α_g and black stars correspond to α. The time averages for both α_… and ⟨Q⟩ were taken over an interval of $120{\mathrm{\Omega }}_0^{-1}$ in the fully developed turbulent state of each simulation. The MHD simulations sg-mhd-1 and sg-mhd-2 that will be discussed in later sections are shown as diamonds and crosses. The red dashed line is a guide for the eye (linear fit) for the trend in α_g.

Fig. 4. Time-averaged dimensionless turbulent stresses (Eqs. (8) and (9)) as function of the volume and time-averaged Toomre parameter ⟨⟨Q⟩⟩t. Each data point corresponds to one simulation. The black stars show the sum total of Reynolds and gravitational stress α = αr + αg for the hydrodynamical simulations (GI058, GI072, GI087, GI095). The corresponding contributions from the gravitational stress αg are hown as red dots. Depicted as black diamonds and red crosses are α and αg for the MHD simulations sg-mhd-1 and sg-mhd-2 (discussed in the following sections). In MHD cases, one has α = αr + αg + αm, with Maxwell contribution αm. The red dashed line is a guide to the eye (linear fit) for the trend in αg.

The energy balance is also tested in the nonlinear state. The only way the system can lose or gain energy is by extracting kinetic energy from the background shear flow (Kepler) or by radiative cooling and irradiation heating (Eq. (5)). Minor losses over the vertical boundaries may also occur. In a stationary state, the kinetic energy extracted from the Kepler shear (related to α) is balanced by net cooling and irradiation (related to $⟨\rho \dot{q}⟩$). The balance can be used to derive a relation between the cooling rate and the dimensionless turbulent stress α (see Gammie 2001; Rice et al. 2011). For the cooling-model used here, this is

$$\begin{array}{c}\hfill \alpha =\frac{4}{9\gamma (\gamma -1){\tau }_{\mathrm{c}}{\mathrm{\Omega }}_0}(1-\frac{⟨\rho ⟩c_{\mathrm{s},0}^2}{\gamma (\gamma -1)⟨E_{\mathrm{th}}⟩}),\end{array}$$(17)

with P = (γ − 1)E_th. The values α_p calculated via Eq. (17) are also shown in Table 2. As can be seen, for all simulations the α_p are in good agreement with the values α obtained via Eq. (8), proving that the energy balance between turbulent dissipation and cooling is statistically satisfied with a high level of accuracy.

5. MHD simulations including self-gravity and irradiation heating

5.1. Initial conditions

To investigate interactions of the gravito-turbulent state with magnetic fields, a simulation with fully developed gravito-turbulence but no magnetic field is restarted introducing a magnetic seed field (see also the method used in Riols & Latter 2018a, 2019). The initial field is seeded into the turbulent state of GI072, $40{\mathrm{\Omega }}_0^{-1}$ after the start of the simulation. This corresponds to the density slices shown in Fig. 2. It is noted that, unless stated otherwise, the time axes in all subsequent plots start at tΩ₀ = 120, which is due to the fact that GI072 is restarted from an initial simulation (see Appendix B) and tΩ₀ = 120 corresponds to the absolute time point when the seed field is introduced.

For simplicity, the aforementioned combined simulation is referred to as sg-mhd-1 in the following. The magnetic seed field is chosen such that β ∼ 10¹⁰, to make sure that no artefacts from the initial conditions arise. The initial seed-field configuration is of zero-net-flux type, with the field pointing in the z direction:

$$\begin{array}{cc}& {\mathbf{B}}_0=B_{0}sin(2\pi x/L_x){\mathbf{e}}_z\hfill \end{array}$$(18)

$$\begin{array}{cc}\hfill & B_0={10}^{-5}\leftrightarrow ⟨\beta ⟩=2⟨P⟩/⟨B^2⟩\approx 2.66\times {10}^{10}.\hfill \end{array}$$(19)

All hydrodynamic simulation parameters of GI072 are passed on to sg-mhd-1. This includes the radiative cooling model of Eq. (5), with a cooling-time of τ_cΩ₀ = 10.

5.2. Observed time evolution

This section is used to report on some prominent features that are observed in the self-gravitating MHD simulations. More detailed discussions of some of the aspects are provided in subsequent sections.

The time evolution of the box-averaged magnetic energy density for sg-mhd-1 can be seen in the first image of Fig. 5 as the dashed and dotted lines, corresponding to the left abscissa. Most of the energy is in the y component ($⟨B_y^2⟩$, dotted line), followed by minor contributions from the x component ($⟨B_x^2⟩$, dashed line) and the z component ($⟨B_z^2⟩$, dash-dotted line). The system enters a growth phase immediately after initialisation, with the magnetic energy density increasing between five and eight orders of magnitude within 120 ≲ tΩ₀ ≲ 400. The growth then stagnates for tΩ₀ ≳ 400. In the latter phase, the volume-averaged magnetic field components are oscillating (see e.g., ⟨B_y⟩ as the black solid line in the first image of Fig. 5). The simulation sg-mhd-2 (red solid line) is discussed in more detail in the last paragraph of this section. ⟨B_y⟩ yields significantly larger oscillation amplitudes than the other components, and, during maxima, ∼30% of the total magnetic energy is assigned to ⟨B_y⟩²/2.

Fig. 5. Time evolution for a selection of volume-averaged quantities, whereby the first image shows quantities related to the magnetic field and the second image depicts the turbulent stresses. First image: left axis, dashed and dotted lines: box-averaged magnetic energy densities for Bx, By, and Bz. Right axis, solid lines: box-averaged toroidal field component ⟨By⟩ for simulations sg-mhd-1(black) and sg-mhd-2(red). Simulation sg-mhd-2 is discussed in the last paragraph of this section. Second image: turbulent stresses normalised by pressure sorted for αr (Reynolds), αg (gravity), and αm (Maxwell). The stresses (αr, αg) are smoothed using a Gaussian convolute (standard deviation $\sigma =3{\mathrm{\Omega }}_0^{-1}$).

In order to asses whether MRI can be active, one can calculate the Q_MRI, j parameters (Eq. (16)) for each field component, B_j, using the time evolution of the magnetic energy densities. First considered is the initial phase, 120 ≤ tΩ₀ ≲ 400. Similarly to the energy densities, the Q_MRI, j also increase with time, and it is found that Q_MRI, y ≥ Q_MRI, x ≥ Q_MRI, z. Time averages over three different intervals in the initial phase are shown in Table 3.

As can be seen, Q_MRI, z <  1, except for the last interval. The z contribution is important as the Alfven component v_A, z is connected to the most unstable MRI mode (see e.g., Sano et al. 2004). It was pointed out in Sano et al. (2004) that Q_MRI, z <  6 can lead to the decay of turbulence or significantly lower saturation levels (see also Sect. 1). In Simon et al. (2011), stratified, resistive, and viscous simulations are studied, and it is found that certain parameter regimes lead to time intervals of low turbulent activity. It is also mentioned that those intervals coincide with periods of unresolved poloidal magnetic fields, and turbulence sets back in once the shear has generated enough toroidal field strength B_y to restart turbulence. Hence, one concludes that MRI is, at maximum, marginally resolved in the early phases of sg-mhd-1. Especially in the interval 120 ≤ tΩ₀ ≤ 200, all Q_MRI, j are significantly lower than one, and MRI is almost certainly absent. Further hints for this are also given in Sect. 6. Nevertheless, significant field amplifications are observed in that time period, as can be seen from the energy densities in Fig. 5. However, in the last image of Fig. 6 it is also observed that the ratio

$$\begin{array}{c}\hfill {\langle \frac{F_{x,y}^B}{F_{x,y}^{\mathrm{SG}}}\rangle }_{\mathit{xy}}={\langle \frac{\sqrt{{(\mathbf{B}·\mathrm{\nabla }{B}_x)}^2+{(\mathbf{B}·\mathrm{\nabla }{B}_y)}^2}}{\sqrt{{(\rho {\partial }_x\mathrm{\Phi })}^2+{(\rho {\partial }_y\mathrm{\Phi })}^2}}\rangle }_{\mathit{xy}},\end{array}$$(20)

Fig. 6. (z, t) diagrams for different horizontally averaged quantities in sg-mhd-1. The first two images show the magnetic field components ⟨Bx⟩xy and ⟨By⟩xy, respectively. The third image shows ${⟨B_{\mathit{y}}⟩}_{x\mathit{y}}^{\mathrm{norm}}={⟨B_{\mathit{y}}⟩}_{x\mathit{y}}/\sqrt{⟨B^2⟩}$. The last image shows the ratio of field-line bending forces to gravitational forces in the horizontal direction (see text and Sect. 6 for more details).

relating horizontal forces due to the field line being bent to horizontal forces due to self-gravity, is significantly lower than one for tΩ₀ ≲ 400. Here, the values shown in Fig. 6 are horizontal averages. It is noted that the field-line bending forces are also an important part of MRI. Weak magnetic pressure forces are guaranteed by the high initial thermal-to-magnetic pressure ratio β ∼ 10¹⁰. Hence, one concludes that MRI is most likely absent from the initial phase, and Lorentz forces in general are negligible compared to self-gravity forces. Consequently, one has to assume that the system would be stable without GI and the initial field amplifications must then be caused by kinematic effects of GI turbulence. This suggests that GI provides a possible dynamo mechanism. Additionally, in contrast to the low initial Q_MRI, j values, the saturation level of the Q_MRI, j is roughly two-to-three times larger than for the pure MRI simulations. For the toroidal component, this yields Q_MRI, y ∼ 50−90 for tΩ₀ ≳ 600.

The second image in Fig. 5 shows the turbulent stresses normalised by the volume-averaged pressure, following Eqs. 8 and 9. It is noted that both α_r and α_m significantly fluctuate with time and are therefore smoothed by folding the time series with a Gaussian distribution (standard deviation $\sigma =3{\mathrm{\Omega }}_0^{-1}$).

Similarly to the magnetic energy, the Maxwell contribution to the stress (α_m) is also increasing significantly and even becomes the dominant stress contribution for tΩ₀ ≳ 700. This is in contrast to Riols & Latter (2018a), in which α_g is the dominant contribution.

The volume-averaged toroidal field component ⟨B_y⟩ for sg-mhd-1 is shown in the first image of Fig. 5 as the solid black curve. For tΩ₀ ≳ 400, clear oscillations with growing amplitude can be observed. We only show the y component, though the same oscillation occurs for the ⟨B_x⟩ component as well. This can also be seen in the horizontally averaged field components ⟨B_x⟩_xy and ⟨B_y⟩_xy, which are shown in the first two images of Fig. 6. The oscillations appear as a clearly regular pattern that closely resembles that of a butterfly diagram, similar to MRI turbulence (see e.g., Gressel 2010; Davis et al. 2010; Simon et al. 2011). Field reversals are also found in Riols & Latter (2018a), although they are much less regular and more concentrated to the mid-plane, especially for cooling times tΩ₀ ≲ 20. The pattern obtained here appears to be more akin to that found in the global simulations of Deng et al. (2020).

The first of the bottom two images in Fig. 6 shows ⟨B_y⟩_xy, except that the averaged field components are normalised by the root-mean-squared field, which is averaged over the entire box volume at that particular time point: ${⟨B_{x,\mathit{y}}⟩}_{x\mathit{y}}/\sqrt{⟨B^2⟩}$. From that, one can see that the initial phase tΩ₀ ≲ 400 is qualitatively different from later times. One field reversal near the mid-plane is observed shortly before tΩ₀ = 400, though this may already be a first manifestation of the emerging butterfly diagram. Also, the magnetic activity seems to be concentrated closer to the mid-plane, which is clearly in contrast to the butterfly-diagram that develops later. This initial state appears to be qualitatively closer to the state observed in Riols & Latter (2018a), where most of the activity is near the mid-plane and field reversals occur less regularly.

The last image in Fig. 6 shows the ratio of horizontal field-line bending forces, $\sqrt{{(\mathbf{B}·\mathrm{\nabla }{B}_x)}^2+{(\mathbf{B}·\mathrm{\nabla }{B}_{\mathit{y}})}^2}$, to horizontal self-gravity forces, $\sqrt{{(\rho {\partial }_x\mathrm{\Phi })}^2+{(\rho {\partial }_{\mathit{y}}\mathrm{\Phi })}^2}$, whereby the base-ten logarithm of the values is depicted. For tΩ₀ ≲ 400, self-gravity clearly dominates (see also the discussion above), whereas for later times, forces due to field-line bending can significantly exceed self-gravity above and below the mid-plane (see also Sect. 6).

The total mass in the simulation volume (M = ⟨ρ⟩ L_x L_y L_z) starts to decrease once the magnetic field becomes more dominant tΩ₀ ≳ 400. This trend is shown for the volume-averaged mass density ⟨ρ⟩ by the blue declining curve in Fig. 7. We interpret this as an increased wind activity at the vertical box boundaries. Something similar is also observed for the global simulations in Deng et al. (2020). One can take advantage of this by restarting sg-mhd-1 at a given time point (here tΩ₀ = 600) and holding the total mass in the box volume constant from then on. This restarted simulation is referred to as sg-mhd-2 in the following. Consequently, sg-mhd-1 and sg-mhd-2 can really only be compared for tΩ₀ >  600. Simulation sg-mhd-1 for tΩ₀ >  600 and simulation sg-mhd-2 are exactly equal, with the only difference being that the total mass contained in the box volume is held constant in sg-mhd-2. Hence, the volume-averaged mass density ⟨ρ⟩ stays constant in sg-mhd-2, whereas it slowly decreases in sg-mhd-1 (see the dashed blue lines in Fig. 7). Clearly, at tΩ₀ = 600 the values for ⟨ρ⟩ are exactly equal in sg-mhd-1 and sg-mhd-2. This reflects on the Toomre values ⟨⟨Q⟩⟩_t and the latter split at tΩ₀ = 600 with ⟨Q⟩ increasing for sg-mhd-1 and ⟨Q⟩ staying roughly constant for sg-mhd-2 (see the solid black lines in Fig. 7). This allows one to study the dependence of various quantities on the Toomre parameter in the presence of magnetic fields. For example, one now has two data points for MHD cases in Fig. 4 instead of only one. The last point is elaborated on in more detail in the next section. The replenishing is achieved by calculating the total box mass at each moment and comparing it to the original value (here at tΩ₀ = 600). Once the mass deviates by more than 1%, the latter is restored to its original value by adding a density correction δρ at each grid point, such that the integral yields exactly the missing mass, ∫δρ dV = △M. Thereby, the density correction is distributed homogeneously over (x, y) and has a Gaussian profile in the z direction. Here, $\delta \rho =(▵Mexp(-0.5{(z/h)}^2)/(\sqrt{2\pi }h{L}_x{L}_{\mathit{y}}))$ with h = 0.5 is chosen. Between sg-mhd-1 and sg-mhd-2, no qualitative changes are observed, with the butterfly diagram, for instance, being unchanged. The box-averaged toroidal field component ⟨B_y⟩ for sg-mhd-2 is shown as the red, solid curve in Fig. 5, and it can be seen that no qualitative changes occur, except slight changes in frequency and amplitude. However, quantitative changes in the time-averaged stresses can be observed between sg-mhd-1 and sg-mhd-2. This is discussed in more detail in the next section.

Fig. 7. Saturated Toomre parameter ⟨Q⟩ (solid black, left axis) and box-averaged mass density ⟨ρ⟩ (dashed blue, right axis). At tΩ0 = 600, sg-mhd-1 was restarted, with ⟨ρ⟩ held constant from then on (sg-mhd-2). Therefore, the values of ⟨ρ⟩ start to deviate between sg-mhd-1 and sg-mhd-2 at tΩ0 = 600. Similarly, the values of ⟨Q⟩ also begin to deviate at tΩ0 = 600. For sg-mhd-1, ⟨Q⟩ slightly increases, and for sg-mhd-2, ⟨Q⟩ stays almost constant.

6. GI-MRI coexistence

6.1. Testing the presence of MRI

Table 4 lists important quantities averaged over both volume and time for simulations sg-mhd-1 and sg-mhd-2. For sg-mhd-1, two different time intervals are considered for the average, 120 ≤ tΩ₀ ≤ 320, left column, and 700 ≤ tΩ₀ ≤ 1000, middle column. The first time interval corresponds to the previously discussed initial phase, where the most unstable MRI mode is significantly under-resolved and where MRI is therefore most likely absent. This is further backed up below. Listed in the third column are time averages for sg-mhd-2 (mass held constant) using the second time interval (700 ≤ tΩ₀ ≤ 1000).

In sg-mhd-1, the time-averaged α_m is dominant for 700 ≤ tΩ₀ ≤ 1000, whereas α_g has decreased by roughly 50% compared to the 120 ≤ tΩ₀ ≤ 320 interval. The Reynolds contribution α_r is also reduced by ∼20%, which is significantly less than the reduction of α_g. The total stress (α = α_r + α_g + α_m) increases from the first to the second time interval by ∼14%, which can be attributed to the emergence of α_m for tΩ₀ ≳ 400. Also shown in Table 4 is the predicted total stress, α_p, calculated from the energy balance between turbulent production and radiative cooling (heating) (Eq. (17)). In the low stress interval (120 ≤ tΩ₀ ≤ 320) of sg-mhd-1, both α and α_p are in good agreement, indicating that most of the energy extracted from the background is eventually lost through radiative cooling. In the second interval (700 ≤ tΩ₀ ≤ 1000) of sg-mhd-1, the α_p-value did not change, but the total stress α increased by ∼14%. This agrees with the observed loss of mass, as the generation of a wind will also lead to an additional energy flux that is not accounted for by radiative cooling. Hence, the most dominant effect of introducing magnetic fields is that the gravitational stress α_g is decreasing, whereas the Maxwell stress α_m is increasing and also overcompensating the reductions of both α_g and α_r.

As it turns out, the reduction of α_g in sg-mhd-1 from the first to the second time interval is related to an increase in the saturated Toomre parameter, ⟨Q⟩. The reason for changes in ⟨Q⟩ is twofold. The dominant reason in sg-mhd-1 is that, as discussed in Sect. 4, winds at the vertical boundaries lead to a sinking mass in the box volume, thereby increasing both the irradiation-Toomre Q₀ as well as the saturated value ⟨Q⟩ (see Table 4). The second reason includes changes in the saturated thermal energy density. The irradiation value Q₀ does only react to changes in mass as the irradiation temperature is held constant (c_s, 0 = const.), but the saturated value ⟨Q⟩ can change due to thermal energy changes as it also depends on the saturated sound-speed ⟨c_s⟩. Hence, the ⟨⟨Q⟩⟩_t-difference between sg-mhd-2 and sg-mhd-1 is in part due to the higher mass in sg-mhd-2 (mass held constant), but also due to the slightly different saturation levels of the thermal energy density (see Table 4).

Regardless of the exact changing mechanism, the difference of ⟨⟨Q⟩⟩_t between sg-mhd-1 and sg-mhd-2 allows one to investigate the dependence of α_g on the saturated Toomre-parameter ⟨⟨Q⟩⟩_t for the MHD cases. It turns out that the reduction of the time-averaged α_g from the first to the second time interval is only 33% for sg-mhd-2 compared to 50% for sg-mhd-1 (see Table 4). A possible trend can be directly analysed using Fig. 4. The dimensionless gravitational stresses α_g are shown for sg-mhd-1,2 (700 ≤ tΩ₀ ≤ 1000) via the two red crosses in Fig. 4, whereby the lower α_g value corresponds to sg-mhd-1 (see also Table 4). These values can then be compared to α_g, obtained from the purely hydrodynamical GI-simulations and shown as red dots. One can clearly see a trend for α_g as a function of ⟨⟨Q⟩⟩_t. The data point obtained for sg-mhd-2 is especially close to a hydrodynamical simulation (GI087). It is stressed at this point that sg-mhd-1 (and therefore sg-mhd-2) is restarted initially from GI072 by introducing a seed field. Hence, it appears that the newly formed turbulent states with magnetic fields lead to significantly different values in both ⟨⟨Q⟩⟩_t and α_g; though, the dependence α_g(⟨⟨Q⟩⟩_t) remains similar to the hydrodynamical states. This is clearly not the case for the total stress α = α_r + α_g + α_m. From Fig. 4, one can see that simulations sg-mhd-1,2 (black diamonds) deviate significantly from the trend observed for the purely hydrodynamical cases (black stars). This is of course attributed to the emergence of a strong magnetic stress contribution α_m (see also Table 4). Despite α_g being larger for sg-mhd-2 compared with sg-mhd-1, one can see from Table 4 that the total stress (α = α_r + α_g + α_m) for sg-mhd-2 is exactly equal to that found for sg-mhd-1. Consequently, α_m is found to be lower in sg-mhd-2 compared to sg-mhd-1. What this suggests is that GI with magnetic fields is qualitatively similar to purely hydrodynamic GI; though, the presence of magnetic fields may have an influence on the Toomre-parameter and thereby on the absolute strength of GI, as measured in α_g.

It then remains to be tested whether the magnetic stress-contribution α_m originates from MRI-turbulence. As shown by the previous considerations, the introduction of magnetic fields did not change the turbulence characteristics of GI, although the absolute values ⟨⟨Q⟩⟩_t and α_g can change. Therefore, a simple superposition may be assumed and the turbulent stress is separated into two contributions: α = α(GI)+α(MRI). The magnetic part α_m then originates from MRI and the gravitational stress α_g originates from GI. The Reynolds part has contributions from both MRI and GI. Hence, the following separation of α is suggested:

$$\begin{array}{cc}\hfill \alpha & =\alpha (\mathrm{GI})+\alpha (\mathrm{MRI})\hfill \\ \hfill & ={\alpha }_{\mathrm{g}}+{\alpha }_{\mathrm{r}}(\mathrm{GI})+{\alpha }_{\mathrm{r}}(\mathrm{MRI})+{\alpha }_{\mathrm{m}}.\hfill \end{array}$$(21)

Therein, the Reynolds stress has the contributions α_r = α_r(GI)+α_r(MRI). One then determines the relations of (α_m with α_r) and (α_g with α_r) for pure MRI and pure-GI states, respectively. More precisely, the ratios F_m = α_m/α_r for MRI and F_g = α_g/α_r for GI are defined. It is noted that this assumes the ratios F_g and F_m to be set independently from each other. Hence, the presence of MRI does not affect F_g and GI does not affect F_m. It is again noted that this is supported by the previous considerations, showing that the hydrodynamical trend for α_g can be extended to the MHD cases as well. To see whether the combined states, sg-mhd-1 and 2, are consistent with the independently set values of F_m and F_g, one first determines the ratio F_g from GI simulations and then calculates the resulting F_m. Simulations of magneto-rotational instability suggest typical values of 3.6 ≲ F_m ≲ 4.0 (see Sect. 4.1 and Sect. 8 below). Values for F_g are obtained from the pure GI simulations discussed in Sect. 4 and are listed in the last row of Table 2. One could have started with F_m and derived F_g from that, but F_m is rather a fixed (typical) value, whereas F_g can depend on the mass content measured by Q₀, as can be seen in Table 2. As Table 2 provides F_g for simulations with different total box masses and therefore different values of ⟨Q₀⟩_t, one can estimate the F_g for sg-mhd-1 (reduced mass) and sg-mhd-2 (constant mass) by interpolating between the values in Table 2. For sg-mhd-1 one thereby finds ⟨Q₀⟩_t = 0.86 leading to F_g ≈ 2.0 and for sg-mhd-2 it is ⟨Q₀⟩_t = 0.75 leading to F_g ≈ 1.85.

In the combined state, one has F_g = α_g/(α_r(GI)) and as α_r = α_r(MRI)+α_r(GI), one finds α_r(MRI)=α_r − α_g/F_g. This is then substituted in F_m = α_m/(α_r(MRI)), yielding

$$\begin{array}{c}\hfill F_{\mathrm{m}}=\frac{F_{\mathrm{g}}{\alpha }_{\mathrm{m}}}{F_{\mathrm{g}}{\alpha }_{\mathrm{r}}-{\alpha }_{\mathrm{g}}}.\end{array}$$(22)

Applying Eq. (22) to the values for sg-mhd-1 and sg-mhd-2 (second and third column in Table 4), one finds

$$\begin{array}{c}\hfill F_{\mathrm{m}}=\{\begin{array}{c}3.85\mathrm{For}\mathrm{sg}-\mathrm{mhd}-1\hfill \\ 4.09\mathrm{For}\mathrm{sg}-\mathrm{mhd}-2.\hfill \end{array}\end{array}$$(23)

Hence, the values are in very good agreement with the expectation for MRI turbulence. This also demonstrates that the combined state is consistent with a coexistence of GI and MRI, whereby the ratios F_g and F_m are set independently from each other. For comparison, one can also evaluate F_m for the initial phase 120 ≤ tΩ₀ ≤ 320 of simulation sg-mhd-1 (first column in Table 4). There, one finds ⟨Q₀⟩_t ∼ 0.72, which is expected, as sg-mhd-1 is restarted from GI072 and the initial magnetic fields are weak. From Table 2 then follows F_g = 1.8. Using the time averages in Table 4 and applying Eq. (22), one obtains F_m = 0.223 ≪ 4. This is further proof supporting the argument in Sect. 5.2 that MRI is absent from the initial phase of sg-mhd-1. It also shows that, even when GI is corrected for, the kinematic stress contributions dominate the magnetic contribution and the fields are likely passive in the initial phase. This further strengthens the picture that the initial field amplification is caused by GI turbulence. This is of course significantly different from the saturated state primarily studied in this section.

Similarly to the pure MRI simulations, one can also check for the ratio of Maxwell stress to magnetic pressure r_sp = ⟨⟨−B_xB_y⟩/⟨B²/2⟩⟩_t. For sg-mhd-1, we find r_sp = 0.32, and for sg-mhd-2 r_sp = 0.34, whereby the time averages are evaluated for 700 ≤ tΩ₀ ≤ 1000. Hence, the ratio also agrees well with the ratios obtained for the pure MRI simulations (see Sect. 4.1).

The coexistence hypothesis is further backed up by the last image in Fig. 6, depicting the ratio of horizontal field-line bending forces to horizontal forces of self-gravity ${⟨F_{x,\mathit{y}}^B/F_{x,\mathit{y}}^{\mathrm{SG}}⟩}_{x\mathit{y}}$, evaluated for sg-mhd-1. As can be observed in Fig. 6, values of ${⟨F_{x,\mathit{y}}^B/F_{x,\mathit{y}}^{\mathrm{SG}}⟩}_{x\mathit{y}}>1$ occur for tΩ₀ ≳ 400. The ratio increases with |z|, and for higher altitudes one even finds values up to ∼100 and the bending forces clearly dominate. Only near the mid-plane (|z|≲1H) does one find self-gravity forces to be significantly larger than field-line bending. Since the instability mechanism of MRI relies on field-line bending forces (see e.g., Balbus & Hawley 1998), this is a strong indication that the latter is active, at least for higher altitudes.

6.2. Discussion and comparison to previous studies

At this point, it is useful to present a short discussion about the interpretation of the previously obtained results and also to give a comparison to the findings for the non-irradiated cases studied in Riols & Latter (2018a). First, an overview of similarities and differences to the setup of the SGMRI simulations in Riols & Latter (2018a) is given. The similarities are discussed first. One of them is the horizontal box size, which is 20H in both the x and the y direction. The horizontal resolutions are also similar, whereby Riols & Latter (2018a) used ∼26 points per H, and here, 22 points per H are used. This is true for the vertical resolution as well, which is ∼21 point per H in both cases. The boundary conditions are also similar for both setups. These are outflow-type boundaries for the velocity and vertical-field boundary conditions for the magnetic field. An important difference is the cooling model. In addition to a pure cooling function, the model used here also includes a heating term (see the second term in Eq. (5)). The latter is not used in Riols & Latter (2018a). Additionally, the vertical box size used here (−4H <  z <  4H) is slightly larger than the size used in Riols & Latter (2018a) (−3H <  z <  3H). The comparison provided here especially refers to the SGMRI-(...) simulations in Riols & Latter (2018a), as they also introduced a magnetic seed field into a gravito-turbulent state. It is noted that Riols & Latter (2018a) also provided simulations starting with an MRI state and introducing self-gravity, whereby no significant differences are observed.

Now, differences in the simulation outcomes are elucidated. One difference is the saturation level of the turbulent stresses. Here, it is found that α_m is the dominant contribution in both simulations sg-mhd-1 and 2 and the stresses α_g and α_r are of the same order of magnitude. This is in contrast to SGMRI-10 (τ_cΩ₀ = 10) in Riols & Latter (2018a). There, α_m and α_r are of the same order of magnitude and α_g is the dominant contribution. Moreover, the relative strength of α_m is even weaker for higher cooling times. Differences are also observed in both the time evolution and the vertical structure of the magnetic field. Our simulations sg-mhg-1 and 2 both develop butterfly diagrams with clearly distinguishable field oscillations (see Fig. 6). Furthermore, magnetic activity extends to higher altitudes. This is a distinguishing feature of Riols & Latter (2018a), in which more irregular field reversals were found, and for cooling times τ_cΩ₀ <  100 the magnetic activity is mainly confined to the mid-plane region. It is noted that the concentration of magnetic activity near the mid-plane is also found here, though this is only in the early phases tΩ₀ ≲ 400 where MRI is not resolved or fully developed. We agree with Riols & Latter (2018a) that GI provides a dynamo (see the initial field amplification between $120{\mathrm{\Omega }}_0^{-1}$ and $400{\mathrm{\Omega }}_0^{-1}$). However, it is found that this does not imply suppression of MRI for efficient cooling (τ_cΩ₀ = 10), but rather that the saturated state is characterised by a coexistence of both.

It is noted at this point that the previously discussed simulations do include the effect of irradiation heating (second term in Eq. (5)), which is not the case in Riols & Latter (2018a). Whether the differences are due to the influence of irradiation is elucidated in Sect. 7 below. The boundary conditions used here are comparable to those used in Riols & Latter (2018a); hence, the differences are unlikely to arise from boundary artefacts.

From the previously derived results, one concludes that the turbulent states of simulations sg-mhd-1 and 2 consist of a coexistence between MRI and GI. However, this does not mean that the stresses can simply be added, as is also pointed out in Riols & Latter (2018a). Equation (26) may seem to do that, but it is noted that the saturation value of the Toomre-parameter ⟨Q⟩ also changes with the introduction of magnetic fields. Hence, this is different than adding MRI stresses to the GI stress before magnetic fields have been introduced, as the presence of magnetic fields has an influence on the overall energy balance. Here, the changes in ⟨Q⟩ are in part due to mass changes, but even with constant mass, the additional Maxwell stress −B_xB_y provides additional means to extract energy from the background flow and thereby to change the saturation level of thermal energy.

Finally, it is noted that the results obtained here are more akin to the findings in the global simulations of Deng et al. (2020). There, hints for field-oscillations reminiscent of a butterfly diagram are found as well. Furthermore, the ‘grvmhd1’ simulation in Deng et al. (2020) develops a Maxwell stress in excess of the other stress contributions, which is also a common feature of the simulations presented here.

7. The influence of irradiation-heating

All previous simulations include irradiation heating (second term in Eq. (10)), which is not the case for the simulations in Riols & Latter (2018a). As notable differences arise in the previous chapters, the influence of irradiation will be elucidated in more detail here. For that purpose, restarts are conducted from sg-mhd-1 at tΩ₀ = 920. For the restarts, irradiation is turned off and the total mass held constant. Two different restarts are run, with the first changing the cooling time to τ_cΩ₀ = 20 (sg-mhd-noirrad-20) (a more detailed discussion is provided for that cooling time in Riols & Latter 2018a) and the second keeping the original cooling time τ_cΩ₀ = 10 (sg-mhd-noirrad-10).

The restart of sg-mhd-noirrad-20 is done in two steps. At tΩ₀ = 920, irradiation is turned off, and later, at tΩ₀ = 1060, the cooling time is set to τ_cΩ₀ = 20.

The horizontally averaged toroidal field component ⟨B_y⟩_xy is depicted for sg-mhd-noirrad-20 as functions of z and t in the first image of Fig. 8. As can be seen, the butterfly diagram is at first disrupted but eventually returns. The resulting state is qualitatively different from that found for SGMRI-20 in Riols & Latter (2018a). It is noted that for sg-mhd-noirrad-10, the field oscillations become more irregular near the mid-plane (|z|≲1H), but the butterfly is still visible for higher altitudes. The second image of Fig. 8 depicts the dimensionless turbulent stresses α_r, α_g, and α_m as a function of time. Averaging the stresses over time, starting at tΩ₀ = 1100 yields the values listed in the first line of Table 5. The corresponding values found in Riols & Latter (2018a) (SGMRI-20) are shown in the second line. The same is also shown for sg-mhd-noirrad-10 (third line) and the respective SGMRI-10 simulation of Riols & Latter (2018a) (fourth line).

Fig. 8. Time evolutions of the vertical magnetic field structure and the dimensionless turbulent stresses for simulation sg-mhd-noirrad-20. First image: horizontally averaged field component ⟨By⟩xy as function of (z, t). Second image: dimensionless turbulent stresses (αr, αg, αm). Stresses (αr, αg) were smoothed using a Gaussian convolute (standard deviation $\sigma =3{\mathrm{\Omega }}_0^{-1}$).

The stresses for SGMRI-20 and SGMRI-10 were obtained by applying Eqs. (8) and (9) to the values provided in Riols & Latter (2018a) (see also Appendix D). In addition to the stresses, the ratio of Maxwell stress to magnetic pressure r_sp = ⟨⟨−B_xB_y⟩/⟨B²/2⟩⟩_t is provided in the last column.

As one can see, all stress contributions (α_r, α_g and α_m) in sg-mhd-noirrad-10 are comparable to those obtained in SGMRI-10. Regarding the ratio r_sp, both simulations yield values from 0.3−0.4, which is also typical of MRI. The situation is different for sg-mhd-noirrad-20 and SGMRI-20, with the stresses differing significantly. The Reynolds-contribution α_r is mostly similar in both cases, but the Maxwell and gravitational contributions switch places. In sg-mhd-noirrad-20, the Maxwell part α_m is dominant whereas in SGMRI-20 the gravitational contribution α_g is dominant. Also, the magnetic stress-to-pressure ratio for sg-mhd-noirrad-20 is r_sp ∼ 0.32, which is also MRI typical, but for SGMRI-20 it is found to be r_sp ∼ 0.28, which is considerably lower than all other values. We also analysed the sum total stress α for both cooling times. One first notices that the cases of Riols & Latter (2018a) and our simulations agree in that their α values differ by roughly a factor of two between τ_cΩ₀ = 10 and τ_cΩ₀ = 20. This finding is reasonable, considering that the total stress α is a measure of the rate at which energy is drawn from the Kepler flow and that the cooling timescale is a measure of how fast the system can lose thermal energy. Assuming a steady state, the extracted energy must be balanced by the losses. Reducing the τ_cΩ₀ by a factor of two means doubling the energy-loss rate, which correlates with a doubling of turbulent energy production α. However, considering each cooling time separately, the α values found in Riols & Latter (2018a) and the values found here differ by ∼10%. One likely reason for the differences in α might be different rates of wind cooling. It may also be speculated that this is related to the other differences observed for τ_cΩ₀ = 20 as well. This is also supported by the observed butterfly pattern, as the latter is characterised by magnetic activity rising to higher altitudes. Hence, this is likely an interesting route for future investigations.

One concludes that for τ_cΩ₀ = 10 we can partly reproduce the results obtained in Riols & Latter (2018a), especially that the turbulent stresses are comparable. Although, for τ_cΩ₀ = 20, significant differences arise. Although we currently have no clear explanation for the latter differences, one might speculate that they are related to wind activity at higher altitudes. Furthermore, we come to a different conclusion regarding the coexistence of GI and MRI. Similarly to the irradiated cases, the ratio of Maxwell stress to magnetic pressure is typical of MRI in both cases, and for τ_cΩ₀ = 20 one finds a clearly visible butterfly pattern. Hence, a GI–MRI coexistent state can also be sustained without irradiation, and MRI is not suppressed.

8. Turning off self-gravity

To elaborate on the role of self-gravity, the simulation sg-mhd-1 is restarted at tΩ₀ = 920, with the influence of self-gravity being effectively removed. For that purpose, the gravitational constant is rescaled by a factor of 100 (G = 1 → G = 0.01). That simulation is run until tΩ₀ = 1520 and is referred to as sg-mhd-G001. Time evolutions of the ⟨B_y⟩-field component and the magnetic energy densities are shown in the first image of Fig. 9. The vertical dashed line marks the time point of rescaling, tΩ₀ = 920. The energy densities are mapped to the left abscissa, with curves for $⟨B_x^2⟩/2$ (black, dashed), $⟨B_y^2⟩/2$/2 (blue, dotted), and $⟨B_x^2⟩/2$ (red, dash-dotted). The ⟨B_y⟩-component is mapped to the right abscissa as the black solid curve. Shown in the second image of Fig. 9 is the toroidal field component, averaged horizontally ⟨B_y⟩_xy and as function of (z, t), whereby t starts at $920{\mathrm{\Omega }}_0^{-1}$.

Fig. 9. Time evolution for a selection of volume-averaged quantities, related to the magnetic field, and the evolution of the vertical magnetic field structure, for simulation sg-mhd-G001. First image: volume-averaged magnetic energy densities separated into components $⟨B_i^2⟩/2$ (left abscissa) and volume-averaged toroidal field ⟨By⟩ (right abscissa). Second image: horizontally averaged toroidal field component ⟨By⟩xy as a function of (z, t).

Some observations can immediately be made. The volume-averaged magnetic energy densities decrease significantly. Comparing the time intervals (500 ≤ tΩ₀ ≤ 920) and (1200 ≤ tΩ₀ ≤ 1520), one finds that $⟨⟨B_x^2/2⟩{⟩}_t$ decreases by 92%, $⟨⟨B_y^2/2⟩{⟩}_t$ by 84%, and $⟨⟨B_z^2/2⟩{⟩}_t$ by 93%. The oscillations in the volume-averaged ⟨B_y⟩ component also have significantly lower amplitudes and appear to be less regular without self-gravity. The oscillation frequency stays roughly equal. The butterfly diagram remains, although the mid-plane activity is significantly reduced (see the second image of Fig. 9).

As the magnetic field energy drops by roughly an order of magnitude, one concludes that the GI dynamo dominates in generating the magnetic field. One also concludes that the influence of GI is especially pronounced near the mid-plane, which is expected, as self-gravity has the most influence there. This also agrees with the last image in Fig. 6 (sg-mhd-1), indicating that Lorentz-forces can exceed self-gravity forces, especially for higher altitudes. The potential of GI to act as a dynamo is also consistent with the findings in Riols & Latter (2018a, 2019). The image that emerges is that MRI is operating on the field generated by the GI-dynamo, thereby leading to nonlinear saturation in the butterfly pattern. It is noted that the lack of coherence in the oscillations in the top image of Fig. 9 probably originates from a decoupling of the regions above and below the mid-plane, with the latter two regions in turn oscillating slightly out of phase. This then causes the volume-averaged values to oscillate more irregularly. The less regular oscillations are also present in the smaller box MRI simulations of Sect. 4.1.

Finally, the fact that one obtains a turbulent state for 1200 ≤ tΩ₀ ≤ 1520 is further proof that in the original simulations one does have both GI and MRI. The stresses found in the new turbulent state are α_m = 0.0024 and α_r = 0.00058, yielding a ratio of α_m/α_r ≈ 4.1, which, as expected, agrees with the typical value for MRI.

9. Dynamo properties

Decomposing the velocity (v = (3Ω₀/2)xe_y + δv) into background shear and perturbation δv, one can analytically show that the induction equation (Eq. (2d)) is of the following form:

$$\begin{array}{c}\hfill {\partial }_t\mathbf{B}=\frac{3{\mathrm{\Omega }}_0}{2}x{\partial }_y\mathbf{B}-q{B}_x{\mathbf{e}}_y+\mathrm{\nabla }\times \mathcal{E},\end{array}$$(24)

with electromotive force ℰ = δv × B (see also Vishniac & Brandenburg 1997; Guan & Gammie 2011; Simon et al. 2011). A horizontal average is applied, yielding for the x and y components, respectively (see also Riols & Latter 2019),

$$\begin{array}{cc}& {\partial }_t{⟨B_x⟩}_{\mathit{xy}}=-{\partial }_z{⟨{\mathcal{E}}_y⟩}_{\mathit{xy}},\hfill \end{array}$$(25)

$$\begin{array}{cc}\hfill & {\partial }_t{⟨B_y⟩}_{\mathit{xy}}=-\frac{3{\mathrm{\Omega }}_0}{2}{⟨B_x⟩}_{\mathit{xy}}+{\partial }_z{⟨{\mathcal{E}}_x⟩}_{\mathit{xy}}.\hfill \end{array}$$(26)

It is then assumed that the connection between turbulent electromotive force and averaged field is of the form ${\mathcal{E}}_i={\widehat{\alpha }}_{\mathit{ij}}{⟨B_j⟩}_{x\mathit{y}}$, with coefficients ${\widehat{\alpha }}_{\mathit{ij}}$ and i, j ∈ {x, y}. The z derivatives in the induction equation are estimated to be of the order of ∂_z ∼ 1/H. Substitution into Eqs. (25)-(26) and averaging over z yields (see also Guan & Gammie 2011)

$$\begin{array}{cc}& {\partial }_t⟨B_x⟩=\frac{{\alpha }_{\mathit{xx}}}{H}⟨B_x⟩+\frac{{\alpha }_{\mathit{xy}}}{H}⟨B_y⟩,\hfill \end{array}$$(27a)

$$\begin{array}{cc}\hfill & {\partial }_t⟨B_y⟩=(-\frac{3{\mathrm{\Omega }}_0}{2}+\frac{{\alpha }_{\mathit{yx}}}{H})⟨B_x⟩+\frac{{\alpha }_{\mathit{yy}}}{H}⟨B_y⟩.\hfill \end{array}$$(27b)

It is noted that in the two equations above, the more convenient coefficients ${\alpha }_{\mathit{xx}}=-{\widehat{\alpha }}_{\mathit{y}x}$, ${\alpha }_{x\mathit{y}}=-{\widehat{\alpha }}_{\mathit{yy}}$, ${\alpha }_{\mathit{y}x}={\widehat{\alpha }}_{xx,}$ and ${\alpha }_{\mathit{yy}}={\widehat{\alpha }}_{x\mathit{y}}$ are defined and are also used throughout the paper.

Depending on the exact choice of the α_ij, oscillations in ⟨B_x, y⟩ will occur. As Eqs. (27a)–(27b) are linear in the field components, one can analytically solve the model with an ansatz of the form ⟨B_j⟩=A_jexp(s t), whereby the both s and the A_j are complex numbers. Elementary analytical calculations then yield

$$\begin{array}{cc}\hfill s_{1,2}& =\frac{\mathit{Tr}}{2}\pm \sqrt{{\left(\frac{\mathit{Tr}}{2}\right)}^2-D},\hfill \end{array}$$(28a)

$$\begin{array}{cc}\hfill \mathit{Tr}& ={\alpha }_{\mathit{xx}}+{\alpha }_{\mathit{yy}},\hfill \end{array}$$(28b)

$$\begin{array}{cc}\hfill D& ={\alpha }_{\mathit{xx}}{\alpha }_{\mathit{yy}}-{\alpha }_{\mathit{xy}}({\alpha }_{\mathit{yx}}-1.5{\mathrm{\Omega }}_0).\hfill \end{array}$$(28c)

The second term on the right hand side of Eq. (28a) is important for oscillations. For (Tr/2)² − D <  0, the solutions oscillate with a frequency of $\omega =\sqrt{D-{(Tr/2)}^2}$. In the simpler case of vanishing diagonal elements, α_xx, α_yy, the condition for oscillation reduces to α_xy (α_yx − 1.5Ω₀)< 0. In case oscillations are present, the growth or decay of the latter is set by whether Tr/2 is positive or negative.

All simulations carried out in this study show field oscillations in the form of butterfly diagrams. In the early phases of sg-mhd-1 (tΩ₀ <  400), sign reversals of both ⟨B_x⟩ and ⟨B_y⟩ occasionally occur. This is in contrast to Riols & Latter (2018a), in which irregularly occurring field reversals were found (more akin to the initial phase of sg-mhd-1), and also to Riols & Latter (2019), with field reversals being absent completely.

To test how the simulations relate to the simple dynamo-model of Eqs. (27a)–(27b), correlations between the time derivatives ∂_t⟨B_i⟩ and the components themselves ⟨B_j⟩ were studied. This is visualised in the scatter plots of Fig. 10, with ∂_t⟨B_i⟩ mapped to the y-axis and ⟨B_j⟩ mapped to the x-axis. Two subsequent time points in the scatter plots are separated by $\delta t=0.1{\mathrm{\Omega }}_0^{-1}$, and the time derivatives are explicitly calculated from that. The set of all dots corresponds to the time interval (700 ≤ tΩ₀ ≤ 1000) of sg-mhd-2. Shown in the legends are the slopes and y incidents, which were determined via linear regression. The appearance of ellipsoidal shapes in some plots is an indication that non-zero diagonal elements α_xx and α_yy also contribute. As H = 1 in code units, the slopes can be identified with the coefficients α_ij (see Eqs. (27a)–(27b)). The resulting coefficients are summarised in Table 6.

Fig. 10. Time derivatives of the volume-averaged magnetic field components ∂t⟨Bi⟩ as functions of the field components themselves ⟨Bj⟩. Each dot corresponds to a given time point in the interval (700 ≤ tΩ0 ≤ 1000) of sg-mhd-2, separated by $\delta t=0.1{\mathrm{\Omega }}_0^{-1}$. The partial time derivatives are explicitly calculated from the corresponding time series.

A peculiarity arises for the coefficient α_yx, as determining the latter also requires subtracting the shear (−1.5Ω₀), that is, α_yx = slope + 1.5. Directly reading off the slope from the bottom left image of Fig. 10, one finds the slope ≈ − 1.3. From this, one concludes that the effect of shearing (often called the Ω effect) is the dominant source for ⟨B_y⟩ production.

The α_yx can be calculated more directly by evaluating the term ⟨∂_z⟨ℰ_x⟩_xy⟩_z in Eq. (26) as function of ⟨B_x⟩, thereby avoiding to calculate the time derivative. Generating a similar scatter-plot with this term instead, yields α_yx ≈ 0.197. That is in very good agreement with the value derived by explicitly calculating the time derivative ∂_t⟨B_y⟩, providing a direct check.

The analytical dynamo model in Eqs. (28a)–(28c) is then applied to the values in Table 6, yielding

$$\begin{array}{c}\hfill s=0.000118\pm i·0.0273.\end{array}$$(29)

This implies that oscillations should be present in agreement with the observed butterfly diagram. From ω = 2π/T, one finds a time period $T\approx 230{\mathrm{\Omega }}_0^{-1}$. The actual time evolution of ⟨B_y⟩, evaluated for sg-mhd-2, is shown via the red solid curve in Fig. 5. The dominant frequency obtained by calculating the Fourier series for the interval (700 ≤ tΩ₀ ≤ 1000) yields a period of $T_{\mathrm{s}\mathrm{i}\mathrm{m}}\approx 150{\mathrm{\Omega }}_0^{-1}$. Hence, the measured period is of the same order of magnitude as the period predicted by the model, though they are not exactly equal. The predicted growth rate (Tr/2) is found to be two orders of magnitude lower than the frequency part, which is considered an indication that the simulation is in a quasi-steady state, as amplitude growth stagnates. Doing the same analysis for (400 ≤ tΩ₀ ≤ 700), one finds s ≈ 0.0022 ± i ⋅ 0.031. The frequency is roughly equal in both intervals, though the growth rate is one order of magnitude larger in the earlier time interval. This is expected, as the amplitude growth can clearly be observed for (400 ≤ tΩ₀ ≤ 700) (see Fig. 5).

Riols & Latter (2019) proposed a possible mechanism for the GI-dynamo, which involves horizontal rolls (see Riols & Latter 2018b) above and below the density waves in the mid-plane. These rolls are indeed observed in the simulations presented in this paper (see Fig. 3). Also seen in Fig. 3 is that neighbouring rolls often rotate with opposite orientation. Furthermore, as the rolls are connected to density waves, the rolls are often not exactly parallel to the y-axis. Thus, a magnetic field with only a B_y component can be bent by neighbouring rolls to locally produce a B_x component (see the mechanism proposed in Riols & Latter 2019).

As the horizontal rolls appear slightly above and below the mid-plane (see Fig. 3), they tend to separate the upper atmosphere regions from the mid-plane. In order to demonstrate this, the horizontally averaged enstrophy ${⟨{\omega }_{\mathit{y}}^2⟩}_{x\mathit{y}}^{\mathrm{norm}}={⟨{\omega }_{\mathit{y}}^2⟩}_{x\mathit{y}}/{⟨\delta {\mathit{v}}^2⟩}_{x\mathit{y}}$ is calculated, whereby only the y vorticity ω_y = ∂_zv_x − ∂_xv_z is considered. The enstrophy value is normalised by the mean-squared velocity perturbation ⟨δv²⟩_xy in order to isolate changes in the vortical structure from changes in the magnitude of the velocity. The vertical profile, ${⟨{⟨{\omega }_{\mathit{y}}^2⟩}_{x\mathit{y}}^{\mathrm{norm}}⟩}_t$, obtained by averaging over the time interval (700 ≤ tΩ₀ ≤ 1000), is depicted in Fig. 11, right abscissa. It is noted that the time average is obtained by evaluating ${⟨{\omega }_{\mathit{y}}^2⟩}_{x\mathit{y}}^{\mathrm{norm}}$ for ten equidistant time points in the considered interval. As can be seen from Fig. 11, the enstrophy has pronounced maxima at |z|≈1H, separating the inner regions from regions in the upper atmosphere. Shown in black (left y-axis) are the two contributions to the electromotive force (ℰ_y = v_zB_x − v_xB_z), ⟨⟨v_xB_y⟩_xy⟩_t (dashed) and ⟨⟨v_zB_x⟩_xy⟩_t (dotted). The horizontally averaged magnetic field components periodically change sign, and so do the terms ⟨v_iB_j⟩_xy. For that reason, the ⟨v_iB_j⟩_xy are multiplied with the sign of the volume-averaged field component ⟨B_x⟩ at that particular time point. An average over the time interval (700 ≤ tΩ₀ ≤ 1000) with ten equidistant time points then yields the black curves in Fig. 11. It can be seen that both ⟨⟨v_xB_z⟩_xy⟩_t and ⟨⟨v_zB_x⟩_xy⟩_t change behaviour at |z|≈1H, and on average ∂_z(v_xB_z) and ∂_z(v_zB_x) change sign at |z|≈1H. This is a strong hint that the rolls indeed play an important role in the dynamo mechanism (Riols & Latter 2019).

Fig. 11. Vertical profiles for a selection of quantities (averaged over (x, y), and over time). Left (vertical) axis: Vertical, turbulent shearing (⟨vxBz⟩xy, black dashed line), and vertical compression (⟨vzBx⟩xy, black dotted line) as a function of z. Right (vertical) axis: Toroidal enstrophy normalised by the mean-squared velocity perturbations (${⟨{\omega }_{\mathit{y}}^2⟩}_{x\mathit{y}}^{\mathrm{norm}}={⟨{\omega }_{\mathit{y}}^2⟩}_{x\mathit{y}}/{⟨\delta v^2⟩}_{x\mathit{y}}$, blue solid line) as a function of z.

This is further addressed in Table 7, which shows the α_ij one obtains if the field components are not averaged over the entire box height, but only over given height intervals. Two intervals are chosen, separating the mid-plane region |z|< 1H (interval I₁) and the region above the mid-plane |z|> 1H (interval I₂). This also coincides with the regions below and above the rolls, respectively. The values are calculated by reconstructing the ∂_t⟨B_i⟩ from ∂_z(⟨ℰ_j⟩_xy) and averaging over the corresponding height interval in z. The frequencies can also be calculated for intervals I₁ and I₂, yielding s₁ = 0.000251 ± i ⋅ 0.0505 and s₂ = 0.00535 ± 0.0101, respectively. In terms of time periods, that is $T=2\pi /\omega =124{\mathrm{\Omega }}_0^{-1}$ for I₁ and $T=624{\mathrm{\Omega }}_0^{-1}$ for I₂. The mid-plane value is by far the closest to the observed butterfly frequency of $150{\mathrm{\Omega }}_0^{-1}$, indicating that the frequency is predominantly set near the mid-plane.

The values in Table 7 demonstrate that the dynamo properties do change as a function of z. At higher altitudes, the Ω-effect is less dominant over the turbulent contributions compared to the volume near the mid-plane. α_xy decreases towards higher altitudes. Though, the diagonal elements both increase significantly from the mid-plane to the higher regions. This is in agreement with the situation depicted in Fig. 11 suggesting that the dynamo changes properties at |z|≈1H.

The simulations presented here are hence consistent with the dynamo mechanism based on rolls, as proposed in Riols & Latter (2019); yet, we come to the opposite conclusion in terms of oscillations. In Riols & Latter (2019), the average field grows without sign reversals; although, it has to be noted that Riols & Latter (2019) also included a finite resistivity. The ideal runs in Riols & Latter (2018a) also show field reversals, though they are much less regular than the butterfly pattern found here. In our understanding, the model relying on the rolls connected with the density waves naturally predicts field reversals. Consider a purely toroidal magnetic field line with B_y >  0. The mechanism in Riols & Latter (2019) then suggests that the field line is lifted upwards between two neighbouring rolls. As the rolls are not perfectly aligned with the y-axis, some parts of the field line come to rest on the first roll, and other parts on the second roll. As the rolls rotate in opposite directions, the field line is bent between the two rolls, producing a positive B_x component. This is in agreement with the observation that α_xy is positive. However, the thus produced B_x >  0 is subjected to shear, producing a δB_y <  0. Hence, the newly produced y field has opposite sign to the original field. That will lead to the originally positive B_y being reduced and eventually crossing B_y = 0.

10. Summary

We first investigated the direct interaction of gravitational turbulence with magnetic fields by introducing a magnetic seed field into a gravito-turbulent state (simulation sg-mhd-1, with cooling time τ_cΩ₀ = 10). Significant amplifications in the magnetic energy density and the volume-averaged fields can be observed in the initial phase, shortly after introducing a magnetic seed field. After $\sim 300{\mathrm{\Omega }}_0^{-1}$, the system changes state and a butterfly diagram develops that also extends towards higher altitudes in the z direction. This is in contrast to earlier studies of Riols & Latter (2018a) finding field reversals to occur much less regularly and magnetic activity to be more localised in the mid-plane region. It is noted that a concentration near the mid-plane is found in simulation sg-mhd-1 as well; though, this only occurs in the aforementioned initial phase after introducing the seed-field.

In accordance with the growing magnetic field strength, a Maxwell contribution α_m to the total stress also emerges. As α_m grows, the gravitational stress α_g is found to decrease, whereas the sum total of all stresses slightly increases. At saturation, α_m is the dominant stress contribution, or at least it is of the same order of magnitude as the remaining stress contributions. The reduction in α_g is found to be consistent with an increase in the saturated Toomre parameter ⟨Q⟩. In simulation sg-mhd-1, this can mostly be attributed to a reduction of the total mass in the box volume due to vertical wind activity initiated by the magnetic fields. A second simulation (sg-mhd-2) is restarted from the early turbulent state of sg-mhd-1, whereby the sum total mass is then held constant. Hence, the mass is slightly higher in sg-mhd-2, leading to slightly lower ⟨Q⟩ values, and α_g indeed saturates at a slightly higher level in sg-mhd-2 compared to sg-mhd-1. The total stress α is equal in both cases (sg-mhd-1 and 2), consequently, the saturation value α_m is lower in sg-mhd-2 compared to sg-mhd-1. In addition to that, a method of separating possible GI and MRI stress contributions (α = α(GI)+α(MRI)) is suggested, and it is found that α_m/α_r(MRI)=3.85,  4.09 for sg-mhd-1 and 2, respectively. For the ratio of Maxwell stress to magnetic pressure, one obtains values of ⟨⟨−B_xB_y⟩/⟨B²/2⟩⟩_t = 0.32,  0.34 for sg-mhd-1 and 2, respectively.

The influence of irradiation is demonstrated by using two restarts of sg-mhd-1, with irradiation being turned off at the restart point. One restart also changes the cooling time from τ_cΩ₀ = 10 to τ_cΩ₀ = 20. Both cases still develop butterfly diagrams, whereby the pattern is more pronounced for τ_cΩ₀ = 20. For τ_cΩ₀ = 20 especially, the clearly visible butterfly pattern is different from the findings for SGMRI-20 in Riols & Latter (2018a). The ratio of Maxwell stress to magnetic pressure ⟨⟨−B_xB_y⟩/⟨B²/2⟩⟩_t is found to yield values similar to those obtained with irradiation. The averaged stresses for τ_cΩ₀ = 10 are closest to the values found in Riols & Latter (2018a) (SGMRI-10), though for τ_cΩ₀ = 20 significant differences remain. Most notably, the Maxwell contribution to the dimensionless stress α_m is dominant for the cases presented here, whereas in (SGMRI-20) of Riols & Latter (2018a), the gravitational contribution is dominant α_g.

The role of self-gravity is elucidated by reducing the gravitational constant G in Poisson’s equation to one hundredth of its original value. With self-gravity being effectively removed, the butterfly diagram remains, and we found α_m/α_r ≈ 4.1. The oscillation amplitude and the magnetic energy in the box volume are significantly reduced (one order of magnitude) once self-gravity is weakened. Significantly larger field strengths are observed with full self-gravity being enabled, especially near the mid-plane.

Finally, the saturated state is shown to be consistent with an α − Ω dynamo of the from shown, for example, in Gressel (2010), Guan & Gammie (2011), and Simon et al. (2011). The measured dynamo coefficients (α_xx, α_xy, α_yx, α_yy) agree with the field oscillations apparent in the butterfly diagrams. The order of magnitude of the expected oscillation period, as derived using the simple α − Ω model, is in agreement with the time period observed in the simulations. Taking the intervals |z|> 1H and |z|< 1H separately, it is found that the dynamo coefficients differ in both intervals. This is backed up by the observation that the toroidal vorticity ${⟨{\omega }_{\mathit{y}}^2⟩}_{x\mathit{y}}^{\mathrm{norm}}$ takes on a maximum at |z|≈1H. The vorticity-maximum appears to be connected with rolls aligned with the density waves; the importance of these rolls is also pointed out in Riols & Latter (2018b, 2019).

11. Conclusion

Here, the conclusions that can be drawn from the aforementioned findings are summarised. Essentially, there are two important points; gravitational turbulence provides the means to act as a dynamo, and the saturated state of the GI-dynamo is a coexistence between gravitational turbulence and magneto-rotational turbulence. The first point regarding the GI-dynamo agrees with the findings in Riols & Latter (2018a, 2019), whereas differences arise regarding the coexistence of GI and MRI, which does not seem to be the case for the corresponding simulations in Riols & Latter (2018a) (zero-net-flux, with τ_cΩ₀ ∼ 10−20).

The possibility of a GI-dynamo is most prominently demonstrated in the initial phases, shortly after introducing the magnetic seed field. Significant field amplifications are observed, whereas it is shown that MRI is likely not present as the weak initial field renders it highly unresolved. This is further backed up by the fact that field line bending forces are significantly weaker than self-gravity in the early phases. Hence, the field amplifications can only originate from the kinematic effects of GI turbulence. The GI dynamo is also supported by the significant reduction of magnetic energy once self-gravity is effectively removed. As stated in Sect. 9, vortical motions in the (x, z) plane (rolls) seem to play an important role in the dynamo characteristics, which is in agreement with the findings of Riols & Latter (2018b, 2019). Although, in contrast to Riols & Latter (2019), field oscillations and a butterfly pattern are found that might be linked to the pattern observed in the global simulations of Deng et al. (2020). We thus point out that the roll dynamo must not be in opposition to the field oscillations.

Regarding the second point, there are strong hints that the saturated state is comprised of a coexistence between GI and MRI. One direct hint is the previously mentioned oscillating pattern, which is reminiscent of a butterfly diagram for pure MRI. It is also found that the turbulent stresses (α_r, α_m, α_g) in the nonlinear state agree with a superposition of GI and MRI turbulence. More precisely, one can self-consistently decompose the stresses into contributions of GI and MRI. Thereby, the ratios α_m/α_r and α_g/α_r are assumed to be independent of each other and are only determined by pure MRI and pure GI states, respectively. The presence of MRI is also supported by the ratio of Maxwell stress to magnetic pressure ⟨⟨−B_xB_y⟩/⟨B²/2⟩⟩_t, it being in the MRI typical interval of 0.3−0.4. Hence, in contrast to Riols & Latter (2018a), we come to the conclusion that the turbulent state is consistent with GI–MRI coexistence. Though, care must be taken, as one can not simply add the MRI stresses to the GI stresses, as these are evaluated before the magnetic field is introduced. This is due to the fact that the growing magnetic field strength reduces the saturated Toomre parameter ⟨Q⟩ and, therefore, also the turbulent stresses associated with GI. Cases without irradiation heating are found to be consistent with GI–MRI coexistence as well, and a butterfly diagram also emerges. Riols & Latter (2018a) did not include irradiation, though, and significant differences remain, especially for τ_cΩ₀ = 20. At this point, we have no clear explanation for the discrepancies, although one might speculate that they are related to possible wind activity at the vertical boundaries. Future work is to be done to clarify how MRI might cause the GI dynamo to saturate. Moreover, the influence of non-ideal effects and winds is to be elucidated in more detail.

References

Appendix A: Notes on the Poisson solver

Here, the implementation of the Poisson solver described in Sect. 3.4 is outlined in more detail. The method closely follows that suggested in Koyama & Ostriker (2009). In the following, (x_n1, y_n2, z_n3) with n₁, n₂, n₃ ∈ ℕ refers to a position in the real space grid, and (k_mn = (2πm/L_x, 2πn/L_y)) with m, n ∈ ℕ refers to a point in the two-dimensional Fourier grid. Wave numbers with m, n >  N_x, y/2 are interpreted as negative wave numbers (m − N_x, n − N_y). The absolute value of the wave-vector is k_mn (see also Eq. 12). Green’s function (Eq. 14) is here abbreviated as G_mn(z_n3 − z_j)=G(k_mn, z_n3 − z_j). The expression for the potential (Eq. 15) then reads

$$\begin{array}{cc}& \mathrm{\Phi }(x_{n_1},y_{n_2},z_{n_3})=\frac{\delta z}{N_x{N}_y}{\displaystyle \sum _{m=0}^{N_x-1}\sum _{n=0}^{N_y-1}\sum _{j=0}^{N_z-1}}...\hfill \\ \hfill & ...exp\left(2\pi i[\frac{n_{1}m}{N_x}+\frac{n_{2}n}{N_y}]\right){\rho }_{\mathit{mn}}(z_j)G_{\mathit{mn}}(z_{n_3}-z_j),\hfill \end{array}$$(A.1)

whereby δz = L_z/N_z is the vertical discretisation length and

$$\begin{array}{c}\hfill {\rho }_{\mathit{mn}}={\displaystyle \sum _{n_1=0}^{N_x-1}\sum _{n_2=0}^{N_y-1}}\rho (x_{n_1},y_{n_2},z_{n_3})exp(-2\pi i[\frac{n_{1}m}{N_x}+\frac{n_{2}n}{N_y}])\end{array}$$(A.2)

is the two-dimensional Fourier transform of the mass density. It is also noted that the mass density is assumed to vanish outside the vertical boundaries. Then, G_mn(z_n3 − z_j) is expanded as a Fourier series, whereby z_n3 − z_j ∈ [ − L_z, L_z] (see also Koyama & Ostriker 2009)

$$\begin{array}{cc}\hfill G_{\mathit{mn}}(z_{n_3}-z_j)& =\frac{1}{2N_z}{\displaystyle \sum _{l=0}^{2N_z-1}}exp\left(2\pi i\frac{(n_3-j)l}{2N_z}\right)G_{\mathit{mnl}}\hfill \\ \hfill G_{\mathit{mnl}}& ={\displaystyle \sum _{q=-N_z}^{N_z-1}}exp(-2\pi i\frac{\mathit{ql}}{2N_z})G_{\mathit{mn}}(q).\hfill \end{array}$$(A.3)

Substituting G_mn(z_n3 − z_j) from Eq. 14 and evaluating the sum yields

$$\begin{array}{c}\hfill G_{\mathit{mnl}}=4\pi G\{\begin{array}{c}\frac{1-e^{i\pi l}{e}^{-k_{\mathit{mn}}{L}_z}}{2k_{\mathit{mn}}}\frac{\sinh (k_{\mathit{mn}}\delta z)}{cos(\frac{\pi l}{N_z})-\cosh (k_{\mathit{mn}}\delta z)}\mathrm{for}{k}_{\mathit{mn}}\ne 0\hfill \\ \frac{\delta z}{2}\frac{1+e^{i\pi l}}{cos(k_{\mathit{mn}}\delta z)-1}\mathrm{for}{k}_{\mathit{mn}}=0\hfill \end{array}.\end{array}$$(A.4)

It is noted that this expression for G_mnl is slightly different from that given in Koyama & Ostriker (2009). Though, by expanding the hyperbolic functions and the cosine function for small k_mnδz and lπ/N_z , respectively, one recovers the result given in Koyama & Ostriker (2009).

Then, substituting the expression obtained in Eq. A.4 into Eq. A.1 yields

$$\begin{array}{cc}& \mathrm{\Phi }(x_{n_1},y_{n_2},z_{n_3})=\frac{\delta z}{N_x{N}_y(2N_z)}{\displaystyle \sum _{m=0}^{N_x-1}\sum _{n=0}^{N_y-1}\sum _{j=0}^{N_z-1}\sum _{l=0}^{2N_z-1}}T,\hfill \\ \hfill & T=G_{\mathit{mnl}}{\rho }_{\mathit{mn}}(z_j)exp\left(2\pi i[\frac{n_{1}m}{N_x}+\frac{n_{2}n}{N_y}+\frac{(n_3-j)l}{2N_z}]\right).\hfill \end{array}$$(A.5)

The j summation arises from the convolution and the l summation arises from the substitution of G_mnl. The j summation over ρ_mn(z_j)⋅exp(−2π i jl/(2N_z)) can then be carried out. The summation result is abbreviated by ${\rho }_{\mathit{mnl}}^{\mathit{corr}}$. The expression for the potential then becomes

$$\begin{array}{cc}\hfill \mathrm{\Phi }(x_{n_1},y_{n_2},z_{n_3})& =\frac{\delta z}{N_x{N}_y(2N_z)}{\displaystyle \sum _{m=0}^{N_x-1}\sum _{n=0}^{N_y-1}\sum _{l=0}^{2N_z-1}}{G}_{\mathit{mnl}}{\rho }_{\mathit{mnl}}^{\mathit{corr}}\hfill \\ \hfill & exp\left(2\pi i[\frac{n_{1}m}{N_x}+\frac{n_{2}n}{N_y}+\frac{n_{3}l}{2N_z}]\right),\hfill \end{array}$$(A.6)

which represents a three-dimensional, inverse Fourier transform of size (N_x × N_y × (2N_z)). The inverse FFT yields 2N_z values in the z direction, but the index n₃ is only defined for N_z values. Hence, the potential is extracted by keeping the first N_z indices of the inverse FFT array.

Analytically evaluating the j-summation in ${\rho }_{mnl,}^{\mathit{corr}}$ one can show that ${\rho }_{\mathit{mnl}}^{\mathit{corr}}$ can be calculated using two forward-Fourier transforms of size N_x × N_y × N_z:

$$\begin{array}{c}\hfill {\rho }_{\mathit{mnl}}^{\mathit{corr}}=\{\begin{array}{cc}{\mathcal{FT}}_3(\rho (x_{n_1},y_{n_2},z_{n_3}))(m,n,q)\hfill & \mathrm{for}l\mathrm{even}\\ {\mathcal{FT}}_3(\rho (x_{n_1},y_{n_2},z_{n_3})exp(-i\pi \frac{n_3}{N_z}))(m,n,q)\hfill & \mathrm{for}l\mathrm{odd},\end{array}\end{array}$$(A.7)

whereby the Fourier coordinate q relates to l as l = 2q for even l and l = 2q + 1 for odd l. It is noted that ${\rho }_{\mathit{mnl}}^{\mathit{corr}}$ may also be defined differently by calculating two-dimensional FFTs for each z_j instead of two three-dimensional FFTs. To use the latter method, G_mnl also needs to be modified (see e.g., Koyama & Ostriker 2009). We used three-dimensional FFTs as outlined above, and it is pointed out that our method may deviate from that used in Koyama & Ostriker (2009) at that point.

Appendix B: GI benchmark

As noted in Sect. 4.2, the simulations listed in Table 2 are restarted from an initial simulation that has been brought into the nonlinear state. The initial simulation (initGI18) is started with adiabatic index γ = 1.8, cooling-time τ_cΩ₀ = 10, and initial Toomre-parameter Q₀(t = 0)=0.6 = ⟨Q⟩. A CFL number of 0.06 is used. The choice γ = 1.8 is only applied to initGI18 and has turned out to be necessary, as the linear-nonlinear transition is rather disruptive and directly choosing γ = 1.64 has always lead to a clumping state. Occasionally, the simulation can also halt at the transition point. Choosing a time step that is too large (CFL >  0.06) turns out to also produce clumping states after entering the nonlinear phase. The linear-nonlinear transition occurs roughly $40{\mathrm{\Omega }}_0^{-1}$ after the start of the simulation. After a total of $80{\mathrm{\Omega }}_0^{-1}$, initGI18 is restarted with γ = 1.64 (GI072) and run for $160{\mathrm{\Omega }}_0^{-1}$. Simulations GI058, GI087, and GI095 are restarted from GI072 at tΩ₀ = 120 ($40{\mathrm{\Omega }}_0^{-1}$ after initial restart from initGI18) by rescaling the density on every grid point with a fixed factor to yield different box masses (⇒ different Q₀). However, the exact value of Q₀ may later depend on time due to losses over the vertical boundaries.