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

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.


Introduction
Fluid-based magnetic dynamos can be thought of as coming in two flavours: smallscale 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 lengthscale, 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 ≫ 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 lengthscales 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 et al. (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 = 0) are not captured (see also, for example, the discussion in 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 to those produced in real systems (cf.Tobias & Cattaneo 2013).
In Section 2 we briefly describe the relevant background to accretion disc physics and the relevant observed properties of such discs.In Section 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 socalled anomalous viscosity.In these objects, if this viscosity is magneto-hydrodynamic in origin, then the observations imply that the magnetic fields are strong (plasma β ∼ 1).In Section 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 (β ≫ 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 Section 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 Section 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 Section 7.

Accretion Discs
The theory of accretion discs is set out by Shakura & Sunyaev (1973, see also Pringle 1981;Frank et al. 2002).The basic disc flow, in cylindrical polar (R, φ, z) coordinates, is in the azimuthal direction, where M is the central gravitating mass.The disc thickness or scaleheight, H, in the z-direction is given by where c s is the local sound speed.For the usual thin disc approximation we require that H/R ≪ 1.Thus the azimuthal motion is supersonic, with Mach number ∼ R/H.Inflow (i.e.accretion) through the disc requires the action of a so-called "anomalous viscosity" which taps the azimuthal (Rφ) shear and transfers angular momentum outwards and mass inwards.The viscosity is generally assumed to be caused by small-scale, l 0 H, magneto-hydrodynamic turbulence within the disc, and it is the maintenance of this turbulence that requires dynamo action.Shakura & Sunyaev (1973) introduced a parameter α which is a dimensionless measure of the size of the anomalous viscosityessentially a dimensionless measure of the z-averaged Rφ-stress.Thus the (z-averaged) effective kinematic viscosity of the disc may be written as where is the angular velocity at radius R. For MHD turbulence, they note that (in appropriate units), where ρ is the disc density.They also introduced physical arguments as to why we might expect α 1 (see also the discussion in Martin et al. 2019).
In terms of α, the radial accretion flow velocity is (where the minus sign indicates inward flow), and is subsonic.

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 lowmass, solar-type star which is losing its outer layers to its companion.The companion is a more massive, but more compact star, about 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 ionized, with low magnetic diffusivity, and (ii) a dim quiescent state in which the accretion rate is low, the disc ionization is low and the magnetic diffusivity relatively high.

Hot, highly-ionized 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, ν, that is, proportional to the Shakura-Sunyaev parameter α (see, for example, Lynden-Bell & Pringle 1974;Pringle 1981).During the decay from a dwarf nova outburst, the mass drains from the disc on to the central white dwarf.The timescale for this decay depends directly on the value of α.The decay timescale of the outburst allows for a measurement of α in the hot state from modeling 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 α ≈ 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 α.The models imply that α ≈ 0.2 − 0.3 (e.g.Pringle et al. 1986;Smak 1998Smak , 1999;;Buat-Ménard et al. 2001;Cannizzo 2001a,b;Schreiber et al. 2003Schreiber et al. , 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 α deduced from time-dependent disc behaviour in other systems with highly ionized accretion discs, for example the soft X-ray transients and the Be stars (see the discussion in Martin et al. 2019).
If, as we assume here, the dominant stresses that give rise to α are magnetic, this implies (see equation 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φ) shear, we expect that B 2 φ ≫ B 2 R , it is evident that we require B 2 φ /ρc 2 s 1.In other words, the observations seem to imply that the β parameter of the disc plasma is such that β 1.
From the observed disc properties, we may estimate the physical value of the magnetic diffusivity, η, and hence the magnetic Reynolds number defined as 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 = 1M ⊙ = 2 × 10 33 g, and consider a typical radius in the disc, R = 10 10 cm.For a highly ionized disc, in the bright state of a dwarf nova, we take a typical accretion rate of Ṁ = 10 18 g/s, and we take the dimensionless viscosity parameter to be α = 0.3 in line with observations.We evaluate Ω 0 as For the magnetic diffusivity we assume the usual Spitzer value of η = 3.5 × 10 12 T −3/2 cm 2 /s (Potter & Balbus 2014;Spitzer 1962). From Frank et al. (2002) we find that the temperature in the plane of the disc is T c = 7.1 × 10 4 K. Thus we have that η = 1.9 × 10 5 cm 2 /s.Similarly we find that the disc scaleheight is given by H = 3.8 × 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 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) magneto-hydrodynamic, the mechanism that produces it is able to drive the fluid motions only up to, or close to, the sound speed.
The fact that α 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 β ≈ 1), not only do the turbulent velocities become trans-sonic, but, in addition, the timescale for the Parker instability (leading to loss of magnetic flux from the disc) becomes comparable with the shearing timescale (growth timescale for magnetic flux) ≈ 1/Ω (cf.Tout & Pringle 1992).
The corollary of this basic finding is that numerical simulations of disc turbulence (for highly ionized 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 α 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).

Cool discs
The value of α 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 α 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 α is less than that found in the outburst disc by at least a factor of ≈ 10 (Meyer & Meyer-Hofmeister 1983;Pringle et al. 1986;Cannizzo 1993Cannizzo , 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 R m 10 3 − 10 4 (Hawley et al. 1996;Stone et al. 1996;Fleming et al. 2000;Davis et al. 2010).These results strengthen the argument that the anomalous viscosity (at least in dwarf novae) is an MHD effect (Martin et al. 2019).

Local simulations -the shearing box
Most simulations of disc dynamos make use of the shearing box approximation (Hawley et al. 1995).Here the computational domain is a Cartesian box, of size ∼ H ≪ R co-moving with the fluid at some fixed radius, R 0 , where the angular velocity is 2 Ω 0 , and the box rotates with angular velocity Ω 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 ∝ sin kx.The early results are summarized 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φ magnetic stress was around α ∼ 0.01, about an order of magnitude smaller than that required by observations.Balbus & Hawley (1998) conceded that while the dynamo saturation with α ≈ 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 characterized by unstable growth continuously balancing dissipation scale losses.This leads to subthermal magnetic fields and a dimensionless stress tensor of order α ≈ 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 open question." 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 α at most around α ≈ 0.01 and often less.(2016) using a spectral code find that with no net imposed flux dynamo activity decays, so that α = 0, and Mamatsashvili et al. (2020) using a spectral code find α ≈ 0.006 − 0.01.It is clear that such low values are not in agreement with observations.However, there is a discussion of this tension in Shi et al. (2016) who present grid-based shearing box simulations.There they find that higher values (α ≈ 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 α ≈ 0.01 or less, much smaller than the values of α ≈ 0.2 − 0.3 required to account for the observational data.So far, there is no explanation of why the value of α ≈ 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.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 R/H ≫ 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 β ∼ 1, magnetic buoyancy plays an important role in the conversion of toroidal field to poloidal and acts at toroidal wavelengths λ φ ≫ 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 α that are obtained in such simulations.We therefore need to consider global simulations.

Global simulations
We consider two recent global, grid-based, numerical simulations, presented by Pjanka & Stone (2020) and Oyang et al. (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 ∼ H) of magnetic flux, of magnitude β ≫ 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 ≈ 0.02 − 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 α ≈ 0.003.The disc in Oyang et al. (2021) has H/R = 0.044, and they find values of α ≈ 0.01.Thus, neither of these global disc dynamo simulations is capable of accounting for the observed values of α.

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φ shear acting on the radial B R field and converting it to azimuthal field, B φ .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;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 timesteps, and the orders of the numerical code in both space and time (see, for example, Rembiasz et al. 2017).But 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 ≈ 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, η 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 B reverses is ≈ H.After a time t, the angle an initially radial field makes with the azimuthal direction is ≈ (U ′ t) −1 .If the radial grid cell has size ∆R, then after a time t max , where U ′ t max = H/(2∆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 ≈ U ′ t 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) The codes used by Pjanka & Stone (2020) and Oyang et al. (2021) employ, respectively, adaptive mesh refinement (AMR) and static mesh refinement (SMR), and thus have variable cell size.In Pjanka & Stone (2020) for the thinner disc case (which has disc thickness closer to that required for modelling dwarf nova discs), they take H = 0.1R, The initial cells have ∆R/R = 0.1, and there are up to 5 levels of cell refinement.Thus, in principle, ∆R might be reduced by a factor of up to 32.In their model this would imply ∆R/R ≈ 0.003.From this we deduce that the maximum achievable enhancement of magnetic field energy is ≈ 16 2 .In Oyang et al. (2021) there are 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 ∆R/R ≈ 0.012, to be compared with the disc scaleheight H/R = 0.044.Thus here the maximum field energy enhancement is by a factor of ≈ 4.
In order 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.

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φ−shear flow, is the conversion (stretching) of radial magnetic field to form azimuthal field.In Pjanka & Stone (2020) and Oyang et al. (2021) small loops of magnetic flux are used 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φ−plane.Such loops would then be stretched by the underlying Rφ−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 u = (U ′ y, 0). with U ′ 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.

Straight (radial) field lines
We take an initial magnetic field, at time t = 0 of the form B = B 0 (0, cos(kx)).This can be thought of as an approximation to a magnetic loop of size l = π/k in a shearing (but not rotating) box of size ∼ l.We assume the fluid to be incompressible, and to obey the standard MHD equations with a magnetic diffusivity η.
In this case we can define the magnetic flux function A(x, y, t) such that (Davidson 2001 The flux function is such that the magnetic field lines are given by A = const.At time t = 0 we have that If we have ideal MHD, so that η = 0, we expect the field lines to be simply advected by the shear flow, so that (6.4) In this case, the magnetic field grows indefinitely in a linear fashion, viz., For non-ideal MHD, with magnetic diffusivity η > 0, we can see that the solution is separable, and of the form Substituting this into Equation 6.2 we find that

Thus we have that
and, therefore, that where R m = U ′ /(k 2 η) is the relevant magnetic Reynolds number.Thus the spatially averaged (e.g. over a box of size −l x, y l) value of the magnetic energy (E B ) is of the form (cf. Pringle et al. 2017) We plot this in Figure 1, for various values of R m .From this we can see that the evolution is as follows: (i) Initially the magnetic field energy starts to decay at a rate of 2U ′ /R m due to magnetic diffusivity.This initial decay phase is brief.
(ii) Provided that R m is large enough, the shear is strong enough to provide an enhancement of the magnetic field strength.For large R m , field growth occurs after a time U ′ t ≈ R −1 m .In this phase the magnetic energy grows linearly and this is a direct consequence of the shear.
(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(−ηk 2 t), which is what happens if there is no shear, because the combination of shear and diffusivity hastens the process.
The maximum possible enhancement of the magnetic field energy, E B (t max )/E B (0), depends on R m .We plot this in Figure 2.For R m ≫ 1, the maximum enhancement occurs when U ′ t ≈ R 1/3 m , and is equal to E B (t max )/E B (0) ≈ (R m /e) 2/3 .In Figure 3 we plot the geometric evolution of the magnetic field lines at different The maximum growth factor of the magnetic field energy plotted as a function of the magnetic Reynolds number, Rm, for the case of initial radial field lines (i.e. the maximum value attained from equation 6.10).The red-dashed line indicates the prediction appropriate to the limit of large Rm, which is accurate for Rm 10.For Rm 3.8 the magnetic energy never grows back above the original value.
times as an illustration; note that the value of R m does not affect the geometry and only affects the field strength in this case.From left to right the panels correspond to times of U ′ t = 0, 0.1, 1, 10.

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 A(x, y, t = 0) = sin k x x sin k y y . (6.11) In the region 0 k x x, k y y π the magnetic field takes the form of loops around the point (π/2, π/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 Thus with zero diffusivity we have simple advection of the field by the shear so that Because we have have a combination of two Fourier modes, for η > 0 we now look for solutions of the form which, as before, can be substituted into the evolution equation (6.2) and solve 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 For clarity, we also plot 5 field lines, which at U ′ t = 0 correspond to A = −k −1 B0 sin(kγ) with k = π, B0 = 1, and γ = −2/3, −1/3, 0, 1/3, 2/3 (the solid black lines represent those field lines with A 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 Rm (see Figure 1).Here we have plotted the case where where η → 0 (corresponding to Rm → ∞), and as such the range of values for A remains fixed (cf.equation 6.4).and 6.16) Putting this all together we have In this case the evolution of the spatially averaged magnetic field energy (e.g. over a Rm 32 with the magnetic energy on a linear scale.Right panel shows values of Rm up to 10 7 with the magnetic energy on a log scale.For small Rm the energy decays rapidly, while for large Rm the field initially decays before exhibiting growth and final decay.The behaviour is similar, particularly at large Rm ≫ 1 with larger Rm 14.7 required to exhibit field energy growth.

box of size
. (6.18) As before, we define the magnetic Reynolds number as R m = U ′ /(k 2 x η), 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 (equation 6.18) in Figure 4.The result is similar to that for initial straight field.For small R m the field energy decays rapidly.For large R m there is a period of significant growth of the field energy before the field later decays.For loops the critical R m at which the field exhibits any growth is significantly larger than for the straight field lines, requiring R m 14.7 for loops rather than just R m 3.8 for lines.
In Figure 5 we plot the maximum growth of the magnetic energy against R m .In this figure the black solid line corresponds to maximum growth from initial field loops, while the red-dashed lines is 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 R m .However, we can also see that both cases exhibit the same functional dependence of the maximum growth rate on R m ; for field loops with k x = k y = 1 the maximum growth in magnetic energy is E B (t max )/E B (0) ≈ (1/2)(R m /e) 2/3 .This is to be expected as at late times, where the maximum is reached for large 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 ′ t ≫ 1 the field structure is similar to the case with initial lines of field as the loops become strongly stretched in the x-direction. .The maximum growth factor of the magnetic field energy plotted as a function of the magnetic Reynolds number, Rm, for the case of initial loops of magnetic field (i.e. the maximum value attained from equation 6.18).The red-dashed line shows the solution presented in Figure 2 for the case of initial radial field lines.Here, for Rm 14.7 the magnetic energy never grows back above the original value, which contrasts with the value of Rm 3.8 in the initial radial field line case.

Discussion & 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 ionized 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 −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 around ∼ 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 that the reason is 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 (Section 5) that the maximum factor by which the field strength of a small magnetic loop can be enhanced in these codes is ∼ 16, which corresponds to a factor of 16 2 in the magnetic field energy.From Section 6 we see that E B (t max )/E B (0) ≈ (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 R m ∼ 10 4 (for Pjanka & Stone 2020, and less for Oyang et al. 2021).This implies that the numerical magnetic diffusivities, η N , of these codes, compared to the values physically expected in such discs, η, are too high by factors of at least 6 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 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 0 and the dashed lines represent those with A < 0) with the exception of the final panel 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 Rm (see Figure 4).Here we have plotted the case where where η → 0 (corresponding to Rm → ∞), and as such the range of values for A remains fixed (cf.equation 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 ′ t ≫ 1. parameter, α.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 about 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 β ∼ 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 R m ≈ 10 10 in order that global disc dynamo simulations can achieve the required values of β ∼ 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.

Figure 1 .
Figure1.Evolution of the magnetic energy, scaled to the initial value, with time for several values of the magnetic Reynolds number, Rm, for the case of initial radial field lines (equation 6.10).Left panel shows values of 0.1 Rm 10 with the magnetic energy on a linear scale.Right panel shows values of Rm up to 10 7 with the magnetic energy on a log scale.For small Rm the energy decays rapidly, while for large Rm the field initially decays before exhibiting growth and final decay.

Figure 3 .
Figure3.Evolution of the magnetic field lines from an initially linear radial (y-direction) configuration.From left to right, and then top to bottom, time increases from U ′ t = 0 (upper left panel) up to U ′ t = 10 (lower right panel).The colour bar denotes the value of A, with lines of constant A delineating the field lines.For clarity, we also plot 5 field lines, which at U ′ t = 0 correspond to A = −k −1 B0 sin(kγ) with k = π, B0 = 1, and γ = −2/3, −1/3, 0, 1/3, 2/3 (the solid black lines represent those field lines with A 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 Rm (see Figure1).Here we have plotted the case where where η → 0 (corresponding to Rm → ∞), and as such the range of values for A remains fixed (cf.equation 6.4).

Figure 4 .
Figure 4. Evolution of the magnetic energy, scaled to the initial value, with time for several values of the magnetic Reynolds number, Rm, for the case of initial loops of magnetic field (equation 6.18).Left panel shows values of 0.1Rm 32 with the magnetic energy on a linear scale.Right panel shows values of Rm up to 10 7 with the magnetic energy on a log scale.For small Rm the energy decays rapidly, while for large Rm the field initially decays before exhibiting growth and final decay.The behaviour is similar, particularly at large Rm ≫ 1 with larger Rm 14.7 required to exhibit field energy growth.
Figure5.The maximum growth factor of the magnetic field energy plotted as a function of the magnetic Reynolds number, Rm, for the case of initial loops of magnetic field (i.e. the maximum value attained from equation 6.18).The red-dashed line shows the solution presented in Figure2for the case of initial radial field lines.Here, for Rm 14.7 the magnetic energy never grows back above the original value, which contrasts with the value of Rm 3.8 in the initial radial field line case.

Figure 6 .
Figure6.Evolution of the magnetic field lines from initial loops of magnetic field with kx = ky = π.From left to right, and then top to bottom, time increases from U ′ t = 0 (upper left panel) up to U ′ t = 10 (lower right panel).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 0 and the dashed lines represent those with A < 0) with the exception of the final panel 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 Rm (see Figure4).Here we have plotted the case where where η → 0 (corresponding to Rm → ∞), and as such the range of values for A remains fixed (cf.equation 6.4).As with the initial linear field case (see Figure3), 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 ′ t ≫ 1.