# 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 nonlinear 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 (MRI). It also exhibits non-linear saturation behaviour, and we use this to extract the asymptotic scaling of various transport coefficients in physically interesting limits.

###### keywords:

convection - Stars; rotation - Stars; interiors - Stars: evolution^{1}

^{2}

## 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 which 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 timescales 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 timescales of minutes and nuclear timescales 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 et al., 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 et al. (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 nonetheless distinct from attempts to simulate turbulence in all its detail. Some may introduce new dynamical fields, as in the 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 linearised 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 et al., 1997; Lohse & Xia, 2010). With some notable exceptions (Sukoriansky et al., 2007) this cascade begins at a large scale set by the overall structure of the fluid flow and ends at an extremely small scale 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 for a given turbulent flow which fully characterises the turbulence as seen by measurements performed over length scales . This is the modern interpretation and justification of the original mixing length hypothesis, which asserts that turbulent fluctuations on scales are not dynamically coupled to the large-scale () 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

(1) |

where is the velocity, denotes the outer product, denotes the time-averaged expectation, is the amplitude of the Fourier mode with wavevector and is the tensor specifying how different components of the same mode are correlated with one another. It is crucial to notice that the quantity is also the Reynolds stress of mode . This, and several other closely related quantities, are ultimately what we seek. These two-point correlation functions suffice to characterise not only the stresses but also all higher-order correlations through Wick’s theorem and perturbation theory (Wick, 1950; Isserlis, 1918).

To determine we begin by writing the linearised equations of motion as

(2) |

where is a linear operator of its first argument and 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 may be treated as translation-invariant and so we may compute a Fourier transform in without coupling different modes. This gives

(3) |

The modes are decoupled in this regime so can be represented by a matrix , and we write

(4) |

When is independent of equation (4) is straightforward to solve and gives us

(5) |

where are the initial mode amplitudes and and are respectively the normalised right eigenvectors and eigenvalues of . The vectors 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

