# A Multi-dimensional Code for Isothermal Magnetohydrodynamic Flows in AstrophysicsTo appear in the Astrophysical Journal

## 1 INTRODUCTION

Over the last two decades, conservative upwind differencing schemes have proven to be very efficient for solving adiabatic hydrodynamic and magnetohydrodynamic (MHD) equations. These methods generally depend on the calculated estimates of mass, momentum and energy fluxes as well as magnetic field flux across cell boundaries based on the so-called “Riemann” solutions from the basic conservation laws. Examples for hydrodynamics include extensions of Godunov’s scheme (Godunov 1959), such as the MUSCL scheme (Van Leer 1979) and the PPM scheme (Colella & Woodward 1984), as well as those based on approximate flow eigenstates such as the Roe’s scheme (Roe 1981), the TVD scheme (Harten 1983) and the ENO scheme (Harten et al. 1987). Works for MHD include Brio & Wu (1988), Zachary & Colella (1992), Zachary, Malagoli, & Colella (1994), Dai & Woodward (1994a,1994b), Ryu & Jones (1995) (RJ, hereafter), Ryu, Jones, & Frank (1995) (RJF, hereafter), Powell et al. (1995), and Roe & Balsara (1996). Brio & Wu applied the Roe’s approach to the MHD equations. Zachary and collaborators used the BCT scheme to estimate fluxes for the MHD conservation equations. Dai & Woodward adapted the PPM scheme to MHD. Ryu and collaborators extended the Harten’s TVD scheme to MHD. Powell and collaborators developed a Roe-type Riemann solver with an eight-wave structure for MHD. Roe & Balsara constructed one variety of linearized Riemann solutions for MHD. The upwind schemes generally share an ability to sharply and cleanly define fluid discontinuities, especially shocks, and exhibit a robustness that makes them broadly applicable.

The assumption of adiabatic flows holds in the limit where cooling is negligible or the cooling time scale is much larger than the dynamical time scale. However, in the other limit where the cooling time scale is much shorter than the dynamical time scale, the assumption of isothermal flows becomes physically more valid (e.g., Draine & McKee 1993 and references therein). Of course, if cooling time scale is comparable to dynamical time scale, cooling should be considered explicitly.

Usually, numerical simulations of isothermal flows are made with adiabatic codes by setting the adiabatic index, $`\gamma `$, close to unity. Truelove et al. (1998) showed that with $`\gamma `$ as close as to unity as $`1.001`$, their adiabatic hydrodynamic code can follow isothermal collapse without any significant deterioration of accuracy. We also observed that with $`\gamma =1.001`$ the adiabatic TVD MHD code (RJ) captures structures in isothermal magnetohydrodynamic (IMHD) shock tubes without noticeable error. Yet, it is desirable to build codes specifically for isothermal flows, since those codes are simpler and faster than adiabatic ones. That is because the energy conservation equation need not to be solved in isothermal codes. As a result, the entropy mode, which carries the contact discontinuity, need not to be considered. In the current paper we describe an IMHD code based on Harten’s TVD scheme. It is the same scheme that was used for the adiabatic MHD code in RJ and RJF. Balsara (1998b) developed an IMHD code also based on an upwind scheme, but his scheme is different from ours.

In §2, we give recipes for the development of one- and multi-dimensional IMHD codes. In §3, we present the results of tests that include one-dimensional and two-dimensional shock tube problems, the decay of an Alfvén wave, and the nonlinear evolution of the Parker instability under a uniform gravity. Conclusions follow in §4.

## 2 NUMERICAL SCHEME

### 2.1 The Equations for Isothermal Magnetohydrodynamics

MHDs describes the behavior of the combined system of a conducting fluid and magnetic fields in the limit that the displacement current and the separation between ions and electrons are neglected. So the MHD equations represent coupling of the equations of fluid dynamics with Maxwell’s equations of electrodynamics. By ignoring the effects of electrical resistivity, viscosity, and thermal conductivity, and imposing isothermality on the conducting fluid, we get the following IMHD equations:

$$\frac{\rho }{t}+(\rho v)=0,$$
(2-1)   
$$\frac{v}{t}+vv+\frac{1}{\rho }(a^2\rho )\frac{1}{\rho }(\times B)\times B=0,$$
(2-2)   
$$\frac{B}{t}\times (v\times B)=0,$$
(2-3)   
with an additional constraint

$$B=0,$$
(2-4)   
for the absence of magnetic monopoles. Here $`a`$ is an isothermal sound speed, and other notations have their usual meanings. We incorporate a factor of $`1/\sqrt{4\pi }`$ into the definition of $`B`$ so that the factor of $`4\pi `$ does not appear in Eq. (2-2).

In Cartesian coordinates, the above equations are written in a conservative form as

$$\frac{q}{t}+\frac{F_x}{x}+\frac{F_y}{y}+\frac{F_z}{z}=0,$$
(2-5)   
$$q=\left(\begin{array}{c}\rho \\ \rho v_x\\ \rho v_y\\ \rho v_z\\ B_x\\ B_y\\ B_z\end{array}\right),F_x=\left(\begin{array}{c}\rho v_x\\ \rho v_x^2+a^2\rho +(B_y^2+B_z^2B_x^2)/2\\ \rho v_xv_yB_xB_y\\ \rho v_xv_zB_xB_z\\ 0\\ v_xB_yv_yB_x\\ v_xB_zv_zB_x\end{array}\right),$$
(2-6)   
with $`F_y`$ and $`F_z`$ obtained by properly permuting indices. With the state vector $`q`$ and the flux functions $`F_x(q)`$, $`F_y(q)`$, and $`F_z(q)`$, the Jacobian matrices, $`A_x(q)=F_x/q`$, $`A_y(q)=F_y/q`$, and $`A_z(q)=F_z/q`$ are formed. A system is called hyperbolic if all the eigenvalues of the Jacobian matrices are real and distinct and the corresponding set of right eigenvectors is complete (Jeffrey & Taniuti 1964). The system of the ideal, adiabatic MHD equations is known as non-strictly hyperbolic, since some eigenvalues coincide at some points (Brio & Wu 1988; Roe & Balsara 1996). The eigen-structure of the IMHD equations, which is presented in the next subsection, is very similar to that of the adiabatic ones. It is easy to show that the IMHD equations also form a non-strictly hyperbolic system.

