text:
# Constraining URCA cooling of neutron stars from the neutron radius of 208Pb

## I Introduction

Neutron stars are created with very high temperatures in supernova explosions. Indeed, neutrinos observed from SN1987A indicate a neutrinosphere temperature that could be as high as 5 MeV raf . Neutrons stars then cool, primarily by neutrino emission Pe92 . In the standard scenario, the modified URCA reaction,

$$n+nn+p+e^{}+\overline{\nu }_e,$$
(1)
emits neutrinos from the volume of the star. This process, however, is relatively slow as a second nucleon is necessary to conserve momentum.

Recent X-ray observations of the neutron star in 3C58 3c58 , Vela vela , and Geminga geminga  indicate low surface temperatures. Moreover, the low quiescent luminosity in the transiently accreting binaries KS 1731-260 wijnands  and Cen X-4 page  suggest rapid cooling. As X-ray observatories progress and our knowledge of neutron-star atmospheres and ages improve, additional “cold” neutron stars may be discovered. Such low surface temperatures appear to require enhanced cooling from reactions that proceed faster than the modified URCA process of Eq. (1). Enhanced cooling may occur via the weak decay of additional hadrons such as pion or kaon condensates kaon , hyperons Pr92 , or quark matter quarks . Yet perhaps the most conservative enhanced-cooling mechanism is the direct URCA process urca ; leinson  of neutron beta decay followed by electron capture:

$`np+e^{}+\overline{\nu }_e,`$ (2a)
$`e^{}+pn+\nu _e.`$ (2b)
This mechanism is not “exotic” as it only requires protons, neutrons, and electrons—constituents known to be present in dense matter. However, to conserve momentum in Eq. (2a) the sum of the Fermi momenta of the protons plus that of the electrons must be greater than (or equal to) the neutron Fermi momentum. This requires a relatively large proton fraction.

Yakovlev and collaborators yakovlev  are able to reproduce measured neutron-star temperatures using a relativistic mean-field equation of state that allows direct URCA for neutron stars with masses above $`1.358M_{}`$ ($`M_{}=`$ solar mass). In contrast, Tsuruta and collaborators t5  rely on pion condensation to reproduce the measured temperatures. They argue that microscopic calculations of neutron-rich matter fp  using nonrelativistic nucleon-nucleon interactions yield too small a proton fraction for the URCA process to operate. Unfortunately, these microscopic calculations depend on a poorly known three-nucleon force and on relativistic effects that could end up increasing the proton fraction at high densities.

Superconductivity and superfluidity can greatly influence neutron-star cooling Ya99 ; Pa00 . For temperatures much lower than the pairing gap, pairing correlations suppress exponentially the rate of many cooling reactions. Yet for temperatures of the order of the pairing gap, the thermal breaking and subsequent reformation of nucleon “Cooper” pairs promotes an additional neutrino-emission mechanism that rapidly cools the star pairing . However, it has been argued in Ref. yakovlev  that this mechanism alone is unlikely to explain the low temperature of some neutron stars. This is because for a large enough neutron-pairing gap, pair breaking would rapidly cool all neutron stars at a rate almost independent of the mass of the star. This would disagree with observations of some warm neutron stars. Tsuruta and collaborators have claimed that microscopic calculations with a high proton concentration show a small proton pairing gap tpairing . If so, a direct URCA process (one not controlled by pairing correlations) will cool a star so quickly that thermal radiation would become invisible t5 . However, we caution that drawing definitive conclusions from microscopic calculations of pairing gaps may be premature, as significant uncertainties remain in the interactions, equation of state, composition, and phases of high-density matter.

Although the precise mechanism remains unknown, some kind of enhanced cooling appears to be required to explain the recent observations of cold neutron stars. While the need for exotic phases of matter is appealing, more conventional cooling scenarios, such as the direct URCA process, can not be dismissed on purely theoretical grounds. Moreover, neutron-star observations alone may not be able to resolve the detailed mechanism of enhanced cooling. Thus, we consider complementary laboratory experiments that could help us confirm (or possibly dismiss) the direct URCA process. This can be achieved by constraining the symmetry energy of dense matter. The symmetry energy describes how the energy of (asymmetric) nuclear matter increases as one departs from equal numbers of neutrons and protons. The proton fraction $`Y_p=Z/A`$ of nuclear matter in beta-equilibrium is sensitive to the symmetry energy urca . A large symmetry energy imposes a stiff penalty on the system for upsetting the $`N=Z`$ balance hereby forcing it to retain a large proton fraction.

Energetic heavy-ion collisions probe the symmetry energy at high nuclear densities li . Possible observables include the ratio of $`\pi ^{}`$-to-$`\pi ^+`$ production and the neutron-proton differential collective flow. However, these reactions may suffer from important uncertainties associated with the complex strong interactions of the heavy-ion collisions. Thus, we rely on a purely electroweak reaction that can be unambiguously interpreted. Parity violating elastic electron scattering from a heavy nucleus is sensitive to the neutron density. This is because the weak charge of a neutron is much larger than the weak charge of a proton. The Parity Radius Experiment at the Jefferson Laboratory aims to measure the neutron radius in 208Pb to a 1% accuracy ($`\pm 0.05`$ fm) prex . This measurement can be both accurate and model independent bigpaper . In Ref. brown  Brown showed that the neutron radius of 208Pb determines the pressure of neutron-rich matter at normal densities which, in turn, is related to the density dependence of the symmetry energy prakash . In an earlier work we showed how the neutron radius of 208Pb determines properties of the neutron-star surface, such as the transition density from a solid crust to a liquid interior prl . Furthermore, we argued that by comparing the neutron radius of 208Pb (a low-density observable) to the radius of a neutron star (a high and low-density observable) evidence may be provided in support of a phase transition in dense matter radii .

In the present work we show how the neutron radius of a heavy nucleus (such as 208Pb) controls the density dependence of the symmetry energy. Unfortunately, the density dependence of the symmetry energy ($`da_{\mathrm{sym}}/d\rho `$) is poorly known. Thus a measurement of the neutron radius of 208Pb seems vital, as it will constrain the density dependence of the symmetry energy at low density. This, in turn, will allow a more reliable extrapolation of the symmetry energy, and thus a more reliable determination of the proton fraction at the higher densities required in the study of neutron-star structure. While in principle collective modes of nuclei, such as the giant-dipole or isovector-monopole resonances, are sensitive to $`da_{\mathrm{sym}}/d\rho `$, in practice this sensitivity is small. Moreover, the parameter sets used in the calculations (see various tables) have been adjusted so that well known ground-state properties remain fixed while changing the neutron radius. This shows that existing ground-state information, such as charge densities or binding energies, do not determine the neutron radius uniquely. Thus the need for a new measurement—such as the neutron radius in 208Pb—that will provide important information on $`da_{\mathrm{sym}}/d\rho `$.

The paper has been organized as follows. In Sec. II, relativistic effective-field theories for both dense matter and finite nuclei are discussed. A large number of parameter sets are considered so that the density-dependence of the symmetry energy may be changed while reproducing existing ground-state data. In Sec. III, results for the equilibrium proton fraction as a function of baryon density are presented using interactions that predict different neutron radii in 208Pb. Our summary and conclusions are offered in Sec. IV. In particular, we conclude that for models with a large neutron skin in 208Pb ($`R_nR_p0.25`$ fm) the symmetry energy rises rapidly with density and the direct URCA cooling of a $`1.4M_{}`$ neutron star is likely. Conversely, if the neutron radius is small ($`R_nR_p0.20`$ fm) it is unlikely that the direct URCA process occurs. In this case, the enhanced cooling of neutron stars may indeed require the presence of exotic states of matter, such as meson condensates, hyperonic, and/or quark matter.

## II Formalism

Our starting point will be the relativistic effective-field theory of Ref. horst  supplemented with additional couplings between the isoscalar and the isovector mesons. This allows us to correlate nuclear observables sensitive to the density dependence of the symmetry energy, such as the neutron radius of 208Pb, with neutron-star properties, such as the threshold mass for URCA cooling. As the density dependence of the symmetry energy is poorly known, uncertainties in these correlations will be explored by considering a wide range of model parameters. The interacting Lagrangian density is thus given by prl ; horst

$`_{\mathrm{int}}`$ $`=`$ $`\overline{\psi }\left[g_\mathrm{s}\varphi \left(g_\mathrm{v}V_\mu +{\displaystyle \frac{g_\rho }{2}}𝝉𝐛_\mu +{\displaystyle \frac{e}{2}}(1+\tau _3)A_\mu \right)\gamma ^\mu \right]\psi `$ (3)
$``$ $`{\displaystyle \frac{\kappa }{3!}}(g_\mathrm{s}\varphi )^3{\displaystyle \frac{\lambda }{4!}}(g_\mathrm{s}\varphi )^4+{\displaystyle \frac{\zeta }{4!}}g_\mathrm{v}^4(V_\mu V^\mu )^2`$
$`+`$ $`g_\rho ^2𝐛_\mu 𝐛^\mu \left[\mathrm{\Lambda }_\mathrm{s}g_\mathrm{s}^2\varphi ^2+\mathrm{\Lambda }_\mathrm{v}g_\mathrm{v}^2V_\mu V^\mu \right].`$
The model contains an isodoublet nucleon field ($`\psi `$) interacting via the exchange of two isoscalar mesons, the scalar sigma ($`\varphi `$) and the vector omega ($`V^\mu `$), one isovector meson, the rho ($`𝐛^\mu `$), and the photon ($`A^\mu `$). In addition to meson-nucleon interactions the Lagrangian density includes scalar and vector self-interactions. The scalar-meson self-interactions ($`\kappa `$ and $`\lambda `$) soften the equation of state (EOS) of symmetric nuclear matter at (and near) saturation density while the $`\omega `$-meson self-interactions ($`\zeta `$) soften the high-density component of the EOS. Finally, the nonlinear couplings ($`\mathrm{\Lambda }_\mathrm{s}`$ and $`\mathrm{\Lambda }_\mathrm{v}`$) are included to modify the density dependence of the symmetry energy prl ; radii .

