Hostname: page-component-848d4c4894-pftt2 Total loading time: 0 Render date: 2024-05-06T22:42:02.035Z Has data issue: false hasContentIssue false

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

Published online by Cambridge University Press:  05 January 2024

C.J. Nixon*
Affiliation:
School of Physics and Astronomy, University of Leeds, Sir William Henry Bragg Building, Woodhouse Ln., Leeds LS2 9JT, UK School of Physics and Astronomy, University of Leicester, University Road, Leicester LE1 7RH, UK
C.C.T. Pringle
Affiliation:
Centre for Fluid and Complex Systems, Coventry University, Coventry CV1 5FB, UK
J.E. Pringle
Affiliation:
Institute of Astronomy, University of Cambridge, Madingley Road, Cambridge CB3 0HA, UK
*
Email address for correspondence: c.j.nixon@leeds.ac.uk
Rights & Permissions [Opens in a new window]

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.

Type
Research Article
Creative Commons
Creative Common License - CCCreative Common License - BY
This is an Open Access article, distributed under the terms of the Creative Commons Attribution licence (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted re-use, distribution and reproduction, provided the original article is properly cited.
Copyright
Copyright © The Author(s), 2024. Published by Cambridge University Press

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 Reference Brandenburg and Subramanian2005; Rincon Reference Rincon2019; Schekochihin Reference Schekochihin2022). 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. Reference Padoan, Federrath, Chabrier, Evans II, Johnstone, Jørgensen, McKee, Nordlund, Beuther, Klessen, Dullemond and Henning2014; Federrath Reference Federrath2016) or in cosmological simulations (e.g. Martin-Alvarez et al. Reference Martin-Alvarez, Katz, Sijacki, Devriendt and Slyz2021. 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. Reference Brandenburg, Nordlund, Stein and Torkelsson1995; Balbus & Hawley Reference Balbus and Hawley1998), perhaps enhanced by magnetic buoyancy (Tout & Pringle Reference Tout and Pringle1992; Johansen & Levin Reference Johansen and Levin2008).

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 (Reference King, Pringle and Livio2007) 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 Reference Parkin and Bicknell2013). There also remain questions regarding the convergence of such models when they are extended to very high resolution (Bodo et al. Reference Bodo, Cattaneo, Mignone and Rossi2014). 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 Reference Tobias and Cattaneo2013).

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 (Reference Shakura and Sunyaev1973, see also Pringle Reference Pringle1981; Frank, King & Raine Reference Frank, King and Raine2002). 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 (Reference Shakura and Sunyaev1973) 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. Reference Martin, Nixon, Pringle and Livio2019). 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 Reference Warner1995). Because of angular momentum considerations, the flow takes the form of an accretion disc (see, for example, Pringle & Wade Reference Pringle and Wade1985). 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 Reference Lynden-Bell and Pringle1974; Pringle Reference Pringle1981). 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 (Reference Bath and Pringle1981) 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 Reference Pringle, Verbunt and Wade1986; Smak Reference Smak1998, Reference Smak1999; Buat-Ménard, Hameury & Lasota Reference Buat-Ménard, Hameury and Lasota2001; Cannizzo Reference Cannizzo2001a,Reference Cannizzob; Schreiber, Hameury & Lasota Reference Schreiber, Hameury and Lasota2003, Reference Schreiber, Hameury and Lasota2004; Balman & Revnivtsev Reference Balman and Revnivtsev2012; Kotko & Lasota Reference Kotko and Lasota2012; Coleman et al. Reference Coleman, Kotko, Blaes, Lasota and Hirose2016).

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. Reference Martin, Nixon, Pringle and Livio2019).

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 Reference Spitzer1962; Potter & Balbus Reference Potter and Balbus2014). From Frank et al. (Reference Frank, King and Raine2002), 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 (Reference Shakura and Sunyaev1973), 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 (Reference Shakura and Sunyaev1973, Reference Shakura and Sunyaev1976), 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 Reference Tout and Pringle1992).

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. Reference King, Pringle and Livio2007).

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 Reference Meyer and Meyer-Hofmeister1983; Pringle et al. Reference Pringle, Verbunt and Wade1986; Cannizzo Reference Cannizzo1993, Reference Cannizzo2001b; Coleman et al. Reference Coleman, Kotko, Blaes, Lasota and Hirose2016).

