On the role of numerical diffusivity in MHD simulations of global accretion disc dynamos

Abstract

Observations, mainly of outbursts in dwarf novae, imply that the anomalous viscosity in highly ionized accretion discs is magnetic in origin and requires that the plasma ${\beta \sim 1}$. Until now, most simulations of the magnetic dynamo in accretion discs have used a local approximation (known as the shearing box). While these simulations demonstrate the possibility of a self-sustaining dynamo, the magnetic activity generated in these models saturates at $\beta \gg 1$. This long-standing discrepancy has previously been attributed to the local approximation itself. There have been recent attempts at simulating magnetic activity in global accretion discs with parameters relevant to the dwarf novae. These too find values of $\beta \gg 1$. We speculate that the tension between these models and the observations may be caused by numerical magnetic diffusivity. As a pedagogical example, we present exact time-dependent solutions for the evolution of weak magnetic fields in an incompressible fluid subject to linear shear and magnetic diffusivity. We find that the maximum factor by which the initial magnetic energy can be increased depends on the magnetic Reynolds number as ${\mathcal {R}}_{m}^{2/3}$. We estimate that current global numerical simulations of dwarf nova discs have numerical magnetic Reynolds numbers around six orders of magnitude less than the physical value found in dwarf nova discs of ${\mathcal {R}}_{m} \sim 10^{10}$. We suggest that, given the current limitations on computing power, expecting to be able to compute realistic dynamo action in observable accretion discs using numerical MHD is, for the time being, a step too far.

1. Introduction

