\documentstyle[12pt]{article}
\title{Self-similar dissipation of three-dimensional magnetohydrodynamic 
       turbulence}
\author{A.Brandenburg$^{1,2)}$ and V.Urpin$^{1,3)}$\\
     1) NORDITA, Blegdamsvej 17, \\
        DK-2100  Copenhagen \O, Denmark \\
     2) Department of Mathematics \\
        University of Newcastle \\
        Newcastle upon Tyne NE1 7RU, UK \\
     3) A.F.Ioffe Institute of Physics and Technology \\
        194021 St.Petersburg, Russia}

\textwidth=16cm
\textheight=22cm

\begin{document}

\baselineskip15pt

\maketitle

%\newpage
\begin{abstract}

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx

{\bf Key words}: xxxxxxxxxxxxxxxxxx 

\end{abstract}

\newpage
\section{Introduction}

In many cases, especially in astrophysics, plasma is characterized by
turbulent magnetic fields. For instance, star-forming clouds of
interstellar gas can maintain sub- or trans-Alfvenic turbulence which is
thought to decay relatively slowly (Arons \& Max 1975, Gammie \&
Ostriker 1996). The decay of turbulence determines the most important
characteristics of these clouds such as their lifetime (Mac Low, Klessen
\& Burkert 1998). Turbulent magnetic fields can also play an important
role in observational manifestations of various types of stars. For 
example, recent models relate a superstrong magnetic field measured in 
some neutron stars ($\sim 10^{14}-10^{15}$ G) with the field generated
by turbulent dynamo action during the very early evolutionary stage
(the first $\sim20$ seconds!)
of these objects when they are convective (Thompson \& Duncan 1993). The 
decay of MHD turbulence after the end of the convective phase is here the
crucial issue for the further evolution of such neutron stars.    

Usually, the conductivity of astrophysical plasmas, $\sigma$, is rather 
high and therefore the corresponding magnetic decay time scale, $\sim
\sigma L^{2}/c^{2}$ where $L$ is the characteristic length scale, is very 
long. However, in many cases the viscosity, $\nu$, can be much larger than 
the magnetic diffusivity, $\eta = c^{2}/ 4 \pi \sigma$, providing therefore
an important
mechanism of dissipation of magnetohydrodynamic (MHD) turbulence. The decay 
of turbulent magnetic fields for large magnetic Prandtl number
$\mbox{Pr}_{\rm M}=\nu/\eta$ may have 
a number of qualitative features, and the present Letter is addressed this
problem.

The decay of MHD turbulence has been the subject of computational modelling 
in many papers. Numerical studies have mostly been concentrated on 
two-dimensional simulations, where high Reynolds numbers can be reached
more readily (see Kinney, McWilliams \& Tajima 1995, Galtier, Politano \& 
Pouquet 1997). Computational models suggest that within a wide range of
parameters the energy of MHD turbulence can decay in accordance with a 
power law, $ \propto t^{-\gamma}$ (see, e.g., Biskamp 1994). For 
incompressible MHD turbulence, the decay behavior seems to be consistent 
with $2/3 \leq \gamma \leq 1$ (Biskamp \& Welter 1984, Hatory 1984, 
Kinney, McWilliams \& Tajima 1995). Recent 3D computations of supersonic 
and super-Alfvenic turbulence also indicates that the decay can represented
by a power law with a somewhat larger exponent, $1.2 \geq \gamma \geq 0.85$
(Mac Low et al. 1998). The detailed numerical analysis of freely decaying 
incompressible 3D turbulence done by Biskamp \& M\"uller (1999) show that 
the decay law may be more complex and can depend on the magnetic helicity. 
However, this paper analyses the decay only for one value of the magnetic 
Prandl number, $\mbox{Pr}_{\rm M}=1$.

In the present Letter we present the results of a study of freely 
decaying 3D MHD turbulence for magnetic Prandl numbers ranging from
0.1 to 200. We discuss the decay 
laws of the integral quantities, in particular, the magnetic and kinetic
energies, $E_{\rm M}$ and $E_{\rm K}$, respectively. The energy decay is found to
follow a simple law which is almost independent on the initial conditions.  


\section{Equations and numerical results}