These findings are in line with the idea (suggested for the dwarf novae quiescent discs by Gammie & Menou Reference Gammie and Menou1998) 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 Reference Hawley, Gammie and Balbus1996; Stone et al. Reference Stone, Hawley, Gammie and Balbus1996; Fleming, Stone & Hawley Reference Fleming, Stone and Hawley2000; Davis, Stone & Pessah Reference Davis, Stone and Pessah2010). These results strengthen the argument that the anomalous viscosity (at least in dwarf novae) is an MHD effect (Martin et al. Reference Martin, Nixon, Pringle and Livio2019).

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 Reference Hawley, Gammie and Balbus1995). 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. (Reference Brandenburg, Nordlund, Stein and Torkelsson1995) and Hawley et al. (Reference Hawley, Gammie and Balbus1996) 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 (Reference Balbus and Hawley1998). 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 (Reference Balbus and Hawley1998) conceded that while the dynamo saturation with $\alpha \approx 1$ postulated by Shakura & Sunyaev (Reference Shakura and Sunyaev1973) 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. (Reference Fromang, Papaloizou, Lesur and Heinemann2007) using both grid-based and spectral codes obtained $\alpha \approx 0.001$ in their best resolved simulations, Heinemann & Papaloizou (Reference Heinemann and Papaloizou2009) use a spectral code and find $\alpha \approx 0.005$, Davis et al. (Reference Davis, Stone and Pessah2010) using a grid-based code find $\alpha \approx 0.01$, Salvesen et al. (Reference Salvesen, Simon, Armitage and Begelman2016) using a grid-based code find $\alpha \approx 0.01$, Walker, Lesur & Boldyrev (Reference Walker, Lesur and Boldyrev2016) using a spectral code find that with no net imposed flux dynamo, the activity decays, so that $\alpha = 0$, and Mamatsashvili et al. (Reference Mamatsashvili, Chagelishvili, Pessah, Stefani and Bodo2020) 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 (Reference Shi, Stone and Huang2016) 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. Reference Ryan, Gammie, Fromang and Kestener2017).

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. (Reference King, Pringle and Livio2007) 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 (Reference Tout and Pringle1992) 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 (Reference Pjanka and Stone2020) and Oyang, Jiang & Blaes (Reference Oyang, Jiang and Blaes2021), 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 Reference Livio and Pringle1994). Oyang et al. (Reference Oyang, Jiang and Blaes2021) 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. Reference Frank, King and Raine2002). The thinner disc in the simulations of Pjanka & Stone (Reference Pjanka and Stone2020) 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. (Reference Oyang, Jiang and Blaes2021) 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 Reference Pjanka and Stone2020; Oyang et al. Reference Oyang, Jiang and Blaes2021). (Note that, in contrast, e.g. Fromang et al. (Reference Fromang, Papaloizou, Lesur and Heinemann2007), Walker et al. (Reference Walker, Lesur and Boldyrev2016), Gogichaishvili et al. (Reference Gogichaishvili, Mamatsashvili, Horton, Chagelishvili and Bodo2017) and Mamatsashvili et al. (Reference Mamatsashvili, Chagelishvili, Pessah, Stefani and Bodo2020) 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. (Reference Fromang, Papaloizou, Lesur and Heinemann2007) 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. Reference Rembiasz, Obergaulinger, Cerdá-Durán, Aloy and Müller2017). 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. (Reference Rembiasz, Obergaulinger, Cerdá-Durán, Aloy and Müller2017). 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 (Reference Pjanka and Stone2020) and Oyang et al. (Reference Oyang, Jiang and Blaes2021) respectively employ adaptive mesh refinement (AMR) and static mesh refinement (SMR), and thus have variable cell size. Pjanka & Stone (Reference Pjanka and Stone2020), 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. (Reference Oyang, Jiang and Blaes2021) 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 (Reference Pjanka and Stone2020) and Oyang et al. (Reference Oyang, Jiang and Blaes2021) 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 Reference Davidson2001; Yeates & Hornig Reference Yeates and Hornig2011)