(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 and integrating over gives us

(7) |

Generally some 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 (Pessah & Goodman, 2009; Lesaffre et al., 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 . That is,

(8) |

where depends on the large scale properties of the flow but is independent of , is the number of modes per wavevector and is the index of the turbulence. Following Kolmogorov (1941a) we choose in our model. Appendix A contains a detailed discussion of this choice.

The wavenumber is just that of the characteristic scale, and is given by

(9) |

Replacing the divergent expression in equation (7) with this amplitude we find

(10) |

It only remains to determine . To do this we note that there is one characteristic length scale and one characteristic timescale, the growth rate of the mode. Because has dimensions of velocity squared we find

(11) |

where 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 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 only acts 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 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 .

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

(12) |

This series can be truncated to produce an approximation of which is accurate in a certain window around .

We may likewise write the velocity at a given wavevector as

(13) |

This suggests defining a new vector

(14) |

which, in principle, encodes the full time evolution of the velocity field. This vector evolves according to

(15) |

where is formed of blocks given by

(16) |

By definition though we also have

(17) |

where , 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 . Doing so makes the assumption that the behaviour of the system at all greater is known. Inspired by the solution for time-independent , we try an exponential behaviour. This truncates equation (17) such that it applies only to and means that we are searching for vectors with

(18) |

and

(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 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

(20) |

where and are matrices acting only on the -block. For example, in the case where , our equations are

(21) | ||||

and | ||||

(22) |

which may be put in the form of equation (20) with

(23) | ||||

and | ||||

(24) |

The eigenvectors of this system are solutions of the original equation (4) because if is such an eigenvector then

(25) |

solves

(26) |

over the time window for which is well-approximated at -th order. As a result we say that are the instantaneous modes of the system at -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

(27) |

This approximation is controlled in the sense that so long as converges as grows, so does the inferred velocity history. In this work we present results with so that involves both and . We leave the exploration of larger 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^{3}

(28) |

In a fixed coordinate system differential rotation is difficult to analyze so we make two reference frame changes. First we switch from an inertial frame to one rotating at

(29) |

Secondly we make a formal change of coordinates

(30) |

without altering the corresponding unit vectors. Under this last change the gradient transforms as

(31) |

Because the operator is most easily expressed in Fourier space we define the transformed wavevector as

(32) |

With this the transformed MHD and Navier-Stokes equations may be written as

(33) |

and

(34) |

where is the specific entropy and

(35) |

is the pressure tensor with 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 . 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

(36) |

so in Fourier space

(37) |

Combining this with equation (33) and the Boussinesq approximation (see Appendix B) we find

(38) |

Note that as did Balbus & Schaan (2012) we take 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 to eliminate and find

(39) |

where the notation denotes the component perpendicular to .

To construct the matrix version 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 we choose the unit vectors

(40) | ||||

and | ||||

(41) |

where is any unit vector with . This choice of basis ensures that our vectors are perpendicular to the wavevector.

A choice of particular convenience for is

(42) |

With this choice is time-independent, because the component of perpendicular to is time-independent, and so we may write

(43) |

and

(44) |

Note that there is a removeable singularity when . The matrix is then given by computing the relation between and . 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

(45) |

Combining the linearised equations of motion with our closure scheme we can compute the correlation function

(46) |

where the index ranges over eigenvectors. This function contains all of the usual stresses and transport functions. For instance, the Reynolds stress is

(47) |

Likewise up to a dimensionless constant of order unity the turbulent diffusivity is

(48) |

and the turbulent viscosity is

(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 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 et al., 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 , and 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ä et al. (2004). In order to match the units for the rotation rates we put everything in terms of the coriolis number

(50) |

where, following the convention of Käpylä et al. (2004), was computed for a non-rotating system.

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 , which is seen in these and other simulations (Rüdiger et al., 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 , 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.

Despite these difficulties some trends are clear and sustained between both sets of data. For instance in the northern hemisphere (), , while in both hemispheres , in keeping with predictions and simulations by Rüdiger et al. (2005b). Likewise we find that 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.

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 and , which are the correlation functions controlling this transport.

Both correlators vary at second order in in the slow rotation limit as expected (Lesaffre et al., 2013; Kitchatinov, 2013). In the rapid rotation limit on the other hand they exhibit clear 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, stabilising modes.

The peak of each correlator is of order unity and occurs when . In fact for the stress the maximum is while for the diffusivity it is , 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

(51) |

Next we consider the effect of rotation on the 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 and as a function of the rotation rate.

In the slow-rotation regime both quantities scale as , while in the rapid rotation limit they scale as . The peak is of order unity and occurs near . This gives rise to the approximation

(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 and stress and diffusivity correlations as a function of latitude. The 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 stabilise modes.

By contrast the 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 unit vector is ill-defined. The sign change between the northern and southern hemispheres occurs because has the same sign everywhere while 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 and . 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 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 which couple strongly to the Coriolis effect which are stabilised 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).

### 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 is at most of order unity.

Fig. 8 shows the and velocity and diffusivity correlation functions as a function of differential rotation for a situation where is at an angle of 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) are helpful because their model was implemented in a two-dimensional solar model which compared well with helioseismic observations.

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 greatly reduces near , 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 nonlinear behaviour despite being carried out to linear order in . 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 . Fig. 9 shows these derivatives as a function of . The stress derivative is constant in . This means that the stress scales as . 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

(53) |

By contrast the derivatives of the correlations as well as the diffusivity all diverge in the limit as . In particular, the velocity correlation diverges as , the diffusivity correlation diverges as and the diffusivity diverges as . These divergences are signatures of symmetry breaking. They indicate that the direction in which the limit is approached matters. That is, this limit can be approached by first letting and then differentiating or by differentiating and then taking and the divergence we find in the latter approach indicates that the order matters.

When and there is a symmetry between and between . As a result both the and terms vanish in this limit. When these symmetries are broken by the rotation and we know from Figs. 5 and 7 that this occurs at first order for and second order for . In the opposing limit the situation is different because in the time evolution described by equation (39) is independent of when . There is, however, a dependence on through the time-dependence of . This breaks the symmetry because is proportional to and hence is sensitive to . It does not, however, break the symmetry, because is symmetric with respect to changing the signs of both and . It follows then that we should find divergences in the correlation derivatives owing to the path-dependence of the zero-rotation limit and that we should find the derivatives to be generally well-behaved.

The curious divergence is then that in the diffusivity, because this correlation function does not suffer from a symmetry-derived path-dependence. This arises because the differential rotation means that 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 and their growth may proceed in the azimuthal direction until bounded by the Coriolis effect at a time . As a result these modes contribute to the diffusivity as and hence lead to a diverging derivative in as .

### 6.3 Differential Rotation and Stable Stratification

Stably stratified regions are those with

(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

(55) |

is satisfied (Zahn, 1993). Here is the velocity and 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 non-vanishing stress components in a rotating stably stratified zone with latitudinal rotational shear. All 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,

(56) |

where is some function of the angular velocity. Fig. 11 shows the dependence of the stress components on for fixed . The components all scale as in both regimes. Thus, for instance, because the viscosity is the derivative of the stress with respect to the shear, and hence

(57) |

The scaling in equation (57) arises owing to the centrifugal term, which has a destabilising effect when increases with . When this effect is not present so the system is stable but introducing a small differential rotation produces an acceleration proportional to and hence

(58) |

which means that the stress scales as and thence the viscosity scales as .

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 , 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 non-vanishing stresses are quadratic in the shear, 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 , because the mode would need to be an eigenvector of . 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 stabilising 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 relative to the vertical. These angles represent a compromise between maximising the magnitude of the centrifugal acceleration and maximising its projection on to the velocity, both subject to the Boussinesq condition that motion be in the plane perpendicular to .

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 maximise the centrifugal term while the velocity can be chosen to maximise the projection of the acceleration on to the velocity.

### 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 and 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 .

By contrast the correlations behave approximately as , and vanishes when . 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 . 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 in a mildly rotating () convection zone as a function of for three polarisations; radial (), latitudinal (), and longitudinal (). As the field increases the stress falls off. This is because the field quenches the turbulence by providing a stabilising 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.

By contrast consider , 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 polarisation 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 polarisation relative to the others reflects the fact that motion is damped perpendicular to the magnetic field so the polarisation damps motion in both directions involved in this component of the stress whereas the and polarisations only damp motion in one of the two directions.

### 6.6 Magnetorotational Instability

As a final example we consider the magnetorotational instability (Chandrasekhar, 1960). This instability arises in magnetised fluids undergoing Keplerian orbital motion.

Fig. 19 shows the and Reynolds and Maxwell stresses for an accretion disc with a vertical magnetic field. Contrary to predictions (Chandrasekhar, 1960) not all 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 . 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 . This is to be expected because they are proportional to . Additionally, the - stresses vanish for all because the system is symmetric under reversing both and , and similarly the - stresses vanish because the system is symmetric under the simultaneous reversal of and .

As the magnetic field increases the Reynolds stress changes sign. This indicates the onset of MRI modes, which have the opposite sign to the zero-field correlations. This effect saturates when , where is the scale height of the disc. The total stress saturates at roughly , 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 et al., 2006). Furthermore there are challenges with the -disk prescription which make the specific stress components more difficult to interpret (Pessah et al., 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 timescales.

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 IOA, ENS and 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 STFC Grant ST/L003910/1 for CPU cycles which aided in this work.

## 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 (Kolmogorov, 1941a). There is more debate over the index to use for convection, with answers ranging from (Benzi et al., 1994) to (Procaccia & Zeitak, 1989) and (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 et al., 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 that it is divergent for and not for . This favours the scenario of Goldreich & Sridhar (1995), who argue that in the strongly-magnetised limit the index ought to be .

In order to consistently treat both the non-magnetic and the strongly-magnetised limits, we choose a simple prescription in which when one of , or exceeds and use otherwise. This means that there is a critical wavenumber

(59) |

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

(60) |

Taking the time derivative of both sides we see that

(61) |

As a result

(62) |

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 visualisation of the data was with the Python modules Numpy (van der Walt et al., 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 et al. (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 only supports single-threaded operation, though it may be used inside parallelised 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 optimised for convecting systems for which achieving accuracy better than relative and absolute typically requires between and 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 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

(63) |

where in the calculations presented in this work. This corresponds to setting an upper bound on the length scale of the displacements , namely

(64) |

which means that 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 cutoff . For example, figure 20 shows the and correlations as functions of for a stably stratified differentially rotating system. The results are constant over many orders of magnitude so long as , which is easily satisfied by our default.

### Footnotes

- pubyear: 2017
- pagerange: Turbulence Closure for Mixing Length Theories–C
- It would not be difficult, however, to incorporate them into this framework at a later date.

### References

- Ashkenazi S., Steinberg V., 1999, Physical Review Letters, 83, 4760
- Balbus S. A., Hawley J. F., 1991, ApJ, 376, 214
- Balbus S. A., Schaan E., 2012, MNRAS, 426, 1546
- Benzi R., Tripiccione R., Massaioli F., Succi S., Ciliberto S., 1994, EPL (Europhysics Letters), 25, 341
- Berntsen J., Espelid T. O., Genz A., 1991, ACM Trans. Math. Softw., 17, 437
- Böhm-Vitense E., 1958, Z. Astrophys., 46, 108
- Canuto V. M., 1994, ApJ, 428, 729
- Canuto V. M., 1997, ApJ, 482, 827
- Canuto V. M., Hartke G. J., 1986, A&A, 168, 89
- Canuto V. M., Cheng Y., Howard A. M., Esau I. N., 2008, Journal of Atmospheric Sciences, 65, 2437
- Carati D., 1990, Phys. Rev. A, 41, 3129
- Chan K. L., 2001, ApJ, 548, 1102
- Chandrasekhar S., 1960, Proceedings of the National Academy of Science, 46, 253
- Dobrowolny M., Mangeney A., Veltri P., 1980, Phys. Rev. Lett., 45, 144
- Floquet G., 1883, Annales scientifiques de l\textquotesingleÉcole normale supérieure, 12, 47
- Garaud P., Ogilvie G. I., Miller N., Stellmach S., 2010, MNRAS, 407, 2451
- Garaud P., Gagnier D., Verhoeven J., 2017, The Astrophysical Journal, 837, 133
- Genz A., Malik A., 1980, Journal of Computational and Applied Mathematics, 6, 295
- Goldreich P., Sridhar S., 1995, ApJ, 438, 763
- Gough D., 1977, in Spiegel E. A., Zahn J.-P., eds, Lecture Notes in Physics, Berlin Springer Verlag Vol. 71, Problems of Stellar Convection. pp 15–56, doi:10.1007/3-540-08532-7_31
- Gough D. O., 2012, ISRN Astronomy and Astrophysics, 2012, 987275
- Guennebaud G., Jacob B., et al., 2010, Eigen v3, http://eigen.tuxfamily.org
- Hunter J. D., 2007, Computing in Science Engineering, 9, 90
- Isserlis L., 1918, Biometrika, 12, 134
- Käpylä P. J., Korpi M. J., Tuominen I., 2004, A&A, 422, 793
- Kichatinov L. L., 1986, Geophysical and Astrophysical Fluid Dynamics, 35, 93
- Kichatinov L. L., 1987, Geophysical and Astrophysical Fluid Dynamics, 38, 273
- Kichatinov L. L., Rudiger G., 1993, A&A, 276, 96
- Killworth P. D., 1980, Dynamics of Atmospheres and Oceans, 4, 143
- King A. R., Pringle J. E., Livio M., 2007, MNRAS, 376, 1740
- Kitchatinov L. L., 2013, in Kosovichev A. G., de Gouveia Dal Pino E., Yan Y., eds, IAU Symposium Vol. 294, Solar and Astrophysical Dynamos and Magnetic Activity. pp 399–410 (arXiv:1210.7041), doi:10.1017/S1743921313002834
- Kolmogorov A., 1941a, Akademiia Nauk SSSR Doklady, 30, 301
- Kolmogorov A. N., 1941b, Akademiia Nauk SSSR Doklady, 32, 16
- Launder B., Spalding D., 1974, Computer Methods in Applied Mechanics and Engineering, 3, 269
- Lebovitz N. R., 1965, ApJ, 142, 1257
- Lee D., 2013, Journal of Computational Physics, 243, 269
- Lesaffre P., Balbus S. A., Latter H., 2009, MNRAS, 396, 779
- Lesaffre P., Chitre S. M., Potter A. T., Tout C. A., 2013, MNRAS, 431, 2200
- Lohse D., Xia K.-Q., 2010, Annual Review Of Fluid Mechanics, 42, 335
- Luschgy H., Pagès G., 2006, Moment estimates for LÃ©vy Processes (arXiv:math/0607282)
- Maeder A., 1995, A&A, 299, 84
- McKinney J. C., Tchekhovskoy A., Sadowski A., Narayan R., 2014, MNRAS, 441, 3177
- Murphy G. C., Pessah M. E., 2015, ApJ, 802, 139
- Pessah M. E., Goodman J., 2009, ApJ, 698, L72
- Pessah M. E., Chan C.-K., Psaltis D., 2006, MNRAS, 372, 183
- Pessah M. E., Chan C.-K., Psaltis D., 2008, MNRAS, 383, 683
- Procaccia I., Zeitak R., 1989, Physical Review Letters, 62, 2128
- Rüdiger G., Egorov P., Ziegler U., 2005a, Astronomische Nachrichten, 326, 315
- Rüdiger G., Egorov P., Kitchatinov L. L., Küker M., 2005b, A&A, 431, 345
- Salvesen G., Simon J. B., Armitage P. J., Begelman M. C., 2016, MNRAS, 457, 857
- Schou J., et al., 1998, ApJ, 505, 390
- Smolec R., Houdek G., Gough D., 2011, in Brummell N. H., Brun A. S., Miesch M. S., Ponty Y., eds, IAU Symposium Vol. 271, Astrophysical Dynamics: From Stars to Galaxies. pp 397–398 (arXiv:1011.3813), doi:10.1017/S1743921311017972
- Sridhar S., Goldreich P., 1994, ApJ, 432, 612
- Starling R. L. C., Siemiginowska A., Uttley P., Soria R., 2004, Monthly Notices of the Royal Astronomical Society, 347, 67
- Sukoriansky S., Dikovskaya N., Galperin B., 2007, Journal of Atmospheric Sciences, 64, 3312
- Wick G. C., 1950, Phys. Rev., 80, 268
- Yakhot V., Orszag S. A., 1986, Journal of Scientific Computing, 1, 3
- Zahn J.-P., 1993, Space Sci. Rev., 66, 285
- Zhou Y., McComb D. W., Vahala G., 1997
- van der Walt S., Colbert S. C., Varoquaux G., 2011, Computing in Science Engineering, 13, 22