A model of nonlinear evolution and saturation of the turbulent MHD dynamo

The growth and saturation of magnetic field in conducting turbulent media with large magnetic Prandtl numbers are investigated. This regime is very common in low-density hot astrophysical plasmas. During the early (kinematic) stage, weak magnetic fluctuations grow exponentially and concentrate at the resistive scale, which lies far below the hydrodynamic viscous scale. The evolution becomes nonlinear when the magnetic energy is comparable to the kinetic energy of the viscous-scale eddies. A physical picture of the ensuing nonlinear evolution of the MHD dynamo is proposed. Phenomenological considerations are supplemented with a simple Fokker--Planck model of the nonlinear evolution of the magnetic-energy spectrum. It is found that, while the shift of the bulk of the magnetic energy from the subviscous scales to the velocity scales may be possible, it occurs very slowly -- at the resistive, rather than dynamical, time scale (for galaxies, this means that generation of large-scale magnetic fields cannot be explained by this mechanism). The role of Alfvenic motions and the implications for the fully developed isotropic MHD turbulence are discussed.


Introduction
It has long been appreciated [1] that a generic three-dimensional turbulent flow of a conducting fluid is likely to be a dynamo, i.e., it amplifies an initial weak seed magnetic field provided the magnetic Reynolds number is above a certain instability threshold (usually between 10 and 100). The amplification is exponentially fast and occurs due to stretching of the magnetic field lines by the random velocity shear associated with the turbulent eddies. The same mechanism leads to an equally rapid decrease of the field's spatial coherence scale, as the stretching and folding by the ambient flow brings the field lines with oppositely directed fields ever closer together [2,3,4].
If the magnetic Prandtl number (P r = ν/η, the ratio of the viscosity and the magnetic diffusivity of the fluid) is large and the hydrodynamic turbulence has Kolmogorov form, the scale of the magnetic field can decrease by a factor of P r 1/2 below the viscous scale of the fluid. This situation is known as the Batchelor regime [5] and is characterized by the magnetic-energy concentration at subviscous scales. This MHD regime is very common in astrophysical plasmas, which tend to have very low density and high temperature. Examples include warm interstellar medium, some accretion discs, jets, protogalaxies, intracluster gas in galaxy clusters, early Universe, etc. (a representative set of recent references is [6,7,8,9,10,11,12,13,14,15]). The object of our principal interest are magnetic fields of galaxies, often thought to have been generated by the turbulent dynamo in the interstellar medium. The challenge is to reconcile the preponderance of small-scale magnetic fluctuations resulting from the weak-field (kinematic) stage of the dynamo with the observed galactic fields coherent at scales of approximately 1 kpc, or 10 times larger than the outer scale of the interstellar turbulence, which is forced by supernova explosions at scales of ∼ 100 pc ( [16,17,18] are some of the recent reviews of relevant observations).
The galactic magnetic fields generally have energies comparable to the energy of the turbulent motions of the interstellar medium, and, therefore, cannot be considered weak. Indeed, the growth of the magnetic energy established for the kinematic regime inevitably leads to the field becoming strong enough to resist further amplification and to modify the properties of the ambient turbulence by exerting a force (Lorentz tension) on the fluid. This marks the transition from the kinematic (linear) stage of the dynamo to the nonlinear regime. In the astrophysical context, the main issue is what happens to the coherence scale of the field during the nonlinear stage and, specifically, whether a mechanism could be envisioned that would effectively transfer the magnetic energy from small (subviscous) to velocity scales: perhaps as large as the outer scale of the turbulence and above.
A large amount of work has been devoted to this issue during the last 50 years (a small subset of the relevant references is [1,19,20,21,22,23,24,25,26]). Since most phenomenological theories of the steady-state MHD turbulence envision a state of eventual detailed scale-by-scale equipartition between magnetic and velocity energies [27,28,29], it has been widely expected that the dynamo would converge to such a state via some form of inverse cascade of magnetic energy. The implication would be an eventual dominant magnetic-energy concentration at the outer scale of the turbulencea convenient starting point for the operation of the helical mean-field dynamo [30,31] or of the helical inverse cascade [32,19]. Note that we believe the inverse cascade of magnetic helicity that takes place in helical MHD turbulence [32,19] to be unrelated to the small-scale processes we discuss here: we are interested in ways to shift the magnetic energy from subviscous scales to the outer (forcing) scale of the turbulence, while the helical inverse cascade transfers magnetic helicity from the forcing scale to even larger scales. In the galactic dynamo theory, helicity-related effects should not be relevant for the small-scale turbulence because they operate at the time scale associated with the overall rotation of the galaxy, which is 10 times larger than the turnover time of the largest turbulent eddies (10 8 and 10 7 years, respectively [6]).
To prevent confusion, let us spell out that by turbulent dynamo we mean the turbulent amplification of the energy of magnetic fluctuations, not of the mean field. Fluctuation-dynamo theories describe magnetic fields at the scales of the turbulence, mean-field-dynamo theories refer to much larger scales (the mean field is usually the field averaged over the flow scales). While fast mean-field dynamo is not possible in an isotropic system without net helicity (see, e.g., [30]), the fluctuation dynamo is [1,35,6]. Moreover, fluctuation dynamos in three-dimensional turbulent flows are usually fast, i.e., proceed at a finite rate in the limit of arbitrarily small magnetic diffusivity.
Thus, in this work, we study the possibility of the nonhelical, nonlinear, large-P r fluctuation dynamo saturating with magnetic energy concentrated at the outer scale of the turbulence. Using a simple spectral Fokker-Planck model inspired by the quantitative theory of the kinematic stage of the dynamo and by a physical model of the nonlinear stage, we establish a possible mechanism for shifting the bulk of the magnetic energy from the subviscous scales to the velocity scales. The energy transfer mechanism proposed here involves selective resistive decay of small-scale fields combined with continued amplification of the field at larger scales. Therefore, the time needed for this process to complete itself turns out to be the resistive time associated with the velocity scales of the turbulence, not a dynamical time. In large-P r systems, this time is enormously long and, in fact, can easily exceed the age of the Universe. Thus, a steady state with magnetic energy at velocity scales, even if attainable in principle, is not relevant for such systems. However, the scenario of the evolution and saturation of the nonlinear dynamo proposed here has fundamental implications for our understanding of the structure of the forced isotropic MHD turbulence and of the results of the recent numerical experiments.
The plan of the rest of this paper is as follows. In section 2, we review the kinematic theory. It is worked out in some detail because it constitutes the mathematical framework of the subsequent nonlinear extension presented in section 4. In section 3, a phenomenological theory of the nonlinear dynamo is proposed. This section contains our main physical argument and is, therefore, the part of this paper that a reader not interested in details should peruse before jumping to the discussion section at the end. In section 4 we proceed to formulate a Fokker-Planck model of the evolution of the magneticenergy spectrum. The model is then solved numerically and analytically. Self-similar solutions are obtained. Both the qualitative physics of section 3 and the analytical solutions of section 4 exhibit a resistively controlled effective transfer of the magnetic energy to the velocity scales. The resulting long-time asymptotic state of the fully developed isotropic MHD turbulence is discussed in section 3.4. Section 5 contains a discussion of our findings and of the issues left open in our approach.

