text:
# TVD Scheme for the Numerical Simulation of the Axisymmetrical Selfgravitating MHD Flows

## 1 Introduction

Magnetic field plays an important role in the various problems of plasma physics and astrophysics. It is well known, for example, that the density of the magnetic energy in the interstellar clouds is compared with the density of the turbulent energy and can exceed the density of the thermal energy (see, e.g., ). There are a wide class of problems dealing with the dynamics of the selfgravitating MHD flows. Note among others the gravitational collapse of the interstellar clouds, protostellar clusters and galaxies; origin of the molecular bipolar outflows and jets from young stars and galaxy cores; dynamics of the accretion flows and circumstellar disks.

A variety of finite-difference schemes for the solution of the gasdynamical equations have been developed during last five decades. Corresponding 1D, 2D, and 3D numerical codes have been applied to the simulation of various gasdynamical flows. The direct extension of the these schemes to the MHD, however, is quite complex problem due to the specific properties of the MHD equations (anisotropy, vanishing divergence of the magnetic field, etc.).

Among the non-monotonic generalization of gasdynamical schemes for numerical simulation of the MHD flows the Lax-Vendroff method (and its modifications) is used (see, e.g., ). This method is subjected to the non-physical oscillations on strong discontinuities . Another prevailing method for the simulation of the multidimensional selfgravitating MHD flows is the MHD extension of the SPH method .

A finite-difference scheme should approximate correctly the physical oscillations of the flows containing shock waves, contact discontinuities and rarefaction waves. The nonlinear monotonic methods of high resolution were suggested for the solution of these problems. Historically the first monotonic method was FCT method . It was developed for of linear advection equation and for gasdynamical equations. At the present time FCT schemes are used for the simulations of the MHD flows as well .

There are two basic approaches that are applied for the construction of the monotonic finite-difference schemes for the numerical solution of the MHD problems. Nonlinear Riemann problem is solved in the original Godunov’s methods at every mesh interface and every time step by various iteration procedures (see, e.g., ). On the basis of nonlinear Riemann solver, Dai and Woodward  extended the PPM method  there. Powell et al.  constructed an eight-wave eigensystem for the approximate Riemann solver. Jiang and Wu  applied a high-order WENO interpolation scheme to MHD equations. Another approach is splitting the hyperbolicity matrix onto positive and negative parts and following solution of the linear Riemann problem .

We have suggested recently a simple TVD scheme for the numerical solution of the MHD equations that has high resolution in the smooth regions of the solution (see ). This scheme does not demand the solution of full eigenvalues problem (i.e. a calculation of all eigenvectors and eigenvalues of the hyperbolicity matrix), since it uses only the spectral radius. On the basis of developed TVD scheme we elaborated 2D numerical method for the solution of the selfgravitational MHD equations in cylindrical coordinates . This method can be used for the simulation of isothermal and adiabatic plasma flows. We have applied it to simulation of the collapse of magnetized rotating protostellar clouds.

In outline this paper proceeds as follows. In Section 2 we consider the mathematical problem statement for the simulation of selfgravitating MHD flows in the frame of axisymmetrical approximation. In Section 3 the TVD scheme for the numerical solution of the 1D MHD system is constructed. In Section 4 the developed 1D TVD scheme is extended to the 2D case and numerical algorithm for simulation of the 2D axisymmetrical selfgravitating MHD flows is described. Section 5 gives the description of 2D MHD numerical code and the results of few test problems. The results of the numerical simulations of the magnetized protostellar cloud collapse is presented in Section 6. In the conclusion (Section 7) we discuss briefly the main properties of the presented numerical code and general perspectives of its development in future.

## 2 Problem statement

### 2.1 Basic equations

The system of the gravitational MHD reads:

$$\frac{\rho }{t}+_k(\rho v^k)=0,$$
(1)
$$\frac{}{t}\left(\rho v^i\right)+_k\left(\rho v^iv^k+g^{ik}P\sigma ^{ik}\right)=\rho ^i\mathrm{\Phi },$$
(2)
$$\frac{B^i}{t}[,[𝐯,𝐁]]^i=0,$$
(3)
$`{\displaystyle \frac{}{t}}\left(\rho \epsilon +\rho {\displaystyle \frac{𝐯^2}{2}}+{\displaystyle \frac{𝐁^2}{8\pi }}\right)+_k\left\{\rho v^k\left(\epsilon +{\displaystyle \frac{P}{\rho }}+{\displaystyle \frac{𝐯^2}{2}}\right)+[𝐁,[𝐯,𝐁]]^k\right\}`$
$`=\rho v^k_k\mathrm{\Phi },`$ (4)
$$^2\mathrm{\Phi }=4\pi G\rho ,$$
(5)
where

$$\sigma ^{ik}=\frac{B^iB^k}{4\pi }g^{ik}\frac{𝐁^2}{8\pi }$$
(6)
is the Maxwellian tension tensor, $`g^{ik}`$ is the metric tensor. The other variables in this system (and further) have its usual physical sense. The system (16) is written in covariant form. Therefore it can be rewritten easy for the case of arbitrary curvilinear coordinates.

To close the system (16) we use the equation of state of a perfect gas $`P=(\gamma 1)\rho \epsilon `$, where $`\gamma `$ is the adiabatic index. The equation of energy (4) in the case of isothermal plasma should be excluded and the pressure should be determined from the relation $`P=c_T^2\rho `$, $`c_T^2=RT/\mu `$, where $`c_T`$ is the isothermal sound speed, $`\mu `$ is the molecular weight.

The collapse of a rotating magnetized protostellar cloud can be studied in the frame of axisymmetric approximation if the initial uniform magnetic field $`𝐁`$ is parallel to the vector of angular velocity $`𝛀`$. In this case we can use the cylindrical coordinates $`(r,\phi ,z)`$, axis $`z`$ being parallel to the vectors $`𝐁`$ and $`𝛀`$. The variables in (16) do not depend on the azimuthal coordinate $`\phi `$, so equations (14) can be written then in divergent form as follows:

$$\frac{𝐮}{t}+\frac{𝐅}{r}+\frac{𝐆}{z}=𝐑,$$
(7)
where the vectors of the conservative variables $`𝐮`$, fluxes $`𝐅`$, $`𝐆`$ and sources $`𝐑`$ are determined by expressions:

$`𝐮`$ $`=`$ $`\{r\rho ,r\rho v_r,r^2\rho v_\phi ,r\rho v_z,`$ (8)
$`B_r,B_\phi ,rB_z,r(\rho \epsilon +\rho {\displaystyle \frac{𝐯^2}{2}}+{\displaystyle \frac{𝐁^2}{8\pi }})\}^T,`$
$`𝐅`$ $`=`$ $`\{r\rho v_r,r(\rho v_r^2+P+{\displaystyle \frac{𝐁^2}{8\pi }}{\displaystyle \frac{B_r^2}{4\pi }}),r^2(\rho v_rv_\phi {\displaystyle \frac{B_rB_\phi }{4\pi }}),`$ (9)
$`r\left(\rho v_rv_z{\displaystyle \frac{B_rB_z}{4\pi }}\right),0,v_rB_\phi v_\phi B_z,r\left(v_rB_zv_zB_r\right),`$
$`r\rho v_r(\epsilon +{\displaystyle \frac{P}{\rho }}+{\displaystyle \frac{𝐯^2}{2}})+rv_r{\displaystyle \frac{𝐁^2}{4\pi }}r{\displaystyle \frac{B_r}{4\pi }}(𝐯𝐁)\}^T,`$
$`𝐆`$ $`=`$ $`\{r\rho v_z,r(\rho v_zv_r{\displaystyle \frac{B_zB_r}{4\pi }}),r^2(\rho v_\phi v_z{\displaystyle \frac{B_\phi B_z}{4\pi }}),`$ (10)
$`r\left(\rho v_z^2+P+{\displaystyle \frac{𝐁^2}{8\pi }}{\displaystyle \frac{B_z^2}{4\pi }}\right),v_zB_rv_rB_z,v_zB_\phi v_\phi B_z,0,`$
$`r\rho v_z(\epsilon +{\displaystyle \frac{P}{\rho }}+{\displaystyle \frac{𝐯^2}{2}})+rv_z{\displaystyle \frac{𝐁^2}{4\pi }}r{\displaystyle \frac{B_z}{4\pi }}(𝐯𝐁)\}^T,`$
$`𝐑`$ $`=`$ $`\{0,\rho v_\phi ^2+P+{\displaystyle \frac{𝐁^2}{8\pi }}{\displaystyle \frac{B_\phi ^2}{4\pi }}r\rho {\displaystyle \frac{\mathrm{\Phi }}{r}},0,r\rho {\displaystyle \frac{\mathrm{\Phi }}{z}},`$ (11)
$`0,0,0,r\rho v_r{\displaystyle \frac{\mathrm{\Phi }}{r}}r\rho v_z{\displaystyle \frac{\mathrm{\Phi }}{z}}\}^T.`$
Here index $`T`$ denotes the operation of the transposition.

