# Turbulence closure for mixing length theories

Turbulence closure for mixing length theories Abstract We present an approach to turbulence closure based on mixing length theory with three-dimensional fluctuations against a two-dimensional background. This model is intended to be rapidly computable for implementation in stellar evolution software and to capture a wide range of relevant phenomena with just a single free parameter, namely the mixing length. We incorporate magnetic, rotational, baroclinic, and buoyancy effects exactly within the formalism of linear growth theories with non-linear decay. We treat differential rotation effects perturbatively in the corotating frame using a novel controlled approximation, which matches the time evolution of the reference frame to arbitrary order. We then implement this model in an efficient open source code and discuss the resulting turbulent stresses and transport coefficients. We demonstrate that this model exhibits convective, baroclinic, and shear instabilities as well as the magnetorotational instability. It also exhibits non-linear saturation behaviour, and we use this to extract the asymptotic scaling of various transport coefficients in physically interesting limits. convection, stars: evolution, stars: interiors, stars: rotation 1 INTRODUCTION An understanding of turbulent transport and stresses remains one of the major outstanding problems in the astrophysics of fluids. While many pieces of this puzzle are understood in broad strokes, the nature of this problem is such that the details are almost as important as the big picture. The magnetorotational instability (MRI), for instance, is understood conceptually but making predictions that match observed accretion discs is a persistent problem (Murphy & Pessah 2015). Similarly, the solar differential rotation is understood to arise from turbulent stresses but precisely how this works and in balance with what other forces remains uncertain (Schou et al. 1998). Significant progress has indeed been made with three-dimensional turbulence simulations (for examples, see Lee 2013; McKinney et al. 2014; Salvesen et al. 2016) but these are generally relevant only on short time-scales and in small volumes. Performing so-called global simulations over large times and distances requires a turbulence closure model to substitute for resolution at small scales (Launder & Spalding 1974; Canuto 1994). At the other extreme models of stellar evolution generally assume extremely simple analytical transport coefficients to overcome the tremendous gap between turbulent time-scales of minutes and nuclear time-scales of millions of years (Maeder 1995). A variety of such approaches have been developed. For instance, the mixing length theory of Böhm-Vitense (1958) provided a closure of convection. This was then put on firmer theoretical ground by Gough (1977, 2012) and extended to include additional phenomena (Smolec, Houdek & Gough 2011; Lesaffre et al. 2013). Kichatinov (1986) introduced an entirely different closure formalism, arriving at an expression for the so-called Λ-effect (Kichatinov 1987), and later incorporating it under the α–Λ formalism with Rudiger (Kichatinov & Rudiger 1993). What these formalisms have in common is a minimal set of free parameters: the mixing length formalism has just the mixing length, and the formalism of Kichatinov & Rudiger (1993) has just the anisotropy parameter. Another set of models has arisen, which aims to reproduce higher order moments of the turbulent fields. This increases the number of free parameters and a number of approaches have been developed to deal with this. For instance, Garaud, Gagnier & Verhoeven (2017) and Garaud et al. (2010) fit their free parameters against small-scale simulations, while Canuto (1997) fits his against experimental results. In addition, there are models, such as that of Canuto (1994), which fix at least some free parameters by introducing new assumptions, in this case regarding the various relevant time-scales. Regardless of the details of how they close the equations of turbulent moments, models of this sort generally take the form of physically motivated analytic expressions, which provide ready access to scaling laws. Their free parameters then serve to better their agreement with data, at the cost of being less straightforwardly interpreted and extended. The availability of growing computational resources in recent years has provided a new niche in this landscape in the form of computational closure models. These are models which do not seek analytic solutions but which are none the less distinct from attempts to simulate turbulence in all its detail. Some may introduce new dynamical fields, as in the k − ε model (Launder & Spalding 1974), while others invoke effective theories of small-scale motion (Canuto & Hartke 1986). The latter kind are essentially renormalized theories, which accept the cost of having to numerically accommodate complex behaviour in exchange for more precision over a wider variety of phenomena. Combined with perturbation theory, this approach represents a tunable middle-ground between expensive simulations and simple analytic models, allowing the computational cost to be traded off against fidelity to suit the problem at hand. The model we present here is in this spirit. We construct a mixing-length theory, which incorporates three-dimensional fluctuations against a two-dimensional axisymmetric background. This is done by treating each mode as growing with its linear growth rate before saturating at an amplitude set by the turbulent cascade (Lesaffre et al. 2013). Beyond this, the motion in each mode is taken to be uncorrelated. We treat the geometry of the flow in full generality, allowing for baroclinic effects as well as magnetism and rotational shears. To incorporate differential rotation, we use a time-dependent sheared coordinate system (Balbus & Schaan 2012). In this frame, there is a continual flow of modes across Fourier space, lending a time dependence to growth rates. Corrections to saturation amplitudes owing to this flow are incorporated perturbatively with the time derivatives of the growth rate. In Section 2, we describe our closure framework in more detail, paying particular attention to the choice of mixing length. We then develop a perturbative approach for correcting the saturation amplitude in Section 3. In section 4, we introduce the sheared coordinate system and the linearized equations of motion. Finally, in Section 6, we show results from our theory, including calculations for the solar convection zone and accretion discs. The software implementing our model is open source and available under a GPLv3 license. Details of the implementation are given in Appendix C. Tabulated transport coefficients produced by the code are also available under the same license and both may be found at github.com/adamjermyn/Mixer. 2 CLOSURE FORMALISM Turbulent phenomena generically exhibit a cascade of energy between large and small scales (Zhou, McComb & Vahala 1997; Lohse & Xia 2010). With some notable exceptions (Sukoriansky, Dikovskaya & Galperin 2007), this cascade begins at a large scale L0 set by the overall structure of the fluid flow and ends at an extremely small scale Lν related to the microscopic viscosity. Between these scales, yet far from each of them, lies the so-called inertial range where the fluid flow is scale-free (Kolmogorov 1941b). In this range, all correlations of the turbulent motion obey simple power laws. This statement was originally proved by Kolmogorov (1941b) for isotropic turbulence. It was later found to be a broader consequence of the renormalizability of the Navier–Stokes equation (Yakhot & Orszag 1986; Carati 1990) and consequently holds quite generally. This means that there is a single relevant scale L0 for a given turbulent flow, which fully characterises the turbulence as seen by measurements performed over length scales L ≫ L0. This is the modern interpretation and justification of the original mixing length hypothesis, which asserts that turbulent fluctuations on scales L ≪ L0 are not dynamically coupled to the large-scale (L ≫ L0) flow properties (Böhm-Vitense 1958). The scale-free nature of turbulence in the inertial range means that modes of significantly different wavevectors are uncorrelated. A natural extension of this is to assume that all modes of distinct wavevectors are at least approximately uncorrelated. That is, we assume that   \begin{eqnarray} \langle \tilde{\boldsymbol {v}}_{\boldsymbol {k}} \otimes \tilde{\boldsymbol {v}}^*_{\boldsymbol {k}^{\prime }} \rangle &= (2{\pi} )^3 \delta ^3 (\boldsymbol {k}-\boldsymbol {k}^{\prime }) \mathsf {V}_{\boldsymbol {k}}, \end{eqnarray} (1)where v is the velocity, ⊗ denotes the outer product, 〈⋅⋅⋅〉 denotes the time-averaged expectation, $$\tilde{\boldsymbol {v}}_{\boldsymbol {k}}$$ is the amplitude of the Fourier mode with wavevector $$\boldsymbol {k}$$ and $$\mathsf {V}_{\boldsymbol {k}}$$ is the tensor specifying how different components of the same mode are correlated with one another. It is crucial to notice that the quantity $$\mathsf {V}_{\boldsymbol {k}}$$ is also the Reynolds stress of mode $$\boldsymbol {k}$$. This, and several other closely related quantities, are ultimately what we seek. These two-point correlation functions suffice to characterize not only the stresses but also all higher order correlations through Wick’s theorem and perturbation theory (Wick 1950; Isserlis 1918). To determine $$\mathsf {V}_{\boldsymbol {k}}$$, we begin by writing the linearized equations of motion as   \begin{eqnarray} \partial _t \boldsymbol {v}(\boldsymbol {r}) = \mathcal {L}\left[\boldsymbol {v}(\boldsymbol {r}), \partial _i \boldsymbol {v}, \partial _i \partial _j \boldsymbol {v},\ldots , \boldsymbol {r}, t\right], \end{eqnarray} (2)where $$\mathcal {L}$$ is a linear operator of its first argument and $$\boldsymbol {v}$$ is the fluctuating part of the velocity field. In principle, we can work with this operator, though the derivatives of the velocity field make it highly inconvenient. Fortunately, at short length scales, the operator $$\mathcal {L}$$ may be treated as translation-invariant and so we may compute a Fourier transform in $$\boldsymbol {r}$$ without coupling different modes. This gives   \begin{eqnarray} \frac{{\rm d}\tilde{\boldsymbol {v}}_{\boldsymbol {k}}}{{\rm d}t} = \tilde{\mathcal {L}}\left[\tilde{\boldsymbol {v}}_{\boldsymbol {k}}, \boldsymbol {k}, t\right]. \end{eqnarray} (3)The modes are decoupled in this regime so $$\tilde{\mathcal {L}}$$ can be represented by a matrix $$\mathsf {L}$$, and we write   \begin{eqnarray} \frac{{\rm d}\tilde{\boldsymbol {v}}_{\boldsymbol {k}}}{{\rm d}t} = \mathsf {L}(\boldsymbol {k},t) \boldsymbol {v}_{\boldsymbol {k}}. \end{eqnarray} (4) When $$\mathsf {L}$$ is independent of t equation (4) is straightforward to solve and gives us   \begin{eqnarray} \frac{{\rm d}\tilde{\boldsymbol {v}}_{\boldsymbol {k}}}{{\rm d}t} = \sum _i v_{0,i} \hat{v}_{\boldsymbol {k},i} e^{\lambda _i t}, \end{eqnarray} (5)where v0, i are the initial mode amplitudes and $$\hat{v}_{\boldsymbol {k},i}$$ and λi are the normalized right eigenvectors and eigenvalues of $$\mathsf {L}$$, respectively. The vectors $$\hat{v}_{\boldsymbol {k},i}$$ then specify the modes of the system at a given wavevector. If the eigenvalues are not precisely degenerate then modes which begin in phase rapidly become uncorrelated and we may extend equation (1) to the modes at each wavevector and write   \begin{eqnarray} \langle \tilde{\boldsymbol {v}}_{\boldsymbol {k},i} \otimes \tilde{\boldsymbol {v}}^*_{\boldsymbol {k}^{\prime },j} \rangle &= (2{\pi} )^3 \delta ^3 (\boldsymbol {k}-\boldsymbol {k}^{\prime }) \delta _{ij} \mathsf {V}_{\boldsymbol {k},i}. \end{eqnarray} (6)This result holds even when modes are degenerate. Because the time evolution of the Navier–Stokes equation is deterministic, the expectation 〈⋅⋅⋅〉 represents a sum over initial conditions. In this sum all relative phases between the modes are explored, so even degenerate modes become uncorrelated. Inserting equation (5) into equation (6) and summing over j and integrating over $$\boldsymbol {k}$$ gives us   \begin{eqnarray} \mathsf {V}_{\boldsymbol {k},i} &= \hat{v}_{\boldsymbol {k},i} \otimes \hat{v}_{\boldsymbol {k},i} \langle |v_{0,i}|^2 \exp \left[2 t \Re \left[\lambda _i\right]\right] \rangle . \end{eqnarray} (7)Generally, some λi have positive real parts and so in a long-term expectation this exponential diverges. Indeed, it turns out that these growing modes are precisely those which matter! What happens of course is just that these modes eventually reach amplitudes where the linear approximation fails. By assumption, the system is stable over long times relative to the turbulent scale so this must result in these modes saturating. This has been variously described as mode crashing or the action of parasitic modes (Lesaffre, Balbus & Latter 2009; Pessah & Goodman 2009) but, regardless of the mechanism, it simply means that these modes exit the linear regime and find their growth impeded. To complete the closure, we must find the saturation amplitude. Relying again on the scale-free nature of turbulence, we note that this must be a power law in k. That is   \begin{eqnarray} \langle \tilde{v}_{\boldsymbol {k}, i}^2 \rangle = \mathrm{Tr}\left[\mathsf {V}_{\boldsymbol {k}, i}\right] = \frac{A}{M} \left(\frac{k_0}{k}\right)^n, \end{eqnarray} (8)where A depends on the large-scale properties of the flow but is independent of $$\boldsymbol {k}$$, M is the number of modes per wavevector and n is the index of the turbulence. Following Kolmogorov (1941a), we choose n = 11/6 in our model. Appendix A contains a detailed discussion of this choice. The wavenumber k0 is just that of the characteristic scale, and is given by   \begin{eqnarray} k_0 = \frac{2{\pi} }{L_0}. \end{eqnarray} (9)Replacing the divergent expression in equation (7) with this amplitude, we find   \begin{eqnarray} \mathsf {V}_{\boldsymbol {k},i} = \frac{A}{M} \left(\frac{k_0}{k}\right)^n \hat{v}_{\boldsymbol {k},i} \otimes \hat{v}_{\boldsymbol {k},i}. \end{eqnarray} (10) It only remains to determine A. To do this, we note that there is one characteristic length scale L0 and one characteristic time scale, the growth rate ℜ[λi] of the mode. Because A has dimensions of velocity squared, we find   \begin{eqnarray} \mathsf {V}_{\boldsymbol {k},i} = \frac{c}{M} L_0^2 \Re \left[\lambda _i\right]^2 \left(\frac{k_0}{k}\right)^n \hat{v}_{\boldsymbol {k},i} \otimes \hat{v}_{\boldsymbol {k},i}, \end{eqnarray} (11)where c is a dimensionless constant of order unity. This constant, known as the mixing length parameter, varies from theory to theory, so for clarity we set c = 1 in this work but this degree of freedom is important to note when comparing between models. In effect, what we have done is incorporate the non-linearity of turbulence by means of the spectrum while using linear growth rates to set the characteristic scale. In practice, the spectrum acts only to provide a convergent measure over modes (see Appendix A for further discussion), and it is the growth rate and the modes themselves that yield the anisotropies and other phenomena of interest. This is closely related to the approaches of Lesaffre et al. (2013) and Canuto & Hartke (1986). This prescription is easily extended in cases where there are additional dynamical fields, such as the turbulent displacement or a fluctuating magnetic field. The additional fields are simply incorporated into the vector describing the state and M is increased accordingly. We can continue to use equation (8) to fix the amplitude of the entire mode against that of the velocity as long as we know the turbulent index n. Up to this point this prescription is mathematically identical to that of Lesaffre et al. (2013), with the exception that we define the mixing wave vector as in equation (9), while they use π/L0 instead. In the next section, we introduce perturbative corrections to this model to capture a wider variety of phenomena. 3 PERTURBATIVE CORRECTIONS Now consider the case where the matrix $$\mathsf {L}$$ is time-dependent. Most of our reasoning about the behaviour of modes from the previous section still holds but, because the eigenvectors are time-dependent, we no longer have a well-defined notion of a mode as a long-running solution to the equations of motion. When the time dependence is periodic Floquet theory applies (Floquet 1883), but in the cases of interest the time dependence is aperiodic. To recover modes when the time evolution matrix itself evolves and does so aperiodically, we begin by expanding as   \begin{eqnarray} \mathsf {L}(t) = \mathsf {L}(0) + t \frac{{\rm d}\mathsf {L}}{{\rm d}t} + \frac{1}{2}t^2 \frac{{\rm d}^2\mathsf {L}}{{\rm d}t^2} +\cdots \,. \end{eqnarray} (12)This series can be truncated to produce an approximation of $$\mathsf {L}$$, which is accurate in a certain window around t = 0. We may likewise write the velocity at a given wavevector as   \begin{eqnarray} \tilde{\boldsymbol {v}}_{\boldsymbol {k}}(t) = \tilde{\boldsymbol {v}}_{\boldsymbol {k}}(0) + t\left. \frac{{\rm d}\tilde{\boldsymbol {v}}_{\boldsymbol {k}}}{{\rm d}t}\right|_0 + \frac{1}{2}t^2\left. \frac{{\rm d}^2\tilde{\boldsymbol {v}}_{\boldsymbol {k}}(t)}{{\rm d}t^2}\right|_0 +\cdots \,. \end{eqnarray} (13)This suggests defining a new vector   \begin{eqnarray} \Phi _{\boldsymbol {k}}(t) \equiv \left\lbrace \tilde{\boldsymbol {v}}_{\boldsymbol {k}}, \frac{{\rm d}\tilde{\boldsymbol {v}}_{\boldsymbol {k}}}{{\rm d}t}, \frac{{\rm d}^2\tilde{\boldsymbol {v}}_{\boldsymbol {k}}}{{\rm d}t^2},\ldots \right\rbrace , \end{eqnarray} (14)which, in principle, encodes the full time evolution of the velocity field. This vector evolves according to   \begin{eqnarray} \frac{{\rm d}\Phi _{\boldsymbol {k}}}{{\rm d}t} = \mathsf {A} \Phi _{\boldsymbol {k}}, \end{eqnarray} (15)where $$\mathsf {A}$$ is formed of blocks given by   \begin{eqnarray} \mathsf {A}_{ij} = {i \atopwithdelims ()j}\frac{{\rm d}^{i-j}}{{\rm d}t^{i-j}} \mathsf {L}. \end{eqnarray} (16)By definition though we also have   \begin{eqnarray} \frac{{\rm d}\Phi _{\boldsymbol {k},i}}{{\rm d}t} = \Phi _{\boldsymbol {k},i+1}, \end{eqnarray} (17)where $$\Phi _{\boldsymbol {k},0} = \tilde{\boldsymbol {v}}_{\boldsymbol {k}}$$, $$\Phi _{\boldsymbol {k},1} = {\rm d}\tilde{\boldsymbol {v}}_{\boldsymbol {k}}/{\rm d}t$$ and so on. Thus, we are searching for a simultaneous solution of equations (15) and (17). In order to close the system, we must truncate it at some finite order N. Doing so makes the assumption that the behaviour of the system at all greater N is known. Inspired by the solution for time-independent $$\mathsf {L}$$, we try an exponential behaviour. This truncates equation (17) such that it applies only to i < N − 1 and means that we are searching for vectors with   \begin{eqnarray} \left(A \Phi _{\boldsymbol {k}}\right)_{N-1} = \lambda \Phi _{\boldsymbol {k},N-1} \end{eqnarray} (18)and   \begin{eqnarray} \Phi _{\boldsymbol {k},i+1} = (\mathsf {A} \Phi _{\boldsymbol {k}})_i, i < N-1. \end{eqnarray} (19)These equations are most straightforwardly written as a general eigensystem and this has the advantage of restricting the dimension of the linear space to just those states obeying the constraint. This is possible because both $$\mathsf {A}$$ and the constraint are lower triangular in the same basis, and so each row may be substituted into the next, leading to an eigenproblem of the form   \begin{eqnarray} \mathsf {Q} \Phi _{\boldsymbol {k,0}} = \lambda \mathsf {W} \Phi _{\boldsymbol {k,0}}, \end{eqnarray} (20)where $$\mathsf {Q}$$ and $$\mathsf {W}$$ are matrices acting only on the 0-block. For example, in the case where N = 2, our equations are   \begin{eqnarray} \Phi _{\boldsymbol {k},1} &= \mathsf {M} \Phi _{\boldsymbol {k},0} \end{eqnarray} (21)and   \begin{eqnarray} \mathsf {M} \Phi _{\boldsymbol {k},1} + \dot{\mathsf {M}} \Phi _{\boldsymbol {k},0} &= \lambda \Phi _{\boldsymbol {k},1}, \end{eqnarray} (22)which may be put in the form of equation (20) with   \begin{eqnarray} \mathsf {Q} &= \mathsf {M}^2 + \dot{\mathsf {M}} \end{eqnarray} (23)and   \begin{eqnarray} \mathsf {W} &= \mathsf {M}. \end{eqnarray} (24)The eigenvectors of this system are solutions of the original equation (4) because if $$\psi ^i_{\boldsymbol {k}}$$ is such an eigenvector then   \begin{eqnarray} \tilde{\boldsymbol {v}}_{\boldsymbol {k},i}(t) \equiv \sum _{j=0}^{N} \frac{t^j}{j!}\psi ^i_{j} \end{eqnarray} (25)solves   \begin{eqnarray} \frac{{\rm d}\tilde{\boldsymbol {v}}_{\boldsymbol {k},i}}{{\rm d}t} = \mathsf {L}(t) \tilde{\boldsymbol {v}}_{\boldsymbol {k},i}(t) \end{eqnarray} (26)over the time window for which $$\mathsf {L}$$ is well-approximated at N-th order. As a result, we say that ϕi(t) are the instantaneous modes of the system at N-th order and use them and in equation (11). In place of the eigenvalue, we use the instantaneous growth rate of the velocity, which is given by   \begin{eqnarray} g \equiv \frac{1}{2}\frac{{\rm d} v^2}{{\rm d}t} = \frac{\Re \left(\Phi _{\boldsymbol {k},0} \cdot \Phi _{\boldsymbol {k},1}\right)}{|\Phi _{\boldsymbol {k},0}|^2}. \end{eqnarray} (27)This approximation is controlled in the sense that so long as $$\mathsf {L}(t)$$ converges as N grows, so does the inferred velocity history. In this work, we present results with N = 2 so that $$\mathsf {A}$$ involves both $$\mathsf {L}$$ and $$\dot{\mathsf {L}}$$. We leave the exploration of larger N to later work. 4 EQUATIONS OF MOTION We now specialise to the case of an ideal gas obeying the ideal MHD equations. This section largely follows the derivation of Balbus & Schaan (2012) so we present only the pieces necessary to understand later parts of this work as well as the few places where our derivation diverges from theirs. We take the background to be axisymmetric, the fluctuations to be adiabatic and we work in cylindrical coordinates. We neglect both the microscopic viscosity and the microscopic thermal diffusivity because these are both negligible in most circumstances in stellar physics.1 Because our closure model treats turbulent properties as local, we compute all background quantities at a reference point $$\boldsymbol {r}_0$$. Relative to this point, we define the Lagrangian separation $$\boldsymbol {\delta r}$$ and velocity $$\boldsymbol {\delta v}$$ equivalent to $$\boldsymbol {\xi }$$ and $${\rm D}\boldsymbol {\xi }/{\rm D}t$$ of Balbus & Schaan (2012). In addition, we take the Boussinesq approximation that density variations are ignored except in terms involving gravitational acceleration. With the above definitions, the continuity equation may be written as   \begin{eqnarray} \nabla \cdot \boldsymbol {\delta r} = 0. \end{eqnarray} (28) In a fixed coordinate system differential rotation is difficult to analyse so we make two reference frame changes. First, we switch from an inertial frame to one rotating at   \begin{eqnarray} \Omega _0 \equiv \Omega (\boldsymbol {r}_0). \end{eqnarray} (29)Secondly, we make a formal change of coordinates   \begin{eqnarray} \phi \rightarrow \phi - t\delta \boldsymbol {r}\cdot \nabla \Omega , \end{eqnarray} (30)without altering the corresponding unit vectors. Under this last change the gradient transforms as   \begin{eqnarray} \nabla &\rightarrow \nabla - t(\nabla \Omega )\partial _{\phi }. \end{eqnarray} (31)Because the operator $$\mathcal {L}$$ is most easily expressed in Fourier space, we define the transformed wavevector as   \begin{eqnarray} \boldsymbol {q} \equiv \boldsymbol {k} - t k_\phi R \nabla \Omega . \end{eqnarray} (32)With this the transformed MHD and Navier–Stokes equations may be written as   \begin{eqnarray} \boldsymbol {\delta \tilde{B}} = \boldsymbol {B} \cdot \boldsymbol {q} \boldsymbol {\delta \tilde{r}} \end{eqnarray} (33)and   \begin{eqnarray} {\partial _t \delta \tilde{\boldsymbol {v}} + 2\boldsymbol {\Omega }\times \delta \tilde{\boldsymbol {v}}+ \hat{R} R \delta \tilde{\boldsymbol {r}}\cdot \nabla \Omega ^2}\nonumber\\ {\quad- \frac{1}{\gamma \rho }\left(\boldsymbol {\delta \tilde{r}}\cdot \nabla \sigma \right)\nabla \cdot \boldsymbol {\Pi } + \frac{i}{\rho }\boldsymbol {q}\cdot \boldsymbol {\delta \tilde{\Pi }}=0,} \end{eqnarray} (34)where σ is the specific entropy and   \begin{eqnarray} \boldsymbol {\Pi } \equiv p \mathsf {I} - \frac{1}{\mu _0}\left(\boldsymbol {B}\otimes \boldsymbol {B} - \frac{1}{2}B^2 \mathsf {I}\right) \end{eqnarray} (35)is the pressure tensor with $$\mathsf {I}$$ the identity matrix. All quantities prefixed with δ are fluctuating, a tilde denotes the Fourier transformed function, and all other quantities are background fields evaluated at $$\boldsymbol {r}_0$$. It is straightforward to see that this is the same equation as that derived by Balbus & Schaan (2012) once the appropriate relations for the pressure and magnetic force are substituted. The fluctuation in the pressure tensor may be written as   \begin{eqnarray} \boldsymbol {\delta \Pi } = \delta p \mathsf {I} - \frac{1}{\mu _0}\left(\boldsymbol {B}\otimes \boldsymbol {\delta B} + \boldsymbol {\delta B}\otimes \boldsymbol {B} - \mathsf {I} \boldsymbol {B}\cdot \boldsymbol {\delta B} \right), \end{eqnarray} (36)so in Fourier space   \begin{eqnarray} \boldsymbol {\delta \tilde{\Pi }} = \delta \tilde{p} \mathsf {I} - \frac{1}{\mu _0}\left(\boldsymbol {B}\otimes \boldsymbol {\delta \tilde{B}} + \boldsymbol {\delta \tilde{B}}\otimes \boldsymbol {B} - \mathsf {I} \boldsymbol {B}\cdot \boldsymbol {\delta \tilde{B}} \right). \end{eqnarray} (37)Combining this with equation (33) and the Boussinesq approximation (see Appendix B), we find   \begin{eqnarray} \boldsymbol {q}\cdot \boldsymbol {\delta \tilde{\Pi }} = \boldsymbol {q} \delta p - \frac{i}{\mu _0} (\boldsymbol {B}\cdot \boldsymbol {q})^2 \boldsymbol {\delta \tilde{r}}. \end{eqnarray} (38)Note that as did Balbus & Schaan (2012), we take $$\boldsymbol {B}\cdot \boldsymbol {q}$$ to be constant in time as implied by the Boussinesq and ideal-MHD conditions. We now depart from prior work and use this equation along with equation (34) taking the component perpendicular to $$\boldsymbol {q}$$ to eliminate δp and find   \begin{eqnarray} 0 &=&\left(\partial _t \delta \tilde{\boldsymbol {v}} + 2\boldsymbol {\Omega }\times \delta \tilde{\boldsymbol {v}}+ \hat{R} R \delta \tilde{\boldsymbol {r}}\cdot \nabla \Omega ^2\right. \nonumber \\ &&- \left. \frac{1}{\gamma \rho }\left(\boldsymbol {\delta \tilde{r}}\cdot \nabla \sigma \right)\nabla \cdot \boldsymbol {\Pi } + \frac{1}{\mu _0 \rho } (\boldsymbol {B}\cdot \boldsymbol {q})^2 \boldsymbol {\delta \tilde{r}}\right)_{\perp \boldsymbol {q}}, \end{eqnarray} (39)where the notation $$\left(... \right)_{\perp \boldsymbol {q}}$$ denotes the component perpendicular to $$\boldsymbol {q}$$. To construct the matrix version $$\mathsf {L}$$ of these equations, we must choose a coordinate system. Both because of the constraint (28) and because equation (39) is written in the plane perpendicular to $$\boldsymbol {q}$$ we choose the unit vectors   \begin{eqnarray} \hat{\boldsymbol {a}} \equiv \frac{\hat{\boldsymbol {q}}\times \hat{\boldsymbol {w}}}{\sqrt{1-(\hat{\boldsymbol {q}}\,\cdot\, \hat{\boldsymbol {w}})^2}} \end{eqnarray} (40)and   \begin{eqnarray} \hat{\boldsymbol {b}} &\equiv \hat{\boldsymbol {q}}\times \hat{\boldsymbol {a}}, \end{eqnarray} (41)where $$\hat{w}$$ is any unit vector with $$\hat{\boldsymbol {w}} \cdot \hat{\boldsymbol {q}} \ne 1$$. This choice of basis ensures that our vectors are perpendicular to the wavevector. A choice of particular convenience for $$\hat{w}$$ is   \begin{eqnarray} \hat{\boldsymbol {w}} = \frac{\nabla \Omega }{|\nabla \Omega |}. \end{eqnarray} (42)With this choice $$\hat{a}$$ is time-independent, because the component of $$\boldsymbol {q}$$ perpendicular to $$\boldsymbol {w}$$ is time-independent, and so we may write   \begin{eqnarray} \boldsymbol {\delta \tilde{r}} = \alpha \hat{\boldsymbol {a}} + \beta \hat{\boldsymbol {b}} \end{eqnarray} (43)and   \begin{eqnarray} \boldsymbol {\delta \tilde{v}} = \dot{\alpha } \hat{\boldsymbol {a}} + \dot{\beta } \hat{\boldsymbol {b}} + \beta \partial _t \hat{\boldsymbol {b}}. \end{eqnarray} (44)Note that there is a removable singularity when $$\hat{w} \parallel \hat{q}$$. The matrix $$\mathsf {L}$$ is then given by computing the relation between $$\partial _t \lbrace \alpha , \beta , \dot{\alpha }, \dot{\beta }\rbrace$$ and $$\lbrace \alpha , \beta , \dot{\alpha }, \dot{\beta }\rbrace$$. The result is quite unwieldy so we do not present it here but note that it is fully documented in the software in which we implement these equations. 5 STRESSES AND TRANSPORT The equations of motion contain the position and the velocity, so our expanded vector space is   \begin{eqnarray} \Phi = \left\lbrace \delta \boldsymbol {r}, \delta \boldsymbol {v}, \partial _t \delta \boldsymbol {v},\ldots , \right\rbrace . \end{eqnarray} (45)Combining the linearized equations of motion with our closure scheme, we can compute the correlation function   \begin{eqnarray} \langle \Phi \otimes \Phi \rangle = \int \frac{{\rm d}^3 \boldsymbol {k}}{(2{\pi} )^3} \sum _i \langle \Phi ^{i}_{\boldsymbol {k}}\otimes \Phi ^{i*}_{\boldsymbol {k}}\rangle , \end{eqnarray} (46)where the index i ranges over eigenvectors. This function contains all of the usual stresses and transport functions. For instance, the Reynolds stress is   \begin{eqnarray} R \equiv \langle \delta \boldsymbol {v}\otimes \delta \boldsymbol {v}\rangle = \langle \Phi _1 \otimes \Phi _1 \rangle . \end{eqnarray} (47)Likewise up to a dimensionless constant of order unity the turbulent diffusivity is   \begin{eqnarray} d \equiv \langle \delta \boldsymbol {v}\otimes \delta \boldsymbol {r}\rangle = \langle \Phi _1 \otimes \Phi _0 \rangle . \end{eqnarray} (48)and the turbulent viscosity is   \begin{eqnarray} Q \equiv \langle \delta \boldsymbol {v}\otimes \delta \boldsymbol {r}\rangle +\langle \delta \boldsymbol {r}\otimes \delta \boldsymbol {v}\rangle = \langle \Phi _1 \otimes \Phi _0 \rangle +\langle \Phi _0 \otimes \Phi _1 \rangle . \end{eqnarray} (49)Similar expressions hold for the dynamo effect, the transport of magnetic fields, and material diffusion. 6 RESULTS In this section, we exhibit a number of results which come from applying our model to a wide variety of astronomically and physically relevant circumstances. We also compare with the results of Lesaffre et al. (2013) and Kichatinov & Rudiger (1993). We modify the former to use the convention in equation (9) to avoid spurious differences in scale. We likewise assume that our L0 is equal to three times the mixing length of Kichatinov & Rudiger (1993), as this is an inherent freedom in the formalism and resolves an otherwise-persistent scale difference between our model and theirs. These models have been well-tested against a variety of data, most notably helioseismic results, and so provide a useful reference for our work. We have also included more direct comparisons but, because direct experiments are extremely difficult to perform under most circumstances relevant to astrophysics, we have instead included comparisons with simulations and observations where available and applicable. Simulations are often the most useful comparison for stellar phenomena, because a variety of processes, including meridional circulation, can mask the effects of turbulent transport (Kitchatinov 2013). In accretion discs, however, there are several observable quantities, which are thought to correlate closely with the underlying turbulence and these provide very helpful constraints (King, Pringle & Livio 2007). These comparisons and calculations are not intended to be a complete collection of the results our model can produce, nor have we exhaustively explored the circumstances and dependencies of each result. Rather, it is our hope to demonstrate that there is a great deal of interesting physics in this model, that our perturbative corrections give rise to realistic results and reproduce many known results, and that there is much to warrant further exploration along these lines. 6.1 Rotating convection We begin with the effect of rotation on convection in the case of a rotating system with radial pressure and entropy gradients. It is useful to start by comparing our results with those from simulations. Fig. 1 shows the ratios $$\sqrt{\langle \delta v_r^2\rangle /\langle \delta v^2\rangle }$$, $$\sqrt{\langle \delta v_\theta ^2\rangle /\langle \delta v^2\rangle }$$ and $$\sqrt{\langle \delta v_\phi ^2\rangle /\langle \delta v^2\rangle }$$ for several rotation rates as a function of latitude. The positive latitudes come from table 2 of Chan (2001), while the negative are from table 2 of Käpylä, Korpi & Tuominen (2004). In order to match the units for the rotation rates, we put everything in terms of the coriolis number   \begin{eqnarray} \mathrm{Co} \equiv \frac{\Omega h}{\langle \delta v^2\rangle ^{1/2}}, \end{eqnarray} (50)where, following the convention of Käpylä et al. (2004), 〈δv2〉1/2 was computed for a non-rotating system. Figure 1. View largeDownload slide The ratios $$\sqrt{\langle \delta v_r^2\rangle /\langle \delta v^2\rangle }$$ (blue), $$\sqrt{\langle \delta v_\theta ^2\rangle /\langle \delta v^2\rangle }$$ (red), and $$\sqrt{\langle \delta v_\phi ^2\rangle /\langle \delta v^2\rangle }$$ (purple) are shown for our model (solid) and for simulations by (Käpylä et al. 2004, dots, negative latitude) and (Chan 2001, dots, positive latitude) for a wide range of rotation rates as a function of latitude. The rotation rate is captured by the Coriolis number Co = Ωh/〈δv2〉1/2. Our model general overestimates the anisotropy but captures its variation well. Figure 1. View largeDownload slide The ratios $$\sqrt{\langle \delta v_r^2\rangle /\langle \delta v^2\rangle }$$ (blue), $$\sqrt{\langle \delta v_\theta ^2\rangle /\langle \delta v^2\rangle }$$ (red), and $$\sqrt{\langle \delta v_\phi ^2\rangle /\langle \delta v^2\rangle }$$ (purple) are shown for our model (solid) and for simulations by (Käpylä et al. 2004, dots, negative latitude) and (Chan 2001, dots, positive latitude) for a wide range of rotation rates as a function of latitude. The rotation rate is captured by the Coriolis number Co = Ωh/〈δv2〉1/2. Our model general overestimates the anisotropy but captures its variation well. Our model overestimates the anisotropy of the turbulence but captures its symmetries and trends well. For instance, we find that near the poles and in non-rotating systems the θ and ϕ components of the velocity fluctuations have identical magnitudes, in line with the simulations. We reproduce the trend of decreasing anisotropy towards the equator and decreasing anisotropy with increasing rotation, and, in cases where there are differences between the θ and ϕ velocities, we reproduce both their sign and magnitude. In particular, we find that $$\langle \delta v_r^2\rangle \ge \langle \delta v_\theta ^2 \rangle \ge \langle \delta v_\phi ^2\rangle$$, which is seen in these and other simulations (Rüdiger, Egorov & Ziegler 2005a). Likewise, we find that radial motion makes up a greater fraction of the total velocity near the poles than at the equator, and that as the Coriolis number increases $$\langle \delta v_r^2 - \delta v_\theta ^2 - \delta v_\phi ^2 \rangle \rightarrow 0$$, all of which is in agreement with the predictions of Rüdiger et al. (2005b). Our overestimate of the anisotropy may be due to our model incorporating the large-scale fields on all scales, as noted by Lesaffre et al. (2013). This suggests that a future refinement might be to use estimates of the large-scale modes to compute the environment of those at smaller scales, but we do not treat such complications for now. As a further comparison, we consider the off-diagonal Reynolds stresses of both Chan (2001) and Käpylä et al. (2004). These numbers were extracted from table 3 of the former and also table 3 of the latter and are shown along with our predictions in Fig. 2. In the former, they were straightforward to analyse but in the latter they do not provide a precise test because the simulations included a bulk shear. To correct for this, we used a linear expansion to subtract results across simulations, which were identical in all conditions other than the rotation and thereby determine the effect of the rotation alone. As we will see in Section 6.2 this procedure is problematic because the shear may interact non-linearly with the rotation. Furthermore, because these corrections are of the same order as the terms themselves some care must be taken in interpreting the results. Figure 2. View largeDownload slide The ratios $$\sqrt{\langle \delta v_r \delta v_\theta \rangle /\langle \delta v^2\rangle }$$ (red), $$\sqrt{\langle \delta v_\theta \delta v_\phi \rangle /\langle \delta v^2\rangle }$$ (purple), and $$\sqrt{\langle \delta v_r \delta v_\phi \rangle /\langle \delta v^2\rangle }$$ (blue) are shown from our model (solid) and from simulations by (Käpylä et al. 2004, dots, negative latitude) and (Chan 2001, dots, positive latitude) as a function of latitude. Note that Käpylä et al. (2004) cautions that the moderate rotation simulations had difficulty converging, and these results arise as the difference between two simulations, so it is not clear how significant this test is. Our model generally overestimates these stresses, and suggests a different symmetry for the variation (going as sin θ rather than sin (2θ)). Figure 2. View largeDownload slide The ratios $$\sqrt{\langle \delta v_r \delta v_\theta \rangle /\langle \delta v^2\rangle }$$ (red), $$\sqrt{\langle \delta v_\theta \delta v_\phi \rangle /\langle \delta v^2\rangle }$$ (purple), and $$\sqrt{\langle \delta v_r \delta v_\phi \rangle /\langle \delta v^2\rangle }$$ (blue) are shown from our model (solid) and from simulations by (Käpylä et al. 2004, dots, negative latitude) and (Chan 2001, dots, positive latitude) as a function of latitude. Note that Käpylä et al. (2004) cautions that the moderate rotation simulations had difficulty converging, and these results arise as the difference between two simulations, so it is not clear how significant this test is. Our model generally overestimates these stresses, and suggests a different symmetry for the variation (going as sin θ rather than sin (2θ)). Despite these difficulties some trends are clear and sustained between both sets of data. For instance, in the northern hemisphere (θ > 0), 〈δvrδvθ〉 < 0, while in both hemispheres 〈δvrδvϕ〉 < 0, in keeping with predictions and simulations by Rüdiger et al. (2005b). Likewise, we find that 〈vθvϕ〉 > 0 in the northern hemisphere, in agreement with the findings of Rüdiger et al. (2005a). Once more, however, our model overestimates these anisotropic terms by an amount, which is largely invariant as a function of rotation. This suggests that this overestimate is a systematic offset rather than an error in scaling. We also have some difficulty to reproduce the signs of some of the stresses, particularly in the results of Käpylä et al. (2004), though this could simply be a subtraction difficulty. This is supported by the fact that the simulations themselves do not agree on the signs of these terms and highlights the challenges of making comparisons of terms, which are small in magnitude relative to the scale of the turbulence. To better understand which trends are significant and which are artefacts, we have placed data from comparable rotation rates for the two sets of simulations side-by-side in Fig. 3. The top five panels show the same data as in Fig. 1, while the bottom three show the data from Fig. 2. In general, there is good agreement in the top five panels. The data of Käpylä et al. (2004) gives systematically larger anisotropies and the two sets of simulations occasionally differ on the relative magnitudes of the velocity components (i.e. their ordering), but otherwise the two are in good agreement. By contrast, the bottom three panels paint two very divergent pictures. Neither ordering, trends nor signs are consistent between the two sets of simulations. Only the magnitudes agree in these cases. Thus, the two sets of simulations agree that our model systematically overestimates anisotropies and that, beyond that, our model agrees with them to the extent that they agree with one another. Figure 3. View largeDownload slide The functions shown in Figs 1 and 2 are shown from our model (solid), simulations by (Käpylä et al. 2004, dots, negative latitude) and Chan (2001) (crosses, positive latitude) as a function of latitude. The most comparable pairs of rotation rates were placed side-by-side for each function. A solid black line is shown along the equator where the latitude is zero. There is reasonable agreement on the distribution of velocities in direction but not on the correlations between different velocity directions.. Figure 3. View largeDownload slide The functions shown in Figs 1 and 2 are shown from our model (solid), simulations by (Käpylä et al. 2004, dots, negative latitude) and Chan (2001) (crosses, positive latitude) as a function of latitude. The most comparable pairs of rotation rates were placed side-by-side for each function. A solid black line is shown along the equator where the latitude is zero. There is reasonable agreement on the distribution of velocities in direction but not on the correlations between different velocity directions.. Having compared in detail with these simulations, we now consider predictions which go beyond the domain where simulations are possible. In convection with radial gradients, the leading order effect is to transport heat and material radially. Fig. 4 shows 〈δvrδvr〉 and 〈δvrδrr〉, which are the correlation functions controlling this transport. Figure 4. View largeDownload slide The radial velocity correlation function 〈δvrδvr〉 (red) and the radial diffusivity 〈δvrδrr〉 (blue) are shown in linear scale for Ω < |N| (left-hand panel) and log-log scale for Ω > |N| (right-hand panel). These results are for uniform rotation at a latitude of π/4 with no magnetic field. On this and all subsequent figures $$v_r v_r/L_0^2|N|^2$$ should be read as $$\langle \delta v_r \delta v_r \rangle / L_0^2|N|^2$$ and similarly for other correlations. Shown in purple (*, dashed) for comparison is the result of Kichatinov & Rudiger (1993) with an anisotropy factor of 2, which agrees in sign, scale, and variation. The bumps in our results reflect parameter values where the numerical integration was more difficult. All quantities are given in units of the mixing length and the Brünt–Väisälä frequency. Figure 4. View largeDownload slide The radial velocity correlation function 〈δvrδvr〉 (red) and the radial diffusivity 〈δvrδrr〉 (blue) are shown in linear scale for Ω < |N| (left-hand panel) and log-log scale for Ω > |N| (right-hand panel). These results are for uniform rotation at a latitude of π/4 with no magnetic field. On this and all subsequent figures $$v_r v_r/L_0^2|N|^2$$ should be read as $$\langle \delta v_r \delta v_r \rangle / L_0^2|N|^2$$ and similarly for other correlations. Shown in purple (*, dashed) for comparison is the result of Kichatinov & Rudiger (1993) with an anisotropy factor of 2, which agrees in sign, scale, and variation. The bumps in our results reflect parameter values where the numerical integration was more difficult. All quantities are given in units of the mixing length and the Brünt–Väisälä frequency. Both correlators vary at second order in Ω in the slow rotation limit as expected (Kitchatinov 2013; Lesaffre et al. 2013). In the rapid rotation limit on the other hand, they exhibit clear Ω−1 scaling, consistent with what is seen in other closure models and in simulations (Garaud et al. 2010). The quenching of turbulence in this limit arises because the Coriolis effect acts as a restoring force, stabilizing modes. The peak of each correlator is of order unity and occurs when Ω = 0. In fact, for the stress, the maximum is 0.254647 while for the diffusivity it is 0.28125, both of which are consistent to this precision with Lesaffre et al. (2013), noting that we used the definition in equation (9) for their mixing length. This is because our model is precisely the same as theirs in this limit. Based on this and the observed scalings a good approximation is   \begin{eqnarray} \langle \delta v_r \delta r_r \rangle \approx \langle \delta v_r \delta v_r \rangle \approx \frac{1-(\Omega /|N|)^2}{1 - (\Omega /|N|)^3}. \end{eqnarray} (51) Next, we consider the effect of rotation on the r − θ correlation functions. These functions are responsible for latitudinal transport of heat, mass, and momentum and vanish as a result of spherical symmetry in the non-rotating limit. Fig. 5 shows 〈δvrδvθ〉 and 〈δvrδrθ〉 as a function of the rotation rate. Figure 5. View largeDownload slide The absolute value of the r − θ velocity correlation function 〈δvrδvθ〉 (red) and corresponding diffusivity 〈δvrδrr〉 (blue) are shown in log-log scale against rotation rate. These results are for uniform rotation at a latitude of π/4 with no magnetic field. Shown in purple (*, dashed) for comparison is the result of Kichatinov & Rudiger (1993) with an anisotropy factor of 2, which agrees in sign, scale, and variation. All quantities are given in units of the mixing length and the Brünt–Väisälä frequency. Figure 5. View largeDownload slide The absolute value of the r − θ velocity correlation function 〈δvrδvθ〉 (red) and corresponding diffusivity 〈δvrδrr〉 (blue) are shown in log-log scale against rotation rate. These results are for uniform rotation at a latitude of π/4 with no magnetic field. Shown in purple (*, dashed) for comparison is the result of Kichatinov & Rudiger (1993) with an anisotropy factor of 2, which agrees in sign, scale, and variation. All quantities are given in units of the mixing length and the Brünt–Väisälä frequency. In the slow-rotation regime both quantities scale as Ω2, while in the rapid rotation limit they scale as Ω−1. The peak is of order unity and occurs near Ω = |N|. This gives rise to the approximation   \begin{eqnarray} \langle v_r r_\theta \rangle \approx \langle v_r v_\theta \rangle \approx \frac{(\Omega /|N|)^2}{1 + (\Omega /|N|)^3}. \end{eqnarray} (52)These scalings may be interpreted as a competition between symmetry breaking and quenching: the correlation function rises as rotation breaks symmetries but excessive rotation stabilises the system and quenches the turbulent motions. The symmetry is broken quadratically because, at first order, the Coriolis effect only couples radial and azimuthal motions. The properties of turbulence vary with latitude in a rotating system because the rotation axis picks out a preferred direction. Fig. 6 shows the r − r and r − θ stress and diffusivity correlations as a function of latitude. The r − r correlations vary similarly to one another, exhibiting a minimum at the equator and maxima on-axis. On-axis the rotation drops out of the equations and so the on-axis functions are just those for non-rotating convection. The effect of rotation is then largest at the equator, where the convective motion is predominantly perpendicular to the rotation axis. The correlation functions are smallest where the rotation has the largest effect because rotation primarily acts to stabilize modes. Figure 6. View largeDownload slide Various correlation functions are shown as a function of the angle θ from the rotation axis. The functions are the r − r (left-hand panel) and r − θ (right-hand panel) velocity (red) and diffusivity (blue) correlation functions. These results are for uniform rotation at Ω = 0.2|N| (top panel), Ω = |N| (middle panel), and Ω = 5|N| (bottom panel). Shown in purple (*, dashed) for comparison is the KR result, which agrees in sign and variation but not scale. For slow rotation the scale of the variation is generally smaller than we predict, while for fast rotation the variation is somewhat larger. All quantities are given in units of the mixing length and the Brünt–Väisälä frequency. Figure 6. View largeDownload slide Various correlation functions are shown as a function of the angle θ from the rotation axis. The functions are the r − r (left-hand panel) and r − θ (right-hand panel) velocity (red) and diffusivity (blue) correlation functions. These results are for uniform rotation at Ω = 0.2|N| (top panel), Ω = |N| (middle panel), and Ω = 5|N| (bottom panel). Shown in purple (*, dashed) for comparison is the KR result, which agrees in sign and variation but not scale. For slow rotation the scale of the variation is generally smaller than we predict, while for fast rotation the variation is somewhat larger. All quantities are given in units of the mixing length and the Brünt–Väisälä frequency. By contrast, the r − θ correlator is largest in magnitude at mid-latitudes, vanishing both on-axis and at the equator. On-axis this correlation function must vanish because the $$\hat{\theta }$$ unit vector is ill-defined. The sign change between the northern and southern hemispheres occurs because $$(\hat{\boldsymbol {r}} \times \boldsymbol {\Omega })_\phi$$ has the same sign everywhere while $$(\hat{\boldsymbol {\theta }} \times \boldsymbol {\Omega })_\phi$$ changes sign between the hemispheres. This also explains the vanishing correlation at the equator. The quantities of particular interest for studying the origins of differential rotation are the radial-azimuthal correlation functions 〈δvrδvϕ〉 and 〈δvrδrϕ〉. The former provides a stress coupling the angular momentum to radial motions known as the Λ-effect, while the latter provides a viscosity coupling radial shears to azimuthal motion and so acts as a proxy for the α-effect (Kichatinov & Rudiger 1993). Fig. 7 shows these quantities as a function of the rotation rate. In the slow-rotation limit, both scale as Ω before peaking near unity and falling off as Ω−2 in the rapid-rotation limit. The linear scaling at slow rotation rates is a consequence of the Coriolis effect directly coupling radial and azimuthal motions. These quantities fall off more rapidly than the others in the case of rapid rotation because it is preferentially the modes that couple strongly to the Coriolis effect, which are stabilized the most. The absolute scale of our Λ-effect is approximately what is seen in simulations, slightly overestimating relative to Käpylä et al. (2004) and similar to other theoretical predictions (Kitchatinov 2013; Gough 2012). Figure 7. View largeDownload slide The absolute value of the r − ϕ velocity correlation function 〈δvrδvϕ〉 (red) and corresponding diffusivity 〈δvrδrϕ〉 (blue) are shown in log-log scale versus rotation rate. These results are for uniform rotation at a latitude of π/4 with no magnetic field. Shown in purple (*, dashed) for comparison is the result of Kichatinov & Rudiger (1993) with an anisotropy factor of 2, which agrees in sign, variation, and scale up until Ω = |N|, at which point the behaviour differs significantly. Shown in grey (**, dotted) for comparison is 〈δvrδvϕ〉 from that of Lesaffre et al. (2013). This agrees precisely in the Ω → 0 limit and the agreement is good even near Ω ≈ 0.5|N|. All quantities are given in units of the mixing length and the Brünt–Väisälä frequency. Figure 7. View largeDownload slide The absolute value of the r − ϕ velocity correlation function 〈δvrδvϕ〉 (red) and corresponding diffusivity 〈δvrδrϕ〉 (blue) are shown in log-log scale versus rotation rate. These results are for uniform rotation at a latitude of π/4 with no magnetic field. Shown in purple (*, dashed) for comparison is the result of Kichatinov & Rudiger (1993) with an anisotropy factor of 2, which agrees in sign, variation, and scale up until Ω = |N|, at which point the behaviour differs significantly. Shown in grey (**, dotted) for comparison is 〈δvrδvϕ〉 from that of Lesaffre et al. (2013). This agrees precisely in the Ω → 0 limit and the agreement is good even near Ω ≈ 0.5|N|. All quantities are given in units of the mixing length and the Brünt–Väisälä frequency. 6.2 Differential rotation and convection We now turn to the dependence of convective transport coefficients on differential rotation. We expand our closure model to linear order in the shear and so restrict this analysis to cases where the dimensionless shear |R∇ln Ω| is at most of order unity. Fig. 8 shows the r − θ and r − ϕ velocity and diffusivity correlation functions as a function of differential rotation for a situation where ∇Ω is at an angle of π/4 relative to the pressure gradient. All four functions behave linearly near the origin, with intercept set by the stress and diffusivity in the uniform rotation limit. This is precisely as expected: the intercept is non-zero, giving rise to the Λ-effect, while the slope is non-zero, giving rise to the α-effect (Kichatinov & Rudiger 1993). Note that the favourable comparison of our results with those of Kitchatinov (2013) is helpful because their model was implemented in a two-dimensional solar model, which compared well with helioseismic observations. Figure 8. View largeDownload slide The r − θ (left-hand panel) and r − ϕ (right-hand panel) velocity (red) and diffusivity (blue) correlation functions are shown in log-scale versus the differential rotation. These results are for a convecting region with differential rotation in the cylindrical radial direction, Ω = 0.1|N| and no magnetic field at a latitude of π/4. Shown in purple (*, dashed) for comparison is the result of Kichatinov & Rudiger (1993) with an anisotropy factor of 2. This disagrees in sign and on the magnitude of the slope but agrees in the sign of the slope. Shown in grey (**, dotted) for comparison is 〈δvrδvϕ〉 of Lesaffre et al. (2013). This generally predicts smaller stresses though with the same sign and slope sign as our model. All quantities are given in units of the mixing length and the Brünt–Väisälä frequency. Figure 8. View largeDownload slide The r − θ (left-hand panel) and r − ϕ (right-hand panel) velocity (red) and diffusivity (blue) correlation functions are shown in log-scale versus the differential rotation. These results are for a convecting region with differential rotation in the cylindrical radial direction, Ω = 0.1|N| and no magnetic field at a latitude of π/4. Shown in purple (*, dashed) for comparison is the result of Kichatinov & Rudiger (1993) with an anisotropy factor of 2. This disagrees in sign and on the magnitude of the slope but agrees in the sign of the slope. Shown in grey (**, dotted) for comparison is 〈δvrδvϕ〉 of Lesaffre et al. (2013). This generally predicts smaller stresses though with the same sign and slope sign as our model. All quantities are given in units of the mixing length and the Brünt–Väisälä frequency. A key difference between our work and what we compare with in Fig. 8 is that, while we predict the same sign and comparable magnitude for the α-effect in the zero-shear limit, the effect changes sign near |R∇ln Ω| ≈ 0.5, indicating that, at least for this configuration, this is the point at which non-linear effects become important. This does not represent a particularly severe shear and highlights a key point that the correlation functions we find are generally non-linear in all of the small parameters in which one might wish to expand. Our model captures this non-linear behaviour despite being carried out to linear order in |R∇ln Ω|. This is because, in our expansion, the time evolution operator is what is expanded linearly. The resulting eigenvalues and eigenvectors are generally non-linear functions of this operator. This caution aside, there is a significant regime where the α–Λ expansion is valid and, in this regime, key quantities of interest are the derivatives of the various correlation functions with respect to the shear |R∇Ω|. Fig. 9 shows these derivatives as a function of Ω. The r − ϕ stress derivative is constant in Ω. This means that the stress scales as R∇Ω. This is as expected (see, e.g. equation 79 of Lesaffre et al. 2013) and indicates that there is a well-defined effective viscosity transporting angular momentum. This viscosity is given by   \begin{eqnarray} \nu _{r\phi } \approx L_0^2 |N|. \end{eqnarray} (53) Figure 9. View largeDownload slide The derivatives of various correlation functions with respect to |R∇Ω| are shown as a function of Ω, with both axes log-scaled. The functions are the r − θ (left-hand panel) and r − ϕ (right-hand panel) velocity (red) and diffusivity (blue) correlation functions. These results are for a convecting region with differential rotation in the cylindrical radial direction and no magnetic field at a latitude of π/4. All quantities are given in units of the mixing length and the Brünt–Väisälä frequency. Figure 9. View largeDownload slide The derivatives of various correlation functions with respect to |R∇Ω| are shown as a function of Ω, with both axes log-scaled. The functions are the r − θ (left-hand panel) and r − ϕ (right-hand panel) velocity (red) and diffusivity (blue) correlation functions. These results are for a convecting region with differential rotation in the cylindrical radial direction and no magnetic field at a latitude of π/4. All quantities are given in units of the mixing length and the Brünt–Väisälä frequency. By contrast, the derivatives of the r − θ correlations as well as the r − ϕ diffusivity all diverge in the limit as Ω → 0. In particular, the r − θ correlations diverge as Ω−1, while the r − ϕ diffusivity diverges as Ω−2. These divergences are signatures of symmetry breaking. They indicate that the direction in which the R∇Ω → 0 limit is approached matters. That is, this limit can be approached by first letting Ω → 0 and then differentiating or by differentiating and then taking Ω → 0 and the divergence we find in the latter approach indicates that the order matters. When Ω = 0 and |R∇Ω| = 0, there is a symmetry between ±θ and between ±ϕ. As a result both the r − θ and r − ϕ terms vanish in this limit. When Ω ≠ 0 these symmetries are broken by the rotation and we know from Figs 5 and 7 that this occurs at first order for r − ϕ and second order for r − θ. In the opposing limit, the situation is different because in the time evolution described by equation (39) $$\mathsf {L}$$ is independent of |R∇Ω| when Ω = 0. There is, however, a dependence on |R∇Ω| through the time-dependence of $$\boldsymbol {q}$$. This breaks the ϕ symmetry because $$\partial _t \boldsymbol {q}$$ is proportional to qϕR∇Ω and hence is sensitive to ϕ. It does not, however, break the θ symmetry, because qϕR∇Ω is symmetric with respect to changing the signs of both θ and $$\boldsymbol {q}$$. It follows then that we should find divergences in the r − θ correlation derivatives owing to the path-dependence of the zero-rotation limit and that we should find the r − ϕ derivatives to be generally well-behaved. The curious divergence is then that in the r − ϕ diffusivity, because this correlation function does not suffer from a symmetry-derived path-dependence. This arises because the differential rotation means that $$\mathsf {L}$$ is time-dependent. This introduces polynomial corrections to the usual exponential growth, as discussed in Section 3. This formalism captures the fact that the differential rotation turns vertical displacement into ϕ displacements, which vary as polynomials in time. There are, therefore, modes with very small radial velocities, which nevertheless have large azimuthal displacements and these dominate the diffusivity derivative. These modes grow proportional to |R∇Ω| and their growth may proceed in the azimuthal direction until bounded by the Coriolis effect at a time Ω−1. As a result these modes contribute to the diffusivity as |R∇ln Ω| and hence lead to a diverging derivative in |R∇Ω| as Ω → 0. 6.3 Differential rotation and stable stratification Stably stratified regions are those with   \begin{eqnarray} N^2 > 0, \end{eqnarray} (54)such that buoyancy acts to counter perturbations in the vertical direction. This tends to damp turbulence. In the presence of such damping, there can still be turbulence if there is also a shear. The classic example of this is the Kelvin–Helmholtz phenomenon, which can occur in such a system if the Richardson criterion   \begin{eqnarray} \frac{|{\rm d}u/{\rm d}z|^2}{|N|^2} > \frac{1}{4} \end{eqnarray} (55)is satisfied (Zahn 1993). Here, u is the velocity and z is the coordinate parallel to the stratification. Even when this criterion is not satisfied, latitudinal shear can still generate turbulence (Canuto et al. 2008). These motions are suppressed in vertical extent by the stratification and hence are primarily confined to the plane perpendicular to the stratification direction. Fig. 10 shows the dependence on shear strength of all six stress components in a rotating stably stratified zone with latitudinal rotational shear. All six exhibit linear scaling with the shear strength. This is unusual in an otherwise-stable zone because it implies a viscosity which, to leading order, does not depend on the shear. That is,   \begin{eqnarray} \nu _{ij} \approx L_0^2 N f_{ij}\left(\frac{\Omega }{|N|}\right), \end{eqnarray} (56)where fij is some function of the angular velocity. Fig. 11 shows the dependence of the stress components on Ω/|N| for fixed |R∇ln Ω| = 0.1. The r − θ and θ − ϕ stresses vary as Ω3 in the slow-rotation regime and as Ω2 for rapid rotation. The other components all scale as Ω2 in both regimes. Thus, for instance, frϕ = Ω/|N| because the viscosity is the derivative of the stress with respect to the shear, and hence   \begin{eqnarray} \nu _{r\phi } \approx 10^{-5} L_0^2 \Omega . \end{eqnarray} (57) Figure 10. View largeDownload slide The absolute value of the r − r (red), r − θ (blue) and r − ϕ (purple) velocity correlation functions are shown as a function of |R∇ln Ω|, with both axes log-scaled. These results are for a stably stratified region with differential rotation in the radial direction, Ω = 0.1|N| and no magnetic field. The data is computed for a point on the equator with radial differential rotation. All quantities are given in units of the mixing length and the Brünt–Väisälä frequency. Figure 10. View largeDownload slide The absolute value of the r − r (red), r − θ (blue) and r − ϕ (purple) velocity correlation functions are shown as a function of |R∇ln Ω|, with both axes log-scaled. These results are for a stably stratified region with differential rotation in the radial direction, Ω = 0.1|N| and no magnetic field. The data is computed for a point on the equator with radial differential rotation. All quantities are given in units of the mixing length and the Brünt–Väisälä frequency. Figure 11. View largeDownload slide The absolute value of the r − r (red), r − θ (blue), and r − ϕ (purple) velocity correlation functions are shown as a function of Ω/|N| for fixed |R∇ln Ω| = 0.1, with both axes log-scaled. These results are for a stably stratified region with differential rotation in the radial direction and no magnetic field. The data is computed for a point on the equator with radial differential rotation. All quantities are given in units of the mixing length and the Brünt–Väisälä frequency. Figure 11. View largeDownload slide The absolute value of the r − r (red), r − θ (blue), and r − ϕ (purple) velocity correlation functions are shown as a function of Ω/|N| for fixed |R∇ln Ω| = 0.1, with both axes log-scaled. These results are for a stably stratified region with differential rotation in the radial direction and no magnetic field. The data is computed for a point on the equator with radial differential rotation. All quantities are given in units of the mixing length and the Brünt–Väisälä frequency. The scaling in equation (57) arises owing to the centrifugal term, which has a destabilising effect when Ω increases with $$\hat{R}$$. When |R∇Ω| = 0 this effect is not present so the system is stable but introducing a small differential rotation produces an acceleration proportional to $$\Omega R \boldsymbol {\delta r}\cdot \nabla \Omega$$ and hence   \begin{eqnarray} \partial _t^2 \boldsymbol {\delta r} \approx g^2 \boldsymbol {\delta r} \propto \hat{R}\Omega R \boldsymbol {\delta r}\cdot \nabla \Omega , \end{eqnarray} (58)which means that the stress scales as Ω∇Ω and thence the viscosity scales as Ω. Because the magnitudes of the stresses are always ordered in the same way, the same terms are always the most significant. From largest to smallest, the stresses are r − r, r − ϕ, ϕ − ϕ, θ − θ, θ − ϕ, and r − θ. This group is nearly separated into diagonal stresses, which are larger, and off-diagonal stresses, which are smaller. The exception to this rule is the r − ϕ stress, which is special because it is the term which directly couples to the shear. The ordering of the remaining terms is not surprising because the off-diagonal stresses are typically mediated by a coupling between different directions, whereas the on-diagonal stresses require no such coupling. To better understand the effect of our perturbative corrections, we computed the same results without them. This produced stresses, which were zero to within numerical precision in all cases, indicating that the entire contribution in this case is coming from the perturbation. However with a different angle of differential rotation, we obtained non-zero results. It is instructive then to compare Fig. 12 with Fig. 13. These show the same correlation functions as each other in the same physical scenario, with differential rotation this time at an angle of π/4, but the former uses the first order perturbative expansion while the latter only expands to zeroth order. The difference between the two calculations is striking: many of the correlation functions have fundamentally different scalings when the perturbative corrections are taken into account. In particular, the r − θ and r − r stresses are both quadratic in the shear and the r − ϕ and θ − ϕ stresses both vary as the shear to the 3/2 power, whereas they are all linear in the shear in the expanded calculation. This difference relates in part to the centrifugal term, which couples the displacement to the acceleration. Without expanding the equations of motion, we would have δr ∝ δv, because the mode would need to be an eigenvector of $$\mathsf {M}$$. The modes which couple to the centrifugal term would still grow according to equation (58) but, for most modes, arranging for the displacement to couple to this term requires coupling to the stabilizing buoyant term too. To make this clearer, in Fig. 14, we have computed the growth rate as a function of wave-vector orientation without using the perturbative expansion. There are several rapidly growing regions, oriented at angles of ± π/4 relative to the vertical. These angles represent a compromise between maximizing the magnitude of the centrifugal acceleration and maximizing its projection on to the velocity, both subject to the Boussinesq condition that motion be in the plane perpendicular to $$\boldsymbol {q}$$. Figure 12. View largeDownload slide The absolute value of the r − r (red), r − θ (blue), and r − ϕ (purple) velocity correlation functions are shown as a function of |R∇ln Ω|, with both axes log-scaled. The correlation functions are evaluated at first order in the perturbative expansion rather than first order. These results are for a stably stratified region with differential rotation in the radial direction, Ω = 0.1|N| and no magnetic field. The data are computed for a point on the equator with differential rotation at an angle of π/4. All quantities are given in units of the mixing length and the Brünt–Väisälä frequency. Figure 12. View largeDownload slide The absolute value of the r − r (red), r − θ (blue), and r − ϕ (purple) velocity correlation functions are shown as a function of |R∇ln Ω|, with both axes log-scaled. The correlation functions are evaluated at first order in the perturbative expansion rather than first order. These results are for a stably stratified region with differential rotation in the radial direction, Ω = 0.1|N| and no magnetic field. The data are computed for a point on the equator with differential rotation at an angle of π/4. All quantities are given in units of the mixing length and the Brünt–Väisälä frequency. Figure 13. View largeDownload slide The absolute value of the r − r (red), r − θ (blue) and r − ϕ (purple) velocity correlation functions are shown as a function of |R∇ln Ω|, with both axes log-scaled. The correlation functions are evaluated at zeroth order in the perturbative expansion rather than first order. These results are for a stably stratified region with differential rotation in the radial direction, Ω = 0.1|N| and no magnetic field. The data are computed for a point on the equator with differential rotation at an angle of π/4. All quantities are given in units of the mixing length and the Brünt–Väisälä frequency. Figure 13. View largeDownload slide The absolute value of the r − r (red), r − θ (blue) and r − ϕ (purple) velocity correlation functions are shown as a function of |R∇ln Ω|, with both axes log-scaled. The correlation functions are evaluated at zeroth order in the perturbative expansion rather than first order. These results are for a stably stratified region with differential rotation in the radial direction, Ω = 0.1|N| and no magnetic field. The data are computed for a point on the equator with differential rotation at an angle of π/4. All quantities are given in units of the mixing length and the Brünt–Väisälä frequency. Figure 14. View largeDownload slide The square of the growth rate is shown as a function of wave-vector orientation on a logarithmic colour scale. The wave-vector is specified by a magnitude and two angles, θ(q) and ϕ(q), which are spherical angles relative to the $$\hat{z}$$ direction. These rates were computed with a zeroth-order expansion. Regions with squared growth rates below 10−16 are shown in white. All quantities are given in units of the mixing length and the Brünt–Väisälä frequency. Figure 14. View largeDownload slide The square of the growth rate is shown as a function of wave-vector orientation on a logarithmic colour scale. The wave-vector is specified by a magnitude and two angles, θ(q) and ϕ(q), which are spherical angles relative to the $$\hat{z}$$ direction. These rates were computed with a zeroth-order expansion. Regions with squared growth rates below 10−16 are shown in white. All quantities are given in units of the mixing length and the Brünt–Väisälä frequency. By contrast, the growth rates in the expanded system, shown in Fig. 15, are significant over a much wider swath of parameter space. This is because, in the expanded system, the displacement and velocity need not be parallel so the displacement can be chosen to maximize the centrifugal term while the velocity can be chosen to maximize the projection of the acceleration on to the velocity. Figure 15. View largeDownload slide The square of the growth rate is shown as a function of wave-vector orientation on a logarithmic colour scale. The wave-vector is specified by a magnitude and two angles, θ(q) and ϕ(q), which are spherical angles relative to the $$\hat{z}$$ direction. These rates were computed with a first-order expansion. Regions with squared growth rates below 10−16 are shown in white. All quantities are given in units of the mixing length and the Brünt–Väisälä frequency. Figure 15. View largeDownload slide The square of the growth rate is shown as a function of wave-vector orientation on a logarithmic colour scale. The wave-vector is specified by a magnitude and two angles, θ(q) and ϕ(q), which are spherical angles relative to the $$\hat{z}$$ direction. These rates were computed with a first-order expansion. Regions with squared growth rates below 10−16 are shown in white. All quantities are given in units of the mixing length and the Brünt–Väisälä frequency. 6.4 Baroclinic instability The baroclinic instability arises in otherwise stably stratified fluids when the entropy gradient is not parallel to the pressure gradient (Killworth 1980). In fact, this is part of a family of instabilities, which includes the convective instability (Lebovitz 1965). This family provides a continuous connection between the unstable convective and stably stratified limits. To explore, it consider Fig. 16 which shows the variation of r − r and r − θ correlation functions against the angle δ between the entropy gradient and the pressure gradient. The radial correlations peak when the two gradients are aligned. This is the convective limit. These correlations fall to zero in the opposing limit where the two gradients are anti-aligned, which is the stably stratified limit. In between these limits, the behaviour is approximately that of cos 2δ. Figure 16. View largeDownload slide Various correlation functions are shown as a function of the angle δ between the entropy and pressure gradients. The functions are the r − r (left-hand panel) and r − θ (right-hand panel) velocity (red) and diffusivity (blue) correlation functions. These results are for a non-rotating convective region on the equator and with no magnetic field. Shown in purple (*, dashed) for comparison is the KR result. This agrees in sign, and for r − r agrees in scale, but their r − θ prediction is considerably larger. Notably, this comparison is precisely as cos (δ) (left-hand panel) and sin (δ) (right-hand panel) and crosses zero at non-extremal angles. This is most likely because their theory is not designed for nearly stable regions with extreme baroclinicity. All quantities are given in units of the mixing length and the Brünt–Väisälä frequency. Figure 16. View largeDownload slide Various correlation functions are shown as a function of the angle δ between the entropy and pressure gradients. The functions are the r − r (left-hand panel) and r − θ (right-hand panel) velocity (red) and diffusivity (blue) correlation functions. These results are for a non-rotating convective region on the equator and with no magnetic field. Shown in purple (*, dashed) for comparison is the KR result. This agrees in sign, and for r − r agrees in scale, but their r − θ prediction is considerably larger. Notably, this comparison is precisely as cos (δ) (left-hand panel) and sin (δ) (right-hand panel) and crosses zero at non-extremal angles. This is most likely because their theory is not designed for nearly stable regions with extreme baroclinicity. All quantities are given in units of the mixing length and the Brünt–Väisälä frequency. By contrast, the r − θ correlations behave approximately as sin δ, and vanishes when δ = 0. This is because both the aligned and the anti-aligned limits are spherically symmetric and so must have this correlation function vanish. Deviations from the convective limit give rise to linear scaling so the convective baroclinic instability transports heat and momentum at first order in the baroclinicity. This is an entirely distinct phenomenon from the thermal wind balance, which is a large-scale effect while this results from integrating out the small-scale turbulent modes. In the stable limit perturbations arise quadratically, a deviation from the behaviour of sin δ. This is because there are no existing turbulent motions to perturb, and so each position and velocity component is linear in δ and gives rise to a quadratic two-point correlation function. 6.5 Stellar magnetism We now turn to the impact of the magnetic field on convective turbulence in stars. Fig. 17 shows 〈δvrδvr〉 in a mildly rotating (Ω = 0.1|N|) convection zone as a function of B for three polarisations; radial ($$\boldsymbol {B} \parallel \hat{\boldsymbol {r}}$$), latitudinal ($$\boldsymbol {B} \parallel \hat{\boldsymbol {\theta }}$$), and longitudinal ($$\boldsymbol {B} \parallel \hat{\boldsymbol {\phi }}$$). As the field increases, the stress falls off. This is because the field quenches the turbulence by providing a stabilizing restoring force, and is in general agreement with Canuto & Hartke (1986). Interestingly, the only significant differences are between the radial and angular field polarisations! The θ and ϕ polarisations show precisely the same behaviour out to very strong fields. This is a result of symmetry, because the radial stress is not sensitive to rotation about the radial direction. The deviation seen with strong fields is a numeric artefact and decreases with increasing integration time. Figure 17. View largeDownload slide The stress 〈δvrδvr〉 is shown as a function of magnetic field strength. The magnetic field is polarized radially (red), longitudinally (purple), and latitudinally (blue). The system is rigidly rotating at Ω = 0.1|N| at a latitude of π/4. All quantities are given in units of the mixing length and Brünt–Väisälä frequency. Figure 17. View largeDownload slide The stress 〈δvrδvr〉 is shown as a function of magnetic field strength. The magnetic field is polarized radially (red), longitudinally (purple), and latitudinally (blue). The system is rigidly rotating at Ω = 0.1|N| at a latitude of π/4. All quantities are given in units of the mixing length and Brünt–Väisälä frequency. By contrast, consider 〈δvrδvϕ〉, shown in Fig. 18. This component, along with the corresponding Maxwell stress, is responsible for transporting angular momentum. Interestingly, it shows differences amongst all polarisations, with the strongest difference between the θ polarization and the others. This is because the stress is mixed between different directions and so is sensitive to all variations in the magnetic field direction. The large difference of the θ polarization relative to the others reflects the fact that motion is damped perpendicular to the magnetic field so the θ polarization damps motion in both directions involved in this component of the stress, whereas the r and ϕ polarizations only damp motion in one of the two directions. Figure 18. View largeDownload slide The stress 〈δvrδvϕ〉 is shown as a function of magnetic field strength. The magnetic field is polarized radially (red), longitudinally (purple) and latitudinally (blue). The system is rigidly rotating at Ω = 0.1|N| at a latitude of π/4. All quantities are given in units of the mixing length and Brünt–Väisälä frequency. Figure 18. View largeDownload slide The stress 〈δvrδvϕ〉 is shown as a function of magnetic field strength. The magnetic field is polarized radially (red), longitudinally (purple) and latitudinally (blue). The system is rigidly rotating at Ω = 0.1|N| at a latitude of π/4. All quantities are given in units of the mixing length and Brünt–Väisälä frequency. 6.6 Magnetorotational instability As a final example, we consider the MRI (Chandrasekhar 1960). This instability arises in magnetized fluids undergoing Keplerian orbital motion. Fig. 19 shows the r − r, θ − ϕ, r − ϕ, and r − θ Reynolds and Maxwell stresses for an accretion disc with a vertical magnetic field. Contrary to predictions (Chandrasekhar 1960) none of the Reynolds stresses vanish in the zero-field limit. This is because the linear system supports short-term growing modes but, while they only grow in the short-time limit, our numerical methods are not sensitive to that effect at this order. In principle, at higher order, this phenomenon should become evident and so this may be interpreted as an artefact associated with our expanding to low order in |R∇Ω| > 1. Despite this, it is likely that other non-magnetic processes can destabilise these modes even in the long term and so we feel it is appropriate to at least consider them (cf. Luschgy & Pagès 2006). The Maxwell stresses by contrast do vanish as B → 0. This is to be expected because they are proportional to B2. Figure 19. View largeDownload slide From the top to bottom panel are the r − ϕ, r − r, r − θ, and θ − ϕ stresses. The Reynolds (velocity) stresses are in red and the Maxwell (magnetic) stresses are in blue. Note that it is the negative r − r Maxwell stress, which is shown to make the comparison with the Reynolds stresses clearer. In all cases, they are shown as a function of B for a Keplerian disc. The magnetic field is taken parallel to $$\hat{z}$$. The system is taken to be stably stratified in the vertical direction with |N| = Ω and hence L0 = h = R. All quantities are given in units of the mixing length and Ω. Figure 19. View largeDownload slide From the top to bottom panel are the r − ϕ, r − r, r − θ, and θ − ϕ stresses. The Reynolds (velocity) stresses are in red and the Maxwell (magnetic) stresses are in blue. Note that it is the negative r − r Maxwell stress, which is shown to make the comparison with the Reynolds stresses clearer. In all cases, they are shown as a function of B for a Keplerian disc. The magnetic field is taken parallel to $$\hat{z}$$. The system is taken to be stably stratified in the vertical direction with |N| = Ω and hence L0 = h = R. All quantities are given in units of the mixing length and Ω. As the magnetic field increases, the r − ϕ and θ − ϕ Reynolds stresses change sign. This indicates the onset of MRI modes, which have the opposite sign to the zero-field correlations. This effect saturates when vA ≈ Ωh, where h is the scale height of the disc. The total r − ϕ stress saturates at roughly 10−2(hΩ)2, which lies between those typically found in simulations and those inferred from observations (Starling et al. 2004; King et al. 2007). Note that at the saturation point the Maxwell and Reynolds stresses are comparable, and beyond this point the Maxwell stress increases while the Reynolds stress falls off. Above the saturation point the Reynolds stresses drop off as the magnetic field quenches the turbulence. This is precisely what is expected for the MRI (Balbus & Hawley 1991). The Maxwell stresses, however, continue to grow, again in line with expectations. Some care is required to interpret these results because they were computed for a fixed field and that field may or may not be stable under the action of the turbulence it generates (Pessah, Chan & Psaltis 2006). Furthermore, there are challenges with the α-disc prescription, which make the specific stress components more difficult to interpret (Pessah, Chan & Psaltis 2008). Nevertheless, it is encouraging that what we see matches well with both observations and simulations. 7 CONCLUSION We have derived a turbulent closure model, which incorporates shear, rotation, and magnetism as well as a full three-dimensional spectrum of fluctuations. We have also presented a new perturbative approach to incorporate time-dependence in the evolution equations. This model, which is implemented in an open source numerical software package, fully reproduces many known phenomena such as the MRI, baroclinic instability, rotational quenching, and more classic shear instabilities. Using this model, we have determined the asymptotic behaviour of a wide variety of correlation functions and transport coefficients under a wide range of circumstances, many of which do not appear in the literature. We have further explored the behaviour of turbulent transport coefficients in intermediate regimes where no single phenomenon dominates, such as in the critical MRI. In these cases, the behaviour is generally complex and does not separate easily into components associated with the different pieces of input physics. The closure formalism developed here fills a new niche in the landscape of solutions to turbulent transport, covering enough phenomena to be useful to understand those operating in stars, planets, and accretion discs, while being rapid enough to be incorporated into stellar evolution codes on nuclear time-scales. In the future, we hope to provide further refinements and comparisons with direct numerical simulations as well as experiments. In addition, it would be interesting to explore the results of this model to higher order in the shear and, even at this order, there are many results which deserve more analysis than we have given here. Acknowledgements ASJ acknowledges financial support from a Marshall Scholarship as well as support from the Institute of Astronomy, École Normale Supérieure (ENS), and Centre for Excellence in Basic Sciences (CEBS) to work at ENS Paris and CEBS in Mumbai. PL acknowledges travel support from the french PNPS (Programme National de Physique Stellaire) and from CEBS. CAT thanks Churchill College for his fellowship. SMC is grateful to the IOA for support and hopsitality and thanks the Cambridge-Hamied exchange program for financial support. The authors also thank Rob Izzard and Science and Technology Facilities Council Grant ST/L003910/1 for CPU cycles, which aided in this work. Footnotes 1 It would not be difficult, however, to incorporate them into this framework at a later date. REFERENCES Ashkenazi S., Steinberg V., 1999, Phys. Rev. Lett. , 83, 4760 CrossRef Search ADS   Balbus S. A., Hawley J. F., 1991, ApJ , 376, 214 CrossRef Search ADS   Balbus S. A., Schaan E., 2012, MNRAS , 426, 1546 CrossRef Search ADS   Benzi R., Tripiccione R., Massaioli F., Succi S., Ciliberto S., 1994, Europhys. Lett. , 25, 341 CrossRef Search ADS   Berntsen J., Espelid T. O., Genz A., 1991, ACM Trans. Math. Softw. , 17, 437 CrossRef Search ADS   Böhm-Vitense E., 1958, ZAp , 46, 108 Canuto V. M., 1994, ApJ , 428, 729 CrossRef Search ADS   Canuto V. M., 1997, ApJ , 482, 827 CrossRef Search ADS   Canuto V. M., Hartke G. J., 1986, A&A , 168, 89 Canuto V. M., Cheng Y., Howard A. M., Esau I. N., 2008, J. Atmos. Sci. , 65, 2437 CrossRef Search ADS   Carati D., 1990, Phys. Rev. A , 41, 3129 CrossRef Search ADS PubMed  Chan K. L., 2001, ApJ , 548, 1102 CrossRef Search ADS   Chandrasekhar S., 1960, Proc. Natl. Acad. Sci. , 46, 253 CrossRef Search ADS   Dobrowolny M., Mangeney A., Veltri P., 1980, Phys. Rev. Lett. , 45, 144 CrossRef Search ADS   Floquet G., 1883, Annales scientifiques de l’École normale supérieure , 12, 47 CrossRef Search ADS   Garaud P., Ogilvie G. I., Miller N., Stellmach S., 2010, MNRAS , 407, 2451 CrossRef Search ADS   Garaud P., Gagnier D., Verhoeven J., 2017, ApJ , 837, 133 CrossRef Search ADS   Genz A., Malik A., 1980, J. Comput. Appl. Math. , 6, 295 CrossRef Search ADS   Goldreich P., Sridhar S., 1995, ApJ , 438, 763 CrossRef Search ADS   Gough D., 1977, in Spiegel E. A., Zahn J.-P., eds, Astrophysics and Space Science Library, Vol. 71, Problems of Stellar Convection . Springer-Verlag, Berlin p. 15 Gough D. O., 2012, A&A , 2012, 987275 Guennebaud G. et al.  , 2010. Available at http://eigen.tuxfamily.org Hunter J. D., 2007, Comput. Sci. Eng. , 9, 90 CrossRef Search ADS   Isserlis L., 1918, Biometrika , 12, 134 CrossRef Search ADS   Käpylä P. J., Korpi M. J., Tuominen I., 2004, A&A , 422, 793 CrossRef Search ADS   Kichatinov L. L., 1986, Geophys. Astrophys. Fluid Dyn. , 35, 93 CrossRef Search ADS   Kichatinov L. L., 1987, Geophys. Astrophys. Fluid Dyn. , 38, 273 CrossRef Search ADS   Kichatinov L. L., Rudiger G., 1993, A&A , 276, 96 Killworth P. D., 1980, Dyn. Atmos. Oceans , 4, 143 CrossRef Search ADS   King A. R., Pringle J. E., Livio M., 2007, MNRAS , 376, 1740 CrossRef Search ADS   Kitchatinov L. L., 2013, in Kosovichev A. G., de Gouveia Dal Pino E., Yan Y., eds, Proc. IAU Symp. 294, Solar and Astrophysical Dynamos and Magnetic Activity . Kluwer, Dordrecht, p. 399 Kolmogorov A., 1941a, Akademiia Nauk SSSR Doklady , 30, 301 Kolmogorov A. N., 1941b, Akademiia Nauk SSSR Doklady , 32, 16 Launder B., Spalding D., 1974, Comput. Methods Appl. Mech. Eng. , 3, 269 CrossRef Search ADS   Lebovitz N. R., 1965, ApJ , 142, 1257 CrossRef Search ADS   Lee D., 2013, J. Comput. Phys. , 243, 269 CrossRef Search ADS   Lesaffre P., Balbus S. A., Latter H., 2009, MNRAS , 396, 779 CrossRef Search ADS   Lesaffre P., Chitre S. M., Potter A. T., Tout C. A., 2013, MNRAS , 431, 2200 CrossRef Search ADS   Lohse D., Xia K.-Q., 2010, Annu. Rev. Fluid Mech. , 42, 335 CrossRef Search ADS   Luschgy H., Pagès G., 2006, preprint (arXiv:math/0607282) Maeder A., 1995, A&A , 299, 84 McKinney J. C., Tchekhovskoy A., Sadowski A., Narayan R., 2014, MNRAS , 441, 3177 CrossRef Search ADS   Murphy G. C., Pessah M. E., 2015, ApJ , 802, 139 CrossRef Search ADS   Pessah M. E., Goodman J., 2009, ApJ , 698, L72 CrossRef Search ADS   Pessah M. E., Chan C.-K., Psaltis D., 2006, MNRAS , 372, 183 CrossRef Search ADS   Pessah M. E., Chan C.-K., Psaltis D., 2008, MNRAS , 383, 683 CrossRef Search ADS   Procaccia I., Zeitak R., 1989, Phys. Rev. Lett. , 62, 2128 CrossRef Search ADS PubMed  Rüdiger G., Egorov P., Ziegler U., 2005a, Astron. Nachr. , 326, 315 CrossRef Search ADS   Rüdiger G., Egorov P., Kitchatinov L. L., Küker M., 2005b, A&A , 431, 345 CrossRef Search ADS   Salvesen G., Simon J. B., Armitage P. J., Begelman M. C., 2016, MNRAS , 457, 857 CrossRef Search ADS   Schou J. et al.  , 1998, ApJ , 505, 390 CrossRef Search ADS   Smolec R., Houdek G., Gough D., 2011, in Brummell N. H., Brun A. S., Miesch M. S., Ponty Y., eds, Proc. IAU Symp. 271, Astrophysical Dynamics: From Stars to Galaxies . Kluwer, Dordrecht, p. 397 Sridhar S., Goldreich P., 1994, ApJ , 432, 612 CrossRef Search ADS   Starling R. L. C., Siemiginowska A., Uttley P., Soria R., 2004, MNRAS , 347, 67 CrossRef Search ADS   Sukoriansky S., Dikovskaya N., Galperin B., 2007, J. Atmos. Sci. , 64, 3312 CrossRef Search ADS   van der Walt S., Colbert S. C., Varoquaux G., 2011, Comput. Sci. Eng. , 13, 22 CrossRef Search ADS   Wick G. C., 1950, Phys. Rev. , 80, 268 CrossRef Search ADS   Yakhot V., Orszag S. A., 1986, J. Sci. Comput. , 1, 3 CrossRef Search ADS   Zahn J.-P., 1993, Space Sci. Rev. , 66, 285 CrossRef Search ADS   Zhou Y., McComb D. W., Vahala G., 1997, NASA Contractor Report 201718  APPENDIX A: TURBULENT INDEX The general question of which turbulent index to use and under what circumstances remains open though many specific cases are well understood. In the case of isotropic incompressible turbulence the Kolmogorov index is well-known to be n = 11/6 (Kolmogorov 1941a). There is more debate over the index to use for convection, with answers ranging from n = 5/2 (Benzi et al. 1994) to n = 21/10 (Procaccia & Zeitak 1989) and n = 2.4 ± 0.2 (Ashkenazi & Steinberg 1999). There has also been work attempting to determine the spectrum in a context-sensitive manner through energy balance arguments (Yakhot & Orszag 1986). In the magnetised case sources differ even more, with some suggesting that this range still applies (Dobrowolny, Mangeney & Veltri 1980), some arguing for a Kolmogorov-like spectrum (Goldreich & Sridhar 1995) and others giving a range of indices depending on geometry and the direction of the wavevector (Sridhar & Goldreich 1994). From numerical experiments with our closure model, we have found that the magnetic stress scales sufficiently rapidly with k that it is divergent for n = 11/6 and not for n = 8/3. This favours the scenario of Goldreich & Sridhar (1995), who argue that in the strongly magnetized limit the index ought to be n = 8/3. In order to consistently treat both the non-magnetic and the strongly magnetized limits, we choose a simple prescription in which n = 11/6 when one of |N|, or |R∇Ω| exceeds kvA and use n = 8/3 otherwise. This means that there is a critical wavenumber   \begin{eqnarray} k_{\rm c} \equiv \frac{\max \left(|N|, |R\nabla \Omega |\right)}{v_A} \end{eqnarray} (A1)at which the spectrum changes. In the non-magnetic case, the evolution matrix is independent of the magnitude of the wavevector and so altering the index just alters the correlation coefficients by a multiplicative factor. In the magnetic case, the potential for error is larger because the magnitude of the wavevector is relevant but there appears to be no consensus on the best prescription and so we make do with what is available. APPENDIX B: BOUSSINESQ ODDITIES In this work, we have taken the Boussinesq approximation. In Fourier space this is   \begin{eqnarray} \boldsymbol {q}\cdot \tilde{\delta \boldsymbol {r}}=0. \end{eqnarray} (B1)Taking the time derivative of both sides we see that   \begin{eqnarray} \partial _t\left(\boldsymbol {q}\cdot \delta \tilde{\boldsymbol {r}}\right) = \boldsymbol {q}\cdot \delta \tilde{\boldsymbol {v}} + \delta \tilde{\boldsymbol {r}}\cdot \partial _t\boldsymbol {q} = 0. \end{eqnarray} (B2)As a result   \begin{eqnarray} \delta \tilde{\boldsymbol {v}}\cdot \boldsymbol {q} = - \delta \tilde{\boldsymbol {r}}\cdot \partial _t\boldsymbol {q} \ne 0. \end{eqnarray} (B3)This is quite peculiar, but is just an artefact of our coordinate system. Because the wavevectors are time-dependent, maintaining the volume of a fluid parcel requires that the displacement be orthogonal to the wavevector, which actually means that the velocity is generally not orthogonal to the wavevector. APPENDIX C: SOFTWARE DETAILS The software used for this work is Mixer version 1, which we have released under a GPLv3 license at github.com/adamjermyn/Mixer. All data produced for this work are available at the same location as HDF5 tables with attributes documenting the physical inputs. Post-processing and visualization of the data was with the Python modules Numpy (van der Walt, Colbert & Varoquaux 2011) and Matplotlib (Hunter 2007) and the relevant scripts for this are included with Mixer. The core of Mixer is written in C++, for performance reasons, and the code is supplied with a Makefile, which supports compilation on both Linux and MacOS. Mixer makes use of the Eigen library (Guennebaud et al. 2010) for linear algebra. Mixer also uses the Cubature library for numerical integration. This library is an implementation of the algorithms by Genz & Malik (1980) and Berntsen, Espelid & Genz (1991). These integration routines are supplemented by a Python integration routine tailored for integrands with small support regions. The details will be explored in later work. In addition, many routines provide a Python interface. Currently Mixer supports only single-threaded operation, though it may be used inside parallelized scripts through the Python wrapper. The version of Mixer used to generate the data in this work was compiled against Cubature version 1.0.2 and Eigen version 3.3.3, though the code does not use any features which require recent versions, so many likely suffice. Mixer is optimized for convecting systems for which achieving accuracy better than 10−5 relative and absolute typically requires between 1 ms and 1 s on a single core of a 2016 Intel CPU. This is further improved when the differential rotation is minimal, in which case the perturbative expansion may be turned off to save a factor of several in runtime. In stably stratified zones and those with magnetic fields up to 103s may be required to achieve good convergence. In cases where the code has more difficulty, it is quite likely that Mixer becomes the bottleneck in simulations and so, under these circumstances, we recommend tabulating results in advance. This is still considerably more performant than direct numerical simulation, and the results can generally be guaranteed to converge at much higher precision, so that derivatives may be extracted as well. At various points in the software, we must divide by the magnitude of the velocity of an eigenmode. This may approach zero in some cases. To avoid dividing by zero in these cases, we place a lower bound on this magnitude, such that   \begin{eqnarray} |\delta v|^2 \ge \epsilon , \end{eqnarray} (C1)where $$\epsilon = 10^{-20} L_0^2 |N|^2$$ in the calculations presented in this work. This corresponds to setting an upper bound on the length scale d of the displacements $$\delta \boldsymbol {r}$$, namely   \begin{eqnarray} |\delta r|^2 \le L_0^3 |N| \epsilon ^{-1/2}, \end{eqnarray} (C2)which means that d = 1010L0 in this work. To verify that this numerical fix does not impact our results, we have examined the correlation functions in several scenarios as a function of this numerical cut off L. For example, Fig. C1 shows the r − θ and r − ϕ correlations as functions of d for a stably stratified differentially rotating system. The results are constant over many orders of magnitude so long as d > 104L0, which is easily satisfied by our default. Figure C1. View largeDownload slide The absolute values of 〈δvrδvθ〉 (left-hand panel) and 〈δvrδvϕ〉 (right-hand panel) are shown as functions of d, with both axes log-scaled. These results are for a stably stratified region with differential rotation in the radial direction with |R∇ln Ω| = 10−3, Ω = 0.1|N| and no magnetic field. The data is computed for a point on the equator with differential rotation at an angle of π/4. All quantities are given in units of the mixing length and the Brünt–Väisälä frequency. Figure C1. View largeDownload slide The absolute values of 〈δvrδvθ〉 (left-hand panel) and 〈δvrδvϕ〉 (right-hand panel) are shown as functions of d, with both axes log-scaled. These results are for a stably stratified region with differential rotation in the radial direction with |R∇ln Ω| = 10−3, Ω = 0.1|N| and no magnetic field. The data is computed for a point on the equator with differential rotation at an angle of π/4. All quantities are given in units of the mixing length and the Brünt–Väisälä frequency. © 2018 The Author(s) Published by Oxford University Press on behalf of the Royal Astronomical Society http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Monthly Notices of the Royal Astronomical Society Oxford University Press