Fluid-based magnetic dynamos can be thought of as coming in two flavours: small scale and large scale. Both types can be found in different astrophysical contexts (see the reviews by Brandenburg & Subramanian 2005; Rincon 2019; Schekochihin 2022). Small-scale dynamos tend to produce magnetic fields that are correlated on the length scale, $l_0$, (or smaller) of the driving turbulence, such as is seen in models of the internal dynamics of molecular clouds (e.g. Padoan et al. 2014; Federrath 2016) or in cosmological simulations (e.g. Martin-Alvarez et al. 2021. Large-scale dynamos typically comprise small-scale turbulence (scale $l _0$) set within a large scale shear flow (scale $L \gg l_0$) and show large-scale spatial coherence. Obvious examples of these are the dynamo responsible for the solar cycle and the dynamo thought to be responsible for driving the accretion flow within accretion discs. In the former, the large-scale flow is the differential rotation of the sun and the turbulence is driven by convective motions. In the latter, the large-scale flow is the Keplerian differential flow of the disc and the small-scale turbulence is due to the magneto-rotational instability (MRI; Brandenburg et al. 1995; Balbus & Hawley 1998), perhaps enhanced by magnetic buoyancy (Tout & Pringle 1992; Johansen & Levin 2008).

Due to the complexity of the resulting dynamics, particularly in the nonlinear, turbulent phase, it is routine to resort to numerical simulations to gain an understanding of the flows generated from such dynamo action. For many years, it was not straightforward to resolve the length scales associated with the growth of the MRI in global accretion disc simulations, and it is therefore necessary to employ a local, or ‘shearing box’, approach to study the dynamics in a simplified framework with less computational power required to achieve high resolution. However, the shearing box approach does not reach the level of dynamo activity that is required to explain the observations. King, Pringle & Livio (2007) note that the shearing box approach necessarily restricts the available modes that can be produced in the flow and, in particular, low-$m$ modes (with $m \ne 0$) are not captured (see also, for example, the discussion by Parkin & Bicknell 2013). There also remain questions regarding the convergence of such models when they are extended to very high resolution (Bodo et al. 2014). We do not pursue such arguments here.

The discrepancy in the strength of the dynamo activity in shearing box models and observed accretion discs motivates the development of global MHD models of accretion discs. However, the computational demands of global models mean that the spatial resolution, typically measured as a fraction of the local disc scale height or the wavelengths of the fastest growing modes, may not be sufficient. In particular, insufficient resolution may imply strong levels of numerical magnetic diffusivity compared with the physical diffusivity expected in the simulated systems. If the numerical diffusivity is far greater than the expected physical diffusivity, then it seems reasonable to expect that this will have implications for the types (strengths and geometries) of fields that are produced in the simulations compared with those produced in real systems (cf. Tobias & Cattaneo 2013).

In § 2, we briefly describe the relevant background to accretion disc physics and the relevant observed properties of such discs. In § 3, we present what is known about the properties of the disc dynamo as can be deduced from observations of the discs in dwarf novae (a particularly well-studied subclass of cataclysmic variable stars in which a white dwarf accretes from a companion star). The discs in the outburst state of these objects are those for which we have the best understanding of the properties of the so-called anomalous viscosity. In these objects, if this viscosity is magnetohydrodynamic in origin, then the observations imply that the magnetic fields are strong (plasma $\beta \sim 1$). In § 4, we discuss the numerical simulations of the MHD properties of accretion discs. We note that in the usual shearing box approximation, it has long been known that the predicted magnetic field strengths are too low ($\beta \gg 1$) to satisfy the observations. It is possible that here the problem lies with the shearing box approximation itself, although we must remark that there is no evidence that this is the case. We therefore consider some recent numerical simulations of global accretion discs, but we note that in these too, the field strengths that are found are similarly small. In § 5, we speculate that for the global disc simulations, the numerical magnetic diffusivity may be too large to allow the required growth of the field. To illustrate this, in § 6, we present calculations of the evolution of magnetic fields embedded in a fluid subject to an inexorable shear. We vary the level of diffusivity to illuminate the effect of it on the maximum magnetic field growth that is achievable. We highlight the limitations if excessive diffusion is included either explicitly or through numerical effects. Finally, we discuss our results in the context of some of the global simulations presented in the literature and draw conclusions in § 7.

2. Accretion discs

The theory of accretion discs is set out by Shakura & Sunyaev (1973, see also Pringle 1981; Frank, King & Raine 2002). The basic disc flow, in cylindrical polar $(R, \phi, z)$ coordinates, is in the azimuthal direction,

(2.1)\begin{equation} u_\phi = \sqrt{GM/R}, \end{equation}

where $M$ is the central gravitating mass. The disc thickness or scale height, $H$, in the $z$-direction is given by

(2.2)\begin{equation} H/R \approx c_s/u_\phi, \end{equation}

where $c_s$ is the local sound speed. For the usual thin disc approximation, we require that $H/R \ll 1$. Thus, the azimuthal motion is supersonic, with Mach number $\sim R/H$. Inflow (i.e. accretion) through the disc requires the action of a so-called ‘anomalous viscosity’ which taps the azimuthal ($R \phi$) shear and transfers angular momentum outwards and mass inwards. The viscosity is generally assumed to be caused by small-scale, $l_0 \le H$, magnetohydrodynamic turbulence within the disc, and it is the maintenance of this turbulence that requires dynamo action. Shakura & Sunyaev (1973) introduced a parameter $\alpha$ which is a dimensionless measure of the size of the anomalous viscosity – essentially a dimensionless measure of the $z$-averaged $R\phi$-stress. Thus, the ($z$-averaged) effective kinematic viscosity of the disc may be written as

(2.3)\begin{equation} \nu \approx \alpha c_s H \approx \alpha H^2 \varOmega, \end{equation}

where

(2.4)\begin{equation} \varOmega = u_\phi/R = \sqrt{GM/R^3} \end{equation}

is the angular velocity at radius $R$. For MHD turbulence, they note that

(2.5)\begin{equation} \alpha \approx \langle B_R B_\phi \rangle / \rho c_s^2 \end{equation}

(in appropriate units), where $\rho$ is the disc density. They also introduced physical arguments as to why we might expect $\alpha \le 1$ (see also the discussion by Martin et al. 2019). In terms of $\alpha$, the radial accretion flow velocity is

(2.6)\begin{equation} u_R \sim{-} \alpha (H/R) c_s \end{equation}

(where the minus sign indicates inward flow) and is subsonic.

3. The disc dynamo – observable properties

The stars for which we have the most comprehensive understanding of the properties (thermal, magnetic) of the dynamo fluid, and for which we have the best handle on the global properties of the dynamo itself, are the subclass of the cataclysmic variable stars, known as dwarf novae. The cataclysmic variables are binary stars, consisting of a low-mass, solar-type star which is losing its outer layers to its companion. The companion is a more massive, but more compact star, approximately the size of the Earth, known as a white dwarf (Warner 1995). Because of angular momentum considerations, the flow takes the form of an accretion disc (see, for example, Pringle & Wade 1985). The dwarf novae show two states: (i) a bright outburst state in which the mass accretion rate onto the white dwarf is high and the disc is highly ionised, with low magnetic diffusivity; and (ii) a dim quiescent state in which the accretion rate is low, the disc ionisation is low and the magnetic diffusivity relatively high.

3.1. Hot, highly ionised discs

The evolution of the surface density of an accretion disc is described by a diffusion equation, with the diffusion parameter proportional to the disc kinematic viscosity, $\nu$, that is, proportional to the Shakura–Sunyaev parameter $\alpha$ (see, for example, Lynden-Bell & Pringle 1974; Pringle 1981). During the decay from a dwarf nova outburst, the mass drains from the disc onto the central white dwarf. The time scale for this decay depends directly on the value of $\alpha$. The decay time scale of the outburst allows for a measurement of $\alpha$ in the hot state from modelling the outburst light curve. The disc size is known from the properties of the system and the disc temperature is obtained from the spectra. A simple calculation by Bath & Pringle (1981) suggested that $\alpha \approx 1$. Since then, there has been extensive and more detailed modelling of dwarf nova outburst behaviour, and all models point to relatively large values of $\alpha$. The models imply that $\alpha \approx 0.2\unicode{x2013}0.3$ (e.g. Pringle, Verbunt & Wade 1986; Smak 1998, 1999; Buat-Ménard, Hameury & Lasota 2001; Cannizzo 2001a,b; Schreiber, Hameury & Lasota 2003, 2004; Balman & Revnivtsev 2012; Kotko & Lasota 2012; Coleman et al. 2016).

It should be noted that these measurements are in line with the values of $\alpha$ deduced from time-dependent disc behaviour in other systems with highly ionised accretion discs, for example, the soft X-ray transients and the Be stars (see the discussion by Martin et al. 2019).

If, as we assume here, the dominant stresses that give rise to $\alpha$ are magnetic, this implies (see (2.5)) that the magnetic pressure in the disc is comparable to the gas pressure. Indeed, formally, since for the MHD dynamo driven by ($R \phi$) shear we expect that $\langle B_\phi ^2 \rangle \gg \langle B_R^2 \rangle$, it is evident that we require $\langle B_\phi ^2 \rangle / \rho c_s^2 \gtrsim 1$. In other words, the observations seem to imply that the $\beta$ parameter of the disc plasma is such that $\beta \lesssim 1$.

From the observed disc properties, we may estimate the physical value of the magnetic diffusivity, $\eta$, and hence the magnetic Reynolds number defined as

(3.1)\begin{equation} {\mathcal{R}}_{m} = c_s H /\eta, \end{equation}

at a typical point in the plane of a dwarf nova accretion disc in outburst.

We take the central star to have a mass $M = 1 M_\odot = 2 \times 10^{33}$ g, and consider a typical radius in the disc, $R = 10^{10}$ cm. For a highly ionised disc, in the bright state of a dwarf nova, we take a typical accretion rate of ${\dot M} = 10^{18}$ g s$^{-1}$ and we take the dimensionless viscosity parameter to be $\alpha = 0.3$ in line with observations.

We evaluate $\varOmega _0$ as

(3.2)\begin{equation} \varOmega_0 = \left(\frac{GM}{R^3} \right)^{1/2} = 1.2 \times 10^{{-}2}\,{\rm radian}\,{\rm s}^{{-}1}. \end{equation}

For the magnetic diffusivity, we assume the usual Spitzer value of $\eta = 3.5 \times 10^{12}\,T^{-3/2}$ cm$^2$ s$^{-1}$ (Spitzer 1962; Potter & Balbus 2014). From Frank et al. (2002), we find that the temperature in the plane of the disc is $T_{c} = 7.1 \times 10^4$ K. Thus, we have that $\eta = 1.9 \times 10^5$ cm$^2$ s$^{-1}$. Similarly, we find that the disc scale height is given by $H = 3.8 \times 10^8$ cm, so that the disc opening angle is $H/R = 0.038$.

From these, we are able to estimate a typical magnetic Reynolds number as

(3.3)\begin{equation} {\mathcal{R}}_{m} = 9.1 \times 10^{9}. \end{equation}

In summary, we note that, in these discs, whatever the origin of the turbulent behaviour within the disc that gives rise to the observed effective viscosity, whether it is purely hydrodynamic or (as is generally believed) magnetohydrodynamic, the mechanism that produces it is able to drive the fluid motions only up to, or close to, the sound speed. The fact that $\alpha$ is always found to be close to this limit (for these discs) implies that whatever instability might give rise to the driving mechanism for the turbulence, the turbulent velocities grow until they become trans-sonic.

Thus, in agreement with the original conjecture of Shakura & Sunyaev (1973), saturation of the turbulence is achieved once the motions become trans-sonic. This must be the result of the fact that once the motions approach the sound speed, the nature of the turbulence changes in a fundamental fashion. In line with the ideas of Shakura & Sunyaev (1973, 1976), it is evident that the change in the nature of the turbulence might occur for one, or both, of two physical reasons. First, in the case of hydrodynamic turbulence, as the turbulence becomes trans-sonic, shocks begin to dominate the dissipative process. Second, in the case of MHD turbulence, once the Alfvén speed approaches the sound speed (i.e. once $\beta \approx 1$), not only do the turbulent velocities become trans-sonic, but, in addition, the time scale for the Parker instability (leading to loss of magnetic flux from the disc) becomes comparable with the shearing time scale (growth time scale for magnetic flux) $\approx 1/\varOmega$ (cf. Tout & Pringle 1992).

The corollary of this basic finding is that numerical simulations of disc turbulence (for highly ionised discs) which do not find that the strength of the turbulence grows until limited by the sound speed (and which therefore do not find the large values of $\alpha$ implied by the observational data) cannot provide an adequate description of observed accretion disc dynamos. It seems likely that such models must be missing some fundamental physics (King et al. 2007).

3.2. Cool discs

The value of $\alpha$ in the low state dwarf nova accretion disc is difficult to determine, as the disc in that state shows little in the way of time evolution, and what time evolution there is appears to be dominated by the continued addition of mass to the disc from the companion star. In addition, estimates of $\alpha$ in the quiescent disc require modelling of the complete outburst cycle. However, all models of dwarf nova outburst cycles seem to require that in the quiescent disc, the value of $\alpha$ is less than that found in the outburst disc by at least a factor of ${\approx }10$ (Meyer & Meyer-Hofmeister 1983; Pringle et al. 1986; Cannizzo 1993, 2001b; Coleman et al. 2016).

These findings are in line with the idea (suggested for the dwarf novae quiescent discs by Gammie & Menou 1998) that the driving force for MHD disc turbulence, namely the MRI, is suppressed once the magnetic diffusivity becomes too large. A number of (local, shearing box) simulations suggest that the MRI becomes inoperative once ${\mathcal {R}}_{m} \lesssim 10^3\unicode{x2013}10^4$ (Hawley, Gammie & Balbus 1996; Stone et al. 1996; Fleming, Stone & Hawley 2000; Davis, Stone & Pessah 2010). These results strengthen the argument that the anomalous viscosity (at least in dwarf novae) is an MHD effect (Martin et al. 2019).

4. The disc dynamo – numerical simulations

4.1. Local simulations – the shearing box

Most simulations of disc dynamos make use of the shearing box approximation (Hawley, Gammie & Balbus 1995). Here the computational domain is a Cartesian box of size ${\sim }H \ll R$ co-moving with the fluid at some fixed radius, $R_0$, where the angular velocity is $\varOmega _0$. Thus, $(R, \phi, z) \rightarrow (x = R - R_0, y = \phi - \varOmega _0 t, z)$. The base flow in the box is a linear shear $u_y = U^\prime x$, where in a Keplerian disc, $U^\prime = \frac {3}{2} \varOmega _0$ and the box rotates with angular velocity $\varOmega _0$. In the simulations, the disc may, or may not, be stratified in the $z-$direction.

Typically, the simulations start with a small initial field. For example, Brandenburg et al. (1995) and Hawley et al. (1996) take an initial poloidal field, with zero net flux, of the form $B_z \propto \sin k x$. The early results are summarised in a review by Balbus & Hawley (1998). The most important finding was that with small initial seed fields, a steady, turbulent MHD disc dynamo could occur. However, for seed fields which had no externally imposed net field (i.e. for those simulations relevant to astrophysical discs), the magnitude of the $R \phi$ magnetic stress was around $\alpha \sim 0.01$, approximately an order of magnitude smaller than that required by observations. Balbus & Hawley (1998) conceded that while the dynamo saturation with $\alpha \approx 1$ postulated by Shakura & Sunyaev (1973) was an attractive physical possibility, this was not what was found in the simulations. They concluded, ‘It appears likely, therefore, that there is a dynamo regime that is characterised by unstable growth continuously balancing dissipation scale losses. This leads to subthermal magnetic fields and a dimensionless stress tensor of order $\alpha \approx 0.01$. Whether there are other modes of dynamo operation that arise naturally in accretion disks – at different magnetic Prandtl numbers, for example – is a fascinating and completely openquestion.’

In the 25 years or so since then, there have been many shearing box simulations of the accretion disc dynamo, using a variety of codes both grid-based and spectral, and a variety of physical parameters, with steadily increasing sophistication and resolution. The net result has been a much deeper understanding of the inner workings of the shearing box dynamo process, but the conclusions are unchanged. The simulations that assume zero externally imposed net magnetic flux typically find values of $\alpha$ at most around $\alpha \approx 0.01$ and often less. For example, Fromang et al. (2007) using both grid-based and spectral codes obtained $\alpha \approx 0.001$ in their best resolved simulations, Heinemann & Papaloizou (2009) use a spectral code and find $\alpha \approx 0.005$, Davis et al. (2010) using a grid-based code find $\alpha \approx 0.01$, Salvesen et al. (2016) using a grid-based code find $\alpha \approx 0.01$, Walker, Lesur & Boldyrev (2016) using a spectral code find that with no net imposed flux dynamo, the activity decays, so that $\alpha = 0$, and Mamatsashvili et al. (2020) using a spectral code find $\alpha \approx 0.006\unicode{x2013}0.01$. It is clear that such low values are not in agreement with observations. However, there is a discussion of this tension by Shi, Stone & Huang (2016) who present grid-based shearing box simulations. There they find that higher values ($\alpha \approx 0.1$) can be obtained in an unstratified shearing box with an unusually large height-to-width ratio. Whether these larger values of alpha in tall, unstratified boxes carries over to the more realistic stratified case has been questioned (e.g. Ryan et al. 2017).

In summary, while such numerical simulations have successfully demonstrated the existence of a self-sustaining disc dynamo, they have not been able to produce magnetic fields of sufficient magnitude to agree with the observational data. Typically, they produce values of $\alpha \approx 0.01$ or less, much smaller than the values of $\alpha \approx 0.2\unicode{x2013}0.3$ required to account for the observational data. So far, there is no explanation of why the value of $\alpha \approx 0.01$ emerges from the majority of the shearing box dynamo simulations. A physical understanding of what gives rise to the saturation of the dynamo process in numerical simulations would be of great value in understanding this difference between simulated and observed accretion discs.

4.2. Global simulations

King et al. (2007) drew attention to the discrepancy between the magnitudes of the disc magnetic fields that are required by observations and those that can be achieved in numerical, shearing box simulations. They discussed various reasons as to why the numerical simulations at that time might be inadequate. They noted that the shearing box approximation limits azimuthal structures to azimuthal wavenumbers $m=0$ or $m \gtrsim R/H \gg 1$, and suggested that this might play a role in suppressing large-scale azimuthal fields. For example, Tout & Pringle (1992) note that in their semi-analytic dynamo model, which produces fields with $\beta \sim 1$, magnetic buoyancy plays an important role in the conversion of toroidal field to poloidal and acts at toroidal wavelengths $\lambda _\phi \gg H$. Thus, it may be that it is the localisation of dynamo activity in the shearing box approximation that is responsible for the small values of $\alpha$ that are obtained in such simulations. We therefore need to consider global simulations.

We consider two recent global, grid-based, numerical simulations presented by Pjanka & Stone (2020) and Oyang, Jiang & Blaes (2021), which are specifically designed to simulate the accretion discs in cataclysmic variables. Our main interest here is the numerical diffusivity present in such simulations. In both of these simulations, material is introduced at the outer edge of the disc (to model the stream of material from the low mass star) and is taken to contain small loops (size ${\sim }H$) of magnetic flux, of magnitude $\beta \gg 1$ which can then initialise dynamo action. This seems reasonable, since the low mass star is sun-like and so presumably has surface magnetic fields (Livio & Pringle 1994). Oyang et al. (2021) also experiment with initially introducing small loops of flux at all radii in the disc.

The disc thicknesses in dwarf novae are typically such that $H/R \approx 0.02\unicode{x2013}0.05$ (see, for example, Frank et al. 2002). The thinner disc in the simulations of Pjanka & Stone (2020) has $H/R = 0.1$ and for that disc, they find in the body of the disc that $\alpha \approx 0.003$. The disc of Oyang et al. (2021) has $H/R = 0.044$ and they find values of $\alpha \approx 0.01$. Thus, neither of these global disc dynamo simulations is capable of accounting for the observed values of $\alpha$.

5. What is missing?

We have seen that the shearing box approach does not seem able to produce a saturation level to dynamo action which involves strong enough magnetic fields to agree with the observational constraints. We have noted that one possibility for this might be that the shearing box approximation itself is too small scale (particularly in the azimuthal direction) to be able to provide an accurate description of dynamo activity in a global accretion disc.

However, we have also noted that those global disc simulations presented so far are also unable to produce the necessary saturation level to explain the observed dynamo action. We speculate here that in the global models, this might come about for a different reason. The dominant physical process by which shear energy is converted into magnetic energy in an accretion disc is by the $R \phi$ shear acting on the radial $B_R$ field and converting it to an azimuthal field, $B_\phi$. It is this process that ultimately drives the dynamo. If that physical process is artificially restricted, for example, by too high a magnetic diffusivity, then the saturation level of dynamo activity might also be restricted. (We thank one of the referees for pointing out that this consideration may also be an issue for shearing box simulations.)

Numerical simulations are, of course, an approximation to reality. One way in which numerical simulations differ from reality is the inescapable inclusion of numerical dissipation, e.g. numerical viscosity and numerical magnetic diffusivity that are present in all numerical calculations (either explicitly as additional terms in the equations of motion or implicitly through e.g. truncation error at the grid scale). Because of these effects, it is not possible to provide a solution to the equations of inviscid hydrodynamics with a numerical approach, as the calculation will necessarily involve some level of numerical viscosity. Similarly, it is not possible to solve the equations of ideal MHD with a numerical approach, as the calculation will necessarily involve numerical magnetic diffusivity. However, it seems to be typical to ‘solve the equations of ideal MHD’ with a numerical code without providing discussion on, or estimates of, the magnitude of these numerical effects (see, for example, Pjanka & Stone 2020; Oyang et al. 2021). (Note that, in contrast, e.g. Fromang et al. (2007), Walker et al. (2016), Gogichaishvili et al. (2017) and Mamatsashvili et al. (2020) use spectral codes in which they explicitly include viscosity and magnetic diffusivity of sufficient magnitude that their effects can be resolved numerically. However, such spectral codes are not readily applicable to global disc simulations. See also, e.g. Fromang et al. (2007) for the use of explicit dissipation in non-spectral approaches.) It is of course possible to make estimates of the numerical viscosity and magnetic diffusivity. The magnitudes of these depend typically on the sizes of the grid cells, the time steps, and the orders of the numerical code in both space and time (see, for example, Rembiasz et al. 2017). However, it is rare that sufficient information is presented for such estimates to be made by the reader.

As an illustrative example, we consider the evolution of a loop of magnetic flux of size ${\approx }H$ in the numerical codes employed in these global simulations. The simple question we ask is: what is the maximum possible flux amplification that can be achieved in these codes, ignoring any instabilities (such as MRI, buoyancy etc.) which might limit the effect? Ideally, one would make use of the detailed code quantities, together with formulae for numerical magnetic diffusivity, $\eta _N$, such as those proposed by Rembiasz et al. (2017). However, not enough information is provided in these papers to undertake such an analysis. For this reason, we shall content ourselves with the following approximation.

For the loop we are considering, the distance in the azimuthal direction over which the radial component of $\boldsymbol B$ reverses is ${\approx }H$. After a time $t$, the angle an initially radial field makes with the azimuthal direction is $\approx (U^\prime t)^{-1}$. If the radial grid cell has size ${\rm \Delta} R$, then after a time $t_{\rm max}$, where $U^\prime t_{\rm max} = H/(2 {\rm \Delta} R)$, the grid cell should, in principle, contain magnetic fluxes of opposing signs. We may assume that at this point, numerical magnetic diffusivity ensures that this does not happen. We also note that at that time, the initial magnetic field has been amplified by a factor of ${\approx }U^\prime t_{\rm max}$. Thus, as a first approximation, we may assume that in an Eulerian numerical code, the maximum amplification of the field energy is given approximately by

(5.1)\begin{equation} {\mathcal{E}}_B(t_{\rm max})/{\mathcal{E}}_B(0) = (B(t_{\rm max}) / B_0)^2 \approx \left(\frac{H}{2 {\rm \Delta} R}\right)^2. \end{equation}

The codes used by Pjanka & Stone (2020) and Oyang et al. (2021) respectively employ adaptive mesh refinement (AMR) and static mesh refinement (SMR), and thus have variable cell size. Pjanka & Stone (2020), for the thinner disc case (which has disc thickness closer to that required for modelling dwarf nova discs), take $H = 0.1 R$, the initial cells have ${\rm \Delta} R/R = 0.1$ and there are up to five levels of cell refinement. Thus, in principle, ${\rm \Delta} R$ might be reduced by a factor of up to 32. In their model, this would imply ${\rm \Delta} R/R \approx 0.003$. From this, we deduce that the maximum achievable enhancement of magnetic field energy is ${\approx }16^2$. Oyang et al. (2021) use up to 256 cells logarithmically spaced between the inner radius $r=1$ and the outer radius $r=23.4$. This gives a smallest radial cell size of ${\rm \Delta} R/R \approx 0.012$, to be compared with the disc scale height $H/R = 0.044$. Thus, here the maximum field energy enhancement is by a factor of ${\approx }4$.

To relate these numbers to the actual numerical viscosities present in the codes, we now present a formal computation of the evolution of such field loops in a simple linear shear flow.

6. Formal computation of loop evolution

To illustrate the effect of magnetic diffusivity on magnetic field amplification, we consider a simple two-dimensional flow with differing initial magnetic field configurations. We have noted that the fundamental mechanism by which magnetic energy is increased, by tapping the energy available in the background $R\phi$-shear flow, is the conversion (stretching) of the radial magnetic field to form an azimuthal field. Pjanka & Stone (2020) and Oyang et al. (2021) used small loops of magnetic flux to seed the radial field. Alternatively, one may think of the effect of the MRI on an initially azimuthal field as producing radial undulations of the field, which can be thought of as loops of flux in the $R\phi$-plane. Such loops would then be stretched by the underlying $R\phi$-shear flow.

To contrast with the usual shearing box nomenclature (in which rotation, i.e. coriolis forces, are included), we shall consider two dimensions, but take the $x$-coordinate to correspond to the ‘azimuthal’ direction and the $y$-coordinate to the ‘radial’ direction. Thus, in the $xy$-plane, we assume an inexorable linear shear flow of the form ${\boldsymbol u} = (U^\prime y, 0)$ with $U^\prime$ a constant. First we consider the simple case where the initial field has only a component in the $y$-direction. We then consider the evolution of initial magnetic field loops. The exact general solution for arbitrary initial conditions can be found in Appendix A.

6.1. Straight (radial) field lines

We take an initial magnetic field, at time $t=0$, of the form ${\boldsymbol B} = B_0 (0, \cos (kx))$. This can be thought of as an approximation to a magnetic loop of size $l = {\rm \pi}/k$ in a shearing (but not rotating) box of size ${\sim }l$. We assume the fluid to be incompressible and to obey the standard MHD equations with a magnetic diffusivity $\eta$.

In this case, we can define the magnetic flux function $A(x,y, t)$ such that (Davidson 2001; Yeates & Hornig 2011)

(6.1)\begin{equation} {\boldsymbol B} = \left( \frac{\partial A}{\partial y}, - \frac{\partial A}{\partial x} \right), \end{equation}

which, in this case, obeys the equation

(6.2)\begin{equation} \frac{\partial A}{\partial t} + U^\prime y \frac{\partial A}{\partial x} = \eta \nabla^2 A. \end{equation}

The flux function is such that the magnetic field lines are given by $A =$ const.

At time $t=0$, we have that

(6.3)\begin{equation} A(x,y,t=0) ={-} k^{{-}1} B_0 \sin (kx). \end{equation}

If we have ideal MHD, so that $\eta = 0$, we expect the field lines to be simply advected by the shear flow, so that

(6.4)\begin{equation} A(x,y,t) ={-} k^{{-}1} B_0 \sin k(x - U^\prime t y). \end{equation}

In this case, the magnetic field grows indefinitely in a linear fashion, viz.

(6.5)\begin{equation} {\boldsymbol B} = B_0 (U^\prime t, 1) \cos k (x - U^\prime ty). \end{equation}

For non-ideal MHD, with magnetic diffusivity $\eta > 0$, we can see that the solution is separable, and of the form

(6.6)\begin{equation} A = f(t) \sin k(x - U^\prime ty). \end{equation}

Substituting this into (6.2), we find that

(6.7)\begin{equation} \frac{{\rm d} f}{{\rm d} t} ={-} \eta f k^2 [1 + (U^\prime)^2 t^2]. \end{equation}

Thus, we have that

(6.8)\begin{equation} A(x,y,t) ={-} k^{{-}1} B_0 \sin [k(x - U^\prime ty)] \exp \left\{- \eta k^2 \left[t + \tfrac{1}{3} (U^\prime)^2 t^3\right] \right\}, \end{equation}

and, therefore, that

(6.9)\begin{equation} {\boldsymbol B} = B_0 (U^\prime t, 1) \cos [k(x-U^\prime ty)] \exp \left\{- {\mathcal{R}}_{m}^{{-}1} \left[U^\prime t + \tfrac{1}{3} (U^\prime t)^3\right] \right\}, \end{equation}

where ${\mathcal {R}}_{m} = U^\prime /(k^2 \eta )$ is the relevant magnetic Reynolds number.

Thus, the spatially averaged (e.g. over a box of size $-l \le x,y \le l$) value of the magnetic energy (${\mathcal {E}}_B$) is of the form (cf. Pringle, McMillan & Teaca 2017)

(6.10)\begin{equation} {\mathcal{E}}_B(t) / {\mathcal{E}}_B(0) = [1 + (U^\prime t)^2] \exp \left\{- 2{\mathcal{R}}_{m}^{{-}1} \left[U^\prime t + \tfrac{1}{3} (U^\prime t)^3\right]\right\}. \end{equation}

We plot this in figure 1 for various values of ${\mathcal {R}}_{m}$. From this, we can see that the evolution is as follows.

1. (i) Initially, the magnetic field energy starts to decay at a rate of $2U^\prime /{\mathcal {R}}_{m}$ due to magnetic diffusivity. This initial decay phase is brief.

2. (ii) Provided that ${\mathcal {R}}_{m}$ is large enough, the shear is strong enough to provide an enhancement of the magnetic field strength. For large ${\mathcal {R}}_{m}$, field growth occurs after a time $U^\prime t \approx {\mathcal {R}}_{m}^{-1}$. In this phase, the magnetic energy grows linearly and this is a direct consequence of the shear.

3. (iii) Eventually, the shear decreases the spatial length scale of the magnetic field sufficiently that diffusivity wins, and, as is well known, the field ultimately decays (cf. Weiss 1966; Moffatt 1978). The exponential decay is quicker than the simple $\exp (- \eta k^2 t)$, which is what happens if there is no shear, because the combination of shear and diffusivity hastens the process.

Figure 1. Evolution of the magnetic energy, scaled to the initial value, with time for several values of the magnetic Reynolds number, ${\mathcal {R}}_{m}$, for the case of initial radial field lines (6.10).(a) Values of $0.1 \le {\mathcal {R}}_{m} \le 10$ with the magnetic energy on a linear scale. (b) Values of ${\mathcal {R}}_{m}$ up to $10^7$ with the magnetic energy on a log scale. For small ${\mathcal {R}}_{m}$, the energy decays rapidly, while for large ${\mathcal {R}}_{m}$, the field initially decays before exhibiting growth and final decay.

The maximum possible enhancement of the magnetic field energy, ${\mathcal {E}}_B(t_{\rm max})/{\mathcal {E}}_B(0)$, depends on ${\mathcal {R}}_{m}$. We plot this in figure 2. For ${\mathcal {R}}_{m} \gg 1$, the maximum enhancement occurs when $U^\prime t \approx {\mathcal {R}}_{m}^{1/3}$, and is equal to ${\mathcal {E}}_B(t_{\rm max})/{\mathcal {E}}_B(0) \approx ({\mathcal {R}}_{m}/e)^{2/3}$.

Figure 2. The maximum growth factor of the magnetic field energy plotted as a function of the magnetic Reynolds number, ${\mathcal {R}}_{m}$, for the case of initial radial field lines (i.e. the maximum value attained from (6.10)). The red-dashed line indicates the prediction appropriate to the limit of large ${\mathcal {R}}_{m}$, which is accurate for ${\mathcal {R}}_{m} \gtrsim 10$. For ${\mathcal {R}}_{m} \lesssim 3.8$, the magnetic energy never grows back above the original value.

In figure 3, we plot the geometric evolution of the magnetic field lines at different times as an illustration; note that the value of ${\mathcal {R}}_{m}$ does not affect the geometry and only affects the field strength in this case. Panels (a)–(d) correspond to times of $U^\prime t = 0,0.1,1,10$, respectively.

Figure 3. Evolution of the magnetic field lines from an initially linear radial ($y$-direction) configuration. Time increases as $U^\prime t$ is: (a) 0; (b) 0.1; (c) 1 and (d) 10. The colour bar denotes the value of $A$, with lines of constant $A$ delineating the field lines. For clarity, we also plot five field lines, which at $U^\prime t = 0$ correspond to $A = -k^{-1}B_0 \sin (k\gamma )$ with $k={\rm \pi}$, $B_0=1$, and $\gamma = -2/3, -1/3, 0, 1/3, 2/3$ (the solid black lines represent those field lines with $A \ge 0$ and the dashed lines represent those with $A < 0$). In subsequent panels, only these field lines are plotted, with the apparent number increased due to stretching of the field lines in the $x$-direction. This figure illustrates the conversion of radial field ($y$-direction) into azimuthal field ($x$-direction). The strength of the field varies in time and the evolution of the strength depends on ${\mathcal {R}}_{m}$ (see figure 1). Here we have plotted the case where $\eta \rightarrow 0$ (corresponding to ${\mathcal {R}}_{m} \rightarrow \infty$) and, as such, the range of values for $A$ remains fixed (cf. (6.4)).

6.2. Loops of magnetic field

We now explore the behaviour of loops of magnetic flux. To stay with one Fourier mode in each direction, we take

(6.11)\begin{equation} A(x,y,t=0) = \sin k_x x \sin k_y y. \end{equation}

In the region $0 \le k_x x,\ k_y y \le {\rm \pi}$, the magnetic field takes the form of loops around the point $({\rm \pi} /2, {\rm \pi}/2)$ with the field pointing in the anti-clockwise direction. The rest of space is populated with similar rectangular tiles containing loops of magnetic field of alternating chirality.

We note that we may rewrite the flux function as

(6.12)\begin{equation} A(x,y,t=0) = \tfrac{1}{2} [ \cos (k_x x- k_y y) - \cos (k_x x+k_y y) ]. \end{equation}

Thus, with zero diffusivity, we have simple advection of the field by the shear so that

(6.13)\begin{equation} A(x,y,t) = \tfrac{1}{2} [ \cos (k_x x - U^\prime t k_x y - k_y y) - \cos (k_x x - U^\prime t k_x y + k_y y)]. \end{equation}

Because we have have a combination of two Fourier modes, for $\eta > 0$, we now look for solutions of the form

(6.14)\begin{equation} A(x,y,t) = \tfrac{1}{2} [ f_1(t) \cos (k_x x - U^\prime t k_x y - k_y y) - f_2(t) \cos (k_x x - U^\prime t k_x y + k_y y)], \end{equation}

which, as before, can be substituted into the evolution equation (6.2) and solved for $f_1(t)$ and $f_2(t)$. Using the initial conditions that at $t=0$, $f_1(0) = f_2(0) = 1$, we have

(6.15)\begin{equation} f_1(t) = \exp \left\{ - \eta \left[ (k_x^2 + k_y^2)t + k_x k_y U't^2 + \tfrac{1}{3}k_x^2 (U')^2 t^3 \right] \right\}, \end{equation}

