2.1. The long-wavelength stability problem

The basic state of our problem comprises statistically steady and homogeneous incompressible MHD turbulence, driven by a prescribed forcing $\boldsymbol {F}(\boldsymbol {x},t)$. It is supposed that the velocity field $\boldsymbol {U}$ and the magnetic field $\boldsymbol {B}$ are periodic in a cuboidal domain with typical size $L$. After a standard non-dimensionalisation, the momentum and induction equations take the form

(2.1)$$\begin{gather} \frac{\partial \boldsymbol{U}}{\partial t} +\boldsymbol{U}\boldsymbol{\cdot }\boldsymbol{\nabla}\boldsymbol{U} ={-}\boldsymbol{\nabla} P + \boldsymbol{B} \boldsymbol{\cdot }\boldsymbol{\nabla}\boldsymbol{B} + Re^{{-}1} \nabla^2 \boldsymbol{U} +\boldsymbol{F}, \end{gather}$$

(2.2)$$\begin{gather}\frac{\partial \boldsymbol{B}}{\partial t} +\boldsymbol{U}\boldsymbol{\cdot }\boldsymbol{\nabla}\boldsymbol{B} = \boldsymbol{B}\boldsymbol{\cdot } \boldsymbol{\nabla}\boldsymbol{U} + Rm^{{-}1} \nabla^2 \boldsymbol{B} . \end{gather}$$

The Reynolds number $Re$ and the magnetic Reynolds number $Rm$ are defined as usual by $Re=\hat {U}L/\nu$ and $Rm=\hat {U}L/\eta$, where $\hat {U}$ is a typical velocity, $\nu$ is the kinematic viscosity and $\eta$ is the magnetic diffusivity. Here we prescribe the forcing and not the velocity, and so $\hat {U}^2/L=\hat {F}$, where $\hat {F}$ is a typical magnitude of the forcing. We are interested in finding the growth rate of very long-wavelength disturbances to this basic state. We assume that these disturbances are small, so that the perturbation equations can be linearised. We start by writing down the linearised equations for the perturbation quantities, denoted by $\boldsymbol {u}$, $p$ and $\boldsymbol {b}$:

(2.3)$$\begin{gather} \frac{\partial \boldsymbol{u}}{\partial t} +\boldsymbol{U}\boldsymbol{\cdot }\boldsymbol{\nabla}\boldsymbol{u} +\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{\nabla}\boldsymbol{U} ={-}\boldsymbol{\nabla} p + \boldsymbol{B}\boldsymbol{\cdot }\boldsymbol{\nabla}\boldsymbol{b} + \boldsymbol{b}\boldsymbol{\cdot }\boldsymbol{\nabla}\boldsymbol{B}+ Re^{{-}1} \nabla^2 \boldsymbol{u}, \end{gather}$$

(2.4)$$\begin{gather}\frac{\partial \boldsymbol{b}}{\partial t} +\boldsymbol{U}\boldsymbol{\cdot }\boldsymbol{\nabla}\boldsymbol{b} +\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{\nabla}\boldsymbol{B} = \boldsymbol{B}\boldsymbol{\cdot }\boldsymbol{\nabla}\boldsymbol{u} +\boldsymbol{b}\boldsymbol{\cdot }\boldsymbol{\nabla}\boldsymbol{U} + Rm^{{-}1} \nabla^2 \boldsymbol{b} . \end{gather}$$

We now assume that these quantities evolve on two disparate length and time scales. By averaging over intermediate scales, we can therefore write $\boldsymbol {u}=\bar {\boldsymbol {u}}+\boldsymbol {u}'$, etc., where the barred (mean) variables vary only on length and time scales, denoted by $\boldsymbol {X}$ and $T$, respectively, that are very long compared with the scales $\boldsymbol {x}$ and $t$ of the basic state. Equations (2.3) and (2.4) may then be decomposed into mean and fluctuating parts. The mean equations take the form

(2.5)$$\begin{gather} \frac{\partial \bar{\boldsymbol{u}}}{\partial T} +\frac{\partial}{\partial X_j}(\overline{U_j \boldsymbol{u}'} + \overline{u'_j \boldsymbol{U}}) ={-}\boldsymbol{\nabla}_{\!\!{\tiny X}} \bar{p} +\frac{\partial}{\partial X_j}(\overline{B_j \boldsymbol{b}'} + \overline{b'_j \boldsymbol{B}})+ Re^{{-}1} {\nabla}_{\!\!{\tiny X}}^2 \bar{\boldsymbol{u}}, \end{gather}$$

(2.6)$$\begin{gather}\frac{\partial \bar{\boldsymbol{b}}}{\partial T} = \boldsymbol{\nabla}_{\!\!{\tiny X}}\times( \overline{\boldsymbol{U}\times\boldsymbol{b}'} + \overline{\boldsymbol{u}'\times\boldsymbol{B}}) + Rm^{{-}1} {\nabla}_{\!\!{\tiny X}}^2 \bar{\boldsymbol{b}} , \end{gather}$$

where $(\boldsymbol {\nabla }_{\!\!{\tiny X}})_i=\partial /\partial X_i$. To close these equations we need to express the averaged quantities in terms of $\bar {\boldsymbol {u}}$ and $\bar {\boldsymbol {b}}$. If we wish only to find terms including first derivatives of $\bar {\boldsymbol {u}}$ and $\bar {\boldsymbol {b}}$ with respect to $\boldsymbol {X}$ (since these terms will dominate at sufficiently long scales), we can self-consistently neglect derivatives of $\bar {\boldsymbol {u}}$, $\bar {\boldsymbol {b}}$ in the equations for the fluctuations. With this approximation, these then take the form