The energy of neutron-rich matter may be written in terms of the energy of symmetric nuclear matter ($`\rho _p=\rho _n`$) and the symmetry energy $`a_{\mathrm{sym}}(\rho )`$. That is,

$$\frac{E}{A}(\rho ,t)=\frac{E}{A}(\rho ,t=0)+t^2a_{\mathrm{sym}}(\rho )+𝒪(t^4),$$
(4)
where the neutron excess has been defined as

$$t\frac{\rho _n\rho _p}{\rho _n+\rho _p}.$$
(5)
Here $`\rho _n`$ is the neutron and $`\rho _p`$ the proton density, and

$$\rho =\rho _p+\rho _n\frac{2k_\mathrm{F}^3}{3\pi ^2}.$$
(6)
The symmetry energy describes how the energy of the system increases as one moves away from $`\rho _p=\rho _n`$. It is discussed in Ref. radii  where it is shown that it is given by

$$a_{\mathrm{sym}}(\rho )=\frac{k_F^2}{6E_F^{}}+\frac{g_\rho ^2}{12\pi ^2}\frac{k_F^3}{m_\rho ^2},$$
(7)
where $`k_\mathrm{F}`$ is the Fermi momentum, $`E_F^{}=\sqrt{k_\mathrm{F}^2+M^2}`$, and $`M^{}=Mg_s\varphi _0`$ is the effective nucleon mass. Further, the effective rho-meson mass has been defined as follows:

$$m_\rho ^2=m_\rho ^2+2g_\rho ^2\left(\mathrm{\Lambda }_\mathrm{s}g_\mathrm{s}^2\varphi _0^2+\mathrm{\Lambda }_\mathrm{v}g_\mathrm{v}^2V_0^2\right).$$
(8)
The symmetry energy is given as a sum of two contributions. The first term in Eq. (7) represents the increase in the kinetic energy of the system due to the displacement of the Fermi levels of the two species (neutrons and protons). This contribution has been fixed by the properties of symmetric nuclear matter as it only depends on the nucleon effective mass $`M^{}`$. By itself, it leads to an unrealistically low value for the symmetry energy; for example, at saturation density this contribution yields $`15`$ MeV, rather than the most realistic value of $`37`$ MeV. The second contribution is due to the coupling of the rho meson to an isovector-vector current that no longer vanishes in the $`NZ`$ system. It is by adjusting the strength of the $`NN\rho `$ coupling constant that one can now fit the empirical value of the symmetry energy at saturation density. However, the symmetry energy at saturation is not well constrained experimentally. Yet an average of the symmetry energy at saturation density and the surface symmetry energy is constrained by the binding energy of nuclei. Thus, the following prescription is adopted: the value of the $`NN\rho `$ coupling constant is adjusted so that all parameter sets have a symmetry energy of $`a_{\mathrm{sym}}=25.67`$ MeV at $`k_F=1.15`$ fm-1($`\rho =0.10`$ fm-3prl . That is,

$$g_\rho ^2=\frac{m_\rho ^2\mathrm{\Delta }a_{\mathrm{sym}}}{{\displaystyle \frac{k_F^3}{12\pi ^2}}2\left(\mathrm{\Lambda }_\mathrm{s}g_\mathrm{s}^2\varphi _0^2+\mathrm{\Lambda }_\mathrm{v}g_\mathrm{v}^2V_0^2\right)\mathrm{\Delta }a_{\mathrm{sym}}},$$
(9)
where $`\mathrm{\Delta }a_{\mathrm{sym}}(a_{\mathrm{sym}}k_F^2/6E_F^{})`$. This prescription insures accurate binding energies for heavy nuclei, such as 208Pb. Following this prescription the symmetry energy at saturation density is predicted (for $`\mathrm{\Lambda }_\mathrm{s}=\mathrm{\Lambda }_\mathrm{v}=0`$) to be $`37.3`$, $`36.6`$, and $`36.3`$ MeV for the three families of parameter sets considered in this work, namely NL3 NL3 , S271 prl , and Z271prl , respectively (see various tables). Moreover, all these parameter sets reproduce the following properties of symmetric nuclear matter: i) nuclear saturation at a Fermi momentum of $`k_\mathrm{F}=1.30`$ fm-1 with ii) a binding energy per nucleon of $`16.24`$ MeV, and iii) a nuclear incompressibility of $`K=271`$ MeV. These values follow from the successful parametrization of Ref. NL3  and, thus, have been adopted for the other sets (S271 and Z271) as well. Yet the various parameter sets differ in i) their values for the effective nucleon mass at saturation density, ii) the value of the $`\omega `$-meson quartic coupling ($`\zeta 0`$ for set Z271 but vanishes for the NL3 and S271 sets) and iii) the nonlinear couplings $`\mathrm{\Lambda }_\mathrm{s}`$ and $`\mathrm{\Lambda }_\mathrm{v}`$ (see various tables). Note that changing $`\mathrm{\Lambda }_\mathrm{s}`$ or $`\mathrm{\Lambda }_\mathrm{v}`$ modifies the density dependence of the symmetry energy through a change in the effective rho-meson mass \[see Eq. (8)\]. In general, increasing either $`\mathrm{\Lambda }_\mathrm{s}`$ or $`\mathrm{\Lambda }_\mathrm{v}`$ causes the symmetry energy to soften, that is, to grow slower with increasing density. This, in turn, allows for a larger neutron-proton mismatch or equivalently, for a lower equilibrium proton fraction at high density.

The neutron radius of 208Pb also depends on the density dependence of the symmetry energy. A stiff density dependence (i.e., pressure) for neutron matter pushes neutrons out against surface tension, leading to a large neutron radius. The pressure of neutron-rich matter depends on the derivative of the energy of symmetric matter with respect to the density ($`dE(\rho ,t=0)/d\rho `$) and on the derivative of the symmetry energy ($`da_{\mathrm{sym}}/d\rho `$). While the former is well known, at least in the vicinity of the saturation density, the latter is not. Hence, by changing the values of $`\mathrm{\Lambda }_\mathrm{s}`$ or $`\mathrm{\Lambda }_\mathrm{v}`$ one can adjust the density dependence of the symmetry energy $`da_{\mathrm{sym}}/d\rho `$, while keeping a variety of well-known ground-state properties (such as the proton radius and the binding energy of 208Pb) unchanged. Note that parameter sets with a large “pressure”, $`da_{\mathrm{sym}}/d\rho `$, yield a large neutron radius in 208Pb.

## III Results

In this section results are presented for various observables that have been computed using an equation of state for matter composed of neutrons, protons, electrons, and muons in beta equilibrium:

$`np+e^{}+\overline{\nu }_e,`$ (10a)
$`e^{}\mu ^{}+\nu _e+\overline{\nu }_\mu .`$ (10b)
These reactions demand that the chemical potential ($`\mu `$) of the various constituents be related by the following equation:

$$\mu _n\mu _p=\mu _e=\mu _\mu ,$$
(11)
where $`\mu _n`$, $`\mu _p`$, $`\mu _e`$, and $`\mu _\mu `$ represent the chemical potentials of neutrons, protons, electrons, and muons, respectively. Neglecting the rest mass of the electron, Eq. (10b) is equivalent (for $`k_\mathrm{F}^em_\mu `$) to the following equation expressed in terms of the Fermi momenta of the electron and the muon:

$$k_\mathrm{F}^e=\sqrt{\left(k_\mathrm{F}^\mu \right)^2+m_\mu ^2}.$$
(12)
Finally, charge neutrality imposes the following constraint on the system: $`\rho _p=\rho _e+\rho _\mu `$, or equivalently,

$$\left(k_\mathrm{F}^p\right)^3=\left(k_\mathrm{F}^e\right)^3+\left(k_\mathrm{F}^\mu \right)^3.$$
(13)
For a given proton Fermi momentum, the corresponding Fermi momenta for the electrons and the muons are readily obtained by solving Eqs. (12) and (13). With these in hand, Eq. (11) determines the equilibrium neutron ($`Y_n=N/A`$) and proton ($`Y_p=Z/A`$) fractions in the system.