The growth of magnetic energy
The kinematic dynamo is described by the induction equation where B is the magnetic field, η is the magnetic diffusivity, and u is the random incompressible velocity field, whose statistics are externally prescribed. The following exact evolution law for the magnetic energy W = 1 2 B 2 is an immediate corollary of (1): where, by definition, and M (t, k) = 1 2 dΩ k |B(t, k)| 2 is the k-shell-integrated magnetic-energy spectrum. The angle brackets signify ensemble averages, but, the system being homogeneous, can also be thought of as volume integrals. Batchelor [1] argued that, in a turbulent flow, γ(t) would tend to some constantγ, thus implying exponential growth of the magnetic energy. This is essentially a consequence of the line-stretching property of turbulence [33,34].
The growth of the magnetic energy is accompanied by the decrease of the characteristic scale of the magnetic fluctuations, which also proceeds exponentially fast in time at a rate ∼γ [1]. In a large-P r medium, this process quickly transfers the bulk of the magnetic energy to scales below the hydrodynamic viscous cutoff. In the Kolmogorov turbulence, the fastest eddies are the viscous-scale ones, so it is their turnover time that determines the time scale at which the small-scale kinematic dynamo operates:γ ∼ k ν u ν , where k ν and u ν are the characteristic wave number and velocity of these eddies [6].

The magnetic-energy spectrum
Since the subviscous-scale magnetic fields "see" the ambient flow as a random linear shear, one might conjecture that their statistics should not be very sensitive to the specific structure of the flow, provided the latter is sufficiently random in space and time. This gives a physical rationale for modelling the ambient turbulence with some synthetic random velocity field whose statistics would be simple enough to render the problem analytically solvable. The most eminently successful such scheme has been the Kazantsev-Kraichnan model of passive advection [35,36] which uses a Gaussian white-in-time velocity as a stand-in for the real turbulence: On the basis of this model, it is possible to show that the magnetic-energy spectrum at subviscous scales (k ≫ k ν ) satisfies the following equation [35,37,38,6,39,8] , which is consistent with its physical meaning as the characteristic eddy-turnover rate of the turbulence [40,8]. Note that, in order for equation (5) to be an adequate description of the spectral properties of the magnetic field, it must be a conservative approximation, i.e., be consistent with the exact energyevolution law (2). Indeed, we can rewrite equation (5) in a conservative Fokker-Planck form where D(k) =γk 2 /5 is the diffusion coefficient and V (k) = 4γk/5 the "drift velocity" in k space (V > 0 implies the spectrum moving towards higher k). We see that Remarks on the validity of the Fokker-Planck model. Equation (5) is valid in the limit that the flow velocity has a zero correlation time as given by (4), which implies an infinite kinetic-energy density. In order to calculate the relevant parameters of this model for comparison with the more realistic situation where the flow has a finite correlation time, one may take the following steps. In (4) where τ c is the correlation time of the velocity. Doing a calculation parallel to that in [40], one gets where E(k) is the energy spectrum of the fluid turbulence. [For a more refined model, one might allow for the decorrelation rate τ −1 c to vary with k and remain inside the integral in (7), or use a three-mode decorrelation rate instead to account for a finite ratio of the velocity and magnetic-field correlation times. But as long as one uses a value for τ c characteristic of the fastest (highest-k) eddies, which dominate the integral in (7), and τ c does not become large compared to the magnetic-field correlation time, the main scaling is the same.] Note that the integral in (7) is proportional to the square of the characteristic eddy-turnover rate τ −1 eddy of the flow. It is clear that, in a fully developed turbulence, τ c ∼ τ eddy , so we may rewrite (7) as where c 1 is a coefficient of order unity.
The other main assumption in (5) is that we are focusing on magnetic excitations at k ≫ k ν , which raises the question of the appropriate boundary conditions or modifications to the theory at low k. The effect the large-scale structure of the velocity field may have on the statistics of the small-scale magnetic field is incompletely understood. However, when the bulk of the magnetic energy is at the subviscous scales, it is probably safe to assume that the spectrum is not significantly affected by the details of what happens at large scales (see discussion in [8]). For P r ≫ 1, we shall find that k rms ≫ k ν , except perhaps at very large t. In this regime, our results are insensitive to the precise form of the boundary condition at small k, and we can simply specify a suitably low boundary k * ≪ k rms across which there is no flux of magnetic energy.
At k k ν , there is no scale separation between the velocity and magnetic field, the local (in k space) approximation (5) breaks down, and one must, strictly speaking, use a more general integrodifferential equation for M (t, k) [35,6,8], but simplifications may be possible. In the Kolmogorov picture of isotropic turbulence, the dominant stretching of eddies at a particular scale is due to interactions with other eddies at comparable scales. This suggests that the effect of small-scale velocity fields on larger-scale magnetic fields can be neglected, i.e., that the "turbulent diffusivity" acting on the magnetic field at some k < k ν is dominated by the contribution from the velocities at the same k.
As an approximation, one would then just continue to use the Fokker-Planck equation (6) but with the expression (8) forγ modified to account for the fact that velocities at a given scale are not very effective in causing net stretching of magnetic fields at larger scales: For k < k ν , the dominant contribution to this integral comes from k ′ ∼ k, i.e., from eddies of sizes comparable to the scale of the field that is being stretched. Note thatγ(k) in (9) has the same scaling as the k-dependent eddy-damping rate in the popular phenomenological EDQNM closure [41], where setting c 1 = 0.36 results in a good fit to the Kolmogorov constant [19]. The shearing rateγ(k) drops to zero for k below the forcing wavenumber k f of the turbulence, so the "infrared" (IR) cutoff k * just needs to be low enough thatγ(k * ) [and thus the flux of magnetic energy in k space, the expression in square brackets in equation (6)] also vanishes. In the analytic solutions presented below, we shall neglect any k dependence ofγ for simplicity and write the results in terms of a general IR cutoff k * . One could set k * ∼ k f as a crude model of the argument above, but the precise value of k * will not affect any of our results until section 4.5, where the choice and the meaning of k * will be revisited.