(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 Reference Pringle, McMillan and Teaca2017)

(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 Reference Weiss1966; Moffatt Reference Moffatt1978). 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. Reference King, Pringle and Livio2007; Martin et al. Reference Martin, Nixon, Pringle and Livio2019). 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. Reference Fromang, Papaloizou, Lesur and Heinemann2007; Walker et al. Reference Walker, Lesur and Boldyrev2016; Mamatsashvili et al. Reference Mamatsashvili, Chagelishvili, Pessah, Stefani and Bodo2020). 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 (Reference Pjanka and Stone2020) and Oyang et al. (Reference Oyang, Jiang and Blaes2021). However, of necessity, the numerical methods used by Pjanka & Stone (Reference Pjanka and Stone2020) and Oyang et al. (Reference Oyang, Jiang and Blaes2021) 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 (Reference Pjanka and Stone2020), and less for Oyang et al. Reference Oyang, Jiang and Blaes2021). 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 (Reference Pjanka and Stone2020) 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. (Reference Oyang, Jiang and Blaes2021) 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 Reference Tout and Pringle1992). Balbus & Hawley (Reference Balbus and Hawley1998) 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.

Acknowledgements

We thank both Referees and the Editor for useful comments and discussion.

Editor Alex Schekochihin thanks the referees for their advice in evaluating this article.

Funding

This work was supported by the Science and Technology Facilities Council (C.J.N., grant number ST/Y000544/1); and the Leverhulme Trust (C.J.N., grant number RPG-2021-380).

Declaration of interests

The authors report no conflict of interest.

Appendix A. Solution for arbitrary initial field structure

As in the main text, we are considering a flow in the $xy$-plane with an inexorable linear shear flow of the form ${\boldsymbol u} = (U^\prime y, 0)$, with $U^\prime$ a constant. We assume the fluid to be incompressible and to obey the standard MHD equations with a magnetic diffusivity $\eta$, so that the evolution equation for the magnetic flux function $A(x,y, t)$ is

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

We consider the evolution of $A$ on the domain $L_x \times L_y$. At time $t=0$, the flux function $A(x,y,0)$ can be expanded as a Fourier series in the form

(A2)\begin{equation} A(x,y,0) = \sum_n \sum_k a_{n,k} \exp ({\rm i} \alpha n y) \exp ({\rm i} \beta k x), \end{equation}

where $\alpha = 2 {\rm \pi}/ L_y$ and $\beta = 2 {\rm \pi}/ L_x$.

At later times we expect the form of $A$ to be