and

(6.16)\begin{equation} f_2(t) = \exp \left\{ - \eta \left[ (k_x^2 + k_y^2)t - k_x k_y U't^2 + \tfrac{1}{3}k_x^2 (U')^2 t^3\right] \right\}. \end{equation}

Putting this all together, we have

(6.17)\begin{align} A(x,y,t) & = [ \sin k_x(x - U'ty) \sin k_yy \cosh ( \eta k_x k_y U't^2) \nonumber\\ & \quad -\cos k_x (x - U'ty) \cos k_y y \sinh ( \eta k_x k_y U't^2) ] \nonumber\\ & \quad \times \exp \left\{ - \eta \left[(k_x^2 + k_y^2)t + \tfrac{1}{3} k_x^2 (U')^2 t^3\right] \right \}. \end{align}

In this case, the evolution of the spatially averaged magnetic field energy (e.g. over a box of size $-{\rm \pi} \le k_x x, k_y y \le {\rm \pi}$) is

(6.18)\begin{equation} \frac{{\mathcal{E}}_B(t)}{{\mathcal{E}}_B(0)} = \frac{[k_x^2 + (k_y + k_x U't)^2] f^2_1(t) + [ k_x^2 + (k_y - k_x U't)^2] f^2_2(t)}{2 [k_x^2 + k_y^2]}. \end{equation}

As before, we define the magnetic Reynolds number as ${\mathcal {R}}_{m} = U^\prime /(k_x^2\eta )$, but here, a secondary number ($k_y/k_x$) is required to define the flow. For illustration, we make the following figures with $k_y/k_x = 1$. First we plot the evolution of the magnetic energy (6.18) in figure 4. The result is similar to that for an initial straight field. For small ${\mathcal {R}}_{m}$, the field energy decays rapidly. For large ${\mathcal {R}}_{m}$, there is a period of significant growth of the field energy before the field later decays. For loops, the critical ${\mathcal {R}}_{m}$ at which the field exhibits any growth is significantly larger than for the straight field lines, requiring ${\mathcal {R}}_{m} \gtrsim 14.7$ for loops rather than just ${\mathcal {R}}_{m} \gtrsim 3.8$ for lines.

Figure 4. Evolution of the magnetic energy, scaled to the initial value, with time for several values of the magnetic Reynolds number, ${\mathcal {R}}_{m}$, for the case of initial loops of magnetic field (6.18). (a) Values of $0.1 \le {\mathcal {R}}_{m} \le 32$ with the magnetic energy on a linear scale. (b) Values of ${\mathcal {R}}_{m}$ up to $10^7$ with the magnetic energy on a log scale. For small ${\mathcal {R}}_{m}$, the energy decays rapidly, while for large ${\mathcal {R}}_{m}$, the field initially decays before exhibiting growth and final decay. The behaviour is similar, particularly at large ${\mathcal {R}}_{m} \gg 1$, with larger ${\mathcal {R}}_{m} \gtrsim 14.7$ required to exhibit field energy growth.

In figure 5, we plot the maximum growth of the magnetic energy against ${\mathcal {R}}_{m}$. In this figure, the black solid line corresponds to maximum growth from initial field loops, while the red-dashed lines are for the case of field lines plotted in figure 2. Here we can see that the maximum growth is weaker for loops of field compared to field lines for the same value of ${\mathcal {R}}_{m}$. However, we can also see that both cases exhibit the same functional dependence of the maximum growth rate on ${\mathcal {R}}_{m}$; for field loops with $k_x = k_y = 1$, the maximum growth in magnetic energy is ${\mathcal {E}}_B(t_{\rm max})/{\mathcal {E}}_B(0) \approx (1/2)({\mathcal {R}}_{m}/e)^{2/3}$. This is to be expected as at late times, where the maximum is reached for large ${\mathcal {R}}_{m}$, the initial field loops have already been sheared out into lines of field. To illustrate this, we plot the magnetic field lines resulting from the initial loops of field in figure 6. This figure again shows the conversion of radial to azimuthal field, and that for $U^\prime t \gg 1$, the field structure is similar to the case with initial lines of field as the loops become strongly stretched in the $x$-direction.

Figure 5. The maximum growth factor of the magnetic field energy plotted as a function of the magnetic Reynolds number, ${\mathcal {R}}_{m}$, for the case of initial loops of magnetic field (i.e. the maximum value attained from (6.18)). The red-dashed line shows the solution presented in figure 2 for the case of initial radial field lines. Here, for ${\mathcal {R}}_{m} \lesssim 14.7$, the magnetic energy never grows back above the original value, which contrasts with the value of ${\mathcal {R}}_{m} \lesssim 3.8$ in the initial radial field line case.

Figure 6. Evolution of the magnetic field lines from initial loops of magnetic field with $k_x = k_y = {\rm \pi}$. Time increases as $U^\prime t$ is: (a) 0; (b) 0.1; (c) 1 and (d) 10. The colour bar denotes the value of $A$, with lines of constant $A$ delineating the field lines. The regions with $A > 0$ (yellow) correspond to regions where the magnetic field lines are oriented in the counterclockwise direction, while for regions with $A < 0$ (blue), the magnetic field lines are oriented with the opposite chirality. For clarity, we also plot field lines corresponding to $A = -3/4, -1/2, -1/4, 1/4, 1/2, 3/4$ (the solid black lines represent those field lines with $A \ge 0$ and the dashed lines represent those with $A < 0$) with the exception of panel (d) where only field lines with $A=-1/2$ and $1/2$ are plotted for clarity. The strength of the field varies in time and the evolution of the strength depends on ${\mathcal {R}}_{m}$ (see figure 4). Here we have plotted the case where $\eta \rightarrow 0$ (corresponding to ${\mathcal {R}}_{m} \rightarrow \infty$) and, as such, the range of values for $A$ remains fixed (cf. (6.4)). As with the initial linear field case (see figure 3), the field lines are stretched due to the shear. As time proceeds, the solution for initial loops of field takes on a similar geometry to the case with initial lines of field for $U^\prime t \gg 1$.

7. Discussion and conclusions

It has long been known that shearing box numerical simulations of the dynamo dynamics to be found in accretion discs do not provide an explanation of the observational data. In particular, for the highly ionised discs in outbursting dwarf novae, the dimensionless viscosity parameter exceeds that obtained in numerical simulations by more than an order of magnitude (King et al. 2007; Martin et al. 2019). We note that in the spectral code implementation of the shearing box, the magnetic diffusivity is controlled and magnetic Reynolds numbers in the range $10^4\unicode{x2013}10^5$ can be achieved (Fromang et al. 2007; Walker et al. 2016; Mamatsashvili et al. 2020). As we have seen, the physical value of the magnetic Reynolds number in the observed dwarf nova discs is ${\sim }10^{10}$. Thus, while the lack of agreement between the shearing box simulations and the observational data may yet be due to a discrepancy in magnetic diffusivity, it may also be due to the shearing box approximation itself. If so, then the next step is to undertake global disc simulations.

Such simulations have been recently presented by Pjanka & Stone (2020) and Oyang et al. (2021). However, of necessity, the numerical methods used by Pjanka & Stone (2020) and Oyang et al. (2021) have intrinsic numerical magnetic diffusivities. We estimate (§ 5) that the maximum factor by which the field strength of a small magnetic loop can be enhanced in these codes is $\sim 16$, which corresponds to a factor of $16^2$ in the magnetic field energy. From § 6, we see that ${\mathcal {E}}_B(t_{\rm max})/{\mathcal {E}}_B(0) \approx ({\mathcal {R}}_{m}/e)^{2/3}$, suggesting that growth of the magnetic field energy by a factor of $16^2$ would correspond to a magnetic Reynolds number of ${\mathcal {R}}_{m} \sim 10^4$ (for Pjanka & Stone (2020), and less for Oyang et al. 2021). This implies that the numerical magnetic diffusivities, $\eta _N$, of these codes, compared to the values physically expected in such discs, $\eta$, are too high by factors of at least six orders of magnitude.

Pjanka & Stone (2020) provide a brief, preliminary discussion of the difficulty of accessing the physical parameter space using current numerical techniques, limited by current numerical processing power. They conclude that, nevertheless, ‘the global structure and behavior of these models should be reflected properly’. This is hard to square with the large discrepancy between the models and the reality of the fundamental measured disc parameter, $\alpha$. In contrast, Oyang et al. (2021) note that the size of viscosity produced by MHD processes in their simulation is too small to be able to account for the superhump phenomenon that they are investigating. It is only once they artificially increase the viscosity by approximately an order of magnitude that they are able to account for the phenomenon.

Furthermore, the fact that the observed values of the viscosity parameter imply values of the plasma parameter $\beta \sim 1$ means that magnetic buoyancy effects may play a far greater role than is found in the simulations (Tout & Pringle 1992). Balbus & Hawley (1998) noted that the simulations vastly overestimate the scale at which resistive losses occur. We agree with this assessment. In view of all this, we suggest that the large discrepancy mentioned above may well imply that the nature of the dynamos found in the simulations is fundamentally different from that which occurs in real accretion discs.

We have speculated that in the global numerical simulations, it may be the large numerical diffusivities in those simulations that present the problem. If that is the case, it will be important to discover how close the numerical magnetic Reynolds number needs to be to the physical value of ${\mathcal {R}}_{m} \approx 10^{10}$ in order that global disc dynamo simulations can achieve the required values of $\beta \sim 1$. Given the current limitations on computing power, it may be that expecting to be able to compute realistic dynamo action in observable accretion discs using numerical MHD codes is, for the time being, a step too far.