(2.7) \begin{align} \frac{\partial \boldsymbol{u}'}{\partial t} +(\boldsymbol{U}\boldsymbol{\cdot }\boldsymbol{\nabla}\boldsymbol{u}' +\boldsymbol{u}'\boldsymbol{\cdot }\boldsymbol{\nabla}\boldsymbol{U})' +\bar{\boldsymbol{u}}\boldsymbol{\cdot }\boldsymbol{\nabla}\boldsymbol{U} &={-}\boldsymbol{\nabla} p'+(\boldsymbol{B}\boldsymbol{\cdot }\boldsymbol{\nabla}\boldsymbol{b}'+\boldsymbol{b}' \boldsymbol{\cdot }\boldsymbol{\nabla}\boldsymbol{B})'\nonumber\\ &\quad +\bar{\boldsymbol{b}}\boldsymbol{\cdot }\boldsymbol{\nabla}\boldsymbol{B}+ Re^{{-}1} \nabla^2 \boldsymbol{u}', \end{align}

(2.8) \begin{align} \frac{\partial \boldsymbol{b}'}{\partial t} + ( \boldsymbol{U} \boldsymbol{\cdot } \boldsymbol{\nabla} \boldsymbol{b}' + \boldsymbol{u}' \boldsymbol{\cdot } \boldsymbol{\nabla} \boldsymbol{B} )' + \bar{\boldsymbol{u}}\boldsymbol{\cdot }\boldsymbol{\nabla} \boldsymbol{B} &= ( \boldsymbol{B} \boldsymbol{\cdot } \boldsymbol{\nabla} \boldsymbol{u}' + \boldsymbol{b}' \boldsymbol{\cdot } \boldsymbol{\nabla}\boldsymbol{U} )' \nonumber\\ &\quad + \bar{\boldsymbol{b}} \boldsymbol{\cdot }\boldsymbol{\nabla} \boldsymbol{U} + Rm^{{-}1} \nabla^2 \boldsymbol{b}' . \end{align}

It can be seen that each of $\boldsymbol {u}'$ and $\boldsymbol {b}'$ is subject to forcing from both $\bar {\boldsymbol {u}}$ and $\bar {\boldsymbol {b}}$. It is assumed that the long-time solutions of (2.7) and (2.8) are (statistically) unique, consistent with the basic state being a (statistical) attractor. If (2.7) and (2.8) are formally solved, and the results substituted into (2.5) and (2.6), then we obtain at leading order, after dropping the diffusive terms, which are subdominant at small wavenumbers:

(2.9)$$\begin{gather} \frac{\partial \bar{u}_i}{\partial T} +\frac{\partial}{\partial X_j}(\varGamma^{(1)}_{ijl}\bar{u}_l+\varGamma^{(2)}_{ijl}\bar{b}_l) ={-}\frac{\partial\bar{p}}{\partial X_i} , \end{gather}$$

(2.10)$$\begin{gather}\frac{\partial \bar{{b}}_i}{\partial T} = \epsilon_{ijk}\frac{\partial}{\partial X_j}( \alpha^{(1)}_{kl} \bar{{b}}_l +\alpha^{(2)}_{kl}\bar{u}_l ) . \end{gather}$$

The quantities $\boldsymbol {\varGamma }^{(i)}$, $\boldsymbol {\alpha }^{(i)}$ are tensors, depending only on the properties of the basic state $\boldsymbol {U}$, $\boldsymbol {B}$ and so ultimately on those of $\boldsymbol {F}$ (to be more precise, $\boldsymbol {\varGamma }^{(1)}$ and $\boldsymbol {\alpha }^{(2)}$ are true tensors, whereas $\boldsymbol {\varGamma }^{(2)}$ and $\boldsymbol {\alpha }^{(1)}$ are pseudo-tensors). The $\boldsymbol {\varGamma }^{(i)}$ are always symmetric in their first two indices, whereas the $\boldsymbol {\alpha }^{(i)}$ have no symmetries in general. It is important to note that since we are considering forced turbulence in a fixed frame of reference, our system is not Galilean invariant; in a Galilean-invariant system, the tensors $\boldsymbol {\varGamma }^{(1)}$ and $\boldsymbol {\alpha }^{(2)}$ would vanish.

Two of the mean quantities appearing in (2.9) and (2.10) are familiar: $\boldsymbol {\alpha }^{(1)}$ is the well-known dynamo $\alpha$-effect, which has been extensively discussed in the literature; $\boldsymbol {\varGamma }^{(1)}$ can be recognised as the AKA effect introduced by Frisch et al. (Reference Frisch, She and Sulem1987). The other terms, which couple the two mean field equations, and were first introduced by Courvoisier et al. (Reference Courvoisier, Hughes and Proctor2010a), have not been widely studied. It is interesting to note that reflectionally symmetric turbulence leads to constraints on $\boldsymbol {\alpha }^{(1)}$ and $\boldsymbol {\varGamma }^{(2)}$, being pseudo-tensors (e.g. the symmetric part of $\boldsymbol {\alpha }^{(1)}$ vanishes), whereas there are no such constraints on $\boldsymbol {\alpha }^{(2)}$ and $\boldsymbol {\varGamma }^{(1)}$.

As noted in the introduction, it is of interest to compare our formulation with other approaches to obtaining mean field tensors in turbulent MHD. Yoshizawa (Reference Yoshizawa1990) (see also Yokoi Reference Yokoi2023) considered a two-scale approach to MHD turbulence, adopting the two-scale direct-interaction approximation (TSDIA), retaining the effects both of the mean quantities and their derivatives, thus capturing effects such as eddy diffusivities, which we do not consider here. Yoshizawa (Reference Yoshizawa1990) considered a system that was Galilean invariant and, as such, the tensors $\boldsymbol {\varGamma }^{(1)}$ and $\boldsymbol {\alpha }^{(2)}$ vanish. Moreover, under the TSDIA, $\boldsymbol {\varGamma }^{(2)}$ also vanishes. Schrinner et al. (Reference Schrinner, Rädler, Schmitt, Rheinhardt and Christensen2005) adopted a very different approach. They extracted the fluctuating velocity from a stationary MHD turbulent state – possibly one with a mean magnetic field – and then solved the fluctuating induction equation after the imposition of so-called ‘test fields’, which may be functions of position. They then calculated the resulting electromotive force (e.m.f.), which they related to the $\alpha _{ij}$ and $\beta _{ijk}$ tensors of mean field MHD. Although their starting point is a saturated MHD state, they did not consider the equation for the mean velocity.

2.2. Determining the tensors ${\alpha}^{(i)}$ and ${\varGamma}^{(i)}$ under the ‘short-sudden’ approximation

Obtaining analytic expressions for the mean quantities $\bar {\boldsymbol {u}}$, $\bar {\boldsymbol {b}}$ is not, in general, possible. However, there are two special cases in which progress can be made. One is when both $Re$ and $Rm$ are small, which has already been considered in Rädler & Brandenburg (Reference Rädler and Brandenburg2010) and Courvoisier et al. (Reference Courvoisier, Hughes and Proctor2010b). The other, which we consider here, is when the correlation time of the turbulence is short and one may adopt the ‘short-sudden’ approximation (see e.g. Krause & Rädler Reference Krause and Rädler1980). In this case, (2.7) and (2.8) are approximated at leading order by

(2.11a,b)\begin{equation} \boldsymbol{u}'=\tau_1(\bar{\boldsymbol{b}}\boldsymbol{\cdot }\boldsymbol{\nabla}\boldsymbol{B}-\bar{\boldsymbol{u}}\boldsymbol{\cdot }\boldsymbol{\nabla}\boldsymbol{U}),\quad \boldsymbol{b}'=\tau_2(\bar{\boldsymbol{b}}\boldsymbol{\cdot }\boldsymbol{\nabla}\boldsymbol{U}-\bar{\boldsymbol{u}}\boldsymbol{\cdot }\boldsymbol{\nabla}\boldsymbol{B}), \end{equation}

where $\tau _1$ and $\tau _2$ are, respectively, the correlation times for the fluctuating flow and field. Note that $\boldsymbol {u}'$ is automatically solenoidal, and so the pressure gradient does not enter this approximation. Substituting for $\boldsymbol {u}'$ and $\boldsymbol {b}'$ into (2.5) and (2.6) allows us to calculate the mean field tensors. It is straightforward to show that $\boldsymbol {\varGamma }^{(1)}$ vanishes in this special case, since

(2.12)\begin{align} \varGamma_{ijl}^{(1)} &={-}\tau_1 \left( \overline{U_j \frac{\partial U_i}{\partial x_l}} + \overline{\frac{\partial U_j}{\partial x_l} U_i }\right) + \tau_2 \left( \overline{B_j \frac{\partial B_i}{\partial x_l}} + \overline{\frac{\partial B_j}{\partial x_l} B_i }\right) \end{align}

(2.13)\begin{align} &={-}\tau_1 \overline{\frac{\partial \left(U_i U_j \right)}{\partial x_l}} + \tau_2 \overline{\frac{\partial \left( B_i B_j \right)} {\partial x_l}} =0. \end{align}

The elements of the tensors $\boldsymbol {\alpha }^{(1)}$, $\boldsymbol {\alpha }^{(2)}$ and $\boldsymbol {\varGamma }^{(2)}$ are given by

(2.14)$$\begin{gather} \alpha_{il}^{(1)} = \epsilon_{ipq} \left( \tau_2\overline{U_p \frac{\partial U_q}{\partial x_l}} - \tau_1 \overline{B_p \frac{\partial B_q}{\partial x_l}} \right), \end{gather}$$

(2.15)$$\begin{gather}\alpha^{(2)}_{il}=(\tau_1-\tau_2)\epsilon_{ipq}\gamma_{pql}, \end{gather}$$

(2.16)$$\begin{gather}\varGamma^{(2)}_{ijl}= (\tau_2+\tau_1)(\gamma_{ijl}+\gamma_{jil}), \end{gather}$$

where

(2.17)\begin{equation} \gamma_{ijl} = \overline{U_i\frac{\partial B_j}{\partial x_l}} ={-} \overline{B_j\frac{\partial U_i}{\partial x_l}} = \frac{1}{2} \left( \overline{U_i\frac{\partial B_j}{\partial x_l}} - \overline{B_j\frac{\partial U_i}{\partial x_l}} \right). \end{equation}

It should be noted that under the short-sudden approximation, the equations for the fluctuating flow and field, (2.11a,b), take a simplified form. As a consequence, $\boldsymbol {\alpha }^{(1)}$ contains only products of $\boldsymbol {U}$ with itself and $\boldsymbol {B}$ with itself – in general, however, there is no reason to expect that $\boldsymbol {\alpha }^{(1)}$ will not contain products of $\boldsymbol {U}$ and $\boldsymbol {B}$. In a similar vein, when $\tau _1=\tau _2$, the symmetry of (2.11a,b) leads to the vanishing of the $\boldsymbol {\alpha }^{(2)}$ tensor, although $\tau _1$ and $\tau _2$ cannot be calculated a priori in general. There is an analogous situation in the approximated system introduced by Rädler & Brandenburg (Reference Rädler and Brandenburg2010), in which $\boldsymbol {\alpha }^{(2)}$ vanishes when $Re=Rm$.

The expression for $\boldsymbol {\alpha }^{(1)}$, which first appeared in Pouquet, Frisch & Léorat (Reference Pouquet, Frisch and Léorat1976), is the extension of the classical $\alpha$-effect from hydrodynamic to magnetohydrodynamic basic states; the velocity and magnetic fields are clearly on an equal footing, as indicated in the introduction. It is worth noting that this form of $\boldsymbol {\alpha }^{(1)}$ is sometimes interpreted as a ‘quenching’ mechanism for the $\alpha$-effect as a dynamo evolves to finite amplitude; there are, however, problems with such an interpretation, as noted by Proctor (Reference Proctor2003). The two terms that couple the basic state velocity and magnetic fields ($\boldsymbol {\alpha }^{(2)}$ and $\boldsymbol {\varGamma }^{(2)}$) rely for their existence on different parts of $\gamma _{ijl}$. In particular, since $\boldsymbol {\varGamma }^{(2)}$ is symmetric in its first two arguments, it will vanish if the statistics of the basic state are isotropic. The coefficients of the $\gamma _{ijl}$ tensor are a priori almost unconstrained, although since $\boldsymbol {U}$ and $\boldsymbol {B}$ are solenoidal, we do have $\gamma _{iji} = \gamma _{ijj} = 0$. It is straightforward to show that for isotropic turbulence, $\alpha ^{(2)}_{il} = \tilde \alpha \delta _{il}$, where

(2.18)\begin{equation} \tilde{\alpha}={-}\left( \frac{\tau_1-\tau_2}{3} \right) \overline{\boldsymbol{U} \boldsymbol{\cdot } \boldsymbol{\nabla} \times \boldsymbol{B}} ={-}\left( \frac{\tau_1-\tau_2}{3} \right) \overline{\boldsymbol{B} \boldsymbol{\cdot } \boldsymbol{\nabla} \times \boldsymbol{U}}. \end{equation}

It is of interest to note that even if the traditional $\alpha$-effect ($\boldsymbol {\alpha }^{(1)}$) vanishes, it still seems that there can be a long-wavelength instability induced by the coupling terms. Seeking solutions proportional, for example, to $\exp ({\rm i}(KZ + \omega T))$, so that $\bar {u}_3=\bar {{b}}_3=0$, we have the coupled two-dimensional algebraic equations (where $G^{(1)}_{ij}=\gamma _{i3j}$, $G^{(2)}_{ij}=\gamma _{3ij}$),

(2.19a,b)\begin{equation} \omega \bar{\boldsymbol{u}}={-}(\tau_1+\tau_2) K({\boldsymbol{G}^{(1)}+ \boldsymbol{G}^{(2)})}\bar{\boldsymbol{b}}, \quad \omega \bar{\boldsymbol{b}}= (\tau_1 - \tau_2) K({\boldsymbol{G}^{(1)}- \boldsymbol{G}^{(2)}})\bar{\boldsymbol{u}},\end{equation}

and so $\omega ^2$ is an eigenvalue of the matrix

(2.20)\begin{equation} (\tau_2^2 -\tau_1^2)K^2 ({\boldsymbol{G}^{(1)}+\boldsymbol{G}^{(2)}})({\boldsymbol{G}^{(1)}-\boldsymbol{G}^{(2)}}). \end{equation}

Unless both eigenvalues $\omega ^2$ are real and positive, the basic state will be unstable. It therefore appears that there is a new generic mechanism for long-wavelength instability, relying on coupling between the mean momentum and induction equations and on anisotropy in the basic flow statistics.

2.3. From formalism to calculation

In theory, numerical calculation of the $\boldsymbol {\varGamma }^{(i)}$ and $\boldsymbol {\alpha }^{(i)}$ tensors in (2.9) and (2.10) is straightforward, albeit computationally expensive. Following the standard mean field prescription, the fluctuating fields $\boldsymbol {u}'$ and $\boldsymbol {b}'$ are calculated from the linear equations (2.7) and (2.8), solved in concert with the nonlinear equations (2.1) and (2.2), following the imposition of uniform mean velocity and magnetic fields, $\bar {\boldsymbol {u}}$ and $\bar {\boldsymbol {b}}$. From the fluctuating fields, we construct the linearised mean total stress and the linearised mean e.m.f., defined by

(2.21)$$\begin{gather} R_{ij} = \overline{(U_j u_i'} + \overline{u'_j U_i)} - \overline{(B_j b_i'} + \overline{b'_j B_i)}, \end{gather}$$

(2.22)$$\begin{gather}{\boldsymbol{\mathcal{E}}} = \overline{\boldsymbol{U}\times\boldsymbol{b}'} + \overline{\boldsymbol{u}'\times\boldsymbol{B}}. \end{gather}$$

Varying the directions of $\bar {\boldsymbol {u}}$ and $\bar {\boldsymbol {b}}$ allows calculation of all the tensorial coefficients, after spatial and temporal averaging of the various quadratic interactions, with the proviso that the requisite averages exist.

By contrast, when the averages do not exist, it is necessary, as discussed in the introduction, to consider fully nonlinear responses to the imposition of mean fields and flows. We restore the nonlinear terms into the perturbation equations (cf. (2.3) and (2.4)), which now take the form

(2.23)$$\begin{gather} \frac{\partial \boldsymbol{u}}{\partial t} +\boldsymbol{U}\boldsymbol{\cdot }\boldsymbol{\nabla}\boldsymbol{u} +\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{\nabla}\boldsymbol{U} +\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{\nabla}\boldsymbol{u} ={-}\boldsymbol{\nabla} p + \boldsymbol{B}\boldsymbol{\cdot }\boldsymbol{\nabla}\boldsymbol{b} + \boldsymbol{b}\boldsymbol{\cdot }\boldsymbol{\nabla}\boldsymbol{B} + \boldsymbol{b}\boldsymbol{\cdot }\boldsymbol{\nabla}\boldsymbol{b} + Re^{{-}1} \nabla^2 \boldsymbol{u}, \end{gather}$$

(2.24)$$\begin{gather}\frac{\partial \boldsymbol{b}}{\partial t} +\boldsymbol{U}\boldsymbol{\cdot }\boldsymbol{\nabla}\boldsymbol{b} +\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{\nabla}\boldsymbol{B} +\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{\nabla}\boldsymbol{b} = \boldsymbol{B}\boldsymbol{\cdot }\boldsymbol{\nabla}\boldsymbol{u} +\boldsymbol{b}\boldsymbol{\cdot } \boldsymbol{\nabla}\boldsymbol{U} + \boldsymbol{b}\boldsymbol{\cdot }\boldsymbol{\nabla}\boldsymbol{u} + Rm^{{-}1} \nabla^2 \boldsymbol{b} , \end{gather}$$

where, again, $\boldsymbol {u}=\bar {\boldsymbol {u}}+\boldsymbol {u}'$, etc. Equations (2.23) and (2.24) are solved in concert with (2.1) and (2.2). The mean total stress and the mean e.m.f. now take the form

(2.25)$$\begin{gather} R_{ij} = \overline{(U_j u_i'} + \overline{u'_j U_i} + \overline{u_i' u_j')} - \overline{(B_j b_i'} + \overline{b'_j B_i} + \overline{b_i' b_j')}, \end{gather}$$

(2.26)$$\begin{gather}{\boldsymbol{\mathcal{E}}} = \overline{\boldsymbol{U}\times\boldsymbol{b}'} + \overline{\boldsymbol{u}'\times\boldsymbol{B}} + \overline{\boldsymbol{u}' \times \boldsymbol{b}'}. \end{gather}$$

The MHD equations are solved on a cubic, $2{\rm \pi}$-periodic domain, using a parallelised, dealiased, pseudospectral code. Time advancement of the diffusive terms is carried out exactly using an integrating factor, and the remaining terms are treated using a third-order Runge–Kutta scheme. For a detailed description of the numerical method, see Cattaneo, Emonet & Weiss (Reference Cattaneo, Emonet and Weiss2003). The simulations with the largest Reynolds numbers have a spatial resolution of $128^3$ grid points. Very long time evolutions are required in order to obtain meaningful averages; a single simulation can consume on the order of $10^4$ core hours. A series of simulations with different orientations and strengths of the mean fields is required in order to measure $\boldsymbol {\alpha }^{(i)}$ and $\boldsymbol {\varGamma }^{(i)}$.

3.1. When the linear theory works well

In this subsection, we investigate long-wavelength instabilities to the basic state generated by forcings $\boldsymbol {F}_1$ and $\boldsymbol {F}_2$ at low values of $Re$ and fairly low values of $Rm$ (though high enough to support dynamo action, so that the basic states are magnetohydrodynamic).

First we consider forcing $\boldsymbol {F}_1$, with $Re=Rm=12$; the basic MHD state for these particular parameters has been studied in some detail by Galanti, Sulem & Pouquet (Reference Galanti, Sulem and Pouquet1992). In the absence of magnetic field, the 1:1:1 ABC flow is hydrodynamically stable; it does though support dynamo action for $Rm=12$, and so an initially weak field is amplified until a stationary MHD state is attained in which the flow varies from its hydrodynamic state. For these parameter values, the basic MHD state has regular periodic oscillations, with mean values of $\langle U^2 \rangle \approx 1.67$ and $\langle B^2 \rangle \approx 0.057$, where $\langle {\cdot } \rangle$ denotes a volume and time average. This new state is not isotropic, but has a preferred direction, which depends on initial conditions; this state corresponds to the long-time solution shown in figure 3 of Galanti et al. (Reference Galanti, Sulem and Pouquet1992). Representative two-dimensional slices of the flow and field of the basic MHD state are shown in figure 1; here the preferred direction is along the $x$-axis. Interestingly, although the flow differs from its kinematic (hydrodynamic) form, it is, nonetheless, almost maximally helical, i.e. with $\langle\boldsymbol {U} \boldsymbol{\cdot } \boldsymbol {\nabla } \times \boldsymbol {U}\rangle/\sqrt {\langle U^2 \rangle \langle |\boldsymbol {\nabla } \times \boldsymbol {U}|^2 \rangle } \approx 1$.

Figure 1. Snapshots of representative slices of the basic state flow and field for forcing $\boldsymbol {F}_1$, with $Re=Rm=12$. From left to right: (a,b,c) $U_x(y,z)$, $U_y(x,z)$, $U_z(x,y)$; (d,e,f) $B_x(y,z)$, $B_y(x,z)$, $B_z(x,y)$.

The flow is perturbed, in turn, by the imposition of kinematic uniform fields and flows in the $x$, $y$ and $z$ directions, allowing us to calculate $R_{ij}$ and ${\boldsymbol{\mathcal{E}}}$ given by (2.21) and (2.22). By imposing magnetic fields, we can determine $\boldsymbol {\alpha }^{(1)}$ and $\boldsymbol {\varGamma }^{(2)}$, whereas from the imposed flows we can determine $\boldsymbol {\alpha }^{(2)}$ and $\boldsymbol {\varGamma }^{(1)}$. Here, the linear perturbations are also periodic in time and hence the quadratic quantities needed to calculate the mean field tensors must have well-defined averages. When the mean values are non-zero, convergence to the mean is rapid; however, for certain cases where the mean value is zero, some components of the Reynolds stress converge in amplitude faster than others. After imposing mean magnetic fields, the tensor $\boldsymbol {\alpha }^{(1)}$ is calculated as

(3.3)\begin{equation} \boldsymbol{\alpha}^{(1)} = \left( \begin{array}{@{}c@{\quad}c@{\quad}c@{}} -1.22 & 0 & 0 \\ 0 & -0.043 & 0 \\ 0 & 0 & -0.043 \end{array} \right), \end{equation}

where the symmetries of the tensor indicate that the basic state for this run is isotropic in the $yz$-plane. The symmetries of the flow and field are such that all components of the tensor $\boldsymbol {\varGamma }^{(2)}$ are zero. On imposing a mean flow, all components of the tensor $\boldsymbol {\alpha }^{(2)}$ are similarly zero; the only non-zero components of $\boldsymbol {\varGamma }^{(1)}$ are given by

(3.4)\begin{equation} \varGamma^{(1)}_{132} = \varGamma^{(1)}_{312} ={-}\varGamma^{(1)}_{123} ={-}\varGamma^{(1)}_{213} =0.058. \end{equation}

Since both the tensors $\boldsymbol {\alpha }^{(2)}$ and $\boldsymbol {\varGamma }^{(2)}$ are zero, as a consequence of the symmetries of the basic state, the mean field equations (2.9) and (2.10) are decoupled. Any instability that results is therefore either purely magnetic or purely hydrodynamic, although the basic flow of course depends crucially on both the flow and field.

To obtain insight into the magnetic instability, we consider a scaled form of $\boldsymbol {\alpha }^{(1)}$ given by (3.3), namely

(3.5)\begin{equation} \boldsymbol{\alpha}^{(1)} = \left( \begin{array}{@{}c@{\quad}c@{\quad}c@{}} -1 & 0 & 0 \\ 0 & -\varepsilon & 0 \\ 0 & 0 & -\varepsilon \end{array} \right). \end{equation}

On seeking solutions of the form $\exp ({\rm i} \boldsymbol {K} \boldsymbol{\cdot } \boldsymbol {X} + S T)$, (2.10) leads to the following expression for the growth rate $S$:

(3.6)\begin{equation} S ={\pm} \left( \varepsilon K_H^2 + \varepsilon^2 K_X^2 \right)^{1/2}, \end{equation}

where $K_H^2 = K_Y^2 + K_Z^2$. The next order correction to (3.6) is $O(|\boldsymbol {K}|^2)$ and comprises molecular diffusion together with eddy effects that would result from inclusion of the first derivatives of $\bar {\boldsymbol {u}}$ and $\bar {\boldsymbol {b}}$ in the above analysis. For illustrative purposes, if we just retain the molecular terms, the extension of (3.6) becomes

(3.7)\begin{equation} S ={\pm} \left( \varepsilon K_H^2 + \varepsilon^2 K_X^2 \right)^{1/2} - Rm^{{-}1} \left( K_H^2 + K_X^2\right). \end{equation}

For $\varepsilon < 1$, the maximum growth rate is $\varepsilon Rm /4$, when $K_X=0$ and $K_H=\sqrt {\varepsilon } Rm/2$. Note that this dynamo mechanism is very different from that resulting from two-dimensional flows such as first studied by Roberts (Reference Roberts1972), which may be thought of as the case of very large $\varepsilon$.

To analyse the flow instability, we consider a scaled form of $\boldsymbol {\varGamma }^{(1)}$ given by (3.4), with the following non-zero entries:

(3.8)\begin{equation} \varGamma^{(1)}_{132} = \varGamma^{(1)}_{312} ={-}\varGamma^{(1)}_{123} ={-}\varGamma^{(1)}_{213} = 1. \end{equation}

The growth rate resulting from (2.9), again with the inclusion of the dissipative term for illustrative purposes, is then given by

(3.9)\begin{equation} S ={\pm} K_X \left( \frac{K_X^2 - K_H^2}{K_H^2 + K_X^2}\right)^{1/2} - Re^{{-}1} \left( K_H^2 + K_X^2\right). \end{equation}

Note that there is instability provided that $K_X^2 > K_H^2$ and the wavenumbers are sufficiently small that dissipative effects can be ignored.

We have also considered the nonlinear theory for these parameters and verified, as expected, that there is good agreement with the linear theory, provided that the imposed fields and flows are sufficiently small. This point is illustrated in figure 2, which shows the non-zero components of the e.m.f. ${\boldsymbol{\mathcal{E}}}$, now calculated via the nonlinear prescription (2.26), for various orientations and strengths of the imposed mean field $\bar {\boldsymbol {b}}$. The figure shows that the e.m.f. depends linearly on $\bar {\boldsymbol {b}}$ until $|\bar {\boldsymbol {b}}|$ reaches a value of $O(Rm^{-1})$.

Figure 2. The nonlinear e.m.f.s for different orientations of the mean field (along (a) $\hat {\boldsymbol {x}}$, (b) $\hat {\boldsymbol {y}}$, (c) $\hat {\boldsymbol {z}}$) for the nonlinear perturbations for forcing $F_1$, with $Re=Rm=12$. The dashed straight line shows the linear prediction.

To explore a less symmetrical situation than that provided by the forcing $\boldsymbol {F}_1$, with the aim of coupling the mean field equations (2.9) and (2.10), we now consider the forcing $\boldsymbol {F}_2$ with $Re=1$ and $Rm=12$. Representative two-dimensional slices of the flow and field of the basic state are shown in figure 3. For this case, the mean field tensors take the form:

(3.10a,b) \begin{align} \boldsymbol{\alpha}^{(1)} &= \left(\begin{array}{@{}c@{\quad}c@{\quad}c@{}} -0.31 & 0 & 0.339 \\ 0 & 0.418 & 0 \\ 0.71 & 0 & -0.403 \end{array} \right), \quad \boldsymbol{\alpha}^{(2)} = \left(\begin{array}{@{}c@{\quad}c@{\quad}c@{}} 0 & -0.083 & 0 \\ 0.0208 & 0 & 0.187 \\ 0 & 0.074 & 0 \end{array} \right) , \end{align}

(3.11a,b) \begin{align} \boldsymbol{\varGamma}_{ij1}^{(1)} &= \left( \begin{array}{@{}c@{\quad}c@{\quad}c@{}} 0 & -0.20 & 0 \\ -0.20 & 0 & -0.27 \\ 0 & -0.27 & 0 \end{array} \right), \quad \boldsymbol{\varGamma}_{ij1}^{(2)} = \left( \begin{array}{@{}c@{\quad}c@{\quad}c@{}} 2.27 & 0 & 1.64 \\ 0 & 3.41 & 0 \\ 1.64 & 0 & 2.84 \end{array} \right), \end{align}

(3.12a,b) \begin{align} \boldsymbol{\varGamma}_{ij2}^{(1)} &= \left( \begin{array}{@{}c@{\quad}c@{\quad}c@{}} -0.124 & 0 & 0.024 \\ 0 & -0.190 & 0 \\ 0.024 & 0 & 0.188 \end{array} \right), \quad \boldsymbol{\varGamma}_{ij2}^{(2)} = \left( \begin{array}{ccc} 0 & -0.128 & 0 \\ -0.128 & 0 & 0.081 \\ 0 & 0.081 & 0 \end{array} \right), \end{align}

(3.13a,b) \begin{align} \boldsymbol{\varGamma}_{ij3}^{(1)} &= \left( \begin{array}{@{}c@{\quad}c@{\quad}c@{}} 0 & 0.614 & 0 \\ 0.614 & 0 & 0.254 \\ 0 & 0.254 & 0 \end{array} \right), \quad \boldsymbol{\varGamma}_{ij3}^{(2)} = \left( \begin{array}{@{}c@{\quad}c@{\quad}c@{}} -2.24 & 0 & -0.738 \\ 0 & -2.92 & 0 \\ -0.738 & 0 & -1.70 \end{array} \right) . \end{align}

Equations (2.9) and (2.10) are now coupled, with all mean field tensors being non-zero. We can investigate the linear stability of the basic state by seeking solutions to (2.9) and (2.10) proportional to $\exp ({\rm i}(\boldsymbol {K} \boldsymbol{\cdot } \boldsymbol {X} + \omega T))$. On eliminating the pressure, these take the form

(3.14)$$\begin{gather} \omega \bar{u}_i ={-} K_j \left( \delta_{ip} - \frac{ K_i K_p}{K^2} \right) \left( \varGamma^{(1)}_{pjl}\bar{u}_l + \varGamma^{(2)}_{pjl}\bar{b}_l \right) \equiv Q^{(1)}_{il} \bar{u}_l + Q^{(2)}_{il} \bar{{b}}_l, \end{gather}$$

(3.15)$$\begin{gather}\omega \bar{{b}}_i = \epsilon_{ijp} K_j \left( \alpha^{(1)}_{pl} \bar{{b}}_l + \alpha^{(2)}_{pl}\bar{u}_l \right) \equiv R^{(1)}_{il} \bar{{b}}_l + R^{(2)}_{il} \bar{u}_l. \end{gather}$$

Thus, $\omega$ is an eigenvalue of the $6 \times 6$ real matrix

(3.16) \begin{equation} \left[\!\begin{array}{c|c} \boldsymbol{Q}^{(1)} & \boldsymbol{Q}^{(2)} \\ \hline \boldsymbol{R}^{(2)} & \boldsymbol{R}^{(1)^{\vphantom{A}}} \\ \end{array}\!\right]. \end{equation}

Figure 3. Snapshots of slices of the basic state flow and field for forcing $\boldsymbol {F}_2$ with $Re=1$, $Rm=12$. From left to right: (a,b,c) $U_x(y,z)$, $U_y(x,z)$, $U_z(x,y)$; (d,e,f) $B_x(y,z)$, $B_y(x,z)$, $B_z(x,y)$.

Since diffusion becomes negligible when $|\boldsymbol {K}|$ is sufficiently small, we see that in this case there will be instability for any chosen $\boldsymbol {K}$ if there is a complex eigenvalue of the above matrix. By contrast, if all the eigenvalues are real, then there is stability owing to the effect of diffusion. Writing $\boldsymbol {K} = |\boldsymbol {K}| (\sin \theta \cos \phi, \sin \theta \sin \phi, \cos \theta )$, we can find the eigenvalues $\omega$ for each $\theta$ and $\phi$. Figure 4(a) is a Mollweide projection of the largest value of $| \textrm {Im} (\omega )|/|\boldsymbol {K}|$ as a function of $\theta$ and $\phi$, showing that there is instability for a wide range of wavenumber directions. We can compare this figure with those that would arise by considering the stability problem when retaining the $\boldsymbol {\alpha }^{(1)}$ and $\boldsymbol {\varGamma }^{(1)}$ tensors in isolation, as shown in figure 4(b,c). It is of interest to note that the full problem has a wider range of instability than the union of the unstable regions for the matrices $\boldsymbol {Q}^{(1)}$ and $\boldsymbol {R}^{(1)}$, thus showing that the cross terms ($\boldsymbol {Q}^{(2)}$ and $\boldsymbol {R}^{(2)}$) can indeed lead to additional instabilities.

Figure 4. Mollweide $(\theta, \phi )$ projections of the scaled growth rate (the largest value of $| \textrm {Im} (\omega )|/|\boldsymbol {K}|$), for the system described by (3.14), (3.15); longitude is in the range $-{\rm \pi} < \phi < {\rm \pi}$, latitude in the range $0 < \theta < {\rm \pi}$ (the equator is $\theta ={\rm \pi} /2$): (a) growth rate for the full system, governed by the matrix (3.16); (b) growth rate for the system in which only the $\boldsymbol {\alpha }^{(1)}$ effect operates ($\boldsymbol {Q}^{(1)}$, $\boldsymbol {Q}^{(2)}$, $\boldsymbol {R}^{(2)}$ all zero); (c) growth rate for the system in which only the $\boldsymbol {\varGamma }^{(1)}$ effect operates ($\boldsymbol {Q}^{(2)}$, $\boldsymbol {R}^{(1)}$, $\boldsymbol {R}^{(2)}$ all zero). In each plot, the maximum of the scaled growth rate is normalised to unity; the numerical values of the maximum scaled growth rate are (a) $0.338$, (b) $0.267$ and (c) $0.338$.

3.2. When the linear theory fails

The flows considered in § 3.1, with reasonably small values of $Re$ and $Rm$, are fairly simple in space and time, with well-defined values of the mean field tensors when nonlinear terms are neglected. In this subsection, we again consider flows driven by the forcing $\boldsymbol {F}_1$, but with higher values of $Re$ and $Rm$; specifically we consider the case of $Re = Rm = 100$. Figure 5 shows $\langle U^2 \rangle$ and $\langle B^2 \rangle$ vs time for the basic state; the flows and fields are now disordered, in both space and time.

Figure 5. Plots of (a) $\langle U^2 \rangle$ and (b) $\langle B^2 \rangle$ vs time for forcing $\boldsymbol {F}_1$ with $Re = Rm = 100$.

On imposing uniform fields and flows, and in stark contrast to the flows of § 3.1, we find that the solutions to the linear equations (2.7) and (2.8) increase without bound. This is demonstrated by figure 6, which shows the exponential growth, on average, of the e.m.f.s and combined Reynolds and Maxwell stresses calculated from expressions (2.21) and (2.22). As such, the linear theory can give no guide to the mean field coefficients. We thus have to consider the full nonlinear MHD equations (2.23) and (2.24) – for which perturbations are guaranteed to be bounded – and to see whether the average response depends linearly on the imposed mean fields and flows when these are sufficiently small. To provide some insight into the nature of the response of a complicated system to an applied perturbation, we delay to § 5 the discussion of fully nonlinear perturbations in mean field MHD, in order to consider an ostensibly much simpler problem, namely the nature of the response to small symmetry-breaking perturbations of one-dimensional maps.

Figure 6. Components of the mean e.m.f. (a–c) and mean total stress (d–i), calculated from the linear theory, with forcing $\boldsymbol {F}_1$ and $Re=Rm=100$, after the imposition of a magnetic field of unit magnitude in the $y$-direction.

4.1. The invariant measure

The invariant measure (or probability density) $\mu(x)$ of a map $f(x)$ is specified in the following way. First, we need to know the effect of the map on intervals in the domain. For every point $x$ there is at least one pre-image $\tilde {x}$ such that ${f(\tilde {x})=x}$. The number of pre-images depends on the map, but for a chaotic map there are at least two. An interval $\mathrm {d}\kern0.7pt\tilde {x}$ at $\tilde {x}$ is mapped into an interval ${{\rm d}\kern0.7pt x}=|\,f'(\tilde {x})| \,\mathrm {d}\kern0.7pt\tilde {x}$ at $x$. Thinking of $\mu (x)$ as a (relative) ‘density’ of points in the domain, we can see that the points in the interval ${{\rm d}\kern0.7pt x}$ at $x$ derive from the sum of the quantities $(\mu (\tilde {x}_i)/|\,f'(\tilde {x}_i)|) \,\mathrm {d}\kern0.7pt\tilde {x}$ over the pre-images $\tilde {x}_i$. Since the measure or density is invariant under the map, we have the following functional equation (the Frobenius–Perron equation) determining $\mu (x)$, where the map is defined on the range $-1\leqslant x\leqslant 1$:

(4.1a,b)\begin{equation} \mu(x)=\mathop{\sum}\limits_i \frac{\mu(\tilde{x}_i)}{|\,f(\tilde{x}_i)|}\quad{\rm and}\quad\int_{{-}1}^1 \mu(x)\, {{\rm d}\kern0.7pt x}=1. \end{equation}

4.2. Maps with linear response (i): the ‘cubic tent map’

In general, $\mu (x)$ is very complicated and non-differentiable in nature; see, for example, the plots of $\mu (x)$ for the logistic map in Gottwald, Wormell & Wouters (Reference Gottwald, Wormell and Wouters2016). However, there are some special problems for which both $\mu$ and its dependence on the parameters of the map can be calculated explicitly. As such an example, here we consider the one-parameter family $f(x,a)$ defined on $-1\leqslant x\leqslant 1$ (see figure 7):

(4.2)\begin{align} \left.\begin{array}{@{}l@{\quad}c@{}} f(x,a) ={-}2-3x & -1 \leqslant x \leqslant -1/3, \\ \phantom{f(x,a)}= 3x & -1/3 \leqslant x \leqslant 0, \\ \phantom{f(x,a)}= ax & 0 \leqslant x \leqslant 1/a, \\ \phantom{f(x,a)}=\dfrac{1+a-2ax}{a-1} & 1/a \leqslant x \leqslant 1. \end{array}\right\} \end{align}