Consider the evolution of the magnetic field in a fluid with uniform
dynamical viscosity, $\mu$, and magnetic diffusivity, $\eta$. The 
continuity, momentum and induction equations are
\begin{equation}
\frac{\partial \rho}{\partial t} + \nabla \cdot (\rho {\bf v}) = 0,
\end{equation}
\begin{equation}
\frac{\partial {\bf v}}{\partial t} + ({\bf v} \cdot \nabla) {\bf v} =
-\frac{\nabla p}{\rho} + \frac{1}{4 \pi \rho} (\nabla \times {\bf B})
\times {\bf B} + {\mu\over\rho}\left(\Delta {\bf v}
+{\textstyle{1\over3}}\nabla\nabla\cdot{\bf v}\right),
\end{equation}
\begin{equation}
\frac{\partial {\bf B}}{\partial t} = \nabla\times({\bf v} \times {\bf B})
+ \eta \Delta {\bf B},
\end{equation}
where $p$ and $\rho$ are pressure and density, respectively. An isothermal
dependence of $p$ on $\rho$ is adopted, $p = c_{\rm s}^2\rho$.
We assume that initially, at $t=0$, there exists a magnetic field with 
the initial distribution ${\bf B}_{0}({\bf r}) \neq 0$.
% maintained by the corresponding distribution of currents.
%AB: the above is misleading. It's not really maintained, its decaying.
For the sake of simplicity, we
also assume that there are no fluid motions at $t=0$, but fluid motion
can be generated by the Lorentz force at $t>0$. We are interested in the 
behaviour of the magnetic and kinetic energies, $E_{\rm M}= \int (B^{2}/8 \pi) dV$
and $E_{\rm K}= \int (\rho v^{2}/2) dV$, respectively; integration is over the volume.   

To gain a basic idea about the behaviour of the magnetic field we first
solve Equations (1)-(3) numerically with an explicit code using sixth
order centered differences to compute spatial derivatives and a third
order Runge-Kutta scheme for the timestep. Further details can be found
in Brandenburg (2000). The computational domain is $0 < (x, y, z) < L$,
where $L=2\pi$, with either $120^3$ or $240^3$ meshpoints. We adopt
random helical initial conditions for ${\bf B}=\nabla\times{\bf A}$
by calculating the magnetic vector potential as
\begin{equation}
{\bf A}={\bf\tilde A}+h_{\rm M}\nabla\times{\bf\tilde A},
\end{equation}
where $h_{\rm M}=0.3$, unless specified otherwise, and ${\bf\tilde A}$
is gaussian distributed vector field that is smoothed such that the
peak in the initial powerspectrum of ${\bf B}$ is at wavenumber $k=20$
(for $120^3$ meshpoints). For smaller wavenumbers the magnetic
powerspectrum has a $k^4$ powerlaw.

The results of computations are plotted in Figures 1 and 2. In Fig.1, we 
show the time dependence of $E_{\rm M}(t)/E_{\rm M}(0)$ and $E_{\rm K}(t)/E_{\rm M}(0)$.
These dependences are shown for different values of $\nu$ and $\eta$; we normalize 
$\nu$ and $\eta$ to $c_{A} L$ where $c_{A}= \sqrt{2 E_{\rm M}(0) / \rho L^{3}}$ 
is the characteristic Alfven velocity at $t=0$. The decay curves show a 
remarkably similar behaviour for a wide range of parameters. The Lorentz 
force leads to a very fast generation of hydrodynamic motions which occurs 
on a hydrodynamic time scale, $\sim \ell / c_{A}$, where $\ell$ is the 
characteristic length scale of a random initial magnetic field. After a 
relatively short initial phase which lasts typically not longer than $\sim 
0.1 L/c_{A}$, both $b_{rms}$ and $u_{rms}$ reach a self-similar regime and 
evolve further coherently. During this phase, the ratio $E_{\rm K}/E_{\rm M}$ 
is approximately constant, and the decay follows a power law with a high 
accuracy. The exponent depends weakly on $\nu$ and $\eta$ being larger for 
larger dissipative coefficients. The fact that $b_{rms}$ and $u_{rms}$ are 
approximately parallel for given values of $\nu$ and $\eta$ implies that 
there exists a strong coupling between the rate of dissipation of the 
magnetic field and the rate of generation of hydrodynamic motions and,
hence, the magnetic field is not complitely force free. In the absense of 
such coupling the field should decay on a time scale characterized by the 
magnetic diffusivity whereas dissipation of motions is determined by the 
viscous time scale. The power law dependence of $b_{rms}$ and $u_{rms}$ 
means that the magnetic and kinetic energies are coupled by
\begin{equation}
\frac{\dot{E}_{\rm M}}{E_{\rm M}} \approx \frac{\dot{E}_{\rm K}}{E_{\rm K}}
\end{equation}  
for each calculated model. Note that the behaviour shown in Fig.1 seems to 
be rather general. To check this we have performed computations for 
time-dependent viscosity and magnetic diffusivity. These computations also
clearly indicate the existence of a self-similar regime with a coherent
decay of $b_{rms}$ and $u_{rms}$. Departures from the simple expression (4)
are significant only at very large $t$, and the reason of these departures 
is seen from the next figure. 