### 2.2 Initial and boundary conditions

Let us consider the problem of initial and boundary condition on the example of the interstellar (protostellar) clouds collapse. Under the condition of axial and equatorial symmetries we can considered the 2D computational domain in the cylindrical coordinate system: $`D=(0rR,0zZ)`$ with the characteristic sizes $`R`$ and $`Z`$. We suppose that contracting cloud is located inside the computational domain and cloud boundary does not coincide with the grid boundary.

The conditions on the external boundaries of the computational domain are defined by cloud contraction from the given volume. In this case the normal components of the velocity $`𝐯_n`$ on the external boundary are equal to zero. For density $`\rho `$, internal energy $`\epsilon `$, magnetic field $`𝐁`$ and tangential components of the velocity $`𝐯_\tau `$ we can write boundary conditions as $`/𝐧=0`$, where $`𝐧`$ is the outer normal unit vector. The internal boundary conditions correspond to the axial and equatorial symmetries ($`v_r=v_\phi =0`$ on the axe and $`/z=0`$, $`v_z=B_\phi =B_r=0`$ on the equator).

Initial conditions for the collapse of the rotating magnetic cloud are determined by values of the mass $`M_0`$, temperature $`T_0`$ and parameters $`ϵ_t`$, $`ϵ_m`$, $`ϵ_\omega `$. These parameters are ratios of the internal, magnetic and kinetic energies of rotation to the absolute value of gravitational energy, respectively. The parameters $`ϵ_t`$, $`ϵ_m`$, $`ϵ_\omega `$ allow us to obtain the following similarity relations for the variables characterizing the initial cloud state:

$$R_0=8.911110^{16}ϵ_t\left(\frac{M_0}{M_{}}\right)\left(\frac{T_0}{10\mathrm{K}}\right)^1\mathrm{cm},$$
(12)
$$\rho _0=6.747510^{19}ϵ_t^3\left(\frac{M_0}{M_{}}\right)^2\left(\frac{T_0}{10\mathrm{K}}\right)^3\mathrm{g}\mathrm{cm}^3,$$
(13)
$$\mathrm{\Omega }_0=7.520410^{13}ϵ_t^{3/2}ϵ_\omega ^{1/2}\left(\frac{M_0}{M_{}}\right)^1\left(\frac{T_0}{10\mathrm{K}}\right)^{3/2}\mathrm{s}^1,$$
(14)
$$B_0=1.234210^4ϵ_t^2ϵ_m^{1/2}\left(\frac{M_0}{M_{}}\right)^1\left(\frac{T_0}{10\mathrm{K}}\right)^2\mathrm{G}.$$
(15)
These similarity relations allow to generalize the results of simulation with the specified values $`ϵ_t`$, $`ϵ_m`$, $`ϵ_\omega `$ on the contracting clouds with other values of mass and temperature.

The minimal initial values of the parameters $`ϵ_t`$, $`ϵ_m`$, $`ϵ_\omega `$ depend on the critical mass $`M_{}`$. We can obtain the value of critical mass from virial theorem as follows:

$$\frac{M_{}}{M_0}=ϵ_\omega +\frac{3}{2}(\gamma 1)ϵ_t+\sqrt{(ϵ_\omega +\frac{3}{2}(\gamma 1)ϵ_t)^2+ϵ_m}.$$

## 3 Finite-difference scheme

### 3.1 TVD schemes

Let us construct the finite-difference scheme for adiabatical and isothermal MHD flows using the total variation diminishing (TVD) principle. TVD principle for schemes approximating the linear hyperbolic system is an extension of monotonic principle for the schemes approximating the linear advection equation.

The general idea of this principle consists in the numerical implementation of the total variation diminishing for the finite-difference Riemann invariants . An efficient way of the monotonic schemes construction was proposed by Vyaznikov et al. . It consists in the following steps:

1. The selection of a monotonic ‘base’ scheme of the first order of approximation.
2. The transformation of the ‘base’ scheme to the high resolution scheme by addition of the antidiffusional terms.
3. The construction of the antidiffusional coefficients for monotonicity guarantee.

The constructed TVD scheme should have a high resolution in the smooth regions of the solution but the order of approximation can decrease in the regions of large gradients.

### 3.2 Advection equation

To investigate the general ideology of the TVD scheme construction and its basic properties let us consider the linear advection equation:

$$\frac{\rho }{t}+v\frac{\rho }{x}=0,v=\mathrm{const},$$
(16)
that is defined in the domain $`D=\{t>0`$, $`\mathrm{}<x<+\mathrm{}\}`$ with the initial values $`\rho (x,0)=\rho _0(x)`$. We introduce an uniform spatial grid $`\{x_i=ih`$, $`i=0,\pm 1,\pm 2,\mathrm{}\}`$ with the fixed stepsize $`h`$ and consider the following class of ‘base’ schemes for the equation (16):

$`{\displaystyle \frac{\rho _i^{n+1}\rho _i^n}{\tau }}`$ $`+`$ $`v{\displaystyle \frac{\rho _{i+1}^n\rho _{i1}^n}{2h}}{\displaystyle \frac{w}{2h}}\left(\rho _{i+1}^n2\rho _i^n+\rho _{i1}^n\right)`$ (17)
$``$ $`{\displaystyle \frac{vw}{4}}\left(\alpha _{i+1/2}^{}{\displaystyle \frac{\rho _{i+1}^n\rho _i^n}{h}}\alpha _{i1/2}^{}{\displaystyle \frac{\rho _i^n\rho _{i1}^n}{h}}\right)`$
$`+`$ $`{\displaystyle \frac{v+w}{4}}\left(\alpha _{i+1/2}^+{\displaystyle \frac{\rho _{i+1}^n\rho _i^n}{h}}\alpha _{i1/2}^+{\displaystyle \frac{\rho _i^n\rho _{i1}^n}{h}}\right)=0.`$
The parameter $`w`$ is a coefficient of numerical viscosity, while the coefficients $`\alpha _{i+1/2}0`$ determine the antidiffusional terms.

If all $`\alpha _{i+1/2}=0`$ the scheme (17) gives the ‘base’ scheme of the first order of approximation in space and time. The condition of monotonicity of this scheme is $`w|v|`$. In the case $`w=|v|`$ the ‘base’ scheme transforms into the classical donor cell scheme .

To obtain the conditions of monotonicity in the case of $`\alpha _{i+1/2}0`$ we have to use the principle of maximum , that for the scheme (17) gives the inequalities:

$`1{\displaystyle \frac{1}{2}}\alpha _{i+1/2}^{}+{\displaystyle \frac{1}{2}}{\displaystyle \frac{\alpha _{i1/2}^{}}{R_{i1/2}^{}}}0,`$
$`1{\displaystyle \frac{1}{2}}\alpha _{i1/2}^++{\displaystyle \frac{1}{2}}{\displaystyle \frac{\alpha _{i+1/2}^+}{R_{i+1/2}^+}}0,`$ (18)
$`{\displaystyle \frac{wv}{2h}}\tau \left(1{\displaystyle \frac{1}{2}}\alpha _{i+1/2}^{}+{\displaystyle \frac{1}{2}}{\displaystyle \frac{\alpha _{i1/2}^{}}{R_{i1/2}^{}}}\right)+{\displaystyle \frac{w+v}{2h}}\tau \left(1{\displaystyle \frac{1}{2}}\alpha _{i1/2}^++{\displaystyle \frac{1}{2}}{\displaystyle \frac{\alpha _{i+1/2}^+}{R_{i+1/2}^+}}\right)`$
$`1,`$
where

$$R_{i+1/2}^+=\frac{\rho _i^n\rho _{i1}^n}{\rho _{i+1}^n\rho _i^n},R_{i1/2}^{}=\frac{\rho _{i+1}^n\rho _i^n}{\rho _i^n\rho _{i1}^n}$$