When $a=3$ the map is antisymmetric, with the pre-images of $x$ being $x/3$, $(2-x)/3$ and $- (2+x)/3$. Since the gradients of $f(x,a=3)$ at each pre-image have modulus $3$, (4.1a,b) is solved by the constant value $\mu (x)=1/2$. Thus $\langle x\rangle \equiv \int _{-1}^1 x\mu (x)\, {{\rm d} \kern0.7pt x}=0$.

Figure 7. The map $f(x,a)$ (the cubic tent map, defined by (4.2)), for $a=3$ (red) and $a=5$ (green); the blue line denotes the invariant measure $\mu (x)$ for $a=5$.

For $a \ne 3$, the (anti-)symmetry of the map is broken, as shown in figure 7. There are now two different (absolute) gradients of the map for positive $x$, and yet a third for $x<0$, and so $\mu$ is no longer constant. However, it turns out that $\mu$ is still piecewise constant, with $\mu =\mu _+$, $x> 0$; $\mu =\mu _-$, $x<0$. For each $x$, there are three pre-images $\tilde {x}_{i}$: for $x>0$, $\tilde {x}_{1}<0$ and $0<\tilde {x}_{2}< a^{-1}<\tilde {x}_{3}<1$, whereas for $x<0$, $\tilde {x}_{1,2}<0$ and $\tilde {x}_{3}>a^{-1}$. The (moduli of the) gradients in $x<0$, $0< x< a_{-1}$, $a_{-1}< x<1$ are $3$, $a$, $2a/(a-1)$, respectively. The assumption of piecewise constancy then leads to the two equations