The diffusion-free (ideal-MHD) regime
Let us discuss how a magnetic excitation initially concentrated at some wave numbers k ′ such that k ν ≪ k ′ ≪ k η will evolve with time. Neglecting the resistive-diffusion term for the time being (η = 0), we immediately find the Green's-function solution of (5): where M 0 (k ′ ) is the initial spectrum. This solution shows that (i) the amplitude of each Fourier mode grows exponentially in time at the rate (3/4)γ, (ii) the number of excited modes [i.e., the width of the lognormal envelope in (10)] also grows exponentially at the rate (4/5)γ, (iii) the peak of the excitation originating from each initially present mode k ′ moves towards larger k: k peak = k ′ exp [(3/5)γt], leaving a power spectrum ∼ k 3/2 behind. The cumulative effect is the total magnetic-energy growth at the rate 2γ. The solution (10) becomes invalid when the magnetic excitation reaches the velocity scale k −1 ν and/or the resistive scale k −1 η . In the former case, including the k dependence ofγ as discussed at the end of section 2.2 would be the minimal adjustment necessary to model large-scale effects. However, the bulk of the magnetic energy is at subviscous scales and is not significantly affected by such modifications. We shall, therefore, continue to rely on (5) and ask what happens when the resistive scale is reached.

The resistive regime
The asymptotic form of the magnetic-energy spectrum in the resistive regime can be determined by solving an eigenvalue problem. Seeking the solution in the form [6] we get The solutions of (12) are Bessel functions. Demanding that Φ(ξ → ∞) = 0, we find where K ν is the Macdonald function.
The eigenvalue λ must be determined from the boundary condition at small k. This can be done in several different ways leading to the same result in the asymptotic case P r ≫ 1 (see, e.g., [39,8] and the exhaustive reference list in [8]). Here we choose some finite IR cutoff k * and impose a zero-flux boundary condition at k = k * . This gives [see equation (6)] Substituting (13) yields the following transcendental equation for λ = λ(ξ * ): It is immediately clear that (15) has no solutions with real ν(λ), i.e., with λ ≥ 3/4. Indeed, such solutions could only exist if the first term were positive: ν(λ) > 5/2, which would imply λ > 2. Since the total energy cannot grow faster than at the rate 2γ, λ > 2 is impossible. Now, the Macdonald function K ν(λ) (ξ) of imaginary order ν(λ) = iν(λ) oscillates as it approaches ξ = 0. Since the spectrum must be a nonnegative function, all zeros of K iν(λ) (ξ) must lie to the left of the IR cutoff ξ * : As we are considering a large-P r case, we must take ξ * = k * /k η ∼ P r −1/2 ≪ 1. The condition (16) then impliesν(λ) ≪ 1. In these limits, (15) reduces to Thus, in the large-P r limit, the zero-flux boundary condition (14) is equivalent to requiring that the spectrum should vanish at the IR cutoff: Φ(ξ * ) = 0 (cf. [8]). Taking the first (rightmost) zero in (17), so that (16) is observed, we find the desired eigenvalue For ξ ≫ ξ * , the spectrum is well described by (13) with ν(λ) = 0. Thus, in the kinematic regime with diffusion, the magnetic-energy spectrum has a constant profile with a k 3/2 scaling in the range k ν ≪ k ≪ k η [see (13)], and all modes grow exponentially at the same rate ≃ (3/4)γ [6]. Most importantly, the bulk of the magnetic energy is concentrated at the resistive scales ∼ k −1 η .

The onset of nonlinearity
When the energy of the magnetic field reaches values such that the back reaction on the flow cannot be neglected, the dynamo becomes nonlinear and the velocity u can no longer be considered as given. It has to be determined self-consistently from the Navier-Stokes equation: where ν is the fluid viscosity, f is some random large-scale forcing, p is the pressure determined from the incompressibility condition ∇ · u = 0, and p and B have been normalized by ρ and (4πρ) 1/2 , respectively, ρ being the (constant) density of the medium. From equation (19), it is evident that the nonlinear regime is characterized by at least partial balancing of the inertial and magnetic-tension terms: u · ∇u ∼ B · ∇B. In order to correctly estimate the tension force B · ∇B, it is not enough to know the energy of the field and its characteristic scale as determined from the isotropic spectra in section 2. Indeed, the gradient in B · ∇B is locally aligned with the field direction, so one must invoke the geometrical structure of the field lines. As was shown in our earlier studies [3,4,26], the small-scale fields generated by the kinematic dynamo are organized in folds: the smallness of their characteristic scale is due to rapid spatial-direction reversals, while the field lines remain mostly straight up the scale of the stretching eddy. It follows that, during the kinematic stage, B · ∇B ∼ k ν B 2 . Since the inertial term is u · ∇u ∼ k ν u 2 ν , the nonlinearity becomes important when the small-scale magnetic energy equalizes with the energy of the viscous-scale eddies: [3]. Note that if the folded structure of the field had been ignored, then estimating B · ∇B ∼ k η B 2 would have led to the incorrect conclusion that nonlinearities become important at the much lower magnetic energies: B 2 ∼ u 2 ν /P r 1/2 .