are smoothness analyzers. Note that last inequality (18) determines also the stability condition for the scheme (17).

To monotonize the scheme (17) a number of approaches were suggested (see, e.g., ). Vyaznikov et al.  proposed to consider the coefficients $`\alpha `$ as a piecewise-linear functions of the analyzers $`R`$. These functions must satisfy to the following conditions: 1) $`\alpha =0`$, $`R0`$; 2) $`0\alpha 2`$, $`R0`$. Let us take $`\alpha (R)`$ in the form of two-parametric set of piecewise-linear functions proposed by Chakravarthy and Osher :

$$\alpha (R)=\frac{1+\psi }{2}\beta minmod(1/\beta ,R)+\frac{1\psi }{2}minmod(\beta ,R),$$

where

$$minmod(x,y)=\frac{1}{2}(sign(x)+sign(y))\mathrm{min}(|x|,|y|).$$

The function $`minmod(x,y)`$ equals to zero if $`x`$ and $`y`$ have the different signs and to argument with the minimal modulus otherwise. The function $`\alpha (R)`$ is linear in the region of the smooth solution: $`1/\beta R\beta `$,

$$\alpha (R)=\frac{1+\psi }{2}+\frac{1\psi }{2}R$$
(19)
and increases from value $`\alpha _1=((1+\psi )\beta +(1\psi ))/(2\beta )`$ to $`\alpha _2=((1+\psi )+(1\psi )\beta )/2`$.

The scheme (17) has the second or third order of approximation on space in the regions of smooth solution. This result can be obtained by of expanding of the scheme (17) to the Taylor series in the neighborhood of point $`M=(t=t^n,x=x_i)`$. The corresponding first difference approximation for the scheme (17) is:

$$\frac{\rho }{t}+v\frac{\rho }{x}+\frac{v}{4}\left(\frac{1}{3}\psi \right)\frac{^3\rho }{x^3}h^2+O(\tau +h^3),$$

Therefore at the value $`\psi =1/3`$ in the region of smooth solution the scheme (17) attains the third order of approximation in space. For the other values of $`\psi `$ the scheme (17) has the second order of approximation.

The parameters $`\alpha _1`$ and $`\alpha _2`$ determine the range of values of the function $`\alpha (R)`$ in the smoothness region (where $`R1`$). The optimal value of the parameter $`\beta `$ attains at $`\alpha _2=2`$ when $`\beta =\beta _{max}=(3\psi )/(1\psi )`$. In this case we can find that $`\alpha _1=2/(3\psi )`$, $`\alpha _2=2`$. The stability condition is determined by last inequality (18). It can be rewritten in the following form:

$`{\displaystyle \frac{w\tau }{h}}`$ $`+`$ $`{\displaystyle \frac{(wv)\tau }{4h}}\underset{R}{\mathrm{max}}\left\{{\displaystyle \frac{\alpha }{R}}\right\}{\displaystyle \frac{(wv)\tau }{4h}}\underset{R}{\mathrm{min}}\left\{\alpha \right\}`$
$`+`$ $`{\displaystyle \frac{(w+v)\tau }{4h}}\underset{R}{\mathrm{max}}\left\{{\displaystyle \frac{\alpha }{R}}\right\}{\displaystyle \frac{(w+v)\tau }{4h}}\underset{R}{\mathrm{min}}\left\{\alpha \right\}1.`$
From this inequality and from the definition of $`\alpha (R)`$ we can obtain the stability condition of Courant, Friedrichs and Lewy (CFL condition, ):

$$\frac{w\tau }{h}\frac{4}{5\psi +\beta (1+\psi )}.$$
(20)
In the case of optimal values of the parameters $`\psi =1/3`$ and $`\beta =4`$ we can rewrite (20) by the following way:

$$\frac{w\tau }{h}0.4.$$

The developed scheme can be rewritten more briefly in the conservative form:

$$\frac{\rho _i^{n+1}\rho _i^n}{\tau }+\frac{f_{i+1/2}f_{i1/2}}{h}=0,$$
(21)
where

$`f_{i+1/2}`$ $`=`$ $`f_{i+1/2}^0{\displaystyle \frac{1\psi }{4}}minmod(f_{i+3/2}^{},\beta f_{i+1/2}^{})`$ (22)
$``$ $`{\displaystyle \frac{1+\psi }{4}}minmod(f_{i+1/2}^{},\beta f_{i+3/2}^{})`$
$`+`$ $`{\displaystyle \frac{1+\psi }{4}}minmod(f_{i+1/2}^+,\beta f_{i1/2}^+)`$
$`+`$ $`{\displaystyle \frac{1\psi }{4}}minmod(f_{i1/2}^+,\beta f_{i+1/2}^+),`$
$$f_{i+1/2}^0=\frac{f_{i+1}+f_i}{2}\frac{w}{2}\left(\rho _{i+1}^n\rho _i^n\right),$$

$$f_{i+1/2}^{}=f_{i+1/2}^0f_i,f_{i+1/2}^+=f_{i+1}f_{i+1/2}^0$$

and fluxes $`f_i=v\rho _i^n`$.

### 3.3 System of hyperbolic equations

The developed methodology of difference scheme construction for the linear advection equation can be generalized onto the system of the linear hyperbolic equations:

$$\frac{𝐮}{t}+𝖠\frac{𝐮}{x}=0,$$
(23)
where $`𝐮`$ is the vector of the unknown variables and $`𝖠`$ is the matrix of the constant coefficients (hyperbolicity matrix). This system can be rewritten in the conservative form:

$$\frac{𝐮}{t}+\frac{𝐅}{x}=0$$

with the vector of fluxes: $`𝐅=𝖠𝐮`$.

The system of equations (23) can be represented in the equivalent form of the equations for the Riemann invariants

$$S^\alpha =(𝐥^\alpha 𝐮)$$
(24)
as

$$\frac{S^\alpha }{t}+\lambda _\alpha \frac{S^\alpha }{x}=0,$$
(25)
where $`𝐥^\alpha `$ are the left eigenvectors of the matrix $`𝖠`$ and $`\lambda _\alpha `$ are its eigenvalues. The original variables $`𝐮`$ can be obtained using the inverse transformation:

$$𝐮=\underset{\alpha =1}{\overset{N}{}}𝐫_\alpha S^\alpha ,$$
(26)
where $`𝐫_\alpha `$ are the right eigenvectors of the matrix $`𝖠`$.

We approximate each equation in the system (25) by the difference scheme (21):

$$\frac{S_i^{\alpha ,n+1}S_i^{\alpha ,n}}{\tau }+\frac{H_{i+1/2}^\alpha H_{i1/2}^\alpha }{h}=0,$$
(27)
where fluxes

$`H_{i+1/2}^\alpha `$ $`=`$ $`H_{i+1/2}^{\alpha ,0}{\displaystyle \frac{1\psi }{4}}minmod(H_{i+3/2}^\alpha ,\beta H_{i+1/2}^\alpha )`$
$``$ $`{\displaystyle \frac{1+\psi }{4}}minmod(H_{i+1/2}^\alpha ,\beta H_{i+3/2}^\alpha )`$
$`+`$ $`{\displaystyle \frac{1+\psi }{4}}minmod(H_{i+1/2}^{\alpha +},\beta H_{i1/2}^{\alpha +})`$
$`+`$ $`{\displaystyle \frac{1\psi }{4}}minmod(H_{i1/2}^{\alpha +},\beta H_{i+1/2}^{\alpha +}),`$
$$H_{i+1/2}^{\alpha ,0}=\frac{H_{i+1}^\alpha +H_i^\alpha }{2}\frac{w_\alpha }{2}\left(S_{i+1}^{\alpha ,n}S_i^{\alpha ,n}\right),H_i^\alpha =\lambda _\alpha S_i^{\alpha ,n},w_\alpha |\lambda _\alpha |,$$

$$H_{i+1/2}^\alpha =H_{i+1/2}^{\alpha ,0}H_i^\alpha ,H_{i+1/2}^{\alpha +}=H_{i+1}^\alpha H_{i+1/2}^{\alpha ,0}.$$

To construct the scheme for the variables $`𝐮`$ we have to perform the transformation (24) in (27). Finally the scheme reads:

