{'fraction_non_alphanumeric': 0.031607731240870955, 'fraction_numerical': 0.011838037976833452, 'mean_word_length': 4.636511583381199, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Numerical Simulations of Magnetic Fields in Astrophysical Turbulence

The generation and evolution of astrophysical magnetic fields occurs largely through the action of turbulence. In many situations, the magnetic field is strong enough to influence many important properties of turbulence itself. Numerical simulation of magnetized turbulence is especially challenging in the astrophysical regime because of the high magnetic Reynolds numbers involved, but some aspects of this difficulty can be avoided in weakly ionized systems.


Introduction
The interaction of magnetic fields with turbulence is a basic feature of many astrophysical systems. Important, basic MHD turbulence problems common to many fields of astrophysics include the nature of MHD turbulence itself, the dynamo problem, the effects of magnetic fields on turbulent transport and turbulent mixing, the effects of small scale turbulence on large scale dynamics, and the formation and evolution of current sheets.
Although progress on all of these problems has been made analytically, numerical simulations are an increasingly powerful means of approaching them. Accurate simulation of astrophysical magnetic fields under turbulent conditions presents extreme challenges of its own. The main reason is the smallness of the magnetic diffusivity, which leads naturally to the formation of thin current layers. Most of the Ohmic dissipation, and most of the change in magnetic topology, takes place in these layers. It is important to understand these diffusive effects in order to answer questions such as: How do dynamos amplify magnetic fields on large scales but not small scales? What determines the ratio of magnetic flux to mass in self gravitating regions? How is magnetic energy dissipated in a turbulent plasma?
The purpose of this review is to isolate some important problems related to magnetized turbulence in astrophysics, with emphasis on their computational aspects. The subject is vast, and we make no claim to be comprehensive. We cite literature up to early 2002.
In §2, we write down the magnetic induction equation, derive some of its basic properties, and discuss the use of conservation laws in testing MHD codes.
As an illustration, we carry out such tests on the ZEUS code. In §3, we discuss the parameter space for astrophysical MHD. In §4, we discuss some results on the dynamo problem, the turbulence problem, the formation of singularities, and the role of turbulence in large scale dynamics. In §5, we discuss lightly ionized media, for which ambipolar drift yields an effective diffusivity which is much higher than the Ohmic value. Although the two forms of diffusion are not equivalent, ambipolar drift is a partial solution to the problem of small diffusivity. Section 6 is a summary and discussion of future prospects.

Induction Equation and Consequences
According to Faraday's law, also known as the magnetic induction equation Throughout this paper, we will assume E is given by where σ is the electrical conductivity, λ Ω ≡ c 2 /4πσ is the magnetic diffusivity, and in the last step we have used Ampere's law, neglecting the displacement current.
The first and second terms on the RHS of eqn. (2) represent inductive and resistive effects, respectively. In order to compare these terms, we introduce a characteristic speed V 0 and a characteristic lengthscale L 0 , and write v and r in terms of dimensionless velocities and coordinates; v ≡ V 0 u; r ≡ L 0 s. Equation (2) can then be written as where the dimensionless parameter R m is defined by In some problems it is convenient to set V 0 to a typical Alfvén speed v A ≡ B/(4πρ) 1/2 , in which case R m is known as the Lundquist number and denoted by S. The case λ Ω ∝ S −1 ≡ 0 is called ideal MHD. It is clear from eqns. (1) and (2) that the limit S → ∞ is a singular limit, in the sense that at this limit the order of the magnetic induction equation drops from second to first. We therefore expect that as S → ∞, thin boundary layers, or current sheets, will form. This should be anticipated when choosing numerical schemes ( §3.3).
The high S limit describes most astrophysical problems. The ideal form of the magnetic induction equation can be solved exactly in terms of fluid trajectories.
Suppose that the fluid at position r at time t was at position r 0 at time 0. Define the deformation matrix D by It can then be shown that the magnetic field B(r, t) is related to the initial field B 0 (r 0 (r, t), 0) by where |D| is the determinant of D. Equation (6) is known as the Cauchy solution of the magnetic induction equation. It is the basis for the numerical technique used by Kinney et al. [53] to simulate 2D MHD turbulence. Conservation laws can be used to test MHD codes. The Cauchy solution embodies magnetic flux conservation. This basic property is in practice difficult to test, because in turbulent flow, fluid elements follow complex paths. The Cauchy solution only makes sense as long as fluid trajectories do not intersect, but this cannot be guaranteed in a finite difference code, which is always somewhat diffusive. Therefore, we turn to a globally conserved quantity: magnetic helicity.
The helicity H of a fixed volume V of fluid with magnetic vector potential A and magnetic field B = ∇ × A is Uncurling eqn. (1) leads to an evolution equation for A where φ is a free gauge function. According to eqns. (1) and (8), the rate of change of H is where we have integrated once by parts and used Gauss's theorem. We now take periodic boundary conditions on V , so that the surface integral vanishes, and use eqn. (2). Equation (9) then becomes Equation (10) shows that helicity is conserved in an ideal medium. The rate at which helicity varies is a global measure of the magnetic diffusivity λ Ω . Diffusion can cause helicity growth as well as decay. We obtain a slightly different conservation law if we assume that V is comoving instead of fixed. In this case, we find that in an ideal medium, H is conserved within V as long as B · dS ≡ 0. However, in view of the difficulty of following comoving volumes in a turbulent fluid, it is more useful to treat V as fixed, and to take it as the computational domain.
Finally, we derive an equation for magnetic energy W Taking the scalar product of eqn. (1) with B, integrating over space, and assuming periodic boundary conditions yields The first term on the RHS of eqn. (12) represents the work done by the flow on the field, and appears with opposite sign in the evolution equation for kinetic energy. The second term represents energy loss by Ohmic decay.

Helicity Conservation in the ZEUS Code
The ZEUS-3D code [99,100] solves the equations of ideal, compressible MHD using a finite difference scheme and a von Neumann artificial viscosity to capture shocks. The MHD induction equation is followed using the method of consistent transport along characteristics [46]. ZEUS-3D is publically available, and has been of great service in the astrophysical community. If there were no numerical dissipation in the ZEUS code, magnetic helicity would be strictly conserved (see eqn. (10)). Here, we investigate helicity conservation in the ZEUS code in two applications.
The first problem is the evolution of a twisted magnetic flux tube which is unstable to the kink mode. Initially, the helicity of the tube is entirely in the twist of the field about the axis. Theory predicts that the instability drives the system to a new equilibrium state in which some of the helicity is carried by a writhing deformation of the tube axis, and the total magnetic energy is reduced while the total helicity is fixed.
The simulations confirm the broad outline of this picture. The growth rate of the kink during the linear phase agrees with an analytical calculation, and significant motion occurs only during the kinking phase, while magnetic energy is being released. However, both magnetic energy and magnetic helicity decline steadily once the system reaches equilibrium, as shown in Figure (1) for computations with 80 3 and 160 3 gridpoints. During the dynamical kink phase, the energy declines much faster than the helicity, confirming the importance of both dynamical and resistive processes in the evolution of magnetic energy.
We have used eqn. (10) to estimate the mean pseudo-magnetic diffusivity λ Ω , i.e. the mean numerical magnetic diffusivity, by writing the integral on the right hand side as which can be regarded as the definition of λ Ω . The result is shown in Figure  (2), where λ Ω is plotted for the two simulations in units of ∆xv A0 , where ∆x is the grid scale. The near coincidence of the two curves shows that the numerical diffusion is linear in ∆x. During the dynamical phase, λ Ω is enhanced over its "quiet" value ǫ 0 ∆xv A0 by about a factor of 3, perhaps suggesting that numerical resistivity, like numerical viscosity, is proportional to velocity. The results shown in Figure (2) allow us to estimate the Lundquist number S. According to eqn. (4), S can be written in terms of ǫ 0 , the number of gridpoints N , and the magnetic lengthscale L B in units of the box size L as In these computations, the tube has an initially Gaussian magnetic profile with L B /L = 0.1. Taking ǫ 0 ∼ 0.01, which is representative of the "quiet" value shown in Figure (2) and N 1/3 ∼ 100, we see that S 10 3 . We have also computed the evolution of H in models of molecular clouds [47,48], with initially uniform magnetic fields, stirred by supersonic turbulence. In these models H is initially zero, and should remain so, although the helicity density A · B can vary arbitrarily between positive and negative values from point to point. Figure (3) shows H relative to E m L for two different runs at various times, measured in units of the sound crossing time. In these units, the errors in H are at the 1% level, and, reassuringly, there is no evidence that H drifts steadily away from zero. The variance of H around zero is smaller for the run at Alfven Mach number M A = 0.7 than for the run at M A = 6.7, presumably because the field becomes less tangled if it is relatively strong. Figure (4) shows mass density and helicity density, for the weak field model, in a plane perpendicular to the mean magnetic field. The figure shows that helicity density varies strongly with position, and that it is well correlated with mass density (the correlation is not as good in the strong field case).
The density compressions in these models are associated with shock fronts. Helicity density can also increase across a shock. For example, consider a locally sheared, force free magnetic field of the form The vector potential is so the helicity density is −k −1 0 B 2 0 . Suppose the fluid is compressed in theẑ direction, so that the initial and final coordinates z 0 and z are related by z = αz 0 , with α < 1. It can be shown from eqn. (6), or from intuitive arguments, that B still has the form (15), but with B 0 → α −1 B 0 and k 0 → α −1 k 0 . The helicity density then increases by a factor α −1 , which is the same factor by which the mass density increases.

Parameter Space
The dimensionless parameters of astrophysical turbulence define the scope of the problem and should influence the choice of numerical technique. Here we define and give quantitative expressions for some important parameters.

Macroscopic Parameters
The ratio of gas pressure to magnetic pressure is usually denoted by β ≡ 8πP/B 2 , which is closely related to the ratio of the sound speed v S to the Alfvén speed v A ; v 2 S /v 2 A = γβ/2, where γ is the ratio of specific heats. If β is either very small or very large, the acoustic and Alfvén frequencies are very different, with consequences for choice of simulation technique.
The ratios of the mean flow speed v to v S and v A , respectively, are denoted by the sonic and Alfvén Mach numbers M and M A . When M ≫ 1, compressibility effects are important and the flow is pervaded by shocks, requiring an algorithm which treats them accurately.
The importance of self gravity is measured by the Jeans length λ J . Turbulent pressure effectively increases λ J . This is captured by a heuristic expression for λ J in a turbulent medium where M 2 (λ J ) is the Mach number of the turbulence at wavelengths less that λ J ; turbulence at longer wavelengths does not contribute to pressure support [22,12,90]. Equation (17) reverts to the usual definition of the Jeans length as M → 0. Truelove et al. [106] have shown that unless the size of the grid, ∆x, remains less than about 0.25 of the local λ J , the system is subject to spurious local fragmentation.

Microscopic Parameters
Complete expressions for plasma transport coefficients are given by [13]. We are primarily concerned here with the magnetic and viscous diffusivities, which are conveniently given in terms of the electron and ion collision times τ e and τ i τ e = 2.9 × 10 −2 (Λ/10) T 3/2 n e Z s; τ i = 1.7 (Λ/10) In eqn. (18), T is given in degrees K, Λ is the Coulomb logarithm (Λ = 9.4 − 1.15 log n e + 3.45 log T for T < 5.8 × 10 5 K and Λ = 15.9 − 1.15 log n e + 2.3 log T for T > 5.8 × 10 5 K), and A and Z are the ionic atomic number and charge, respectively. The electron density n e is expressed in cgs units; cm −3 . The magnetic diffusivity λ Ω introduced in eqn. (2) can also be expressed in terms of the electron skin depth δ e ≡ c(m e /4πn e e 2 ) 1/2 = 5.4 × 10 5 n −1/2 e cm and τ e as λ Ω = δ 2 e τ e = 9.9 × 10 12 (Λ/10) This numerical expression allows us to estimate the Lundquist number S where the last expression holds in a hydrogen plasma. Under astrophysical conditions, S is enormous: for example, in ionized interstellar gas with T = 10 4 , n e = 1, B = 3 × 10 −6 G, and L B = 1 pc, S = 10 17 . With the exception of extremely dense, cold environments such as protostellar disks, or systems with extremely small lengthscales, such as the outer layers of accreting neutron stars, λ Ω is always much less than the numerical diffusivity arising from discretization. The viscous diffusivity, or kinematic viscosity ν, is The Reynolds number R with respect to an organized flow V on a lengthscale L is defined as LV /ν, and is almost always extremely large, confirming our expectation that astrophysical flow can be structured over a wide range of scales. According to eqns. (19) and (21), the ratio of viscous to magnetic diffusivity, which is known as the magnetic Prandtl number P r , is Typically, P r ≫ 1; in the example of fully ionized interstellar gas introduced above, P r = 10 11 . Therefore, the kinetic energy spectrum should be truncated on much larger scales than the magnetic fluctuation spectrum [55,66,24] The computations by Maron & Cowley [66] show that the essentially constant shear of the velocity field on scales near the resistive scale generates tightly folded, or hairpin-like, magnetic structures. If the ion cyclotron frequency ω ci ∼ 1.8 × 10 4 BZ/A and the ion collision time τ i satisfy the condition ω ci τ i ≫ 1, the viscosity is highly anisotropic with respect to the magnetic field. The viscous force acting on shear flow parallel to the magnetic field is reduced by a factor of (ω ci τ i ) 2 , while the parallel viscosity remains the same (see [13] for a full description). Using eqn. (18), we see that ω ci τ i is indeed large; in our numerical example, it is nearly 10 5 , implying very strong suppression of viscosity perpendicular to the magnetic field. Maron & Cowley [66] have implemented the full tensor viscosity in their simulations of turbulent dynamos.

Implications for Numerical Techniques
In order for a numerical method to include diffusive processes, we would ideally require two things, namely (a) that it can resolve the physically important scales and (b) that it can handle the disparate time scales. Clearly, meeting each of these requirements alone is already a non-trivial task.
Resolving the physically interesting scales means separating the intrinsic numerical diffusion existing in any scheme from its physical counterpart. Simultaneously, it is desirable to cover scales much larger than the diffusive scales as well, in order to embed the problem in a realistic environment. Currently, studies are restricted to one of these regimes, for obvious reasons.
There are various ways to parametrize the small-scale physics in the large scale approach. Shock capturing is a familiar example (for a detailed discussion see e.g. [58]). Many finite-difference codes apply a von-Neumann-Richtmyer artificial viscosity [78], spreading out the shock across several support points and not affecting the solution in smooth regions. Artificial viscosity mimicks physical viscosity in that it satisfies the Rankine-Hugoniot conditions, but on a much larger scale, which renders the numerical scheme very robust at low computational cost. It operates only where the flow is compressive.
Galsgaard & Nordlund [36] implemented an artificial resistivity algorithm similar to artifical viscosity in their simulations of formation and dissipation of current sheets in the solar corona; see also Caunt & Korpi [21].
"Current sheet capturing" is intrinsically more complex than shock capturing. As we discuss in §4.3, the energy dissipated by a shock is independent of the magnitude of the diffusivity, but this does not hold for a current sheet. Moreover, the topological evolution of the magnetic field is determined primarily by processes in current sheets, and this imposes additional requirements on how they are treated. Detailed, small scale studies could be useful for developing a parameterization of the effects of diffusion on energy balance and topology. This is reminiscent of the approach taken in large eddy simulations, which numerically resolve the largest scales while treating the unresolved scales with a subgrid model (e.g. [18,27,93]).
Godunov-type methods [41] avoid broadening the shocks over many zones. They solve the Riemann problem defined by a set of advection equations for the physical variables and the physical states on both sides of the discontinuity. The upstream and downstream states are connected by a sequence of shocks or rarefaction waves, each of which satisfies the Rankine-Hugoniot conditions. As implementation of a full Riemann solver would involve determining the full wave structure and all wave propagation speeds, a variety of approximate Riemann solvers has been designed to simplify the process by disregarding certain wave types or linearizing the problem (see [58] for details). Riemann solvers utilize the wave nature of hyperbolic equations, however, the diffusive terms are parabolic, so that they cannot be included in the Riemann solver directly.
An alternative approach has been taken in the so-called BGK-schemes [7], using the collisionless Boltzmann equation as a model of gas dynamics (see [89], and [104] for MHD).
Needless to say, the P r ≫ 1 regime cannot be captured by numerical codes which solve the ideal fluid equations (η = ν = 0), since the numerical diffusivities for magnetic field and velocity are about the same size and occur at about the same scale, i.e. the grid scale. This is an argument for including viscosity and resistivity in the equations, in which case the physical diffusivities must be larger than their numerical counterparts, and their scales must be well separated. This is the approach taken by Maron & Cowley [66], who used a spectral code to study turbulence with widely separated, but resolved, viscous and resistive scales.
Several efforts have been made to connect the spatial scales by adaptive mesh refinement techniques [5,4] for large eddy simulations [27], general hyperbolic systems [6], incompressible non-ideal MHD [34] and for compressible MHD [110,2]. However, turbulent flows by definition connect scales over the whole domain, which complicates finding a suitable refinement criterion and may in the end yield only modest savings of time.
We have seen in §3.1 that timescales in astrophysical problems can differ by several orders of magnitude. Including diffusive processes aggravates the problem, as their time scales are generally orders of magnitude longer than the dynamical times. Any disparity in timescales leads to a stiff problem. Explicit methods -often the first choice because of their low time and memory needsadvance the solution by a time step ∆t which is given by the so-called Courant -Friedrichs -Levy (CFL) condition [28] ∆t ≤ ∆x c (23) with the distance between two support points ∆x and a characteristic propagation speed c. Disparate characteristic speeds then lead to severe time step restrictions. Moreover, the slowly varying components may introduce numerical errors [50]. One solution to the problem is to choose a short timestep for the fast processes, and update the variables controlled by slow processes less frequently. Mac Low et al. [63] used this so-called subcycling in their implementation of ambipolar diffusion in the ZEUS code. Depending on the problem, the critical processes often can be broken out and treated implicitly, while the remaining equations can be treated explicitly [64]. Fully implicit schemes (e.g. [51]) can be unconditionally stable and independent of the time step chosen. Their main limitation lies in the tremendous computational needs in the multi-dimensional case. For a discussion of implicit methods and for a comparison between explicit and implicit shock capturing methods see [50,109]. For a discussion of higher order finite difference schemes and their application to MHD turbulence simulations see [14].
As this brief discussion shows, a great variety of techniques can be brought to bear on MHD turbulence problems. Which technique is best depends very much on the problem at hand.

Results from Theories and Simulations
Theory offers an interpretive framework for simulations, while numerical experiments test turbulence theory, and can inspire new developments. Here, we discuss amplification of a field by turbulence in the context of dynamo theory, magnetic fluctuations, the formation of current sheets, and the dynamical effects of turbulence on large scales. The first three of these problems are closely related. The fourth is topical in view of recent simulations of star forming regions.

Turbulent Amplification of a Weak Field
Amplification of a weak magnetic field by turbulence is one of the main components of dynamo theory. A successful astrophysical dynamo theory must demonstrate that a dominant, large scale magnetic field can be generated by small scale turbulence, on a timescale which is virtually independent of the Lundquist number S in the limit S → ∞. This is a highly nonlinear problem, which involves multiple scales. It probably cannot be solved without numerical simulations, but in view of the extreme parameters involved it is unlikely to yield to brute force alone.
The growth rate of magnetic energy in a periodic domain is given by eqn. (12). Although the resistive term appears as a sink, some resistivity is necessary to effect irreversible topological change. Amplification occurs through work done by the flow. In general, there is also a surface integral −c dS · E × B on the RHS of eqn. (12), which represents an energy flux through the boundaries.
Equation (12) says nothing about the structure of the field, and, in particular, whether magnetic energy is concentrated at large or small scales. If we formally set λ Ω ≡ 0, the ratio of fieldstrength to line length is preserved by incompressible motions. The only way to lengthen fieldlines in a finite volume is to wind or tangle them. Therefore, we expect that most of the energy will be in the rms field, rather than in the mean field. For example, if the protogalactic magnetic field were 10 −17 G [40], the fieldlines would have been lengthened by more than a factor of 10 11 in the course of amplifying the field to its present strength. Yet, the large scale and small scale components of the galactic magnetic field are observed to be roughly equal [113]. It is not enough to argue that magnetic forces will prevent the field from becoming tangled on small scales. If the lines aren't somehow lengthened, the field won't be strengthened.
Of course, we are interested in λ Ω → 0, not λ Ω ≡ 0. The fieldlines are allowed to break, and the question is whether they do so in a way which prevents power from piling up at the resistive scale and peaking instead at the large scale. Because S is so large, topological constraints are strong, and we expect resistive boundary layers to form.
Let's see how theories and simulations of dynamo theory address this problem. The most influential and detailed dynamo theory is the so-called mean field theory, which was proposed by Parker [83], and extended and formalized by Steenbeck, Krause & Rädler [97]. Mean field theory is based on the idea that a large scale magnetic field B can be generated from the average inductive electric field E associated with small scale, helical velocity and magnetic fluctuations δv, δB. The dynamo property of small scale, helical fluctuations is known as the α effect. There is also diffusive transport of B , which is known as the β effect. Although α and β are in general tensors, it suffices for our purposes to write them as scalars. The induction equation for the mean field is where In eqn. (25), V is the large scale velocity field, and the "..." represent additional terms involving moments of the turbulent fluctuations and successively higher derivatives of B (see Moffatt [74]). Note the correspondence between eqns. (24) and (25) and eqns. (1) and (2).
If α and β are independent of B , eqn (24) is linear in B . Physically, this corresponds to neglecting the back reaction of the magnetic field on the turbulent velocity, so the linear case is sometimes called the kinematic case.
The linear version of eqn. (24) has solutions which vary exponentially in time. At large shear and low wavenumber k, the growth rate is approximately the geometric mean of the shear rate V ′ and the turbulent frequency kα.
Standard mean field dynamo theory does not directly address the evolution of the small scale field, and the role of resistivity, on which the fate of the small scale field depends, is left somewhat implicit. Usually λ Ω is dropped from eqn. (24), since β is assumed to be much larger. A calculation by Moffatt [74] of the α effect due to helical Alfvén waves demonstrates that α ≡ 0 unless there is a phase difference between δv and δB. In Moffatt's model, the phase lag arises from magnetic diffusivity. At small λ Ω , α ∝ λ Ω , implying slow dynamo action at large R m .
In more general calculations of α, the diffusive damping time (k 2 λ Ω ) −1 is replaced by the correlation time τ of the turbulence. At large R m , this is much shorter than the resistive time, and implies large growth rates for the dynamo. But β does not distinguish between converting large scale field to small scale field and actually destroying field, so the behavior of the small scale fields are obscured by this treatment. Kulsrud & Anderson [55] used the tools of mean field theory itself to follow the growth of the small scale field. They showed that at large R m , the rms field B 2 1/2 grows much faster than the mean field B , and is dominated by small scale fluctuations with a k 3/2 power spectrum.
A completely different perspective on kinematic dynamo theory emerges from solutions of the full magnetic induction equation for prescribed, chaotic flow. These studies, which are fully described up to 1995 in the volume by Childress & Gilbert [23], use a combination of numerical simulations at large R m and analytical theory to establish connections between the properties of the flow and the properties of the magnetic field. Although there are known examples of flows which amplify the magnetic field at a rate independent of R m , as originally shown by Galloway & Proctor [35], the fields produced by these so-called fast dynamos are highly intermittent in space and fluctuate in sign, confirming at least qualitatively that small scale fields are dominant in kinematic dynamos at high R m .
The kinematic theories are modified by dynamical effects. It was expected on general grounds that the dynamo would saturate as the magnetic field approached equipartition with the kinetic energy. A saturation process known as "Alfvénization" was first described by (Pouquet et al [88]). Saturation occurs because the small scale fluctuations increasingly resemble Alfvén waves as the mean field grows. In an ideal medium δB → ±δv in Alfvén units, so their cross product also tends to zero, eliminating the α effect (also shown in the calculation by Moffatt). Alfvénization was first identified by Pouquet et al. [88] using a spectral closure scheme, and later derived using quasilinear theory by Gruzinov & Diamond [44]. In the model of Pouquet et al., kinetic helicity injected at the forcing scale generates magnetic helicity, which undergoes both a direct cascade, to high k, where it is resistively dissipated, and an inverse cascade, to low k. Alfvénization saturates the dynamo once the large scale field is roughly in equipartition with the kinetic energy. In the model of Gruzinov & Diamond, saturation by Alfvénization takes place when the ratio of kinetic to large scale magnetic energy is O(R m ), confirming the dominance of small scale fields at large R m . The difference between these calculations may lie in their treatments of diffusion. Gruzinov & Diamond explicitly include λ Ω , while Pouquet et al. took all the diffusivities to be eddy diffusivities, which are much larger than molecular diffusivities.
The role of kinetic helicity injection has been tested in simulations. Computations by Maron & Cowley [66], without kinetic helicity injection, show a peak in spectral power at the resistive scale, as predicted by Kulsrud & Anderson [55]. Simulations by Maron & Blackman [65] show that a large scale field is generated when kinetic helicity is injected, but that the growth rate is proportional to the resistive time. Nonlinear simulations with helical forcing by Brummell et al. [17] show that the magnetic lengthscale is proportional to R −1/2 m , while kinematic investigations testing the sensitivity of fast dynamo action to kinetic helicity have not found an inverse cascade [49].
Thus, we see that much hinges on the resistive time, or, more broadly, the mechanism by which small scale fields are dissipated. Turbulent diffusion of B , parameterized by the β effect, merely stands for spectral transfer out of the large scale. The possibility remains that turbulence could mix the field to the small scale, efficiently destroying it. However, numerical models by Cattaneo et al. [20], and subsequent analytical work by Kim [52] suggest that the mixing is self limiting due to Lorentz forces.
The possibility that the small scale fields might escape through an open boundary rather than being destroyed in situ, leading to growth of the large scale field on the convective timescale rather than the resistive timescale, was raised by Parker [86] and Blackman & Field [11] Up to now, neither analytical studies [92] nor numerical models [15] have confirmed this idea.
Astrophysical dynamo theory is at a critical juncture. The fundamental premise of mean field dynamo theory -that kinetic helicity injected at small scales can drive the growth of a magnetic field at large scales -is supported by a variety of analytical treatments and numerical investigations. Mean field theory is a felicitous outcome of turbulence theory in the sense that it can be applied just by calculating a few parameters (or functions): α, β, and the large scale shear. However, neither theory nor simulations have fully come to grips with the large value of R m and its implications for the small scale fields. From a numerical point of view, computations at large R m which can capture current sheets as well as accommodating a variety of boundary conditions would seem to be a prerequisite for either validating or superseding mean field theory.

Magnetic Fluctuations
Steady hydrodynamic turbulence, uninfluenced by boundaries, has been characterized in two regimes. Subsonic turbulence is incompressible, and follows the Kolmogorov scaling (E(k) ∝ k −5/3 ). Supersonic hydrodynamic turbulence is compressible, and better described as Burgers turbulence (E(k) ∝ k −2 ). (These power law spectra, and all subsequent ones, describe only the so-called inertial range, in which there is no driving or microscopic dissipation, only nonlinear spectral transfer.) Both these limiting forms of turbulence are spatially intermittent in the sense that the variances of quantities such as the rate of energy dissipation undergo substantial fluctuations from point to point. In Burgers turbulence the basic structures are shock waves, and the large fluctuations and gradients occur principally in the shock fronts. In Kolmogorov turbulence, the basic units are eddies, and intermittency is represented by the concentration of vorticity into small structures. Intermittency does not affect the energy spectrum of Burgers turbulence (the Fourier spectrum of a shock itself is k −2 ), but is thought to steepen the Kolmogorov spectrum. This occurs through the local enhancement of the dissipation rate at small scales [57].
Hydrodynamic turbulence in unstratified, nonrotating systems is isotropic. A uniform component of magnetic field causes turbulence to be anisotropic.
An ideal, adiabatic fluid with a uniform component of magnetic field has three linear modes: the shear Alfvén mode, and the fast and slow magnetosonic modes. The velocity field δv of the Alfvén mode, is solenoidal, and this mode is driven purely by magnetic tension. The dispersion relation is ω = k v A ≡ ω A , where k is the component of k parallel to the large scale magnetic field. The fast and slow magnetosonic modes are compressive, and driven by a combination of gas pressure, magnetic pressure, and magnetic tension.
Since the shear Alfvén mode is noncompressive, it is reasonable to ask whether a strong (β < 1) magnetic field permits a regime of pure shear, supersonic turbulence. The answer is no.
First, suppose k ⊥ ≡ 0. Then, the nonlinear magnetic pressure gradient associated with the fluctuating transverse magnetic field drives a compressive parallel flow, causing an Alfvén wavetrain to steepen nonlinearly into a train of weak shocks [26], much as a sound wave steepens. The only exception is the infinitely long, circularly polarized wave, which is an exact solution of the ideal MHD equations. The outcome of the evolution of an ensemble of parallel propagating Alfvén waves is a series of shocks, similar to Burgers turbulence, with k −2 kinetic and magnetic energy spectra [39]. The steepening time depends on amplitude as (δv/v A ) −2 , while the steepening time of an acoustic wave is proportional to (δv/v S ) −1 . In this sense, shear waves survive longer in a low β plasma than compressive waves of the same amplitude, by a factor of order β −1/2 (v A /δv), but the outcome in either case is a compressive flow dominated by shocks, and the difference in timescales is only significant for small wave amplitude.
Although Burgers turbulence was originally based on 1D flow, 3D simulations of supersonic MHD turbulence show similar rapid evolution to a shock dominated flow. The distribution of shock strengths in simulations of driven and decaying turbulence, with and without magnetic fields, has been studied by Smith et al. [95,94].
The steepening of Alfvén waves, which can be viewed as a cascade in k , is suppressed if β ≥ 1. Under the assumption of incompressibility, and in contrast to the highly compressible low β case, interactions occur only between oppositely directed wave packets. In particular there is no self interaction of the kind which leads to steepening.
The incompressible shear Alfvén wave cascade is highly anisotropic, with power transferred much more rapidly in k ⊥ than in k . In physical space, the correlation length transverse to the mean field is much shorter than the correlation length along it. These anisotropic states can be described by so-called reduced MHD [101], which is a quasi-2D approximation describing elongated magnetic structures. It is often applied to (and, indeed, was derived for), low β plasmas, but it neglects the parallel steepening which leads to shocks. This, as we said above, is a good approximation for long parallel Alfvén transit times and small turbulent amplitudes.
There is not yet complete consensus on the cascade itself. Weak MHD turbulence, in which an individual Alfvén wave packet survives for many periods, can be viewed as weakly perturbative resonant interactions between multiple waves. Sridhar & Goldreich [96] argued that the dominant interaction is a four wave interaction, which results in a spectrum E(k) ∝ k −7/3 . Others [79,8,38] claim that the dominant interactions are three wave interactions, and predict In strong MHD turbulence, wave packets survive for only about one wave period. Goldreich & Sridhar [43] developed the concept of a critically balanced cascade in which the wave frequency ω A is the same as the nonlinear frequency k ⊥ v ⊥ . This, together with the requirement of constant spectral energy flux argument k ⊥ v 3 ⊥ , leads to a Kolmogorov spectrum E(k) ∝ k −5/3 and wavenumber anisotropy k /k ⊥ ∝ k −1/3 ⊥ . Other theories of strong MHD turbulence predict E(k) ∝ k −3/2 , including the original isotropic Iroshnikov-Kraichnan theory [77].
Numerical simulations, if free of confounding computational effects, could play a role in resolving these disagreements. Biskamp & Müller [10] and Cho & Vishniac [25] find E(k) ∝ k −5/3 over about one decade in k space. Maron & Goldreich [67] slightly extended the inertial range and found E(k) ∝ k −3/2 . These authors argue that the spectrum is flattened because of intermittency. The crucial difference between the HD and MHD cases is that energy cascades in MHD only when oppositely directed wave packets collide; the small filling factor associated with intermittency reduces the collision rate, more than offsetting the enhanced dissipation associated with small structures.
Larger computations, with an extended inertial range, should shortly become available. The differences in the spectra obtained with different codes and under different forms of driving (Biskamp & Müller studied decaying or isotropically forced turbulence, Cho & Vishniac studied isotropically forced turbulence, and Maron & Goldreich studied anisotropically forced turbulence) are at this point comparable to the theoretical disagreements. A joint exercise in which different groups simulated identical problems and implemented identical diagnostics might be quite enlightening.
A variety of processes can terminate the cascade at short wavelengths (see [59] for a recent discussion). These include fluid effects, such as ion-neutral friction, and kinetic effects such as ion gyroresonance absorption and electron Landau damping [91]. The extension of the magnetic spectrum to scales below the viscous cutoff in a high Prandtl number plasma is discussed by [66] and [24].

Intermittency: Current Sheets and Flux Tubes
Simulations of MHD turbulence show strong intermittency in the distribution of magnetic field gradients, or current. In the weak field, high Prandtl number simulations of [66], the fluctuating field is tightly folded. These results suggest that we look at current sheets.
Current sheets and filaments are sites of intense local heating, and possibly particle acceleration, which makes them interesting in their own right. It is quite likely that in most astrophysical systems, significant departures from flux freezing occur only in these regions, so their existence is important in the evolution of magnetic field topology, and for dynamo action, (see §4.1). Reconnection at magnetic X-points drives strong, small scale jets. Simulations by Galsgaard & Nordlund [37] show highly time variable dissipation rates associated with the formation and disruption of large scale current sheets.
Current sheets are a prerequisite of, and accompany, magnetic reconnection. Reconnection is a resistive process which occurs on a timescale intermediate between the Ohmic time and the Alfvén time. The reconnection time τ rec is often parameterized by the Lundquist number S; τ rec scales as S p with 0 ≤ p < 1. If p = 0, the reconnection is said to be fast (in correspondence with the definition of a fast dynamo). In Sweet-Parker reconnection, p = 1/2 [85], and resistive tearing modes in slabs and cylinders have p = 3/5, and p = 1/3, respectively [3]. The formation and behavior of current sheets has been studied for many years. Most of this work has focussed on plasmas in or near magnetostatic equilibrium, which can be expected if β ≪ 1 and the plasma is either allowed to relax or is driven at frequencies far below its characteristic Alfvén frequency. Such conditions hold, for example, in stellar coronae, in which the magnetic fieldlines are tied to the underlying photosphere, and evolve quasistatically in response to slow photospheric motion.
Current sheets form as a plasma seeks equilibrium while obeying the topological constraints imposed by flux freezing [75]. Topological constraints can arise from features such as null points [102] or separatrices [62,115,60]. Sweet [103] and Parker [84] suggested that shearing and winding the footpoints of an initially uniform field will create current sheets. Even a simple sinusoidal deformation of the boundaries of a sheared magnetic slab can induce the formation of a current sheet [45].
As a simple example of how current sheets can form, consider the magnetic field of a pair of neighboring, uniformly magnetized superconducting spheres, with antiparallel magnetic dipole moments, surrounded by a perfectly conducting, zero pressure plasma. Some fieldlines have both endpoints on the same sphere, and others have one endpoint on each sphere. Now suppose that one sphere rotates on its magnetic axis by an angle φ. The fieldlines which are not connected to the other sphere also rotate by φ, but the fieldlines which are connected to the other sphere cannot rotate uniformly because one endpoint remains fixed. The separatrix surface which divides the two domains of magnetic connectivity becomes a current sheet: the field outside it is sheared while the field inside it is not.
There is still no comprehensive theory for exactly when current sheets form, and how their properties reflect the underlying magnetic topology. Nevertheless, analytical arguments [108,61] and numerical experiments [73,61,36] show that shear, or equivalently current density, grows exponentially rapidly in line tied magnetic fields with randomly moving footpoints. This has some of the same consequences as an exact singularity, which cannot form in any case in even a slightly resistive plasma.
These results are not fully applicable to MHD turbulence, in which the magnetic field is not in equilibrium. Nevertheless, current sheets form in flows with magnetic null points [87]. The underlying mechanism for the growth of the current is the random stretching of neighboring fieldlines, which leads to large crossfield shear.
What are the consequences of current sheet formation for dissipation of magnetic energy? That is, if energy per unit volumeė is being added to a system by a driving process, under what circumstances can this energy can be dissipated in current sheets?
Consider a current sheet of area L 2 and width δ. Assume the current sheet structure is given by the Sweet-Parker theory, so that δ = LS −1/2 = (Lλ Ω /v A ) 1/2 . If N (L)dL is the number density of current sheets with L between L and L + dL, then the rate at which energy is dissipated by current sheets in this size range is The rate of energy dissipation by all current sheets is obtained by integrating eqn. (26) over L:ẇ Ifė is independent of λ Ω as λ Ω → 0, then dLN (L)L 3/2 must scale as λ −1/2 Ω in this limit. This would imply an unbounded number of current sheets as λ Ω → 0, which is virtually equivalent to a pileup of magnetic energy in fluctuations at the resistive scale. 1 There are alternative possibilities. One is that the current sheets are not structured according to the Sweet-Parker picture. If δ scaled as λ Ω , as it does (up to a logarithmic factor) in Petschek's model of fast reconnection (reviewed in [85]), then dLN (L)L 3/2 could be independent of λ Ω . Note, however, that while in the Sweet-Parker theory, magnetic energy is converted to kinetic energy and to heat at about the same rate, in Petschek's theory magnetic energy is converted predominantly to kinetic energy, and that the realizability of Petschek's model has recently been questioned [9,107].
The second possibility is that the magnetic and velocity fields adjust themselves such that the driving just balances the damping, as they do in mechanical systems such as linear harmonic oscillators. Thus, the energy input rateė would depend on λ Ω , while the applied force itself, presumably would not.
Dedicated simulations of current sheets and reconnection regions themselves allow us to study the small scale aspects of the problem. It is equally important to understand how current sheets are embedded in the overall flow. This has already been done to some extent for laboratory experiments [107]. Ultimately, it may be possible to parameterize magnetic reconnection in simulations on larger scales.
Finally, we mention an even more extreme form of intermittency: the segregation of the magnetic field into thin tubes. This is observed in the solar photosphere, where magnetic flux tubes, or sheets, just a few hundred km thick are seen at the borders of convective cells. Flux tubes can form only in high β plasma, so that gas pressure can balance magnetic pressure at the tube walls.
Early attempts to explain photospheric flux tubes relied on a two stage process, proposed by Parker, for concentrating a diffuse field into organized structures. According to this scenario, thermal convection would sweep the field to the borders of convective cells, concentrating the field up to equipartition with the flow. Then reduced convective heat transport would cool the gas, "collapsing" the tube and bringing it to equipartition with the thermal plasma.
An alternative possibility is that intermittent fields are generated in situ. It has been shown that small scale turbulent thermal convection operates as a dynamo, and generates a magnetic field with a broad tail of high energy features [19].
It is unclear whether a dynamo operating in a high β plasma should always produce flux tubes. In the solar case, the turbulence is organized by stratification and the thermal gradient, and is subsonic. If the turbulence were supersonic and the field were amplified to equipartition, then flux tubes could not persist.

Dynamical Effects of Turbulence
The stresses associated with turbulence can exert dynamical forces. In MHD turbulence, these forces include both magnetic stresses and the Reynolds stress. Observations suggest that these stresses are important in the interstellar medium, and especially in molecular clouds.
It was realized early that supersonic motions are nearly ubiquitous in molecular gas, and pose a problem. If the motions are associated with gravitational collapse, the implied star formation rate would be very high [111], but if the motions are turbulent, they should be dissipated very quickly [42]. This conundrum, together with the expectation that reasonably strong magnetic fields should be present, prompted analytical studies of Alfvénic turbulence under molecular cloud conditions [1,114,32]. This work demonstrated that even when M A ≪ M and the waves are purely transverse, ion-neutral friction and nonlinear steepening are strong damping mechanisms, which limit the lifetimes of Alfvén waves to a few 10 6 yr or less under molecular cloud conditions. Even so, under the assumption that the waves might be replenished by some energy source, various aspects of the theory of self gravitating clouds supported by Alfvén wave stresses were developed and applied to model clouds [90,33,69,68,70]. According to this theory, which is based on the the weak turbulence approach pioneered by Dewar [29], Alfvén waves have an isotropic stress tensor P = P w I. The waves pressure -density relation is P w ∝ ρ 1/2 in a stratified, static medium, while P w ∝ ρ 3/2 under slow, spatially uniform changes in density with time. The negative polytropic index (n = −1/2) of the P −ρ relation leads to a large center to surface pressure contrast, in accord with observations (see [70] for a good discussion).
The theory of wave supported clouds is most relevant to objects which are slightly magnetically supercritical, that is, the ratio of mass to magnetic flux is slightly too large for magnetostatic support. If the clouds were subcritical they would be supported by the DC magnetic field. If they were highly supercritical, the amplitude of the turbulence required to sustain them would be so large that the weak turbulence theory would probably not apply, although they could still be turbulently supported. ¿From an observational viewpoint, the most likely candidates for wave support are dense cores with sonic or mildly supersonic linewidths [76].
The picture of strongly supersonic turbulence which emerges from numerical simulations of molecular clouds is quite different from this analytical picture, primarily because, as we discussed in §4.2, generation of compressive disturbances is entirely unavoidable. These compressive disturbances, or shocks, sweep up most of the mass into thin layers in which the local Jeans length is small. These layers, and the even denser structures formed where layers intersect, collapse if the domain is magnetically supercritical, whether the turbulence is freely decaying or maintained in a steady state by driving [47,81]. The driven models are globally stable in the sense that λ J computed using the total turbulent energy and mean density is larger than the length of the domain. The small, high density structures collapse because the turbulent power at wavelengths less than the local Jeans length is relatively small (see eqn. (17)), and because the intensity of turbulence inside the cores is comparable with the density outside. Nevertheless, the bound structures are small, and the role of numerical diffusion within them, which damps the turbulence and removes the support by the DC magnetic field, has not been quantified. Therefore, while the formation of turbulently supported clumps has not been found in the numerical models, their existence cannot be ruled out based on the present models.

Ambipolar Drift
As we have discussed in the preceding sections, both theoretical and computational aspects of astrophysical magnetohydrodynamics are dogged persistently by the enormous value of R m . There is one situation, however, in which the magnetic diffusivity is large whether or not it is enhanced by turbulence: in weakly ionized gas, the magnetic field and plasma drift with respect to the neutrals. This so-called ambipolar drift is not entirely equivalent to resistive diffusion, because it preserves magnetic topology (the field is frozen to the plasma), but it does change the ratio of magnetic flux to mass, as originally pointed out in [72], and dissipates fluctuations on small scales.
Full treatment of ambipolar drift requires solving the equations of MHD for the plasma and neutrals treated as distinct species coupled by collisions, ionization, and recombination [30]. At the low ionization fractions expected in molecular clouds, the ion Alfvén speed v Ai ≡ B/ √ 4πρ i can be several orders of magnitude larger than the other characteristic speeds in the problem, reducing the maximum possible timestep by a similar factor. This and other numerical issues related to implementation of ambipolar drift are discussed by [63,105,64,98].
On timescales longer than the ion -neutral collision time τ in , and at low ionized mass fractions, the ion-neutral drift v D is determined by balancing the Lorentz force on the ions against the frictional drag by neutrals wheref is a unit vector in the direction of the Lorentz force. The second equality in eqn. (28) can be taken as a definition of the magnetic lengthscale L B . The quantity v 2 Ai τ in has units of diffusivity, and from now on we refer to it as λ AD . With this definition, v D = λ AD /L B . Note that ρ i /τ in = ρ n /τ ni and ρ n ≈ ρ, so v 2 Ai τ in = v 2 A τ ni , and λ AD is often expressed in the latter terms. The ambipolar Reynolds number R AD is defined as the ratio of the bulk velocity v to the drift velocity v D and measures how well the magnetic field is frozen to the bulk fluid on scale L B [112]. Setting R AD = 1 correctly predicts the thickness of shocks in which the main dissipation mechanism is ion-neutral friction [31], the wavelength at which Alfvén waves are critically damped [56], and the minimum size of an eddy which can wind up a magnetic field [112].
Simulations of supersonic MHD turbulence [82] with varying strengths of ambipolar drift shows that, on average, the drift smooths out the current, and reduces the rms Lorentz force.
We expect (with one caveat discussed at the end of this section) that the magnetic field should show very little structure below the scale L min ≡ λ AD /v at which R AD = 1. Consider a volume L 3 of turbulent gas with mean density ρ, turbulent velocity v, and mean magnetic field B. If we introduce the crossing time τ d ≡ L/v and the Alfvén Mach number M A ≡ v/v A , then we can write Expressing v in units of km s −1 as v 5 , and L in pc as L pc , τ d ∼ 3×10 13 L pc /v 5 s. For the particle species expected in molecular clouds, τ ni ∼ 5 × 10 8 n −1 i s. With the ionization equilibrium relation n i ∼ 10 −5 n Comparing eqn. (31) with the grid spacing ∆x/L = 1/N , we see that a moderately large numerical simulation (say 256 3 ), which is not highly super-Alfvénic or extremely dense, should resolve almost all the magnetic structure associated with turbulence. Equation (30) can also be written in terms of the critical magnetic field B c ≡ 2πG 1/2 ρL and the free fall time τ f f ≡ (4πGρ) −1/2 ∼ 5.5 × 10 14 n −1/2 n s as where the last step holds for virial equilibrium; τ f f ∼ τ d . Of course, eqn. (30) is entirely independent of G.
The modest value of R AD in molecular clouds has been commented upon and exploited by [54,80]. Ambipolar drift is essentially diffusive when the magnetic field has a well ordered, well combed structure. Because the ambipolar diffusivity is proportional to B 2 , the diffusion is nonlinear. Sharp fronts can be generated in the vicinity of magnetic nulls, and minima are steepened [16,63,112,64].

Summary and Future Agenda
The theme of this paper is how numerical simulations can best cope with the enormous magnetic Reynolds numbers encountered in most astrophysical problems. We argued that, because R m (and the Lundquist number S) are so large, resistive effects occur primarily in thin sheets or filaments, which must be treated accurately in order to follow the topological evolution of the magnetic field.
In §2, we reviewed the magnetic induction equation, giving the solution in the ideal limit, and derived equations for magnetic helicity and magnetic energy. We tested helicity conservation in the ZEUS code by simulating the evolution of the ideal kink mode, and studied fluctuations of helicity in a simulation of a turbulent cloud. The first of these problems is especially suitable for benchmarking numerical codes. That is, it yields a quantitative estimate of the numerical resistivity in the code, from which one can compute the resistive timescale (at least in smooth regions, on large scales) and assess whether it is much longer than other timescales of interest. One can also compare the performance of different codes.
In §3, we discussed some important dimensionless parameters in astrophysical MHD problems. We verified the large size of the Lundquist number S and magnetic Reynolds number R m under most conditions, and showed that the magnetic Prandtl number, the ratio of viscous to magnetic diffusivity, is also usually large. This implies that magnetic structure will generally extend to smaller scales than velocity structure. We briefly summarized a variety of numerical techniques from the perspective of the large S limit.
In §4, we discussed theoretical and numerical results on several basic problems. In §4.1 we discussed the amplification of magnetic fields by turbulence, and argued that the successful operation of a large scale dynamo depends on fast diffusion of the magnetic field. In §4.2 we discussed MHD turbulence itself. In §4.3 we discussed current sheets, which may provide the fast diffusion necessary for a dynamo, and may be produced in turbulence as a form of intermittency. In §4.4, we discussed the dynamical effects of turbulence, and its importance in molecular clouds.
In §5, we discussed ambipolar drift as a mechanism for increasing the magnetic diffusivity in low density, weakly ionized gas. It is now possible to simulate turbulent molecular clouds at a resolution which accounts for all magnetic structure down to the ambipolar diffusion scale. Under certain conditions, however, ambipolar drift can mediate the formation of current sheets, possibly leading to rapid magnetic reconnection.
It is fortunate that a number of different groups are simulating MHD turbulence under astrophysical conditions, using a variety of codes and techniques. There is substantial agreement on a number of issues, such as the short decay time of supersonic turbulence, but disagreement on others, such as the spectrum of strong MHD turbulence. A dedicated effort at benchmarking, in which different groups simulated identical problems and implemented identical diagnostics, would provide some perspective. We suggested in §4.2 that such an exercise be carried out for strong MHD turbulence.
It is unlikely that it will ever be possible to adequately resolve current sheets and global dynamics in a single simulation. Local studies of current sheets and magnetic reconnection are essential in capturing the physics of these layers. It is equally important to understand how they are embedded in the overall flow, as this constrains their properties, sets the boundary conditions for reconnection, and is necessary for understanding how to parameterize the effects of current sheets in global simulations. The inability to deal adequately with this small scale structure is the greatest present challenge to understanding magnetic field evolution in turbulent flows.