(4.3a,b)\begin{equation} \mu_+=\frac{\mu_-}{3}+\left(\frac{a+1}{2a}\right)\mu_+, \quad \mu_-=\frac{2\mu_-}{3}+\frac{a-1}{2a}\mu_+, \end{equation}

and hence

(4.4a,b)\begin{equation} \mu_+=\frac{2a}{5a-3}, \quad \mu_-=\frac{3(a-1)}{5a-3}. \end{equation}

Therefore,

(4.5)\begin{equation} \langle x\rangle=\mu_-\int_{{-}1}^{0} x\,{{\rm d}\kern0.7pt x}+\mu_+\int_{0}^{1}x\, {{\rm d}\kern0.7pt x}=\frac{1}{2}(\mu_+-\mu_-)=\frac{3-a}{2(5a-3)}. \end{equation}

Now suppose that $a=3+\delta$, where $\delta \ll 1$; then, to leading order, $\mu _+\approx (1-\delta /12)/2$, $\mu _-\approx (1+\delta /12)/2$ and $\langle x\rangle \approx -\delta /24$, confirming the linear response. This result can be reproduced by means of a direct perturbation theory.

It should be noted that the perturbation investigated is not the most general. Were we to change the map so that the two peaks were at different heights, then the invariant measure would become non-differentiable, and the situation would be much like that for the cubic logistic map discussed in § 4.4.