$$\frac{𝐮_i^{n+1}𝐮_i^n}{\tau }+\frac{𝐅_{i+1/2}𝐅_{i1/2}}{h}=0,$$
(28)
where fluxes are:

$`𝐅_{i+1/2}`$ $`=`$ $`𝐅_{i+1/2}^0{\displaystyle \frac{1\psi }{4}}{\displaystyle \underset{\alpha }{}}minmod(𝐅_{i+3/2}^\alpha ,\beta 𝐅_{i+1/2}^\alpha )`$ (29)
$``$ $`{\displaystyle \frac{1+\psi }{4}}{\displaystyle \underset{\alpha }{}}minmod(𝐅_{i+1/2}^\alpha ,\beta 𝐅_{i+3/2}^\alpha )`$
$`+`$ $`{\displaystyle \frac{1+\psi }{4}}{\displaystyle \underset{\alpha }{}}minmod(𝐅_{i+1/2}^{\alpha +},\beta 𝐅_{i1/2}^{\alpha +})`$
$`+`$ $`{\displaystyle \frac{1\psi }{4}}{\displaystyle \underset{\alpha }{}}minmod(𝐅_{i1/2}^{\alpha +},\beta 𝐅_{i+1/2}^{\alpha +}),`$
$$𝐅_{i+1/2}^0=\frac{𝐅_{i+1}+𝐅_i}{2}\frac{1}{2}\mathrm{𝖱𝖶𝖫}\left(𝐮_{i+1}^n𝐮_i^n\right),$$

$$𝖶=diag(w_1,w_2,\mathrm{},w_N),$$

$$𝐅_{i+1/2}^{\alpha \pm }=\frac{1}{2}\left(\lambda _\alpha \pm w_\alpha \right)\left(𝐥_\alpha \left(𝐮_{i+1}^n𝐮_i^n\right)\right)𝐫^\alpha ,$$

where $`𝖫`$, $`𝖱`$ – matrixes of the left and of the right eigenvectors of the matrix $`𝖠`$ respectively.

The difference scheme (28,29) satisfies to the total variation diminishing (TVD) principle , since it is constructed on the basis of the monotonic schemes for the Riemann invariants. In the present paper we consider the Lax–Friedrichs schemes with the coefficients $`w_\alpha =w\underset{\alpha }{\mathrm{max}}\{|\lambda _\alpha |\}`$. For this case the fluxes in the base scheme reads :

$$𝐅_{i+1/2}^0=\frac{𝐅_{i+1}+𝐅_i}{2}\frac{w}{2}\left(𝐮_{i+1}^n𝐮_i^n\right).$$
(30)
During the calculation of fluxes (29) the numerical procedure of the $`minmod`$-evaluation should be called very often. To decrease the computation time we consider a more simple scheme in which the additional fluxes are summed up. Finally the scheme (28) is determined by the fluxes:

$`𝐅_{i+1/2}`$ $`=`$ $`𝐅_{i+1/2}^0{\displaystyle \frac{1\psi }{4}}minmod(𝐅_{i+3/2}^{},\beta 𝐅_{i+1/2}^{})`$ (31)
$``$ $`{\displaystyle \frac{1+\psi }{4}}minmod(𝐅_{i+1/2}^{},\beta 𝐅_{i+3/2}^{})`$
$`+`$ $`{\displaystyle \frac{1+\psi }{4}}minmod(𝐅_{i+1/2}^+,\beta 𝐅_{i1/2}^+)`$
$`+`$ $`{\displaystyle \frac{1\psi }{4}}minmod(𝐅_{i1/2}^+,\beta 𝐅_{i+1/2}^+),`$
where

$$𝐅_{i+1/2}^{}=𝐅_{i+1/2}^0𝐅_i,𝐅_{i+1/2}^+=𝐅_{i+1}𝐅_{i+1/2}^0$$

and the fluxes $`𝐅_{i+1/2}^0`$ are determined by the expression (30).

This methodology can be generalized easy onto the system of the nonlinear equations:

$$\frac{𝐮}{t}+\frac{𝐅(𝐮)}{x}=0.$$
(32)
The finite-difference scheme (28) for the system (32) is determined by the same fluxes (30,31), but the coefficients $`w`$ we should take in the form:

$$w=w_{i+1/2}=\varphi \underset{\alpha }{\mathrm{max}}\{|\lambda _{\alpha ,i}|,|\lambda _{\alpha ,i+1}|\},\varphi 1,$$
(33)
where $`\lambda _\alpha `$ are the eigenvalues of the hyperbolicity matrix $`𝖠=𝐅/𝐮`$.

### 3.4 MHD modification of the finite-difference scheme

We construct now the finite-difference TVD scheme for the isothermal or adiabatical 1D MHD in Cartesian coordinates. We assume that the magnetic field and the velocity have all three components: $`𝐁=\{B_x,B_y,B_z\}`$, $`𝐯=\{v_x,v_y,v_z\}`$.

The MHD equations in the conservative form,

$$\frac{𝐮}{t}+\frac{𝐅}{x}=0$$

are determined by the vector of the conservative variables:

$$𝐮=\{\rho ,\rho v_x,\rho v_y,\rho v_z,B_x,B_y,B_z,\rho \epsilon +\rho \frac{𝐯^2}{2}+\frac{𝐁^2}{8\pi }\}^T$$

and by the vector of the fluxes:

$`𝐅`$ $`=`$ $`\{\rho v_x,P+\rho v_x^2+{\displaystyle \frac{𝐁^2}{8\pi }}{\displaystyle \frac{B_x^2}{4\pi }},\rho v_xv_y{\displaystyle \frac{B_xB_y}{4\pi }},\rho v_xv_z{\displaystyle \frac{B_xB_z}{4\pi }},`$
$`0,v_xB_yv_yB_x,v_xB_zv_zB_x,`$
$`\rho v_x(\epsilon +{\displaystyle \frac{P}{\rho }}+{\displaystyle \frac{𝐯^2}{2}}+{\displaystyle \frac{𝐁^2}{4\pi }}){\displaystyle \frac{B_x}{4\pi }}(v_xB_x+v_yB_y+v_zB_z)\}^T.`$
For the adiabatic case the hyperbolicity matrix $`𝖠`$ is

$$\left(\begin{array}{cccccccc}0& 1& 0& 0& 0& 0& 0& 0\\ & & & & & & & \\ A_{21}& (3\gamma )v_x& (1\gamma )v_y& (1\gamma )v_z& \gamma \frac{B_x}{4\pi }& A_{26}& A_{27}& \gamma 1\\ & & & & & & & \\ v_xv_y& v_y& v_x& 0& \frac{B_y}{4\pi }& \frac{B_x}{4\pi }& 0& 0\\ & & & & & & & \\ v_xv_z& v_z& 0& v_x& \frac{B_z}{4\pi }& 0& \frac{B_x}{4\pi }& 0\\ & & & & & & & \\ 0& 0& 0& 0& 0& 0& 0& 0\\ & & & & & & & \\ A_{61}& \frac{B_y}{\rho }& \frac{B_x}{\rho }& 0& v_y& v_x& 0& 0\\ & & & & & & & \\ A_{71}& \frac{B_z}{\rho }& 0& \frac{B_x}{\rho }& v_z& 0& v_x& 0\\ & & & & & & & \\ A_{81}& A_{82}& A_{83}& A_{84}& A_{85}& A_{86}& A_{87}& \gamma v_x\end{array}\right)$$
(34)
where:

$$A_{21}=\frac{\gamma 1}{2}𝐯^2v_x^2,A_{26}=(2\gamma )\frac{B_y}{4\pi },A_{27}=(2\gamma )\frac{B_z}{4\pi },$$

$$A_{61}=\frac{v_yB_xv_xB_y}{\rho },A_{71}=\frac{v_zB_xv_xB_z}{\rho },$$

$$A_{81}=v_x\left(\frac{c^2}{\gamma 1}+\left(1\frac{\gamma }{2}\right)𝐯^2+𝐚^2\right)+a_x(𝐚𝐯),$$

$$A_{82}=\frac{c^2}{\gamma 1}+\frac{𝐯^2}{2}+𝐚^2a_x^2(\gamma 1)v_x^2,$$

$$A_{83}=(1\gamma )v_xv_ya_xa_y,A_{84}=(1\gamma )v_xv_za_xa_z,$$

