1 Introduction
The magnetorotational instability (MRI) (Velikhov Reference Velikhov1959) is one of the most important processes in modern astrophysics, thought to be responsible for the outward transport of angular momentum in magnetised accretion disks (Balbus & Hawley Reference Balbus and Hawley1991), and has been the subject of numerous publications in recent years (see the review by Julien & Knobloch (Reference Julien and Knobloch2010)). It is driven by the interaction of large-scale Keplerian rotation and mean magnetic fields, and the nonlinear saturation of the instability can lead to turbulent dynamics and efficient angular momentum transport. The presence of an axial magnetic field enables the instability to proceed for fluids with angular momentum profiles that increase with radius; for the same configurations hydrodynamic flows are stable via the Rayleigh criterion.
Direct numerical simulations of this instability (Bai & Stone Reference Bai and Stone2014; Suzuki & Inutsuka Reference Suzuki and Inutsuka2014; Gressel & Pessah Reference Gressel and Pessah2015; Meheut et al. Reference Meheut, Fromang, Lesur, Joos and Longaretti2015), although illuminating for the underlying physics, are performed at parameter values that are far from those pertaining to astrophysical disks. For this reason, approaches other than direct solution of the equations have been brought to bear on this problem. The first of these is the utilisation of laboratory experiments in parameter regimes that may not be accessible to numerical simulations. This approach has revealed much about the interaction of rotating fluids, turbulence and magnetic fields as described below, but its direct applicability to the problem of angular momentum transport in disks remains problematic. The second strategy is to utilise direct statistical simulation (DSS) for astrophysical flows (Tobias, Dagon & Marston Reference Tobias, Dagon and Marston2011), where the statistics of such flows are obtained directly rather than emerging from averaging of the dynamics from DNS.
As noted above, an important tool to understand such an instability is laboratory experiments. The first of these, by Ji, Goodman & Kageyama (Reference Ji, Goodman and Kageyama2001) and Rüdiger & Zhang (Reference Rüdiger and Zhang2001) simultaneously, but independently, proposed an experimental set-up with a liquid metal in a cylindrical Taylor–Couette geometry with an imposed axial magnetic field. Instability occurs when the non-dimensional magnetic Reynolds number
$Rm={\it\Omega}_{i}r_{i}^{2}/{\it\eta}\geqslant O(10)$
. Together with the fact that liquid metals have magnetic Prandtl number
$Pm\sim O(10^{-6})$
, this translates into a hydrodynamic Reynolds number
$Re=Rm/Pm\geqslant O(10^{7})$
. Such large values present severe technical difficulties (Schartman, Ji & Burin Reference Schartman, Ji and Burin2009), especially when it comes to conditions at the end plates, which have a strong tendency to dominate the flow (Hollerbach & Fournier Reference Hollerbach and Fournier2004; Avila Reference Avila2012; Gissinger, Goodman & Ji Reference Gissinger, Goodman and Ji2012). As a result, this ‘classical’ set-up has not yet succeeded in obtaining the standard MRI (SMRI) – although the Princeton MRI experiment has yielded many other interesting results (Nornberg et al.
Reference Nornberg, Ji, Schartman, Roach and Goodman2010; Roach et al.
Reference Roach, Spence, Gissinger, Edlund, Sloboda, Goodman and Ji2012). Very recently it has also been suggested (Flanagan et al.
Reference Flanagan, Clark, Collins, Cooper, Khalzov, Wallace and Forest2015) that the SMRI could be obtained using plasma rather than liquid metal.
Alternatively, Hollerbach & Rüdiger (Reference Hollerbach and Rüdiger2005) suggested adding an azimuthal field
$B_{{\it\phi}}\sim r^{-1}$
, imposed by an electric current running along the central axis of the cylinders. This has a surprisingly large influence on the onset of instability, switching the criterion from
$Rm\geqslant O(10)$
to
$Re\geqslant O(10^{3})$
– that is, reducing the critical Reynolds number by some four orders of magnitude. Despite this dramatic reduction, and indeed fundamentally different scaling with
$Pm$
, this new instability, now known as the helical MRI (HMRI), is continuously connected to the SMRI (although a number of subtleties are involved (Kirillov & Stefani Reference Kirillov and Stefani2010; Kirillov, Stefani & Fukumoto Reference Kirillov, Stefani and Fukumoto2014)). Following its theoretical prediction, the HMRI was quickly observed in a series of experiments in the PROMISE facility (see Stefani et al. (Reference Stefani, Gerbeth, Gundrum, Hollerbach, Priede, Rüdiger and Szklarski2009) and earlier references therein). The HMRI also continues to be the subject of further theoretical analysis (Liu et al.
Reference Liu, Goodman, Herron and Ji2006; Priede, Grants & Gerbeth Reference Priede, Grants and Gerbeth2007; Priede Reference Priede2011).
The main purpose of this paper, however, is to determine the regime of applicability of the GQL approximation for a system of wall bounded, electrically conducting flows. The HMRI is an excellent system to elucidate the strengths and weaknesses of this approach as it shows turbulent dynamics even in the axisymmetric configuration described above. This enables the calculation of the statistics of turbulent states at computationally accessible resolutions, which is ideal for evaluating the effectiveness of the GQL approximation – enabling the analysis of systems far from criticality where the quasilinear approximation is expected to break down. The GQL representation has been examined in the context of two-dimensional driven turbulence on a spherical surface and
${\it\beta}$
-plane and shown to be more effective than the regular quasilinear (QL) approximation in reproducing both the dynamics and statistics of these flows away from equilibrium (Marston, Chini & Tobias Reference Marston, Chini and Tobias2016). The importance of establishing the regime of applicability of various approximations arises owing to their use in deriving statistical theories.
As noted by Tobias et al. (Reference Tobias, Dagon and Marston2011), the enormous range of spatial and temporal scales present in many astrophysical phenomena poses a problem for DNS in that each scale must be properly resolved; this is unlikely, even with today’s powerful parallel computers – and indeed those of the foreseeable future. The problem is that in most astrophysical applications, the large-scale dynamics is influenced by, and in turn influences, the small-scale interactions.
A complementary approach, that has been afforded much attention recently, is that of direct statistical simulation (DSS), in which the low-order statistics are obtained directly from a hierarchy of cumulant equations, and a closure is applied to this hierarchy in order to provide tractable calculations. DSS has a number of advantages over DNS, chiefly that the statistics vary slowly in time and are smoother in space and can be explored via the evolution of relatively few modes on a slow manifold. There have been a number of examples of direct statistical simulation in recent works, most notably involving various planetary jets; Tobias et al. (Reference Tobias, Dagon and Marston2011) derived a formulation for flow on rotating spherical surface and compared with the DNS solutions, showing considerable accuracy. Indeed, DSS is shown to be useful in analysing the underlying physics of flows in that it in this case demonstrates that jets may be formed without an inverse cascade. Moreover, a quasilinear statistical framework has allowed Squire & Bhattacharjee (Reference Squire and Bhattacharjee2015) to probe the MRI dynamo in a shearing box in ways inaccessible in conventional DNS. All of the above investigations used the most simple form of DSS (known as CE2) – this approximation is the statistical representation of QL. At this quasilinear level of approximation DSS is formally equivalent to the stochastic structural stability theory (SSST or S3T) of Farrell, Ioannou and collaborators (Farrell & Ioannou Reference Farrell and Ioannou2007; Farrell & Ioannou Reference Farrell and Ioannou2008; Constantinou, Farrell & Ioannou Reference Constantinou, Farrell and Ioannou2013), which has been extensively used to model a number of physical systems and can be justified for systems near equilibrium for which there is a separation of time-scales (Bouchet, Nardini & Tangarife Reference Bouchet, Nardini and Tangarife2013).
Of course, this is not to imply that quasilinear DSS (i.e. CE2) is universally applicable; there are a number of examples of inconsistencies with fully nonlinear DNS (NL DNS), evident for example in out-of-equilibrium jets studied with DSS for a number of parameter combinations. Far from equilibrium, CE2 ceases to be an adequate representation of the system, as shown by Tobias & Marston (Reference Tobias and Marston2013). In that paper it was demonstrated that when the system is in statistical equilibrium, and thus dominated by strong jets, CE2 does in fact do a good job of reproducing both the number and strength of jets in the system. However, the further away from equilibrium, the worse CE2 is seen to perform, resulting in both inaccurate jet strengths and inaccurate numbers of jets (see also Srinivasan & Young Reference Srinivasan and Young2012). Two options have been suggested to remedy the situation. The first is to include eddy–eddy scattering in the truncation. This leads to the CE3 (or CE2.5) approximation, which has been shown to improve performance of DSS away from equilibrium (Marston, Qi & Tobias Reference Marston, Qi and Tobias2014).
A second possibility is to generalise the quasilinear approximation (Marston et al. Reference Marston, Chini and Tobias2016) and derive a corresponding statistical theory (GCE2). We explore the accuracy of such DSS of the HMRI at GCE2 by performing a series of generalised quasilinear DNS (QL DNS) in the inductionless limit. First we formulate the model and give details of the methods utilised for DNS. We continue by defining the GQL approximation and compare its efficacy with that of QL theory by examining how well it reproduces a number of diagnostic quantities, such as the energy spectra and first and second cumulants of the flow. We conclude by briefly discussing the implications of our results to possible statistical simulations.
2 Formulation and numerical method
We consider axisymmetric cylindrical Taylor–Couette flow, with radii and angular velocities at the inner and outer cylinders given by
$r_{i},r_{o}$
and
${\it\Omega}_{i}$
and
${\it\Omega}_{o}$
respectively. We set
$r_{i}/r_{o}=0.5$
to coincide with the aspect ratio of the PROMISE experiment. Unlike an experiment however, with its inevitable end plates, we take the cylinders to be axially periodic with a basic periodicity of 40 times the gap width. This is sufficiently long to allow turbulent structures to develop in the axial direction; it is the dynamics and interactions of these structures that we then wish to explore in this work.
The flow and magnetic field are described by the Navier–Stokes and induction equations,
$$\begin{eqnarray}\displaystyle & \displaystyle Re\frac{\partial \boldsymbol{U}}{\partial t}-{\rm\nabla}^{2}\boldsymbol{U}=-\boldsymbol{{\rm\nabla}}p-Re\boldsymbol{U}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}\boldsymbol{U}+Ha^{2}Rm^{-1}(\boldsymbol{{\rm\nabla}}\times \boldsymbol{B})\times \boldsymbol{B}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\partial \boldsymbol{B}}{\partial t}=\boldsymbol{{\rm\nabla}}\times (\boldsymbol{U}\times \boldsymbol{B})+Rm^{-1}{\rm\nabla}^{2}\boldsymbol{B}, & \displaystyle\end{eqnarray}$$
where we have non-dimensionalised as follows: we scale lengths with the gap width
$(r_{o}-r_{i})\hat{r}$
, time on the rotational time-scale
$t={\it\Omega}_{i}^{-1}\hat{t}$
,
$U$
is proportional to length over time,
$U=(r_{o}-r_{i}){\it\Omega}_{i}\hat{U}$
and
$B=B_{z}\hat{B}$
, with
$B_{z}$
the strength of the imposed axial magnetic field. The non-dimensional quantities
$$\begin{eqnarray}Rm=(r_{o}-r_{i}){\it\Omega}_{i}/{\it\eta},\quad Re=(r_{o}-r_{i}){\it\Omega}_{i}/{\it\nu},\quad Ha=B_{z}(r_{o}-r_{i})/\sqrt{{\it\eta}{\it\nu}{\it\mu}_{0}{\it\rho}}\end{eqnarray}$$
are the magnetic and hydrodynamic Reynolds numbers and the Hartmann number respectively. Here
${\it\eta}$
is the magnetic diffusivity,
${\it\nu}$
is the kinematic viscosity,
${\it\mu}_{0}$
is the magnetic permeability of free space and
${\it\rho}$
is the fluid density.
As noted above, for the SMRI it is crucial that
$Rm$
be retained, and exceed
$O(10)$
. In contrast, the HMRI continues to exist even in the so-called inductionless limit
$Rm\rightarrow 0$
. Splitting the magnetic field into imposed and induced fields as
$\boldsymbol{B}=\boldsymbol{B}_{\mathbf{0}}+Rm\,\boldsymbol{b}$
and then letting
$Rm\rightarrow 0$
yields
$$\begin{eqnarray}\displaystyle & \displaystyle Re\frac{\partial \boldsymbol{U}}{\partial t}-{\rm\nabla}^{2}\boldsymbol{U}=-\boldsymbol{{\rm\nabla}}p-Re\boldsymbol{U}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}\boldsymbol{U}+Ha^{2}(\boldsymbol{{\rm\nabla}}\times \boldsymbol{b})\times \boldsymbol{B}_{\boldsymbol{ 0}}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle 0={\rm\nabla}^{2}\boldsymbol{b}+\boldsymbol{{\rm\nabla}}\times (\boldsymbol{U}\times \boldsymbol{B}_{\mathbf{0}}). & \displaystyle\end{eqnarray}$$
Note also that in this limit the only nonlinearity arises from the inertial term
$\boldsymbol{U}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}\boldsymbol{U}$
.
The imposed field is given by
$\boldsymbol{B}_{\mathbf{0}}=\boldsymbol{e}_{z}+{\it\beta}r^{-1}\boldsymbol{e}_{{\it\phi}}$
; the HMRI operates when
${\it\beta}\geqslant O(1)$
(Hollerbach & Rüdiger Reference Hollerbach and Rüdiger2005). The rotation ratio is fixed at
${\it\Omega}_{o}/{\it\Omega}_{i}=0.27$
, slightly beyond the Rayleigh value of 0.25. That is, for
${\it\Omega}_{o}/{\it\Omega}_{i}<0.25$
one would obtain purely hydrodynamic instabilities, which we are not interested in here, but for
${\it\Omega}_{o}/{\it\Omega}_{i}>0.25$
one obtains only the HMRI, for appropriate values of the three governing parameters
$Re$
,
$Ha$
and
${\it\beta}$
.
The boundary conditions on
$\boldsymbol{U}$
are no-slip as they must be for any Couette flow. For
$\boldsymbol{b}$
the HMRI exists in much the same form for either insulating (Hollerbach & Rüdiger Reference Hollerbach and Rüdiger2005) or perfectly conducting (Rüdiger et al.
Reference Rüdiger, Hollerbach, Stefani, Gundrum, Gerbeth and Rosner2006) boundaries. Since figure 2 of Rüdiger et al. (Reference Rüdiger, Hollerbach, Stefani, Gundrum, Gerbeth and Rosner2006) mapped out a particularly thorough set of linear onset curves for perfectly conducting boundary conditions, we will here consider only perfectly conducting boundaries.
Because the HMRI is an intrinsically axisymmetric phenomenon (just like the SMRI), we restrict attention to axisymmetric solutions, and further decompose as
$$\begin{eqnarray}\boldsymbol{U}=v\boldsymbol{e}_{{\it\phi}}+\boldsymbol{{\rm\nabla}}\times ({\it\psi}\boldsymbol{e}_{{\it\phi}}),\quad \boldsymbol{b}=a\boldsymbol{e}_{{\it\phi}}+\boldsymbol{{\rm\nabla}}\times (b\boldsymbol{e}_{{\it\phi}}),\end{eqnarray}$$
thereby automatically satisfying
$\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\boldsymbol{U}=\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\boldsymbol{b}=0$
.
The full details of the numerical solution for
$v$
,
${\it\psi}$
,
$a$
and
$b$
are presented in Hollerbach (Reference Hollerbach2008). Here we merely note that the radial structure is expanded in 50 Chebyshev polynomials, and the axial structure in 600 Fourier modes of the form
$\exp (\text{i}{\it\kappa}z)$
, where
${\it\kappa}=2{\rm\pi}k/z_{0}$
, with
$z_{0}=40$
the imposed basic periodicity. The essence of the GQL again is to study interactions between different Fourier modes, so such a decomposition is convenient both numerically and theoretically.
2.1 Generalised quasilinear approximation
In this section we outline how to apply the generalised quasilinear approximation introduced by Marston et al. (Reference Marston, Chini and Tobias2016) to this axisymmetric HMRI system. For simplicity, let us consider the velocity
$\boldsymbol{U}$
, although similar considerations are applied to the perturbed magnetic field
$\boldsymbol{b}$
. In general the fields may be split via a standard Reynolds decomposition into mean and fluctuation (eddy) components, for example for the azimuthal velocity
$$\begin{eqnarray}v=\overline{v}+v^{\prime },\end{eqnarray}$$
where overbar denotes a mean and prime denotes fluctuations about that mean. In previous papers the averaging procedure has involved taking a zonal mean, but in this axisymmetric geometry a natural definition of the mean is that part that is independent of the
$z$
-coordinate (i.e. the
${\it\kappa}=0$
mode). This standard definition for the mean and fluctuations (eddies) then naturally leads to the quasilinear approximation. In this approximation certain mode interactions are discarded from the dynamics. In particular whilst mean/eddy
$\rightarrow$
eddy and eddy/eddy
$\rightarrow$
mean interactions are kept, eddy/eddy
$\rightarrow$
eddy interactions are discarded. In this quasilinear approximation therefore only non-local (in wavenumber) elements of cascades or inverse cascades are possible (Tobias et al.
Reference Tobias, Dagon and Marston2011).
A similar decomposition is used for GQL. However here the velocity is decomposed into large- and small-scale modes, where large and small scale are defined by the spectral wavenumber. That is we set
$$\begin{eqnarray}v(r,z)=\mathop{\sum }_{k=0}^{{\it\Lambda}}v_{lk}(r)\text{e}^{\text{i}{\it\kappa}z}+\mathop{\sum }_{k={\it\Lambda}+1}^{N_{k}}v_{hk}(r)\text{e}^{\text{i}{\it\kappa}z},\end{eqnarray}$$
where the
$v_{l}$
and
$v_{h}$
are the large-scale (denoted
$l$
for ‘low’ modes) and smaller-scale (denoted
$h$
for ‘high’ modes) modes respectively, and
${\it\Lambda}$
is the mode cutoff point between the two. It is clear then that when
${\it\Lambda}=0$
this expansion is the same as that for QL given in (2.7). Hence the GQL formalism reverts to QL when
${\it\Lambda}=0$
as noted by Marston et al. (Reference Marston, Chini and Tobias2016).
When
${\it\Lambda}\neq 0$
and the generalised quasilinear approximation differs from QL, the triad interactions that are retained and discarded are selected so as to satisfy the relevant conservation laws and enable closure (Marston et al.
Reference Marston, Chini and Tobias2016). Briefly this involves retaining low/low
$\rightarrow$
low, high/high
$\rightarrow$
low and low/high
$\rightarrow$
high interactions and discarding all others. Again retaining this set of interactions is consistent with QL when
${\it\Lambda}=0$
. Furthermore as
${\it\Lambda}\rightarrow \infty$
the system consists solely of low modes and is formally the same as direct numerical simulation. The aim of this paper is to determine how well GQL performs in reproducing the NL DNS results as
${\it\Lambda}$
is varied; clearly if it is able to do so for relatively small
${\it\Lambda}$
then this is extremely encouraging for producing direct statistical simulation based on the GQL approximation.
3 Results
3.1 Means of comparison: outputs of interest
As DSS is a statistical theory, we shall evaluate the GQL approximation by accumulating statistics using QL DNS and testing how well GQL reproduces the statistics of the HMRI (Marston et al. (Reference Marston, Chini and Tobias2016) evaluated both the statistics and the dynamics for the problem of the driving of zonal jets). To that end we shall calculate average vertical spectra of the flow and magnetic field as well as the first and second cumulants for the flow. The first cumulant is the mean over
$z$
and time as a function of
$r$
defined as
$$\begin{eqnarray}c_{{\it\psi}}(r)=\langle {\it\psi}(r,z,t)\rangle ,\quad c_{v}(r)=\langle v(r,z,t)\rangle ,\end{eqnarray}$$
where angled brackets indicate averaging over
$z$
and
$t$
. The second cumulants yield information about the (non-local) fluctuation/fluctuation interactions and are defined as
$$\begin{eqnarray}\displaystyle & \displaystyle c_{{\it\psi}{\it\psi}}(r_{1},r_{2},{\it\xi})=\langle {\it\psi}^{\prime }(r_{1},z_{1},t){\it\psi}^{\prime }(r_{2},z_{2},t)\rangle , & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle c_{vv}(r_{1},r_{2},{\it\xi})=\langle v^{\prime }(r_{1},z_{1},t)v^{\prime }(r_{2},z_{2},t)\rangle , & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle c_{v{\it\psi}}(r_{1},r_{2},{\it\xi})=\langle v^{\prime }(r_{1},z_{1},t){\it\psi}^{\prime }(r_{2},z_{2},t)\rangle , & \displaystyle\end{eqnarray}$$
where
${\it\xi}=z_{2}-z_{1}$
. The second cumulant can then be expressed as, for example
$$\begin{eqnarray}c_{{\it\psi}{\it\psi}}(r_{1},r_{2},{\it\xi})=\iint {\it\psi}^{\prime }(r_{1},z_{1}){\it\psi}^{\prime }(r_{2},z_{1}+{\it\xi})\,\text{d}z_{1}\,\text{d}t.\end{eqnarray}$$
In this paper, we choose
$r_{1}$
to be the point halfway between the two cylinders.
3.2 Parameter choice and initialisation
The generalised quasilinear approximation is straightforward to implement in the pseudo-spectral time-stepping code – the various triad interactions discussed by Marston et al. (Reference Marston, Chini and Tobias2016) can be easily isolated in spectral space and discarded or retained as required. We performed a number of simulations under various different parameter regimes for a number of truncations of the large-scale modes
${\it\Lambda}$
.
We initialise each simulation using the following procedure. For each parameter set we run a fully turbulent NL DNS solution until it reaches a statistically stationary state. It is this state that we use as initial conditions for runs at the same parameter values but with different cutoffs for the GQL approximation. We note that the same results are obtained if the GQL calculations are started from rest – although these take longer to equilibrate. For each parameter set we evolve
${\it\Lambda}=0$
(corresponding directly to the QL expansion) and we also examine
${\it\Lambda}=1,2,3,6$
and
${\it\Lambda}=20$
. We evolve each of the equations using a time step
$\text{d}t=5\times 10^{-3}$
, and run each simulation for a minimum of
$8\times 10^{4}$
time steps to make sure that the solutions reached are statistically steady. Any averages (over time) are then calculated using 300 equally spaced (in time) solutions, commencing from the 200 000th step.
We also give here some note about the rationale for our choice of model parameters – care has to be taken to make sure that
$Re$
is sufficiently above
$Re_{c}$
, the critical Reynolds number. Intuitively, one might expect that the supercriticality of the system as measured by
$(Re-Re_{c})$
will have bearing on the effectiveness of the GQL approximations, and so we choose to examine a number of different parameter sets. Based on figure 2 of (Rüdiger et al.
Reference Rüdiger, Hollerbach, Stefani, Gundrum, Gerbeth and Rosner2006), we choose to focus on three specific choices;
${\it\beta}=5$
,
$Ha=16$
,
${\it\beta}=2.5$
,
$Ha=16$
and
${\it\beta}=10$
,
$Ha=8$
, each for
$Re=5000,7000$
and 10 000. These choices ensure that the flow is sufficiently supercritical to generate a turbulent state. Of the three parameters sets,
${\it\beta}=5$
,
$Ha=16$
has
$Re_{c}\approx 1000$
, and is the most supercritical, and
${\it\beta}=2.5$
,
$Ha=16$
has a slightly higher
$Re_{c}\approx 1500$
, and thus will give insight into how the GQL approximation holds for less turbulent flows. The point
${\it\beta}=10$
,
$Ha=8$
is more interesting in that, even though
$Re_{c}<1000$
, on the aforementioned figure the contours of
$Re_{c}$
are closely bunched together. That is, the flow would be viewed as being highly supercritical when the Reynolds parameter alone is examined, but close to marginal stability in terms of the magnetic field parameters. This is of particular interest in that it suggests that this parameter set may be only weakly nonlinear, for which we would expect even the standard QL approximation to give an accurate reproduction of the NL DNS. We claim that these three points and
$Re$
values sufficiently encapsulate the properties of the GQL approximation; simulations for other parameters yield qualitatively similar results.
3.3 NL DNS
We begin by describing the behaviour of the system for NL DNS of the HMRI. In figure 1, we show plots of various quantities – contours of the azimuthal velocity, Hovmöller diagrams, mean azimuthal flows (first cumulants), energy spectra and second cumulants – for each parameter combination at the highest value of the Reynolds number
$Re=10\,000$
. As seen in figure 1(a–c), different parameter choices lead to varying degrees of turbulence, with
${\it\beta}=5$
,
$Ha=16$
(a,d,g,j,m,p) being the strongest, and
${\it\beta}=10$
,
$Ha=8$
(c,f,i,l,o,r) – that which we suspect is only weakly nonlinear – being the weakest.
(a–c) Contour snapshots of the azimuthal velocity in a radial cross-section; (d–f) Hovmöller plots of
${\it\psi}$
as a function of
$z$
and time; (g–i) mean azimuthal flow (first cumulant) as a function of
$r$
; (j–l) vertical energy spectra for the magnetic field (red) and velocity (blue) with poloidal components (solid) and toroidal components (dashed); (m–o) second cumulants
$\langle v^{\prime }v^{\prime }\rangle$
; (p–r) second cumulants
$\langle v^{\prime }{\it\psi}^{\prime }\rangle$
. Shown on (a,d,g,j,m,p) is
${\it\beta}=5$
,
$Ha=16$
, (b,e,h,k,n,q) is
${\it\beta}=2.5$
,
$Ha=16$
and on (c,f,i,l,o,r) is
${\it\beta}=10$
,
$Ha=8$
, each for
$Re=10^{4}$
.
The HMRI is known to be a unidirectionally travelling wave at onset (Hollerbach & Rüdiger Reference Hollerbach and Rüdiger2005; Rüdiger et al.
Reference Rüdiger, Hollerbach, Stefani, Gundrum, Gerbeth and Rosner2006), due to the
$\pm z$
symmetry breaking caused by an imposed field having both axial and azimuthal components (Knobloch Reference Knobloch1996). This drifting nature is conveniently visualised by presenting a number of Hovmöller plots of
${\it\psi}(t,z)$
for a given radius, which allows the identification of the location of Taylor vortices. The streamfunction
${\it\psi}$
undergoes a reversal in sign when transitioning from one Taylor vortex to another. Therefore, simply plotting the sign of
${\it\psi}$
is equivalent to tracking the position of the vortices over time. It can be seen clearly that in DNS the vortices propagate downward in the axial direction as a travelling wave. Also apparent on these plots is the number of dislocations present in the flow; these are an effective measure of the strength of turbulence. It can be seen that in general, as expected, the number of dislocations decreases as
$Re_{c}$
increases for fixed
$Re$
, i.e. as the flow becomes less supercritical. However, at the point identified as only weakly nonlinear few of these nonlinear dislocations occur. The mean azimuthal flows shown in the third row demonstrate the change from the basic state. As expected the role of turbulence is to decrease the gradient of the mean azimuthal flow (removing the source of instability). For the more turbulent flows this change is more significant whilst for the less supercritical flows the departure from the mean state is much less pronounced, as expected. The resulting outward transport of angular momentum, essential to MRI, is also evident.
Turning our attention to the energy spectra, for which we plot the energy of both the toroidal and poloidal components of the magnetic field (in red) and velocity (in blue), we may note that regardless of the parameter regime these are qualitatively very similar. The magnetic energy is smaller than that for the kinetic energy and has a steeper spectrum. In the weakly nonlinear case the energy peaks at a mode number
$k\approx 13$
, which is similar to the linear mode of maximum growth rate. The energy spectra for
$Re=7000$
and
$Re=5000$
have similar characteristics to these plots. The second cumulants
$\langle v^{\prime }v^{\prime }\rangle$
and
$\langle v^{\prime }{\it\psi}^{\prime }\rangle$
are very much localised in the axial direction, although there is strong radial correlation with the midpoint in radius. The axial localisation is a reflection of the flat energy spectra for the low modes, which leads to this
${\it\delta}$
–function–like axial structure in the second cumulants. For the weaker turbulence (third column), the cumulants are different in that a periodic axial correlation is now apparent, such that flow at a point is correlated not only in the radial direction, but also periodically in each direction along the axis. A wavenumber of
$k\approx 13$
in the second cumulants is observed, which corresponds to the peak observed in the energy spectrum.
These results will be used in the next section as benchmarks to evaluate the effectiveness of the GQL approximation at various truncations.
3.4 GQL DNS
3.4.1 Weak turbulence:
${\it\beta}=10,Ha=8$
We start by examining the effectiveness of GQL DNS in comparison with QL DNS for the weak turbulence case
${\it\beta}=10$
,
$Ha=8$
,
$Re=10\,000$
(for which the DNS results are shown in the third column of figure 1). It is to be expected that the QL approximation will most likely have the greatest success in this parameter regime, as nonlinearities have smaller influence on the instability; recall QL is an exact representation of the dynamics in the formally weakly nonlinear regime. We can see on comparing the Hovmöller plots (figure 2) that this expectation is correct; there are few dislocations in the NL DNS, meaning that even
${\it\Lambda}=0$
is an accurate approximation.
Snapshot contour plots of
$u_{{\it\phi}}$
(a–c) and Hovmöller plots (d–f) for
${\it\beta}=10$
,
$Ha=8$
,
$Re=10^{4}$
. Shown on (a,d) is the NL DNS, (b,e) is
${\it\Lambda}=0$
, and on (c,f) is
${\it\Lambda}=2$
.
We can see from the first cumulants (figure 3
a) that the various
${\it\Lambda}$
approximations are nearly indistinguishable from the NL DNS profile. Even the quasilinear
${\it\Lambda}=0$
gives an excellent approximation to the full solution. Indeed, for the first cumulants there is little advantage in using a GQL approximation over a QL approximation. Again for such a laminar solution this is perhaps to be expected.
First cumulants
$\langle v\rangle$
for
${\it\beta}=10$
$Ha=8$
(a),
${\it\beta}=5$
$Ha=16$
(b) and
${\it\beta}=2.5$
$Ha=16$
(c). Note that for
${\it\beta}=10$
$Ha=16$
, as the flow is only weakly nonlinear the NL DNS and GQL DNS profiles are near indistinguishable.
Moreover, even
${\it\Lambda}=0$
reproduces axial structure the second cumulants fairly well at
$Re=10\,000$
(figure 4). The axial wavelength and magnitudes are marginally inaccurate and the radial correlation is focussed at the midpoint instead of across the full gap width, however they are surprisingly good for such a crude approximation to the nonlinearity. Interestingly, the correlation structure is still finer than that of the NL DNS, due to the peak at
$k\approx 20$
that is present in all QL solutions. This structure becomes finer as
${\it\Lambda}$
is increased, although the magnitudes of the spurious correlations decrease. However, we comment that
${\it\Lambda}=0$
does not give universally adequate approximations even in this weakly turbulent case. If we examine
$Re=5000$
(not shown), the cumulant structure is much less accurate; it is not until
${\it\Lambda}=2$
, that the axial cumulant structure resembles the NL DNS case. So, whilst for much of this weakly nonlinear parameter regime the standard QL approximation produced adequate first and second cumulants, this cannot be guaranteed for all
$Re$
– the
${\it\Lambda}=2$
GQL is required to further guarantee this level of accuracy.
Second cumulants
$\langle v^{\prime }v^{\prime }\rangle$
as a function of
$r_{2}$
and
${\it\xi}$
at
${\it\beta}=10$
$Ha=8$
for NL DNS (a), QL DNS (
${\it\Lambda}=0$
) (b), and GQL DNS with
${\it\Lambda}=2$
(c),
${\it\Lambda}=3$
(d). Recall that here we choose
$r_{1}$
to be the point halfway between the two cylinders.
3.4.2 Strong turbulence:
${\it\beta}=5$
$Ha=16$
We next examine the GQL approximation at the set of parameters
${\it\beta}=5$
,
$Ha=16$
, for which we plot the same diagnostic quantities as for NL DNS. From figure 5, it is clear that the azimuthal velocity profile is quite different in the
${\it\Lambda}=0$
case to that for NL DNS, with the QL approximation failing to reproduce dislocations. Importantly however, we note that the Hovmöller plots show that the travelling wave property is preserved, even by the standard
${\it\Lambda}=0$
QL approximation. The wavelength, although difficult to measure due to dislocations, can be seen to be within an acceptable range of the NL DNS. Note that if the
${\it\Lambda}=2$
GQL approximation is taken, even the nonlinear dislocations are reproduced; a significant improvement over
${\it\Lambda}=0$
.
Snapshot contour plots of
$u_{{\it\phi}}$
(a–c) and Hovmöller plots (d–f) for
${\it\beta}=5$
,
$Ha=16$
,
$Re=10^{4}$
. Shown on (a,d) is the NL DNS, (b,e) is
${\it\Lambda}=0$
, and on (c,f) is
${\it\Lambda}=2$
.
Additionally, we compare the energy spectrum to that of the NL DNS (figure 6). It should be noted that agreement in the lower modes of the spectrum is the most important for the evaluation of the approximation. We observe that the energy in the GQL approximation is located primarily in the lower modes, which peak higher than in the original NL DNS case. Because energy is not being transferred by nonlinear processes the higher modes contain much less, if any, energy and the turbulent cascade is no longer present.
Vertical energy spectra for the various mode truncations at
${\it\beta}=5$
,
$Ha=16$
,
$Re=10\,000$
with magnetic field (red) and velocity (blue) and with poloidal components (solid) and toroidal components (dashed).
The most striking features in these plots are the strong peaks in the spectra for GQL at
${\it\Lambda}=0$
and
${\it\Lambda}=1$
, which correspond to a linear instability and its harmonics. The first peak occurs at a mode number
$k=40$
, and the other peaks appear to be the resonance of this at
$k=80$
and
$k=160$
. A key question one may ask is whether this is due to the initial condition – recall that we started with a fully turbulent state from which the energy could have decayed away in the upper spectra. As noted earlier, even if the simulations are started from a small perturbation about the basic state these peaks still exist in the spectra, therefore the energy must be transferred outwards to them. For
${\it\Lambda}=1$
we still detect these linear peaks, although due to the increased energy scattering off the low modes, some energy has now transferred outwards to the modes between them.
If we turn our attention to
${\it\Lambda}=2$
, it can be seen that the energy in higher modes becomes more like that in the NL DNS, even if it does decay at a much greater rate for
$k>200$
. The linear instability peak is still present, however it is but a small protrusion above the energy stored in adjacent modes, and for the modes immediately higher than
$k=20$
the spectrum begins to resemble the power law from the NL DNS. The improvement in reproduction of the DNS spectrum over
${\it\Lambda}=0$
is clear, yet it can also be seen that increasing
${\it\Lambda}>2$
does not yield similarly great increases in accuracy (
${\it\Lambda}=6$
produces an almost equivalent spectrum). Interestingly there is a sharp decrease in energy around the mode
$k=20$
, which persists even for
${\it\Lambda}=20$
, where the rest of the spectra resembles the NL DNS. For this large
${\it\Lambda}$
, we see that the energy is higher than the NL DNS in the modes slightly lower, and lower in the slightly higher modes, whilst virtually identical to the NL DNS otherwise.
If we turn our attention to the first cumulants
$\langle v\rangle$
in figure 3, we see that for
${\it\Lambda}=0$
(QL) they deviate approximately 10 % from the NL DNS solution. As
${\it\Lambda}$
is increased the cumulants then approach the profile of NL DNS, with a large improvement in agreement shown between
${\it\Lambda}=1$
and
${\it\Lambda}=2$
, in agreement with what is seen in the energy spectra. The deviation gets smaller as
$Re$
is decreased (not shown) but decreasing
$Re$
does not lead to a better approximation; the lower Reynolds number simply means that the cylinder is rotating with lower angular velocity and thus the overall perturbation in the velocity profile is smaller, with the percentage error for QL staying approximately the same.
We also compare plots of the second cumulant,
$\langle v^{\prime }v^{\prime }\rangle$
. Recall we choose to display the cross-correlation with the midpoint of the gap between the cylinders. The choice of
$Re$
has a slight influence on the second cumulants but does not significantly alter their axial structure. Thus, we choose to examine the case when
$Re=7000$
to illustrate the behaviour of the second cumulants in the GQL approximation.
It is immediately clear from figure 7 that the second cumulants for
${\it\Lambda}=0$
and
${\it\Lambda}=1$
are inadequate reproductions of the NL DNS case – for these approximations there is an uncharacteristically large observed axial correlation. Indeed, if we take
${\it\Lambda}=0$
we see that the correlation has a wavenumber of 5, whereas the NL DNS is strongly correlated with one axial location, and is very weakly (or un-) correlated with the rest of the domain. Interestingly,
${\it\Lambda}=1$
has correlation of wavenumber 4, possibly associated with the presence of waves travelling axially (see also Marston et al.
Reference Marston, Chini and Tobias2016). Recall that the localisation of the second cumulant appeared to be associated with the absence of a strong spectral peak. As such, by
${\it\Lambda}=2$
the axial localisation of the second cumulant has been re-established.
Second cumulants
$\langle v^{\prime }v^{\prime }\rangle$
, as in figure 4, at a number of mode truncations for the parameter regime
${\it\beta}=5$
$Ha=16$
,
$Re=7000$
. Shown are; NL DNS (a), QL DNS (
${\it\Lambda}=0$
) (b), and GQL DNS with
${\it\Lambda}=1$
(c),
${\it\Lambda}=2$
(d).
We note similar structure in the cross-cumulants
$\langle v^{\prime }{\it\psi}^{\prime }\rangle$
, which are plotted in figure 8. A similar overcorrelation is seen for
${\it\Lambda}=0$
and
${\it\Lambda}=1$
, although
$v^{\prime }$
is out of phase with
${\it\psi}^{\prime }$
. At
${\it\Lambda}=2$
we see that the structure of the NL DNS is regained, with all but the axial midpoint being very weakly correlated. Note how now the radial correlation is much more consistent across the NL DNS and GQL approximations than previous.
As for figure 7 but for cross-cumulants
$\langle v^{\prime }{\it\psi}^{\prime }\rangle$
.
3.4.3
${\it\beta}=2.5,Ha=16$
: moderate turbulence
As mentioned previously, this parameter set is much less turbulent than
${\it\beta}=5$
,
$Ha=16$
, with fewer dislocations in the NL DNS. Thus, it can inform us of the effect that
$Re_{c}$
has on the viability of the GQL, in a similar fashion to the reduced
$Re$
cases previously. Much like before, it can be found that the fundamental travelling wave property is reproduced even for
${\it\Lambda}=0$
, and it is not until
${\it\Lambda}=2$
that nonlinear dislocations are again found in Hovmöller plots.
If we examine the first cumulants in figure 3 we see that they have profiles that are very close to that of the NL DNS, whatever the degree of approximation – although improvement is seen as
${\it\Lambda}$
is increased. We briefly comment that the energy spectra and second cumulants show similar behaviour to those for the strong turbulence case described above, with GQL beginning to yield a good approximation by
${\it\Lambda}=5$
. Hence we conclude that even if QL appears to give an adequate approximation on examination of the mean flows, GQL may be required for accurate representation of the fluctuations.
4 Conclusion
In this paper we have extensively examined the effectiveness of the generalised quasilinear approximation (GQL) for the case of an axisymmetric HMRI for a number of different parameter regimes, with flows varying from weakly to strongly turbulent. Establishing the approximation’s effectiveness is important not only for determining which interactions play a key role in the dynamics of the system, but also for determining the minimum ingredient of a statistical theory (such as direct statistical simulation) that is necessary to reproduce the low-order statistics of the system.
For weakly turbulent systems we find no fundamental problems with the possibility of DSS and the HMRI – even for the QL case of
${\it\Lambda}=0$
the travelling wave behaviour is reproduced, even if some of the other properties are of varying degrees of accuracy. However when the system is more turbulent we have found that for the GQL approximation at
${\it\Lambda}=2$
, much more accurate results can be obtained, and that increasing
${\it\Lambda}$
further yields diminishing returns, with a much lesser improvement in accuracy per extra mode included. We have seen that
${\it\Lambda}=2$
gives a much improved spectra compared with
${\it\Lambda}=0$
. The various cumulants further reinforce this;
${\it\Lambda}=0$
for the most part yields overcorrelated second cumulants, and it is not until
${\it\Lambda}=2$
that the correlation structures properly resemble the NL DNS.
The choice of parameter regime has a knock–on effect in the accuracy of the GQL approximation, with the degree of supercriticality an important factor. However, for the majority of the parameter space the influence that
${\it\beta}$
,
$Ha$
and
$Re$
have on the approximations has little effect on the relative accuracy of the different
${\it\Lambda}$
approximations; in each regime it was found that increasing to
${\it\Lambda}=2$
offers a marked improvement over the
${\it\Lambda}=0$
case whilst increasing
${\it\Lambda}>2$
yields much smaller increases in accuracy. Special points, such as the weakly nonlinear
${\it\beta}=10$
,
$Ha=8$
may be approximated to a good degree of accuracy using even the
${\it\Lambda}=0$
approximation. However, these are but a small section of the total parameter space.
This leads us to two main conclusions. Firstly, as the HMRI is the subject of much current research, it may be fruitful to perform direct statistical simulations of the instability. Having shown that using a GQL approximation may yield accurate results, the equivalent cumulant expansion to the
${\it\Lambda}=2$
case may be used to give further insight into the instability.
Secondly, and most importantly, this work shows that for some problems it may be advantageous to use cumulant expansions other than the quasilinear CE2 expansion to conduct direct statistical simulation. As mentioned in the introduction, previous research has shown that DSS via the CE2 expansion (equivalent to
${\it\Lambda}=0$
), whilst accurate for problems such as out-of-equilibrium jets, can also prove inaccurate. We have shown that for the HMRI, where taking a standard QL approximation would have been inadequate, even the addition of two more large-scale modes in the
${\it\Lambda}=2$
approximation yields markedly more accurate solutions. It is logical to conclude that a similar
${\it\Lambda}\geqslant 2$
cumulant expansion would allow for accurate DSS in other flows, without compromising the computational advantage of DSS over DNS. These results reinforce those of Marston et al. (Reference Marston, Chini and Tobias2016) who found that the GQL approximation may yield significant enhancement over QL for the hydrodynamic problem of jet formation on a spherical surface or
${\it\beta}$
-plane. In that system the large-scale shear flows emerge as a result of systematic correlations in the turbulence, rather than being imposed through stresses at the boundary. Marston et al. (Reference Marston, Chini and Tobias2016) found that the QL approximation tended to overestimate the long-range correlations and conjectured that this was due to the approximation neglecting the self-interaction of nonlinear Rossby waves and oversimplifying their interaction with the mean flows. This in turn may lead to an incorrect calculation of the Reynolds stress and therefore the mean flows (jets). In this paper, the dynamics is different, with energy being injected at large scales and the instability of a systematic flow leading to turbulence and angular momentum transport. Again the QL approximation neglects the self-interaction of the turbulence and simplifies the response of the turbulence to the means flows. In both cases even including a few modes in the large scales appears to significantly improve the effectiveness of the approximation; we believe that this enables the scattering of energy off the large scales (sometimes termed non-local transfer) which appears to be a crucial factor in reproducing the statistics of the flow.
Finally, we note that while we have here only considered the axisymmetric HMRI, for only slightly larger Hartmann numbers there is also a non-axisymmetric analogue, the so-called azimuthal MRI (AMRI). This mode undergoes a similar switch in scaling from
$Rm$
to
$Re$
as in the SMRI
$\rightarrow$
HMRI transition (Hollerbach, Teeluck & Rüdiger Reference Hollerbach, Teeluck and Rüdiger2010), and has by now also been observed in the PROMISE experiment (Seilmayer et al.
Reference Seilmayer, Galindo, Gerbeth, Gundrum, Stefani, Gellert, Rüdiger, Schultz and Hollerbach2014). Three-dimensional numerical modelling of these modes is only just beginning (Guseva et al.
Reference Guseva, Willis, Hollerbach and Avila2015), but a GQL analysis as presented here would obviously also be of considerable interest, with the separation into different scales then possible in both axial and azimuthal directions. Indeed encouragement for this enterprise is obtained from results for the hydrodynamic rotating Couette problem which indicate that GQL remains an attractive approximation when averaging over two spatial directions (in that case streamwise and spanwise directions).
Acknowledgements
A.C. was supported by an STFC studentship, R.H. and S.T. were supported by STFC grant no. ST/K000853/1, and B.M. was supported by NSF grant no. DMR-1306806.
