A&A 404, 533-543 (2003)
DOI: 10.1051/0004-6361:20030556

Resistive MHD simulations of the Parker instability in galactic disks

G. Kowal1 - M. Hanasz2 - K. Otmianowska-Mazur1

1 - Astronomical Observatory, Jagiellonian University, ul. Orla 171, 30-244 Kraków, Poland
2 - Torun Centre for Astronomy, Nicholas Copernicus University, 87-148 Piwnice/Torun, Poland

Received 30 September 2002 / Accepted 25 March 2003

Parker instability leads to the formation of tangential discontinuities in a magnetic field and subsequent magnetic reconnection due to a numerical and/or an explicit resistivity. In this paper we investigate the role of the uniform, localized and numerical resistivity on the diffusion of magnetic field lines during the growth phase of Parker instability modes. We propose a new method to quantify the diffusion of magnetic field lines which is attributed to the presence of resistivity in ideal and non-ideal MHD codes. The method relies (1) on integration of magnetic lines in between periodic boundaries, (2) on measurements of the dispersion of magnetic field lines with the left and the right periodic boundaries and (3) on a statistical analysis of shifts of a large set of magnetic lines. The proposed method makes it possible to detect topological evolution of magnetic field. We perform a series of resistive MHD simulations of the Parker instability in uniformly rotating galactic disks. We follow the topological evolution of the magnetic field evolving due to the Parker instability and relate it to the ratio of total to uniform magnetic field in galactic disks. We find that after the onset of the Parker instability, the magnetic field becomes first tangled and later on it evolves toward a uniform state due to the presence of resistivity. A similar effect of a varying contribution of a turbulent magnetic field is observed in arms and inter-arm regions of galaxies.

Key words: ISM: magnetic diffusion - galaxies: structure - magnetohydrodynamics (MHD)

1 Introduction

The issue of stability of magnetized astrophysical disks is of primary importance for their dynamical evolution. Parker (1966, 1967) discovered that galactic disks containing a large scale horizontal magnetic field and vertical disk gravity are unstable against the buoyancy of magnetic field and cosmic rays. When the contribution of the weightless magnetic field and cosmic rays to the total pressure is significant, small perturbations of the system grow in amplitude, form wavelike perturbations of magnetic field lines and promote sliding of the heavy gas to magnetic valleys. The typical horizontal wavelength of unstable modes is of the order of a few hundreds parsecs up to a few kiloparsecs. This kind of instability transports the magnetic field out of the disks. On the other hand Parker instability is supposed to contribute to the magnetic field amplification through the process of fast galactic dynamo introduced originally by Parker (1992) and developed subsequently in a series of papers by Hanasz & Lesch (1993, 1998, 2000), Hanasz et al. (2002, hereafter referred to as Paper I), Brandenburg & Schmitt (1998) and Moss et al. (1999).

Another relevant MHD instability is the magneto-rotational instability, referred to as Balbus-Hawley instability (1991, 1992), which relies on the angular momentum exchange, by means of perturbed magnetic field lines, between different orbits of differentially rotating magnetized gaseous disks. The development of this instability initiates and sustains turbulence and leads to the amplification of the small scale magnetic field up to a certain level of saturation. The common feature of both instabilities is the development of non-uniformities and contact discontinuities on the initially smooth magnetic fields.

According to Ampere's law, magnetic field discontinuities correspond to current sheets. Currents are dissipated according to Ohm's law and magnetic energy is converted into the thermal energy. As a consequence magnetic field is no longer frozen in the gas flows and the magnetic field topology can evolve due to a nonidealness of the medium. Resistive processes may influence the development of MHD instabilities in astrophysical conditions. The issue has been undertaken in both cases of Balbus-Hawley (Fleming et al. 2000) and Parker instabilities (Paper I).

A significant part of studies of MHD instabilities in astrophysics have been done with the aid of numerical simulations. The intrinsic property of numerical codes solving MHD equations is the presence of a code-dependent numerical diffusivity. Therefore, even if the ideal set of equations is actually being solved, the numerical diffusion of the magnetic field changes the magnetic field topology. In view of the above facts, the main goal of the present project is to analyze the influence of the different kinds of resistivity on the magnetic field topology in numerical simulations of the Parker instability.

We apply three types of a resistivity: the uniform one, the localized resistivity operating only above some critical value of the current density, and the numerical one. Our model study is based on the general setup similar to that applied in Paper I, i.e. we consider a simple equilibrium state in a uniform gravity, assuming that magnetic pressure is initially proportional to the gas pressure. The initial magnetic field is vertically stratified, and is parallel to the azimuthal direction. We superpose a perturbation of velocity and simulate the Parker instability in a local Cartesian box with coordinates (x,y,z)parallel locally to the directions of r, $\phi$ and z coordinates. We take into account the disk rotation introducing the Coriolis force. We apply the resistivity model following Paper I (and references therein).