$$A_{85}=\gamma v_x\frac{B_x}{4\pi }v_y\frac{B_y}{4\pi }v_z\frac{B_z}{4\pi },A_{86}=\gamma v_x\frac{B_y}{4\pi }+v_x\frac{B_y}{2\pi }v_y\frac{B_x}{2\pi },$$

$$A_{87}=\gamma v_x\frac{B_z}{4\pi }+v_x\frac{B_z}{2\pi }v_z\frac{B_x}{2\pi }.$$

We use here the sound speed, $`c=\sqrt{\gamma P/\rho }`$ and $`𝐚=𝐁/\sqrt{4\pi \rho }`$. The calculations give the following expressions for the eigenvalues of the hyperbolicity matrix (34):

$$\lambda _1=0,\lambda _2=v_x,\lambda _{3,4}=v_x\pm a_x,$$

$$\lambda _{5,6,7,8}=v_x\pm \sqrt{\frac{c^2+𝐚^2}{2}\pm \frac{1}{2}\sqrt{(c^2+𝐚^2)^24a_x^2c^2}}.$$
(35)
The zero value of $`\lambda _1`$ is the consequence of zero value of the flux $`𝐅`$ component corresponding to $`B_x`$. The value of $`\lambda _2`$ corresponds to the entropic wave, $`\lambda _{3,4}`$ correspond to the Alfvén waves, and $`\lambda _{5,6,7,8}`$ correspond to the fast and slow magnetosonic waves. In the isothermal case the entropic wave disappears and instead adiabatical sound speed it should be used the isothermal sound speed $`c_T`$.

We should construct now the calculation rule for the coefficients $`w`$ in (33). The values of these coefficients should be not less than the maximal modulus of all eigenvalues of the hyperbolicity matrix. It is clear that this condition can be satisfied by the following way:

$$w_{i+1/2}=\varphi \mathrm{max}\{|v_{x,i}|+\mathrm{\Psi }_{x,i},|v_{x,i+1}|+\mathrm{\Psi }_{x,i+1}\},\varphi 1,$$
(36)
where

$$\mathrm{\Psi }_x=\sqrt{\frac{c^2+𝐚^2}{2}+\frac{1}{2}\sqrt{(c^2+𝐚^2)^24a_x^2c^2}}.$$

Therefore we can keep in the relation (33) only eigenvalues corresponding to the fast MHD waves. This allow us to adapt the scheme (28,30,31,36) to the simulation of the 1D MHD flows. If the magnetic field is absent the fast magnetosonic speed $`\mathrm{\Psi }_x`$ equals the sound speed $`c`$, therefore the transition to the gasdynamics occurs correctly. With the help of the splitting technique we can construct the scheme for multidimensional (2D and 3D) MHD flows. One of this scheme for the axisymmetric case we will describe in the next section.

## 4 Numerical method

### 4.1 General scheme

The computational grid is determined by points $`(r_i,z_j)`$ indexed by integer indexes $`i=0,1,\mathrm{},I`$; $`j=0,1,\mathrm{},J`$. These points denote the centers of the grid cells. The coordinates of the cell boundaries are defined as the arithmetic mean of corresponding coordinates of the cell centers and have half-integer indexes. The computational variables are related to the centers of the cells and the numerical fluxes are related to the cell boundaries. Above we developed the TVD scheme in the frame of 1D MHD approximation. Now we generalize this approach onto the 2D MHD approximation. The ‘base’ difference scheme for equations (7) is the following one:

$$\frac{𝐮_{i,j}^{n+1}𝐮_{i,j}^n}{\mathrm{}t}+\frac{𝐅_{i+1/2,j}𝐅_{i1/2,j}}{\mathrm{}r}+\frac{𝐆_{i,j+1/2}𝐆_{i,j1/2}}{\mathrm{}z}=𝐑_{i,j}.$$

The main problem here is the definition of the calculation rule for the numerical fluxes across the cell boundaries.

Let us consider for simplicity only the radial direction because of the expressions of fluxes for the vertical direction can be derived by similar way. The radial fluxes $`𝐅_{i+1/2,j}`$ and $`𝐅_{i1/2,j}`$ should be calculated with the help of the formulae (30,31,33). Relatively determination of the coefficients $`w`$ we should note the following.

In the definition of the vectors $`𝐮`$, $`𝐅`$ and $`𝐆`$ in the form (810) the radial coordinate $`r`$ is present. Therefore to calculate eigenvalues of the hyperbolicity matrices we can close formally the system (7) by the following obvious equation:

$$\frac{r}{t}=0.$$

It can be shown that this extension of the original MHD system not changes the standard set of the eigenvalues.

In fact we consider the following system of MHD equations:

$$\frac{𝐰}{t}+𝖠\frac{𝐰}{r}+𝖡\frac{𝐰}{z}=𝐑,$$
(37)
with the vector $`𝐰`$ determined as

$$𝐰=\{r,\rho ,\rho v_r,\rho v_\phi ,\rho v_z,B_r,B_\phi ,B_z,\rho \left(\epsilon +\frac{𝐯^2}{2}\right)+\frac{𝐁^2}{8\pi }\}^T.$$

Obviously in Cartesian coordinates the eigenvalues of matrices $`𝖠`$ and $`𝖡`$ are determined by the same expressions (35). Let us consider another vector of the variables $`𝐮`$ that consists of $`r`$ and eight components of the vector (8). Making the transformation of the equation (37) to the new variables $`𝐮`$ we obtain:

$$\frac{𝐮}{t}+𝖠^{}\frac{𝐮}{r}+𝖡^{}\frac{𝐮}{z}=𝐑^{},$$

where

$$𝖠^{}=𝖳^1\mathrm{𝖠𝖳},𝖡^{}=𝖳^1\mathrm{𝖡𝖳},𝐑^{}=𝖳^1𝐑,$$
(38)
and matrix $`𝖳=𝐮/𝐰`$ is Jacobian of transformation. The determinant of this matrix $`det\widehat{𝖳}=1/r^7`$ is not equal to zero. Hence the formulae (38) describes the similarity transformations of matrices $`𝖠`$ and $`𝖡`$ that does not change their eigenvalues.

This discussion allows us to take the values of the coefficients $`w_{i+1/2}`$ in (33) in the following form:

$$w_{i+1/2}=\varphi \mathrm{max}\{|v_{r,i}|+\mathrm{\Psi }_{r,i},|v_{r,i+1}|+\mathrm{\Psi }_{r,i+1}\},\varphi 1,$$
(39)
where

$$\mathrm{\Psi }_r=\sqrt{\frac{c^2+𝐚^2+\sqrt{(c^2+𝐚^2)^24a_r^2c^2}}{2}}.$$
(40)
These coefficients are determined by eigenvalues that correspond to the fast magnetosonic wave. The expressions (39) can be applied both to the adiabatical and to the isothermal plasma flows. The energy equation should be excluded in last case and as $`c`$ in (40) should be used the isothermal sound speed $`c_T`$.

Since the fluxes related to the different space directions are constructed independently, the $`z`$-part of the numerical scheme can be developed by similar way.

### 4.2 Poisson solver

Poisson equation for the gravitational potential is solved by ADI-method (alternating directions implicit method) of Douglas and Rachford :

$$\frac{\overline{\mathrm{\Phi }}\mathrm{\Phi }^p}{\tau _p}=\theta \mathrm{\Lambda }_{rr}\overline{\mathrm{\Phi }}+(1\theta )\mathrm{\Lambda }_{zz}\mathrm{\Phi }^p\frac{\rho }{2}$$

$$\frac{\mathrm{\Phi }^{p+1}\overline{\mathrm{\Phi }}}{\tau _p}=(1\theta )\mathrm{\Lambda }_{rr}\overline{\mathrm{\Phi }}+\theta \mathrm{\Lambda }_{zz}\mathrm{\Phi }^{p+1}\frac{\rho }{2}$$

where $`\theta `$ is the parameter of scheme implicity, $`\mathrm{\Lambda }_{rr}`$, $`\mathrm{\Lambda }_{zz}`$ are second order finite-difference approximations of the corresponding parts of Laplacian operator in $`r`$– and $`z`$–directions, respectively, $`\tau _p`$ is the iteration parameter.

Boundary conditions for the gravitational potential on the inner boundaries are regular: $`\mathrm{\Phi }/𝐧=0`$. Boundary condition on the external boundary is $`\mathrm{\Phi }_{int}=\mathrm{\Phi }_{ext}`$ where $`\mathrm{\Phi }_{int}`$ is the computing potential, $`\mathrm{\Phi }_{ext}`$ is the external potential of the protostellar cloud. It can be obtained from Poisson integral transformation of equation (5) solution:

$$\mathrm{\Phi }_{ext}(r,z)=G\underset{l=0}{\overset{\mathrm{}}{}}\frac{P_l(x)}{(r^2+z^2)^{l+1}/2}\underset{V^{}}{}V^{}\rho (r^{},z^{})(r^2+z^2)^{l/2}P_l(x^{}),$$

where $`x=z(r^2+z^2)^{l/2}`$, and $`P_l(x)`$ is the Legendre polynomial of $`l`$th order. All terms with odd $`l`$ are absent in this representation as a result of the axial symmetry.

### 4.3 Implementation of the $`div𝐁=0`$ condition

Another difficulty of the numerical simulation of the multidimensional MHD flows is connected with the necessity of the numerical implementation of the vanishing divergence condition for the magnetic field

$$div𝐁=0.$$
(41)
Numerical errors due to discretization and finite-difference differential operator splitting can lead to nonzero divergence over the time. This problem can be solved by many ways . In the first approach the 8-wave formulation of the MHD equations is used (see also ). Another way uses the constrained transport suggested by Evans and Hawley . Finally, the third approach is the projection scheme . Our numerical scheme uses the last method (see ).

Denote as $`𝐁^{}`$ the magnetic field obtained by the numerical solution of the induction equation (3). This ‘numerical’ magnetic field contains the numerical error $`𝐛`$, that consist of vortex and potential parts:

$$𝐛=[,𝐚]+\varphi .$$

Therefore exact solution of induction equation is:

$$𝐁=𝐁^{}𝐛.$$
(42)
Vortex part of numerical error $`[,𝐚]`$ does not contributes to divergence of the magnetic field $`𝐁`$. Substitution of (42) into (41) shows that potential $`\varphi `$ must satisfy the Poisson equation

$$^2\varphi =div𝐁^{}.$$
(43)
Resolving this equation numerically for the defined ‘numerical’ magnetic field $`𝐁^{}`$ for each time step we calculate the divergence free magnetic field with the help of formula

$$𝐁^s=𝐁^{}\varphi .$$

Further we will assume that $`𝐁=𝐁^s`$. The algorithm is implemented in our numerical code as follows. Denote the central cell with indexes ($`i`$, $`j`$) by the index $`C`$ for brevity (see. Fig. 1), the neighbor cells by indexes $`W`$, $`N`$, $`E`$, $`S`$ and intermediate coordinates (cell boundaries) by double indexes $`CE`$, $`CS`$ etc.

Then we can approximate the divergence of magnetic field by the expression:

$$\left(div𝐁\right)_C=\frac{1}{r_C}\frac{r_{CE}B_{r,CE}r_{CW}B_{r,CW}}{r_{CE}r_{CW}}+\frac{B_{z,CN}B_{z,CS}}{z_{CN}z_{CS}},$$
(44)
where intermediate values are calculated by the formulae

$$B_{r,CE}=\frac{1}{2}\left(B_{r,C}+B_{r,E}\right),B_{r,CS}=\frac{1}{2}\left(B_{r,C}+B_{r,S}\right),\mathrm{etc}.$$

The finite-difference approximation of the Laplacian operator is

$`\left(\mathrm{}\varphi \right)_C`$ $`=`$ $`{\displaystyle \frac{1}{r_C(r_{CE}r_{CW})}}\left(r_{CE}{\displaystyle \frac{\varphi _E\varphi _C}{r_Er_C}}r_{CW}{\displaystyle \frac{\varphi _C\varphi _W}{r_Cr_W}}\right)`$
$`+`$ $`{\displaystyle \frac{1}{(z_{CN}r_{CS})}}\left({\displaystyle \frac{\varphi _N\varphi _C}{z_Nz_C}}{\displaystyle \frac{\varphi _C\varphi _S}{z_Cz_S}}\right).`$
It is easy to see that if the potential $`\varphi `$ is the exact solution of the finite-difference Poisson equation (43) then the magnetic field $`𝐁`$ (for which the finite-difference divergence (44) is equal to zero exactly) can be obtained by the formulae:

$$B_{r,C}=B_{r,C}^{}\frac{1}{2}\left[\frac{\varphi _E\varphi _C}{r_Er_C}+\frac{\varphi _C\varphi _W}{r_Cr_W}\right],$$

$$B_{z,C}=B_{z,C}^{}\frac{1}{2}\left[\frac{\varphi _N\varphi _C}{z_Nz_C}+\frac{\varphi _C\varphi _S}{z_Cz_S}\right].$$

The Poisson equation (43) for the potential $`\varphi `$ is solved in our scheme by the ADI-method that is described in the previous section. As the additional boundary conditions on the inner boundaries ($`r=0`$, $`z=0`$) we use the relations $`\varphi /𝐧=0`$. On the outer boundaries we use the boundary condition $`\varphi =0`$ since the external field has identically vanishing divergence.

## 5 Test computations

### 5.1 Background

The numerical code ‘Moon’ was developed on the basis of the proposed finite-difference scheme. It can simulate the 1D and 2D selfgravitating MHD flows. The code is implemented by the programming language C++ and completely object-oriented. Source text of the program uses only standard C++ elements and does not relate to concrete compiler. The kernel of the code uses intensively the object-oriented library ‘Numerical Tool Box’ (NTB 3.5). This library was developed by A.G.Zhilkin in 1996 for more easy creation of the complex numerical codes. To check the properties of our code we tested it on the number of test problems with known (exact or approximate) analytical solutions.

### 5.2 One-dimensional tests

#### 5.2.1 Test 1. Linear advection

The linear advection of the square density profile was simulated to check the phase-error properties of the numerical scheme. Figure 2 gives the density profile after 100th time step without (triangles) and with flux correction (dots). The computational grid consists of 100 cells. The analytical solution is drawn by the solid line.

It is seen that the 1st-order of approximation method has a large numerical diffusion. This diffusion smears the initial profile rapidly, the width of smearing being increased during calculation. The method of the 3rd-order of approximation produces a profile that is smeared onto 3–4 cells only and conserves its form during a long time. We can say therefore that amplitude and phase errors in the scheme of the 3rd-order of approximation are selfconsistent and its order of approximation is not less actually then the approximation order of the so-called LPE-schemes (little phase error) (see ). In our opinion this circumstance is connected with the successful selection of the antidiffusion limiters (22) in the scheme (21).

#### 5.2.2 Test 2. Decay of arbitrary gasdynamical discontinuity

Numerical simulation of the arbitrary discontinuity decay (Riemann problem) allows us to check the scheme properties relating to the resolution of the shock waves, contact discontinuities and rarefaction waves. At the initial time moment the computational domain along the $`x`$-axis was divided on two subdomains $`A`$ and $`B`$. In the initial state the pressure and density have the following values: in subdomain $`A`$$`\rho _A=1`$, $`P_A=1`$; in subdomain $`B`$$`\rho _B=0.125`$, $`P_B=0.1`$. The velocity in both domains $`v_A=v_B=0`$ and adiabatic index $`\gamma =5/3`$. Figure 3 shows the analytical solution (solid line) and the results of the numerical computation (dots) at the time moment $`t=0.164`$. As a result of the discontinuity decay the rarefaction wave 1 propagates to the left, contact discontinuity 2 and shock wave 3 propagates to the right. The Figure shows that the scheme resolves the shock wave sufficiently well, but gives the small non-physical oscillations on the contact discontinuity. To smooth it is necessary to use additional limiters of antidiffusion fluxes.

#### 5.2.3 Tests 3, 4. Decay of arbitrary MHD discontinuity

The numerical simulation of an arbitrary MHD discontinuity decay is the unique possibility to check the MHD properties of the difference scheme. The problem statement in this test is similar to the problem statement in test 2. The initial parameters for this test and its exact analytical solution were taken from paper of Ryu and Jones  (variant 1a): $`\gamma =5/3`$, $`\rho _A=1`$, $`\rho _B=1`$, $`P_A=20`$, $`P_B=1`$, $`v_{x,A}=10`$, $`v_{x,B}=10`$, $`B_{y,A}=5`$, $`B_{y,B}=5`$, $`B_x=5`$. The values of another variables are taken to be zero. This situation can arise, for example, at the collision of two gaseous masses moving forward each to other at the angle $`\pi /4`$ to the direction of the magnetic field.