The nonlinear-growth stage
The magnetic back reaction must then act to suppress the random shear flows associated with these eddies. However, the next-larger-scale eddies still have energies larger than that of the magnetic field, though the turnover rate of these eddies is smaller. We conjecture that these eddies will continue to amplify the field in the same, essentially kinematic, fashion. The scale of these eddies is ℓ s > k −1 ν , the folds are elongated accordingly (see remark in section 4.3), and the tension force is B · ∇B ∼ B 2 /ℓ s . The corresponding inertial term is u · ∇u ∼ u 2 ℓs /ℓ s , so these eddies become suppressed when B 2 ∼ u 2 ℓs , whereupon it will be the turn of the next-larger eddies to provide the dominant stretching action. At all times, it is the smallest unsuppressed eddy that most rapidly stretches the field. Thus, we define the stretching scale ℓ s (t) at a given time t by In equation (2), γ(t) is now the turnover rate of the eddies of scale ℓ s (t): where χ is some constant of order unity. The physical meaning of this equation is as follows. The turbulent energy injected at the forcing scale cascades hydrodynamically down to the scale ℓ s (t) where a finite fraction χ of it is diverted into the small-scale magnetic fields. Equation (21) implies that the magnetic energy should grow linearly with time during this stage: W (t) ∼ ǫt. Furthermore, comparing (21) and (2), we get where α is a constant order-one coefficient to be determined in section 4.3.
We stress that the field is still organized in folds of characteristic parallel length ℓ s (t) with direction reversals perpendicular to B at the resistive scale. Note that the resistive scale now changes with time: a straightforward estimate gives k η (t) ∼ [γ(t)/η] 1/2 ∼ (ηt) −1/2 . Thus, the resistive cutoff moves to larger scales. Indeed, as the stretching slows down with decreasing γ(t), selective decay eliminates the modes at the high-k end of the spectrum for which the resistive time is now shorter than the stretching time and which, consequently, cannot be sustained anymore. In section 4.3, we shall see that the magnetic-energy spectrum evolves in a self-similar fashion during this stage.
Note that some fluid motions do survive at scales below ℓ s (t) and down to the viscous cutoff. These are the motions that do not amplify the magnetic field. They are waves of Alfvénic nature Figure 1. Sketches of the kinetic and magnetic energy spectra in the nonlinear regime. Left: the nonlinear-growth stage (section 3.2). During the kinematic stage, the velocity spectrum E(k) is Kolmogorov, the magnetic spectrum M (k) is peaked at k η ∼ P r 1/2 k ν . In the nonlinear-growth stage, the magnetic field is amplified by the eddies at the stretching scale ℓ s ∼ k −1 s . They are the eddies whose energy is equal to the total magnetic energy at any given time t. The motions with k < k s (t) are still hydrodynamic, but those with k > k s (t) are Alfvénic. As the magnetic energy grows, both k s and k η decrease. Right: the approach to saturation (section 3.3). Selective decay eliminates high-k modes, k η decreases, the stretching rate γ ∼ ηk 2 η decreases accordingly, so the magnetic energy stays constant.
propagating along the folds of direction-reversing magnetic fields (Appendix A). A finite fraction of the hydrodynamic energy arriving from the large scales is channelled into the turbulence of these waves (figure 1). Since the Alfvénic motions are dissipated viscously, not resistively, we shall assume that they do not have a secular effect on the folded small-scale field, i.e., that they do not significantly change its scale or energy (see Appendix A and further discussion in section 3.4).
The process of scale-by-scale suppression of the stretching motions continues until the magnetic energy becomes comparable to the energy W 0 of the outer-scale eddies (i.e., to the energy of the turbulence) and ℓ s ∼ k −1 f , the forcing scale. The time scale for this to happen is easily estimated to be the turnover time of the outer-scale eddies t 0 ∼ (k f u 0 ) −1 ∼γ −1 Re 1/2 , where u 0 is the velocity of these eddies. The resistive cutoff scale is now larger than during the kinematic stage: we denote the resistive cutoff in the kinematic regime. However, this is not a very significant shift: for example, in the interstellar medium, where Re ∼ 10 4 , it is by a meager factor of 0.1. Recall that k ν /k f ∼ Re 3/4 and k (kin) η /k ν ∼ P r 1/2 , so k η (t 0 )/k f ∼ (Re P r) 1/2 , i.e., the scale separation between the field-generating motions and the resulting magnetic fields is actually increased. The magnetic field is still small scale and arranged in folds of length k −1 f with direction reversals at the resistive scale.