In order to perform a quantitative estimation of the effect of resistivity in a Parker unstable system we introduce a new method of measuring the topological variations of magnetic field lines. The method described in Sect. 3 allows us to examine the effects of different prescriptions of the resistivity and to compare them to the effect of the numerical resistivity. The applied concept is similar to the description of magnetic diffusion with the aid of Liapunov exponents (Barge et al. 1984; Barghouty & Jokipii 1996; Zimbardo 1995). Our method is different from the one applied by Fleming et al. (2000), where the resistivity is measured through the determination of the saturation of the growth of magnetic energy as the result of the Balbus-Hawley instability. In the case of Parker instability such a method does not provide any information about the resistivity since without the action of differential rotation we do not observe any growth of the magnetic energy.

We are primarily interested in the role of the resistivity in uniformization of the magnetic field structure distorted by the onset of the Parker instability. Referring to the picture of the galaxy NGC 6946 presented by Beck & Hoernes (1996) we note that highly unpolarized synchrotron emission comes from the regions of optical arms and a polarized emission comes from inter-arm regions. The reason for the depolarization in arms is obviously a strong agitation of the interstellar medium (ISM) by stellar explosions and winds. These phenomena form a very intense source of perturbations which are supposed to contribute to the excitation of the Parker and other instabilities. Since density waves are interpreted as waves traveling across the galactic disk area (Binney & Tremaine 1987) optical arms were placed in the past in the present location of inter-arm regions. One can ask therefore which mechanism is responsible for the converting of the turbulent medium in the arms in the much more uniform medium in the inter-arm regions. Within the framework of our model, we can suggest that the magnetic reconnection can play an important role in the uniformization of galactic magnetic fields.

2 Method

2.1 Equations

We investigate the Parker instability in galactic disk in the presence of a rigid rotation and a fast magnetic reconnection assuming isothermal evolution of the system following Paper I.

We solve a set of MHD equations including the Coriolis force in the equation of motion and the resistive term in the induction equation:

 \begin{displaymath}\frac{\partial \rho}{\partial t} + \nabla \cdot \left( \rho \vec{v} \right) = 0
~ ,
\end{displaymath} (1)

 \begin{displaymath}\frac{\partial \vec{v}}{\partial t} + \left( \vec{v} \cdot \n...
...{B}}{4 \pi \rho} - 2 \vec{\Omega} \times
\vec{v} + \vec{g} ~ ,
\end{displaymath} (2)

 \begin{displaymath}\frac{\partial \vec{B}}{\partial t} = \nabla \times \left( \vec{v} \times
\vec{B} - \eta ~ \vec{j} \right) ~ ,
\end{displaymath} (3)

where $\vec{\Omega}$ is the angular rotation velocity, $\vec{j} = \nabla \times
\vec{B}$ is the current density, $\vec{g}=g ~ \hat{\vec{e}_{z}}$, $g={\rm const}$is the gravitational acceleration and the other symbols have their usual meaning. Our assumption of uniform gravity is not realistic for actual galaxies, nevertheless it is satisfactory for a theoretical investigations of the Parker instability.

We adopt an isothermal equation of state

\begin{displaymath}p = c_{\rm s}^2 \rho ~ ,
\end{displaymath} (4)

where $c_{\rm s}= {\rm const}$ is the isothermal sound speed.

The following Ansatz for an anomalous resistivity (Ugai 1992; Konz et al. 2000; Tanuma et al. 2001; Paper I) is applied:

 \begin{displaymath}\eta (j) = \eta_1 + \eta_2 \left( j^2 - j_{\rm crit}^2 \right)^{1/2} \Theta
\left( j^2 - j_{\rm crit}^2 \right) ~ ,
\end{displaymath} (5)

where $\eta _1$ and $\eta _2$ are constant coefficients, $j_{\rm crit}$ denotes the critical current density above which the resistivity switches on and $\Theta$ is the Heaviside step function. Such a prescription of the resistivity allows for a useful parametrization of the magnetic reconnection, as discussed in Paper I.

Following Paper I we introduce the Coriolis force corresponding to the rotation angular velocity $\Omega$ which takes a typical galactic value defined in Sect. 2.3. The Coriolis force is responsible for horizontal deflections of magnetic field lines, which are smaller in magnitude than the vertical deflections. The role of the Coriolis force is to twist rising magnetic loops against galactic rotation and to produce, in cooperation with vertical gravity and resistivity, the helical magnetic field structures apparent in simulations of Paper I.

In order to integrate MHD Eqs. (1)-(3) we apply the ZEUS-3D code (Stone & Norman 1992a, 1992b) with our own modifications necessary to introduce the Coriolis force and to incorporate the localized resistivity (see Paper I for more details).

2.2 The method of measuring of the effect of the resistivity

The presence of the resistivity in the induction equation releases the magnetic freezing condition and leads to a topological evolution of magnetic field lines.

Measurements of deviations of a single magnetic field line from the perfect freezing in gas motion is possible through integration of magnetic vector field in a 3D Cartesian volume bounded by periodic boundaries. We propose a method that is illustrated in Fig. 1. We assume that the initial magnetic field is parallel to the Y-axis in the whole volume of the computational box. We apply periodic boundary conditions on the XZ as well as YZ boundaries. We integrate magnetic field lines between the starting XZ boundary denoted by the symbol "S" and the ending XZ boundary denoted by "E". We consider a family of N magnetic field lines. Each member of the family is denoted by the index i. Our choice of initial magnetic field parallel to the Y direction ensures that magnetic field lines are closed by periodic boundaries, i.e. the intersection points  (Xi,Zi)S of the ith line in the S plane are identical as the intersection points (Xi,Zi)E in the E plane.