The results of the numerical computation and the analytical solution of this problem at the time $`t=0.08`$ are shown on the Fig. 4. As a result of the discontinuity decay two fast MHD shock waves (1 and 5) arise and propagate in the opposite directions. The rotational (alfvenic) discontinuity 2, the contact discontinuity 3 and the slow MHD shock wave 4 are formed between this waves. The slow MHD shock wave propagates immediately after the fast MHD shock wave to the region of the gas with the larger pressure.

In another variant of this test we consider the Riemann problem with the initial state vector $`(\rho ,v_x,v_y,v_z,B_y,B_z,P)`$ (see , variant 4c): region $`A`$$`(0.65,0.667,0.257,0,0.55,0,0.5)`$, region $`B`$$`(1,0.4,0.94,0,0,0,0.75)`$ with $`B_x=0.75`$ and $`\gamma =5/3`$. The switch-of slow shock wave is formed in this case. The results of computations at the time $`t=0.15`$ are shown on Fig. 5. A fast MHD shock wave 1, switch-of slow MHD shock wave 2, contact discontinuity 3 and gasdynamical shock wave 4 are formed.

These tests show that the scheme approximation of fast and slow shock waves as well as the rotational (alfvenic) discontinuity is good enough. On the contact discontinuities and on the switch-off shock wave small oscillations occur.

#### 5.2.4 Test 5. Alfvén wave

The problem of propagation of the finite amplitude Alfvén wave (the simple alfvenic wave) is interesting because it has the analytical solution (see ). The comparison of the numerical and the analytical velocity profiles and the magnetic field allows us to make conclusions about behavior of the phase and the amplitude errors in the MHD case. We solve the following problem. At the initial time moment the wave is concentrated in the region $`x_Lxx_R`$, where

$$B_y=10\mathrm{sin}\left(\pi \frac{x_Rx}{x_Rx_L}\frac{\pi }{2}\right),B_z=10\mathrm{cos}\left(\pi \frac{x_Rx}{x_Rx_L}\frac{\pi }{2}\right),$$

$$v_y=\frac{B_y}{\sqrt{4\pi \rho }},v_z=\frac{B_z}{\sqrt{4\pi \rho }}.$$

The wave propagates in the gas with parameters $`\rho =1`$, $`P=1`$, $`v_x=0`$, $`B_x=1`$, $`\gamma =1.4`$ with velocity $`a=B_x/\sqrt{4\pi \rho }`$ without profile distortion. In wave region the tangential component of the magnetic field rotates with the conservation of its absolute value. Figure 6. shows the numerical and analytical distributions $`B_y`$ and $`B_z`$ at the time moment $`t=0.3`$.

The comparison of these curves with each other allows us to make the conclusion that even in the case of full MHD system the essential smearing of the velocity and the magnetic field profiles is absent during a sufficiently long time of the computation. Hence the phase errors balance the amplitude ones with the high order of accuracy.

### 5.3 Two-dimensional tests

#### 5.3.1 Test 6. Two-dimensional advection

To check the diffusion properties of the developed 2D code we obtain the numerical solution of the 2D advection equation (see similar computations of Munz )

$$\frac{\rho }{t}+v_x(y)\frac{\rho }{x}+v_y(x)\frac{\rho }{y}=0$$

in box domain $`(0x,y1)`$. The analytical solution of this equation for the velocity profile

$$v_x=(yy_0)\omega ,v_y=(xx_0)\omega .$$

describes the rotation of the density profile around the point ($`x_0,y_0`$) with the angular velocity $`\omega `$. We use in this computation the following values: $`x_0=y_0=0.5`$, $`\omega =1`$. Hence to the time $`t=2\pi `$ the original density profile should carry out one full rotation and return to its initial position. The computation were done on the uniform grid with the cells number $`N=100\times 100`$. We compare the initial state with the profile that obtained after one full rotation.

Figure 7a shows the original density profile. It is split on two regions. In the first region the density is equal to zero and in the second region the density is equal to $`1`$. The boundary of the regions is circle with the cut-out rectangle. The center of the circle is not coincides with the center of the domain. Figure 7b shows the density profile at the time $`t=2\pi `$ after one full revolution. The profile (as in the 1D case) is smeared onto 3–4 cells. Therefore we can conclude that 2D variant of our scheme has the good selfconsistent amplitude-phase properties as well.

#### 5.3.2 Test 7. Expansion of blast wave

The test problem of blast wave expansion in the uniform medium is very crucial for 2D numerical code due to spherical symmetry of the wave. Simulation the blast wave expansion allows us to check how precisely our 2D axisymmetrical numerical code can resolve spherically-symmetrical problems. We compute, in particular, the expansion of the supernova remnant on the adiabatic stage. This stage is described by Sedov–Taylor solution , the approximation being true while the radiation cooling is small.

We use the following initial parameters. The supernovae explosion energy is equal to $`410^{50}\mathrm{erg}`$, the density of the interstellar medium $`\rho _0210^{24}\mathrm{g}\mathrm{cm}^3`$, the spatial scale of the computational domain $`R_04.0310^{19}\mathrm{cm}`$, the pressure in the central region differs from the pressure in the external medium at the initial time moment approximately by eight order of magnitude. The scales of the computational values are equal to: for distance — $`13\mathrm{pc}`$, for velocity — $`8.9\times 10^8\mathrm{cm}\mathrm{s}^1`$, for pressure — $`1.1\times 10^6\mathrm{erg}\mathrm{cm}^3`$, for time — $`1510\mathrm{Yrs}`$. The expansion of the supernova remnant can be followed numerically over the time $`t=10600\mathrm{Yrs}`$. The radius of the shock wave increases from $`3.6\mathrm{pc}`$ to $`12\mathrm{pc}`$ during this time.

Figure 8 shows the dependencies of the radius $`R`$, velocity $`U`$ and the pressure $`P`$ of the shock wave on the time. All curves correspond to power law: $`(R,U,P)t^k`$. We obtained the following exponents: for radius of the shock wave $`k_R=0.403`$, for velocity $`k_U=0.589`$, for pressure $`k_P=1.139`$. The analytical values are $`k_R=0.4`$, $`k_U=0.6`$ and $`k_P=1.2`$. We see that the numerical values of the exponents are very close to exact Sedov–Taylor’s values. The condition of the spherical symmetry over the computations satisfies within the accuracy of 0.3%.

#### 5.3.3 Test 8. Spherically-symmetrical free-fall collapse

For simulation of free-fall (pressure-free) collapse as an initial state we take the uniform spherically-symmetrical cloud with the mass $`M=1M_{}`$ and with the initial density $`\rho _0=10^{19}\mathrm{g}\mathrm{cm}^3`$. Figure 9 shows the evolution of the numerical density (dots) of the cloud over the time in comparison with the analytical solution (solid line). We can see that the numerical points of density coincide with the high accuracy with values of the analytical solution. The cloud contracts to the point at the free-fall time $`t_{ff}=\sqrt{3\pi /(32G\rho _0)}`$. The computation is carried out to time moment when the cloud size decreases to the size of one grid cell. The difference of numerical and analytical values of density is not greater than 1% at this time and the law of the spherical symmetry conservation satisfies within the accuracy of 0.4%.

#### 5.3.4 Test 9. The contraction of polytropic cloud with $`\gamma =4/3`$

The simulation of the selfgravitating polytropic cloud dynamics with the adiabatic index $`\gamma =4/3`$ allows us to check the quality and the accuracy of the code in the simulation of the selfgravitating MHD flows. The hydrostatic equilibrium of polytropic cloud with $`\gamma =4/3`$ is indifferent one since the total energy, its first and second variations are equal to zero (see ). The density distribution of this cloud can be presented by analytical form as $`\rho (r)=\rho _c\chi ^3(r/R)`$, where $`\rho _c`$ is the central density, $`R`$ is the cloud radius, $`\chi (x)=\theta _3(x/\xi _1)`$ ($`\theta _3(\xi )`$ is the Emden function of index 3 and $`\xi _16.9`$ is its first positive zero).