(A3)\begin{equation} A(x,y,t) = \sum_n \sum_k a_{n,k}(t) \exp ({\rm i} \alpha n y) \exp ({\rm i} \beta k[ x- U'ty]). \end{equation}

Substituting this into (A1), and using the fact that the modes evolve independently, we find that this can be satisfied provided that the coefficients $a_{n,k}(t)$ obey the equations

(A4)\begin{equation} \frac{{\rm d} a_{n,k}}{{\rm d} t} ={-} \eta [ (\alpha n - \beta k U't)^2 + \beta^2 k^2] a_{n,k}. \end{equation}

Thus, we find that

(A5)\begin{equation} a_{n,k}(t) = c_{n,k} \exp \tau_{n,k}(t) \end{equation}

where, for $k \ne 0$, we have

(A6)\begin{equation} \tau_{n,k}(t) = \eta \left[ \frac{1}{3 \beta k U'} (\alpha n - \beta k U't)^3 - \beta^2 k^2 t \right], \end{equation}

and $c_{n,k}$ are constants.

For the $k = 0$ mode (for which the contribution to $\boldsymbol {B}$ is independent of $x$), we find that

(A7)\begin{equation} a_{n,0} = c_{n,0} \exp [ - \eta \alpha^2 n^2 t ]. \end{equation}

Note that, in general, the $c_{n,k}$ can take any values, provided that $c_{n,k} = c^*_{-n,-k}$, where the $^*$ denotes the complex conjugate.

The full solution for $A(x,y,t)$ is therefore

(A8)\begin{equation} A(x,y,t) = \sum_n \sum_k c_{n,k} \exp(-\tau_{n,k}(t)) \exp ({\rm i} \alpha n y) \exp ({\rm i} \beta k[ x- U'ty]), \end{equation}

where

(A9)\begin{equation} \tau_{n,k}(t) = \begin{cases} \eta \alpha^2 n^2 t & {\rm if}\ k = 0, \\ \eta \left[ \dfrac{1}{3 \beta k U'} (\alpha n - \beta k U't)^3 - \beta^2 k^2 t \right] & {\rm if}\ k\ne 0. \end{cases} \end{equation}

Thus, for the particular case $A(x,y,t=0) = \sin (\beta x)$, we take the non-zero values of $c_{n,k}$ to be

(A10)\begin{equation} c_{n, k} = \left\{ \begin{array}{@{}ll@{}} -(1/2)i & n=0,k=1, \\ (1/2)i & n=0, k={-}1, \\ 0 & \mbox{otherwise}. \end{array}\right. \end{equation}

And for the case $A(x,y,t=0) = \sin (\beta x) \sin (\alpha y)$, we take

(A11)\begin{equation} c_{n,k} = \left\{ \begin{array}{@{}ll@{}} -(1/4) \exp[-\tau_{n,k}(0)] & (n,k) ={\pm}(1,1),\\ (1/4) \exp[-\tau_{n,k}(0)] & (n,k) ={\pm}(1,-1), \\ 0 & \mbox{otherwise}. \end{array}\right. \end{equation}

The magnetic field is ${\boldsymbol B} = (\partial _y A, - \partial _x A)$. Thus,

(A12)\begin{equation} {\boldsymbol B} = \sum_n \sum_k ({\rm i} \alpha n - \beta k U't, - {\rm i} \beta k) c_{n,k} \exp ({\rm i}\alpha n y) \exp ({\rm i} \beta k [x - U'ty]) \exp [ \tau_{n,k}(t)]. \end{equation}

Then with the spatially averaged magnetic energy defined to be

(A13)\begin{equation} \langle {\boldsymbol B}(t)^2 \rangle = \int_0^{2 {\rm \pi}/\alpha} \int_0^{2 {\rm \pi}/\beta} | {\boldsymbol B}(t) |^2 \, {{\rm d}x} \, {{\rm d}y}, \end{equation}

we can use Parseval's theorem to deduce that

(A14)\begin{equation} \frac{\langle {\boldsymbol B}(t)^2 \rangle} { \langle {\boldsymbol B}(0)^2 \rangle} = \frac{\displaystyle\sum_n \sum_k [ (\alpha n + \beta k U't)^2 + \beta^2 k^2] |c_{n,k} |^2 \exp [2 \tau_{n,k} (t) ] } { \displaystyle\sum_n \sum_k [ (\alpha^2 n^2 + \beta^2 k^2] |c_{n,k} |^2 \exp [2 \tau_{n,k} (0) ]}. \end{equation}

References

Balbus, S.A. & Hawley, J.F. 1998 Instability, turbulence, and enhanced transport in accretion disks. Rev. Mod. Phys. 70 (1), 153.10.1103/RevModPhys.70.1CrossRefGoogle Scholar
Balman, Ş. & Revnivtsev, M. 2012 X-ray variations in the inner accretion flow of dwarf novae. Astron. Astrophys. 546, A112.CrossRefGoogle Scholar
Bath, G.T. & Pringle, J.E. 1981 The evolution of viscous discs. I. Mass transfer variations. Mon. Not. R. Astron. Soc. 194, 967986.10.1093/mnras/194.4.967CrossRefGoogle Scholar
Bodo, G., Cattaneo, F., Mignone, A. & Rossi, P. 2014 On the convergence of magnetorotational turbulence in stratified isothermal shearing boxes. Astrophys. J. 787 (1), L13.CrossRefGoogle Scholar
Brandenburg, A., Nordlund, A., Stein, R.F. & Torkelsson, U. 1995 Dynamo-generated turbulence and large-scale magnetic fields in a Keplerian shear flow. Astrophys. J. 446, 741.CrossRefGoogle Scholar
Brandenburg, A. & Subramanian, K. 2005 Astrophysical magnetic fields and nonlinear dynamo theory. Phys. Rep. 417 (1-4), 1209.CrossRefGoogle Scholar
Buat-Ménard, V., Hameury, J.M. & Lasota, J.P. 2001 The nature of dwarf nova outbursts. Astron. Astrophys. 366, 612622.CrossRefGoogle Scholar
Cannizzo, J.K. 1993 The accretion disk limit cycle model: toward an understanding of the long-term behavior of SS Cygni. Astrophys. J. 419, 318.CrossRefGoogle Scholar
Cannizzo, J.K. 2001 a The 2001 superoutburst of W.Z. Sagittae: a clue to the dynamics of accretion disks. Astrophys. J. 561 (2), L175L178.CrossRefGoogle Scholar
Cannizzo, J.K. 2001 b The V-EUV delay for dwarf nova outbursts: a case study for VW Hydri, U Geminorum, and SS Cygni. Astrophys. J. 556 (2), 847857.CrossRefGoogle Scholar
Coleman, M.S.B., Kotko, I., Blaes, O., Lasota, J. -P. & Hirose, S. 2016 Dwarf nova outbursts with magnetorotational turbulence. Mon. Not. R. Astron. Soc. 462, 37103726.CrossRefGoogle Scholar
Davidson, P.A. 2001 An Introduction to Magnetohydrodynamics. Cambridge University Press.CrossRefGoogle Scholar
Davis, S.W., Stone, J.M. & Pessah, M.E. 2010 Sustained magnetorotational turbulence in local simulations of stratified disks with zero net magnetic flux. Astrophys. J. 713 (1), 5265.CrossRefGoogle Scholar
Federrath, C. 2016 Magnetic field amplification in turbulent astrophysical plasmas. J. Plasma Phys. 82 (6), 535820601.CrossRefGoogle Scholar
Fleming, T.P., Stone, J.M. & Hawley, J.F. 2000 The effect of resistivity on the nonlinear stage of the magnetorotational instability in accretion disks. Astrophys. J. 530 (1), 464477.CrossRefGoogle Scholar
Frank, J., King, A. & Raine, D.J. 2002 Accretion Power in Astrophysics: Third Edition. Cambridge University Press.CrossRefGoogle Scholar
Fromang, S., Papaloizou, J., Lesur, G. & Heinemann, T. 2007 MHD simulations of the magnetorotational instability in a shearing box with zero net flux. II. The effect of transport coefficients. Astron. Astrophys. 476 (3), 11231132.CrossRefGoogle Scholar
Gammie, C.F. & Menou, K. 1998 On the origin of episodic accretion in dwarf novae. Astrophys. J. 492 (1), L75L78.CrossRefGoogle Scholar
Gogichaishvili, D., Mamatsashvili, G., Horton, W., Chagelishvili, G. & Bodo, G. 2017 Nonlinear transverse cascade and sustenance of MRI turbulence in keplerian disks with an azimuthal magnetic field. Astrophys. J. 845 (1), 70.CrossRefGoogle Scholar
Hawley, J.F., Gammie, C.F. & Balbus, S.A. 1995 Local three-dimensional magnetohydrodynamic simulations of accretion disks. Astrophys. J. 440, 742.10.1086/175311CrossRefGoogle Scholar
Hawley, J.F., Gammie, C.F. & Balbus, S.A. 1996 Local three-dimensional simulations of an accretion disk hydromagnetic dynamo. Astrophys. J. 464, 690.CrossRefGoogle Scholar
Heinemann, T. & Papaloizou, J.C.B. 2009 The excitation of spiral density waves through turbulent fluctuations in accretion discs. II. Numerical simulations with MRI-driven turbulence. Mon. Not. R. Astron. Soc. 397 (1), 6474.CrossRefGoogle Scholar
Johansen, A. & Levin, Y. 2008 High accretion rates in magnetised Keplerian discs mediated by a Parker instability driven dynamo. Astron. Astrophys. 490 (2), 501514.CrossRefGoogle Scholar
King, A.R., Pringle, J.E. & Livio, M. 2007 Accretion disc viscosity: how big is alpha? Mon. Not. R. Astron. Soc. 376, 17401746.CrossRefGoogle Scholar
Kotko, I. & Lasota, J.-P. 2012 The viscosity parameter $\alpha$ and the properties of accretion disc outbursts in close binaries. Astron. Astrophys. 545, A115.CrossRefGoogle Scholar
Livio, M. & Pringle, J.E. 1994 Star spots and the period gap in cataclysmic variables. Astrophys. J. 427, 956.CrossRefGoogle Scholar
Lynden-Bell, D. & Pringle, J.E. 1974 The evolution of viscous discs and the origin of the nebular variables. Mon. Not. R. Astron. Soc. 168, 603637.CrossRefGoogle Scholar
Mamatsashvili, G., Chagelishvili, G., Pessah, M.E., Stefani, F. & Bodo, G. 2020 Zero net flux MRI turbulence in disks: sustenance scheme and magnetic Prandtl number dependence. Astrophys. J. 904 (1), 47.CrossRefGoogle Scholar
Martin, R.G., Nixon, C.J., Pringle, J.E. & Livio, M. 2019 On the physical nature of accretion disc viscosity. New Astron. 70, 711.CrossRefGoogle Scholar
Martin-Alvarez, S., Katz, H., Sijacki, D., Devriendt, J. & Slyz, A. 2021 Unravelling the origin of magnetic fields in galaxies. Mon. Not. R. Astron. Soc. 504 (2), 25172534.10.1093/mnras/stab968CrossRefGoogle Scholar
Meyer, F. & Meyer-Hofmeister, E. 1983 Accretion disks in cataclysmic variables - the influence of the frictional parameter alpha on the structure. Astron. Astrophys. 128 (2), 420425.Google Scholar
Moffatt, H.K. 1978 Magnetic Field Generation in Electrically Conducting Fluids. Cambridge University Press.Google Scholar
Oyang, B., Jiang, Y.-F. & Blaes, O. 2021 Investigating lack of accretion disc eccentricity growth in a global 3D MHD simulation of a superhump system. Mon. Not. R. Astron. Soc. 505 (1), 117.CrossRefGoogle Scholar
Padoan, P., Federrath, C., Chabrier, G., Evans II, N.J., Johnstone, D., Jørgensen, J.K., McKee, C.F. & Nordlund, Å. 2014 The star formation rate of molecular clouds. In Protostars and Planets VI (eds. Beuther, H., Klessen, R.S., Dullemond, C.P. & Henning, Th.). University of Arizona Press. doi:10.2458/azu_uapress_9780816531240-ch004.Google Scholar
Parkin, E.R. & Bicknell, G.V. 2013 Global simulations of magnetorotational turbulence. I. Convergence and the quasi-steady state. Mon. Not. R. Astron. Soc. 435 (3), 22812298.CrossRefGoogle Scholar
Pjanka, P. & Stone, J.M. 2020 Stratified global MHD models of accretion disks in semidetached binaries. Astrophys. J. 904 (2), 90.CrossRefGoogle Scholar
Potter, W.J. & Balbus, S.A. 2014 An accretion disc instability induced by a temperature sensitive $\alpha$ parameter. Mon. Not. R. Astron. Soc. 441, 681689.CrossRefGoogle Scholar
Pringle, C.C.T., McMillan, B.F. & Teaca, B. 2017 A nonlinear approach to transition in subcritical plasmas with sheared flow. Phys. Plasmas 24 (12), 122307.CrossRefGoogle Scholar
Pringle, J.E. 1981 Accretion discs in astrophysics. Annu. Rev. Astron. Astrophys. 19, 137162.10.1146/annurev.aa.19.090181.001033CrossRefGoogle Scholar
Pringle, J.E., Verbunt, F. & Wade, R.A. 1986 Dwarf novae in outburst - modelling the observations. Mon. Not. R. Astron. Soc. 221, 169194.CrossRefGoogle Scholar
Pringle, J.E. & Wade, R.A. 1985 Interacting Binary Stars. Cambridge University Press.Google Scholar
Rembiasz, T., Obergaulinger, M., Cerdá-Durán, P., Aloy, M.-Á. & Müller, E. 2017 On the measurements of numerical viscosity and resistivity in Eulerian MHD codes. Astrophys. J. 230 (2), 18.CrossRefGoogle Scholar
Rincon, F. 2019 Dynamo theories. J. Plasma Phys. 85 (4), 205850401.10.1017/S0022377819000539CrossRefGoogle Scholar
Ryan, B.R., Gammie, C.F., Fromang, S. & Kestener, P. 2017 Resolution dependence of magnetorotational turbulence in the isothermal stratified shearing box. Astrophys. J. 840 (1), 6.CrossRefGoogle Scholar
Salvesen, G., Simon, J.B., Armitage, P.J. & Begelman, M.C. 2016 Accretion disc dynamo activity in local simulations spanning weak-to-strong net vertical magnetic flux regimes. Mon. Not. R. Astron. Soc. 457, 857874.CrossRefGoogle Scholar
Schekochihin, A.A. 2022 MHD turbulence: a biased review. J. Plasma Phys. 88 (5), 155880501.CrossRefGoogle Scholar
Schreiber, M.R., Hameury, J.M. & Lasota, J.P. 2003 Delays in dwarf novae I: the case of SS Cygni. Astron. Astrophys. 410, 239252.CrossRefGoogle Scholar
Schreiber, M.R., Hameury, J.M. & Lasota, J.P. 2004 Delays in dwarf novae II: VW Hyi, the tidal instability and enhanced mass transfer models. Astron. Astrophys. 427, 621635.10.1051/0004-6361:20041148CrossRefGoogle Scholar
Shakura, N.I. & Sunyaev, R.A. 1973 Black holes in binary systems. Observational appearance. Astron. Astrophys. 24, 337355.Google Scholar
Shakura, N.I. & Sunyaev, R.A. 1976 A theory of the instability of disk accretion on to black holes and the variability of binary X-ray sources, galactic nuclei and quasars. Mon. Not. R. Astron. Soc. 175, 613632.CrossRefGoogle Scholar
Shi, J.-M., Stone, J.M. & Huang, C.X. 2016 Saturation of the magnetorotational instability in the unstratified shearing box with zero net flux: convergence in taller boxes. Mon. Not. R. Astron. Soc. 456 (3), 22732289.CrossRefGoogle Scholar
Smak, J. 1999 Dwarf nova outbursts. III. The viscosity parameter alpha. Acta Astron. 49, 391401.Google Scholar
Smak, J.I. 1998 Dwarf nova outbursts. I. The UV delay. Acta Astron. 48, 677693.Google Scholar
Spitzer, L. 1962 Physics of Fully Ionized Gases. Interscience.Google Scholar
Stone, J.M., Hawley, J.F., Gammie, C.F. & Balbus, S.A. 1996 Three-dimensionalmagnetohydrodynamical simulations of vertically stratified accretion disks. Astrophys. J. 463, 656.CrossRefGoogle Scholar
Tobias, S.M. & Cattaneo, F. 2013 Shear-driven dynamo waves at high magnetic Reynolds number. Nature 497 (7450), 463465.CrossRefGoogle ScholarPubMed
Tout, C.A. & Pringle, J.E. 1992 Accretion disc viscosity - a simple model for a magnetic dynamo. Mon. Not. R. Astron. Soc. 259, 604612.CrossRefGoogle Scholar
Walker, J., Lesur, G. & Boldyrev, S. 2016 On the nature of magnetic turbulence in rotating, shearing flows. Mon. Not. R. Astron. Soc. 457 (1), L39L43.CrossRefGoogle Scholar
Warner, B. 1995 Cataclysmic Variable Stars. Cambridge University Press.CrossRefGoogle Scholar
Weiss, N.O. 1966 The expulsion of magnetic flux by Eddies. Proc. R. Soc. Lond. A 293 (1434), 310328.Google Scholar
Yeates, A.R. & Hornig, G. 2011 A generalized flux function for three-dimensional magnetic reconnection. Phys. Plasmas 18 (10), 102118.CrossRefGoogle Scholar
Figure 0

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.

Figure 1

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.

Figure 2

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)).

Figure 3

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.

Figure 4

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 5

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$.