We plot in Fig. 2 the time dependence of the characteristic length scales 
of the magnetic field, $L_{B}$, and velocity, $L_{v}$, normalized to 
$L$. These length scales are determined by
\begin{equation}
L_{B}(t) = \sqrt{ 8 \pi E_{\rm M}(t) \left[ \int (\nabla \times {\bf B})^{2} dV
\right]^{-1}},
\end{equation} 
\begin{equation}
L_{v}(t) = \sqrt{2 E_{\rm K}(t) \left[ \int \rho (\nabla \times {\bf v})^{2} dV
\right]^{-1}}.
\end{equation}
Note that the calculated length scales are approximately equal except for the 
very early and very late phases. The length scales show the same
self-similar behaviour as the energies do. The self-similar regime is 
reached after approximately $0.1 L/ c_{A}$ and lasts $\sim 100-500 L/c_{A}$ 
until $L_{B}$ and $L_{v}$ approach to the length of the computational box. 
At this point the solution loses its self-similarity because the boundary 
conditions becomes important for the further evolution. In the 
self-similar regime, both length scales evolve $\propto t^{1/2}$, which
implies that the dissipation scales increase due to diffusion (either magnetic 
or viscous or both) rather than due to a nonlinear fluid interaction.  


\section{Simple analytic decay model}

The results of our computations can easily be interpreted in terms of
a simple analytic model. For the sake of simplicity, instead of solving
the general equations (1)-(3), consider incompressible fluid with
\begin{equation}
\nabla \cdot {\bf v} = 0
\end{equation} 
and with the magnetic field and velocity governed by Equations (2) and (3).
The results of our analysis should be in agreement with the
numerical modelling if the fluid velocity is less than the sound speed,
which is indeed always the case.
Again, assume that initially there exists the magnetic field in a fluid
which vanishes at infinity. The energy equation for the magnetic field
can be obtained from Equation (3) by multiplying it by ${\bf B}/4 \pi$
and integrating over the volume. After some transformations, we have
\begin{equation}
\frac{\partial E_{\rm M}}{\partial t} = - W_{\rm L} + \dot{\varepsilon}_{\eta} +
\frac{1}{4 \pi} \int \nabla \cdot [ ({\bf v} \times {\bf B}) \times
{\bf B}] dV - \frac{\eta}{4 \pi} \int \nabla \cdot [ (\nabla \times 
{\bf B}) \times {\bf B}] dV,
\end{equation}
where
\begin{equation}
W_{\rm L}= \frac{1}{4 \pi} \int {\bf v} \cdot [ (\nabla \times {\bf B}) \times
{\bf B}] dV
\end{equation}   
is the work done by the Lorentz force per unit time, and
\begin{equation}
\dot{\varepsilon}_{\eta} = \frac{\eta}{4 \pi} \int (\nabla \times
{\bf B})^{2} dV
\end{equation}
is the rate of Joule heating. The last two integrals can be transformed
into surface integrals which are vanishing if ${\bf B} = 0$ at infinity,
or for periodic boundary conditions that have been adopted in the present work.
Then, the magnetic energy obeys
\begin{equation}
\frac{\partial E_{\rm M}}{\partial t} = -W_{\rm L} - \dot{\varepsilon}_{\eta}.
\end{equation}

The equation for turbulent kinetic energy can be obtained in a similar
way by multiplying Equation (2) by ${\bf v}$ and integrating over the
volume. We have
\begin{equation}
\frac{\partial E_{\rm K}}{\partial t} = W_{\rm L} - \dot{\varepsilon}_{\nu} -
\int \nabla (p {\bf v}) dV  - \frac{1}{2} \int \nabla (v^{2} {\bf v}) dV -
\int \nabla \cdot [(\nabla \times {\bf v}) \times {\bf v}],
\end{equation}  
where
\begin{equation}
\dot{\varepsilon}_{\nu} = \rho \nu \int (\nabla \times {\bf v})^{2} dV
\end{equation}
is the rate of heat production due to viscosity. The last three integrals
can be transformed into surface integrals and vanish. Therefore, we obtain
\begin{equation}
\frac{\partial E_{\rm K}}{\partial t} = W_{\rm L} - \dot{\varepsilon}_{\nu}.
\end{equation}
The work done by the Lorentz force characterizes a departure of the magnetic
configuration from a force-free one because we have $W_{\rm L}=0$ for a force-free
field. 