In the case of ideal MHD the magnetic freezing condition means that a given element of the fluid will remain on the same magnetic flux tube for an arbitrarily long period of time. This means that if (Xi,Zi)S = (Xi,Zi)Eat t=0 then (Xi,Zi)S = (Xi,Zi)E at any later time (instant), no matter what the evolution of the box interior is. If, on the other hand, the system is non-ideal i.e. resistivity is present, then magnetic field lines become open $(X_i,Z_i)_{S} \neq (X_i,Z_i)_{E}$ at t > 0. We define deviations of the intersection points of magnetic lines with the S and E boundaries, as $(\Delta
X_i,\Delta Z_i) = (X_i,Z_i)_{E} - (X_i,Z_i)_{S}$ (see Fig. 1).

\end{figure} Figure 1: Figure describing method of measuring the dispersion of magnetic field lines due to resistivity. S and E are boundary plane at beginning and ending of y coordinate. PS and PE are points on the XZ plane where magnetic line starts and ends. $\Delta {x_i}$, $\Delta {z_i}$, $\Delta {r_i}$ is a deviation between the starting and the ending points of the ith line. The right panel presents the projection of the deviation on the XZ plane.
Open with DEXTER

As soon as $(\Delta X_i,\Delta Z_i)$ are known as functions of time one can perform a statistical analysis through the computation of the mean and the standard dispersion of deviations of the E and S intersection points.

           $\displaystyle \overline{\Delta x}$ = $\displaystyle \frac{1}{N} \sum_{i=1}^N \Delta x_{i},$ (6)
$\displaystyle \overline{\Delta z}$ = $\displaystyle \frac{1}{N} \sum_{i=1}^N \Delta z_{i},$ (7)

                           $\displaystyle \sigma_{x}$ = $\displaystyle \left( \frac{1}{N-1} \sum_{i=1}^N{(\Delta x_{i}
- \overline{\Delta x})^2}\right)^{1/2} ~ ,$ (8)
$\displaystyle \sigma_{z}$ = $\displaystyle \left( \frac{1}{N-1} \sum_{i=1}^N{(\Delta z_{i}
- \overline{\Delta z})^2}\right)^{1/2} ~ .$ (9)

In the rest of the paper, for the sake of brevity, this kind of dispersion will be referred to as the dispersion of magnetic field lines. In actual implementation of the method of determination of the magnetic diffusion the integration of magnetic field lines is performed with the aid of the Runge-Kutta method of 4th order. This procedure is applied for a large sample of staring points PS on the first XZ boundary for different timesteps of evolution of the Parker instability.

The time evolution of the statistical moments depends primarily on values of the resistivity coefficients $\eta _1$ and $\eta _2$, as well as on the grid resolution adopted for the MHD simulations. The presented method makes it possible to determine the intrinsic numerical resistivity of the MHD code by means of comparison of the dispersion measure of magnetic field lines for low resolution runs with the vanishing resistivity and high resolution runs with different values of the resistivity coefficient $\eta _1$. Results of this procedure will be presented in Sect. 3. Finally we would like to stress out that the periodic boundary conditions are essential for our method because we measure the departures of magnetic field lines from the closure (by periodic boundaries). This measure indicates topological changes of field lines and therefore the presence of non-ideal effects.

2.3 Initial conditions

Our initial equilibrium state is an exponentially stratified isothermal disk (Parker 1966; Paper I). The initial equilibrium is characterized by an uniform vertical gravity with the vertical gravitational acceleration $g = -2
\times 10^{-9}~{\rm cm ~ s^{-1}}$, the midplane gas density $n_0= 1~{\rm cm}^{-1}$ and the sound speed $c_{\rm s} = 7~{\rm km ~ s^{-1} }$, constant across the disk. The fixed ratio of the magnetic pressure to the gas pressure is given by $\alpha = p_{\rm mag} / p_{\rm gas}$, so that the total pressure is given by $p_{\rm tot} = (1+\alpha) p_{\rm gas}$. We assume that $\alpha = 1.0$, corresponding to the magnetic field strength of $4.5~{\rm\mu
G}$ at the galactic midplane. The initial magnetic field is parallel to the azimuthal Y-axis. The dependence of equilibrium quantities on z is given by

 \begin{displaymath}\frac{n_0(z)}{n_0(0)} = \frac{p_0(z)}{p_0(0)} =
\frac{B_0^2(z)}{B_0^2(0)} = \exp \left( - \frac{z}{H} \right) ~ ,
\end{displaymath} (10)

where $H = \left( 1 + \alpha \right) ~ c_{\rm s}^2 / \vert g\vert \approx 159 ~ {\rm
pc}$ the vertical scale-height. We apply the angular velocity of the galactic rotation $\Omega = 0.025 ~ {\rm Myr}^{-1}$ corresponding to that of the Solar orbit.

In order to excite the Parker instability we apply initial perturbations of vertical velocity component in the following form