4.3. Maps with linear response (ii): the Lorenz map

In the example above, it was possible to compute the invariant measure $\mu$ explicitly; in general, however, $\mu$ can only be estimated numerically. Here we give an example of a chaotic map, where $\mu$ has to be calculated numerically, and for which the response to small perturbations is indeed linear. The ultra-high resolution computations leading to the results for this map and the following logistic map have been kindly provided by A.M. Rucklidge, and are described in the Appendix.

The Lorenz map, defined by

(4.6)\begin{equation} g(x,\mu_0) = \mu_0 + \mathrm{sgn}(x)\left({-}1 + 1.5 \sqrt{|x|} \right), \end{equation}

is chaotic and depends on the parameter $\mu _0$ in such a way that when $\mu _0=0$, the map is odd and the symmetric chaotic attractor has zero mean. The map, which is shown in figure 8(a), maps the interval $[-1.9,1.9]$ to itself. Having computed the average $\langle {x}\rangle$ for varying $\mu _0$ according to the methods in the Appendix, we can see from figure 9(a) that $\langle {x}\rangle$ is proportional to $\mu _0$ down to $\mu _0=10^{-6}$, below which $\langle {x}\rangle$ is indistinguishable from zero. This linear dependence is in conformity with the results of Bahsoun & Galatolo (Reference Bahsoun and Galatolo2024), who showed that there is a linear response in a tent-like family of maps with a cusp at the turning point.

