1 Introduction
It is well known in the field of plasma physics – both laboratory and astrophysical – that many magnetic configurations are susceptible to instability. Indeed, ideal instabilities driven by current or pressure gradients may provide an operational limit for the magnetic configurations in plasma devices. Challenges in fusion research include both identifying stable plasma configurations and controlling the nonlinear development of instabilities.
In an astrophysical context, magnetic fields may also become dynamically unstable. One example is the loss of stability of magnetic fields in stellar atmospheres that may lead to huge releases of energy (and indeed flaring behaviour) (see e.g. Cowley et al. Reference Cowley, Wilson, Hurricane and Fong2003). The description of the loss of stability in the astrophysical context draws heavily on the pioneering theory from fusion plasmas such as Furth, Killeen & Rosenbluth (Reference Furth, Killeen and Rosenbluth1963), Taylor (Reference Taylor1968), Taylor & Newton (Reference Taylor and Newton2015).
A second example is the instability of magnetic configurations in stellar radiative zones (Tayler Reference Tayler1973) that potentially leads to the generation of turbulence there and even dynamo action (Dikpati & Gilman Reference Dikpati and Gilman2001). These instabilities are related to the current-driven instabilities of Tayler (Reference Tayler1973, Reference Tayler1980), Pitts & Tayler (Reference Pitts and Tayler1985). These authors examined the stability of toroidal magnetic fields to non-axisymmetric perturbations, both in cylindrical and spherical geometries (see also the extensive discussion in Spruit Reference Spruit1999). Current-driven instabilities use the magnetic field as their energy source, with a strong magnetic field required for the instability to proceed; the role of rotation is simply to mediate the rate at which energy can be extracted.
A related set of instabilities, which are most relevant to this report, have magnetic configurations that are stable in isolation but can be destabilized by the presence of a differential rotation (see e.g. Gilman & Fox Reference Gilman and Fox1997; Gilman & Dikpati Reference Gilman and Dikpati2002; Cally, Dikpati & Gilman Reference Cally, Dikpati and Gilman2003, Reference Cally, Dikpati and Gilman2008; Hollerbach & Cally Reference Hollerbach and Cally2009). These joint instabilities occur for relatively weak magnetic fields provided that the differential rotation is sufficiently strong. Here, the axisymmetric differential rotation and magnetic field, which in isolation are linearly stable, are together jointly unstable. The toroidal magnetic field therefore acts as a conduit to allow the extraction of energy from the differential rotation, though some energy may also be extracted from the current. In recent years, significant attention has focused on these instabilities, as it is believed that they may be important in the solar tachocline (Tobias Reference Tobias, Soward, Jones, Hughes and Weiss2005).
The tachocline is the thin layer of radial and latitudinal shear in the Sun at the base of the solar convection zone (Spiegel & Zahn Reference Spiegel and Zahn1992; Tobias Reference Tobias, Soward, Jones, Hughes and Weiss2005). Estimates place the tachocline at
$0.693\pm 0.002$
solar radii with a width of
$0.039\pm 0.013$
solar radii (Charbonneau et al.
Reference Charbonneau, Christensen-Dalsgaard, Henning, Larsen, Schou, Thompson and Tomczyk1999). Global magnetohydrodynamic (MHD) instabilities in this layer have been invoked in models to drive turbulence below the base of the convection zone and have also been used to provide an electromotive force that may contribute to the solar dynamo (Dikpati & Gilman Reference Dikpati and Gilman2001).
Investigations into these instabilities began following Watson’s linear study on purely hydrodynamic instabilities in two-dimensional shells with differential rotation, as is found in the tachocline layer. With a simple quadratic angular frequency, a two-dimensional shell was found to be hydrodynamically stable at parameter values appropriate to the tachocline (Watson Reference Watson1980). In a landmark paper, Gilman & Fox (Reference Gilman and Fox1997) extended Watson’s analysis by adding a toroidal field profile, stable in the absence of rotation, to the model. Their analysis determined that the MHD shell may be unstable to a non-axisymmetric joint instability if a large enough magnetic field is present. For plausible magnetic field configurations, even the differential rotation of the tachocline may provide a conduit to instability. It is also reasonable to expect that such instabilities arise in other solar-type stars with radiative interiors and convective outer envelopes. This instability can be found both for ideal MHD configurations as well as in models with diffusion, though in the latter case some forcing is required to maintain the basic state.
We comment here on the physics of the non-axisymmetric joint instability. The presence of the magnetic field allows the instability to proceed, and enables the subsequent extraction of energy from the differential rotation. This is clearly reminiscent of the physics of the magnetorotational instability (MRI) (Velikhov Reference Velikhov1959; Chandrasekhar Reference Chandrasekhar1960; Balbus & Hawley Reference Balbus and Hawley1991), where coupling induced by even a weak magnetic field can trigger an instability in a hydrodynamically stable configuration. However, the instability we consider here is more akin to the MRI with a toroidal magnetic field as studied by Ogilvie & Pringle (Reference Ogilvie and Pringle1996); in that case a non-axisymmetric instability arises that allows the system to access a turbulent lower energy state. As for that case, it is conjectured that, for the non-dissipative system, the joint instability acts for arbitrarily weak magnetic fields.
More recently, computational advances have made thorough nonlinear studies of MHD instabilities possible. Notably, Cally (Reference Cally2001) examined a variety of magnetic field profiles, identifying a clamshell instability in a two-dimensional model with a broad toroidal field. Using a three-dimensional thin-shell model with an additional poloidal field, Miesch (Reference Miesch2007) was able to reach a statistical steady state in which the clamshell instability could operate in a sustained manner. In this state, he observed quasiperiodic solutions in which energy was exchanged between magnetic, kinetic and potential energies in the nonlinear regime. Similarly, a ‘seasonal’ quasiperiodic exchange of energies resulting from the joint instability has recently been linked to solar observations by Dikpati et al. (Reference Dikpati, Cally, McIntosh and Heifetz2017).
In this paper, we describe striking and, to the best of our knowledge, novel behaviour of the joint instability in a two-dimensional nonlinear model without an imposed poloidal field. We explore the emergence of a multidecadal cycle in which energy is transformed from kinetic energy to magnetic potential energy, with sharp transitions between two distinct modes.
We also use this joint instability as a paradigm problem to investigate the effectiveness of a range of methods in describing the statistics of the nonlinear saturation in a turbulent state. In addition to linear stability analysis and fully nonlinear direct numerical simulation, a number of other methods for understanding the behaviour of a system of partial differential equations exist. In this paper, we will also apply generalized quasi-linear approximations, as well as a class of methods known as direct statistical simulation (DSS). DSS can be complementary to direct numerical simulation, as it solves for the evolution of the statistics of the system, rather than integrating the equations of motion themselves forward in time. This provides a number of advantages, both in terms of computational efficiency and physical insight. Tobias, Dagon & Marston (Reference Tobias, Dagon and Marston2011) and Constantinou & Parker (Reference Constantinou and Parker2018) fruitfully used DSS to study idealized, stochastically driven, astrophysical flows, calling for a more thorough characterization of DSS methods in these systems.
In the next section, we introduce the model of the two-dimensional instability. We proceed in § 3 to describe the range of tools we employ to investigate the nonlinear development of the instability. We describe the transitions that occur and how well these are captured in § 4, before drawing conclusions in § 5.
2 Description of the model
We consider the simplest possible model of the joint instability of differential rotation and magnetic fields. The evolution of the system, in the absence of forcing, is described by the MHD equations
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}\boldsymbol{u}}{\unicode[STIX]{x2202}t}+(\boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D735})\boldsymbol{u}+2\unicode[STIX]{x1D734}\times \boldsymbol{u}=-\unicode[STIX]{x1D735}p+\boldsymbol{j}\times \boldsymbol{B}+\unicode[STIX]{x1D708}\unicode[STIX]{x1D6FB}^{2}\boldsymbol{u}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}\boldsymbol{B}}{\unicode[STIX]{x2202}t}=\unicode[STIX]{x1D735}\times (\boldsymbol{u}\times \boldsymbol{B})+\unicode[STIX]{x1D702}\unicode[STIX]{x1D6FB}^{2}\boldsymbol{B}, & \displaystyle\end{eqnarray}$$
where
$\boldsymbol{u}$
is the velocity in the frame of the rotating shell,
$\unicode[STIX]{x1D734}$
is the angular rate of rotation of the shell,
$\unicode[STIX]{x1D708}$
is viscosity,
$\unicode[STIX]{x1D702}$
is magnetic diffusivity,
$\boldsymbol{B}$
is the magnetic field and
$\boldsymbol{j}=\unicode[STIX]{x1D735}\times \boldsymbol{B}$
is the current density. Here the density has been set to unity without loss of generality.
We consider evolution on a two-dimensional spherical surface so that the momentum and induction equation can be described by the evolution of two scalar equations, one for the absolute vorticity
$q(\unicode[STIX]{x1D703},\unicode[STIX]{x1D719},t)$
(with
$\unicode[STIX]{x1D703}$
being co-latitude and
$\unicode[STIX]{x1D719}$
longitude), and one for the scalar magnetic potential
$A(\unicode[STIX]{x1D703},\unicode[STIX]{x1D719},t)$
. A background rotation profile and toroidal magnetic field are imposed for all times and we study the behaviour of the fields relative to these prescribed profiles.
The background differential rotation profile can be decomposed into a solid body rotation at frequency
$\unicode[STIX]{x1D6FA}$
, and a latitudinal shear profile. Thus, in an absolute coordinate system, absolute vorticity
$q$
is the sum of the generated vorticity,
$\unicode[STIX]{x1D701}^{\prime }$
, the contribution from the shear profile,
$\unicode[STIX]{x1D701}_{0}$
, and the Coriolis parameter,
$f=2\unicode[STIX]{x1D6FA}\cos \unicode[STIX]{x1D703}$
.
$$\begin{eqnarray}q=\unicode[STIX]{x1D701}^{\prime }+\unicode[STIX]{x1D701}_{0}+2\unicode[STIX]{x1D6FA}\cos \unicode[STIX]{x1D703}.\end{eqnarray}$$
The total magnetic potential,
$A$
, can be similarly decomposed into generated magnetic potential,
$a$
, and the background,
$A_{0}$
$$\begin{eqnarray}A=a+A_{0}.\end{eqnarray}$$
The relative vorticity is defined as
$$\begin{eqnarray}\unicode[STIX]{x1D701}=\unicode[STIX]{x1D701}^{\prime }+\unicode[STIX]{x1D701}_{0}=\unicode[STIX]{x1D6FB}^{2}\unicode[STIX]{x1D713},\end{eqnarray}$$
which can also be represented in terms of the streamfunction
$\unicode[STIX]{x1D713}$
.
The evolution equations for
$q$
and
$A$
can then be written as
$$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}c@{}}{\displaystyle \frac{\unicode[STIX]{x2202}q}{\unicode[STIX]{x2202}t}}=J[q,\unicode[STIX]{x1D713}]+J[A,\unicode[STIX]{x1D6FB}^{2}A]-\unicode[STIX]{x1D705}\unicode[STIX]{x1D701}^{\prime }-\unicode[STIX]{x1D708}_{4}\unicode[STIX]{x1D6FB}^{6}(\unicode[STIX]{x1D6FB}^{2}+2)\unicode[STIX]{x1D701},\\ {\displaystyle \frac{\unicode[STIX]{x2202}A}{\unicode[STIX]{x2202}t}}=J[A,\unicode[STIX]{x1D713}]+\unicode[STIX]{x1D702}\unicode[STIX]{x1D6FB}^{2}a,\end{array}\right\} & & \displaystyle\end{eqnarray}$$
where
$J[A,B]$
gives the Jacobian on the unit sphere:
$$\begin{eqnarray}J[A,B]\equiv \frac{1}{\sin \unicode[STIX]{x1D703}}\left(\frac{\unicode[STIX]{x2202}A}{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}}\frac{\unicode[STIX]{x2202}B}{\unicode[STIX]{x2202}\unicode[STIX]{x1D703}}-\frac{\unicode[STIX]{x2202}A}{\unicode[STIX]{x2202}\unicode[STIX]{x1D703}}\frac{\unicode[STIX]{x2202}B}{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}}\right).\end{eqnarray}$$
Note that, in contrast to previous work such as Tobias et al. (Reference Tobias, Dagon and Marston2011), the Lorentz force term in the momentum equation includes the prescribed field. This modification is necessary given the form of background field used. Friction in the tachocline is parameterized with a frictional coefficient
$\unicode[STIX]{x1D705}$
. To remove enstrophy as it cascades to small scales, hyperviscosity
$\unicode[STIX]{x1D708}_{4}\unicode[STIX]{x1D6FB}^{6}(\unicode[STIX]{x1D6FB}^{2}+2)\unicode[STIX]{x1D701}$
is included in the linear operator of (2.6). The appearance of the operator
$(\unicode[STIX]{x1D6FB}^{2}+2)$
ensures that the hyperviscosity does not change the total angular momentum. The use of hyperviscosity in place of ordinary viscosity has been used to model dynamics in the low magnetic Prandtl number limit appropriate to stellar interiors (and for liquid metals) (Seshasayanan, Dallas & Alexakis Reference Seshasayanan, Dallas and Alexakis2017). Note however that we employ a regular diffusivity in the induction equation to maintain the balance between advection and diffusion there (Miesch et al.
Reference Miesch, Matthaeus, Brandenburg, Petrosyan, Pouquet, Cambon, Jenko, Uzdensky, Stone and Tobias2015).
We work on the unit sphere and in units of time such that
$\unicode[STIX]{x1D6FA}=2\unicode[STIX]{x03C0}$
. A ‘day’ is therefore a unit interval of time. The background vorticity field is chosen to have the form
$\unicode[STIX]{x1D701}_{0}(\unicode[STIX]{x1D703},\unicode[STIX]{x1D719})=\tilde{\unicode[STIX]{x1D701}}_{0}\,Y_{3}^{0}(\cos \unicode[STIX]{x1D703})$
and thus shears the flow symmetrically about the equator. The imposed background toroidal potential is taken to be
$A_{0}(\unicode[STIX]{x1D703},\unicode[STIX]{x1D719})=\tilde{A}_{0}\,Y_{2}^{0}(\cos \unicode[STIX]{x1D703})$
. In all the simulations reported here (except for the linear stability analysis of § 3.2) we set
$\tilde{\unicode[STIX]{x1D701}}_{0}=-3$
and
$\tilde{A}_{0}=1/2$
, in the joint instability regime where neither shear nor the background field by themselves suffice to trigger the instability. We set the friction coefficient to be
$\unicode[STIX]{x1D705}=0.05$
, equivalent to an e-folding time of 20 days. The hyperviscosity coefficient
$\unicode[STIX]{x1D708}_{4}$
is chosen such that the most rapidly dissipating mode decays at a rate of 1.
3 Methods
In this section we give a brief description of the methods we utilize in the paper to analyse the behaviour of the model, both in terms of the dynamics and statistics. These range from conservation laws, linear stability analysis, fully nonlinear simulation (DNS) in both spectral and real space, partially nonlinear simulation (both quasi-linear and generalized quasi-linear) and statistical simulation.
Our aim is to examine how well various approximations and statistical representations capture changes in the behaviour of the model – both in a dynamical and statistical sense.
3.1 Conservation laws
In the absence of damping and driving forces, the equations of motion (EOM) for the vorticity and magnetic potential conserve a number of linear and quadratic quantities. In the pure hydrodynamic case (
$A_{0}=0$
), kinetic energy, enstrophy and angular momentum are conserved. Higher-order Casimirs are also conserved in the continuum limit, but not in the finite-resolution numerical simulations that we employ. For
$A_{0}\neq 0$
the conserved quantities are angular momentum, total energy, mean square potential and cross-helicity, the latter given by the average over the sphere of the product of the absolute vorticity and the scalar magnetic potential,
$$\begin{eqnarray}\displaystyle \frac{1}{4\unicode[STIX]{x03C0}}\int \text{d}^{2}\unicode[STIX]{x1D6FA}qA. & & \displaystyle\end{eqnarray}$$
These invariants are respected by the quasi-linear and generalized quasi-linear approximations (Marston, Chini & Tobias Reference Marston, Chini and Tobias2016), as well as by the second-order cumulant expansion (Marston, Qi & Tobias Reference Marston, Qi and Tobias2019) that we shall utilize, and serve the practical purpose of validating code.
3.2 Linear stability analysis
The joint instability in a two-dimensional incompressible tachocline without viscous or ohmic dissipation was first described using a linear stability analysis (Gilman & Fox Reference Gilman and Fox1997). As described above, in our nonlinear equation solver, GCM,
$\unicode[STIX]{x1D701}_{0}$
and
$A_{0}$
are set to be proportional to the spherical harmonics
$Y_{3}^{0}$
and
$Y_{2}^{0}$
respectively, for the sake of simplicity. Thus,
$\unicode[STIX]{x1D701}_{0}$
in GCM differs slightly from that used by Cally and others (Gilman & Fox Reference Gilman and Fox1997; Cally Reference Cally2001). This is an acceptable substitution for our model, which is highly simplified and not physically precise. We have not been able to detect any qualitative difference in the simulations due to this modification. However, in order to confirm agreement between our work and earlier models, in this subsection we modify GCM to agree with the Gilman & Fox (Reference Gilman and Fox1997) rotation profiles to check that we can reproduce the linear stability analysis.
The equations of motion described in § 2 are rewritten in terms of equilibrium state streamfunctions in an absolute coordinate system. We first address the ideal case, in which friction, magnetic diffusivity and hyperviscosity are set to zero. We group together the Coriolis parameter,
$f$
, and the forcing function,
$\unicode[STIX]{x1D701}_{0}$
, and represent them using a single streamfunction
$\unicode[STIX]{x1D713}_{0}$
.
$$\begin{eqnarray}q=\unicode[STIX]{x1D701}^{\prime }+(\unicode[STIX]{x1D701}_{0}+f)=\unicode[STIX]{x1D6FB}^{2}\unicode[STIX]{x1D713}^{\prime }+\unicode[STIX]{x1D6FB}^{2}\unicode[STIX]{x1D713}_{0}.\end{eqnarray}$$
We define
$\unicode[STIX]{x1D707}=\cos \unicode[STIX]{x1D703}$
, and, assuming our prescribed fields are zonally symmetric, define a rotational angular frequency and Alfvén frequency for the prescribed fields. We set the frequencies to be the profiles given in Gilman & Fox (Reference Gilman and Fox1997).
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D714}_{0}=\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}_{0}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D707}}=\unicode[STIX]{x1D6FA}-\tilde{\unicode[STIX]{x1D701}}_{0}\unicode[STIX]{x1D707}^{2} & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D6FC}_{0}=\frac{\unicode[STIX]{x2202}A_{0}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D707}}=\tilde{A}_{0}\unicode[STIX]{x1D707}. & \displaystyle\end{eqnarray}$$
The equations are linearized in the perturbations
$\unicode[STIX]{x1D713}^{\prime }$
and
$a$
, with each field expanded in terms of associated Legendre polynomials as
$\unicode[STIX]{x1D713}^{\prime }=\sum _{\ell =m}^{\infty }\unicode[STIX]{x1D6F9}_{\ell }^{m}P_{\ell }^{m}(\unicode[STIX]{x1D707})\text{e}^{\text{i}(m\unicode[STIX]{x1D719}-\unicode[STIX]{x1D714}t)}$
and
$a=\sum _{\ell =m}^{\infty }A_{l}^{m}P_{\ell }^{m}(\unicode[STIX]{x1D707})\text{e}^{\text{i}(m\unicode[STIX]{x1D719}-\unicode[STIX]{x1D714}t)}$
where
$\unicode[STIX]{x1D714}$
is the complex-valued angular frequency eigenvalue. The substitutions yield the equations:
$$\begin{eqnarray}(m\unicode[STIX]{x1D714}_{0}-\unicode[STIX]{x1D714})\mathop{\sum }_{\ell =m}^{\infty }A_{\ell }^{m}P_{\ell }^{m}=m\unicode[STIX]{x1D6FC}_{0}\mathop{\sum }_{\ell =m}^{\infty }\unicode[STIX]{x1D6F9}_{\ell }^{m}P_{\ell }^{m}(\unicode[STIX]{x1D707}),\end{eqnarray}$$
and
$$\begin{eqnarray}\displaystyle & & \displaystyle -(m\unicode[STIX]{x1D714}_{0}-\unicode[STIX]{x1D714})\mathop{\sum }_{\ell =m}^{\infty }\ell (\ell +1)\unicode[STIX]{x1D6F9}_{\ell }^{m}P_{\ell }^{m}-m\frac{\text{d}^{2}}{\text{d}\unicode[STIX]{x1D707}^{2}}[(1-\unicode[STIX]{x1D707}^{2})\unicode[STIX]{x1D714}_{0}]\mathop{\sum }_{\ell =m}^{\infty }\unicode[STIX]{x1D6F9}_{\ell }^{m}P_{\ell }^{m}\nonumber\\ \displaystyle & & \displaystyle \quad +\,m\unicode[STIX]{x1D6FC}_{0}\mathop{\sum }_{\ell =m}^{\infty }\ell (\ell +1)A_{\ell }^{m}P_{\ell }^{m}+m\frac{\text{d}^{2}}{\text{d}\unicode[STIX]{x1D707}^{2}}[(1-\unicode[STIX]{x1D707}^{2})\unicode[STIX]{x1D6FC}_{0}]\mathop{\sum }_{\ell =m}^{\infty }A_{\ell }^{m}P_{\ell }^{m}=0.\end{eqnarray}$$
Since the selected profiles for rotation and magnetic field are symmetric and antisymmetric about the equator respectively, either
$\unicode[STIX]{x1D713}^{\prime }$
is symmetric and
$a$
anti-symmetric, or
$a$
is symmetric and
$\unicode[STIX]{x1D713}^{\prime }$
antisymmetric. By truncating
$\ell$
at a finite number
$L_{\text{max}}$
, two
$2L_{\text{max}}\times 2L_{\text{max}}$
matrices, one for the symmetric case and one for the antisymmetric case, can then be diagonalized numerically to find the growth rates of the instabilities.
Using this method, the existence of the joint instability is confirmed in the dissipationless regime. The system continues to exhibit an instability if dissipation is added to our analysis (non-zero friction and/or magnetic diffusivity), provided the dissipation is not too strong. As expected, increasing dissipation decreases the growth rate of the instability (Dikpati, Cally & Gilman Reference Dikpati, Cally and Gilman2004). For all values of parameters examined,
$m=1$
is the most unstable mode, as found previously (Gilman & Fox Reference Gilman and Fox1997; Cally Reference Cally2001). Additionally, GCM with appropriately modified profiles is found via time stepping to reproduce the linear analysis both when nonlinear terms are turned off, and in the limit of small amplitude waves with nonlinear terms included (figure 1).
Figure 1. Comparison between GCM and linear analysis for joint instability parameters with no dissipation. (a) Power in
$\ell =4$
,
$m=1$
spherical harmonic over time in GCM code using prescribed field profiles from Gilman & Fox (Reference Gilman and Fox1997) without nonlinear terms. An exponential fit with a growth rate of 0.027 is shown, as predicted by linear analysis. Inset shows same data on a semilog plot. (b) Power in
$\ell =4$
,
$m=1$
mode over time in GCM, including nonlinear terms. Bolded black points show exponential fit with a growth rate of 0.027 for times before 100 days.
3.3 Direct numerical simulation in real and spectral space
Fully nonlinear (NL) spectral DNS with truncation
$0\leqslant \ell \leqslant L$
and
$|m|\leqslant \min \{\ell ,M\}$
is performed. We set spectral cutoffs
$L=60$
and
$M=15$
and verify, by comparison with high-resolution simulations carried out in real space, that these cutoffs suffice. The pure spectral version of (2.6) is integrated forward in time using a fourth-order accurate Runge–Kutta algorithm with an adaptive time step
$\unicode[STIX]{x0394}t$
. Each time step requires
$O(L^{3}M^{2})$
floating point computations that at high resolutions would be prohibitively expensive compared to a pseudo-spectral algorithm but is feasible here for the moderate resolutions that we study. The computation is made faster by skipping over triads that vanish due to symmetry.
To verify that the full spectral simulation has sufficient resolution, finer-scale DNS of the fluid is also performed in real space on a spherical geodesic grid (Dritschel, Qi & Marston Reference Dritschel, Qi and Marston2015) of
$D=163\,842$
cells; the lattice operators conserve energy and enstrophy. The vorticity evolves forward in time by a second-order accurate leapfrog algorithm, with a Robert–Asselin–Williams filter of 0.001 and
$\unicode[STIX]{x1D6FC}=0.53$
(Williams Reference Williams2009). The time step is fixed at
$\unicode[STIX]{x0394}t=0.003$
.
3.4 Quasi-linear and generalized quasi-linear numerical simulation
The quasi-linear (QL) approximation has its historical origins in the work of Malkus (Reference Malkus1954) and Herring (Reference Herring1963), with perhaps the clearest earliest exposition given in the plasma context by Vedenov, Velikhov & Sagdeev (Reference Vedenov, Velikhov and Sagdeev1962). The approximation is a self-consistent mean-field theory that in this context retains the scattering of eddies (defined as perturbations with respect to the zonal mean) off the mean zonal flow, and processes in which two eddies of equal and opposite zonal wavenumber combine to influence the mean flow. The two retained triadic interactions are depicted in figure 2. Note that the general eddy
$+$
eddy
$\rightarrow$
eddy scattering is dropped.
Figure 2. Triadic interactions retained in the QL approximation. The zonal wavenumber of the eddies is labelled by
$m$
; the zonal mean flow has
$m=0$
and is depicted by the solid line. (a) An eddy scatters off the zonal mean flow. (b) Two eddies of equal but opposite zonal wavenumber combine to modify the zonal mean flow.
A generalization of QL called the generalized quasi-linear approximation (GQL) (Marston et al.
Reference Marston, Chini and Tobias2016) systematically improves upon the QL approximation by including, self-consistently, fully nonlinear interactions among large-scale modes. The zonal mean of the QL approximation is generalized to include modes with low zonal wavenumber at or below a cutoff
$|m|\leqslant \unicode[STIX]{x1D6EC}$
; see figure 3.
$\unicode[STIX]{x1D6EC}=0$
is QL, while the limit
$\unicode[STIX]{x1D6EC}=M$
recovers the full nonlinear EOMs; thus GQL interpolates between QL and NL. Marston et al. (Reference Marston, Chini and Tobias2016) examined GQL in the context of two-dimensional barotropic turbulence on a spherical surface and on a
$\unicode[STIX]{x1D6FD}$
-plane and showed it to be more effective than the QL approximation in reproducing both the dynamics and statistics of these flows away from equilibrium. Remarkably, this remained true even if only a single extra mode was retained in the large scales. Subsequent work demonstrated the effectiveness of the GQL approximation for the case of axisymmetric (two-dimensional) helical magnetorotational instability (HMRI) at low magnetic Reynolds number (Child et al.
Reference Child, Hollerbach, Marston and Tobias2016), and in three-dimensional rotating plane Couette flow (Tobias & Marston Reference Tobias and Marston2016).
Figure 3. Triadic interactions retained and omitted in the GQL approximation. Modes with low zonal wavenumbers
$|m|\leqslant \unicode[STIX]{x1D6EC}$
interact fully nonlinearly (b). The interactions of the low modes with the high wavenumber modes, (a,c), generalizes the quasi-linear interactions shown in figure 2.
3.5 Direct statistical simulations: CE2
3.5.1 CE2 simulations
In order to understand better the physical mechanism behind the dynamics and statistical properties of the model, a direct statistical simulation (DSS) technique based upon expansion of the equations of motion for the equal-time cumulants is also used. We retain only the first and second cumulants (CE2); for a review, see Marston et al. (Reference Marston, Qi and Tobias2019). In the CE2 approximation, the contribution of the third cumulant to the tendency of the second cumulant is neglected, an approximation that is mathematically equivalent to the QL approximation of eliminating the triadic interaction between two eddies that results in a third eddy. The CE2 equations are an exact closure within the QL approximation. DSS is best suited to describing systems with large-scale inhomogeneous and anisotropic flows. The approach was applied to a different, stochastically forced, tachocline model in Tobias et al. (Reference Tobias, Dagon and Marston2011). Here we apply it to the deterministic Cally model.
A program that implements spectral DNS, real-space DNS, QL, GQL and DSS, and also includes all the graphical tools needed to visualize statistics, is freely available.Footnote 1 The Objective-C++ and Swift programming languages are employed. C blocks and Grand Central Dispatch enable the efficient use of multiple CPU cores.
4 Results
We consider the evolution of the instability and determine how the dynamics changes as the magnetic Reynolds number is increased, for fixed profiles of differential rotation and magnetic field and fixed relaxation time. The magnetic Reynolds number is defined as
$$\begin{eqnarray}R_{m}={\displaystyle \frac{UR}{\unicode[STIX]{x1D702}}},\end{eqnarray}$$
where
$U$
is the zonal velocity of the basic state at the equator,
$R$
is the radius and
$\unicode[STIX]{x1D702}$
is the magnetic diffusivity. The results are summarized in the bifurcation diagram of figure 4, where we use the kinetic energy as a fiducial measure of the amplitude of the solution. We note that this diagram has been flipped such that solutions with the kinetic energy of the basic differential rotation appear at the bottom, so that the figure resembles a traditional bifurcation diagram. We begin by describing the transitions in behaviour as
$R_{m}$
is increased as determined by the fully nonlinear simulations. We will then go on to determine how well these transitions are captured by the various approximations (CE2, GQL
$\unicode[STIX]{x1D6EC}=1$
and GQL
$\unicode[STIX]{x1D6EC}=3$
) subsequently.
Figure 4. Bifurcation diagram showing the nature of the solutions obtained by the various methods as a function of magnetic Reynolds number
$R_{m}$
. Note the kinetic energy axis has been inverted so that solutions with the kinetic energy of the basic state appear at the bottom of the diagram. The data have been shifted slightly in
$R_{m}$
for clarity.
4.1 Fully nonlinear dynamics
For low
$R_{m}$
(
${<}90$
), the basic state is stable, as described by the linear theory above. For
$R_{m}>90$
, NL solutions (denoted in purple) are initially located on a primary branch, where the energy in the perturbation is weak and contained almost completely in the
$m=1$
mode, as shown by the power spectrum in figure 5(a). As
$R_{m}$
is increased, the solution remains on this primary branch, with the solution having a characteristic
$m=1$
dependence on longitude and a banded dependence on latitude (as shown in figure 5
b); this solution does have power in higher azimuthal wavenumbers, though these are limited to
$m\lesssim 3$
as shown in figure 5(c).
Figure 5. (a) Power spectrum for
$R_{m}=90$
– energy is largely restricted to the
$m=1$
mode. Here
$l$
is the spherical wavenumber and
$m$
is the zonal wavenumber. (b) Snapshot of the relative vorticity and (c) power spectrum for
$R_{m}\sim 300$
. For this value of
$R_{m}$
, more azimuthal structure is observed.
At
$R_{m}\sim 700$
, our NL simulation is able to find two different types of solution depending on initial conditions and bistability is present. Figure 6 shows the timelines of the zonally averaged relative vorticity for two different initial conditions. It can be seen clearly that there are two different states, one of which, shown in figure 6(a), has a high relative vorticity (as for those at lower
$R_{m}$
) and one (shown in figure 6
b) yields a state with much lower vorticity.
Figure 6. Hysteresis in the fully nonlinear simulations at an intermediate value of the magnetic Reynolds number. Timelines of the zonally averaged relative vorticity for
$R_{m}\sim 700$
show (a) the high vorticity solution and (b) the low vorticity solution.
For the highest value of
$R_{m}$
considered (
$R_{m}=2800$
), interesting relaxation dynamics is found, as shown in figures 7, 8 and 9. From these figures it is clear that the relaxation occurs as an oscillation between the high and low vorticity states described in the hysteretic scenario above. Both of these states are now unstable at this high
$R_{m}$
(presumably having undergone Hopf bifurcations) and the large amplitude relaxation oscillation takes the form of a near-heteroclinic solution oscillating between the two unstable states. The transitions between these two states occur rapidly, on a time scale of approximately 300 days, with a cycle period of approximately 9000 days. Relative vorticity in the high vorticity state is approximately zonally symmetric and power is distributed among a broad range of spherical harmonics. In the low vorticity state, the relative vorticity snapshot shows fine structure and very little power in relative vorticity spherical harmonics (see figure 9).
Figure 7. Abrupt transitions in relative vorticity (
$\unicode[STIX]{x1D701}$
) for a spectral simulation with
$L_{\text{max}}=60$
,
$M_{\text{max}}=15$
and
$R_{m}=2800$
. (a) Fraction of total energy in the form of kinetic energy versus time. (b) Zonally averaged relative vorticity for the same time period. The kinetic energy share tracks the transitions between the two distinct states and highlights the rapidity of the switch.
Figure 8. Abrupt transitions in generated magnetic potential,
$a$
, for a spectral simulation with
$L_{\text{max}}=60$
,
$M_{\text{max}}=15$
and
$R_{m}=2800$
. (a) Fraction of total energy in the form of magnetic potential energy versus time. (b) Zonally averaged generated magnetic potential for the same time period.
Figure 9. Snapshots of both high and low vorticity states for
$R_{m}=2800$
. In the top timeline, the double bars at 18 000–19 000 and 24 000–25 000 days show the period over which time averaging is performed for figure 14. The two columns of figures are representative snapshots taken within these two windows respectively (not time averaged). (a,d) Power spectra of vorticity modes. (b,e) Relative vorticity on a cylindrical projection. (c,f) Second vorticity cumulant on a cylindrical projection.
We note that this behaviour is robust. It appears in both spectral and geodesic simulations. The transitions do not change qualitatively as resolution in real space is increased, and, as seen in figure 10, the transitions are reproducible using both a fully spectral simulation in a basis of spherical harmonics with
$L_{\text{max}}=60$
and
$M_{\text{max}}=15$
, as well as a real space simulation using six (pictured) and seven geodesics. The transitions are consistent and regular even in our longest, highest resolution runs. The addition of small stochastic forcing of the vorticity also does not alter the behaviour, demonstrating that the transitions are not triggered by noise.
Perhaps it is most informative to describe this relaxation oscillation in the ‘potential energy–kinetic energy’ phase plane. In figure 11, the abrupt transitions appear as a near-elliptical path. The low vorticity state occurs along the
$y$
-axis, and the high vorticity state occurs along the rest of the path. In one period of the relaxation oscillation, the system moves slowly down the
$y$
-axis and then much more quickly along the rest of the ellipse, such that the time spent in both states is comparable, as we see in figure 7.
The physics of this relaxation oscillation can be elucidated by performing additional numerical simulations for varying applied magnetic field strengths
$B$
(at fixed
$R_{m}$
). Figure 12 shows both the birth and death of the relaxation oscillation as
$B$
is varied. For small
$B=0.5$
or 0.9, the solution takes the form of a high vorticity state (with large kinetic energy and small kinetic energy), whilst for high
$B=2$
it takes the form of a low vorticity state. For intermediate values of
$B$
the solutions has the character of a relaxation oscillation. These additional simulations also show that the period of the oscillation is controlled by the near-heteroclinic nature of the solution – as
$B$
is increased the period increases as the solution spends increasing amount of time near the invariant low vorticity state. This behaviour is reminiscent of other ‘limit cycle oscillations’ (or LCOs) in plasma such as edge harmonic oscillation (EHO) that is present owing to the action of an instability (the kink-peeling instability) in the presence of a sheared rotation (in this case the
$\boldsymbol{E}\times \boldsymbol{B}$
rotation), see, for example, Wilks et al. (Reference Wilks, Garofalo, Diamond, Guo, Hughes, Burrell and Chen2018).
Figure 10. A real space simulation using the same parameters as the spectral space simulation of relative vorticity in figure 7 (
$R_{m}=2800$
). The behaviour is qualitatively similar.
Figure 11. Abrupt transitions in kinetic energy–potential energy space for
$R_{m}=2800$
. (a) Real space simulations. (b) Spectral space simulations. Arrows show the direction of time. Both simulations reach a steady cycle after a brief initialization period.
Figure 12. Abrupt transitions in kinetic energy–potential energy space for
$R_{m}=2800$
and varying applied field
$B$
. (a) Kinetic energy versus time. (b) Potential energy versus time.
4.2 CE2 and GQL
We now turn to how well the DSS method at CE2 and the DNS methods utilizing the GQL approximations with various severe truncations (
$\unicode[STIX]{x1D6EC}=1,3$
) perform. Figure 4 gives the average kinetic energy of the solutions found by these methods as a function of
$R_{m}$
. We note that the DSS and GQL runs were performed at the same
$R_{m}$
as the corresponding NL solution though the results have been shifted slightly in the bifurcation diagram for ease of viewing.
We begin by examining the efficacy of the most severe quasi-linear approximation CE2. Figure 4 shows that CE2 captures the linear instability (as it must) and the amplitude of the primary branch, with the orange circles tracking the purple circles for
$R_{m}<200$
. This is encouraging for quasi-linear theory. However by
$R_{m}\sim 300$
discrepancies occur between CE2 and NL. Here CE2 finds an oscillatory low vorticity state (whereas NL found a steady high vorticity state as described above). Interestingly, subsequent increase of
$R_{m}$
to
${\sim}700$
increases the accuracy of CE2. It is remarkable that it is able to reproduce the hysteresis found by NL (though it finds an oscillatory rather than a steady low vorticity state). However, perhaps unsurprisingly, for the highest
$R_{m}$
considered CE2 is unable to reproduce the relaxation oscillations shown by NL (figure 13). When initialized with sufficient magnetic potential energy, CE2 passes close to a state with very low vorticity and then makes a rapid transition to an periodic solution that is statistically similar to the high vorticity mode seen in NL trials (figure 13
a). After this first transition, the system remains in this state. When only the vorticity field is perturbed during initialization, CE2 immediately reaches the approximate fixed point (figure 13
b). Interestingly, this transition breaks north–south symmetry imposed at initialization when it reaches the high vorticity mode.
Time averaging is performed in both the low vorticity transient state and for the final oscillatory state. The zonally averaged fields are compared with time averaged results from spectral NL for high and low vorticity states, as shown in figure 14. The NL low vorticity state and the CE2 low vorticity state have statistically similar field profiles, while the NL high vorticity state and the CE2 high vorticity fixed point also have statistically similar field profiles. Thus, CE2 is capable of accessing both of the states that the NL simulation transits between, but is incapable of reproducing the periodic transitions between them.
This interpretation agrees well with the time series viewed in energy space. In contrast to figure 11, a cycle does not form. In figure 15(a), the system seems to reach the same low vorticity state as the nonlinear simulation during its transient, but then becomes stuck in a small region of energy space that has a higher potential energy than any of the points in the nonlinear elliptical cycle. In figure 15(b), the system immediately finds itself stuck in this same region of high potential and kinetic energy, and performs a small cycle there.
Figure 13. Zonally averaged relative vorticity versus time for CE2 simulations at
$R_{m}=2800$
. (a) CE2 initialized with perturbations in magnetic potential. (b) CE2 initialized with perturbations in vorticity. Time averaging is performed for the 1000 day period between the double bars in both figures for figure 14.
Figure 14. Comparison between time averaged field profiles for
$R_{m}=2800$
. Thick lines correspond to fully nonlinear spectral DNS, and thin lines correspond to CE2. The thick solid line is the first time averaged period (high vorticity) in figure 9, and the thick dashed line is the second period (low vorticity) in figure 9. The thin dashed line is the time averaged portion of figure 13(a), and the thin solid line is the time averaged portion of figure 13(b). Each panel compares these four time averaged and zonally averaged fields, with (a) comparing relative vorticity, (b) comparing zonal velocity (not including solid body rotation), (c) comparing generated magnetic fields and (d) comparing toroidal fields. Fields are plotted as a function of latitude.
Figure 15. CE2 simulations in kinetic energy–potential energy space for
$R_{m}=2800$
. Arrows show the direction of time, and the insets show the long-time behaviour. (a) CE2 initialized with perturbations in magnetic potential, as in figure 13(a). (b) CE2 initialized with perturbations in vorticity, as in figure 13(b).
The inability to transition back to the high magnetic potential state in CE2 shows that these transitions are dependent on processes that cannot be captured by a quasi-linear approximation. It may be that they owe their origin to eddy–eddy
$\rightarrow$
eddy interactions, since these interactions are forbidden in CE2. The single transition that occurs may be the result of lingering initialization energy, likely aided by the fact that the high vorticity mode’s power is concentrated in low zonal modes, while the low vorticity mode’s power is more spread out, making a transition to the high vorticity mode possible without eddy–eddy interactions. We may hope that the problem is remedied by employing a GQL approximation.
Figure 16. GQL simulations at
$R_{m}=2800$
. (a) GQL simulation with
$\unicode[STIX]{x1D6EC}=1$
. (b) GQL simulation with
$\unicode[STIX]{x1D6EC}=3$
. Resolution remains
$L_{\text{max}}=60$
and
$M_{\text{max}}=15$
.
The bifurcation diagram in figure 4 shows that GQL with
$\unicode[STIX]{x1D6EC}=1$
or
$\unicode[STIX]{x1D6EC}=3$
indeed performs better than CE2 in reproducing the NL results for all
$R_{m}$
. Figures 16 and 17 show the relative vorticity timelines and the corresponding phase plane analysis for both of these thresholds
$\unicode[STIX]{x1D6EC}$
at
$R_{m}=2800$
. These timelines illustrate that GQL with
$\unicode[STIX]{x1D6EC}=1$
performs similarly to CE2 at this choice of parameter, yielding a weakly oscillatory state. However, increasing the threshold to
$\unicode[STIX]{x1D6EC}=3$
allows the GQL approximation to reproduce the relaxation oscillation. Remarkably, GQL reproduces both the amplitude and period of the oscillation at this threshold, which is a testament to the effectiveness of this approximation.
5 Conclusions
We have shown that a simple, two-dimensional tachocline model can reproduce periodic transitions qualitatively similar to the seasonal exchange between toroidal field and wave-like dynamics observed in the solar record. The model undergoes complicated nonlinear transitions as the magnetic Reynolds number,
$R_{m}$
, is increased, with both hysteresis and relaxation oscillations emerging naturally. This near-heteroclinic behaviour has, to the authors’ knowledge, not yet been seen in joint instability simulations, though oscillations have also been seen in sustained magnetoshear instabilities (Miesch Reference Miesch2007) and quasiperiodic energy exchange between toroidal field and magnetic Rossby waves have also recently been observed (Dikpati et al.
Reference Dikpati, Cally, McIntosh and Heifetz2017) both in models and in solar observations. Our simple model is less expensive than three-dimensional thin-shell models or shallow water models that are often studied (Gilman & Dikpati Reference Gilman and Dikpati2002; Miesch Reference Miesch2007), which allows us to examine long time scale behaviour at the cost of restricting instabilities to two dimensions.
We have also assessed the efficacy of the quasi-linear (QL) and generalized quasi-linear (GQL) approximations in reproducing this dynamics and that of direct statistical simulation (DSS at CE2) in reproducing the statistics. We find that CE2 is capable of reproducing many of the nonlinear states (and indeed the hysteresis between states), but for the most complicated relaxation oscillation states, CE2 cannot reproduce the transitions between states. This is another example of the limitation of quasi-linear models, that may become less relevant as the system moves away from statistical equilibrium (Tobias & Marston Reference Tobias and Marston2013). We speculate that the addition of noise to DSS might enable this statistical formalism to reproduce the transitions. On the other hand, GQL appears to contain enough nonlinearity to represent the complicated nonlinear relaxation oscillation. This result is encouraging for a program of direct statistical simulation based on GQL or upon DSS formulations predicated on more sophisticated averaging procedures than simple zonal averaging (Bakas & Ioannou Reference Bakas and Ioannou2013, Reference Bakas and Ioannou2014; Allawala, Tobias & Marston Reference Allawala, Tobias and Marston2017).
Acknowledgements
A.P. was supported in part by a Brown University LINK award. S.M.T. was supported partially by a Leverhulme Fellowship and partially by funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement no. D5S-DLV-786780).