= \left\{\begin{array}{lll}
v_0 \cos \left( \frac{2 \...
... z\vert \leq H\\
0 & {\rm elsewhere} \\
\end{array} \right.
\end{displaymath} (11)

where nx, ny are the numbers of the harmonic components of the X, Y directions, respectively. Lx, Ly are the physical sizes of the simulation domain and the initial velocity amplitude is $v_0 = 1.0~{\rm km~s^{-1}}$. Following Paper I the initial velocity perturbations are nonvanishing only up to the height H above the disk midplane in order to take into account the fact that the perturbations originate from the disk activity, i.e SN explosions, stellar winds, etc. For the sake of simplicity we perturb only the vertical velocity component.

2.4 Numerical setup and input parameters

We analyze a set of numerical simulations in a domain of physical sizes of  $L_{x} = 600 ~ {\rm pc}$, $L_{y} = 1800 ~ {\rm pc}$, $L_{z}
= 1200 ~ {\rm pc}$ in the Cartesian reference frame x, y, z corresponding locally to the radial, azimuthal and vertical galactic coordinates. Given the fact that the scaleheight of the disk is 159 pc, the density contrast between the bottom and the top of the computational box is about 1900.

The initial perturbations follow the formula (11) with $n_{x}
\times n_{y} = 3 \times 3$ in all cases. The periodic boundary conditions are applied at all boundaries perpendicular to the galactic plane (at the XZ planes and the YZ planes). At the lower and upper boundary parallel to the galactic plane (the XY planes) the reflecting boundary conditions are applied. We note that the reflecting boundary conditions imposed on horizontal boundaries are formally not compatible with the vertical stratification of By. However, this choice of boundary conditions allows for a stable evolution of the system on long timescales without any significant contribution of numerical artifacts. Moreover, the integrated magnetic field lines are located in the middle range of the full height of the computational box, i.e. the starting points of magnetic field lines span the range $240~{\rm pc} \leq z \leq
960~{\rm pc}$ out of the full height of the computational box.

In our study we adopt several different values of the resistivity parameters $\eta _1$, $\eta _2$ and two different grid resolutions: $30\times 90 \times 60$ (low) and  $60\times 180\times 120$ (high) as listed in Table 1. We assume  ${\rm pc^2 ~Myr^{-1}}$ as the unit of $\eta$.

Table 1: Input parameters for computed models.

2.5 Statistical analysis of the simulation results

For each model displayed in Table 1 we compute a time evolution of the following six quantities:

the maximum of vertical velocity $v_{z ~ {\rm max}}$ representing the amplitude of unstable Parker modes;

the ratio of mean energy densities of total and uniform magnetic fields in the disk ($z\leq H$),

\begin{displaymath}U_{ B_{\rm t}}/U_{ B_{\rm u}}=~<B^2>/<B>^2.

The averaging denoted by <...> is performed over a volume between z=0 and z=H;

the maximum resistivity $\eta_{\rm max}$ in the computational volume;

the maximum of current density $j_{\rm max}$ in the computational volume;

the dispersion $\sigma_x$ of magnetic field lines in the X direction;
the dispersion $\sigma_z$ of magnetic field lines in the Z direction.

3 Results

All simulations presented in this paper start from the same initial state as described in Sect. 2.3. The basic difference between different considered models listed in Table 1 follows from the choice of different values of resistivity parameters $\eta _1$ and $\eta _2$. The identical initial conditions lead to initially identical evolution, however the differences become apparent as soon as resistivity starts to play a role.

\end{figure} Figure 2: Cross-sectional slices of the computational domain at t = 700 Myr corresponding to the end of the exponential growth phase. Panel  a) shows magnetic vectors and the greyplot of density perturbation $\Delta \rho /\rho _0$ in the YZ plane, panel  b) shows velocity vectors and the density perturbation in the XZ plane and panel  c) shows magnetic vectors and current density in the same plane.
Open with DEXTER

An example of cross-sectional slices at the end of exponential growth phase for the "ideal" model p0 ( $\eta_1 = \eta_2=0$) is shown in Fig. 2. Panel a displays magnetic field vectors on the background of density perturbation  $\Delta \rho /\rho _0$, where $\rho_0$ is the density distribution at t=0. The initially horizontal magnetic field lines are deformed in the familiar way by an undulatory mode of the Parker instability. Gas condensation of magnitude  $\Delta \sim \rho_0$ is forming in the magnetic valley. Panel b shows the density perturbation and velocity vectors in the XZ plane. It is remarkable that the falling gas forms dense vertical sheets between rising magnetic lobes. Panel c displays the XZ-components of magnetic vectors superposed onto the current density greyplot. As it is apparent, the current sheets form between rising magnetic lobes and coincide with the dense vertical gas sheets. The cross-sectional slices for other models (except the cases of large values of uniform resistivity) remain very similar to the presented one at time T = 700 Myr.