The approach to saturation
At this point, the hydrodynamic energy stirred up at the forcing scale is transfered directly into the small-scale folded magnetic field (i.e., in (2), γ ∼ k f u 0 ). As there are no scales in the system larger than k −1 f , there can be no further growth of the magnetic energy. Indeed, the tension in the folds is now B · ∇B ∼ k f W and the force balance in (19) requires B · ∇B ∼ u · ∇u ∼ k f W 0 , so W ∼ W 0 . Note that setting dW/dt = 0 in (2) gives γ(t) = ηk 2 rms (t). This is a reflection of the fact that the resistive cutoff k η ∼ k rms is now determined by the zero-sum balance between the field amplification (stretching) and the resistive dissipation: However, at scales larger than the resistive cutoff, i.e., for k < k η , the field-amplification rate γ is larger than the resistive-dissipation rate ηk 2 ( figure 1, right panel). Since the total magnetic energy is not allowed to grow anymore, the growth of modes with k < k η must be compensated by the decrease in the value of k η , i.e., by further selective decay of the high-k modes. This, in turn, should lead to further decrease of γ in accordance with the resistive balance (23). We shall see in section 4.4 that the evolution of the magnetic-energy spectrum in this regime is likely to be self-similar with γ(t) ∼ 1/t and k η (t) ∼ (ηt) −1/2 . Thus, we expect that this last stage of the nonlinear dynamo should be characterized by γ dropping below the turnover rate k f u 0 of the outer-scale eddies. Therefore, the rate of the energy transfer into the small-scale magnetic field decreases, as this channel of energy dissipation becomes inefficient. Instead, we expect the energy injected by the forcing to be increasingly diverted into the Alfvénic turbulence that is left throughout the inertial range in the wake of the suppression of the stretching motions.

The fully developed isotropic MHD turbulence
We saw in the preceding sections that k η decreases both in the nonlinear-growth stage and during the subsequent slower approach to saturation. It is very important to understand how far the decrease of k η can proceed. We recall that, in our arguments so far, we have disregarded the Alfvénic component of the turbulence. It is not, however, justified to do so after the small-scale magnetic energy reaches the velocity scales. Indeed, the decrease of k η (t) is basically a consequence of the balance between the field-amplification and the resistive-dissipation terms in the energy equation (2): as γ(t) drops, so does k η (t) ∼ k rms (t); once the characteristic scale of the small-scale fields reaches the viscous scale k −1 ν , the Alfvénic turbulence will start to affect k rms in an essential way: since the waves are damped at the viscous scale, k rms must stabilize at k rms ∼ k ν . The resulting turbulent state features folded magnetic fields reversing directions at the viscous scale plus Alfvén waves in the inertial range propagating along the folds. There are then two possibilities: either (i) this represents the final steady state of the isotropic MHD turbulence, or (ii) further evolution will lead to unwinding of the folds and continued energy transfer to larger scales, so that the spectrum will eventually peak at the outer scale and have an Alfvénic power tail extending through the inertial range. The turbulence in the inertial range would in the latter case be of the usual Alfvén-wave kind [27,28,29], where the large-scale magnetic fluctuations provide a mean field along which the inertial-range Alfvén waves propagate. A generic MHD turbulent steady state would thus emerge at the end of the nonlinear-dynamo evolution. This dichotomy cannot be resolved at the level of the present model and requires further study. However, in either case, k rms ∼ k ν , so the statistical steady state of the fully developed isotropic MHD turbulence is characterized by the equalization of the resistive and viscous scales.
The time scale at which such an equalization is brought about is the resistive time associated with the viscous scale of the turbulence: t η (k ν ) ∼ ηk 2 ν −1 . We immediately notice that, in order for the slow approach to saturation described in section 3.3 to be realizable, this time scale must be longer than the turnover time of the outer-scale eddies: t η (k ν ) ≫ t 0 , which requires P r ≫ Re 1/2 . In other words, the resistive scale at the end of the nonlinear-growth stage must still be smaller than the viscous scale: k η (t 0 )/k ν ∼ P r 1/2 Re −1/4 ≫ 1. This constitutes the "true large-P r regime," which is the one relevant for astrophysical plasmas (thus, in the interstellar medium, Re ∼ 10 4 , P r ∼ 10 14 ). In this regime, the magnetic energy saturates at the equipartition level, W ∼ W 0 . From the energy balance (2), we get an estimate of the amount of turbulent power that is still dissipated resistively: χǫ ∼ γW ∼ ηk 2 ν W 0 ∼ ǫ Re 1/2 P r −1 . Thus, χ ≪ 1, i.e., most of the injected power now goes into the Alfvénic motions (which are dissipated viscously).
The other possibility is P r Re 1/2 . In this case, the nonlinear-growth stage described in section 3.2 is curtailed with magnetic energy at a subequipartition value: and resistivity continuing to take a significant (χ is order one) part in the dissipation of the turbulent energy. If no further evolution takes place [possibility (i) above], (24) gives an estimate of the saturation energy of the magnetic component of the turbulence. This regime is characterized by Prandtl numbers which may still be much larger than unity, but are not large enough to capture all of the physics of large-P r dynamo. Most of the extant numerical simulations appear to be in this regime (see section 5).
We should now like to examine the implications of the proposed scenario in a somewhat more quantitative fashion, assuming from now on P r ≫ Re 1/2 ≫ 1. A busy reader with no time for these slow-paced developments may skip to section 5.