### 2.2 One-Dimensional Code

Our strategy for developing a one-dimensional IMHD code is based the TVD scheme (Harten 1983) which was devised to improve the first-order-accurate Roe’s upwind scheme (Roe 1981) into a second-order-accurate one. For it, we derive the eigenvalues and eigenvectors of the system of the IMHD equations, which are given below. With the eigenvalues and eigenvectors, it is straightforward to apply the construction procedure for the one-dimensional adiabatic MHD code (RJ) to an isothermal analogue. Even though the procedure is described in RJ, we here repeat it to make this paper self-contained. Special attention is given to the orthonormal eigenvectors of the system of the IMHD equations.

We consider, as an example, plane-symmetric, one-dimensional flows exhibiting variation along the $`x`$-direction. Then $`y`$\- and $`z`$-derivatives in Eq. (2-5) are zero, and we have the one-dimensional IMHD equations

$$\frac{q}{t}+\frac{F_x}{x}=0,$$
(2-7)   
where $`q`$ and $`F_x`$ are defined in Eq. (2-6). The fifth equation in the system of Eqs. (2-7) is $`\frac{}{t}B_x=0`$, and the constraint in Eq. (2-4) is $`\frac{}{x}B_x=0`$. These imply that initially $`B_x`$ should be spatially constant and be kept constant during the evolution of the flow. So we need not include the equation for $`B_x`$ in a one-dimensional code.

The Jacobian matrix $`F_x/q`$ of the system of Eqs. (2-7) is given by

$$A_x=\left(\begin{array}{cccccc}0& 1& 0& 0& 0& 0\\ a^2v_x^2& 2v_x& 0& 0& b_y\sqrt{\rho }& b_z\sqrt{\rho }\\ v_xv_y& v_y& v_x& 0& b_x\sqrt{\rho }& 0\\ v_xv_z& v_z& 0& v_x& 0& b_x\sqrt{\rho }\\ v_x\frac{b_y}{\sqrt{\rho }}+v_y\frac{b_x}{\sqrt{\rho }}& \frac{b_y}{\sqrt{\rho }}& \frac{b_x}{\sqrt{\rho }}& 0& v_x& 0\\ v_x\frac{b_z}{\sqrt{\rho }}+v_z\frac{b_x}{\sqrt{\rho }}& \frac{b_z}{\sqrt{\rho }}& 0& \frac{b_x}{\sqrt{\rho }}& 0& v_x\end{array}\right),$$
(2-8)   
where $`b_{x,y,z}=B_{x,y,z}/\sqrt{\rho }`$. The six eigenvalues in non-decreasing order are

$$a_1=v_xc_f,$$
(2-9)   
$$a_2=v_xc_a,$$
(2-10)   
$$a_3=v_xc_s,$$
(2-11)   
$$a_4=v_x+c_s,$$
(2-12)   
$$a_5=v_x+c_a,$$
(2-13)   
$$a_6=v_x+c_f,$$
(2-14)   
where $`c_f`$, $`c_a`$, $`c_s`$ are the fast, Alfvén, and slow characteristic speeds, respectively. There is no entropy mode for the IMHD equations. The quantities $`a_1,\mathrm{},a_6`$ represent the six speeds with which information is propagated locally by three MHD wave families. The three characteristic speeds are expressed as

$$c_f=\left\{\frac{1}{2}\left[a^2+b_x^2+b_y^2+b_z^2+\sqrt{(a^2+b_x^2+b_y^2+b_z^2)^24a^2b_x^2}\right]\right\}^{1/2},$$
(2-15)   
$$c_a=|b_x|,$$
(2-16)   
$$c_s=\left\{\frac{1}{2}\left[a^2+b_x^2+b_y^2+b_z^2\sqrt{(a^2+b_x^2+b_y^2+b_z^2)^24a^2b_x^2}\right]\right\}^{1/2}.$$
(2-17)   
The six right eigenvectors corresponding to the six eigenvalues are

$$R_{v_x\pm c_f}=\left(\begin{array}{c}1\\ v_x\pm c_f\\ v_y\frac{c_fb_xb_y}{c_f^2b_x^2}\\ v_z\frac{c_fb_xb_z}{c_f^2b_x^2}\\ \frac{c_f^2b_y}{(c_f^2b_x^2)\sqrt{\rho }}\\ \frac{c_f^2b_z}{(c_f^2b_x^2)\sqrt{\rho }}\end{array}\right),$$
(2-18)   
$$R_{v_x\pm c_a}=\left(\begin{array}{c}0\\ 0\\ b_z\mathrm{sign}(b_x)\\ \pm b_y\mathrm{sign}(b_x)\\ \frac{b_z}{\sqrt{\rho }}\\ \frac{b_y}{\sqrt{\rho }}\end{array}\right),$$
(2-19)   
$$R_{v_x\pm c_s}=\left(\begin{array}{c}1\\ v_x\pm c_s\\ v_y\frac{c_sb_xb_y}{c_s^2b_x^2}\\ v_z\frac{c_sb_xb_z}{c_s^2b_x^2}\\ \frac{c_s^2b_y}{(c_s^2b_x^2)\sqrt{\rho }}\\ \frac{c_s^2b_z}{(c_s^2b_x^2)\sqrt{\rho }}\end{array}\right).$$
(2-20)   
Near the point where either $`b_x=0`$ or $`b_y=b_z=0`$, the above right eigenvectors are not well defined, with some elements becoming singular. By re-normalizing the eigenvectors, the singularities can be removed. The renormalized right eigenvectors are