# Turbulence closure for mixing length theories

, Volume 476 (1) – May 1, 2018
17 pages

/lp/ou_press/turbulence-closure-for-mixing-length-theories-GTXDAuD4Yo
Publisher
The Royal Astronomical Society
© 2018 The Author(s) Published by Oxford University Press on behalf of the Royal Astronomical Society
ISSN
0035-8711
eISSN
1365-2966
D.O.I.
10.1093/mnras/sty255
Publisher site
See Article on Publisher Site

### Abstract

Abstract We present an approach to turbulence closure based on mixing length theory with three-dimensional fluctuations against a two-dimensional background. This model is intended to be rapidly computable for implementation in stellar evolution software and to capture a wide range of relevant phenomena with just a single free parameter, namely the mixing length. We incorporate magnetic, rotational, baroclinic, and buoyancy effects exactly within the formalism of linear growth theories with non-linear decay. We treat differential rotation effects perturbatively in the corotating frame using a novel controlled approximation, which matches the time evolution of the reference frame to arbitrary order. We then implement this model in an efficient open source code and discuss the resulting turbulent stresses and transport coefficients. We demonstrate that this model exhibits convective, baroclinic, and shear instabilities as well as the magnetorotational instability. It also exhibits non-linear saturation behaviour, and we use this to extract the asymptotic scaling of various transport coefficients in physically interesting limits. convection, stars: evolution, stars: interiors, stars: rotation 1 INTRODUCTION An understanding of turbulent transport and stresses remains one of the major outstanding problems in the astrophysics of fluids. While many pieces of this puzzle are understood in broad strokes, the nature of this problem is such that the details are almost as important as the big picture. The magnetorotational instability (MRI), for instance, is understood conceptually but making predictions that match observed accretion discs is a persistent problem (Murphy & Pessah 2015). Similarly, the solar differential rotation is understood to arise from turbulent stresses but precisely how this works and in balance with what other forces remains uncertain (Schou et al. 1998). Significant progress has indeed been made with three-dimensional turbulence simulations (for examples, see Lee 2013; McKinney et al. 2014; Salvesen et al. 2016) but these are generally relevant only on short time-scales and in small volumes. Performing so-called global simulations over large times and distances requires a turbulence closure model to substitute for resolution at small scales (Launder & Spalding 1974; Canuto 1994). At the other extreme models of stellar evolution generally assume extremely simple analytical transport coefficients to overcome the tremendous gap between turbulent time-scales of minutes and nuclear time-scales of millions of years (Maeder 1995). A variety of such approaches have been developed. For instance, the mixing length theory of Böhm-Vitense (1958) provided a closure of convection. This was then put on firmer theoretical ground by Gough (1977, 2012) and extended to include additional phenomena (Smolec, Houdek & Gough 2011; Lesaffre et al. 2013). Kichatinov (1986) introduced an entirely different closure formalism, arriving at an expression for the so-called Λ-effect (Kichatinov 1987), and later incorporating it under the α–Λ formalism with Rudiger (Kichatinov & Rudiger 1993). What these formalisms have in common is a minimal set of free parameters: the mixing length formalism has just the mixing length, and the formalism of Kichatinov & Rudiger (1993) has just the anisotropy parameter. Another set of models has arisen, which aims to reproduce higher order moments of the turbulent fields. This increases the number of free parameters and a number of approaches have been developed to deal with this. For instance, Garaud, Gagnier & Verhoeven (2017) and Garaud et al. (2010) fit their free parameters against small-scale simulations, while Canuto (1997) fits his against experimental results. In addition, there are models, such as that of Canuto (1994), which fix at least some free parameters by introducing new assumptions, in this case regarding the various relevant time-scales. Regardless of the details of how they close the equations of turbulent moments, models of this sort generally take the form of physically motivated analytic expressions, which provide ready access to scaling laws. Their free parameters then serve to better their agreement with data, at the cost of being less straightforwardly interpreted and extended. The availability of growing computational resources in recent years has provided a new niche in this landscape in the form of computational closure models. These are models which do not seek analytic solutions but which are none the less distinct from attempts to simulate turbulence in all its detail. Some may introduce new dynamical fields, as in the k − ε model (Launder & Spalding 1974), while others invoke effective theories of small-scale motion (Canuto & Hartke 1986). The latter kind are essentially renormalized theories, which accept the cost of having to numerically accommodate complex behaviour in exchange for more precision over a wider variety of phenomena. Combined with perturbation theory, this approach represents a tunable middle-ground between expensive simulations and simple analytic models, allowing the computational cost to be traded off against fidelity to suit the problem at hand. The model we present here is in this spirit. We construct a mixing-length theory, which incorporates three-dimensional fluctuations against a two-dimensional axisymmetric background. This is done by treating each mode as growing with its linear growth rate before saturating at an amplitude set by the turbulent cascade (Lesaffre et al. 2013). Beyond this, the motion in each mode is taken to be uncorrelated. We treat the geometry of the flow in full generality, allowing for baroclinic effects as well as magnetism and rotational shears. To incorporate differential rotation, we use a time-dependent sheared coordinate system (Balbus & Schaan 2012). In this frame, there is a continual flow of modes across Fourier space, lending a time dependence to growth rates. Corrections to saturation amplitudes owing to this flow are incorporated perturbatively with the time derivatives of the growth rate. In Section 2, we describe our closure framework in more detail, paying particular attention to the choice of mixing length. We then develop a perturbative approach for correcting the saturation amplitude in Section 3. In section 4, we introduce the sheared coordinate system and the linearized equations of motion. Finally, in Section 6, we show results from our theory, including calculations for the solar convection zone and accretion discs. The software implementing our model is open source and available under a GPLv3 license. Details of the implementation are given in Appendix C. Tabulated transport coefficients produced by the code are also available under the same license and both may be found at github.com/adamjermyn/Mixer. 2 CLOSURE FORMALISM Turbulent phenomena generically exhibit a cascade of energy between large and small scales (Zhou, McComb & Vahala 1997; Lohse & Xia 2010). With some notable exceptions (Sukoriansky, Dikovskaya & Galperin 2007), this cascade begins at a large scale L0 set by the overall structure of the fluid flow and ends at an extremely small scale Lν related to the microscopic viscosity. Between these scales, yet far from each of them, lies the so-called inertial range where the fluid flow is scale-free (Kolmogorov 1941b). In this range, all correlations of the turbulent motion obey simple power laws. This statement was originally proved by Kolmogorov (1941b) for isotropic turbulence. It was later found to be a broader consequence of the renormalizability of the Navier–Stokes equation (Yakhot & Orszag 1986; Carati 1990) and consequently holds quite generally. This means that there is a single relevant scale L0 for a given turbulent flow, which fully characterises the turbulence as seen by measurements performed over length scales L ≫ L0. This is the modern interpretation and justification of the original mixing length hypothesis, which asserts that turbulent fluctuations on scales L ≪ L0 are not dynamically coupled to the large-scale (L ≫ L0) flow properties (Böhm-Vitense 1958). The scale-free nature of turbulence in the inertial range means that modes of significantly different wavevectors are uncorrelated. A natural extension of this is to assume that all modes of distinct wavevectors are at least approximately uncorrelated. That is, we assume that   \begin{eqnarray} \langle \tilde{\boldsymbol {v}}_{\boldsymbol {k}} \otimes \tilde{\boldsymbol {v}}^*_{\boldsymbol {k}^{\prime }} \rangle &= (2{\pi} )^3 \delta ^3 (\boldsymbol {k}-\boldsymbol {k}^{\prime }) \mathsf {V}_{\boldsymbol {k}}, \end{eqnarray} (1)where v is the velocity, ⊗ denotes the outer product, 〈⋅⋅⋅〉 denotes the time-averaged expectation, $$\tilde{\boldsymbol {v}}_{\boldsymbol {k}}$$ is the amplitude of the Fourier mode with wavevector $$\boldsymbol {k}$$ and $$\mathsf {V}_{\boldsymbol {k}}$$ is the tensor specifying how different components of the same mode are correlated with one another. It is crucial to notice that the quantity $$\mathsf {V}_{\boldsymbol {k}}$$ is also the Reynolds stress of mode $$\boldsymbol {k}$$. This, and several other closely related quantities, are ultimately what we seek. These two-point correlation functions suffice to characterize not only the stresses but also all higher order correlations through Wick’s theorem and perturbation theory (Wick 1950; Isserlis 1918). To determine $$\mathsf {V}_{\boldsymbol {k}}$$, we begin by writing the linearized equations of motion as   \begin{eqnarray} \partial _t \boldsymbol {v}(\boldsymbol {r}) = \mathcal {L}\left[\boldsymbol {v}(\boldsymbol {r}), \partial _i \boldsymbol {v}, \partial _i \partial _j \boldsymbol {v},\ldots , \boldsymbol {r}, t\right], \end{eqnarray} (2)where $$\mathcal {L}$$ is a linear operator of its first argument and $$\boldsymbol {v}$$ is the fluctuating part of the velocity field. In principle, we can work with this operator, though the derivatives of the velocity field make it highly inconvenient. Fortunately, at short length scales, the operator $$\mathcal {L}$$ may be treated as translation-invariant and so we may compute a Fourier transform in $$\boldsymbol {r}$$ without coupling different modes. This gives   \begin{eqnarray} \frac{{\rm d}\tilde{\boldsymbol {v}}_{\boldsymbol {k}}}{{\rm d}t} = \tilde{\mathcal {L}}\left[\tilde{\boldsymbol {v}}_{\boldsymbol {k}}, \boldsymbol {k}, t\right]. \end{eqnarray} (3)The modes are decoupled in this regime so $$\tilde{\mathcal {L}}$$ can be represented by a matrix $$\mathsf {L}$$, and we write   \begin{eqnarray} \frac{{\rm d}\tilde{\boldsymbol {v}}_{\boldsymbol {k}}}{{\rm d}t} = \mathsf {L}(\boldsymbol {k},t) \boldsymbol {v}_{\boldsymbol {k}}. \end{eqnarray} (4) When $$\mathsf {L}$$ is independent of t equation (4) is straightforward to solve and gives us   \begin{eqnarray} \frac{{\rm d}\tilde{\boldsymbol {v}}_{\boldsymbol {k}}}{{\rm d}t} = \sum _i v_{0,i} \hat{v}_{\boldsymbol {k},i} e^{\lambda _i t}, \end{eqnarray} (5)where v0, i are the initial mode amplitudes and $$\hat{v}_{\boldsymbol {k},i}$$ and λi are the normalized right eigenvectors and eigenvalues of $$\mathsf {L}$$, respectively. The vectors $$\hat{v}_{\boldsymbol {k},i}$$ then specify the modes of the system at a given wavevector. If the eigenvalues are not precisely degenerate then modes which begin in phase rapidly become uncorrelated and we may extend equation (1) to the modes at each wavevector and write   \begin{eqnarray} \langle \tilde{\boldsymbol {v}}_{\boldsymbol {k},i} \otimes \tilde{\boldsymbol {v}}^*_{\boldsymbol {k}^{\prime },j} \rangle &= (2{\pi} )^3 \delta ^3 (\boldsymbol {k}-\boldsymbol {k}^{\prime }) \delta _{ij} \mathsf {V}_{\boldsymbol {k},i}. \end{eqnarray} (6)This result holds even when modes are degenerate. Because the time evolution of the Navier–Stokes equation is deterministic, the expectation 〈⋅⋅⋅〉 represents a sum over initial conditions. In this sum all relative phases between the modes are explored, so even degenerate modes become uncorrelated. Inserting equation (5) into equation (6) and summing over j and integrating over $$\boldsymbol {k}$$ gives us   \begin{eqnarray} \mathsf {V}_{\boldsymbol {k},i} &= \hat{v}_{\boldsymbol {k},i} \otimes \hat{v}_{\boldsymbol {k},i} \langle |v_{0,i}|^2 \exp \left[2 t \Re \left[\lambda _i\right]\right] \rangle . \end{eqnarray} (7)Generally, some λi have positive real parts and so in a long-term expectation this exponential diverges. Indeed, it turns out that these growing modes are precisely those which matter! What happens of course is just that these modes eventually reach amplitudes where the linear approximation fails. By assumption, the system is stable over long times relative to the turbulent scale so this must result in these modes saturating. This has been variously described as mode crashing or the action of parasitic modes (Lesaffre, Balbus & Latter 2009; Pessah & Goodman 2009) but, regardless of the mechanism, it simply means that these modes exit the linear regime and find their growth impeded. To complete the closure, we must find the saturation amplitude. Relying again on the scale-free nature of turbulence, we note that this must be a power law in k. That is   \begin{eqnarray} \langle \tilde{v}_{\boldsymbol {k}, i}^2 \rangle = \mathrm{Tr}\left[\mathsf {V}_{\boldsymbol {k}, i}\right] = \frac{A}{M} \left(\frac{k_0}{k}\right)^n, \end{eqnarray} (8)where A depends on the large-scale properties of the flow but is independent of $$\boldsymbol {k}$$, M is the number of modes per wavevector and n is the index of the turbulence. Following Kolmogorov (1941a), we choose n = 11/6 in our model. Appendix A contains a detailed discussion of this choice. The wavenumber k0 is just that of the characteristic scale, and is given by   \begin{eqnarray} k_0 = \frac{2{\pi} }{L_0}. \end{eqnarray} (9)Replacing the divergent expression in equation (7) with this amplitude, we find   \begin{eqnarray} \mathsf {V}_{\boldsymbol {k},i} = \frac{A}{M} \left(\frac{k_0}{k}\right)^n \hat{v}_{\boldsymbol {k},i} \otimes \hat{v}_{\boldsymbol {k},i}. \end{eqnarray} (10) It only remains to determine A. To do this, we note that there is one characteristic length scale L0 and one characteristic time scale, the growth rate ℜ[λi] of the mode. Because A has dimensions of velocity squared, we find   \begin{eqnarray} \mathsf {V}_{\boldsymbol {k},i} = \frac{c}{M} L_0^2 \Re \left[\lambda _i\right]^2 \left(\frac{k_0}{k}\right)^n \hat{v}_{\boldsymbol {k},i} \otimes \hat{v}_{\boldsymbol {k},i}, \end{eqnarray} (11)where c is a dimensionless constant of order unity. This constant, known as the mixing length parameter, varies from theory to theory, so for clarity we set c = 1 in this work but this degree of freedom is important to note when comparing between models. In effect, what we have done is incorporate the non-linearity of turbulence by means of the spectrum while using linear growth rates to set the characteristic scale. In practice, the spectrum acts only to provide a convergent measure over modes (see Appendix A for further discussion), and it is the growth rate and the modes themselves that yield the anisotropies and other phenomena of interest. This is closely related to the approaches of Lesaffre et al. (2013) and Canuto & Hartke (1986). This prescription is easily extended in cases where there are additional dynamical fields, such as the turbulent displacement or a fluctuating magnetic field. The additional fields are simply incorporated into the vector describing the state and M is increased accordingly. We can continue to use equation (8) to fix the amplitude of the entire mode against that of the velocity as long as we know the turbulent index n. Up to this point this prescription is mathematically identical to that of Lesaffre et al. (2013), with the exception that we define the mixing wave vector as in equation (9), while they use π/L0 instead. In the next section, we introduce perturbative corrections to this model to capture a wider variety of phenomena. 3 PERTURBATIVE CORRECTIONS Now consider the case where the matrix $$\mathsf {L}$$ is time-dependent. Most of our reasoning about the behaviour of modes from the previous section still holds but, because the eigenvectors are time-dependent, we no longer have a well-defined notion of a mode as a long-running solution to the equations of motion. When the time dependence is periodic Floquet theory applies (Floquet 1883), but in the cases of interest the time dependence is aperiodic. To recover modes when the time evolution matrix itself evolves and does so aperiodically, we begin by expanding as   \begin{eqnarray} \mathsf {L}(t) = \mathsf {L}(0) + t \frac{{\rm d}\mathsf {L}}{{\rm d}t} + \frac{1}{2}t^2 \frac{{\rm d}^2\mathsf {L}}{{\rm d}t^2} +\cdots \,. \end{eqnarray} (12)This series can be truncated to produce an approximation of $$\mathsf {L}$$, which is accurate in a certain window around t = 0. We may likewise write the velocity at a given wavevector as   \begin{eqnarray} \tilde{\boldsymbol {v}}_{\boldsymbol {k}}(t) = \tilde{\boldsymbol {v}}_{\boldsymbol {k}}(0) + t\left. \frac{{\rm d}\tilde{\boldsymbol {v}}_{\boldsymbol {k}}}{{\rm d}t}\right|_0 + \frac{1}{2}t^2\left. \frac{{\rm d}^2\tilde{\boldsymbol {v}}_{\boldsymbol {k}}(t)}{{\rm d}t^2}\right|_0 +\cdots \,. \end{eqnarray} (13)This suggests defining a new vector   \begin{eqnarray} \Phi _{\boldsymbol {k}}(t) \equiv \left\lbrace \tilde{\boldsymbol {v}}_{\boldsymbol {k}}, \frac{{\rm d}\tilde{\boldsymbol {v}}_{\boldsymbol {k}}}{{\rm d}t}, \frac{{\rm d}^2\tilde{\boldsymbol {v}}_{\boldsymbol {k}}}{{\rm d}t^2},\ldots \right\rbrace , \end{eqnarray} (14)which, in principle, encodes the full time evolution of the velocity field. This vector evolves according to   \begin{eqnarray} \frac{{\rm d}\Phi _{\boldsymbol {k}}}{{\rm d}t} = \mathsf {A} \Phi _{\boldsymbol {k}}, \end{eqnarray} (15)where $$\mathsf {A}$$ is formed of blocks given by   \begin{eqnarray} \mathsf {A}_{ij} = {i \atopwithdelims ()j}\frac{{\rm d}^{i-j}}{{\rm d}t^{i-j}} \mathsf {L}. \end{eqnarray} (16)By definition though we also have   \begin{eqnarray} \frac{{\rm d}\Phi _{\boldsymbol {k},i}}{{\rm d}t} = \Phi _{\boldsymbol {k},i+1}, \end{eqnarray} (17)where $$\Phi _{\boldsymbol {k},0} = \tilde{\boldsymbol {v}}_{\boldsymbol {k}}$$, $$\Phi _{\boldsymbol {k},1} = {\rm d}\tilde{\boldsymbol {v}}_{\boldsymbol {k}}/{\rm d}t$$ and so on. Thus, we are searching for a simultaneous solution of equations (15) and (17). In order to close the system, we must truncate it at some finite order N. Doing so makes the assumption that the behaviour of the system at all greater N is known. Inspired by the solution for time-independent $$\mathsf {L}$$, we try an exponential behaviour. This truncates equation (17) such that it applies only to i < N − 1 and means that we are searching for vectors with   \begin{eqnarray} \left(A \Phi _{\boldsymbol {k}}\right)_{N-1} = \lambda \Phi _{\boldsymbol {k},N-1} \end{eqnarray} (18)and   \begin{eqnarray} \Phi _{\boldsymbol {k},i+1} = (\mathsf {A} \Phi _{\boldsymbol {k}})_i, i < N-1. \end{eqnarray} (19)These equations are most straightforwardly written as a general eigensystem and this has the advantage of restricting the dimension of the linear space to just those states obeying the constraint. This is possible because both $$\mathsf {A}$$ and the constraint are lower triangular in the same basis, and so each row may be substituted into the next, leading to an eigenproblem of the form   \begin{eqnarray} \mathsf {Q} \Phi _{\boldsymbol {k,0}} = \lambda \mathsf {W} \Phi _{\boldsymbol {k,0}}, \end{eqnarray} (20)where $$\mathsf {Q}$$ and $$\mathsf {W}$$ are matrices acting only on the 0-block. For example, in the case where N = 2, our equations are   \begin{eqnarray} \Phi _{\boldsymbol {k},1} &= \mathsf {M} \Phi _{\boldsymbol {k},0} \end{eqnarray} (21)and   \begin{eqnarray} \mathsf {M} \Phi _{\boldsymbol {k},1} + \dot{\mathsf {M}} \Phi _{\boldsymbol {k},0} &= \lambda \Phi _{\boldsymbol {k},1}, \end{eqnarray} (22)which may be put in the form of equation (20) with   \begin{eqnarray} \mathsf {Q} &= \mathsf {M}^2 + \dot{\mathsf {M}} \end{eqnarray} (23)and   \begin{eqnarray} \mathsf {W} &= \mathsf {M}. \end{eqnarray} (24)The eigenvectors of this system are solutions of the original equation (4) because if $$\psi ^i_{\boldsymbol {k}}$$ is such an eigenvector then   \begin{eqnarray} \tilde{\boldsymbol {v}}_{\boldsymbol {k},i}(t) \equiv \sum _{j=0}^{N} \frac{t^j}{j!}\psi ^i_{j} \end{eqnarray} (25)solves   \begin{eqnarray} \frac{{\rm d}\tilde{\boldsymbol {v}}_{\boldsymbol {k},i}}{{\rm d}t} = \mathsf {L}(t) \tilde{\boldsymbol {v}}_{\boldsymbol {k},i}(t) \end{eqnarray} (26)over the time window for which $$\mathsf {L}$$ is well-approximated at N-th order. As a result, we say that ϕi(t) are the instantaneous modes of the system at N-th order and use them and in equation (11). In place of the eigenvalue, we use the instantaneous growth rate of the velocity, which is given by   \begin{eqnarray} g \equiv \frac{1}{2}\frac{{\rm d} v^2}{{\rm d}t} = \frac{\Re \left(\Phi _{\boldsymbol {k},0} \cdot \Phi _{\boldsymbol {k},1}\right)}{|\Phi _{\boldsymbol {k},0}|^2}. \end{eqnarray} (27)This approximation is controlled in the sense that so long as $$\mathsf {L}(t)$$ converges as N grows, so does the inferred velocity history. In this work, we present results with N = 2 so that $$\mathsf {A}$$ involves both $$\mathsf {L}$$ and $$\dot{\mathsf {L}}$$. We leave the exploration of larger N to later work. 4 EQUATIONS OF MOTION We now specialise to the case of an ideal gas obeying the ideal MHD equations. This section largely follows the derivation of Balbus & Schaan (2012) so we present only the pieces necessary to understand later parts of this work as well as the few places where our derivation diverges from theirs. We take the background to be axisymmetric, the fluctuations to be adiabatic and we work in cylindrical coordinates. We neglect both the microscopic viscosity and the microscopic thermal diffusivity because these are both negligible in most circumstances in stellar physics.1 Because our closure model treats turbulent properties as local, we compute all background quantities at a reference point $$\boldsymbol {r}_0$$. Relative to this point, we define the Lagrangian separation $$\boldsymbol {\delta r}$$ and velocity $$\boldsymbol {\delta v}$$ equivalent to $$\boldsymbol {\xi }$$ and $${\rm D}\boldsymbol {\xi }/{\rm D}t$$ of Balbus & Schaan (2012). In addition, we take the Boussinesq approximation that density variations are ignored except in terms involving gravitational acceleration. With the above definitions, the continuity equation may be written as   \begin{eqnarray} \nabla \cdot \boldsymbol {\delta r} = 0. \end{eqnarray} (28) In a fixed coordinate system differential rotation is difficult to analyse so we make two reference frame changes. First, we switch from an inertial frame to one rotating at   \begin{eqnarray} \Omega _0 \equiv \Omega (\boldsymbol {r}_0). \end{eqnarray} (29)Secondly, we make a formal change of coordinates   \begin{eqnarray} \phi \rightarrow \phi - t\delta \boldsymbol {r}\cdot \nabla \Omega , \end{eqnarray} (30)without altering the corresponding unit vectors. Under this last change the gradient transforms as   \begin{eqnarray} \nabla &\rightarrow \nabla - t(\nabla \Omega )\partial _{\phi }. \end{eqnarray} (31)Because the operator $$\mathcal {L}$$ is most easily expressed in Fourier space, we define the transformed wavevector as   \begin{eqnarray} \boldsymbol {q} \equiv \boldsymbol {k} - t k_\phi R \nabla \Omega . \end{eqnarray} (32)With this the transformed MHD and Navier–Stokes equations may be written as   \begin{eqnarray} \boldsymbol {\delta \tilde{B}} = \boldsymbol {B} \cdot \boldsymbol {q} \boldsymbol {\delta \tilde{r}} \end{eqnarray} (33)and   \begin{eqnarray} {\partial _t \delta \tilde{\boldsymbol {v}} + 2\boldsymbol {\Omega }\times \delta \tilde{\boldsymbol {v}}+ \hat{R} R \delta \tilde{\boldsymbol {r}}\cdot \nabla \Omega ^2}\nonumber\\ {\quad- \frac{1}{\gamma \rho }\left(\boldsymbol {\delta \tilde{r}}\cdot \nabla \sigma \right)\nabla \cdot \boldsymbol {\Pi } + \frac{i}{\rho }\boldsymbol {q}\cdot \boldsymbol {\delta \tilde{\Pi }}=0,} \end{eqnarray} (34)where σ is the specific entropy and   \begin{eqnarray} \boldsymbol {\Pi } \equiv p \mathsf {I} - \frac{1}{\mu _0}\left(\boldsymbol {B}\otimes \boldsymbol {B} - \frac{1}{2}B^2 \mathsf {I}\right) \end{eqnarray} (35)is the pressure tensor with $$\mathsf {I}$$ the identity matrix. All quantities prefixed with δ are fluctuating, a tilde denotes the Fourier transformed function, and all other quantities are background fields evaluated at $$\boldsymbol {r}_0$$. It is straightforward to see that this is the same equation as that derived by Balbus & Schaan (2012) once the appropriate relations for the pressure and magnetic force are substituted. The fluctuation in the pressure tensor may be written as   \begin{eqnarray} \boldsymbol {\delta \Pi } = \delta p \mathsf {I} - \frac{1}{\mu _0}\left(\boldsymbol {B}\otimes \boldsymbol {\delta B} + \boldsymbol {\delta B}\otimes \boldsymbol {B} - \mathsf {I} \boldsymbol {B}\cdot \boldsymbol {\delta B} \right), \end{eqnarray} (36)so in Fourier space   \begin{eqnarray} \boldsymbol {\delta \tilde{\Pi }} = \delta \tilde{p} \mathsf {I} - \frac{1}{\mu _0}\left(\boldsymbol {B}\otimes \boldsymbol {\delta \tilde{B}} + \boldsymbol {\delta \tilde{B}}\otimes \boldsymbol {B} - \mathsf {I} \boldsymbol {B}\cdot \boldsymbol {\delta \tilde{B}} \right). \end{eqnarray} (37)Combining this with equation (33) and the Boussinesq approximation (see Appendix B), we find   \begin{eqnarray} \boldsymbol {q}\cdot \boldsymbol {\delta \tilde{\Pi }} = \boldsymbol {q} \delta p - \frac{i}{\mu _0} (\boldsymbol {B}\cdot \boldsymbol {q})^2 \boldsymbol {\delta \tilde{r}}. \end{eqnarray} (38)Note that as did Balbus & Schaan (2012), we take $$\boldsymbol {B}\cdot \boldsymbol {q}$$ to be constant in time as implied by the Boussinesq and ideal-MHD conditions. We now depart from prior work and use this equation along with equation (34) taking the component perpendicular to $$\boldsymbol {q}$$ to eliminate δp and find   \begin{eqnarray} 0 &=&\left(\partial _t \delta \tilde{\boldsymbol {v}} + 2\boldsymbol {\Omega }\times \delta \tilde{\boldsymbol {v}}+ \hat{R} R \delta \tilde{\boldsymbol {r}}\cdot \nabla \Omega ^2\right. \nonumber \\ &&- \left. \frac{1}{\gamma \rho }\left(\boldsymbol {\delta \tilde{r}}\cdot \nabla \sigma \right)\nabla \cdot \boldsymbol {\Pi } + \frac{1}{\mu _0 \rho } (\boldsymbol {B}\cdot \boldsymbol {q})^2 \boldsymbol {\delta \tilde{r}}\right)_{\perp \boldsymbol {q}}, \end{eqnarray} (39)where the notation $$\left(... \right)_{\perp \boldsymbol {q}}$$ denotes the component perpendicular to $$\boldsymbol {q}$$. To construct the matrix version $$\mathsf {L}$$ of these equations, we must choose a coordinate system. Both because of the constraint (28) and because equation (39) is written in the plane perpendicular to $$\boldsymbol {q}$$ we choose the unit vectors   \begin{eqnarray} \hat{\boldsymbol {a}} \equiv \frac{\hat{\boldsymbol {q}}\times \hat{\boldsymbol {w}}}{\sqrt{1-(\hat{\boldsymbol {q}}\,\cdot\, \hat{\boldsymbol {w}})^2}} \end{eqnarray} (40)and   \begin{eqnarray} \hat{\boldsymbol {b}} &\equiv \hat{\boldsymbol {q}}\times \hat{\boldsymbol {a}}, \end{eqnarray} (41)where $$\hat{w}$$ is any unit vector with $$\hat{\boldsymbol {w}} \cdot \hat{\boldsymbol {q}} \ne 1$$. This choice of basis ensures that our vectors are perpendicular to the wavevector. A choice of particular convenience for $$\hat{w}$$ is   \begin{eqnarray} \hat{\boldsymbol {w}} = \frac{\nabla \Omega }{|\nabla \Omega |}. \end{eqnarray} (42)With this choice $$\hat{a}$$ is time-independent, because the component of $$\boldsymbol {q}$$ perpendicular to $$\boldsymbol {w}$$ is time-independent, and so we may write   \begin{eqnarray} \boldsymbol {\delta \tilde{r}} = \alpha \hat{\boldsymbol {a}} + \beta \hat{\boldsymbol {b}} \end{eqnarray} (43)and   \begin{eqnarray} \boldsymbol {\delta \tilde{v}} = \dot{\alpha } \hat{\boldsymbol {a}} + \dot{\beta } \hat{\boldsymbol {b}} + \beta \partial _t \hat{\boldsymbol {b}}. \end{eqnarray} (44)Note that there is a removable singularity when $$\hat{w} \parallel \hat{q}$$. The matrix $$\mathsf {L}$$ is then given by computing the relation between $$\partial _t \lbrace \alpha , \beta , \dot{\alpha }, \dot{\beta }\rbrace$$ and $$\lbrace \alpha , \beta , \dot{\alpha }, \dot{\beta }\rbrace$$. The result is quite unwieldy so we do not present it here but note that it is fully documented in the software in which we implement these equations. 5 STRESSES AND TRANSPORT The equations of motion contain the position and the velocity, so our expanded vector space is   \begin{eqnarray} \Phi = \left\lbrace \delta \boldsymbol {r}, \delta \boldsymbol {v}, \partial _t \delta \boldsymbol {v},\ldots , \right\rbrace . \end{eqnarray} (45)Combining the linearized equations of motion with our closure scheme, we can compute the correlation function   \begin{eqnarray} \langle \Phi \otimes \Phi \rangle = \int \frac{{\rm d}^3 \boldsymbol {k}}{(2{\pi} )^3} \sum _i \langle \Phi ^{i}_{\boldsymbol {k}}\otimes \Phi ^{i*}_{\boldsymbol {k}}\rangle , \end{eqnarray} (46)where the index i ranges over eigenvectors. This function contains all of the usual stresses and transport functions. For instance, the Reynolds stress is   \begin{eqnarray} R \equiv \langle \delta \boldsymbol {v}\otimes \delta \boldsymbol {v}\rangle = \langle \Phi _1 \otimes \Phi _1 \rangle . \end{eqnarray} (47)Likewise up to a dimensionless constant of order unity the turbulent diffusivity is   \begin{eqnarray} d \equiv \langle \delta \boldsymbol {v}\otimes \delta \boldsymbol {r}\rangle = \langle \Phi _1 \otimes \Phi _0 \rangle . \end{eqnarray} (48)and the turbulent viscosity is   \begin{eqnarray} Q \equiv \langle \delta \boldsymbol {v}\otimes \delta \boldsymbol {r}\rangle +\langle \delta \boldsymbol {r}\otimes \delta \boldsymbol {v}\rangle = \langle \Phi _1 \otimes \Phi _0 \rangle +\langle \Phi _0 \otimes \Phi _1 \rangle . \end{eqnarray} (49)Similar expressions hold for the dynamo effect, the transport of magnetic fields, and material diffusion. 6 RESULTS In this section, we exhibit a number of results which come from applying our model to a wide variety of astronomically and physically relevant circumstances. We also compare with the results of Lesaffre et al. (2013) and Kichatinov & Rudiger (1993). We modify the former to use the convention in equation (9) to avoid spurious differences in scale. We likewise assume that our L0 is equal to three times the mixing length of Kichatinov & Rudiger (1993), as this is an inherent freedom in the formalism and resolves an otherwise-persistent scale difference between our model and theirs. These models have been well-tested against a variety of data, most notably helioseismic results, and so provide a useful reference for our work. We have also included more direct comparisons but, because direct experiments are extremely difficult to perform under most circumstances relevant to astrophysics, we have instead included comparisons with simulations and observations where available and applicable. Simulations are often the most useful comparison for stellar phenomena, because a variety of processes, including meridional circulation, can mask the effects of turbulent transport (Kitchatinov 2013). In accretion discs, however, there are several observable quantities, which are thought to correlate closely with the underlying turbulence and these provide very helpful constraints (King, Pringle & Livio 2007). These comparisons and calculations are not intended to be a complete collection of the results our model can produce, nor have we exhaustively explored the circumstances and dependencies of each result. Rather, it is our hope to demonstrate that there is a great deal of interesting physics in this model, that our perturbative corrections give rise to realistic results and reproduce many known results, and that there is much to warrant further exploration along these lines. 6.1 Rotating convection We begin with the effect of rotation on convection in the case of a rotating system with radial pressure and entropy gradients. It is useful to start by comparing our results with those from simulations. Fig. 1 shows the ratios $$\sqrt{\langle \delta v_r^2\rangle /\langle \delta v^2\rangle }$$, $$\sqrt{\langle \delta v_\theta ^2\rangle /\langle \delta v^2\rangle }$$ and $$\sqrt{\langle \delta v_\phi ^2\rangle /\langle \delta v^2\rangle }$$ for several rotation rates as a function of latitude. The positive latitudes come from table 2 of Chan (2001), while the negative are from table 2 of Käpylä, Korpi & Tuominen (2004). In order to match the units for the rotation rates, we put everything in terms of the coriolis number   \begin{eqnarray} \mathrm{Co} \equiv \frac{\Omega h}{\langle \delta v^2\rangle ^{1/2}}, \end{eqnarray} (50)where, following the convention of Käpylä et al. (2004), 〈δv2〉1/2 was computed for a non-rotating system. Figure 1. View largeDownload slide The ratios $$\sqrt{\langle \delta v_r^2\rangle /\langle \delta v^2\rangle }$$ (blue), $$\sqrt{\langle \delta v_\theta ^2\rangle /\langle \delta v^2\rangle }$$ (red), and $$\sqrt{\langle \delta v_\phi ^2\rangle /\langle \delta v^2\rangle }$$ (purple) are shown for our model (solid) and for simulations by (Käpylä et al. 2004, dots, negative latitude) and (Chan 2001, dots, positive latitude) for a wide range of rotation rates as a function of latitude. The rotation rate is captured by the Coriolis number Co = Ωh/〈δv2〉1/2. Our model general overestimates the anisotropy but captures its variation well. Figure 1. View largeDownload slide The ratios $$\sqrt{\langle \delta v_r^2\rangle /\langle \delta v^2\rangle }$$ (blue), $$\sqrt{\langle \delta v_\theta ^2\rangle /\langle \delta v^2\rangle }$$ (red), and $$\sqrt{\langle \delta v_\phi ^2\rangle /\langle \delta v^2\rangle }$$ (purple) are shown for our model (solid) and for simulations by (Käpylä et al. 2004, dots, negative latitude) and (Chan 2001, dots, positive latitude) for a wide range of rotation rates as a function of latitude. The rotation rate is captured by the Coriolis number Co = Ωh/〈δv2〉1/2. Our model general overestimates the anisotropy but captures its variation well. Our model overestimates the anisotropy of the turbulence but captures its symmetries and trends well. For instance, we find that near the poles and in non-rotating systems the θ and ϕ components of the velocity fluctuations have identical magnitudes, in line with the simulations. We reproduce the trend of decreasing anisotropy towards the equator and decreasing anisotropy with increasing rotation, and, in cases where there are differences between the θ and ϕ velocities, we reproduce both their sign and magnitude. In particular, we find that $$\langle \delta v_r^2\rangle \ge \langle \delta v_\theta ^2 \rangle \ge \langle \delta v_\phi ^2\rangle$$, which is seen in these and other simulations (Rüdiger, Egorov & Ziegler 2005a). Likewise, we find that radial motion makes up a greater fraction of the total velocity near the poles than at the equator, and that as the Coriolis number increases $$\langle \delta v_r^2 - \delta v_\theta ^2 - \delta v_\phi ^2 \rangle \rightarrow 0$$, all of which is in agreement with the predictions of Rüdiger et al. (2005b). Our overestimate of the anisotropy may be due to our model incorporating the large-scale fields on all scales, as noted by Lesaffre et al. (2013). This suggests that a future refinement might be to use estimates of the large-scale modes to compute the environment of those at smaller scales, but we do not treat such complications for now. As a further comparison, we consider the off-diagonal Reynolds stresses of both Chan (2001) and Käpylä et al. (2004). These numbers were extracted from table 3 of the former and also table 3 of the latter and are shown along with our predictions in Fig. 2. In the former, they were straightforward to analyse but in the latter they do not provide a precise test because the simulations included a bulk shear. To correct for this, we used a linear expansion to subtract results across simulations, which were identical in all conditions other than the rotation and thereby determine the effect of the rotation alone. As we will see in Section 6.2 this procedure is problematic because the shear may interact non-linearly with the rotation. Furthermore, because these corrections are of the same order as the terms themselves some care must be taken in interpreting the results. Figure 2. View largeDownload slide The ratios $$\sqrt{\langle \delta v_r \delta v_\theta \rangle /\langle \delta v^2\rangle }$$ (red), $$\sqrt{\langle \delta v_\theta \delta v_\phi \rangle /\langle \delta v^2\rangle }$$ (purple), and $$\sqrt{\langle \delta v_r \delta v_\phi \rangle /\langle \delta v^2\rangle }$$ (blue) are shown from our model (solid) and from simulations by (Käpylä et al. 2004, dots, negative latitude) and (Chan 2001, dots, positive latitude) as a function of latitude. Note that Käpylä et al. (2004) cautions that the moderate rotation simulations had difficulty converging, and these results arise as the difference between two simulations, so it is not clear how significant this test is. Our model generally overestimates these stresses, and suggests a different symmetry for the variation (going as sin θ rather than sin (2θ)). Figure 2. View largeDownload slide The ratios $$\sqrt{\langle \delta v_r \delta v_\theta \rangle /\langle \delta v^2\rangle }$$ (red), $$\sqrt{\langle \delta v_\theta \delta v_\phi \rangle /\langle \delta v^2\rangle }$$ (purple), and $$\sqrt{\langle \delta v_r \delta v_\phi \rangle /\langle \delta v^2\rangle }$$ (blue) are shown from our model (solid) and from simulations by (Käpylä et al. 2004, dots, negative latitude) and (Chan 2001, dots, positive latitude) as a function of latitude. Note that Käpylä et al. (2004) cautions that the moderate rotation simulations had difficulty converging, and these results arise as the difference between two simulations, so it is not clear how significant this test is. Our model generally overestimates these stresses, and suggests a different symmetry for the variation (going as sin θ rather than sin (2θ)). Despite these difficulties some trends are clear and sustained between both sets of data. For instance, in the northern hemisphere (θ > 0), 〈δvrδvθ〉 < 0, while in both hemispheres 〈δvrδvϕ〉 < 0, in keeping with predictions and simulations by Rüdiger et al. (2005b). Likewise, we find that 〈vθvϕ〉 > 0 in the northern hemisphere, in agreement with the findings of Rüdiger et al. (2005a). Once more, however, our model overestimates these anisotropic terms by an amount, which is largely invariant as a function of rotation. This suggests that this overestimate is a systematic offset rather than an error in scaling. We also have some difficulty to reproduce the signs of some of the stresses, particularly in the results of Käpylä et al. (2004), though this could simply be a subtraction difficulty. This is supported by the fact that the simulations themselves do not agree on the signs of these terms and highlights the challenges of making comparisons of terms, which are small in magnitude relative to the scale of the turbulence. To better understand which trends are significant and which are artefacts, we have placed data from comparable rotation rates for the two sets of simulations side-by-side in Fig. 3. The top five panels show the same data as in Fig. 1, while the bottom three show the data from Fig. 2. In general, there is good agreement in the top five panels. The data of Käpylä et al. (2004) gives systematically larger anisotropies and the two sets of simulations occasionally differ on the relative magnitudes of the velocity components (i.e. their ordering), but otherwise the two are in good agreement. By contrast, the bottom three panels paint two very divergent pictures. Neither ordering, trends nor signs are consistent between the two sets of simulations. Only the magnitudes agree in these cases. Thus, the two sets of simulations agree that our model systematically overestimates anisotropies and that, beyond that, our model agrees with them to the extent that they agree with one another. Figure 3. View largeDownload slide The functions shown in Figs 1 and 2 are shown from our model (solid), simulations by (Käpylä et al. 2004, dots, negative latitude) and Chan (2001) (crosses, positive latitude) as a function of latitude. The most comparable pairs of rotation rates were placed side-by-side for each function. A solid black line is shown along the equator where the latitude is zero. There is reasonable agreement on the distribution of velocities in direction but not on the correlations between different velocity directions.. Figure 3. View largeDownload slide The functions shown in Figs 1 and 2 are shown from our model (solid), simulations by (Käpylä et al. 2004, dots, negative latitude) and Chan (2001) (crosses, positive latitude) as a function of latitude. The most comparable pairs of rotation rates were placed side-by-side for each function. A solid black line is shown along the equator where the latitude is zero. There is reasonable agreement on the distribution of velocities in direction but not on the correlations between different velocity directions.. Having compared in detail with these simulations, we now consider predictions which go beyond the domain where simulations are possible. In convection with radial gradients, the leading order effect is to transport heat and material radially. Fig. 4 shows 〈δvrδvr〉 and 〈δvrδrr〉, which are the correlation functions controlling this transport. Figure 4. View largeDownload slide The radial velocity correlation function 〈δvrδvr〉 (red) and the radial diffusivity 〈δvrδrr〉 (blue) are shown in linear scale for Ω < |N| (left-hand panel) and log-log scale for Ω > |N| (right-hand panel). These results are for uniform rotation at a latitude of π/4 with no magnetic field. On this and all subsequent figures $$v_r v_r/L_0^2|N|^2$$ should be read as $$\langle \delta v_r \delta v_r \rangle / L_0^2|N|^2$$ and similarly for other correlations. Shown in purple (*, dashed) for comparison is the result of Kichatinov & Rudiger (1993) with an anisotropy factor of 2, which agrees in sign, scale, and variation. The bumps in our results reflect parameter values where the numerical integration was more difficult. All quantities are given in units of the mixing length and the Brünt–Väisälä frequency. Figure 4. View largeDownload slide The radial velocity correlation function 〈δvrδvr〉 (red) and the radial diffusivity 〈δvrδrr〉 (blue) are shown in linear scale for Ω < |N| (left-hand panel) and log-log scale for Ω > |N| (right-hand panel). These results are for uniform rotation at a latitude of π/4 with no magnetic field. On this and all subsequent figures $$v_r v_r/L_0^2|N|^2$$ should be read as $$\langle \delta v_r \delta v_r \rangle / L_0^2|N|^2$$ and similarly for other correlations. Shown in purple (*, dashed) for comparison is the result of Kichatinov & Rudiger (1993) with an anisotropy factor of 2, which agrees in sign, scale, and variation. The bumps in our results reflect parameter values where the numerical integration was more difficult. All quantities are given in units of the mixing length and the Brünt–Väisälä frequency. Both correlators vary at second order in Ω in the slow rotation limit as expected (Kitchatinov 2013; Lesaffre et al. 2013). In the rapid rotation limit on the other hand, they exhibit clear Ω−1 scaling, consistent with what is seen in other closure models and in simulations (Garaud et al. 2010). The quenching of turbulence in this limit arises because the Coriolis effect acts as a restoring force, stabilizing modes. The peak of each correlator is of order unity and occurs when Ω = 0. In fact, for the stress, the maximum is 0.254647 while for the diffusivity it is 0.28125, both of which are consistent to this precision with Lesaffre et al. (2013), noting that we used the definition in equation (9) for their mixing length. This is because our model is precisely the same as theirs in this limit. Based on this and the observed scalings a good approximation is   \begin{eqnarray} \langle \delta v_r \delta r_r \rangle \approx \langle \delta v_r \delta v_r \rangle \approx \frac{1-(\Omega /|N|)^2}{1 - (\Omega /|N|)^3}. \end{eqnarray} (51) Next, we consider the effect of rotation on the r − θ correlation functions. These functions are responsible for latitudinal transport of heat, mass, and momentum and vanish as a result of spherical symmetry in the non-rotating limit. Fig. 5 shows 〈δvrδvθ〉 and 〈δvrδrθ〉 as a function of the rotation rate. Figure 5. View largeDownload slide The absolute value of the r − θ velocity correlation function 〈δvrδvθ〉 (red) and corresponding diffusivity 〈δvrδrr〉 (blue) are shown in log-log scale against rotation rate. These results are for uniform rotation at a latitude of π/4 with no magnetic field. Shown in purple (*, dashed) for comparison is the result of Kichatinov & Rudiger (1993) with an anisotropy factor of 2, which agrees in sign, scale, and variation. All quantities are given in units of the mixing length and the Brünt–Väisälä frequency. Figure 5. View largeDownload slide The absolute value of the r − θ velocity correlation function 〈δvrδvθ〉 (red) and corresponding diffusivity 〈δvrδrr〉 (blue) are shown in log-log scale against rotation rate. These results are for uniform rotation at a latitude of π/4 with no magnetic field. Shown in purple (*, dashed) for comparison is the result of Kichatinov & Rudiger (1993) with an anisotropy factor of 2, which agrees in sign, scale, and variation. All quantities are given in units of the mixing length and the Brünt–Väisälä frequency. In the slow-rotation regime both quantities scale as Ω2, while in the rapid rotation limit they scale as Ω−1. The peak is of order unity and occurs near Ω = |N|. This gives rise to the approximation   \begin{eqnarray} \langle v_r r_\theta \rangle \approx \langle v_r v_\theta \rangle \approx \frac{(\Omega /|N|)^2}{1 + (\Omega /|N|)^3}. \end{eqnarray} (52)These scalings may be interpreted as a competition between symmetry breaking and quenching: the correlation function rises as rotation breaks symmetries but excessive rotation stabilises the system and quenches the turbulent motions. The symmetry is broken quadratically because, at first order, the Coriolis effect only couples radial and azimuthal motions. The properties of turbulence vary with latitude in a rotating system because the rotation axis picks out a preferred direction. Fig. 6 shows the r − r and r − θ stress and diffusivity correlations as a function of latitude. The r − r correlations vary similarly to one another, exhibiting a minimum at the equator and maxima on-axis. On-axis the rotation drops out of the equations and so the on-axis functions are just those for non-rotating convection. The effect of rotation is then largest at the equator, where the convective motion is predominantly perpendicular to the rotation axis. The correlation functions are smallest where the rotation has the largest effect because rotation primarily acts to stabilize modes. Figure 6. View largeDownload slide Various correlation functions are shown as a function of the angle θ from the rotation axis. The functions are the r − r (left-hand panel) and r − θ (right-hand panel) velocity (red) and diffusivity (blue) correlation functions. These results are for uniform rotation at Ω = 0.2|N| (top panel), Ω = |N| (middle panel), and Ω = 5|N| (bottom panel). Shown in purple (*, dashed) for comparison is the KR result, which agrees in sign and variation but not scale. For slow rotation the scale of the variation is generally smaller than we predict, while for fast rotation the variation is somewhat larger. All quantities are given in units of the mixing length and the Brünt–Väisälä frequency. Figure 6. View largeDownload slide Various correlation functions are shown as a function of the angle θ from the rotation axis. The functions are the r − r (left-hand panel) and r − θ (right-hand panel) velocity (red) and diffusivity (blue) correlation functions. These results are for uniform rotation at Ω = 0.2|N| (top panel), Ω = |N| (middle panel), and Ω = 5|N| (bottom panel). Shown in purple (*, dashed) for comparison is the KR result, which agrees in sign and variation but not scale. For slow rotation the scale of the variation is generally smaller than we predict, while for fast rotation the variation is somewhat larger. All quantities are given in units of the mixing length and the Brünt–Väisälä frequency. By contrast, the r − θ correlator is largest in magnitude at mid-latitudes, vanishing both on-axis and at the equator. On-axis this correlation function must vanish because the $$\hat{\theta }$$ unit vector is ill-defined. The sign change between the northern and southern hemispheres occurs because $$(\hat{\boldsymbol {r}} \times \boldsymbol {\Omega })_\phi$$ has the same sign everywhere while $$(\hat{\boldsymbol {\theta }} \times \boldsymbol {\Omega })_\phi$$ changes sign between the hemispheres. This also explains the vanishing correlation at the equator. The quantities of particular interest for studying the origins of differential rotation are the radial-azimuthal correlation functions 〈δvrδvϕ〉 and 〈δvrδrϕ〉. The former provides a stress coupling the angular momentum to radial motions known as the Λ-effect, while the latter provides a viscosity coupling radial shears to azimuthal motion and so acts as a proxy for the α-effect (Kichatinov & Rudiger 1993). Fig. 7 shows these quantities as a function of the rotation rate. In the slow-rotation limit, both scale as Ω before peaking near unity and falling off as Ω−2 in the rapid-rotation limit. The linear scaling at slow rotation rates is a consequence of the Coriolis effect directly coupling radial and azimuthal motions. These quantities fall off more rapidly than the others in the case of rapid rotation because it is preferentially the modes that couple strongly to the Coriolis effect, which are stabilized the most. The absolute scale of our Λ-effect is approximately what is seen in simulations, slightly overestimating relative to Käpylä et al. (2004) and similar to other theoretical predictions (Kitchatinov 2013; Gough 2012). Figure 7. View largeDownload slide The absolute value of the r − ϕ velocity correlation function 〈δvrδvϕ〉 (red) and corresponding diffusivity 〈δvrδrϕ〉 (blue) are shown in log-log scale versus rotation rate. These results are for uniform rotation at a latitude of π/4 with no magnetic field. Shown in purple (*, dashed) for comparison is the result of Kichatinov & Rudiger (1993) with an anisotropy factor of 2, which agrees in sign, variation, and scale up until Ω = |N|, at which point the behaviour differs significantly. Shown in grey (**, dotted) for comparison is 〈δvrδvϕ〉 from that of Lesaffre et al. (2013). This agrees precisely in the Ω → 0 limit and the agreement is good even near Ω ≈ 0.5|N|. All quantities are given in units of the mixing length and the Brünt–Väisälä frequency. Figure 7. View largeDownload slide The absolute value of the r − ϕ velocity correlation function 〈δvrδvϕ〉 (red) and corresponding diffusivity 〈δvrδrϕ〉 (blue) are shown in log-log scale versus rotation rate. These results are for uniform rotation at a latitude of π/4 with no magnetic field. Shown in purple (*, dashed) for comparison is the result of Kichatinov & Rudiger (1993) with an anisotropy factor of 2, which agrees in sign, variation, and scale up until Ω = |N|, at which point the behaviour differs significantly. Shown in grey (**, dotted) for comparison is 〈δvrδvϕ〉 from that of Lesaffre et al. (2013). This agrees precisely in the Ω → 0 limit and the agreement is good even near Ω ≈ 0.5|N|. All quantities are given in units of the mixing length and the Brünt–Väisälä frequency. 6.2 Differential rotation and convection We now turn to the dependence of convective transport coefficients on differential rotation. We expand our closure model to linear order in the shear and so restrict this analysis to cases where the dimensionless shear |R∇ln Ω| is at most of order unity. Fig. 8 shows the r − θ and r − ϕ velocity and diffusivity correlation functions as a function of differential rotation for a situation where ∇Ω is at an angle of π/4 relative to the pressure gradient. All four functions behave linearly near the origin, with intercept set by the stress and diffusivity in the uniform rotation limit. This is precisely as expected: the intercept is non-zero, giving rise to the Λ-effect, while the slope is non-zero, giving rise to the α-effect (Kichatinov & Rudiger 1993). Note that the favourable comparison of our results with those of Kitchatinov (2013) is helpful because their model was implemented in a two-dimensional solar model, which compared well with helioseismic observations. Figure 8. View largeDownload slide The r − θ (left-hand panel) and r − ϕ (right-hand panel) velocity (red) and diffusivity (blue) correlation functions are shown in log-scale versus the differential rotation. These results are for a convecting region with differential rotation in the cylindrical radial direction, Ω = 0.1|N| and no magnetic field at a latitude of π/4. Shown in purple (*, dashed) for comparison is the result of Kichatinov & Rudiger (1993) with an anisotropy factor of 2. This disagrees in sign and on the magnitude of the slope but agrees in the sign of the slope. Shown in grey (**, dotted) for comparison is 〈δvrδvϕ〉 of Lesaffre et al. (2013). This generally predicts smaller stresses though with the same sign and slope sign as our model. All quantities are given in units of the mixing length and the Brünt–Väisälä frequency. Figure 8. View largeDownload slide The r − θ (left-hand panel) and r − ϕ (right-hand panel) velocity (red) and diffusivity (blue) correlation functions are shown in log-scale versus the differential rotation. These results are for a convecting region with differential rotation in the cylindrical radial direction, Ω = 0.1|N| and no magnetic field at a latitude of π/4. Shown in purple (*, dashed) for comparison is the result of Kichatinov & Rudiger (1993) with an anisotropy factor of 2. This disagrees in sign and on the magnitude of the slope but agrees in the sign of the slope. Shown in grey (**, dotted) for comparison is 〈δvrδvϕ〉 of Lesaffre et al. (2013). This generally predicts smaller stresses though with the same sign and slope sign as our model. All quantities are given in units of the mixing length and the Brünt–Väisälä frequency. A key difference between our work and what we compare with in Fig. 8 is that, while we predict the same sign and comparable magnitude for the α-effect in the zero-shear limit, the effect changes sign near |R∇ln Ω| ≈ 0.5, indicating that, at least for this configuration, this is the point at which non-linear effects become important. This does not represent a particularly severe shear and highlights a key point that the correlation functions we find are generally non-linear in all of the small parameters in which one might wish to expand. Our model captures this non-linear behaviour despite being carried out to linear order in |R∇ln Ω|. This is because, in our expansion, the time evolution operator is what is expanded linearly. The resulting eigenvalues and eigenvectors are generally non-linear functions of this operator. This caution aside, there is a significant regime where the α–Λ expansion is valid and, in this regime, key quantities of interest are the derivatives of the various correlation functions with respect to the shear |R∇Ω|. Fig. 9 shows these derivatives as a function of Ω. The r − ϕ stress derivative is constant in Ω. This means that the stress scales as R∇Ω. This is as expected (see, e.g. equation 79 of Lesaffre et al. 2013) and indicates that there is a well-defined effective viscosity transporting angular momentum. This viscosity is given by   \begin{eqnarray} \nu _{r\phi } \approx L_0^2 |N|. \end{eqnarray} (53) Figure 9. View largeDownload slide The derivatives of various correlation functions with respect to |R∇Ω| are shown as a function of Ω, with both axes log-scaled. The functions are the r − θ (left-hand panel) and r − ϕ (right-hand panel) velocity (red) and diffusivity (blue) correlation functions. These results are for a convecting region with differential rotation in the cylindrical radial direction and no magnetic field at a latitude of π/4. All quantities are given in units of the mixing length and the Brünt–Väisälä frequency. Figure 9. View largeDownload slide The derivatives of various correlation functions with respect to |R∇Ω| are shown as a function of Ω, with both axes log-scaled. The functions are the r − θ (left-hand panel) and r − ϕ (right-hand panel) velocity (red) and diffusivity (blue) correlation functions. These results are for a convecting region with differential rotation in the cylindrical radial direction and no magnetic field at a latitude of π/4. All quantities are given in units of the mixing length and the Brünt–Väisälä frequency. By contrast, the derivatives of the r − θ correlations as well as the r − ϕ diffusivity all diverge in the limit as Ω → 0. In particular, the r − θ correlations diverge as Ω−1, while the r − ϕ diffusivity diverges as Ω−2. These divergences are signatures of symmetry breaking. They indicate that the direction in which the R∇Ω → 0 limit is approached matters. That is, this limit can be approached by first letting Ω → 0 and then differentiating or by differentiating and then taking Ω → 0 and the divergence we find in the latter approach indicates that the order matters. When Ω = 0 and |R∇Ω| = 0, there is a symmetry between ±θ and between ±ϕ. As a result both the r − θ and r − ϕ terms vanish in this limit. When Ω ≠ 0 these symmetries are broken by the rotation and we know from Figs 5 and 7 that this occurs at first order for r − ϕ and second order for r − θ. In the opposing limit, the situation is different because in the time evolution described by equation (39) $$\mathsf {L}$$ is independent of |R∇Ω| when Ω = 0. There is, however, a dependence on |R∇Ω| through the time-dependence of $$\boldsymbol {q}$$. This breaks the ϕ symmetry because $$\partial _t \boldsymbol {q}$$ is proportional to qϕR∇Ω and hence is sensitive to ϕ. It does not, however, break the θ symmetry, because qϕR∇Ω is symmetric with respect to changing the signs of both θ and $$\boldsymbol {q}$$. It follows then that we should find divergences in the r − θ correlation derivatives owing to the path-dependence of the zero-rotation limit and that we should find the r − ϕ derivatives to be generally well-behaved. The curious divergence is then that in the r − ϕ diffusivity, because this correlation function does not suffer from a symmetry-derived path-dependence. This arises because the differential rotation means that $$\mathsf {L}$$ is time-dependent. This introduces polynomial corrections to the usual exponential growth, as discussed in Section 3. This formalism captures the fact that the differential rotation turns vertical displacement into ϕ displacements, which vary as polynomials in time. There are, therefore, modes with very small radial velocities, which nevertheless have large azimuthal displacements and these dominate the diffusivity derivative. These modes grow proportional to |R∇Ω| and their growth may proceed in the azimuthal direction until bounded by the Coriolis effect at a time Ω−1. As a result these modes contribute to the diffusivity as |R∇ln Ω| and hence lead to a diverging derivative in |R∇Ω| as Ω → 0. 6.3 Differential rotation and stable stratification Stably stratified regions are those with   \begin{eqnarray} N^2 > 0, \end{eqnarray} (54)such that buoyancy acts to counter perturbations in the vertical direction. This tends to damp turbulence. In the presence of such damping, there can still be turbulence if there is also a shear. The classic example of this is the Kelvin–Helmholtz phenomenon, which can occur in such a system if the Richardson criterion   \begin{eqnarray} \frac{|{\rm d}u/{\rm d}z|^2}{|N|^2} > \frac{1}{4} \end{eqnarray} (55)is satisfied (Zahn 1993). Here, u is the velocity and z is the coordinate parallel to the stratification. Even when this criterion is not satisfied, latitudinal shear can still generate turbulence (Canuto et al. 2008). These motions are suppressed in vertical extent by the stratification and hence are primarily confined to the plane perpendicular to the stratification direction. Fig. 10 shows the dependence on shear strength of all six stress components in a rotating stably stratified zone with latitudinal rotational shear. All six exhibit linear scaling with the shear strength. This is unusual in an otherwise-stable zone because it implies a viscosity which, to leading order, does not depend on the shear. That is,   \begin{eqnarray} \nu _{ij} \approx L_0^2 N f_{ij}\left(\frac{\Omega }{|N|}\right), \end{eqnarray} (56)where fij is some function of the angular velocity. Fig. 11 shows the dependence of the stress components on Ω/|N| for fixed |R∇ln Ω| = 0.1. The r − θ and θ − ϕ stresses vary as Ω3 in the slow-rotation regime and as Ω2 for rapid rotation. The other components all scale as Ω2 in both regimes. Thus, for instance, frϕ = Ω/|N| because the viscosity is the derivative of the stress with respect to the shear, and hence   \begin{eqnarray} \nu _{r\phi } \approx 10^{-5} L_0^2 \Omega . \end{eqnarray} (57) Figure 10. View largeDownload slide The absolute value of the r − r (red), r − θ (blue) and r − ϕ (purple) velocity correlation functions are shown as a function of |R∇ln Ω|, with both axes log-scaled. These results are for a stably stratified region with differential rotation in the radial direction, Ω = 0.1|N| and no magnetic field. The data is computed for a point on the equator with radial differential rotation. All quantities are given in units of the mixing length and the Brünt–Väisälä frequency. Figure 10. View largeDownload slide The absolute value of the r − r (red), r − θ (blue) and r − ϕ (purple) velocity correlation functions are shown as a function of |R∇ln Ω|, with both axes log-scaled. These results are for a stably stratified region with differential rotation in the radial direction, Ω = 0.1|N| and no magnetic field. The data is computed for a point on the equator with radial differential rotation. All quantities are given in units of the mixing length and the Brünt–Väisälä frequency. Figure 11. View largeDownload slide The absolute value of the r − r (red), r − θ (blue), and r − ϕ (purple) velocity correlation functions are shown as a function of Ω/|N| for fixed |R∇ln Ω| = 0.1, with both axes log-scaled. These results are for a stably stratified region with differential rotation in the radial direction and no magnetic field. The data is computed for a point on the equator with radial differential rotation. All quantities are given in units of the mixing length and the Brünt–Väisälä frequency. Figure 11. View largeDownload slide The absolute value of the r − r (red), r − θ (blue), and r − ϕ (purple) velocity correlation functions are shown as a function of Ω/|N| for fixed |R∇ln Ω| = 0.1, with both axes log-scaled. These results are for a stably stratified region with differential rotation in the radial direction and no magnetic field. The data is computed for a point on the equator with radial differential rotation. All quantities are given in units of the mixing length and the Brünt–Väisälä frequency. The scaling in equation (57) arises owing to the centrifugal term, which has a destabilising effect when Ω increases with $$\hat{R}$$. When |R∇Ω| = 0 this effect is not present so the system is stable but introducing a small differential rotation produces an acceleration proportional to $$\Omega R \boldsymbol {\delta r}\cdot \nabla \Omega$$ and hence   \begin{eqnarray} \partial _t^2 \boldsymbol {\delta r} \approx g^2 \boldsymbol {\delta r} \propto \hat{R}\Omega R \boldsymbol {\delta r}\cdot \nabla \Omega , \end{eqnarray} (58)which means that the stress scales as Ω∇Ω and thence the viscosity scales as Ω. Because the magnitudes of the stresses are always ordered in the same way, the same terms are always the most significant. From largest to smallest, the stresses are r − r, r − ϕ, ϕ − ϕ, θ − θ, θ − ϕ, and r − θ. This group is nearly separated into diagonal stresses, which are larger, and off-diagonal stresses, which are smaller. The exception to this rule is the r − ϕ stress, which is special because it is the term which directly couples to the shear. The ordering of the remaining terms is not surprising because the off-diagonal stresses are typically mediated by a coupling between different directions, whereas the on-diagonal stresses require no such coupling. To better understand the effect of our perturbative corrections, we computed the same results without them. This produced stresses, which were zero to within numerical precision in all cases, indicating that the entire contribution in this case is coming from the perturbation. However with a different angle of differential rotation, we obtained non-zero results. It is instructive then to compare Fig. 12 with Fig. 13. These show the same correlation functions as each other in the same physical scenario, with differential rotation this time at an angle of π/4, but the former uses the first order perturbative expansion while the latter only expands to zeroth order. The difference between the two calculations is striking: many of the correlation functions have fundamentally different scalings when the perturbative corrections are taken into account. In particular, the r − θ and r − r stresses are both quadratic in the shear and the r − ϕ and θ − ϕ stresses both vary as the shear to the 3/2 power, whereas they are all linear in the shear in the expanded calculation. This difference relates in part to the centrifugal term, which couples the displacement to the acceleration. Without expanding the equations of motion, we would have δr ∝ δv, because the mode would need to be an eigenvector of $$\mathsf {M}$$. The modes which couple to the centrifugal term would still grow according to equation (58) but, for most modes, arranging for the displacement to couple to this term requires coupling to the stabilizing buoyant term too. To make this clearer, in Fig. 14, we have computed the growth rate as a function of wave-vector orientation without using the perturbative expansion. There are several rapidly growing regions, oriented at angles of ± π/4 relative to the vertical. These angles represent a compromise between maximizing the magnitude of the centrifugal acceleration and maximizing its projection on to the velocity, both subject to the Boussinesq condition that motion be in the plane perpendicular to $$\boldsymbol {q}$$. Figure 12. View largeDownload slide The absolute value of the r − r (red), r − θ (blue), and r − ϕ (purple) velocity correlation functions are shown as a function of |R∇ln Ω|, with both axes log-scaled. The correlation functions are evaluated at first order in the perturbative expansion rather than first order. These results are for a stably stratified region with differential rotation in the radial direction, Ω = 0.1|N| and no magnetic field. The data are computed for a point on the equator with differential rotation at an angle of π/4. All quantities are given in units of the mixing length and the Brünt–Väisälä frequency. Figure 12. View largeDownload slide The absolute value of the r − r (red), r − θ (blue), and r − ϕ (purple) velocity correlation functions are shown as a function of |R∇ln Ω|, with both axes log-scaled. The correlation functions are evaluated at first order in the perturbative expansion rather than first order. These results are for a stably stratified region with differential rotation in the radial direction, Ω = 0.1|N| and no magnetic field. The data are computed for a point on the equator with differential rotation at an angle of π/4. All quantities are given in units of the mixing length and the Brünt–Väisälä frequency. Figure 13. View largeDownload slide The absolute value of the r − r (red), r − θ (blue) and r − ϕ (purple) velocity correlation functions are shown as a function of |R∇ln Ω|, with both axes log-scaled. The correlation functions are evaluated at zeroth order in the perturbative expansion rather than first order. These results are for a stably stratified region with differential rotation in the radial direction, Ω = 0.1|N| and no magnetic field. The data are computed for a point on the equator with differential rotation at an angle of π/4. All quantities are given in units of the mixing length and the Brünt–Väisälä frequency. Figure 13. View largeDownload slide The absolute value of the r − r (red), r − θ (blue) and r − ϕ (purple) velocity correlation functions are shown as a function of |R∇ln Ω|, with both axes log-scaled. The correlation functions are evaluated at zeroth order in the perturbative expansion rather than first order. These results are for a stably stratified region with differential rotation in the radial direction, Ω = 0.1|N| and no magnetic field. The data are computed for a point on the equator with differential rotation at an angle of π/4. All quantities are given in units of the mixing length and the Brünt–Väisälä frequency. Figure 14. View largeDownload slide The square of the growth rate is shown as a function of wave-vector orientation on a logarithmic colour scale. The wave-vector is specified by a magnitude and two angles, θ(q) and ϕ(q), which are spherical angles relative to the $$\hat{z}$$ direction. These rates were computed with a zeroth-order expansion. Regions with squared growth rates below 10−16 are shown in white. All quantities are given in units of the mixing length and the Brünt–Väisälä frequency. Figure 14. View largeDownload slide The square of the growth rate is shown as a function of wave-vector orientation on a logarithmic colour scale. The wave-vector is specified by a magnitude and two angles, θ(q) and ϕ(q), which are spherical angles relative to the $$\hat{z}$$ direction. These rates were computed with a zeroth-order expansion. Regions with squared growth rates below 10−16 are shown in white. All quantities are given in units of the mixing length and the Brünt–Väisälä frequency. By contrast, the growth rates in the expanded system, shown in Fig. 15, are significant over a much wider swath of parameter space. This is because, in the expanded system, the displacement and velocity need not be parallel so the displacement can be chosen to maximize the centrifugal term while the velocity can be chosen to maximize the projection of the acceleration on to the velocity. Figure 15. View largeDownload slide The square of the growth rate is shown as a function of wave-vector orientation on a logarithmic colour scale. The wave-vector is specified by a magnitude and two angles, θ(q) and ϕ(q), which are spherical angles relative to the $$\hat{z}$$ direction. These rates were computed with a first-order expansion. Regions with squared growth rates below 10−16 are shown in white. All quantities are given in units of the mixing length and the Brünt–Väisälä frequency. Figure 15. View largeDownload slide The square of the growth rate is shown as a function of wave-vector orientation on a logarithmic colour scale. The wave-vector is specified by a magnitude and two angles, θ(q) and ϕ(q), which are spherical angles relative to the $$\hat{z}$$ direction. These rates were computed with a first-order expansion. Regions with squared growth rates below 10−16 are shown in white. All quantities are given in units of the mixing length and the Brünt–Väisälä frequency. 6.4 Baroclinic instability The baroclinic instability arises in otherwise stably stratified fluids when the entropy gradient is not parallel to the pressure gradient (Killworth 1980). In fact, this is part of a family of instabilities, which includes the convective instability (Lebovitz 1965). This family provides a continuous connection between the unstable convective and stably stratified limits. To explore, it consider Fig. 16 which shows the variation of r − r and r − θ correlation functions against the angle δ between the entropy gradient and the pressure gradient. The radial correlations peak when the two gradients are aligned. This is the convective limit. These correlations fall to zero in the opposing limit where the two gradients are anti-aligned, which is the stably stratified limit. In between these limits, the behaviour is approximately that of cos 2δ. Figure 16. View largeDownload slide Various correlation functions are shown as a function of the angle δ between the entropy and pressure gradients. The functions are the r − r (left-hand panel) and r − θ (right-hand panel) velocity (red) and diffusivity (blue) correlation functions. These results are for a non-rotating convective region on the equator and with no magnetic field. Shown in purple (*, dashed) for comparison is the KR result. This agrees in sign, and for r − r agrees in scale, but their r − θ prediction is considerably larger. Notably, this comparison is precisely as cos (δ) (left-hand panel) and sin (δ) (right-hand panel) and crosses zero at non-extremal angles. This is most likely because their theory is not designed for nearly stable regions with extreme baroclinicity. All quantities are given in units of the mixing length and the Brünt–Väisälä frequency. Figure 16. View largeDownload slide Various correlation functions are shown as a function of the angle δ between the entropy and pressure gradients. The functions are the r − r (left-hand panel) and r − θ (right-hand panel) velocity (red) and diffusivity (blue) correlation functions. These results are for a non-rotating convective region on the equator and with no magnetic field. Shown in purple (*, dashed) for comparison is the KR result. This agrees in sign, and for r − r agrees in scale, but their r − θ prediction is considerably larger. Notably, this comparison is precisely as cos (δ) (left-hand panel) and sin (δ) (right-hand panel) and crosses zero at non-extremal angles. This is most likely because their theory is not designed for nearly stable regions with extreme baroclinicity. All quantities are given in units of the mixing length and the Brünt–Väisälä frequency. By contrast, the r − θ correlations behave approximately as sin δ, and vanishes when δ = 0. This is because both the aligned and the anti-aligned limits are spherically symmetric and so must have this correlation function vanish. Deviations from the convective limit give rise to linear scaling so the convective baroclinic instability transports heat and momentum at first order in the baroclinicity. This is an entirely distinct phenomenon from the thermal wind balance, which is a large-scale effect while this results from integrating out the small-scale turbulent modes. In the stable limit perturbations arise quadratically, a deviation from the behaviour of sin δ. This is because there are no existing turbulent motions to perturb, and so each position and velocity component is linear in δ and gives rise to a quadratic two-point correlation function. 6.5 Stellar magnetism We now turn to the impact of the magnetic field on convective turbulence in stars. Fig. 17 shows 〈δvrδvr〉 in a mildly rotating (Ω = 0.1|N|) convection zone as a function of B for three polarisations; radial ($$\boldsymbol {B} \parallel \hat{\boldsymbol {r}}$$), latitudinal ($$\boldsymbol {B} \parallel \hat{\boldsymbol {\theta }}$$), and longitudinal ($$\boldsymbol {B} \parallel \hat{\boldsymbol {\phi }}$$). As the field increases, the stress falls off. This is because the field quenches the turbulence by providing a stabilizing restoring force, and is in general agreement with Canuto & Hartke (1986). Interestingly, the only significant differences are between the radial and angular field polarisations! The θ and ϕ polarisations show precisely the same behaviour out to very strong fields. This is a result of symmetry, because the radial stress is not sensitive to rotation about the radial direction. The deviation seen with strong fields is a numeric artefact and decreases with increasing integration time. Figure 17. View largeDownload slide The stress 〈δvrδvr〉 is shown as a function of magnetic field strength. The magnetic field is polarized radially (red), longitudinally (purple), and latitudinally (blue). The system is rigidly rotating at Ω = 0.1|N| at a latitude of π/4. All quantities are given in units of the mixing length and Brünt–Väisälä frequency. Figure 17. View largeDownload slide The stress 〈δvrδvr〉 is shown as a function of magnetic field strength. The magnetic field is polarized radially (red), longitudinally (purple), and latitudinally (blue). The system is rigidly rotating at Ω = 0.1|N| at a latitude of π/4. All quantities are given in units of the mixing length and Brünt–Väisälä frequency. By contrast, consider 〈δvrδvϕ〉, shown in Fig. 18. This component, along with the corresponding Maxwell stress, is responsible for transporting angular momentum. Interestingly, it shows differences amongst all polarisations, with the strongest difference between the θ polarization and the others. This is because the stress is mixed between different directions and so is sensitive to all variations in the magnetic field direction. The large difference of the θ polarization relative to the others reflects the fact that motion is damped perpendicular to the magnetic field so the θ polarization damps motion in both directions involved in this component of the stress, whereas the r and ϕ polarizations only damp motion in one of the two directions. Figure 18. View largeDownload slide The stress 〈δvrδvϕ〉 is shown as a function of magnetic field strength. The magnetic field is polarized radially (red), longitudinally (purple) and latitudinally (blue). The system is rigidly rotating at Ω = 0.1|N| at a latitude of π/4. All quantities are given in units of the mixing length and Brünt–Väisälä frequency. Figure 18. View largeDownload slide The stress 〈δvrδvϕ〉 is shown as a function of magnetic field strength. The magnetic field is polarized radially (red), longitudinally (purple) and latitudinally (blue). The system is rigidly rotating at Ω = 0.1|N| at a latitude of π/4. All quantities are given in units of the mixing length and Brünt–Väisälä frequency. 6.6 Magnetorotational instability As a final example, we consider the MRI (Chandrasekhar 1960). This instability arises in magnetized fluids undergoing Keplerian orbital motion. Fig. 19 shows the r − r, θ − ϕ, r − ϕ, and r − θ Reynolds and Maxwell stresses for an accretion disc with a vertical magnetic field. Contrary to predictions (Chandrasekhar 1960) none of the Reynolds stresses vanish in the zero-field limit. This is because the linear system supports short-term growing modes but, while they only grow in the short-time limit, our numerical methods are not sensitive to that effect at this order. In principle, at higher order, this phenomenon should become evident and so this may be interpreted as an artefact associated with our expanding to low order in |R∇Ω| > 1. Despite this, it is likely that other non-magnetic processes can destabilise these modes even in the long term and so we feel it is appropriate to at least consider them (cf. Luschgy & Pagès 2006). The Maxwell stresses by contrast do vanish as B → 0. This is to be expected because they are proportional to B2. Figure 19. View largeDownload slide From the top to bottom panel are the r − ϕ, r − r, r − θ, and θ − ϕ stresses. The Reynolds (velocity) stresses are in red and the Maxwell (magnetic) stresses are in blue. Note that it is the negative r − r Maxwell stress, which is shown to make the comparison with the Reynolds stresses clearer. In all cases, they are shown as a function of B for a Keplerian disc. The magnetic field is taken parallel to $$\hat{z}$$. The system is taken to be stably stratified in the vertical direction with |N| = Ω and hence L0 = h = R. All quantities are given in units of the mixing length and Ω. Figure 19. View largeDownload slide From the top to bottom panel are the r − ϕ, r − r, r − θ, and θ − ϕ stresses. The Reynolds (velocity) stresses are in red and the Maxwell (magnetic) stresses are in blue. Note that it is the negative r − r Maxwell stress, which is shown to make the comparison with the Reynolds stresses clearer. In all cases, they are shown as a function of B for a Keplerian disc. The magnetic field is taken parallel to $$\hat{z}$$. The system is taken to be stably stratified in the vertical direction with |N| = Ω and hence L0 = h = R. All quantities are given in units of the mixing length and Ω. As the magnetic field increases, the r − ϕ and θ − ϕ Reynolds stresses change sign. This indicates the onset of MRI modes, which have the opposite sign to the zero-field correlations. This effect saturates when vA ≈ Ωh, where h is the scale height of the disc. The total r − ϕ stress saturates at roughly 10−2(hΩ)2, which lies between those typically found in simulations and those inferred from observations (Starling et al. 2004; King et al. 2007). Note that at the saturation point the Maxwell and Reynolds stresses are comparable, and beyond this point the Maxwell stress increases while the Reynolds stress falls off. Above the saturation point the Reynolds stresses drop off as the magnetic field quenches the turbulence. This is precisely what is expected for the MRI (Balbus & Hawley 1991). The Maxwell stresses, however, continue to grow, again in line with expectations. Some care is required to interpret these results because they were computed for a fixed field and that field may or may not be stable under the action of the turbulence it generates (Pessah, Chan & Psaltis 2006). Furthermore, there are challenges with the α-disc prescription, which make the specific stress components more difficult to interpret (Pessah, Chan & Psaltis 2008). Nevertheless, it is encouraging that what we see matches well with both observations and simulations. 7 CONCLUSION We have derived a turbulent closure model, which incorporates shear, rotation, and magnetism as well as a full three-dimensional spectrum of fluctuations. We have also presented a new perturbative approach to incorporate time-dependence in the evolution equations. This model, which is implemented in an open source numerical software package, fully reproduces many known phenomena such as the MRI, baroclinic instability, rotational quenching, and more classic shear instabilities. Using this model, we have determined the asymptotic behaviour of a wide variety of correlation functions and transport coefficients under a wide range of circumstances, many of which do not appear in the literature. We have further explored the behaviour of turbulent transport coefficients in intermediate regimes where no single phenomenon dominates, such as in the critical MRI. In these cases, the behaviour is generally complex and does not separate easily into components associated with the different pieces of input physics. The closure formalism developed here fills a new niche in the landscape of solutions to turbulent transport, covering enough phenomena to be useful to understand those operating in stars, planets, and accretion discs, while being rapid enough to be incorporated into stellar evolution codes on nuclear time-scales. In the future, we hope to provide further refinements and comparisons with direct numerical simulations as well as experiments. In addition, it would be interesting to explore the results of this model to higher order in the shear and, even at this order, there are many results which deserve more analysis than we have given here. Acknowledgements ASJ acknowledges financial support from a Marshall Scholarship as well as support from the Institute of Astronomy, École Normale Supérieure (ENS), and Centre for Excellence in Basic Sciences (CEBS) to work at ENS Paris and CEBS in Mumbai. PL acknowledges travel support from the french PNPS (Programme National de Physique Stellaire) and from CEBS. CAT thanks Churchill College for his fellowship. SMC is grateful to the IOA for support and hopsitality and thanks the Cambridge-Hamied exchange program for financial support. The authors also thank Rob Izzard and Science and Technology Facilities Council Grant ST/L003910/1 for CPU cycles, which aided in this work. Footnotes 1 It would not be difficult, however, to incorporate them into this framework at a later date. REFERENCES Ashkenazi S., Steinberg V., 1999, Phys. Rev. Lett. , 83, 4760 CrossRef Search ADS   Balbus S. A., Hawley J. F., 1991, ApJ , 376, 214 CrossRef Search ADS   Balbus S. A., Schaan E., 2012, MNRAS , 426, 1546 CrossRef Search ADS   Benzi R., Tripiccione R., Massaioli F., Succi S., Ciliberto S., 1994, Europhys. Lett. , 25, 341 CrossRef Search ADS   Berntsen J., Espelid T. O., Genz A., 1991, ACM Trans. Math. Softw. , 17, 437 CrossRef Search ADS   Böhm-Vitense E., 1958, ZAp , 46, 108 Canuto V. M., 1994, ApJ , 428, 729 CrossRef Search ADS   Canuto V. M., 1997, ApJ , 482, 827 CrossRef Search ADS   Canuto V. M., Hartke G. J., 1986, A&A , 168, 89 Canuto V. M., Cheng Y., Howard A. M., Esau I. N., 2008, J. Atmos. Sci. , 65, 2437 CrossRef Search ADS   Carati D., 1990, Phys. Rev. A , 41, 3129 CrossRef Search ADS PubMed  Chan K. L., 2001, ApJ , 548, 1102 CrossRef Search ADS   Chandrasekhar S., 1960, Proc. Natl. Acad. Sci. , 46, 253 CrossRef Search ADS   Dobrowolny M., Mangeney A., Veltri P., 1980, Phys. Rev. Lett. , 45, 144 CrossRef Search ADS   Floquet G., 1883, Annales scientifiques de l’École normale supérieure , 12, 47 CrossRef Search ADS   Garaud P., Ogilvie G. I., Miller N., Stellmach S., 2010, MNRAS , 407, 2451 CrossRef Search ADS   Garaud P., Gagnier D., Verhoeven J., 2017, ApJ , 837, 133 CrossRef Search ADS   Genz A., Malik A., 1980, J. Comput. Appl. Math. , 6, 295 CrossRef Search ADS   Goldreich P., Sridhar S., 1995, ApJ , 438, 763 CrossRef Search ADS   Gough D., 1977, in Spiegel E. A., Zahn J.-P., eds, Astrophysics and Space Science Library, Vol. 71, Problems of Stellar Convection . Springer-Verlag, Berlin p. 15 Gough D. O., 2012, A&A , 2012, 987275 Guennebaud G. et al.  , 2010. Available at http://eigen.tuxfamily.org Hunter J. D., 2007, Comput. Sci. Eng. , 9, 90 CrossRef Search ADS   Isserlis L., 1918, Biometrika , 12, 134 CrossRef Search ADS   Käpylä P. J., Korpi M. J., Tuominen I., 2004, A&A , 422, 793 CrossRef Search ADS   Kichatinov L. L., 1986, Geophys. Astrophys. Fluid Dyn. , 35, 93 CrossRef Search ADS   Kichatinov L. L., 1987, Geophys. Astrophys. Fluid Dyn. , 38, 273 CrossRef Search ADS   Kichatinov L. L., Rudiger G., 1993, A&A , 276, 96 Killworth P. D., 1980, Dyn. Atmos. Oceans , 4, 143 CrossRef Search ADS   King A. R., Pringle J. E., Livio M., 2007, MNRAS , 376, 1740 CrossRef Search ADS   Kitchatinov L. L., 2013, in Kosovichev A. G., de Gouveia Dal Pino E., Yan Y., eds, Proc. IAU Symp. 294, Solar and Astrophysical Dynamos and Magnetic Activity . Kluwer, Dordrecht, p. 399 Kolmogorov A., 1941a, Akademiia Nauk SSSR Doklady , 30, 301 Kolmogorov A. N., 1941b, Akademiia Nauk SSSR Doklady , 32, 16 Launder B., Spalding D., 1974, Comput. Methods Appl. Mech. Eng. , 3, 269 CrossRef Search ADS   Lebovitz N. R., 1965, ApJ , 142, 1257 CrossRef Search ADS   Lee D., 2013, J. Comput. Phys. , 243, 269 CrossRef Search ADS   Lesaffre P., Balbus S. A., Latter H., 2009, MNRAS , 396, 779 CrossRef Search ADS   Lesaffre P., Chitre S. M., Potter A. T., Tout C. A., 2013, MNRAS , 431, 2200 CrossRef Search ADS   Lohse D., Xia K.-Q., 2010, Annu. Rev. Fluid Mech. , 42, 335 CrossRef Search ADS   Luschgy H., Pagès G., 2006, preprint (arXiv:math/0607282) Maeder A., 1995, A&A , 299, 84 McKinney J. C., Tchekhovskoy A., Sadowski A., Narayan R., 2014, MNRAS , 441, 3177 CrossRef Search ADS   Murphy G. C., Pessah M. E., 2015, ApJ , 802, 139 CrossRef Search ADS   Pessah M. E., Goodman J., 2009, ApJ , 698, L72 CrossRef Search ADS   Pessah M. E., Chan C.-K., Psaltis D., 2006, MNRAS , 372, 183 CrossRef Search ADS   Pessah M. E., Chan C.-K., Psaltis D., 2008, MNRAS , 383, 683 CrossRef Search ADS   Procaccia I., Zeitak R., 1989, Phys. Rev. Lett. , 62, 2128 CrossRef Search ADS PubMed  Rüdiger G., Egorov P., Ziegler U., 2005a, Astron. Nachr. , 326, 315 CrossRef Search ADS   Rüdiger G., Egorov P., Kitchatinov L. L., Küker M., 2005b, A&A , 431, 345 CrossRef Search ADS   Salvesen G., Simon J. B., Armitage P. J., Begelman M. C., 2016, MNRAS , 457, 857 CrossRef Search ADS   Schou J. et al.  , 1998, ApJ , 505, 390 CrossRef Search ADS   Smolec R., Houdek G., Gough D., 2011, in Brummell N. H., Brun A. S., Miesch M. S., Ponty Y., eds, Proc. IAU Symp. 271, Astrophysical Dynamics: From Stars to Galaxies . Kluwer, Dordrecht, p. 397 Sridhar S., Goldreich P., 1994, ApJ , 432, 612 CrossRef Search ADS   Starling R. L. C., Siemiginowska A., Uttley P., Soria R., 2004, MNRAS , 347, 67 CrossRef Search ADS   Sukoriansky S., Dikovskaya N., Galperin B., 2007, J. Atmos. Sci. , 64, 3312 CrossRef Search ADS   van der Walt S., Colbert S. C., Varoquaux G., 2011, Comput. Sci. Eng. , 13, 22 CrossRef Search ADS   Wick G. C., 1950, Phys. Rev. , 80, 268 CrossRef Search ADS   Yakhot V., Orszag S. A., 1986, J. Sci. Comput. , 1, 3 CrossRef Search ADS   Zahn J.-P., 1993, Space Sci. Rev. , 66, 285 CrossRef Search ADS   Zhou Y., McComb D. W., Vahala G., 1997, NASA Contractor Report 201718  APPENDIX A: TURBULENT INDEX The general question of which turbulent index to use and under what circumstances remains open though many specific cases are well understood. In the case of isotropic incompressible turbulence the Kolmogorov index is well-known to be n = 11/6 (Kolmogorov 1941a). There is more debate over the index to use for convection, with answers ranging from n = 5/2 (Benzi et al. 1994) to n = 21/10 (Procaccia & Zeitak 1989) and n = 2.4 ± 0.2 (Ashkenazi & Steinberg 1999). There has also been work attempting to determine the spectrum in a context-sensitive manner through energy balance arguments (Yakhot & Orszag 1986). In the magnetised case sources differ even more, with some suggesting that this range still applies (Dobrowolny, Mangeney & Veltri 1980), some arguing for a Kolmogorov-like spectrum (Goldreich & Sridhar 1995) and others giving a range of indices depending on geometry and the direction of the wavevector (Sridhar & Goldreich 1994). From numerical experiments with our closure model, we have found that the magnetic stress scales sufficiently rapidly with k that it is divergent for n = 11/6 and not for n = 8/3. This favours the scenario of Goldreich & Sridhar (1995), who argue that in the strongly magnetized limit the index ought to be n = 8/3. In order to consistently treat both the non-magnetic and the strongly magnetized limits, we choose a simple prescription in which n = 11/6 when one of |N|, or |R∇Ω| exceeds kvA and use n = 8/3 otherwise. This means that there is a critical wavenumber   \begin{eqnarray} k_{\rm c} \equiv \frac{\max \left(|N|, |R\nabla \Omega |\right)}{v_A} \end{eqnarray} (A1)at which the spectrum changes. In the non-magnetic case, the evolution matrix is independent of the magnitude of the wavevector and so altering the index just alters the correlation coefficients by a multiplicative factor. In the magnetic case, the potential for error is larger because the magnitude of the wavevector is relevant but there appears to be no consensus on the best prescription and so we make do with what is available. APPENDIX B: BOUSSINESQ ODDITIES In this work, we have taken the Boussinesq approximation. In Fourier space this is   \begin{eqnarray} \boldsymbol {q}\cdot \tilde{\delta \boldsymbol {r}}=0. \end{eqnarray} (B1)Taking the time derivative of both sides we see that   \begin{eqnarray} \partial _t\left(\boldsymbol {q}\cdot \delta \tilde{\boldsymbol {r}}\right) = \boldsymbol {q}\cdot \delta \tilde{\boldsymbol {v}} + \delta \tilde{\boldsymbol {r}}\cdot \partial _t\boldsymbol {q} = 0. \end{eqnarray} (B2)As a result   \begin{eqnarray} \delta \tilde{\boldsymbol {v}}\cdot \boldsymbol {q} = - \delta \tilde{\boldsymbol {r}}\cdot \partial _t\boldsymbol {q} \ne 0. \end{eqnarray} (B3)This is quite peculiar, but is just an artefact of our coordinate system. Because the wavevectors are time-dependent, maintaining the volume of a fluid parcel requires that the displacement be orthogonal to the wavevector, which actually means that the velocity is generally not orthogonal to the wavevector. APPENDIX C: SOFTWARE DETAILS The software used for this work is Mixer version 1, which we have released under a GPLv3 license at github.com/adamjermyn/Mixer. All data produced for this work are available at the same location as HDF5 tables with attributes documenting the physical inputs. Post-processing and visualization of the data was with the Python modules Numpy (van der Walt, Colbert & Varoquaux 2011) and Matplotlib (Hunter 2007) and the relevant scripts for this are included with Mixer. The core of Mixer is written in C++, for performance reasons, and the code is supplied with a Makefile, which supports compilation on both Linux and MacOS. Mixer makes use of the Eigen library (Guennebaud et al. 2010) for linear algebra. Mixer also uses the Cubature library for numerical integration. This library is an implementation of the algorithms by Genz & Malik (1980) and Berntsen, Espelid & Genz (1991). These integration routines are supplemented by a Python integration routine tailored for integrands with small support regions. The details will be explored in later work. In addition, many routines provide a Python interface. Currently Mixer supports only single-threaded operation, though it may be used inside parallelized scripts through the Python wrapper. The version of Mixer used to generate the data in this work was compiled against Cubature version 1.0.2 and Eigen version 3.3.3, though the code does not use any features which require recent versions, so many likely suffice. Mixer is optimized for convecting systems for which achieving accuracy better than 10−5 relative and absolute typically requires between 1 ms and 1 s on a single core of a 2016 Intel CPU. This is further improved when the differential rotation is minimal, in which case the perturbative expansion may be turned off to save a factor of several in runtime. In stably stratified zones and those with magnetic fields up to 103s may be required to achieve good convergence. In cases where the code has more difficulty, it is quite likely that Mixer becomes the bottleneck in simulations and so, under these circumstances, we recommend tabulating results in advance. This is still considerably more performant than direct numerical simulation, and the results can generally be guaranteed to converge at much higher precision, so that derivatives may be extracted as well. At various points in the software, we must divide by the magnitude of the velocity of an eigenmode. This may approach zero in some cases. To avoid dividing by zero in these cases, we place a lower bound on this magnitude, such that   \begin{eqnarray} |\delta v|^2 \ge \epsilon , \end{eqnarray} (C1)where $$\epsilon = 10^{-20} L_0^2 |N|^2$$ in the calculations presented in this work. This corresponds to setting an upper bound on the length scale d of the displacements $$\delta \boldsymbol {r}$$, namely   \begin{eqnarray} |\delta r|^2 \le L_0^3 |N| \epsilon ^{-1/2}, \end{eqnarray} (C2)which means that d = 1010L0 in this work. To verify that this numerical fix does not impact our results, we have examined the correlation functions in several scenarios as a function of this numerical cut off L. For example, Fig. C1 shows the r − θ and r − ϕ correlations as functions of d for a stably stratified differentially rotating system. The results are constant over many orders of magnitude so long as d > 104L0, which is easily satisfied by our default. Figure C1. View largeDownload slide The absolute values of 〈δvrδvθ〉 (left-hand panel) and 〈δvrδvϕ〉 (right-hand panel) are shown as functions of d, with both axes log-scaled. These results are for a stably stratified region with differential rotation in the radial direction with |R∇ln Ω| = 10−3, Ω = 0.1|N| and no magnetic field. The data is computed for a point on the equator with differential rotation at an angle of π/4. All quantities are given in units of the mixing length and the Brünt–Väisälä frequency. Figure C1. View largeDownload slide The absolute values of 〈δvrδvθ〉 (left-hand panel) and 〈δvrδvϕ〉 (right-hand panel) are shown as functions of d, with both axes log-scaled. These results are for a stably stratified region with differential rotation in the radial direction with |R∇ln Ω| = 10−3, Ω = 0.1|N| and no magnetic field. The data is computed for a point on the equator with differential rotation at an angle of π/4. All quantities are given in units of the mixing length and the Brünt–Väisälä frequency. © 2018 The Author(s) Published by Oxford University Press on behalf of the Royal Astronomical Society

### Journal

Monthly Notices of the Royal Astronomical SocietyOxford University Press

Published: May 1, 2018

## You’re reading a free preview. Subscribe to read the entire article.

### DeepDyve is your personal research library

It’s your single place to instantly
discover and read the research
that matters to you.

over 18 million articles from more than
15,000 peer-reviewed journals.

All for just $49/month ### Explore the DeepDyve Library ### Search Query the DeepDyve database, plus search all of PubMed and Google Scholar seamlessly ### Organize Save any article or search result from DeepDyve, PubMed, and Google Scholar... all in one place. ### Access Get unlimited, online access to over 18 million full-text articles from more than 15,000 scientific journals. ### Your journals are on DeepDyve Read from thousands of the leading scholarly journals from SpringerNature, Elsevier, Wiley-Blackwell, Oxford University Press and more. All the latest content is available, no embargo periods. DeepDyve ### Freelancer DeepDyve ### Pro Price FREE$49/month
\$360/year

Save searches from
PubMed

Create lists to

Export lists, citations