Formulation of the model
The expression (22) suggests a way to build up the phenomenological picture just laid out into a semiquantitative model of the evolution of the magnetic spectrum throughout the nonlinear regime. Let us keep the general form of equation (5), but replaceγ by a time-dependent function γ(t), for which we stipulate the following ad hoc expression: where c 1 and c 2 are some numerical constants of order unity, E(k) is the hydrodynamic energy spectrum in the absence of magnetic field, k s (t) is defined by the second expression in (25), and W (t) is the total magnetic energy at time t. The Fokker-Planck model we propose for the nonlinear evolution of the magnetic spectrum is then We now assume the following crude model profile of the hydrodynamic energy spectrum of the Kolmogorov turbulence in the absence of magnetic field: where C K ≃ 1.5 is the Kolmogorov constant. The value of k ν is set to enforce consistency of the model profile (27) with the exact relation ǫ = 2ν ∞ 0 dk k 2 E(k): Using (27), it is straightforward to derive from (25) where are of the order of the energies of the viscousand of the outer-scale eddies, respectively. Here,γ = c 1 (ǫ/2ν) 1/2 follows from evaluating (25) for W → 0, k s = k ν . Note that The kinematic regime γ(t) ≃γ (section 2), the nonlinear-growth regime (22) (section 3.2), and the long-time asymptotic γ(t) → 0 (section 3.3) are all contained in the model expression (29).
In the above formulae, E(k) is the hydrodynamic turbulence spectrum in the absence of the magnetic field, and we have effectively assumed that, at each time t during the evolution of the dynamo, E(k) remains unaffected by the magnetic fields at k < k s (t). As was explained in section 3.2, in the suppressed part of the inertial range k > k s (t), Alfvénic turbulence is excited, which draws a finite fraction of the turbulent power and, indeed, in the long run, supplants the small-scale magnetic fields as the main energy-dissipation channel (see section 3.4). This part of the picture is not included in our Fokker-Planck model (26). We believe that it can be considered independently of the evolution of the small-scale magnetic-energy spectrum as long as the bulk of the magnetic energy stays below the viscous scale (see sections 3.4 and 4.5).

Self-similar solutions
The Fokker-Planck model (26) leads to the emergence of two distinct intermediate regimes where the magnetic-energy spectrum evolves in a self-similar fashion. In both regimes, the bulk of the magnetic energy moves from small to larger scales. We stress that the asymptotic self-similar solutions presented below are not sensitive to the particular form (29) of γ(t). The essential features are those discussed in sections 3.2 and 3.3. The expression (29) will only be used to validate our solutions numerically.
Let us now study equation (26). We notice that, if γ(t) can be represented as where α is a numerical constant, equation (26) admits a self-similar solution where the function Φ(ξ) satisfies Clearly, for γ(t) ∝ 1/t, the appearance of self-similarity is inevitable on dimensional grounds because there is no fixed time scale in the problem in this case. The parameters C and µ in (34) must be determined from the assumptions about the physics of the problem. The value of α is fixed by the solvability of equation (35).

The first self-similar regime: nonlinear growth
In section 3.2, we argued that W (t) should grow linearly with time during the intermediate stage of the nonlinear dynamo when the magnetic back reaction gradually eliminates the stretching component of the inertial-range eddies. Thus, (34) should hold with µ = 1. Then, by virtue of (22), γ(t) has the form (31) required for the self-similar solution to exist. The model expression (29) correctly reproduces this regime: substituting (31) into (29), we get Comparing the above asymptotic with equation (21), we note that 1/2α can be interpreted as the fraction of the turbulent power received by the magnetic fields that actually contributes to the growth of their energy. The rest is dissipated resistively. The upper time limit t 0 for this regime is the turnover time of the outer-scale eddies.
In order to determine the value of α, we have to solve the equation (35) with µ = 1 and enforce the right boundary conditions. This is an eigenvalue problem very similar to that which we considered in section 2.4. The solution is again expressible in terms of the Macdonald function: It is now necessary to implement the boundary condition at small k. Again, we should like to pick some finite IR cutoff k * and to impose the zero-flux condition (14) (see discussion in section 4.5 on the choice of k * ). However, strictly speaking, introducing the new dimensional parameter k * breaks the self-similarity of the problem: instead of (32), we must now write where ξ * = k * /k η (t) = k * (10ηα −1 t) 1/2 is a new time-dependent dimensionless variable. Since the non-self-similar solution (39) has a self-similar asymptotic when ξ * → 0, we can still use (37), provided ξ * ≪ 1. This solution remains valid as long as t ≪ αt η , where t η = (10ηk 2 * ) −1 is the resistive time associated with the IR cutoff k * . In large-P r media, this is usually longer than the nonlinear-growth regime itself lasts according to (36) (see sections 3.4 and 4.5 though).
The zero-flux boundary condition for the self-similar solution (32) has the form (14) with ξ * = k * /k η (t). Substituting (37) into (14), we get the following equation for α: As explained above, this equation must be considered in the limit ξ * ≪ 1 and the resulting α must remain finite as ξ * → 0. Again there are no solutions with real ν(α). For imaginary ν(α) = iν(α), the requirement that the spectrum should have no nodes to the right of ξ * ≪ 1 again impliesν(α) ≪ 1, and we get  (26). This is the same run as in figure 2. The spectra are normalized by W (t).
The solution is (see figure 2) Since t ≪ t η , the spectral profile is well described by (37) with ν(α) = 0 and σ(α) ≃ 1.9. Thus, in the first self-similar stage, the scale-by-scale suppression of the inertial-range eddies leads to gradual renormalization of the resistive cutoff k η (t), which is moving leftwards as (ηt) −1/2 in agreement with the estimate in section 3.2 (see figures 2, 3). The turbulent power arriving into the smallest unsuppressed eddies is (partially) transfered into the small-scale magnetic fields. The fraction 1/2α (∼ 16%) of it contributes to the linear-in-time growth of the total magnetic energy, while the rest is dissipated resistively. Since this stage of the dynamo lasts up to t 0 ∼γ −1 Re 1/2 , we recover the physical estimates made in section 3.2: Remark on the elongation of the folds. In section 3.2, we assumed the growth of the effective stretching scale ℓ s (t) during the nonlinear-growth stage to be accompanied by simultaneous elongation of the folds, so that the parallel scale of the magnetic field (characteristic length of the folds) would always be ℓ (t) ∼ ℓ s (t). This assumption can now be supported in the following way. Consider a straight segment of a field line of length ℓ . In the absence of resistivity, volume-and flux-conservation constraints in an incompressible flow imply ℓ /B ∼ const, so d log ℓ /dt = (1/2) d log W/dt = γ(t).
With resistivity included, the growth of the magnetic field is offset by resistive diffusion. As we showed above, this lowers the net amplification rate of the magnetic energy to α −1 γ(t) in the first self-similar regime. However, the length of the field line continues to grow at the rate γ(t), so d log ℓ /dt = α d log W/dt. Using the nonlinear balance (20) and Kolmogorov scaling for velocity gives d log ℓ /d log ℓ s = 2α/3 ≃ 2. By itself, this would indicate that ℓ grows even faster than ℓ s . However, a straight section of a field line remains straight during stretching only as long as it is still short compared to the stretching-eddy size. Once ℓ grows to be comparable to ℓ s , it is limited to ℓ ∼ ℓ s by the curvature in the stretching eddies. The conclusion from this argument is that the length of the folds can quickly adjust to the changing stretching scale. Note that this conclusion does not depend on the precise value of α: it is sufficient that α −1 ≤ 2/3.