$$R_{v_x\pm c_f}=\left(\begin{array}{c}\alpha _f\\ \alpha _f(v_x\pm c_f)\\ \alpha _fv_y\alpha _s\beta _yb_x\\ \alpha _fv_z\alpha _s\beta _zb_x\\ \frac{\alpha _s\beta _yc_f}{\sqrt{\rho }}\\ \frac{\alpha _s\beta _zc_f}{\sqrt{\rho }}\end{array}\right),$$
(2-21)   
$$R_{v_x\pm c_a}=\left(\begin{array}{c}0\\ 0\\ \beta _z\mathrm{sign}(b_x)\\ \pm \beta _y\mathrm{sign}(b_x)\\ \frac{\beta _z}{\sqrt{\rho }}\\ \frac{\beta _y}{\sqrt{\rho }}\end{array}\right),$$
(2-22)   
$$R_{v_x\pm c_s}=\left(\begin{array}{c}\alpha _s\\ \alpha _s(v_x\pm c_s)\\ \alpha _sv_y\pm \alpha _f\beta _ya\mathrm{sign}(b_x)\\ \alpha _sv_z\pm \alpha _f\beta _za\mathrm{sign}(b_x)\\ \frac{\alpha _f\beta _ya^2}{c_f\sqrt{\rho }}\\ \frac{\alpha _f\beta _za^2}{c_f\sqrt{\rho }}\end{array}\right),$$
(2-23)   
where $`\alpha `$’s and $`\beta `$’s are defined by

$$\alpha _f=\frac{\sqrt{c_f^2b_x^2}}{\sqrt{c_f^2c_s^2}},$$
(2-24)   
$$\alpha _s=\frac{\sqrt{c_f^2a^2}}{\sqrt{c_f^2c_s^2}},$$
(2-25)   
$$\beta _y=\frac{b_y}{\sqrt{b_y^2+b_z^2}},$$
(2-26)   
$$\beta _z=\frac{b_z}{\sqrt{b_y^2+b_z^2}}.$$
(2-27)   
At the points where $`b_y=b_z=0`$, $`\beta `$’s are defined as a limiting value, i.e.,

$$\beta _y=\beta _z=\frac{1}{\sqrt{2}}.$$
(2-28)   
Similarly, at the point where $`b_y=b_z=0`$ and $`b_x^2=a^2`$, $`\alpha `$’s are defined as

$$\alpha _f=\alpha _s=1.$$
(2-29)   
The left eigenvectors, which are orthonormal to the right eigenvectors, $`L_lR_m=\delta _{lm}`$, are

$$L_{v_x\pm c_f}=(l_{v_x\pm c_f}^{(1)},l_{v_x\pm c_f}^{(2)},l_{v_x\pm c_f}^{(3)},l_{v_x\pm c_f}^{(4)},l_{v_x\pm c_f}^{(5)},l_{v_x\pm c_f}^{(6)}),$$
(2-30)   
$$l_{v_x\pm c_f}^{(1)}=\frac{1}{\theta _1}\alpha _fa^2\pm \frac{1}{\theta _2}\left[\alpha _fav_x+\alpha _sc_s(\beta _yv_y+\beta _zv_z)\mathrm{sign}(b_x)\right],$$
(2-31)   
$$l_{v_x\pm c_f}^{(2)}=\pm \frac{1}{\theta _2}\alpha _fa,$$
(2-32)   
$$l_{v_x\pm c_f}^{(3)}=\frac{1}{\theta _2}\alpha _s\beta _yc_s\mathrm{sign}(b_x),$$
(2-33)   
$$l_{v_x\pm c_f}^{(4)}=\frac{1}{\theta _2}\alpha _s\beta _zc_s\mathrm{sign}(b_x),$$
(2-34)   
$$l_{v_x\pm c_f}^{(5)}=\frac{1}{\theta _1}\alpha _s\beta _yc_f\sqrt{\rho },$$
(2-35)   
$$l_{v_x\pm c_f}^{(6)}=\frac{1}{\theta _1}\alpha _s\beta _zc_f\sqrt{\rho },$$
(2-36)   
$$L_{v_x\pm c_a}=(l_{v_x\pm c_a}^{(1)},l_{v_x\pm c_a}^{(2)},l_{v_x\pm c_a}^{(3)},l_{v_x\pm c_a}^{(4)},l_{v_x\pm c_a}^{(5)},l_{v_x\pm c_a}^{(6)}),$$
(2-37)   
$$l_{v_x\pm c_a}^{(1)}=\pm \frac{1}{2}(\beta _zv_y\beta _yv_z)\mathrm{sign}(b_x),$$
(2-38)   
$$l_{v_x\pm c_a}^{(2)}=0,$$
(2-39)   
$$l_{v_x\pm c_a}^{(3)}=\frac{1}{2}\beta _z\mathrm{sign}(b_x),$$
(2-40)   
$$l_{v_x\pm c_a}^{(4)}=\pm \frac{1}{2}\beta _y\mathrm{sign}(b_x),$$
(2-41)   
$$l_{v_x\pm c_a}^{(5)}=\frac{1}{2}\beta _z\sqrt{\rho },$$
(2-42)   
$$l_{v_x\pm c_a}^{(6)}=\frac{1}{2}\beta _y\sqrt{\rho },$$
(2-43)   
$$L_{v_x\pm c_s}=(l_{v_x\pm c_s}^{(1)},l_{v_x\pm c_s}^{(2)},l_{v_x\pm c_s}^{(3)},l_{v_x\pm c_s}^{(4)},l_{v_x\pm c_s}^{(5)},l_{v_x\pm c_s}^{(6)}),$$
(2-44)   
$$l_{v_x\pm c_s}^{(1)}=\frac{1}{\theta _1}\alpha _sc_f^2\frac{1}{\theta _2}\left[\alpha _sc_av_x+\alpha _fc_f(\beta _yv_y+\beta _zv_z)\mathrm{sign}(b_x)\right],$$
(2-45)   
$$l_{v_x\pm c_s}^{(2)}=\pm \frac{1}{\theta _2}\alpha _sc_a,$$
(2-46)   
$$l_{v_x\pm c_s}^{(3)}=\pm \frac{1}{\theta _2}\alpha _f\beta _yc_f\mathrm{sign}(b_x),$$
(2-47)   
$$l_{v_x\pm c_s}^{(4)}=\pm \frac{1}{\theta _2}\alpha _f\beta _zc_f\mathrm{sign}(b_x),$$
(2-48)   
$$l_{v_x\pm c_s}^{(5)}=\frac{1}{\theta _1}\alpha _f\beta _yc_f\sqrt{\rho },$$
(2-49)   
$$l_{v_x\pm c_s}^{(6)}=\frac{1}{\theta _1}\alpha _f\beta _zc_f\sqrt{\rho },$$
(2-50)   
where

