{'fraction_non_alphanumeric': 0.026673339171775303, 'fraction_numerical': 0.010158889027899413, 'mean_word_length': 4.62426118772868, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Astrophysical Fluid Dynamics via Direct Statistical Simulation

In this paper we introduce the concept of Direct Statistical Simulation (DSS) for astrophysical flows. This technique may be appropriate for problems in astrophysical fluids where the instantaneous dynamics of the flows are of secondary importance to their statistical properties. We give examples of such problems including mixing and transport in planets, stars and disks. The method is described for a general set of evolution equations, before we consider the specific case of a spectral method optimised for problems on a spherical surface. The method is illustrated for the simplest non-trivial example of hydrodynamics and MHD on a rotating spherical surface. We then discuss possible extensions of the method both in terms of computational methods and the range of astrophysical problems that are of interest.


Introduction
The modeling of astrophysical phenomena is often limited by the huge range of spatial and temporal scales that need to be resolved in order to describe accurately the dynamics. In many cases the large-scale behavior of cosmic bodies depends on interactions at smaller scales that need to be represented properly for a complete understanding of the astrophysical phenomenon in question. The situation is usually complicated by the requirement of including the back-reaction of the large-scale environment on the smaller scale dynamics in a self-consistent manner. These types of problems are ubiquitous in astrophysics, but we list here some important examples. The transport of angular momentum in accretion disks may be mediated by the (magneto-)rotational turbulence that is present in the disk (see e.g. Balbus & Hawley 1998). This turbulence is itself driven and modulated by the large-scale environment of Keplerian rotation and large-scale magnetic fields (see e.g. Jamroz et al. 2008). The differential rotation pattern in stars (including the sun) arises through an interaction of buoyancy-driven turbulence and rotation, with Reynolds stresses at intermediate scales leading to correlations that drive large-scale flows that themselves act back on the turbulence (Ruediger 1989;Brun & Toomre 2002;Rempel 2005;Miesch 2005). This situation is mirrored in planets, where convective processes may create stresses leading to large-scale flows. Such stresses create turbulence in stably stratified outer weather layers (for example in Jupiter and Saturn) that may drive the formation of jets (see e.g. Scott & Polvani 2008, and the references therein).
There are a number of approaches to modeling the fluid interactions in astrophysical objects. The approach taken often depends on whether it is the dynamics or statistics of the system that is of interest. Sometimes information about the dynamics -that is the precise evolution of a particular realisation of a system -is sometimes required for prediction or to compare with observations. More likely is that the statistics -i.e. the average properties of an ensemble of evolutions -is of interest; this may give more insight into the underlying physics of the system.
Theoretically and computationally, a natural first approach is to perform direct numerical simulations (DNS) of the fluid (or MHD) equations for the system. This approach is the most straight-forward and has led to breakthroughs in many branches of astrophysical fluid dynamics. This approach lends itself naturally to determining the dynamics of a given system. However the extreme nature of the astrophysical turbulent environment ensures that not all spatial scales may be faithfully represented even on the most massive parallel computers available today. For this reason the practitioners of DNS must accept that they are not in the correct parameter regime or may claim that the parameters take into account the effects of scales below the grid cut-off via eddy diffusivities (sometimes termed turbulent transport coefficients). These diffusivities are usually chosen in a plausible but ad-hoc manner. Moreover, DNS may not be an efficient algorithm for determining statistics, since the ensemble over a large number of expensive calculations may be required in order to achieve meaningful statistics.
An alternative approach, which is not useful for determining dynamics but may be useful for statistics, is to derive evolution equations for the large-scale dynamics and to formulate closure models for net effects of the dynamics at moderate and small scales. Such models have a long history in astrophysics and have also achieved some measure of success (see e.g. Kitchatinov & Ruediger (1995); Ogilvie (2003); Rempel (2005)). This approach often utilises (either implicitly or explicitly) moment hierarchies (see e.g. Canuto et al. (1994Canuto et al. ( , 2001Farrell & Ioannou (2008); Garaud et al (2010)). In particular it is customary to relate the average of local interactions of the small scale to the local values of the large-scale fields. The weakness of such models is that they usually rely on some ad-hoc assumption to close the model -parameterising the interactions between large and small scales -and often make the assumption of homogeneity or isotropy. They sometimes find it difficult to include self-consistently the effects of the dynamic large-scale environment. Often it is the case in astrophysics that regions of strong transport lie in close proximity to regions of weak or no transport or mixing -for a nice example see the jets in Jupiter -and so closures that rely on homogeneity may lead to misleading large-scale dynamics. Moreover it is often the case that the inclusion of such closure models introduces new adjustable parameters to the problem that can be tuned to fit observations and that little is known about the sensitivity of the large-scale dynamics to changes in the parameterizations of the scales not captured.
In this paper, we present a new approach to the problem of describing astrophysical flows with a range of spatial scales, which we believe will prove useful for a certain class of problems with large-scale inhomogeneous and anisotropic flows. Specifically we describe the development of efficient numerical algorithms to solve truncated hierarchies of cumulant equations, leading directly to the statistical description of astrophysical flows, and we show that this direct statistical simulation (DSS) is able to reproduce qualitatively statistics obtained by time averaging DNS 1 . That DSS has several advantages over DNS has long been recognized, going back at least as far as a seminal monograph by Lorenz (1967). Low-order statistics are smoother in space and stiffer in time than the underlying detailed flows. Statistically stationary fixed points or slowly varying statistics can therefore be described with fewer degrees of freedom and also can be accessed more rapidly. Convergence with increasing resolution can be demonstrated, obviating the need for separate closure models of the subgrid physics, although these may be included in a natural statistical framework. Finally and most importantly, DSS leads more directly to insights, by integrating out fast modes, leaving only the slow modes that contain the physical information of most interest (Lorenz 1967;Marston 2010).
In this paper we develop the techniques illustrated recently in Marston et al. (2008), where the prototypical problem of barotropic flow relaxing toward a point jet was considered and the statistics obtained by DNS were found to be in good qualitative agreement with those found from a second-order cumulant expansion. We begin by examining the general case of constructing a cumulant hierarchy for the evolution of a number of dynamic variables. We describe the derivation and solution of the cumulant equations for the general case, before focussing the discussion on the case of spherical symmetry, where computational efficiencies are available.
Having described the method in general, we illustrate the advantages of the method for a simple model of the interaction of turbulence and mean flows that may be relevant to the generation of zonal flows in stable layers in planets and stars. The model describes the two-dimensional evolution of flows and magnetic fields on a spherical surface. Such an evolution is non-trivial as it is known that for the hydrodynamic problem the Reynolds stresses act to drive inhomogeneous zonal flows; this type of behavior is difficult to parameterize in sub-grid scale closures. These models and their generalizations have been used to describe the dynamics of the outer layers of giant planets such as Jupiter (see e.g. Scott & Polvani (2008) and the references therein) though competing theories for the generation of zonal flows via deep-seated convection (see e.g. Jones & Kuzanyan (2009) and the references therein) are also available. Furthermore there has been much interest in the MHD version of this problem owing to its importance in the dynamics of the solar tachocline (see e.g. ; D. W. Hughes, R. Rosner, & N. O. Weiss (2007)) and potentially in the outer layers of extra-solar planets (Staehling & Cho (2006)). Both of these environments are believed to be turbulent, stably stratified and magnetized.
The tachocline is believed to play a crucial role in the generation of the eleven year solar cycle (see e.g. Tobias & Weiss (2007)) and angular momentum transport through the tachocline may be responsible for spinning down the solar interior. One crucial issue to be resolved is therefore the role of turbulence in transporting angular momentum in the tachocline, and this has been addressed both in the hydrodynamic and magnetohydrodynamic settings (see e.g. Spiegel & Zahn (1992); Gough & McIntyre (1998);McIntyre (2003); ). What has been shown is that whereas anisotropic hydrodynamic two-dimensional turbulence leads to the efficient formation of zonal flows via Reynolds stresses; the addition of a magnetic field leads to Maxwell stresses that can oppose the formation of jets. The suppression of jets is a function of the strength of the large-scale magnetic field and the local magnetic Reynolds number R m .
The paper is organised in the following manner: in the next section we introduce the general method and the computational savings that can be achieved for the case of spherical symmetry in two dimensions. In section 3 the particular model of MHD turbulence on a rotating spherical surface is introduced and a comparison of the large-scale dynamics of DNS and DSS is made 2 . We conclude by discussing extensions to the method and speculating on the range of problems where such a technique may be of use.

Formulation of the Model
In this section we describe the derivation of a general fully spectral algorithm for the direct statistical simulation of astrophysical flows 3 . We develop the method for the typical case of equations with quadratic nonlinearities, before specialising to systems with spherical symmetry in two dimensions.
Consider a system that is represented by partial differential evolution equations (PDEs) for a number r of scalar fields. Typically such a system may be solved directly by discretising the PDEs using a finite-difference, finite volume or finite element method or by deriving equations for the amplitude of modes in a spectral expansion. Formally, this transforms the PDEs into a finite set of ordinary differential equations (ODEs) that may be integrated forward in time. If the discretisation is performed at s discrete points (or for s spectral modes) then the evolution equations can take the formq where 1 ≤ i, j, k ≤ rs and A i , B ij and C ijk are the coefficients. Here the q i are the discretised values of the dependent variables (or the amplitudes of the relevant spectral modes); typically these represent a vector of the values of the fluid properties. We also note here that there is an implicit sum over repeated indices.
Hereinafter, to fix ideas, we shall think of the q i as representing the amplitudes of the spectral modes of a vector of dependent variables -and shall give a concrete example in the next section. The forcing f i (t) can then be interpreted as the statistical forcing of the relevant spectral mode.

Reynolds decompositions, cumulants and moments
One way to formulate the cumulant expansion is by carrying out a Reynolds decomposition of the dynamical variable q i into the sum of a mean value and a fluctuation (or eddy): where we defer, for now, choosing the type of averaging denoted by the angular brackets . Typical choices are temporal or zonal averages, 4 or averages over an ensemble of initial conditions or an ensemble of realisations.
Once the Reynolds decomposition has been implemented, progress is made by defining the first three equal-time cumulants c i , c ij , and c ijk of the combined scalar fields (q i ) as 5 : where m i , m ij , and m ijk are respectively the traditional definitions of the first, second and third moments. We stress here that the second and higher cumulants contain information about correlations that are non-local in space and therefore include interactions that are not included in the simple local moment hierarchies discussed in the introduction. For this reason this approach is more tailored to inhomogeneous problems.

Derivation of the cumulant hierarchy: the Hopf functional approach
The hierarchy of equations of motions for the evolution of the cumulants can be obtained directly be differentiating Eqs. 3 with respect to time and using Eqs. 1, together with repeated back substitution. A more elegant method is to introduce variables p i that are, in analogy to quantum mechanics, conjugate to the q i in the sense that q i = −i∂/∂p i as in Eq. 7 below (Ma & Marston 2005). Then one may define the Hopf generating functional (Frisch 1995): recalling the summation over repeated indices. The Hopf functional obeys a Schrödinger-like equation: with linear operatorĤ given by: as can be verified by combining Eqs. 4, 5, and 6 to reproduce Eq. 1 in the absence of any stochastic forcing.
As Eq. 5 is linear in Ψ, the average Ψ ≡ Ψ[q(t), p] obeys the same equation; however Ψ encapsulates information about the equal-time moments, as can be seen by repeated differentiation of Eq. 4 with respect to p i , followed by averaging: (For the case of time-averaging, the statistics do not vary in time, ∂ ∂t Ψ = 0, and the statistics are obtained from the solution of the time-independent equationĤΨ = 0.) The Hopf functional Ψ may also expressed as the exponential of a power series in p i , the coefficients being the cumulants: as can be checked by use of Eq. 7 to reproduce the moments in terms of the cumulants, Eqs. 3. Stochastic forcing can now be included with the addition of the Γ ij term: Upon substituting Eq. 8 into Eq. 5 and collecting powers of p i one obtains the equations of motion (EOM) for the cumulants that truncated at third order reaḋ where we defer discussion of the parameter µ until later. Here for compactness we have introduced the short-hand notation {} to denote symmetrization over indices that maintains symmetriesċ ij =ċ ji and similarly for the third cumulant.
Truncated at second-order (CE2) the cumulant expansion is realizable (Salmon 1998) 6 and wellbehaved in the sense that the energy density is positive and the second cumulant obeys positivity constraints. Going to third-order (CE3) and beyond introduces difficulties. A phenomenological eddy-damping parameter (Orszag 1977;Andre 1974) µ that models the neglect of the fourth and higher cumulants from the hierarchy is included in the last of Eq. 10 and is required to prevent blow-up. This ad-hoc procedure is somewhat unsatisfactory and more robust methods may be necessary. Indeed determining reliable methods of truncating the hierarchy is a matter of current research.

Symmetry and the derivation of reduced cumulant equations.
In principle, the general set of cumulant equations in Eq. 10 can be solved with enough computational effort. However, efficient algorithms can be developed if the underlying system exhibits further symmetries. This is typically the case for astrophysical systems, which usually exhibit spherical or cylindrical symmetry or a corresponding translational symmetry in a local Cartesian domain. We discuss in detail here the case of cumulants in a sphere.

Equations for fully spectral DNS on a sphere
For systems with an underlying spherical symmetry, the spectral expansion of the dependent variables discussed in section 2 often takes the form where r is spherical radius, θ is co-latitude and φ is longitude. Here the q m (r) are complex functions and the P m are associated Legendre functions. Furthermore on a spherical surface the r-dependence is absent and a fully spectral representation of the equation of motion (equation 1) can be written asq We note here that, because the scalar fields are real-valued in coordinate space, we may focus on the evolution of modes with m ≥ 0 as modes with m < 0 may be obtained by complex conjugation. Moreover, for simplicity in the above and in subsequent equations the index that encodes which state variable is being solved for has been subsumed into the label.
The quadratic nonlinearities have their origin in the Jacobians of Eqs. 16 with coefficients C (+) representing amplitudes for the scattering of two waves with m ≥ 0; C (−) are for waves with m > 0 and m < 0 to scatter. The amplitudes of these coefficients are constructed from the matrix elements of the Jacobian: where m = m 1 ± m 2 . Integrals I (±) are evaluated in a numerically exact manner by Gaussian quadrature.

Equations for fully spectral DSS on a sphere
Similar considerations lead to an efficient representation of the cumulant hierarchy for a spherical shell. These considerations can then be combined with the knowledge of the underlying symmetries of the statistics themselves to derive reduced hierarchies of cumulant equations. These symmetries are preserved whether zonal, temporal or ensemble averages are used. Statistics on the rotating sphere exhibit azimuthal symmetry. The simplest conceptual choice for the averaging operation is therefore the zonal average and we choose that here, and then follow that with a running time average. On symmetry grounds, the first cumulant must be independent of longitude φ and therefore in the spherical harmonic basis only the m = 0 mode c = q ,m=0 is non-zero. Similar symmetry arguments yield the result that the second cumulant depends on the latitudes of the two field points, but only on the difference between their longitudes. It can therefore be written as c 1 2 m = q 1 m q 2 −m − c 1 c 2 δ m0 . Furthermore, zonal averaging then requires that c 1 2 m=0 = 0. Similarly the third cumulant is a function of only 5, not 6, wavenumbers, i.e. it can be written as c ; 1 m 1 ; 2 m 2 . Moreover, because the scalar fields are real-valued in coordinate space, we have c 1 2 m = c * 2 1 m . For models with an imposed north-south reflection symmetry about the equator 7 , the cumulants respect further constraints: c vanishes for all even and c 1 2 m = 0 if 1 is odd and 2 is even, and vice-versa. All of these symmetries therefore lead to a computational saving.
We consider here the simplest non-trivial case where the hierarchy is truncated at second order (CE2), i.e. all higher cumulants are set to zero. The EOM for the cumulants in the basis of spherical harmonics are thenċ where again the convention of summation over repeated indices has been adopted. (There would also be a contribution to the first cumulant from C (+) 1 ; 0; 0 c c but it vanishes for the problems considered here as the Jacobian of two fields with no longitudinal dependence is zero.) These equations may be compared to a coordinate-independent version given by Equations (21) and (22) in Marston et al. (2008). That only the eddy-mean flow interaction is retained in CE2 may be seen by noting that the coupling of the first cumulant with the second involves no mixing of the azimuthal wavenumber m (only a single m appears in Eqs. 15). Eddy-eddy scattering occurs only at third and higher orders 8 .

Turbulent driven MHD on a Spherical Surface: The Model
We consider a simple two-dimensional model of a stably stratified region of hydrodynamic or MHD turbulence. This is the simplest extension of the local β-plane model considered by . We stress again that, although this system is of interest in its own right and the interaction of Reynolds and Maxwell stresses play an important role in the dynamics of the tachocline and other regimes of stably stratified MHD turbulence, in this paper we are utilising this model as a non-trivial example of the utility of direct statistical simulation. We therefore defer discussion of the interaction of the stresses for a subsequent paper.
The behaviour of such a system in two dimensions can be described by the evolution of two scalar fields, namely the relative vorticity ζ(θ, φ, t) (with θ being co-latitude and φ longitude, as before) and the scalar potential for the magnetic field A(θ, φ, t) (cf . When extended to the sphere rotating at angular rate Ω these may be written: where on the unit sphere the Jacobian is given by Here ζ = ∇ 2 ψ, q = ζ + 2Ω cos θ, Hence q is the absolute vorticity. Here κ is a frictional term, ν 2 is a hyperviscosity whilst f (t) is the stochastic forcing. Magnetic diffusion is explicitly included through a magnetic diffusivity η, since we believe it is important to capture this process correctly. We note here that the parameter B 0 measures the strength of a toroidal imposed magnetic field, which is held fixed in time and note that such a field can not be self-consistently maintained in a strictly two-dimensional calculation. The equations have been scaled so that the magnetic field is measured in units of the Alfvén velocity. For purely hydrodynamic simulations we simply set B 0 = A = 0.
These equations may be then be written in the form of equation (1) (with r = 2) by setting the absolute vorticity and magnetic potential scalar fields into two layers labelled by q α with q 1 = q and q 2 = A and discretising the system either to obtain equations for the spectral amplitudes of the form equation (13) or more conveniently for computation a finite difference representation on a spherical geodesic grid. In a similar manner the spectral representation of the cumulant equations (i.e. equations (15)) can simply be derived once the coefficients in equation (13) have been calculated for this model.

DNS
Direct numerical simulation of the two-dimensional system has been implemented using two different techniques. The EOM given by equation (13) may be integrated forward in time in their pure spectral form using a standard fourth-order accurate Runge-Kutta algorithm with an adaptive time step, though in practice it is much faster work directly in real space on a spherical geodesic grid as we do here. The fully spectral code is therefore only used as a validation of the geodesic code below and a useful comparison with the fully spectral direct statistical simulation.
The most efficient numerical integration of the DNS EOM is carried out in real space on a spherical geodesic grid (Heikes & Randall 1995) of D cells with the use of the second-order accurate leapfrog algorithm and a Robert filter. A multigrid algorithm solves Poisson's equation at each time step.

DSS
We take advantage of the stiff nature of the spectral EOM for the cumulants (equations 15). These are integrated forward in time using a semi-implicit backward Euler Full Orthogonal Method (Saad 2003) that is based upon Krylov subspaces and that permits a much longer time step than is possible for explicit integration methods.
We note here that integration of the EOM for CE2, equations (15), requires of order L 3 M operations at each time step, where 0 ≤ ≤ L and 0 ≤ m ≤ min( , M ) define the spectral cutoffs.
A pseudo-spectral implementation of the EOM would require the same order of operations on the sphere and thus offers no advantages over the pure spectral method used here. We find that all c 1 2 m with m greater than the maximum azimuthal wavevector of the stochastic forcing vanish, hence the spectral expansion can be severely truncated by restricting M L without loss of accuracy. This results in substantial speed-up and a reduction in the required memory. Moreover, only a subset of the possible coefficients of the quadratic nonlinearities, C Finally we note that the code implementing both DNS and DSS (via CE2) is written in the Objective-C++ programming language and runs on Apple computers (OS X 10.6) utilizing Cblocks and grand central dispatch (gcd) for efficient SMP parallelism. We stress that the DSS can run an order of magnitude or more faster than DNS.

Conservation Laws, Model Parameters and Initial Conditions
In the absence of damping and driving forces, the EOM for the cumulants, like the EOM for the vorticity and magnetic potential have a number of conservation laws. For example, in the hydrodynamic case, kinetic energy, enstrophy and angular momentum are conserved, whilst for the MHD case the conserved quantities are angular momentum, total energy, cross-helicity, and the mean squared potential. Moreover, for stochastic forcing restricted to wavevectors | | > 0, the case considered here, the angular momentum in the CE2 remains exactly zero, in contrast to DNS.
Just as for direct numerical simulations utilising spherical harmonics there are convenient expressions of the average values of various quantities in terms of the low-order cumulants. For example the mean cross-helicity is given by: where the two layers are labelled explicitly in the final line. Similar expressions are available for the averages of other quadratic quantities.
The models are formulated on the unit sphere with a timescale such that the sphere complete a full rotation in one day of model time. All model parameters may be defined in terms of these length and time scales; for instance Ω = 2π. Friction removes energy at long length scales and is parameterized by rate κ. The hyperviscosity ν 2 that appears in Eq. (16) is included solely to absorb enstrophy at the smallest resolved scales. Consequently it is rescaled with the grid size or spectral cutoff so that ν 2 → ν 2 * (D/4) −2 geodesic, where the maximum eigenvalue of −∇ 2 on a geodesic grid with D cells is approximately D/2. Thus for ν 2 = 1 (the case we consider here) features on the smallest length scales are dissipated on a time scale of order 1 day.
Stochastic forcing is confined to wavevectors L min ≤ ≤ L max and M min < |m| ≤ L max . Within this range of wavevectors the forcing f m that appears in Eq. (13) is given by where gaussian(t) is a complex number randomly drawn, for each value of and m, from a normal distribution of zero mean and unit variance that smoothly transitions from one random number to the next over a time period of ∆. We set ∆ = 0.1 which is large compared with the time step, but small compared with advective time scales. Consequently in Eq. (15) we have Γ m = 2F f or L min ≤ ≤ L max and M min ≤ |m| ≤ 0 else.
In the following we hold fixed κ = 0.02, ν 2 = 1, F = 0.2, and (for the magnetic cases) η = 10 −4 . We study the evolution of the systems (DNS and DSS) for two different choices for the range of the forcing wavevectors {L min , M min , L max }.
We close our description of the set-up of the models by commenting on the choice of initial condition. The DNS integrations are started from rest with zero perturbation to the imposed field. For the DSS, at the start of the CE2 integration we set the first cumulant c = 0 and the second cumulant c 1 2 m = c 2 δ 1 2 which corresponds to initial short-ranged correlations in the vorticity. At low resolutions, the fixed point sometimes has jets that move in directions opposite to those found in DNS; this fixed point, which is an artifact of the spectral truncation, can be avoided by initializing c with small values.

Results
5.1. Small-scale forcing: L min = 8, M min = 8, and L max = 12 We begin by considering the hydrodynamic and magnetohydrodynamic evolutions for the case where the system is forced solely at small scales in the vorticity equation. The DNS of the hydrodynamic case is performed until a statistically steady state is reached and meaningful statistics can be calculated. For this case, this has occurred by t ∼ 1000; after this time a running time-average is performed for another 1000 days. Here the small-scale driving leads to the formation of flows on a range of scales including large-scale jets as shown in Figure 1 (top panels) which show the instantaneous relative vorticity (left) and relative zonal velocity (right). These clearly show the formation of a prograde (westerly) jet at the equator with two retrograde (easterly) jets at high latitudes, with the total angular momentum of the fluid remaining close to zero. As we shall see, these jets are driven by the flows on smaller scales.
The history and statistics of these hydrodynamic jets for DNS is displayed in the timelines in the upper panels of Figure 2. In the left portion of each panel the relative vorticity and relative zonal velocity (averaged over a period of 10 days) is shown as a function of latitude and time. At 1000 days (half-way through the evolution -signified by a vertical line in the figures) temporal averaging is switched on and a running average from that point is displayed in the figures. This running average eventually settles down to show the mean position and strength of the jets. The direct statistical simulation achieves remarkable agreement with the DNS in both the position and strength of the jets. This is confirmed in Figure 3 which demonstrates that the time-averaged zonal mean zonal velocity as calculated by DSS agrees well with DNS except at high latitudes. Moreover, whilst the DSS respects the north-south symmetry as expected, for the DNS the average position of the prograde jet is slightly off-equator, reflecting the finite length of data over which the averages are calculated. Figure 3 also demonstrates that good convergence with increasing resolution is achieved, both for DNS and for DSS.
That DSS and DNS agree well is reflected in the data recorded in Table 1. There we give some non-dimensional ratios that can only be calculated once the kinetic and magnetic energy are in a statistically steady state. 9 These are defined as where we recall that Ω = 2π, L = 1 and η = 10 −4 . That DSS is able to reproduce the jets using a cumulant hierarchy truncated at second order is an interesting result. It is evidence that the forward enstrophy cascade and anisotropic backward energy cascade (Kraichnan & Montgomery 1980;Salmon 1998), which is frequently invoked to explain their existence, is in fact not necessarythere can be no cascade in the absence of eddy-eddy interactions. We note here that the enstrophy cascade argument has also been questioned in the context of planetary atmospheres (Vallis 1992;Schneider & Walker 2006;O'Gorman & Schneider 2007) where shearing and modification of the thermal structure of the atmosphere by eddy fluxes weaken the eddy-eddy interactions. Here it is therefore Reynolds stresses that are primarily responsible for the build-up of the zonal flows. This result is important for our understanding of the driving of zonal flows in planetary atmospheres.
We now examine the effects of including a toroidal magnetic field and examine the dynamics for two imposed field strengths, B 0 . For B 0 = 0.1 the onset of jet formation is delayed but in both DNS and DSS the system eventually settles into a statistically steady state as shown in Figure 4. For both methods, for this relatively weak imposed mean field, there is eventually little suppression of the jets, with slightly more suppression occurring in the DSS. For this choice of parameters, the magnetic energy is small compared with the kinetic energy of the flow (4% for both DNS and DSS (CE2)), and so it is to be expected that the role of the magnetic field will be secondary. Moreover the magnetic field has been expelled to high latitudes by the strong jets and turbulence at low latitudes (see Figure 7). This flux expulsion (Weiss 1966;Tao et al. 1998) leads to separated regions with different dynamics; at low latitudes, where the field is weak, the hydrodynamic evolution continues unimpeded, whilst at high latitudes the magnetic field leads to some suppression of the jets.
At B 0 = 0.5, however, strong qualitative changes are plainly evident as the jets are destroyed by the fluctuations in the magnetic field, both in DNS and in DSS, in agreement with the findings of . Small remnants of the jets persist in DSS at high latitudes, where the imposed toroidal field is weakest (see Fig. 5). DNS results (top panels) show the incoherent nature of the flows -it is this that leads to the suppression of the jets. For this case, the magnetic energy is in approximate equipartition with the kinetic energy, as shown in Table 1. Once again DNS and DSS show a remarkable agreement for this case; see Fig. 6. A comparison of the mean toroidal component of the magnetic field is shown in Fig. 7. Here the situation is reversed from the previous 9 We note that sufficient averaging must be employed in order to obtain meaningful averages in both DSS and DNS. This is easy to achieve for DSS but is problematic for DNS. weaker field case. The magnetic field here is too strong to be expelled by the eddies and the jet never forms at low latitudes. Therefore the field is confined to low latitudes and the (weaker) jets can only be found at high latitudes where the imposed field is weaker. For this strength of imposed field the kinetic energy is reduced (see the values of Rm in Table 1) and the magnetic energy comes into equipartition with the kinetic energy. Clearly the transport of angular momentum by the Reynolds and Maxwell stresses and of magnetic flux by the turbulent advection has acted in a very different manner here. The discussion of these processes and their description via the second cumulants is postponed to a subsequent paper.
5.2. Case L min = 8, M min = 1, and L max = 10 We now consider the effect of reducing the minimum stochastic forcing wavevector in the azimuthal direction, M min , down to wavenumber 1. This brings stochastic effects to larger scales and so presents a more robust challenge for DSS. A comparison of the zonal mean relative vorticity and zonal velocity as calculated by DNS and DSS (CE2) is shown in Fig. 8 for the hydrodynamic problem. In contrast to the previous case, there are two prograde jets at high latitudes, and one equatorial retrograde jet; again the total angular momentum is close to zero. Now, however, the jets are seen to wander significantly in latitude in DNS owing to the continual random forcing at large zonal scales. Once established in DSS, however, they remain fixed in place. As a consequence, the time-averaged zonal means are reduced in magnitude in DNS compared with DSS, as made apparent in the quantitative plot of Fig. 9.
Imposing a relatively weak toroidal field by setting B 0 = 0.1 again has little effect on the zonal mean velocity as shown in Figure 9. However at B 0 = 1 the jets are largely eliminated by the fluctuations in the magnetic field, both in DNS and in DSS. Somewhat larger jet remnants remain at high latitudes, where the imposed toroidal field is weakest, and again stronger jets are found in DSS than in DNS (see Fig. 10).
The toroidal field is tightly confined to latitudes less than roughly 60 • , as Fig. 11 depicts, likely owing to flux expulsion of the field. Again DSS does a reasonably good job of reproducing the mean field.

Discussion
In this paper we have introduced the concept of Direct Statistical Simulation (DSS) for astrophysical fluid dynamics. We have compared the results of Direct Numerical Simulation (DNS) and DSS for the problem of two-dimensional hydrodynamics and magnetohydrodynamics on a spherical surface. Although the set-up of the model is relatively simple, the ensuing dynamics is not. In the hydrodynamic case, non-trivial interactions at moderate scales drive inhomogeneous large-scale zonal flows (jets/winds). With a weak imposed field the jets remain largely unaffected and the magnetic fields are expelled to higher latitudes. With a stronger imposed toroidal magnetic field these winds are suppressed except at high latitudes where the imposed field is weak.
We find that even the simplest formalism of DSS, based upon the truncation of the cumulant hierarchy at second order, is capable of reproducing the driving and suppression of the zonal flows and the flux expulsion of the magnetic fields by the inhomogeneous jets. Because the method includes interactions that are non-local in space it is very well suited to such inhomogeneous problems thatypically arise in astrophysics. Such a truncation is equivalent to keeping the mean/eddy interactions in the eddy equations and the eddy/eddy interactions in the mean equations, whilst suppressing the eddy/eddy interactions in the eddy equations. Thus it is the Reynolds and Maxwell stresses that respectively drive and suppress the jet, not an inverse cascade as is frequently assumed.
The DSS scheme is more numerically efficient than the corresponding DNS. We believe that the results presented here are an encouraging beginning for the concept of DSS in astrophysical fluid dynamics. It is important though to determine the range of validity of such a procedure. Clearly the method is designed to work best when the dynamics leads to the generation of substantial statistical means (e.g. mean flows or magnetic fields) or involves the interactions of prescribed (usually inhomogeneous) mean quantities with smaller scale turbulence. The method is inefficient when the turbulence is dominated completely by small-scales and is largely homogeneous -for example homogeneous, isotropic hydrodynamic turbulence, or the small-scale dynamo problems (see e.g. Tobias et al. 2011). We do believe, however, that many cases of astrophysical interest do fall into the category where DSS techniques may prove useful. Examples currently under consideration include the interaction of mean magnetic fields and shear flows either on a spherical surface (leading to joint instability (see e.g. Cally et al. 2003)) or in a cylindrical domain (leading to magnetorotational instability), the instability and mixing of large-scale shear flows in the presence of a magnetic field, and the driving of zonal flows via convection in a tilted cylindrical annulus (see e.g. Brummell & Hart 1993). It will be interesting to determine how well the techniques described in this paper fare for these problems, and we predict varying degrees of success. We also stress that even utilising DSS will not allow the calculation of statistics at astrophysically realistic values. However it is to be hoped that, whereas the dynamics may be extremely sensitive to the parameters, for a range of problems the statistics may prove less so. We believe that some of these problems may require the inclusion of higher order cumulants in the scheme and are currently engaged in determining efficient numerical procedures for their integration.
Clearly in the longer term, if these techniques prove useful for the simpler problems described above, it will be of interest to apply them to more computationally intensive problems in astrophysical fluids. These include the driving of zonal flows in planets, the mixing of angular momentum and abundances in stellar interiors and the transport by turbulence in accretion disks. JBM thanks P. Kushner and T. Schneider for helpful discussions. The authors are grateful to F. Sabou and W. Strecker-Kellogg for assistance with the semi-implicit Krylov method. The work is supported in part by NSF grant DMR-0605619 (JBM and KD).