4.4.
The second self-similar regime: approach to saturation Thus, after time t 0 ∼γ −1 Re 1/2 , the energy of the small-scale magnetic fields is very close to the equipartition energy W 0 . We observe that in this regime it is still possible to have a self-similar solution in the form (32): the energy is now (nearly) constant, so µ = 0 in (34). Accordingly, substituting (31) into (29), we find We now set µ = 0 in (35) and proceed to solve the resulting eigenvalue problem in exactly the same way as we did in section 4.3. All the same formulae and considerations apply, except that the expressions for α ± in (38) and forν(α) in (41) are now as follows: The new eigenvalue is For t ≪ t η , the spectrum is given by (37) with ν(α) = 0 and σ(α) ≃ 2.2. Thus, in the second self-similar stage, the magnetic energy is nearly constant, but the resistive cutoff k η (t) continues to decrease as (ηt) −1/2 , so the magnetic energy is redistributed from small to larger scales (see figures 2, 3). This is a very slow process, but it can, in principle, go on for a very long time: until self-similarity is broken at t ∼ t η = (10ηk 2 * ) −1 . Note that, according to (43), the magnetic energy changes sufficiently slowly to make the left-hand side of the energy equation (2) subdominant in comparison with the field-amplification and resistive-dissipation terms in the right-hand side. This is consistent with the physical argument of section 3.3.

Steady states at the equipartition energy and below
The self-similar solutions obtained above break down after a finite time ∼ t η due to the presence of the finite IR cutoff k * . Once the resistive cutoff k η becomes comparable to k * , the system arrives at the final steady state. The validity of this state is questionable since our model was derived under the assumption that the magnetic fields were at smaller scales than the velocity field, i.e. k η ≫ k ν . Furthermore, the specifications of the steady-state solution that the model yields clearly will have an order-one dependence on the value of the IR cutoff k * , on the type of the boundary condition imposed, and on the form of the expression (29). However, we proceed to discuss this solution because it highlights several important issues that require further study.
In the steady state, γ(t) → γ(∞) = const, so the corresponding solution of equation (26) is simply the λ = 0 eigenmode of the kinematic equation (5) withγ replaced by γ(∞): using equations (11) and (13), we get The value of γ(∞) must be determined from the boundary condition at k * : This is an equation for ξ * , whose value in turn determines γ(∞) = 10ηk 2 * /ξ 2 * . We have to pick the largest root of (47) to ensure that the spectrum (46) has no nodes. This gives ξ * ≃ 0.59, which is order one, as expected. Substituting γ(∞) into (29), we get the magnetic energy W (∞) in the steady state: which, using W 0 /W ν ∼ Re 1/2 ≫ 1 and γ(∞)/γ = 10ηk 2 * /γξ 2 * ∼ k * /k (kin) η 2 ∼ (k * /k ν ) 2 P r −1 , can be rewritten as follows What is the correct choice of k * ? Of course, as long as the characteristic scale of the field is much smaller than that of all fluid motions, stretching or Alfvénic, i.e., k η ≫ k ν , the precise value of k * does not affect the self-similar asymptotics of sections 4.3 and 4.4 in the leading order. What k * represents is the presence in the problem of a spatial scale which leads to the breakdown of self-similarity in the long run (after time t η ). As we explained in section 3.4, once the resistive and the viscous scales become comparable, the evolution of the small-scale fields cannot be separated from that of the inertial-range Alfvénic turbulence. Clearly, beyond this point, the evolution of the spectrum is not described by the self-similar solutions obtained in sections 4.3 and 4.4. If we set k * ∼ k ν , the end of the selfsimilar regime and the stabilization of k η at a value ∼ k ν are reflected by the steady-state solution of our Fokker-Planck model. Note that, with this choice of k * , equation (49) correctly reproduces both W (∞) ∼ W 0 for P r ≫ Re 1/2 and the physical estimate (24) of the subequipartition saturation level for P r ≪ Re 1/2 . We stress that, since the validity of the Fokker-Planck model itself breaks down together with the self-similarity, the specific functional form of its steady-state solution cannot be expected to have predictive value. The ensuing regime is controlled by the interaction between the folded fields and the Alfvén waves and is left for future study.