Figure 8. The maps (a) $g(x,\mu _0)$ (the Lorenz map, defined by (4.6)) and (b) $h(x,\mu _0)$ (the cubic logistic map, defined by (4.7)) for $\mu _0=0$ (red) and $\mu _0>0$ (green).

Figure 9. Average $\langle {x}\rangle$ as a function of $\mu _0$, averaged over 2000 initial conditions and $10^9$ iterates: (a) Lorenz map; (b,c) cubic logistic map for different ranges of $\mu _0$.

4.4. A map with no linear response: the cubic logistic map

The cubic logistic map, which is shown in figure 8(b), is given by

(4.7)\begin{equation} h(x,\mu_0) = \mu_0 + 2.8x - x^3. \end{equation}

Once more, the term in $\mu _0$ breaks the anti-symmetry.

By comparison with the Lorenz map, the cubic logistic map has an entirely different behaviour as $\mu _0$ is varied. For larger values of $\mu _0$, above about $2\times 10^{-6}$, $\langle {x}\rangle$ increases with $\mu _0$ (figure 9b). However, for smaller values of $\mu _0$, as shown in figure 9(c), $\langle {x}\rangle$ does not even depend monotonically on $\mu _0$; it changes sign, and is significantly different from zero for all $\mu _0$ down to $10^{-11}$. Certainly, there is no obvious linear dependence.

What are the properties of these maps that give rise to the different behaviours? Baladi (Reference Baladi2014) has given a general treatment of this question, building on earlier work of Ershov (Reference Ershov1993), who considered the invariant measure on tent maps not of ‘full height’. The difference between the behaviours of the logistic and Lorenz maps can be principally ascribed to the form of the invariant measure $\mu$. For a linear response, $\mu$ has to be sufficiently smooth (e.g. piecewise constant); this can occur when the map is sufficiently mixing. If there are accumulations of discontinuities, for example due to the existence of arbitrarily long-period stable periodic orbits, then there is no smooth response to a small parameter perturbation. The tent map shown is stretching everywhere and has no stable periodic orbits, and $\mu$ is piecewise constant. The Lorenz map is also stretching everywhere, whereas the cubic logistic map, for which the gradient of the map is bounded, has properties similar to the regular logistic map studied by Gottwald et al. (Reference Gottwald, Wormell and Wouters2016), for which there is no linear response.