We refer the reader to Paper I for more a detailed discussion of different initial configurations, which include cases both with and without the explicit resistivity. As it has been shown in Paper I, in the presence of explicit localized resistivity reconnection starts to operate in current sheets. Before the onset of reconnection, Parker instability deforms the initially horizontal magnetic field into magnetic arches which are additionally twisted by the Coriolis force. The reconnection process transforms the twisted arches into helical magnetic flux tubes. The adjacent tubes undergo next phase of reconnection leading to their coalescence. The final states contain a significant radial magnetic field component, which is produced as a result of Parker instability, the Coriolis force and magnetic reconnection.

In the absence of explicit resistivity our system starts to reconnect as well, at later times, as a result of numerical resistivity. A comparison between cross-sectional slices for two cases, with and without the explicit resistivity is presented in Figs. 2b and 3a of Paper I.

In order to present the results of our statistical analysis of the perturbed configuration we display the six statistical quantities (a)-(e), described in Sect. 2.5, in Figs. 3-6. The plots in Figs. 3 and 4 show the evolution of these quantities for experiments listed in Table 1. i.e. for different values of uniform resistivity coefficient $\eta _1$ and $\eta_2=0$. The results for low resolution runs p0, q0, r0, s0, t0 are presented in Fig. 3 and the results for the high resolution runs P0, Q0, R0, S0, T0 are presented in Fig. 4. Next two figures present results for different values of $\eta _2$ (the localized resistivity coefficient) and $\eta_1=0$. The results for low resolution runs p0, p1, p2, p3, p4 are presented in Fig. 5 and the results for high resolution runs P0, P1, P2, P3, P4 are presented in Fig. 6.

\end{figure} Figure 3: Evolution of statistical quantities for low resolution. Dependence on $\eta _1$. (solid - p0, dotted - q0, dashed - r0, dash dot - s0, dash dot dot - t0).
Open with DEXTER

\end{figure} Figure 4: Evolution of statistical quantities for high resolution. Dependence on $\eta _1$. (solid - P0, dotted - Q0, dashed - R0, dash dot - S0, dash dot dot - T0).
Open with DEXTER

\end{figure} Figure 5: Evolution of statistical quantities for low resolution. Dependence on $\eta _2$. (solid - p0, dotted - p1, dashed - p2, dash dot - p3, dash dot dot - p4).
Open with DEXTER

\end{figure} Figure 6: Evolution of statistical quantities for high resolution. Dependence on $\eta _2$. (solid - P0, dotted - P1, dashed - P2, dash dot - P3, dash dot dot - P4).
Open with DEXTER

3.1 Models with different values of the coefficient $\eta _1$

In general, all statistics are strongly dependent on the value of the resistivity coefficient $\eta _1$ (except the maximum resistivity, which is not representative for the present set of models). The evolution of $v_{z ~ {\rm max}}$ presented in Figs. 3a and 4a demonstrates a common existence of three phases: the initial phase, the exponential growth phase and the saturation one.

During the initial phase, $v_{z ~ {\rm max}}$ decreases because the initial perturbation is a mixture of stable, decaying and unstable modes. The stable modes are quickly radiated out from the perturbation region, the decaying modes are dumped out and then the standing waves of unstable modes start to dominate at the end of the initial phase.

The second phase of exponential growth, belonging to the linear phase of the growth of the Parker instability, starts at about 400 Myr and lasts about 250 Myr. This period is characterized by an exponential rise of the  $v_{z ~ {\rm max}}$ by factor of about 100 (from 0.1 to  $10.0~\rm {km~s^{-1}}$). The inclination of the curve in this phase gives a good qualification of the evolution growth rate. For models p0, q0, r0, s0, P0, Q0, R0, S0 ( $\eta_1 = 0.0
\div 10.0$), the growth rate is about 0.036, but for the cases t0 and T0 with  $\eta_1 = 100.0$ it is only 0.009. The maximum of vz reaches a level of about  $10.0~\rm {km~s^{-1}}$, comparable to the Alfvèn speed, at the beginning of the saturation phase and does not exceed it later on. For the largest value of $\eta _1$ the evolution looks different, however we can still distinguish three phases. Values of  $v_{z ~ {\rm max}}$ in initial phase are larger, the growth rate is smaller, as mentioned earlier, and in the saturation phase, the maximum velocity value decreases faster (but only for the low resolution). In this case the presence of high resistivity apparently slows down the linear phase of growth of the Parker instability.