In Fig. 1 the proton fraction $`Y_p`$ for matter in beta equilibrium is shown as a function of the baryon density for all models discussed in the text (see Table 1). The various curves displayed in each panel are for values of $`\mathrm{\Lambda }_\mathrm{s}`$ or $`\mathrm{\Lambda }_\mathrm{v}`$ that give the indicated values for the neutron skin of 208Pb. Note that the neutron skin of a nucleus is defined as the difference between the neutron ($`R_n`$) and the proton ($`R_p`$) root-mean-square radii. All of the curves yield the same proton fraction at low density ($`k_\mathrm{F}=1.15`$ fm-1 or $`\rho 0.1`$ fm-3) because the symmetry energy has been adjusted to $`a_{\mathrm{sym}}=25.67`$ MeV in order to reproduce the binding energy 208Pb. The proton fraction increases more rapidly with density for those parameter sets that yield larger neutron radii in 208Pb (namely, those with a stiffer symmetry energy). Thus, the neutron radius of 208Pb constrains the slope $`dY_p/d\rho `$ at normal densities. This enables one to make a more reliable extrapolation of $`Y_p`$ to the higher densities where it displays some model dependence.

The direct URCA process is viable only when the proton fraction is large enough to conserve momentum in the neutron beta decay reaction of Eq. (2a). Hence, for this reaction to proceed, the Fermi momenta of neutrons, protons, and electrons must satisfy the following relation:

$$k_\mathrm{F}^nk_\mathrm{F}^p+k_\mathrm{F}^e.$$
(14)
The URCA threshold density $`\rho _{\mathrm{URCA}}`$ is defined as the density at which the equality ($`k_\mathrm{F}^n=k_\mathrm{F}^p+k_\mathrm{F}^e`$) is satisfied. Note that in the simplified case of matter without muons, that is, $`k_\mathrm{F}^e<m_\mu `$ and $`k_\mathrm{F}^p=k_\mathrm{F}^e`$, the proton fraction at the onset of the direct URCA process is $`Y_p=1/90.111`$. In the opposite limit, $`k_\mathrm{F}^em_\mu `$, the threshold proton fraction is $`Y_p0.148`$. Thus, the threshold proton fraction must be contained within these two values for all baryon densities (see Tables 2-5). In Fig. 2 the URCA threshold density is displayed as a function of the neutron skin $`R_nR_p`$ of 208Pb. There is a clear tendency for $`\rho _{\mathrm{URCA}}`$ to decrease with increasing neutron skin. Recall that a large neutron skin implies a stiff symmetry energy and a large proton fraction. Thus the onset for the direct URCA process for models with large neutron skins occurs at low baryon densities. Indeed, parameter sets with neutron skins of $`R_nR_p0.24`$ fm yield a relative low URCA density of $`\rho _{\mathrm{URCA}}0.30`$ fm-3. This density is only a factor of two larger than normal nuclear matter saturation density ($`\rho _00.15`$ fm-3). In contrast, if $`R_nR_p0.21`$ fm, the onset for the URCA process is above $`3\rho _0`$.

The structure of spherical neutron stars in hydrostatic equilibrium is solely determined by the equation of state of neutron-rich matter in beta equilibrium. Having specified the equation of state, we determine the mass of neutron stars that may cool via the direct URCA process by integrating the Tolman-Oppenheimer-Volkoff equations. Our treatment of the low-density crust, where the matter in the star is nonuniform, will be discussed in a later work Ho02 . This region, however, has almost no effect on neutron-star masses. As mentioned earlier, we consider matter composed of neutrons, protons, muons, and (massless) electrons. In Fig. 3 we display, as a function of the neutron skin in 208Pb, the mass of a neutron star whose central density equals the URCA density ($`\rho _{\mathrm{URCA}}`$) of Fig. 2. Neutron stars with larger masses, and thus higher central densities, will cool by the direct URCA process; those with lower masses will not. There is an obvious trend for this threshold mass to decrease with increasing $`R_nR_p`$. Recall that the onset for URCA cooling in models with large neutron skins occurs at low baryon densities, thus lower “URCA masses”. If the neutron skin of 208Pb is less than about $`R_nR_p0.20`$ fm, then all parameter sets considered in this work predict that a neutron star of 1.4 $`M_{}`$ will not undergo URCA cooling. Conversely, if $`R_nR_p0.25`$ fm, then all parameter sets allow URCA cooling for 1.4 $`M_{}`$ neutron star. Note that all well measured neutron stars have masses near 1.4 $`M_{}`$ Th99 .

The threshold neutron star mass for the direct URCA process, $`M_{\mathrm{URCA}}`$, depends on both the “critical” URCA density $`\rho _{\mathrm{URCA}}`$ (of Fig. 2) and on the equation of state at low and high densities. Yet the neutron skin of 208Pb constrains only the low density EOS brown , thereby generating the model dependence displayed by Fig. 3. In particular, the radius of a neutron star, although correlated to $`R_nR_p`$, is not determined uniquely by it radii . This suggests that models with a stiff EOS at high density, such as the NL3 parameter set, will yield neutron stars with relatively large radii and low central densities. This implies, for example, that for a fixed central density of $`\rho _{\mathrm{URCA}}=0.3`$ fm-3, the NL3 set (with large radii) generates an URCA mass of $`M_{\mathrm{URCA}}1.4M_{}`$; in contrast, the softer Z271s set (with small radii) yields an URCA mass of only $`M_{\mathrm{URCA}}1M_{}`$. Equivalently, to make a $`M_{\mathrm{URCA}}=1.4M_{}`$ neutron star the NL3 set requires an interior density of $`\rho _{\mathrm{URCA}}=0.3`$ fm-3, while a central density almost twice as large is needed for the Z271s set to generate the same mass neutron star. These facts suggest that an accurate measurement of neutron-star radii may reduce the model-dependence (i.e., the spread) observed in Fig. 3. Yet, even without further constraints the spread is relatively small and a measurement of $`R_nR_p`$ in 208Pb may still prove decisive.

## IV Summary and Conclusions

Recent X-ray observations suggest that some neutron stars cool quickly. This enhanced cooling could arise from the direct URCA process—that requires a high proton fraction—or from the beta decay of additional hadrons in dense matter, such as pions, kaons, hyperons, or quarks. Yet, it seems unlikely that the X-ray observations alone will determine the origin of the enhanced cooling.

In this work we propose to use a laboratory experiment to constrain the direct URCA process in neutron stars. The Parity Radius Experiment at the Jefferson Laboratory prex ; bigpaper  aims to measure the neutron radius of 208Pb accurately and model independently via parity-violating electron scattering. For the direct URCA process to be realized, the equilibrium proton fraction in the star must be large. The equilibrium proton fraction is determined by the symmetry energy, whose density dependence can be strongly constrained through a measurement of the neutron radius in 208Pb. Such a measurement could provide a reliable extrapolation of the proton fraction to higher densities. Thus, predictions for the neutron radius in 208Pb have been correlated to the proton fraction in dense neutron rich matter by using a wide range of relativistic effective-field theory models. We find that models with a neutron skin in 208Pb of $`R_nR_p0.20`$ fm generate proton fractions that are too small to allow the direct URCA process in 1.4 $`M_{}`$ neutron stars. Conversely, if $`R_nR_p0.25`$ fm, then all models predict the URCA cooling of 1.4 $`M_{}`$ stars.

While this paper has focused on relativistic effective field-theory models, we expect our conclusions to be general and applicable to other approaches, both relativistic and nonrelativistic. For example, the nonrelativistic equation of state of Friedman and Pandharipande fp  predicts too small a proton fraction for URCA cooling to be possible. Moreover, this equation of state yields a neutron skin in 208Pb of only $`R_nR_p=0.16\pm 0.02`$ fm brown . Thus, these results are fully consistent with Fig. 3 that predicts no URCA cooling for such a small value of $`R_nR_p`$.

The equation of state considered in this work consists of matter composed of neutrons, protons, electrons, and muons in beta equilibrium; no exotic component was invoked. Further, no explicit proton or neutron pairing was considered. Nucleon superfluidity is an accepted phenomenon in nuclear physics and superfluid gaps are important for the cooling of neutron stars yakovlev . Thus, the study of pairing gaps in relativistic effective-field theories is an important area of future work; first steps in this direction have been taken in Ref. super . In particular, the proton pairing gap in matter with a high proton concentration must be computed t5 .

In summary, the feasibility of enhanced cooling of neutron stars via the direct URCA process was studied by correlating the proton fraction in dense, neutron-rich matter to the neutron skin of 208Pb. Thus, a measurement of the neutron radius in 208Pb may become vital for confirming (or dismissing) the direct URCA cooling of neutron stars. If direct URCA cooling is ruled out, then observations of enhanced cooling may provide strong evidence in support of exotic states of matter, such as meson condensates and quark matter, at the core of neutron stars.

###### Acknowledgements.

This work was supported by the U.S. DOE under contracts DE-FG02-87ER40365 and DE-FG05-92ER40750.


format: text

text:
# A model of nonlinear evolution and saturation of the turbulent MHD dynamo

## 1 Introduction

It has long been appreciated  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 .

If the magnetic Prandtl number ($`Pr=\nu /\eta `$, 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 $`Pr^{1/2}`$ below the viscous scale of the fluid. This situation is known as the Batchelor regime  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 ). 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 ( 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 ). 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 , 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 turbulence — a convenient starting point for the operation of the helical mean-field dynamo  or of the helical inverse cascade . Note that we believe the inverse cascade of magnetic helicity that takes place in helical MHD turbulence  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 ).

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., ), the fluctuation dynamo is . 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-$`Pr`$ 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-$`Pr`$ 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 magnetic-energy 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.

## 2 The kinematic dynamo

### 2.1 The growth of magnetic energy

The kinematic dynamo is described by the induction equation