Discussion
We have presented a physical scenario of the evolution of the nonlinear large-P r dynamo according to which a magnetic-energy spectrum concentrated at velocity scales does emerge, but the time scale for it is resistive, not dynamical: the resistive scale approaches the viscous scale after t η ∼ ηk 2 ν −1 .
In the interstellar medium, the Spitzer [42] value of the magnetic diffusivity is η ∼ 10 7 cm 2 /s, k ν ∼ 10 −16 cm −1 [8], so t η ∼ 10 17 years, which exceeds the age of the Universe by 7 orders of magnitude! The conclusion is that this mechanism cannot be invoked as a workable feature of the galactic dynamo -at least not if the dissipation of the magnetic field is controlled by the Spitzer magnetic diffusivity.
The statistically steady solutions discussed in section 3.4 represent the possible long-time asymptotic states of the fully developed isotropic large-P r MHD turbulence. Since the approach to these states is controlled by the resistive time scale associated with the velocity spatial scales, they have practically no relevance for astrophysical plasmas, which have large Re and huge P r. Studying these states may, therefore, appear to be a purely academic exercise in fundamental turbulence theory. While the importance of such pursuits ought never to be underestimated, we should like to point out another, more practical angle from which our results could be viewed. The enormous scale separations so often encountered in astrophysical objects are not accessible in the turbulence simulations that constitute the state of the art on Earth, so one usually has to be satisfied with only very modest Re and P r. Most of the numerical studies of the small-scale dynamo have so far adopted the strategy of resolving a reasonable hydrodynamic inertial range while only allowing for Prandtl numbers of at most order 10 [20,21,24,25,26,43]. In this context, our steady-state solution and formula (48) point to an important possibility. When the scale separations in the problem are not very large, the system may converge to a steady state with a subequipartition magnetic energy W (∞) < W 0 . Subequipartion saturated states have indeed been reported in the numerical experiments cited above.
Since the condition for the true asymptotic regime is P r ≫ Re 1/2 (section 3.4), numerical experiments with very high resolution are required to adequately study the large-P r MHD turbulence.
In conclusion, we list some of the unresolved issues that require further study.
• The detailed mechanism of the flow suppression by the small-scale magnetic fields is a longstanding physical problem [44,45,46,47]. In our model, we have conjectured the suppression of the shearing flows but allowed for oscillatory Alfvénic motions, that do not stretch the field (section 3.2). Strictly speaking, only the parallel component of the shear tensor ∇u must be suppressed in order to block the growth of the magnetic energy. The perpendicular components, if unaffected, could perhaps mix the field (without amplifying it) via a quasi-two-dimensional field-line-interchange mechanism and thus prevent any increase in the field's characteristic scale. In fact, this seems to be the physics behind the recent numerical results on the large-P r MHD turbulence with a fixed uniform mean field [48]. However, it is unclear to what extent such motions can remain quasi-two-dimensional and avoid twisting up the magnetic field, which would then resist further mixing -particularly, in our case of a folded field-line structure without an externally imposed mean field.
• Our physical picture of the final approach to saturation (section 3.3) is not based on a solid phenomenological argument such as that for the nonlinear-growth stage, and remains incompletely understood. Indeed, it remains to be proven definitively that the selective decay will continue after the small-scale magnetic energy reaches the energy of the outer-scale eddies. Due to resolution constraints discussed above, no numerical evidence of this regime is available as yet.
• Two possibilities for the long-time behaviour of the isotropic MHD turbulence were identified in section 3.4: saturation in the generic Alfvénic MHD turbulent state and saturation in a state where a significant amount of the magnetic energy remains tied up in the viscous-scale fields. Which one is realized depends on the way the small-scale folded fields interact with the inertial-range Alfvén-wave turbulence. We should like to observe that there is no numerical evidence available at present that would confirm that, the isotropic forced MHD turbulence without externally imposed mean field -at any Prandtl number -attains the Alfvénic state of scale-by-scale equipartition envisioned in Iroshnikov-Kraichnan [27,28] or Goldreich-Sridhar [29] phenomenologies. In fact, medium-resolution numerical investigations [26] rather seem to support the final states with smallscale energy concentration even for P r = 1. This does not mean that the Alfvénic picture is incorrect per se. However, all existing phenomenologies of the Alfvénic turbulence depend on Kraichnan's [28] assumption that the most energetic magnetic fields are those at the outer scale. This is automatically satisfied if a finite mean field is imposed externally [49,50,48]. However, it remains to be seen if such a distribution of energy is set up self-consistently when the turbulence is isotropic. The alternative identified here involves Alfvénic motions superimposed on small-scale folded fields.
• An extension of the MHD description itself that may change the properties of the small-scale magnetic turbulence is the incorporation of the Braginskii [51] tensor viscosity [52,26,53]. Even at magnetic energies small enough for the kinematic approximation to hold, the plasma is already well magnetized and the anisotropy of the viscous stress tensor leads to suppression of the velocity diffusion perpendicular to the local magnetic field [53]. This anisotropy is all the more important in view of the locally anisotropic (folding) structure of the magnetic field itself [3,53].
• Another important extension is to allow compressible motions. The interstellar turbulence is sonic at the outer (supernova) scale but becomes subsonic and predominantly vortical in the inertial range (cf. [54,55]). This is the commonly accepted justification for the use of incompressible MHD in the models of galactic dynamo. Of course, a certain amount of compressive motion is always present in the interstellar medium (as well as in other astrophysical plasmas), and it has been suggested in the literature that the compressibility of the motions and interactions between Alfvén waves and density inhomogeneities can be important in the description of turbulence (see, e.g., [56] Figure A1. Alfvén waves propagate along folded magnetic fields. and many recent references cited therein). The treatment of the kinematic-dynamo problem for a model compressible case can be found in [8,3]. However, it remains to be understood whether and how weak compressibility affects the basic properties of the nonlinear dynamo.
It is a very old observation that any extension of the domain of one's knowledge lengthens the perimeter along which this domain borders on the unknown. The perceptive reader must have realized that this work raises more questions than it gives definitive answers. We shall address these questions in our future investigations.
In (A.3), we must have ℓ −1 s ≪ k ≪ k ν , so the frequency ω is larger than the characteristic time scale of the growth of the small-scale magnetic energy: ω = k v A ≫ v A /ℓ s ∼ W 1/2 /ℓ s ∼ u ℓs /ℓ s ∼ γ. Since the resistive-dissipation rate of the small-scale fields is ηk 2 rms ∼ γ by virtue of (2), the Alfvén waves are mostly dissipated viscously via the MHD turbulent cascade. The velocity spectrum of this Alfvénic turbulence is subject to considerations similar to those arising in the Iroshnikov-Kraichnan [27,28] or the Goldreich-Sridhar [29] phenomenologies. Further investigations of this and related issues will be published elsewhere.