$$\theta _1=2(\alpha _f^2a^2+\alpha _s^2c_f^2),$$
(2-51)   
$$\theta _2=2(\alpha _f^2c_fa+\alpha _s^2c_ac_s).$$
(2-52)   
Some elements in the normalized right and left eigenvectors are not continuous. In order to force them to be continuous,

$$\mathrm{sign}(b_T)=\{\begin{array}{ccccccc}1,& \text{if}& b_y>0& \text{or}& b_y=0& \text{and}& b_z>0\\ 1,& \text{if}& b_y<0& \text{or}& b_y=0& \text{and}& b_z<0,\end{array}$$
(2-53)   
is multiplied to $`R_{v_x\pm c_s}`$ and $`L_{v_x\pm c_s}`$ if $`a^2>c_a^2`$, and to $`R_{v_x\pm c_f}`$ and $`L_{v_x\pm c_f}`$ if $`a^2<c_a^2`$.

Note that our eigenvectors have a different form from those derived in Balsara (1998a). The exact form of the eigenvectors does not matter, once all the singular points are taken care of.

Here, we use the conventional indices. The superscript $`n`$ represents the time step. The subscript $`i`$ indicates quantities at the cell center, while $`i+\frac{1}{2}`$ marks those at the right-hand cell boundary. The subscript $`k`$ represents the characteristic fields, with the order that $`k=1`$ is for the field associated with eigenvalue $`v_xc_f`$, $`k=2`$ for the field with $`v_xc_a`$, $`k=3`$ for the field with $`v_xc_s`$, $`k=4`$ for the field with $`v_x+c_s`$, $`k=5`$ for the field with $`v_x+c_a`$, and finally $`k=6`$ for the field with $`v_x+c_f`$.

An important step in the Roe’s scheme (1981) is to determine a Roe matrix $`\overline{A}_{x,i+1/2}(q_i,q_{i+1})`$ at the cell boundary from the adjacent state vectors, which satisfies the Roe’s suggested properties. One of them is $`F_{x,i+1}F_{x,i}=\overline{A}_{x,i+1/2}(q_{i+1}q_i)`$. For the systems of the adiabatic and isothermal hydrodynamic equations, there exists a Roe matrix evaluated at the $`\sqrt{\rho }`$-weighted average state (Roe 1981; LeVeque 1997). For the system of the adiabatic MHD equations, there is, however, no simple form of the Roe matrix except for the case with an adiabatic index $`\gamma =2`$ (Brio & Wu 1988). We have failed to find a simple form of the Roe matrix for the system of the IMHD equations, too. So we use an arithmetic averaging for the flow quantities at the cell boundary in the IMHD code, which was shown to work well in the adiabatic MHD code (RJ):

$$\rho _{i+1/2}=\frac{\rho _i+\rho _{i+1}}{2},$$
(2-54)   
$$v_{x,i+1/2}=\frac{v_{x,i}+v_{x,i+1}}{2},$$
(2-55)   
$$v_{y,i+1/2}=\frac{v_{y,i}+v_{y,i+1}}{2},$$
(2-56)   
$$v_{z,i+1/2}=\frac{v_{z,i}+v_{z,i+1}}{2},$$
(2-57)   
$$B_{y,i+1/2}=\frac{B_{y,i}+B_{y,i+1}}{2},$$
(2-58)   
$$B_{z,i+1/2}=\frac{B_{z,i}+B_{z,i+1}}{2}.$$
(2-59)   
The state vector $`q`$ at the cell center is updated by calculating the modified fluxes $`\overline{f}_x`$ at the cell boundaries as follows:

$$L_xq_i^{n+1}=q_i^n\frac{\mathrm{\Delta }t^n}{\mathrm{\Delta }x}(\overline{f}_{x,i1/2}\overline{f}_{x,i+1/2}),$$
(2-60)   
$$\overline{f}_{x,i+1/2}=\frac{1}{2}[F_x(q_i^n)+F_x(q_{i+1}^n)]\frac{\mathrm{\Delta }x}{2\mathrm{\Delta }t^n}\underset{k=1}{\overset{6}{}}\beta _{k,i+1/2}R_{k,i+1/2}^n,$$
(2-61)   
$$\beta _{k,i+1/2}=Q_k\left(\frac{\mathrm{\Delta }t^n}{\mathrm{\Delta }x}a_{k,i+1/2}^n+\gamma _{k,i+1/2}\right)\alpha _{k,i+1/2}(g_{k,i}+g_{k,i+1}),$$
(2-62)   
$$\alpha _{k,i+1/2}=L_{k,i+1/2}^n(q_{i+1}^nq_i^n),$$
(2-63)   
$$\gamma _{k,i+1/2}=\{\begin{array}{ccc}\frac{g_{k,i+1}g_{k,i}}{\alpha _{k,i+1/2}}& \text{for}& \alpha _{k,i+1/2}0,\\ 0& \text{for}& \alpha _{k,i+1/2}=0,\end{array}$$
(2-64)   
$$g_{k,i}=\text{sign}(\stackrel{~}{g}_{k,i+1/2})\text{max}[0,\text{min}\{|\stackrel{~}{g}_{k,i+1/2}|,\stackrel{~}{g}_{k,i1/2}\text{sign}(\stackrel{~}{g}_{k,i+1/2})\}]$$
(2-65)   
$$\stackrel{~}{g}_{k,i+1/2}=\frac{1}{2}\left[Q_k\left(\frac{\mathrm{\Delta }t^n}{\mathrm{\Delta }x}a_{k,i+1/2}^n\right)\left(\frac{\mathrm{\Delta }t^n}{\mathrm{\Delta }x}a_{k,i+1/2}^n\right)^2\right]\alpha _{k,i+1/2},$$
(2-66)   
$$Q_k(\chi )=\{\begin{array}{ccc}\frac{\chi ^2}{4ϵ_k}+ϵ_k& \text{for}& |\chi |<2ϵ_k,\\ |\chi |& \text{for}& |\chi |2ϵ_k.\end{array}$$
(2-67)   
Since the use of contact steepener and rotational steepener produces spurious numerical oscillations in the adiabatic MHD code (RJ; RJF), we do not include the rotational steepener in the IMHD code. The time step size $`\mathrm{\Delta }t^n`$ is restricted by the usual Courant condition for stability, $`\mathrm{\Delta }t^n=C_{\mathrm{cour}}\mathrm{\Delta }x/\mathrm{Max}(|v_{x,i+1/2}^n|+c_{f,i+1/2}^n)`$ with $`C_{\mathrm{cour}}<1`$.

### 2.3 Multi-dimensional Code

We extend the one-dimensional IMHD code to more dimensions by using a Strang-type directional splitting (Strang 1968). Here we explain, as an example, the implementation of it in two-dimensional plane-parallel geometry. The two-dimensional IMHD equations written in the conservative form (Eq. \[2-5\]) can be split into

$`L_x(x\text{-sweep}):`$ $`q_t+{\displaystyle \frac{F_x}{x}}=0,`$ (2-68)   
$`L_y(y\text{-sweep}):`$ $`q_t+{\displaystyle \frac{F_y}{y}}=0.`$ (2-69)   
In a time step, we update the state vector $`q(x,y)`$ along the $`x`$-direction with $`y`$ fixed, followed along the $`y`$-direction with $`x`$ fixed,

$$q^{n+1}=L_yL_xq^n.$$
(2-70)   
In order to maintain a second-order accuracy in time, the order of directional sweeps is permuted in the next time step by $`L_xL_y`$. The time step size, $`\mathrm{\Delta }t`$, is calculated at the start of the one complete sequence of $`L_xL_y`$ $`L_yL_x`$ and fixed through the sequence.

In multi-dimensional simulations, numerical solutions may not satisfy $`B=0`$ due to discretization errors. Brackbill & Barnes (1980) pointed out that errors of non-zero $`B`$ appear as a force parallel to the field. Non-zero $`B`$ can be removed, for instances, either by incorporating an explicit divergence cleaning method as described in RJF or by implementing a scheme similar to the constrained transport scheme (Evans & Hawley 1988) which was described in details for the adiabatic MHD code in Ryu et al. (1998). Tests in the next section have been done using the explicit divergence cleaning method, and the next two paragraphs describe it briefly.

At the beginning of MHD simulations, $`B=0`$ is satisfied. The updated magnetic field $`B`$, which is not in general divergence-free, can be decomposed as into two parts,

$$B=\varphi +\times V,$$
(2-71)   
where $`\varphi `$ and $`V`$ are scalar and vector functions respectively. Then the corrected magnetic field defined as $`B^c=B+\varphi `$ becomes divergence-free. So the problem of the divergence-cleaning is reduced to find $`\varphi `$, which is described by the Poisson equation

$$^2\varphi =B.$$
(2-72)   
In two-dimensional Cartesian geometry, for instance, the following finite difference representations

$$B_{x,i,j}^c=B_{x,i,j}+\frac{\varphi _{i+1,j}\varphi _{i1,j}}{2\mathrm{\Delta }x},$$
(2-73)   
$$B_{y,i,j}^c=B_{y,i,j}+\frac{\varphi _{i,j+1}\varphi _{i,j1}}{2\mathrm{\Delta }y},$$
(2-74)   
together with

$$\frac{\varphi _{i+2,j}2\varphi _{i,j}+\varphi _{i2,j}}{(2\mathrm{\Delta }x)^2}+\frac{\varphi _{i,j+2}2\varphi _{i,j}+\varphi _{i,j2}}{(2\mathrm{\Delta }y)^2}=\left(\frac{B_{x,i+1,j}B_{x,i1,j}}{2\mathrm{\Delta }x}+\frac{B_{y,i,j+1}B_{y,i,j1}}{2\mathrm{\Delta }y}\right),$$
(2-75)   
ensures

$$\frac{B_{x,i+1,j}^cB_{x,i1,j}^c}{2\mathrm{\Delta }x}+\frac{B_{y,i,j+1}^cB_{x,i,j1}^c}{2\mathrm{\Delta }y}=0,$$
(2-76)   
within machine round-off error. Extensions to the three-dimension and/or to other geometry are straightforward.

Eq. (2-75) is solved with boundary conditions specific to problems. In the problems with periodic boundaries (as in the decay test of the Alfvén wave in §3.2), a fast Poisson solver based on the fast Fourier transform can be used. In the nonlinear simulation of the Parker instability in §3.3 with reflection boundaries along one direction, the computational domain is doubled to that direction and the resulting boundaries are enforced to be periodic. In the two-dimensional shock tube tests in §3.1.2, doubling the computational domain in both directions also makes the resulting boundaries periodic. Note that in Eq. (2-75) $`\varphi `$’s are coupled with those at every other cell in a column and row. So the (extended) computational domain is divided into four sub-domains, and $`\varphi `$’s are computed in those sub-domains separately.

## 3 TESTS

In this section we present the results of three tests. The first and the second are isothermal versions of MHD shock tubes and decay of an Alfvén wave, respectively (RJ; RJF). The shock tube test shows the ability of the IMHD code to handle all the three MHD wave family structures, while the decay test of an Alfvén wave measures numerical dissipation. The third test is the simulation of a real astrophysical situation, the nonlinear evolution of the Parker instability under a uniform gravity. In all the tests, we set the isothermal speed $`a=1`$.

### 3.1 Shock Tube Tests

Based on the work of RJ, we have devised four shock tube problems which include discontinuities and rarefaction waves of IMHDs. To confirm the validity of our numerical solutions we have compared them to the analytic solutions obtained with an exact, nonlinear MHD Riemann solver described in RJ. That Riemann solver iterates from an initial guess of the solution for the full set of MHD waves based on the given left and right states. Iteration continues until the solutions to the innermost wave zone reached from the two opposite directions agree within some specified limit (in practice, a relative error $`10^5`$). The four shock tube solutions that we have applied in the test described below are listed in Table 1. $`C_{\mathrm{cour}}=0.8`$ and $`ϵ_1=0.3`$ (for fast modes), $`ϵ_2=0.0`$ (for Alfvén modes), $`ϵ_3=0.3`$ (for slow modes) have been used in the shock tube test calculations.

#### 3.1.1 One-Dimensional Shock Tube Tests

The one-dimensional simulations of the shock tube problems have been done with 512 cells in a computational tube bounded by $`x=[0,1]`$. We plot in following figures the resulting $`\rho `$, $`B_y`$, $`B_z`$, $`v_x`$, $`v_y`$, and $`v_z`$ at each cell with open circles and the analytic solutions with lines.

Figure 1a shows the result of the first shock tube problem at $`t=0.1`$ with the initial condition of a left state $`(\rho =1,v_x=0,v_y=0,v_z=0,B_y=5/\sqrt{4\pi },B_z=0)`$, a right state $`(\rho =0.1,v_x=0,v_y=0,v_z=0,B_y=2/\sqrt{4\pi },B_z=0)`$, and $`B_x=3/\sqrt{4\pi }`$. It exhibits the capturing of a fast rarefaction wave, a slow rarefaction wave, a slow shock, and a fast shock whose structures are plotted in the figure from left to right. There is no contact discontinuity. The fast and slow shocks are resolved sharply within several cells. In order to see the capturing rotational discontinuities, we have set up the initial condition of the second shock tube problem as: a left state $`(\rho =1.08,v_x=1.2,v_y=0.01,v_z=0.5,B_y=3.6/\sqrt{4\pi },B_z=2/\sqrt{4\pi })`$, a right state $`(\rho =1,v_x=0,v_y=0,v_z=0,B_y=4/\sqrt{4\pi },B_z=2/\sqrt{4\pi })`$, and $`B_x=2/\sqrt{4\pi }`$. Figure 1b shows the result at $`t=0.2`$. There are two fast shocks propagating outmost, and two slow shocks interior to those. Two rotational discontinuities lie between the fast and slow shocks. Here the strong fast shocks are resolved sharply, but the weak slow shocks and rotational discontinuities spread over more cells. The third shock tube problem has been set up with the initial condition of a left state $`(\rho =0.12,v_x=24,v_y=0,v_z=0,B_y=3/\sqrt{4\pi },B_z=0)`$, a right state $`(\rho =0.3,v_x=15,v_y=0,v_z=0,B_y=0,B_z=3/\sqrt{4\pi })`$, and $`B_x=0`$. Two oppositely moving magnetosonic shocks and a tangential discontinuity at $`t=0.2`$ are shown in Figure 1c. The magnetosonic shocks are again resolved sharply, but the tangential discontinuity spreads over $`20`$ cells. In the fourth shock tube problem, the initial condition has been set up with a left state $`(\rho =1,v_x=1,v_y=0,v_z=0,B_y=1,B_z=0)`$, a right state $`(\rho =1,v_x=1,v_y=0,v_z=0,B_y=1,B_z=0)`$, and $`B_x=0`$. The result at $`t=0.16`$ in Figure 1d shows two oppositely-moving identical magnetosonic rarefactions.

#### 3.1.2 Two-Dimensional Shock Tube Tests

The two-dimensional simulations of the shock tube problems have been done with $`256\times 256`$ cells in a computational domain bounded by $`x=[0,1]`$ and $`y=[0,1]`$. Initially, the domain is divided into two parts by a diagonal line joining the two points (0,1) and (1,0). The left state of the initial conditions for the one-dimensional shock tube problems has been assigned to the lower left part and the right state to the upper right part. The generated structures, including discontinuities and rarefactions, propagate parallel to the other diagonal line joining the two points (0,0) and (1,1).

In Figures 2a and 2b, two-dimensional correspondences of Figures 1a and 1b are plotted. In the figures, the following subscripts are used: $``$ for parallel components of velocity and magnetic field along the diagonal line joining the two points (0,0) and (1,1), $``$ for perpendicular components but still in the computational plane, and $`z`$ for components out of the plane. Although the resolution of the two-dimensional simulations, $`256\times 256`$ cells, is lower than that of the one-dimensional ones, 512 cells, we see all the structures have been captured correctly.

### 3.2 Decay of an Alfvén Wave

RJF carried out a test of the decay of linear waves in order to estimate numerical dissipations (resistivity and viscosity) in their adiabatic MHD code. Following the same idea, the decay of a linear Alfvén wave has been calculated and numerical dissipation in our IMHD code has been estimated. The IMHD equations for viscous and resistive fluid can be written as

$$\frac{\rho }{t}+(\rho v)=0,$$
(3-1)   
$$\frac{v}{t}+vv+\frac{1}{\rho }(a^2\rho )\frac{1}{\rho }(\times B)\times B=\frac{1}{\rho }_k\sigma _{ik},$$
(3-2)   
$$\frac{B}{t}\times (v\times B)=\eta ^2B.$$
(3-3)   
In the momentum equation, the viscosity tensor $`\sigma _{ik}`$ is given by,

$$\sigma _{ik}=\mu (_kv_i+_iv_k\frac{2}{3}\delta _{ik}v)+\zeta \delta _{ik}v,$$
(3-4)   
where $`\mu `$ and $`\zeta `$ are the dynamic shear and bulk viscosity, and $`\eta `$ is the electrical resistivity. Under uniform density, $`\rho _0`$, and uniform magnetic field, $`B=B_0\widehat{x}`$, the complex angular frequency of Alfvén waves is predicted from the linear analysis to be

$$\omega =\frac{i}{2}\left(\frac{\mu }{\rho _0}+\eta \right)k^2\pm c_\mathrm{A}k\left[1\frac{1}{4c_\mathrm{A}^2}\left(\frac{\mu }{\rho _0}\eta \right)^2k^2\right]^{1/2},$$
(3-5)   
where $`c_\mathrm{A}=\sqrt{B_0^2/2\rho _0}`$ is the Alfvén speed along the wave propagation direction and $`k=(k_x^2+k_y^2)^{1/2}`$ is the total wavenumber. Note that the complex angular frequency of the isothermal Alfvén waves is exactly the same as that of the adiabatic ones (RJF). We define the decay rate as

$$\mathrm{\Gamma }_\mathrm{A}=\frac{1}{2}\left(\frac{\mu }{\rho _0}+\eta \right)k^2.$$
(3-6)   
For the decay test of a linear Alfvén wave with the IMHD code, we have set up an initial condition such that, $`\rho _0=1`$, $`\delta v_z=v_{\mathrm{amp}}c_\mathrm{A}\mathrm{sin}(k_xx+k_yy)`$, $`B=1\widehat{x}`$, and all other quantities are equal to zero. The calculations have been done in a square periodic box with size $`L=1`$ using from $`8\times 8`$ cells to $`128\times 128`$ cells by increasing twice the number of cells in each direction. We set $`k_x=k_y=2\pi /L`$. Numerical parameters used are $`ϵ_1=ϵ_2=ϵ_3=0`$ and $`C_{\mathrm{cour}}=0.9`$ (the result is not sensitive to these values). Figure 3a shows the decay of the Alfvén wave calculated with $`32\times 32`$ cells. By fitting the peak points of the decay pattern with respect to time, we have estimated decay rate. In Figure 3b the resulting normalized decay rates as well as Reynolds numbers (see RJF for definition) are shown. Numerical Reynolds numbers scale almost as $`Rn_{\mathrm{cell}}^2`$ indicating the code has a second-order accuracy. Compared to the adiabatic MHD code, the IMHD code has smaller (up to 50%) numerical dissipation. This is partly because the IMHD code has one less mode (entropy mode).

### 3.3 Parker Instability under a Uniform Gravity

Nonlinear development of the Parker instability under a point-mass dominated gravity was simulated by Matsumoto and his collaborators (Matsumoto et al. 1988; Matsumoto et al. 1990; Matsumoto & Shibata 1992). And recently, Basu et al. (1996, 1997) simulated the nonlinear evolution of the Parker instability under a uniform gravity. As the final test of our IMHD code, we have also followed the nonlinear evolution of the Parker instability under the uniform gravity. By comparing our results with those in Basu et al. (1996), the code’s ability to handle a practical problem of astrophysics can be proved.

Since the Parker system is assumed initially to be in an isothermal equilibrium, an IMHD code is a natural tool for simulations. In the IMHD equations the externally given gravity -$`g\widehat{z}`$ is treated as a source term and placed on the right hand side of Eq. (2-5) with the source vector defined by $`S=(0,0,0,gv_z,0,0,0)^T`$. The gravity has the $`z`$-component only, so the source term is evaluated only when the state vector is updated along the $`z`$-axis

$$\frac{q}{t}+\frac{F_z}{z}=S.$$
(3-7)   
Since we use the Strang-type directional splitting in order to reduce multi-dimensional problems to one-dimensional ones, we also use the same technique to split the hyperbolic system with a source term into two parts,

Part A: $`{\displaystyle \frac{q}{t}}+{\displaystyle \frac{F_z}{z}}=0,`$ (3-8)   
Part B: $`{\displaystyle \frac{q}{t}}=S.`$ (3-9)   
Part A is solved by the TVD algorithm, and Part B by a forward-time-difference. To minimize a ‘splitting error’, Part A and Part B are solved by a BAB sequence with time step size 0.5$`\mathrm{\Delta }t`$ for Part B and $`\mathrm{\Delta }t`$ for Part A. The step size $`\mathrm{\Delta }t`$ is determined from the Courant condition.

The Parker system composed of isothermal gas and magnetic field $`\widehat{y}B_0(z)`$ takes under the uniform gravity an equilibrium configuration given by

$$\frac{\rho _0(z)}{\rho _0(0)}=\frac{B_0^2(z)}{B_0^2(0)}=\mathrm{exp}(z/H),$$
(3-10)   
where the gas scale height is $`H(1+\alpha )a^2/g`$ and the initial ratio of the magnetic to the gas pressure $`\alpha (B_0^2/[2\rho _0a^2])`$ is assumed a constant. We have chosen $`\alpha =1.0`$ in the simulation.

The computation domain covers $`0y12H`$ and $`0z12H`$. According to the linear stability analysis $`12H`$ is the horizontal wavelength corresponding to the maximum growth rate (Parker 1966). Periodic condition has been used in the $`y`$-boundaries, while reflecting condition in the $`z`$-boundaries. The density scale height $`H`$, the isothermal sound speed $`a`$, the initial midplane density $`\rho _0(0)`$, and the initial midplane field strength $`B_0(0)`$ have been chosen as the units of length, velocity, density, and magnetic field, respectively.

To initiate the instability we have added random velocity perturbations to the equilibrium profile of Eq. (3-10). Standard deviation of the perturbation velocity is $`10^4a`$ for each of the velocity components. To check whether the system follows, in the initial stage, the prescription of the linear analysis, we plot the logarithmic values of the root-mean-square velocity against time. In Figure 4a the dotted and dashed lines are for the horizontal and vertical components of the velocity, respectively. In the same figure the solid line represents the linear growth of a rate 0.34, which is the maximum growth rate of the system. At the very early stage the system undergoes a transient phase of adjustment, and then quickly develops the Parker instability at the rate predicted by the linear analysis. The linear growth gets saturated near $`t40`$.

The whole development of the Parker instability may be divided into three phases: The linear phase lasts up to $`t`$40, from then on the system undergoes the nonlinear phase until $`t57`$, and finally it reaches the damping oscillatory phase. Iso-contours and grey maps for density (left panels) and magnetic field lines and velocity vectors (right panels), in Figure 4b, present the snap shots of the system at the end of the linear phase ($`t=40`$), at the end of the nonlinear phase ($`t=57`$), and finally at $`t=80`$ of the damping oscillatory phase.

In the linear phase the perturbations grow predominantly in the upper region. In the nonlinear phase the perturbations gradually move towards the midplane. Through the linear and nonlinear phases, our simulation renders features that closely agree with those of Basu et al. (1996). As more matter accumulates, already compressed gas in the valley gets over-compressed. The increased gas pressure bounces the valley matter back to the upper region, and at the same time, the built-up pressure at the valley gets somewhat eased off. This in turn brings the matter back to the valley. The system now enters the oscillatory phase of the Parker instability. As the field lines are pushed deep down to the valley by the weight of over-lying matter, the curvature of the lines becomes small to the degree that magnetic field lines undergo reconnection. Due to the reconnection the matter drops down off the reconnected line, thereby the matter is allowed to move across the magnetic field line. The field line is now relieved from the burden of supporting the gas against the external gravity, and floats upwards. On the other hand, the field line located just below the reconnected one has to support more weight than before. Consequently this line now gets reconnected. In this way a redistribution of matter with respect to the field lines continues to occur, until there is no more reconnection. The system finally settles in an equilibrium. Since the reconnection drives the system to violate the flux-freezing condition, the final configuration of the system is different from that of the Mouschovias equilibrium (1974).

## 4 CONCLUSIONS

We have developed one- and multi-dimensional numerical codes to solve the IMHD equations, which are isothermal analogues of the previous adiabatic codes (RJ; RJF). Both the isothermal and adiabatic codes are based on the same scheme, an explicit finite-difference scheme on an Eulerian grid called TVD, which is a second-order-accurate extension of the Roe-type upwind scheme. The shock tube tests have showed that both codes capture correctly all the structures in MHDs. From the decay test of a linear Alfvén wave, we have found that numerical dissipation of the IMHD code is somewhat smaller than that of the adiabatic MHD code.

The robustness of the adiabatic code has been demonstrated through the simulations of MHD flows such as the Kelvin-Helmholtz instability (Frank et al. 1996; Jones et al. 1997) and jets (Frank et al. 1998; Jones et al. 1998), and that of the isothermal code has been done through the simulation of the Parker instability under the uniform gravity in this paper. Furthermore, both codes are fast enough to simulate multi-dimensional, astrophysical MHD flows using modest computational resources. Both codes run at about 400 MFlops on a Cray C90 processor (RJF), and the isothermal code updates zones about twice as fast as the adiabatic code. Together with the adiabatic code, the isothermal code is a useful tool to study the nonlinear evolution of astrophysical MHD flows.

###### Acknowledgements.

JK was supported by the Ministry of Science and Technology through Korea Astronomy Observatory grant 97-5400-000. The work by DR was supported in part by KOSEF through the 1997 Korea-US Cooperative Science Program 975-0200-006-2. The work by TWJ was supported in part by the NSF through grants AST93-18959, INT95-11654 and AST96-16964, by NASA grant NAG5-5055 and by the University of Minnesota Supercomputing Institute. SSH was supported in part by the Ministry of Education, Basic Science Research Institute grant No. BSRI-97-5411.      

FIGURES