Figures 3b and 4b show the evolution of the ratio of mean energy densities of the total and uniform magnetic fields in the disk $U_{B_{\rm t}}/U_{B_{\rm u}}$ for the low and high resolution respectively. During the period of the first 500 Myr, the ratio is 1.0 in both cases, which means that the turbulent magnetic component is negligible. Then the ratio starts to rise very rapidly to the maximum value of about 1.15 and later on it decreases. We can notice that values of  $U_{B_{\rm t}}/U_{B_{\rm u}}$ are strongly dependent on the uniform resistivity. The experiments with the low and high number of the grid points show that both maximum values and the position of the curve maxima are different for models with different values of the resistivity coefficient $\eta _1$. The evolution of  $U_{B_{\rm t}}/U_{B_{\rm u}}$ in models with low resolution and with the lowest resistivity (p0, q0) have maximum values of 1.18 at the position of 1100 Myr. Increasing $\eta _1$ to 1.0 (r0) and 10.0 (s0) lowers the maximum value to 1.13 but positions of the maxima are different for these two cases: r0 has its maximum at about 1200 Myr, while s0 at 900 Myr. We can see that there is no rule that increasing the resistivity causes a shift of the maximum position to the right or the left side, both are possible. The last experiment with the highest resistivity t0 ( $\eta_1 = 100.0$) has a maximum of about 1.025 more or less at the same time step as for the previous experiment. The high resolution simulations result in a slightly different evolution of the ratio  $U_{B_{\rm t}}/U_{B_{\rm u}}$. Due to the lower value of the numerical resistivity, the maximum for the model Q0 ( $\eta_1=0.1$, see Table 1) is slightly higher than for the case P0 without physical resistivity at all. The maxima for both curves are positioned at about 950 Myr, so shifted to the left in comparison with the low resolution calculations. The run R0 shows two maxima (about 900 and 1100 Myr) and then the evolution is slightly different increasing its value until the end of the evolution. The experiments with the higher resistivity coefficient (S0 and T0) are similar to models with the low number of grid points (s0, t0, Fig. 3b. This means that really high physical resistivity results in that the numerical one is insignificant.

In case of uniform resistivity, plots of Figs. 3c and 4c are presented to ensure of compatibility with the case of the current-dependent resistivity discussed in the next Subsection.

The evolution of $j_{\rm max}$ (Figs. 3d and 4d) shows similar character as the ratio  $U_{B_{\rm t}}/U_{B_{\rm u}}$. In the beginning the maximum of the current density is very small (about 0.03), but after 500 Myr it starts to grow rapidly to the maximum value of about 0.2 and 0.5 for the low and high resolution, respectively. This difference is caused by the construction of the current density calculations, which depends on a size of the grid step. In time less than 200 Myr, $j_{\rm max}$ reaches a maximum, which is very dependent on the uniform resistivity. The evolution of the maximum value of the current density in computations with low grid point number gives a very similar value of the curve maxima (0.2) for all cases, except the case t0 with $\eta_1 = 100.0$ showing only 0.03. The simulations with high resistivity show a larger dispersion of maximum values, three cases (P0, Q0, R0) have maximum value between 0.4 and 0.5, while the experiment S0 ( $\eta_1=10.0$) has 0.27 and T0 ( $\eta_1 = 100.0$) has 0.03. In the saturation phase, the value of  $j_{\rm max}$ decreases. The inclination of all curves in this phase is independent of the uniform resistivity for $\eta_1 < 10.0$ (runs p0, q0, r0), but for the higher one some nonlinear effects may play a significant role and the inclination is different. For example, the case s0 with  $\eta_1=10.0$ (see Fig. 3d) contains three sub-phases in the saturation state: from time 700 Myr, 900 Myr and from 1250 Myr, where the inclination is different.

The deviation of magnetic field lines between opposite XZ boundaries is presented in Figs. 3e, 3f and 4e, 4f for the whole set of runs (see Table 1) corresponding to the constant resistivity. At the end of the exponential growth of perturbations the field lines disperse quickly, but later on the resistivity limits further dispersion. We separated the standard deviation of magnetic field lines into the X direction  $\sigma_x$ (Figs. 3e and 4e) and Z direction  $\sigma_z$(Figs. 3f and 4f). For both resolutions, the magnetic field lines are dispersed more strongly in the Z direction than in X direction, depending on the magnitude of the coefficient $\eta _1$. This can be explained by the dominance of the vertical gravity force over the horizontally directed Coriolis force. This relation between basic forces leads to vertical current sheets as displayed in Fig. 2c. The vertical current sheets lead to the dominating vertical dispersion of reconnected magnetic field lines. The maximum values of  $\sigma_x$ for experiments with zero or the smallest resistivity (p0, q0, P0 and Q0) are about 80 pc for the high resolution models and about 100 pc in the case of the low resolution runs. Values of  $\sigma_z$ are much larger: 250 pc for the high resolution and about 150 pc for the low resolution runs. This behavior is consistent with the expectation that lower resolution of simulations introduces larger numerical resistivity. The maximum values of both dispersions decrease strongly with increasing resistivity, resulting in 50 pc dispersion for the case t0 and T0 (resistivity 100.0, see Table 1). Much smaller displacements  $\sigma_x$ in the horizontal direction are related to the fact that the primary force responsible for the instability is the vertical buoyancy force.

3.2 Models with different values of $\eta _2$

Figures 5 and 6 present the time evolution of our six statistical quantities in simulations with different values of $\eta _2$(the low resolution cases: p0, p1, p2, p3, p4 and the high ones: P0, P1, P2, P3, P4, see Table 1). Generally, the influence of the localized resistivity on the statistical quantities (a)-(f) is much weaker (except  $\eta_{\rm max}$) than the influence of the uniform resistivity.

The time evolution of $v_{z ~ {\rm max}}$ proceeds similar to the case of the uniform resistivity. In the present case, however, the evolution of the exponential growth phase is identical for all runs. This is because the high currents, localized in current sheets, starts to exceed the critical current (may change for $j_{\rm crit}=0.1$) in the middle of the exponential growth phase. In the saturation phase, the maximum of vz falls down less than in the uniform resistivity case.

Comparing to the uniform resistivity (Figs. 3b and 4b) the evolution of  $U_{B_{\rm t}}/U_{B_{\rm u}}$ (Figs. 5b and 6b) looks different especially for the high resolution cases. This quantity exhibits a strong dependence on $\eta _1$ but the effect of localized resistivity is only moderate. Now, the ratio  $U_{B_{\rm t}}/U_{B_{\rm u}}$ reaches comparable values for both grid resolutions, while in the case of the uniform resistivity the resolution influences the maximum values significantly.

For the high resolution cases  $j_{\rm max}$ grows to large values, especially for smaller values of $\eta _1$ (Figs. 5d and 6d). For example, for $\eta_1=0.0$ and  $\eta_2=0.0$ $j_{\rm max}$ reaches the value of 0.2 and 0.5 for the low and high resolution, respectively. This difference we discussed in Sect. 3.1. For a large  $\eta_2=100.0$ the evolution of  $j_{\rm max}$ rises only to the value of 0.2 i.e. exceeding the critical current much less than in the case of small $\eta _2$.

The standard deviations of the magnetic field lines  $\sigma_x$ and  $\sigma_z$ (Figs. 5e, 5f and 6e, 6f) do not exhibit an apparent dependence on the resistivity at the small grid resolution, however at the high resolution the behavior is opposite with respect to the uniform resolution. Up to  $t \sim 900$ Myr all three kinds of the resistivity produce the same dispersion of magnetic field lines. Contrary to the other cases we note a systematic growth of  $\sigma_x$ in the case of $\eta_2=100$ for the high resolution. The dispersion of magnetic field lines in the X direction reaches about 100 pc at t=1500 Myr, while  $\sigma_x$ is only 60 pc for $\eta_2=0$. This property makes a significant difference in the localized resistivity with respect to the uniform resistivity, because it is apparent that the localized one does not tend to eliminate the initial dispersion of magnetic field lines, which appeared in early phases of operation of the resistivity in the period between t=600 and t=900 Myr.

\end{figure} Figure 7: Averaging of standard deviation in Z direction near maximum for models with different values of $\eta _1$. Thick lines indicate region of averaging and computed mean value of  $\bar{\sigma_{z}}$. Boxes around lines indicate errors. (solid - low resolution, dashed - high resolution).
Open with DEXTER

3.3 Determining of the numerical resistivity

In our experiments we try to roughly estimate the value of the numerical resistivity. The Parker instability does not grow infinitely but is suppressed by the resistivity after exceeding some level, manifested by reaching the saturation phase in the evolution of our six statistical quantities (see Sect. 3). The values of these quantities do not rise later on, but after reaching some maximum point, oscillate around it, or decrease exponentially. As an indicator of the total diffusivity in our models, we have chosen the dispersion  $\bar{\sigma_{z}}$ of magnetic field lines defined in Sect. 2.2. We average in time the values of this quantity near maximum, i.e. the averaging is done for a period of about 400 Myr around maximum (see Fig. 7). The computed values of  $\bar{\sigma_{z}}$ are listed together with errors  $\Delta \bar{\sigma_{z}}$ in Table 2.

Table 2: Mean values of maximum of standard deviations in Zdirection for different models averaged from 800 Myr to 1200 Myr.

For the uniform resistivity models p0, q0, and r0 (see Table 2, low resolution), $\bar{\sigma_{z}}$ has the same value of about $150.0 \pm
\Delta \bar{\sigma_{z}}$ pc. Independently of the magnitude of explicit resistivity, the time evolution of  $\sigma_z$ is similar for the runs p0, q0, and r0. This fact indicates that numerical resistivity is large compared to $\eta _1$. For $\eta_1=10.0$ (model s0) the situation is different. The explicit resistivity is now dominant.

Comparing all the above with the high resolution models (P0, Q0, R0, and S0) we conclude that the numerical diffusion in high resolution is much smaller, because  $\bar{\sigma_{z}}$ is very large for the model with no explicit resistivity and its value decreases with growing $\eta _1$.

Although we have very few points, we can roughly determine the numerical resistivity in low resolution through plotting a log-lin graph of dependency $\bar{\sigma_{z}} = A + B \log( \eta_1 )$ for high resolution points (see Fig. 8). We obtain values A=146.3 and B=-31.9. Now we take the value of $\bar{\sigma_{z}}=142.6$ for model p0 (the low resolution, no resistivity) and we project it onto our the curve. We can read the value of $\eta_1 = 1.3$ which corresponds to $\bar{\sigma_{z}}=142.6$. This means that to get evolution of the Parker instability in the high resolution domain similar to model p0 in the low resolution, we have to set $\eta_1 = 1.3$, because this value corresponds to numerical resistivity in low resolution. It is in agreement with our rough estimate above.

We can go further and include the error $\Delta \bar{\sigma_{z}} = 16.2$ in our estimations. By going through the same procedure as described above, but now for values $\bar{\sigma_{z}} + \Delta \bar{\sigma_{z}} = 158.8$ and $\bar{\sigma_{z}} - \Delta \bar{\sigma_{z}} = 126.4$, we estimated a possible range for the value of the numerical resistivity, which is $\eta_{\rm 1
~ numerical} \in (0.4,4.2)$ (see Fig. 8).

The procedure described above does not exclude the numerical resistivity in the high number of grid points experiments. After performing a set of simulations in higher resolutions we should be able to plot a dependency of the numerical resistivity on the resolution. Then it will be possible to determine the value of  $\eta_{\rm 1 ~ numerical}$.

\end{figure} Figure 8: Interpolation of $\eta _{\rm numerical}$.
Open with DEXTER

4 Discussion

The central issue of the present project is to identify the effects of three kinds of resistivity (uniform, current dependent and numerical) on the topological evolution of magnetic field lines after the onset of the Parker instability. We found that the numerical resistivity present in the currently-used code Zeus3D produces dispersion of magnetic field lines, which is present in the action of a weak uniform resistivity.

This finding allowed us to propose a new method of determination of the numerical diffusion in the numerical simulations of the Parker instability. The method described in Sect. 2.2 is based on a statistical analysis of the dispersion of magnetic field lines in a non-ideal medium.

We measure the dispersion of a large set of magnetic field lines and relate this dispersion to the type and the magnitude of the resistivity. This kind of analysis is applicable to MHD simulations in a local approximation involving periodic boundary conditions at least on two planes which are intersected simultaneously by a family of magnetic field lines.

We determine the numerical resistivity for low resolution simulations comparing the mean dispersion of magnetic field lines in the case of the vanishing physical resistivity to the case of high resolution and several values of the uniform resistivity.

The evolution of our statistical quantities shows that for the largest values of the resistivity and higher resolution of simulations the numerical resistivity is negligible. In some cases the extremely large uniform resistivity impedes the Parker instability.

Our considerations imply that the Parker instability contributes to a randomization of galactic magnetic fields. The main sources of perturbations in real galaxies are supernovae supplying locally kinetic, thermal and cosmic ray energy. These perturbations are the most plausible mechanism of the excitation of the Parker instability and the randomization of the magnetic field in spiral arms. The resistivity is then expected to contribute to the annihilation of the random magnetic field component in the inter-arm regions. In galaxies like NGC 6946 and M 51 the magnetic field appears to be much more uniform in the inter-arm regions then in arms (see model of Rohde et al. 1996).

The modeling of the Parker instability in the present simulations is far from being realistic. In order to make our model feasible we applied simplifying assumptions of a plane-parallel initial magnetic field, uniform vertical gravity, isothermal evolution of the thermal gas, a lack of cosmic rays, a lack of selfgravity and differential rotation in the disk. All these circumstances imply that numbers and timescales of different phenomena are not realistic too. Among different effects we can expect that the presence of cosmic rays significantly shortens all the relevant timescales of observed phenomena, however we believe that future incorporation of all the complexity of realistic galactic disks will not change the qualitative picture of the evolution of a magnetic field structure.

5 Summary and conclusions

In the present paper we considered the general problem of the influence of the physical and numerical resistivity on the Parker instability. We focused on the issue of the dispersion of magnetic field lines (which according to Sect. 2.2 is understood in this Paper In the topological sense), resulting from three types i.e. uniform, current-dependent and numerical resistivities.

We performed several experiments applying two different grid resolutions and different types and magnitudes of resistivity. We propose a new statistical method of quantitative estimation of the influence of the explicit and the numerical resistivity on the dispersion of magnetic lines and the amount of a random component of the magnetic field in the considered MHD system.

Our conclusions are:

Except for the most extreme value of the uniform resistivity $\eta_1=100$, the resistivity does not influence the growth rate of the Parker instability. For extremal values of the resistivity, the linear growth rate of Parker instability is significantly diminished.

The time evolution of the dispersion of magnetic field lines and the ratio of the total to uniform magnetic field is strongly dependent on the type and the magnitude of resistivity.

The presence of all kinds of resistivity leads to a rapid dispersion of magnetic field lines at the end of the exponential growth phase of evolution of the Parker instability. After the initial sudden dispersion of magnetic lines the subsequent evolution depends on the type of resistivity:

- uniform resistivity tends to convert the disturbed magnetic field structure to the uniform state (uniformization of the magnetic field). This effect is clearly seen in the plots of the ratio of total to uniform magnetic energies and in the plots of the dispersion of magnetic field lines;

- the localized resistivity diminishes the ratio of the total to uniform magnetic energies less efficiently than the uniform resistivity. On the other hand, for high resolution runs we observe the growth of the horizontal dispersion of magnetic field lines in the course of time, while the vertical dispersion tends to decrease;

- the effect of the numerical resistivity is dependent on the grid resolution and is qualitatively similar to the effect of uniform resistivity. We determine the magnitude of the uniform resistivity equivalent to the numerical resistivity in simulations with lower resolution.

The method can be applied to compare the numerical magnetic diffusivity of different codes and studies of different physical processes e.g. turbulence.


We thank Polish Committee for Scientific Research (KBN) for the support through the grants PB 4264/P03/99/17 and PB 249/P03/2001/21.


Copyright ESO 2003