As it follows from numerical modelling, the rates of decay of the magnetic
and kinetic energy in a self-similar regime are related by Equation (4).
This equation implies that the ratio $\Gamma = E_{\rm K}/E_{\rm M}$ remains constant
that follows from a simple qualitative argument (Biskamp 1993) because
\begin{equation}
\frac{\partial \Gamma}{\partial t}  = \Gamma \left( \frac{\dot{E}_{\rm K}}{E_{\rm K}}
\right)- \frac{\dot{E}_{\rm M}}{E_{\rm M}}.
\end{equation}
Then, solving Equations (11), (14) and (15), we find
\begin{equation}
W_{\rm L} = \frac{1}{E} ( E_{\rm M} \dot{\varepsilon}_{\nu} - 
E_{\rm K} \dot{\varepsilon}_{\eta}),
\end{equation}
\begin{equation}
\dot{E}_{\rm M} = - \frac{E_{\rm M}}{E} (\dot{\varepsilon}_{\nu} +
\dot{\varepsilon}_{\eta}),
\end{equation}
\begin{equation}
\dot{E}_{\rm K} = - \frac{E_{\rm K}}{E} (\dot{\varepsilon}_{\nu} +
\dot{\varepsilon}_{\eta}),
\end{equation}
where $E= E_{\rm M} + E_{\rm K}$. It is convenient to represent 
$\dot{\varepsilon}_{\nu}$ and $\dot{\varepsilon}_{\eta}$ as
\begin{equation}
\dot{\varepsilon}_{\nu} = \frac{2 \nu E_{\rm K}}{L_{\rm K}^{2}},
\dot{\varepsilon}_{\eta} = \frac{2 \eta E_{\rm M}}{L_{B}^{2}}.
\end{equation}
Then, we finally have
\begin{equation}
\frac{\dot{E}_{\rm M}}{E_{\rm M}} = \frac{\dot{E}_{\rm K}}{E_{\rm K}} =
- \frac{2}{1 + \Gamma} \left( \frac{\Gamma \nu}{L_{\rm K}^{2}} +
\frac{\eta}{L_{B}^{2}} \right).
\end{equation}
Both dissipative processes, magnetic diffusivity and viscosity, contribute
to the decay rate of turbulence. From this one may expect that the decay of
$E_{\rm M}$ would be faster if $\nu$ becomes large. Our numerical results
should however that this is not the case. In the following we explain why
increasing $\nu$ does not in fact increase the rate of magnetic energy
dissipation.

\section{additional}

Silk damping

continuous cross-over from mostly helical to less helical fields.
Could correlate energy decay rate with helicity.

Sketch of flow of energy


\section{Conclusion}

\section*{Acknowledgement}



%\newpage
\section*{References}

M. Arons, Max C. 1975, ApJ 196, L77 \\
D. Biskamp, Welter H. 1989, Phys. Fluids B 1, 1964 \\
D. Biskamp 1993, Nonlinear Magnetohydrodynamics, Cambridge Univ. Press, 
Cambridge \\
D. Biskamp, W.-C. M\"uller 1999, Phys. Rev. Lett. 83, 2195 \\
A. Brandenburg 2000, In ..., Ed. Ferriz-Mas, Gordon \& Breach
S. Galtier, Politano H., Pouquet A. 1997, Phys. Rev. Lett. 79, 2807 \\
C. F. Gammie, Ostriker E.C. 1996, ApJ 466, 814 \\
T. Hatory 1984, J. Phys. Soc. Japan 53, 2539 \\
R. Kinney, McWilliams J.C., Tajima T. 1995, Phys. Plasma 2, 3623 \\
M.-M. Mac Low, R. Klessen, A. Burkert, M. Smith 1998, Phys. Rev. Lett. 80, 
2754 \\
C. Thompson, R. Duncan 1993, ApJ 408, 194 \\

%\newpage
\section*{Figure captions}

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx

\end{document}