To check the accuracy of the conservation of mass, momentum, and energy by our scheme we solve numerically the equations of the gravitational gasdynamics with the initial conditions corresponding to the specified hydrostatic state of the cloud. We computed several models with different numbers of cells $`N`$ of the spatial grid. Figure 10 shows the density profiles in the case of $`N=100\times 100`$. From the comparison of numerical (dots) and analytical (solid line) density profiles on this Figure we conclude that the relative error of the numerical solution compared with the analytical one at the time $`t=1.5t_{ff}`$ is less than 2%. The law of total mass conservation over the computations satisfies within the accuracy of 0.003% and the law of the total energy conservation satisfies within the accuracy of 2%.

This problem was solved on the grids with the number of cells $`N=20\times 20`$, $`80\times 80`$ and $`120\times 120`$ as well. We can say that the relative error of the numerical solution obtained by the time $`t=1.5t_{ff}`$ decreases with the increasing of the cells number. In the case of $`N=20\times 20`$ it equals to 20% while in the case of $`N=80\times 80`$ — 7%, and for $`N=120\times 120`$ it is less than one percent. We can conclude therefore that the convergence of the numerical solution to the analytical one with the rise of the grid resolution takes place.

In another variant of this test we specified the deviation from the hydrostatic equilibrium at the initial time moment by the velocity distribution $`v(r,t)=Ur/R`$. The cloud radius in this case evolves in time in accordance with the law $`R(t)=R_0Ut`$ and the evolution of the density and the velocity in the cloud can be described by the following analytical dependencies:

$$\rho (r,t)=\rho _c(0)\left(\frac{R_0}{R(t)}\chi \left(\frac{r}{R(t)}\right)\right)^3,u(r,t)=U\frac{r}{R(t)}.$$

To the time $`t_0=R_0/U`$ the cloud will contract to the point.

Figure 11 shows the numerical distributions of the density (circles) and the velocity (triangles) for different times. The solid lines correspond to analytical curves. From analysis of this Figure we can conclude that the numerical solution with the high accuracy coincides with the analytical one while the grid resolution is sufficient for the adequate spatial approximation of the central region. For the variant with $`N=100\times 100`$ this corresponds approximately to the time moment $`t=0.8t_0`$ (see Fig. 11).

## 6 Collapse of magnetized protostellar clouds

### 6.1 Introduction

Protostellar clouds as the cores of interstellar molecular clouds consist of neutral (atoms, molecules and cosmic dust particles) and charged (electrons, ions and charged dust) components. The equations of similar multicomponent mixture are very complex for the numerical simulations.

The more rough ideal MHD approximation can be used as a first approach. Contracting cloud is transparent to infrared self-radiation on the early stages of collapse. Therefore the initial phases of protostellar collapse can be considered as the isothermal ones with the good approximation. This approximation works while the optical depth of the collapsing cloud core becomes comparable with the unity.

We performed the computation of the isothermal magnetized cloud collapse using the uniform grid with the cells number $`N=120\times 120`$ and Courant number $`C=0.3`$. We use the following values of parameters. The initial ratio of the thermal and kinetic rotational energies to the absolute value of the gravitational energy $`ϵ_t=0.386`$, $`ϵ_\omega =0.0`$, respectively. The mass of the protostellar cloud $`M_0=1.5M_{}`$ and its temperature $`T_0=10\mathrm{K}`$. One the basis of the relations (1215) we can conclude that these parameters correspond to the cloud with the initial density $`\rho _0=5.1810^{18}\mathrm{g}\mathrm{cm}^3`$ and radius $`R_0=5.1710^{16}\mathrm{cm}`$. As a result of computations we calculate the flatness degree (ratio of the cloud thikness to its radius) and the central density up to the free-fall time.

### 6.2 Kinematics

The simplest case of the collapse corresponds to infinitesimally weak magnetic field. If the initial value $`ϵ_m0.510^4`$ then we can neglect the influence of the magnetic field on the collapse dynamics. The magnetic field can be considered as a passive admixture and for the given velocity field and may be determined from the equation of induction (3).

The results of computations show that in this kinematic approach the magnetic field $`B\rho ^k`$ with $`k=2/3`$ since of the cloud contraction is practically spherically-symmetrical one. Note that in the cloud envelope this law fulfills approximately and in the forming core it satisfies exactly . The computations confirm the early obtained conclusion that the magnetic field in the cloud acquires with the time the quasi-radial geometry (see Fig. 12) . This geometry of the field becomes distinct already to the free-fall time.

### 6.3 Dynamics

In the case of intermediate magnetic fields ($`0.1ϵ_m0.5`$) the collapse picture drastically differs from the kinematic statement. The gas particles moving across the magnetic field lines are deviated to the equatorial plane by the electromagnetic forces. As a result the flat (disk-like) structure is formed. Figure 13 shows the typical density profiles and geometry of the magnetic field lines at the moment of time equal to the magnetic free-fall time

$$t_{fm}=\sqrt{\frac{3\pi }{32G(1ϵ_m)\rho _0}}$$

(see ). The magnetic field in the disk is almost quasi-uniform and the ratio $`B_r/B_z`$ is small.

Numerical simulations show that at the free-fall time the flatness degree of the disk is proportional to $`ϵ_m^{1/2}`$ and the central density is inversely proportional to $`ϵ_m`$.

### 6.4 Quasi-steady-state contraction

The growth of the initial magnetic field leads to the increase of the electromagnetic forces that support the cloud from the contraction. Therefore at the large initial values $`ϵ_m`$ the cloud relaxes to the quasi-steady-state equilibrium. Such a cloud should evolve within the diffusion time scale that is determined by the ambipolar diffusion of the magnetic field. The cloud loses the magnetic flux by the ambipolar diffusion and slowly contracts retaining the quasi-hydrostatic equilibrium. As a result this evolution the non-uniform profile of the density is formed. The lifetime of such molecular cloud cores can reach $`10^510^7`$ years.

## 7 Conclusion

The explicit conservative finite-difference TVD scheme of high resolution for the solution of various MHD problems is constructed. The numerical MHD code ‘Moon’ on the basis of this scheme is developed that can simulate the MHD flows in the 1D and 2D approximations. The results of extensive tests show that the code is well adapted to the simulation of many MHD problems of plasma physics and astrophysics.

To check the dissipative properties of the developed scheme the problems of 1D and 2D advection are solved. The results of these tests show that the scheme smears the initial density profile onto 3–4 cells. In the test computation of the blast wave the condition of the spherical symmetry satisfies within the accuracy of 0.3%.

To analyze the convergence of the numerical solution to the analytical one the problem of the hydrostatic equilibrium of the selfgravitating polytropic cloud with the adiabatic index $`\gamma =4/3`$ is computed on the grids with various number of cells. The obtained results allow us to conclude that the errors of the numerical solution decrease to the permissible values on the grids with the number of cells not less than $`N=100\times 100`$. In the case $`N120\times 120`$ the relative error of the numerical solution at the time moment $`t=1.55t_{ff}`$ is less that 1%. The numerical solution in this case coincides with the good accuracy with the analytical one while the grid resolution in the central region is sufficient for acceptable approximation of the grid values.

The numerical computations of the collapse of the magnetized protostellar cloud confirms the earlier obtained results in the framework of 1.5D approximation (see ). The magnetic field in the cloud with the initially small intensity acquires the quasi-radial geometry after some time. The intermediate magnetic field leads to flattening of the collapsing clouds on the late stages of contraction. For the strong initial magnetic field the cloud relaxes to quasi-hydrostatic equilibrium. Such a cloud must evolve within the diffusion time scale.

Note that we can not simulate numerically the advanced stages of the collapse because the described code ‘Moon’ uses the uniform grid only. But investigation of these stages is very important and interesting problem to understand the physics of the formation and further evolution of young stellar objects. We are elaborating now the adaptive mesh refinement (AMR) numerical algorithm for the simulation of advanced stages of rotating magnetized protostellar cloud collapse. The AMR-approach uses the tree-threaded collection of subgrids with the subsequently increasing grid resolution. We are constructing a new 2D MHD numerical code ‘Megalion’ that based on the described in this paper high-resolution TVD scheme for the gravitational compressible MHD flows. This code will allow us to simulate various selfgravitating MHD flows with the large gradients and fast time variations. New numerical code will take into account also the processes of ambipolar diffusion, non-stationary ionization and heating/cooling.

This work is supported partially by grants RFBR N 96–02–19005, N 99–02–16938, N 00–02–17253, and N 00–01–10392. We wish to thank A.V.Koldoba and G.V.Ustyugova for many informative discussions.