$`_t\mathrm{B}\mathbf{+}\mathrm{u}\mathbf{}\mathbf{}\mathrm{B}\mathbf{=}\mathrm{B}\mathbf{}\mathbf{}\mathrm{u}\mathbf{+}𝜼\mathrm{\Delta }\mathrm{B}\mathbf{,}`$ (1)
where $`\mathrm{B}`$ is the magnetic field, $`\eta `$ is the magnetic diffusivity, and $`\mathrm{u}`$ is the random incompressible velocity field, whose statistics are externally prescribed. The following exact evolution law for the magnetic energy $`W=\frac{1}{2}B^2`$ is an immediate corollary of (1):

$`{\displaystyle \frac{}{t}}W=2\gamma \left(t\right)W2\eta k_{\mathrm{rms}}^2\left(t\right)W,`$ (2)
where, by definition,

$`\gamma \left(t\right)={\displaystyle \frac{\mathrm{B}\mathrm{B}\mathbf{:}\mathbf{}\mathrm{u}\mathbf{}}{\mathbf{}\mathrm{B}^2\mathbf{}}},k_{\mathrm{rms}}^2\left(t\right)={\displaystyle \frac{\left(_jB^i\right)\left(_jB^i\right)}{B^2}}={\displaystyle \frac{1}{W}}{\displaystyle _0^{\mathrm{}}}kk^2M(t,k),`$ (3)
and $`M(t,k)=\frac{1}{2}\mathrm{\Omega }k\left|\mathrm{B}\mathbf{(}\mathrm{t}\mathbf{,}\mathrm{k}\mathbf{)}\mathbf{\right|}^2\mathbf{}`$ 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  argued that, in a turbulent flow, $`\gamma \left(t\right)`$ would tend to some constant $`\overline{\gamma }`$, thus implying exponential growth of the magnetic energy. This is essentially a consequence of the line-stretching property of turbulence .

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 $`\overline{\gamma }`$ . In a large-$`Pr`$ 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: $`\overline{\gamma }k_\nu u_\nu `$, where $`k_\nu `$ and $`u_\nu `$ are the characteristic wave number and velocity of these eddies .

### 2.2 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  which uses a Gaussian white-in-time velocity as a stand-in for the real turbulence:

$`u^i(t,\mathrm{x}\mathbf{)}\mathrm{u}^\mathrm{j}\mathbf{(}\mathrm{t}^{\mathbf{}}\mathbf{,}\mathrm{x}^{\mathbf{}}\mathbf{)}=\delta \left(tt^{}\right)\kappa ^{ij}\left(\mathrm{x}\mathbf{}\mathrm{x}^{\mathbf{}}\mathbf{\right)}\mathbf{.}`$ (4)
On the basis of this model, it is possible to show that the magnetic-energy spectrum at subviscous scales ($`kk_\nu `$) satisfies the following equation

$`_tM={\displaystyle \frac{\overline{\gamma }}{5}}\left(k^2{\displaystyle \frac{^2M}{k^2}}2k{\displaystyle \frac{M}{k}}+6M\right)2\eta k^2M,`$ (5)
where $`\overline{\gamma }=\left(1/6\right)\left[\mathrm{\Delta }\kappa ^{ii}\left(\mathrm{y}\mathbf{\right)}\right]_{\mathrm{y}\mathbf{=}0}`$, which is consistent with its physical meaning as the characteristic eddy-turnover rate of the turbulence .

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 energy-evolution law (2). Indeed, we can rewrite equation (5) in a conservative Fokker–Planck form

$`_tM={\displaystyle \frac{}{k}}\left[D\left(k\right){\displaystyle \frac{M}{k}}V\left(k\right)M\right]+2\overline{\gamma }M2\eta k^2M,`$ (6)
where $`D\left(k\right)=\overline{\gamma }k^2/5`$ is the diffusion coefficient and $`V\left(k\right)=4\overline{\gamma }k/5`$ the “drift velocity” in $`k`$ space ($`V>0`$ implies the spectrum moving towards higher $`k`$). We see that $`W\left(t\right)=_0^{\mathrm{}}kM(t,k)`$ satisfies (2) with $`\gamma \left(t\right)=\overline{\gamma }`$.

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), $`\delta \left(tt^{}\right)`$ is replaced by $`\left(1/2\tau _\mathrm{c}\right)\mathrm{exp}\left(\left|tt^{}\right|/\tau _\mathrm{c}\right)`$, where $`\tau _\mathrm{c}`$ is the correlation time of the velocity. Doing a calculation parallel to that in , one gets

$`\overline{\gamma }={\displaystyle \frac{2\tau _\mathrm{c}}{3}}{\displaystyle _0^{\mathrm{}}}kk^2E\left(k\right),`$ (7)
where $`E\left(k\right)`$ is the energy spectrum of the fluid turbulence. \[For a more refined model, one might allow for the decorrelation rate $`\tau _\mathrm{c}^1`$ 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 $`\tau _\mathrm{c}`$ characteristic of the fastest (highest-$`k`$) eddies, which dominate the integral in (7), and $`\tau _\mathrm{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 $`\tau _{\mathrm{eddy}}^1`$ of the flow. It is clear that, in a fully developed turbulence, $`\tau _\mathrm{c}\tau _{\mathrm{eddy}}`$, so we may rewrite (7) as

$`\overline{\gamma }=c_1\left[{\displaystyle _0^{\mathrm{}}}kk^2E\left(k\right)\right]^{1/2},`$ (8)
where $`c_1`$ is a coefficient of order unity.

The other main assumption in (5) is that we are focusing on magnetic excitations at $`kk_\nu `$, 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 ). For $`Pr1`$, we shall find that $`k_{\mathrm{rms}}k_\nu `$, 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_{\mathrm{rms}}`$ across which there is no flux of magnetic energy.

At $`kk_\nu `$, 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)`$ , 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_\nu `$ 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 $`\overline{\gamma }`$ 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:

$`\overline{\gamma }\left(k\right)=c_1\left[{\displaystyle _0^k}k^{}k^2E\left(k^{}\right)\right]^{1/2}.`$ (9)
For $`k<k_\nu `$, 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 $`\overline{\gamma }\left(k\right)`$ in (9) has the same scaling as the $`k`$-dependent eddy-damping rate in the popular phenomenological EDQNM closure , where setting $`c_1=0.36`$ results in a good fit to the Kolmogorov constant . The shearing rate $`\overline{\gamma }\left(k\right)`$ 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 $`\overline{\gamma }\left(k_{}\right)`$ \[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 $`\overline{\gamma }`$ 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.

### 2.3 The diffusion-free (ideal-MHD) regime

Let us discuss how a magnetic excitation initially concentrated at some wave numbers $`k^{}`$ such that $`k_\nu k^{}k_\eta `$ will evolve with time. Neglecting the resistive-diffusion term for the time being ($`\eta =0`$), we immediately find the Green’s-function solution of (5):

$`M(t,k)=\mathrm{}^{\left(3/4\right)\overline{\gamma }t}{\displaystyle _0^{\mathrm{}}}{\displaystyle \frac{k^{}}{k^{}}}M_0\left(k^{}\right)\left({\displaystyle \frac{k}{k^{}}}\right)^{3/2}{\displaystyle \frac{1}{\sqrt{\left(4/5\right)\pi \overline{\gamma }t}}}\mathrm{exp}\left({\displaystyle \frac{\left[\mathrm{ln}\left(k/k^{}\right)\right]^2}{\left(4/5\right)\overline{\gamma }t}}\right),`$ (10)
where $`M_0\left(k^{}\right)`$ is the initial spectrum. This solution shows that (i) the amplitude of each Fourier mode grows exponentially in time at the rate $`\left(3/4\right)\overline{\gamma }`$, (ii) the number of excited modes \[i.e., the width of the lognormal envelope in (10)\] also grows exponentially at the rate $`\left(4/5\right)\overline{\gamma }`$, (iii) the peak of the excitation originating from each initially present mode $`k^{}`$ moves towards larger $`k`$: $`k_{\mathrm{peak}}=k^{}\mathrm{exp}\left[\left(3/5\right)\overline{\gamma }t\right]`$, leaving a power spectrum $`k^{3/2}`$ behind. The cumulative effect is the total magnetic-energy growth at the rate $`2\overline{\gamma }`$.

The solution (10) becomes invalid when the magnetic excitation reaches the velocity scale $`k_\nu ^1`$ and/or the resistive scale $`k_\eta ^1`$. In the former case, including the $`k`$ dependence of $`\overline{\gamma }`$ 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.

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

$`M(t,k)=\mathrm{}^{\lambda \overline{\gamma }t}\mathrm{\Phi }\left(\xi \right),\xi =k/k_\eta ,k_\eta =\left(\overline{\gamma }/10\eta \right)^{1/2},`$ (11)
we get

$`\xi ^2\mathrm{\Phi }^{\prime \prime }2\xi \mathrm{\Phi }^{}+\left(65\lambda \right)\mathrm{\Phi }\xi ^2\mathrm{\Phi }=0.`$ (12)
The solutions of (12) are Bessel functions. Demanding that $`\mathrm{\Phi }\left(\xi \mathrm{}\right)=0`$, we find

$`\mathrm{\Phi }\left(\xi \right)=\mathrm{const}\xi ^{3/2}K_{\nu \left(\lambda \right)}\left(\xi \right),\nu \left(\lambda \right)=\sqrt{5\left(\lambda 3/4\right)},`$ (13)
where $`K_\nu `$ is the Macdonald function.

The eigenvalue $`\lambda `$ 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 $`Pr1`$ (see, e.g.,  and the exhaustive reference list in ). Here we choose some finite IR cutoff $`k_{}`$ and impose a zero-flux boundary condition at $`k=k_{}`$. This gives \[see equation (6)\]

$`\xi _{}\mathrm{\Phi }^{}\left(\xi _{}\right)4\mathrm{\Phi }\left(\xi _{}\right)=0,\xi _{}=k_{}/k_\eta .`$ (14)
Substituting (13) yields the following transcendental equation for $`\lambda =\lambda \left(\xi _{}\right)`$:

$`\xi _{}K_{\nu \left(\lambda \right)}^{}\left(\xi _{}\right){\displaystyle \frac{5}{2}}K_{\nu \left(\lambda \right)}\left(\xi _{}\right)=\left(\nu \left(\lambda \right){\displaystyle \frac{5}{2}}\right)K_{\nu \left(\lambda \right)}\left(\xi _{}\right)\xi _{}K_{\nu \left(\lambda \right)+1}\left(\xi _{}\right)=0.`$ (15)
It is immediately clear that (15) has no solutions with real $`\nu \left(\lambda \right)`$, i.e., with $`\lambda 3/4`$. Indeed, such solutions could only exist if the first term were positive: $`\nu \left(\lambda \right)>5/2`$, which would imply $`\lambda >2`$. Since the total energy cannot grow faster than at the rate $`2\overline{\gamma }`$, $`\lambda >2`$ is impossible. Now, the Macdonald function $`K_{\nu \left(\lambda \right)}\left(\xi \right)`$ of imaginary order $`\nu \left(\lambda \right)=\mathrm{}\stackrel{~}{\nu }\left(\lambda \right)`$ oscillates as it approaches $`\xi =0`$. Since the spectrum must be a nonnegative function, all zeros of $`K_{\mathrm{}\stackrel{~}{\nu }\left(\lambda \right)}\left(\xi \right)`$ must lie to the left of the IR cutoff $`\xi _{}`$:

$`z_{\mathrm{}\stackrel{~}{\nu }\left(\lambda \right),m}\xi _{}.`$ (16)
As we are considering a large-$`Pr`$ case, we must take $`\xi _{}=k_{}/k_\eta Pr^{1/2}1`$. The condition (16) then implies $`\stackrel{~}{\nu }\left(\lambda \right)1`$. In these limits, (15) reduces to

$`K_{\mathrm{}\stackrel{~}{\nu }\left(\lambda \right)}\left(\xi _{}\right){\displaystyle \frac{1}{\stackrel{~}{\nu }\left(\lambda \right)}}\mathrm{sin}\left[\stackrel{~}{\nu }\left(\lambda \right)\mathrm{ln}\left(\xi _{}/2\right)\right]=0,\stackrel{~}{\nu }\left(\lambda \right)=\sqrt{5\left(3/4\lambda \right)}.`$ (17)
Thus, in the large-$`Pr`$ limit, the zero-flux boundary condition (14) is equivalent to requiring that the spectrum should vanish at the IR cutoff: $`\mathrm{\Phi }\left(\xi _{}\right)=0`$ (cf. ). Taking the first (rightmost) zero in (17), so that (16) is observed, we find the desired eigenvalue

$`\lambda {\displaystyle \frac{3}{4}}{\displaystyle \frac{\pi ^2}{5\left[\mathrm{ln}\left(\xi _{}/2\right)\right]^2}}={\displaystyle \frac{3}{4}}\mathrm{O}\left({\displaystyle \frac{1}{\mathrm{ln}^2\left(Pr^{1/2}\right)}}\right).`$ (18)
For $`\xi \xi _{}`$, the spectrum is well described by (13) with $`\nu \left(\lambda \right)=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_\nu kk_\eta `$ \[see (13)\], and all modes grow exponentially at the same rate $`\left(3/4\right)\overline{\gamma }`$ . Most importantly, the bulk of the magnetic energy is concentrated at the resistive scales $`k_\eta ^1`$.

## 3 A phenomenological theory of the nonlinear dynamo

### 3.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 $`\mathrm{u}`$ can no longer be considered as given. It has to be determined self-consistently from the Navier–Stokes equation:

$`_t\mathrm{u}\mathbf{+}\mathrm{u}\mathbf{}\mathbf{}\mathrm{u}\mathbf{=}𝝂\mathrm{\Delta }\mathrm{u}\mathbf{}\mathbf{}\mathrm{p}\mathbf{+}\mathrm{B}\mathbf{}\mathbf{}\mathrm{B}\mathbf{+}\mathrm{f}\mathbf{,}`$ (19)
where $`\nu `$ is the fluid viscosity, $`\mathrm{f}`$ is some random large-scale forcing, $`p`$ is the pressure determined from the incompressibility condition $`\mathbf{}\mathrm{u}\mathbf{=}0`$, and $`p`$ and $`\mathrm{B}`$ have been normalized by $`\rho `$ and $`\left(4\pi \rho \right)^{1/2}`$, respectively, $`\rho `$ 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: $`\mathrm{u}\mathbf{}\mathbf{}\mathrm{u}\mathbf{}\mathrm{B}\mathbf{}\mathbf{}\mathrm{B}`$. In order to correctly estimate the tension force $`\mathrm{B}\mathbf{}\mathbf{}\mathrm{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 $`\mathrm{B}\mathbf{}\mathbf{}\mathrm{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 , 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, $`\mathrm{B}\mathbf{}\mathbf{}\mathrm{B}\mathbf{}\mathrm{k}_𝝂\mathrm{B}^2`$. Since the inertial term is $`\mathrm{u}\mathbf{}\mathbf{}\mathrm{u}\mathbf{}\mathrm{k}_𝝂\mathrm{u}_𝝂^2`$, the nonlinearity becomes important when the small-scale magnetic energy equalizes with the energy of the viscous-scale eddies: $`B^2u_\nu ^2`$, or $`WW_\nu `$ . Note that if the folded structure of the field had been ignored, then estimating $`\mathrm{B}\mathbf{}\mathbf{}\mathrm{B}\mathbf{}\mathrm{k}_𝜼\mathrm{B}^2`$ would have led to the incorrect conclusion that nonlinearities become important at the much lower magnetic energies: $`B^2u_\nu ^2/Pr^{1/2}`$.

### 3.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 $`\mathrm{}_\mathrm{s}>k_\nu ^1`$, the folds are elongated accordingly (see remark in section 4.3), and the tension force is $`\mathrm{B}\mathbf{}\mathbf{}\mathrm{B}\mathbf{}\mathrm{B}^2\mathbf{/}\mathbf{}_𝐬`$. The corresponding inertial term is $`\mathrm{u}\mathbf{}\mathbf{}\mathrm{u}\mathbf{}\mathrm{u}_\mathbf{}_𝐬^2\mathbf{/}\mathbf{}_𝐬`$, so these eddies become suppressed when $`B^2u_\mathrm{}_\mathrm{s}^2`$, 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 $`\mathrm{}_\mathrm{s}\left(t\right)`$ at a given time $`t`$ by

$`W\left(t\right){\displaystyle \frac{1}{2}}u_{\mathrm{}_\mathrm{s}\left(t\right)}^2.`$ (20)
In equation (2), $`\gamma \left(t\right)`$ is now the turnover rate of the eddies of scale $`\mathrm{}_\mathrm{s}\left(t\right)`$: $`\gamma \left(t\right)u_{\mathrm{}_\mathrm{s}\left(t\right)}/\mathrm{}_\mathrm{s}\left(t\right)`$. Then $`\gamma \left(t\right)W\left(t\right)u_{\mathrm{}_\mathrm{s}\left(t\right)}^3/\mathrm{}_\mathrm{s}\left(t\right)ϵ`$, where $`ϵ=\mathrm{const}`$ is the Kolmogorov energy flux. Equation (2) becomes

$`{\displaystyle \frac{}{t}}W\chi ϵ2\eta k_{\mathrm{rms}}^2\left(t\right)W,`$ (21)
where $`\chi `$ 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 $`\mathrm{}_\mathrm{s}\left(t\right)`$ where a finite fraction $`\chi `$ 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\left(t\right)ϵt`$. Furthermore, comparing (21) and (2), we get

$`\gamma \left(t\right){\displaystyle \frac{\chi ϵ}{2W\left(t\right)}}{\displaystyle \frac{\alpha }{t}},`$ (22)
where $`\alpha `$ 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 $`\mathrm{}_\mathrm{s}\left(t\right)`$ with direction reversals perpendicular to $`\mathrm{B}`$ at the resistive scale. Note that the resistive scale now changes with time: a straightforward estimate gives $`k_\eta \left(t\right)\left[\gamma \left(t\right)/\eta \right]^{1/2}\left(\eta t\right)^{1/2}`$. Thus, the resistive cutoff moves to larger scales. Indeed, as the stretching slows down with decreasing $`\gamma \left(t\right)`$, 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 $`\mathrm{}_\mathrm{s}\left(t\right)`$ and down to the viscous cutoff. These are the motions that do not amplify the magnetic field. They are waves of Alfvénic nature propagating along the folds of direction-reversing magnetic fields (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 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 $`\mathrm{}_\mathrm{s}k_f^1`$, 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\left(k_fu_0\right)^1\overline{\gamma }^1Re^{1/2}`$, where $`u_0`$ is the velocity of these eddies. The resistive cutoff scale is now larger than during the kinematic stage: $`k_\eta \left(t_0\right)/k_\eta ^{\left(\mathrm{kin}\right)}\left(t_0\overline{\gamma }\right)^{1/2}Re^{1/4}`$, where by $`k_\eta ^{\left(\mathrm{kin}\right)}`$ we denote the resistive cutoff in the kinematic regime. However, this is not a very significant shift: for example, in the interstellar medium, where $`Re10^4`$, it is by a meager factor of $`0.1`$. Recall that $`k_\nu /k_fRe^{3/4}`$ and $`k_\eta ^{\left(\mathrm{kin}\right)}/k_\nu Pr^{1/2}`$, so $`k_\eta \left(t_0\right)/k_f\left(RePr\right)^{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_f^1`$ with direction reversals at the resistive scale.

### 3.3 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), $`\gamma k_fu_0`$). As there are no scales in the system larger than $`k_f^1`$, there can be no further growth of the magnetic energy. Indeed, the tension in the folds is now $`\mathrm{B}\mathbf{}\mathbf{}\mathrm{B}\mathbf{}\mathrm{k}_\mathrm{f}\mathrm{W}`$ and the force balance in (19) requires $`\mathrm{B}\mathbf{}\mathbf{}\mathrm{B}\mathbf{}\mathrm{u}\mathbf{}\mathbf{}\mathrm{u}\mathbf{}\mathrm{k}_\mathrm{f}\mathrm{W}_0`$, so $`WW_0`$. Note that setting $`W/t=0`$ in (2) gives $`\gamma \left(t\right)=\eta k_{\mathrm{rms}}^2\left(t\right)`$. This is a reflection of the fact that the resistive cutoff $`k_\eta k_{\mathrm{rms}}`$ is now determined by the zero-sum balance between the field amplification (stretching) and the resistive dissipation:

$`\gamma \left(t\right)\eta k_\eta ^2\left(t\right).`$ (23)
However, at scales larger than the resistive cutoff, i.e., for $`k<k_\eta `$, the field-amplification rate $`\gamma `$ is larger than the resistive-dissipation rate $`\eta k^2`$ (figure 1, right panel). Since the total magnetic energy is not allowed to grow anymore, the growth of modes with $`k<k_\eta `$ must be compensated by the decrease in the value of $`k_\eta `$, i.e., by further selective decay of the high-$`k`$ modes. This, in turn, should lead to further decrease of $`\gamma `$ 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 $`\gamma \left(t\right)1/t`$ and $`k_\eta \left(t\right)\left(\eta t\right)^{1/2}`$.

Thus, we expect that this last stage of the nonlinear dynamo should be characterized by $`\gamma `$ dropping below the turnover rate $`k_fu_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.

### 3.4 The fully developed isotropic MHD turbulence

We saw in the preceding sections that $`k_\eta `$ 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_\eta `$ 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_\eta \left(t\right)`$ is basically a consequence of the balance between the field-amplification and the resistive-dissipation terms in the energy equation (2): as $`\gamma \left(t\right)`$ drops, so does $`k_\eta \left(t\right)k_{\mathrm{rms}}\left(t\right)`$; once the characteristic scale of the small-scale fields reaches the viscous scale $`k_\nu ^1`$, the Alfvénic turbulence will start to affect $`k_{\mathrm{rms}}`$ in an essential way: since the waves are damped at the viscous scale, $`k_{\mathrm{rms}}`$ must stabilize at $`k_{\mathrm{rms}}k_\nu `$. 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 , 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_{\mathrm{rms}}k_\nu `$, 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_\eta \left(k_\nu \right)\left(\eta k_\nu ^2\right)^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_\eta \left(k_\nu \right)t_0`$, which requires $`PrRe^{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_\eta \left(t_0\right)/k_\nu Pr^{1/2}Re^{1/4}1`$. This constitutes the “true large-$`Pr`$ regime,” which is the one relevant for astrophysical plasmas (thus, in the interstellar medium, $`Re10^4`$, $`Pr10^{14}`$). In this regime, the magnetic energy saturates at the equipartition level, $`WW_0`$. From the energy balance (2), we get an estimate of the amount of turbulent power that is still dissipated resistively: $`\chi ϵ\gamma W\eta k_\nu ^2W_0ϵRe^{1/2}Pr^1`$. Thus, $`\chi 1`$, i.e., most of the injected power now goes into the Alfvénic motions (which are dissipated viscously).

The other possibility is $`PrRe^{1/2}`$. In this case, the nonlinear-growth stage described in section 3.2 is curtailed with magnetic energy at a subequipartition value:

$`W{\displaystyle \frac{t_\eta \left(k_\nu \right)}{t_0}}W_0{\displaystyle \frac{Pr}{\sqrt{Re}}}W_0,`$ (24)
and resistivity continuing to take a significant ($`\chi `$ 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-$`Pr`$ 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 $`PrRe^{1/2}1`$. A busy reader with no time for these slow-paced developments may skip to section 5.

## 4 The Fokker–Planck model of the dynamo

### 4.1 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 $`\overline{\gamma }`$ by a time-dependent function $`\gamma \left(t\right)`$, for which we stipulate the following ad hoc expression:

$`\gamma \left(t\right)=c_1\left[{\displaystyle _0^{k_\mathrm{s}\left(t\right)}}kk^2E\left(k\right)\right]^{1/2},c_2{\displaystyle _{k_\mathrm{s}\left(t\right)}^{\mathrm{}}}kE\left(k\right)=W\left(t\right),`$ (25)
where $`c_1`$ and $`c_2`$ are some numerical constants of order unity, $`E\left(k\right)`$ is the hydrodynamic energy spectrum in the absence of magnetic field, $`k_\mathrm{s}\left(t\right)`$ is defined by the second expression in (25), and $`W\left(t\right)`$ is the total magnetic energy at time $`t`$. The Fokker–Planck model we propose for the nonlinear evolution of the magnetic spectrum is then

$`_tM={\displaystyle \frac{\gamma \left(t\right)}{5}}\left(k^2{\displaystyle \frac{^2M}{k^2}}2k{\displaystyle \frac{M}{k}}+6M\right)2\eta k^2M.`$ (26)
We now assume the following crude model profile of the hydrodynamic energy spectrum of the Kolmogorov turbulence in the absence of magnetic field:

$`E\left(k\right)=\{\begin{array}{cc}C_\mathrm{K}ϵ^{2/3}k^{5/3}\hfill & \text{for }k[k_f,k_\nu ]\text{,}\hfill \\ 0\hfill & \text{elsewhere,}\hfill \end{array}`$ (27)
where $`C_\mathrm{K}1.5`$ is the Kolmogorov constant. The value of $`k_\nu `$ is set to enforce consistency of the model profile (27) with the exact relation $`ϵ=2\nu _0^{\mathrm{}}kk^2E\left(k\right)`$:

$`k_\nu =k_f\left(1+{\displaystyle \frac{ϵ^{1/3}}{\left(3/2\right)C_\mathrm{K}\nu k_f^{4/3}}}\right)^{3/4}Re^{3/4}k_f.`$ (28)
Using (27), it is straightforward to derive from (25)

$`\gamma \left(t\right)=\overline{\gamma }\left[1{\displaystyle \frac{1}{\left(1+W_0/W_\nu \right)^2}}\right]^{1/2}\left[{\displaystyle \frac{1}{\left(1+W\left(t\right)/W_\nu \right)^2}}{\displaystyle \frac{1}{\left(1+W_0/W_\nu \right)^2}}\right]^{1/2},`$ (29)
where $`W_\nu =\frac{3}{2}c_2C_\mathrm{K}ϵ^{2/3}k_\nu ^{2/3}`$ and $`W_0=c_2_0^{\mathrm{}}kE\left(k\right)`$ are of the order of the energies of the viscous- and of the outer-scale eddies, respectively. Here, $`\overline{\gamma }=c_1\left(ϵ/2\nu \right)^{1/2}`$ follows from evaluating (25) for $`W0`$, $`k_\mathrm{s}=k_\nu `$. Note that

$`{\displaystyle \frac{W_0}{W_\nu }}=\left({\displaystyle \frac{k_\nu }{k_f}}\right)^{2/3}1Re^{1/2}1.`$ (30)
The kinematic regime $`\gamma \left(t\right)\overline{\gamma }`$ (section 2), the nonlinear-growth regime (22) (section 3.2), and the long-time asymptotic $`\gamma \left(t\right)0`$ (section 3.3) are all contained in the model expression (29).

In the above formulae, $`E\left(k\right)`$ 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\left(k\right)`$ remains unaffected by the magnetic fields at $`k<k_\mathrm{s}\left(t\right)`$. As was explained in section 3.2, in the suppressed part of the inertial range $`k>k_\mathrm{s}\left(t\right)`$, 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).

### 4.2 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 $`\gamma \left(t\right)`$. 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 $`\gamma \left(t\right)`$ can be represented as

$`\gamma \left(t\right)={\displaystyle \frac{\alpha }{t}},`$ (31)
where $`\alpha `$ is a numerical constant, equation (26) admits a self-similar solution

$`M(t,k)={\displaystyle \frac{W\left(t\right)}{k_\eta \left(t\right)}}\mathrm{\Phi }\left(\xi \right),\xi =k/k_\eta \left(t\right),{\displaystyle _0^{\mathrm{}}}\xi \mathrm{\Phi }\left(\xi \right)=1,`$ (32)
$`k_\eta \left(t\right)=\left[\gamma \left(t\right)/10\eta \right]^{1/2}=\left(\alpha /10\eta t\right)^{1/2},`$ (33)
$`W\left(t\right)=Ct^\mu ,`$ (34)
where the function $`\mathrm{\Phi }\left(\xi \right)`$ satisfies

$`\xi ^2\mathrm{\Phi }^{\prime \prime }\left(2+{\displaystyle \frac{5}{2}}\alpha ^1\right)\xi \mathrm{\Phi }^{}+\left(65{\displaystyle \frac{2\mu +1}{2}}\alpha ^1\right)\mathrm{\Phi }\xi ^2\mathrm{\Phi }=0.`$ (35)
Clearly, for $`\gamma \left(t\right)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 $`\mu `$ in (34) must be determined from the assumptions about the physics of the problem. The value of $`\alpha `$ is fixed by the solvability of equation (35).

### 4.3 The first self-similar regime: nonlinear growth

In section 3.2, we argued that $`W\left(t\right)`$ 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 $`\mu =1`$. Then, by virtue of (22), $`\gamma \left(t\right)`$ 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

$`W\left(t\right)\alpha ^1\overline{\gamma }W_\nu t{\displaystyle \frac{\chi ϵ}{2\alpha }}t,\overline{\gamma }^1t\overline{\gamma }^1{\displaystyle \frac{W_0}{W_\nu }}\overline{\gamma }^1Re^{1/2}t_0.`$ (36)
Comparing the above asymptotic with equation (21), we note that $`1/2\alpha `$ 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 $`\alpha `$, we have to solve the equation (35) with $`\mu =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:

$`\mathrm{\Phi }\left(\xi \right)=\mathrm{const}\xi ^{\sigma \left(\alpha \right)}K_{\nu \left(\alpha \right)}\left(\xi \right),\sigma \left(\alpha \right)={\displaystyle \frac{3}{2}}+{\displaystyle \frac{5}{4}}\alpha ^1,`$ (37)
$`\nu \left(\alpha \right)={\displaystyle \frac{5}{4}}\left({\displaystyle \frac{1}{\alpha }}{\displaystyle \frac{1}{\alpha _+}}\right)^{1/2}\left({\displaystyle \frac{1}{\alpha }}{\displaystyle \frac{1}{\alpha _{}}}\right)^{1/2},{\displaystyle \frac{1}{\alpha _\pm }}={\displaystyle \frac{2}{5}}\left(\pm 4\sqrt{6}9\right).`$ (38)
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

$`M(t,k)={\displaystyle \frac{W\left(t\right)}{k_\eta \left(t\right)}}\mathrm{\Phi }(\xi ,\xi _{}),{\displaystyle _\xi _{}^{\mathrm{}}}\xi \mathrm{\Phi }(\xi ,\xi _{})=1,`$ (39)
where $`\xi _{}=k_{}/k_\eta \left(t\right)=k_{}\left(10\eta \alpha ^1t\right)^{1/2}`$ is a new time-dependent dimensionless variable. Since the non-self-similar solution (39) has a self-similar asymptotic when $`\xi _{}0`$, we can still use (37), provided $`\xi _{}1`$. This solution remains valid as long as $`t\alpha t_\eta `$, where $`t_\eta =\left(10\eta k_{}^2\right)^1`$ is the resistive time associated with the IR cutoff $`k_{}`$. In large-$`Pr`$ 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 $`\xi _{}=k_{}/k_\eta \left(t\right)`$. Substituting (37) into (14), we get the following equation for $`\alpha `$:

$`\left[\nu \left(\alpha \right)+\sigma \left(\alpha \right)4\right]K_{\nu \left(\alpha \right)}\left(\xi _{}\right)\xi _{}K_{\nu \left(\alpha \right)+1}\left(\xi _{}\right)=0.`$ (40)
As explained above, this equation must be considered in the limit $`\xi _{}1`$ and the resulting $`\alpha `$ must remain finite as $`\xi _{}0`$. Again there are no solutions with real $`\nu \left(\alpha \right)`$. For imaginary $`\nu \left(\alpha \right)=\mathrm{}\stackrel{~}{\nu }\left(\alpha \right)`$, the requirement that the spectrum should have no nodes to the right of $`\xi _{}1`$ again implies $`\stackrel{~}{\nu }\left(\alpha \right)1`$, and we get

$`K_{\mathrm{}\stackrel{~}{\nu }\left(\alpha \right)}\left(\xi _{}\right){\displaystyle \frac{1}{\stackrel{~}{\nu }\left(\alpha \right)}}\mathrm{sin}\left[\stackrel{~}{\nu }\left(\alpha \right)\mathrm{ln}\left(\xi _{}/2\right)\right]=0,\stackrel{~}{\nu }\left(\alpha \right)5^{1/2}6^{1/4}\left({\displaystyle \frac{1}{\alpha _+}}{\displaystyle \frac{1}{\alpha }}\right)^{1/2}.`$ (41)
The solution is (see figure 2)

$`\alpha ^1\alpha _+^1{\displaystyle \frac{\pi ^2}{5\sqrt{6}\left[\mathrm{ln}\left(\xi _{}/2\right)\right]^2}}0.32\mathrm{O}\left({\displaystyle \frac{1}{\mathrm{ln}^2\left(t/t_\eta \right)}}\right).`$ (42)
Since $`tt_\eta `$, the spectral profile is well described by (37) with $`\nu \left(\alpha \right)=0`$ and $`\sigma \left(\alpha \right)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_\eta \left(t\right)`$, which is moving leftwards as $`\left(\eta t\right)^{1/2}`$ in agreement with the estimate in section 3.2 (see figures 23). The turbulent power arriving into the smallest unsuppressed eddies is (partially) transfered into the small-scale magnetic fields. The fraction $`1/2\alpha `$ ($`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\overline{\gamma }^1Re^{1/2}`$, we recover the physical estimates made in section 3.2: $`W\left(t_0\right)Re^{1/2}W_\nu W_0`$ and $`k_\eta \left(t_0\right)Re^{1/4}k_\eta ^{\left(\mathrm{kin}\right)}\left(RePr\right)^{1/2}k_f`$.

Remark on the elongation of the folds. In section 3.2, we assumed the growth of the effective stretching scale $`\mathrm{}_\mathrm{s}\left(t\right)`$ 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 $`\mathrm{}_{}\left(t\right)\mathrm{}_\mathrm{s}\left(t\right)`$. This assumption can now be supported in the following way. Consider a straight segment of a field line of length $`\mathrm{}_{}`$. In the absence of resistivity, volume- and flux-conservation constraints in an incompressible flow imply $`\mathrm{}_{}/B\mathrm{const}`$, so $`\mathrm{log}\mathrm{}_{}/t=\left(1/2\right)\mathrm{log}W/t=\gamma \left(t\right)`$. 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 $`\alpha ^1\gamma \left(t\right)`$ in the first self-similar regime. However, the length of the field line continues to grow at the rate $`\gamma \left(t\right)`$, so $`\mathrm{log}\mathrm{}_{}/t=\alpha \mathrm{log}W/t`$. Using the nonlinear balance (20) and Kolmogorov scaling for velocity gives $`\mathrm{log}\mathrm{}_{}/\mathrm{log}\mathrm{}_\mathrm{s}=2\alpha /32`$. By itself, this would indicate that $`\mathrm{}_{}`$ grows even faster than $`\mathrm{}_\mathrm{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 $`\mathrm{}_{}`$ grows to be comparable to $`\mathrm{}_\mathrm{s}`$, it is limited to $`\mathrm{}_{}\mathrm{}_\mathrm{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 $`\alpha `$: it is sufficient that $`\alpha ^12/3`$.

### 4.4 The second self-similar regime: approach to saturation

Thus, after time $`t_0\overline{\gamma }^1Re^{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 $`\mu =0`$ in (34). Accordingly, substituting (31) into (29), we find

$`W\left(t\right)W_0\left[1{\displaystyle \frac{\alpha ^2}{2}}\left({\displaystyle \frac{t_0}{t}}\right)^2\right],tt_0=\overline{\gamma }^1{\displaystyle \frac{W_0}{W_\nu }}\overline{\gamma }^1Re^{1/2}.`$ (43)
We now set $`\mu =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 $`\alpha _\pm `$ in (38) and for $`\stackrel{~}{\nu }\left(\alpha \right)`$ in (41) are now as follows:

$`{\displaystyle \frac{1}{\alpha _\pm }}=2\left(\pm 2\sqrt{2/5}1\right),\stackrel{~}{\nu }\left(\alpha \right)5^{3/4}2^{1/4}\left({\displaystyle \frac{1}{\alpha _+}}{\displaystyle \frac{1}{\alpha }}\right)^{1/2}.`$ (44)
The new eigenvalue is

$`\alpha ^1\alpha _+^1{\displaystyle \frac{\pi ^2}{5\sqrt{5/2}\left[\mathrm{ln}\left(\xi _{}/2\right)\right]^2}}0.53\mathrm{O}\left({\displaystyle \frac{1}{\mathrm{ln}^2\left(t/t_\eta \right)}}\right).`$ (45)
For $`tt_\eta `$, the spectrum is given by (37) with $`\nu \left(\alpha \right)=0`$ and $`\sigma \left(\alpha \right)2.2`$.

Thus, in the second self-similar stage, the magnetic energy is nearly constant, but the resistive cutoff $`k_\eta \left(t\right)`$ continues to decrease as $`\left(\eta t\right)^{1/2}`$, so the magnetic energy is redistributed from small to larger scales (see figures 23). This is a very slow process, but it can, in principle, go on for a very long time: until self-similarity is broken at $`tt_\eta =\left(10\eta k_{}^2\right)^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.

### 4.5 Steady states at the equipartition energy and below

The self-similar solutions obtained above break down after a finite time $`t_\eta `$ due to the presence of the finite IR cutoff $`k_{}`$. Once the resistive cutoff $`k_\eta `$ 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_\eta k_\nu `$. 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, $`\gamma \left(t\right)\gamma \left(\mathrm{}\right)=\mathrm{const}`$, so the corresponding solution of equation (26) is simply the $`\lambda =0`$ eigenmode of the kinematic equation (5) with $`\overline{\gamma }`$ replaced by $`\gamma \left(\mathrm{}\right)`$: using equations (11) and (13), we get

$`M\left(k\right)=\mathrm{const}\xi ^{3/2}K_{i\sqrt{15}/2}\left(\xi \right),\xi =k/k_\eta \left(\mathrm{}\right),k_\eta \left(\mathrm{}\right)=\left[\gamma \left(\mathrm{}\right)/10\eta \right]^{1/2}.`$ (46)
The value of $`\gamma \left(\mathrm{}\right)`$ must be determined from the boundary condition at $`k_{}`$:

$`\xi _{}K_{i\sqrt{15}/2}^{}\left(\xi _{}\right){\displaystyle \frac{5}{2}}K_{i\sqrt{15}/2}\left(\xi _{}\right)=0,\xi _{}=k_{}/k_\eta \left(\mathrm{}\right).`$ (47)
This is an equation for $`\xi _{}`$, whose value in turn determines $`\gamma \left(\mathrm{}\right)=10\eta k_{}^2/\xi _{}^2`$. We have to pick the largest root of (47) to ensure that the spectrum (46) has no nodes. This gives $`\xi _{}0.59`$, which is order one, as expected. Substituting $`\gamma \left(\mathrm{}\right)`$ into (29), we get the magnetic energy $`W\left(\mathrm{}\right)`$ in the steady state:

$`{\displaystyle \frac{W\left(\mathrm{}\right)}{W_0}}=\left(1+{\displaystyle \frac{W_\nu }{W_0}}\right)\left\{\left[1+\left(\left(1+{\displaystyle \frac{W_0}{W_\nu }}\right)^21\right){\displaystyle \frac{\gamma ^2\left(\mathrm{}\right)}{\overline{\gamma }^2}}\right]^{1/2}\left(1+{\displaystyle \frac{W_0}{W_\nu }}\right)^1\right\},`$ (48)
which, using $`W_0/W_\nu Re^{1/2}1`$ and $`\gamma \left(\mathrm{}\right)/\overline{\gamma }=10\eta k_{}^2/\overline{\gamma }\xi _{}^2\left(k_{}/k_\eta ^{\left(\mathrm{kin}\right)}\right)^2\left(k_{}/k_\nu \right)^2Pr^1`$, can be rewritten as follows

$`{\displaystyle \frac{W\left(\mathrm{}\right)}{W_0}}\left[1+\left({\displaystyle \frac{k_{}}{k_\nu }}\right)^4{\displaystyle \frac{Re}{Pr^2}}\right]^{1/2}Re^{1/2}.`$ (49)
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_\eta k_\nu `$, 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_\eta `$). 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_\nu `$, the end of the self-similar regime and the stabilization of $`k_\eta `$ at a value $`k_\nu `$ 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\left(\mathrm{}\right)W_0`$ for $`PrRe^{1/2}`$ and the physical estimate (24) of the subequipartition saturation level for $`PrRe^{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.

## 5 Discussion

We have presented a physical scenario of the evolution of the nonlinear large-$`Pr`$ 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_\eta \left(\eta k_\nu ^2\right)^1`$. In the interstellar medium, the Spitzer  value of the magnetic diffusivity is $`\eta 10^7\mathrm{cm}^2/\mathrm{s}`$, $`k_\nu 10^{16}\mathrm{cm}^1`$ , so $`t_\eta 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-$`Pr`$ 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 $`Pr`$. 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 $`Pr`$. 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 . 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\left(\mathrm{}\right)<W_0`$. Subequipartion saturated states have indeed been reported in the numerical experiments cited above. Since the condition for the true asymptotic regime is $`PrRe^{1/2}`$ (section 3.4), numerical experiments with very high resolution are required to adequately study the large-$`Pr`$ 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 long-standing physical problem . 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 $`\mathbf{}\mathrm{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-$`Pr`$ MHD turbulence with a fixed uniform mean field . 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  or Goldreich–Sridhar  phenomenologies. In fact, medium-resolution numerical investigations  rather seem to support the final states with small-scale energy concentration even for $`Pr=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  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 . 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  tensor viscosity . 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 . This anisotropy is all the more important in view of the locally anisotropic (folding) structure of the magnetic field itself .
* 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. ). 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.,  and many recent references cited therein). The treatment of the kinematic-dynamo problem for a model compressible case can be found in . 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.

The authors wish to thank B. Chandran for several stimulating discussions. This work was supported by the UKAEA Agreement No. QS06992, the EPSRC Grant No. GR/R55344/01, and by the US-DOE Contract No. DE–AC02–76CHO3073.

## Appendix A Alfvén waves and folded fields

Defining $`H^{ij}=B^iB^j`$, we can write the MHD equations in the following tensor form:

$`_tu^i+u^k_ku^i=_ip+_kH^{ik}+\nu \mathrm{\Delta }u^i+f^i,_iu^i=0,`$ (50)
$`_tH^{ij}+u^k_kH^{ij}=H^{kj}_ku^i+H^{ik}_ku^j+\eta R^{ij},`$ (51)
where $`_k=/x^k`$, and $`R^{ij}=B^i\mathrm{\Delta }B^j+B^j\mathrm{\Delta }B^i`$.

As we explained in section 3, the stretching and folding of the magnetic field by the turbulent eddies of scale $`\mathrm{}_\mathrm{s}`$ lead to the field becoming organized in folds of length $`\mathrm{}_{}\mathrm{}_\mathrm{s}`$ with direction reversals within each fold at the scale $`\mathrm{}_{}k_\eta ^1`$ . Consider now some typical instance during the nonlinear stage of the dynamo when $`k_\nu ^1\mathrm{}_\mathrm{s}<k_f^1`$ (all stretching motions up to the scale $`\mathrm{}_\mathrm{s}`$ are suppressed). There is a clear scale separation in the problem: writing the magnetic-field tensor in the form $`H^{ij}=b^ib^jB^2`$, where $`b^i=B^i/B`$, we see that $`B^2`$ varies at scales $`\mathrm{}_{}k_\eta ^1k_\nu ^1`$ while the tensor $`b^ib^j`$ varies at scales $`\mathrm{}_{}\mathrm{}_\mathrm{s}k_\nu ^1`$ ($`b^i`$ flips its sign at scale $`\mathrm{}_{}`$, but this does not affect the tensor product $`b^ib^j`$). We can, therefore, average all fields over the subviscous scales: $`u^i_{\mathrm{ss}}=u^i`$, $`H^{ij}_{\mathrm{ss}}=b^ib^jB^2_{\mathrm{ss}}`$. In the suppressed part of the inertial range, i.e. at scales $`\mathrm{}`$ such that $`k_\nu ^1\mathrm{}\mathrm{}_\mathrm{s}`$, the tensor $`H^{ij}_{\mathrm{ss}}=b^ib^jB^2_{\mathrm{ss}}`$ can be considered constant. Perturbing around this value and $`u^i=0`$, we find that equations (50)-(51) (without the dissipation terms) admit wave-like solutions with the dispersion relation

$`\omega =\pm k_{}v_\mathrm{A},k_{}=\sqrt{k_ik_jb^ib^j},v_\mathrm{A}=\sqrt{B^2_{\mathrm{ss}}}.`$ (52)
These are Alfvén-like waves propagating along the folds of the direction-reversing magnetic fields, which are straight at scales below $`\mathrm{}_\mathrm{s}`$ (figure 4). These waves are a modification of the previously introduced magnetoelastic waves  with account taken of the folding structure of the small-scale magnetic fields.

In (52), we must have $`\mathrm{}_\mathrm{s}^1k_{}k_\nu `$, so the frequency $`\omega `$ is larger than the characteristic time scale of the growth of the small-scale magnetic energy: $`\omega =k_{}v_\mathrm{A}v_\mathrm{A}/\mathrm{}_\mathrm{s}W^{1/2}/\mathrm{}_\mathrm{s}u_\mathrm{}_\mathrm{s}/\mathrm{}_\mathrm{s}\gamma `$. Since the resistive-dissipation rate of the small-scale fields is $`\eta k_{\mathrm{rms}}^2\gamma `$ 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  or the Goldreich–Sridhar  phenomenologies. Further investigations of this and related issues will be published elsewhere.

## References