A&A 445, 915-929 (2006)
DOI: 10.1051/0004-6361:20053582

Dynamo coefficients in Parker unstable disks with cosmic rays and shear

The new methods of estimation

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


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 7 June 2005 / Accepted 18 July 2005

Abstract
We report a three-dimensional (3D) numerical model of the Parker instability with cosmic rays and shear in a galactic disk and compute the components of the dynamo coefficients $\alpha $ and $\beta $ from the electromotive forces obtained from simulations. For the first time we apply well-known statistical procedures to fit the modeled electromotive forces in terms of the dynamo tensors. We compare our results with kinematic and other methods of dynamo coefficient estimation. Although we do not solve all problems with the physical interpretation of the dynamo tensors we find that the presented methods give us the ability to compute them with the desired statistical quality. The obtained values of $\alpha _{xx}$ (along the radial direction in galaxy) and $\alpha _{zz}$ are of the order of a few km s-1 while the value of $\alpha _{yy}$ (along the azimuthal direction) is one or two orders of magnitude smaller.

Key words: galaxies: ISM - galaxies: magnetic fields - methods: numerical - methods: data analysis

   
1 Introduction

Classical dynamo theory (e.g. Parker 1979; Ruzmaikin et al. 1988; Moffat 1978) holds that the large-scale magnetic field in spiral galaxies is amplified due to the helical turbulent motion of interstellar gas (ISM) and the shear of differentially rotating disks. The crucial element of the theory is the assumption that gaseous velocity and magnetic field can be separated into small and large scales. The theory yields two dynamo coefficients: $\alpha $ and $\beta $, which increase and diffuse the magnetic field, respectively. The dynamo theory explains the observed magnetic field topology in barred and spiral galaxies (e.g. Urbanik et al. 1997; Panesar & Nelson 1992; Moss et al. 2001; Otmianowska-Mazur et al. 2002) quite well but there are still theoretical and observational constraints which are hard to clarify within the framework of classical models (reviewed in Widrow 2002). First, to obtain magnetic field intensity values that are observed in galaxies using classical dynamo theory, we arrive at an enormous amount of time in comparison with the life time of galaxies (e.g. Brandenburg 2001; Widrow 2002). The observation of magnetic fields in objects with high red-shifts (see Beck et al. 1996; Widrow 2002) and in, e.g., irregular galaxies (Chyzy et al. 2000; Otmianowska-Mazur et al. 2000) showed that the amplification of a magnetic field is faster than the time associated with the action of a turbulent dynamo.

The first analytical estimations of $\alpha $ using spectral methods and incorporating the Lorentz force in the absence of the large-scale magnetic field component were done by Frisch et al. (1975) and Pouquet et al. (1976). Later et al. (1999) and Kleeorin et al. (2002, 2003) analyzed anisotropic turbulence and found that the magnetic part of $\alpha $ is determined by the dynamical helicity equation. The authors also found that both algebraic (i.e. quenching of both the $\alpha $ and $\beta $ coefficients) and dynamic nonlinearity limited the growth of the magnetic field at the equipartition level.

Externally forced numerical simulations of a local cube that take into account the back-reaction of a magnetic field on turbulent motion (so called quenching, e.g. Cattaneo & Vainshtein 1991; Vainshtein & Cattaneo 1992; Cattaneo 1994; Cattaneo & Hughes 1996; Ziegler et al. 1996) showed that the Lorentz force strongly suppressed the turbulent dynamo action. One of the aspects of the back reaction problem is related to the conservation of the total magnetic field helicity in media with high magnetic Reynolds numbers $R{_{\rm m}}\gg1$ (e.g. Berger & Field 1984; Brandenburg et al. 2002; Brandenburg & Subramanian 2005). The working dynamo models based on helical motion produce the large-scale magnetic field helicity and prevent the large-scale dynamo action from working on a dynamical time scale (e.g. Brandenburg et al. 2002). There are two ways to solve this problem: constant ejection of the magnetic helicity flux through boundaries (Blackman & Field 2000 and Kleeorin et al. 2000) or creation of helicity with opposite signs at large and small scales (see Brandenburg et al. 2002; Kleeorin et al. 2002). Alternatively, Blackman & Field (2000) analyzed the role of boundary conditions in modeling the magnetic back reaction on turbulent motion. They found that the strong suppression of $\alpha $ obtained by Cattaneo & Hughes (1996) in their numerical simulations was caused by the assumption of periodic boundary conditions. Unfortunately, calculations concerning a local cube with external forcing and open boundaries failed to solve this problem (e.g. Brandenburg & Dobler 2001; Brandenburg & Sokoloff 2002, hereafter BS02). Brandenburg & Dobler (2001) showed that the growth of the large-scale magnetic field in their model was strongly quenched by loss of the magnetic helicity through the open boundary. For these reasons, it seems that calculations with forcing are no longer valid and the crucial matter is to model realistic physical processes (e.g. Parker or MRI instabilities). Such processes produce small-scale motions of gas (e.g. BS02) and provide an opportunity to compute the dynamo coefficients ($\alpha $ and $\beta $) from the electromotive force (EMF, e.g. BS02). Similar calculations were done by Ziegler et al. (1996), however the obtained values of the $\alpha $ coefficients (a few meters per second) were five orders of magnitude smaller than the classical turbulent dynamo coefficient values in galaxies.

On the other hand, Vishniac & Cho (2001) showed that it was possible to build a dynamo model ( $\alpha-\omega$ type) that conserved magnetic helicity. They estimated that the energy of the regular magnetic component in galaxies could grow two orders of magnitude faster than in all other calculations (Brandenburg 2001; Kleeorin et al. 2000). However, Arlt & Brandenburg (2001) did not find any supportive evidence that their dynamo worked. In their next model of the incoherent dynamo in accretion disks, Vishniac & Brandenburg (1997) obtained $\alpha_{\phi\phi}=0$, but the $\alpha-\omega$ dynamo still worked efficiently. The large-scale simulations of Moss et al. (1999) showed that it was possible to obtain a working dynamo in a model by taking into account additional helicity provided by the buoyancy of the large-scale galactic magnetic field component.

In order to solve the mentioned above problems of classical dynamo theory (e.g. too long of a time scale for the large-scale magnetic field growth, limit of the large-scale magnetic field growth by the loss of the magnetic helicity through the open boundary calculations, etc.), the idea of a cosmic-ray driven fast dynamo was first introduced by Parker (1992) and then explored by Hanasz & Lesch (1993, 1997, 1998, 2000). Our 3D numerical simulations of the Parker instability evolution under the influence of cosmic rays (CR) and differential rotation of the underlying disk (Hanasz et al. 2004) showed that amplification of the large scale magnetic field was possible on a timescale of 250 Myr, which was comparable to the period of galactic rotation.

Hanasz et al. (2002) investigated the role of the Parker instability and the process of a magnetic reconnection in uniformly rotating galactic disk. In our next paper (Otmianowska-Mazur 2003) we derived the kinematic (hydrodynamic) dynamo coefficients by integrating their values from small-scale velocity fields obtained in our previous calculations (Hanasz et al. 2002; see also Kowal et al. 2003a,b). The resulting values of $\alpha $ and $\beta $ were of the order of 10 km s-1 and 1025 cm2/s, respectively. Unfortunately, the main component $\alpha_{\phi\phi}$, which is the most important component in the classical $\alpha-\omega$ dynamo, was very small.

The aim of the present work is to compute the components of the dynamo tensors using various methods. We demonstrate that the methods applied by different authors lead to different results. Thus, a natural question arises: which method of computation of the dynamo coefficients is more reliable? We propose the following verification method. The dynamo coefficients are intended to serve as parameters describing the cumulative action of small-scale turbulent motion on the large scale magnetic field. In other words, the information on magnetic field amplification is expected to be contained in a finite set of numbers or functions of one or more coordinates. Therefore, a simple test of the validity of the methods for dynamo coefficient estimation, for numerical MHD experiments showing dynamo action, can be based on the comparison of the observed growth of the mean magnetic field in the MHD (or MHD+CR) experiment to the growth resulting from solution of the dynamo equations. The dynamo coefficients computed in the same experiment are treated as input parameters. In an equivalent formulation, one could compute electromotive forces on the basis of simulation output first, then dynamo coefficients, and finally verify that the dynamo coefficients multiplied by the respective components of the mean magnetic field reproduce the electromotive force. If the electromotive force reconstructed from the dynamo coefficients does not match the directly computed electromotive force, then some essential information on the dynamo process is missing in the derived form of the dynamo coefficients.

We shall demonstrate that the standard methods for computation of the dynamo coefficients fail to reproduce the original electromotive force as well as the amplification rate of the large-scale magnetic field. Therefore we attempt to find new approaches and to modify existing methods of computing the dynamo coefficients. We examine the local and non-local formulations by Ziegler et al. (1996) and Brandenburg & Sokoloff (2002) with and without modification, with the kinematic method following Otmianowska-Mazur (2003), Rädler et al. (1980), Rogachevskii & Kleeorin (2001) and Kleeorin et al. (2003). Subsequently, we propose a new methods based on multidimensional regression (MR) fitting and the Levenberg-Marquardt (LM) method.

The plan of the paper is as follows. In Sect. 2 we present the basic equations. In Sect. 3 we summarize the standard methods of computations of the dynamo coefficients. In Sect. 4 we propose modifications of standard methods as well as new methods of computing the dynamo coefficients. In Sect. 5 we present numerical results. We compare different methods in Sect. 6, discuss our results in Sect. 7, and conclude in Sect. 8.

   
2 Basic equations, the dynamo coefficients and numerical models

2.1 Equations

The induction equation is a starting point for developing dynamo theory. We extend our consideration to the non-ideal case:

 \begin{displaymath}\frac{\partial {\vec{B}}}{\partial {t}} = \nabla \times \left( \vec{V} \times \vec{B} - \eta \nabla \times \vec{B} \right) , \end{displaymath} (1)

where $\vec{B}$ and $\vec{V}$ are the magnetic and velocity fields, respectively, and $\eta$ is the Ohmic diffusion coefficient.

Following Moffat (1978), we separate the velocity field and the magnetic field into mean and fluctuating parts:

 
$\displaystyle \vec{V} = \overline{\vec{V}} + \vec{v},$    $\textstyle \overline{\vec{v}} = 0 ,$   (2)
$\displaystyle \vec{B} = \overline{\vec{B}} + \vec{b},$    $\textstyle \overline{\vec{b}} = 0 .$    

Substituting the decomposed velocity and magnetic field into Eq. (1) and performing algebraic transformations, we obtain the equation for mean magnetic field evolution:

 \begin{displaymath}\frac{\partial {\overline{\vec{B}}}}{\partial {t}} = \nabla \... ...ine{\vec{B}} - \eta \nabla \times \overline{\vec{B}} \right) , \end{displaymath} (3)

where, apart from the terms dependent only on large scale fields, there is also a term describing the effect of small scale fluctuations on large scale electromotive force $\vec{{\cal E}}$. Rädler (1980) suggested representing $\vec{{\cal E}}$ as caused by fluctuating motions:

 \begin{displaymath}\vec{{\cal E}} = \overline{\vec{v} \times \vec{b}} = - \alpha... ...elta \times \left( \nabla \times \overline{\vec{B}} \right) , \end{displaymath} (4)

with the following physical interpretation of the terms: In the above equation we ignored other terms that are proportional to higher derivatives of magnetic field components with respect to spatial coordinates. The linear approximation of the electromotive force (EMF) can be written in the following form (Rädler 1980):

 \begin{displaymath}{{\cal E}}_{i} = a_{ij} \overline{B}_j + b_{ijk} \frac{\parti... ...partial {x_k}} , \qquad {\rm with} \quad i,j,k \in \{x,y,z\} . \end{displaymath} (5)

Where we assume that EMF is determined by $\overline{B}$ and its first derivative (see Rädler 1980). The relation between tensor elements $\alpha_{ij}, \beta_{ij}, \gamma_i, \delta_i$ and aij, bijk is
 
    $\displaystyle \alpha_{ij} = - \frac{1}{2} \left( a_{ij} + a_{ji} \right) , \qquad \gamma_{i} = \frac{1}{2} \epsilon_{ijk} a_{jk} ,$  
    $\displaystyle \beta_{ij} = \frac{1}{4} \left( \epsilon_{ikl} b_{jkl} + \epsilon... ... \right) , \qquad \delta_{i} = - \frac{1}{4} \left( b_{kik} - b_{kki} \right) .$ (6)

We try to derive coefficients aij and bijk from the EMF calculated at each time step in our 3D cosmic-ray driven dynamo numerical experiment (see Sect. 2.2).

In the rest of paper we denote the tensor coefficients calculated directly from the fitting method as aij, bijk (see Eq. (5)). The new coefficients obtained after the tensor operation described by Eq. (6) are denoted as Greek letters, e.g. $\alpha_{ij}$ and  $\beta_{ijk}$.

   
2.2 Numerical model of the cosmic-ray driven dynamo

The calculations of $\alpha $ and $\beta $ involve the velocity field and the magnetic field, gaseous density and EMF obtained in our previous simulations of the cosmic-ray driven dynamo (Hanasz et al. 2004). In that model, we have included the following physical elements: the cosmic ray component described by the diffusion-advection transport equation (see Hanasz & Lesch 2003 for the details of numerical algorithm), cosmic rays diffusing anisotropically along magnetic field lines (Giacalone & Jokipii 1999; and Jokipii 1999), supernova remnants exploding randomly in the disk volume (see Hanasz & Lesch 2000), the localized resistivity of the ISM (see Hanasz et al. 2002; Hanasz & Lesch 2003; Tanuma et al. 2003), and the realistic vertical disk gravity (Ferrière 1998). The sizes of the computational volume are: 500 pc $\times $ 1000 pc $\times $ 1200 pc in X, Y and Z directions extending from z=-600 pc to z=+600 pc with a resolution of $50 \times 100 \times 120$ grid points (see Hanasz et al. 2004). The boundary conditions are periodic in the Y direction, sheared in the X direction (following Hawley et al. 1995) and open in the Z direction. The system of coordinates x,y,z corresponds locally to the global galactic cylindrical system $r, \phi, z$. The disk rotation was defined by the values of the angular velocity $\Omega=0.05{\rm ~ Myr}^{-1}$ and the value of the shearing parameter q=1. In our model (Hanasz et al. 2004) the supernovae explode with the frequency $2 {\rm ~ kpc}^{-2}{\rm ~ Myr}^{-1}$. We assumed that 10% of the $10^{51} {\rm ~ erg}$ kinetic energy output from SN is converted into the cosmic ray energy. The cosmic ray energy is injected instantaneously into the ISM with a Gaussian radial profile ( $r_{\rm SN}=50{\rm ~ pc}$) around the explosion center. The applied value of the CR parallel diffusion coefficient was $K_\parallel = 10^4 {\rm ~ pc}^2{\rm ~ Myr}^{-1}=3 \times 10^{27} {\rm ~ cm}^2~{\rm s}^{-1}$ (i.e. 10% of the realistic value) and the perpendicular one was $K_\perp =10^3 {\rm ~ pc}^2{\rm ~ Myr}^{-1} = 3 \times 10^{26} {\rm ~ cm}^2~ {\rm s}^{-1}$.

   
3 Selected standard methods of determining dynamo coefficients

3.1 Local method

Brandenburg & Sokoloff (2002, hereafter BS02) obtained the dynamo coefficients in tensorial form by solving the following system of equations:

 
                              $\displaystyle \emf_x$ = $\displaystyle a_{xx} \overline{B}_{x} + a_{xy} \overline{B}_{y} + b_{xxz} \overline{B}'_{x} + b_{xyz} \overline{B}'_{y} ,$  
$\displaystyle \emf_y$ = $\displaystyle a_{yx} \overline{B}_{x} + a_{yy} \overline{B}_{y} + b_{yxz} \overline{B}'_{x} + b_{yyz} \overline{B}'_{y} .$ (7)

This system can also be written in the form:

 \begin{displaymath}\vec{E}^{i}(z) = \vec{M}(z) ~ \vec{C}^{i} , \qquad i = x, y \end{displaymath} (8)

where:
                                   $\displaystyle \vec{M}$ = $\displaystyle \left( \begin{array}{cccc} \langle {\overline{B}_x ~ \overline{B}... ...angle {\overline{B}_y' ~ \overline{B}_y'} \rangle \nonumber \end{array}\right),$  
$\displaystyle \vec{E}^i$ = $\displaystyle \left( \begin{array}{c} \langle {{\cal E}_i \overline{B}_x } \ran... ... \begin{array}{c} a_{ix} \\ a_{iy} \\ b_{ixz} \\ b_{iyz} \end{array}\right).$ (9)

The averages $\langle {\ldots} \rangle$ are taken over time. In the above method, only four terms are used. Firstly, due to XY plane averaging, all terms with derivatives over x and y vanish, so in the above equations prime (') denotes a derivative over z. Secondly, in Brandenburg & Sokoloff simulations described in BS02, Bz was zero initially and did not change during the entire evolution. The same method was used by Ziegler et al. (1996), but they calculated only the a-coefficients.

3.2 Non-local method

In order to obtain non-local dynamo coefficients, we use the method applied in BS02. The authors define the EMF as a convolution of dynamo coefficients and magnetic field components

 
                              $\displaystyle \emf_x$ = $\displaystyle a_{xx} \ast B_{x} + a_{xy} * B_{y} + b_{xxz} * B_{x}' + b_{xyz} * B_{y}' ,$ (10)
$\displaystyle \emf_y$ = ayx * Bx + ayy * By + byxz * Bx' + byyz * By' .  

In the above equations, aij(z, z', t) and bijz(z, z', t) are integral kernels and an asterisk (*) denotes convolution. The convolution is written in the integral form as
                 $\displaystyle \left[a_{ij} * B_{j}\right] (z,t)$ $\textstyle \equiv$ $\displaystyle \int_0^{L_z}{a_{ij}(z,z',t) B_{j}(z',t) {\rm d}z'} ,$ (11)
$\displaystyle \left[b_{ijz} * B_{j}'\right] (z,t)$ $\textstyle \equiv$ $\displaystyle \int_0^{L_z}{b_{ijz}(z,z',t) \frac{\partial B_{j}(z',t)}{\partial z'} {\rm d}z'} ,$  

where the integration is taken over the vertical size Lz of domain. They find the Fourier coefficients of the electromotive force z-derivative ( $\hat{{\cal E}}'_{x}$ and $\hat{{\cal E}}'_{y}$) with the aid of a Fourier transform
 
$\displaystyle \hat{{\cal E}}_x' = k \hat{a}_{xx} \tilde{B}_x + k \hat{a}_{xy} \tilde{B}_y - k^2 \hat{b}_{xxz} \hat{B}_x - k^2 \hat{b}_{xyz} \hat{B}_y ,$     (12)
$\displaystyle \hat{{\cal E}}_y' = k \hat{a}_{yx} \tilde{B}_x + k \hat{a}_{yy} \tilde{B}_y - k^2 \hat{b}_{yxz} \hat{B}_x - k^2 \hat{b}_{yyz} \hat{B}_y ,$      

where

 \begin{displaymath}k = k_{\rm n} = (n + 1/2) \pi / L_{z} . \end{displaymath} (13)

They solve Eq. (12) with respect to the components of the $\hat{a}$ and $\hat{b}$ tensors in the same way as the one applied for Eq. (8). The hat ($\hat{~}$) in the above equations denotes cosine transformation coefficients and the tilde ($\tilde{~}$) denotes sine transformation coefficients. They assume integral kernels of the form
 
$\displaystyle a_{ij} (z,z') = (2/L_z) \sin k_0 z \sum_k{\hat{a}_{ij}(k) \sin kz \sin kz'} ,$     (14)
$\displaystyle b_{ijz} (z,z') = (2/L_z) \sum_k{ \hat{b}_{ijz}(k) \sin kz \sin kz'} .$      

The factor $\sin k_0 z$ ensures antisymmetry of the $\alpha $-effect with respect to the equatorial plane. In the rest of paper, the hat ($\hat{~}$) symbol refers to Fourier coefficients.

3.3 Kinematic methods

Another method of computing dynamo coefficients is based on the linear approximation of the $\vec{a}$-tensor, which is calculated from the small-scale velocity and magnetic field in the manner of Otmianowska-Mazur (2003). In the present paper, we compute additionally the magnetic part of the $\vec{a}$-tensor. A conservation of the total magnetic helicity in objects possessing high magnetic Reynolds numbers (e.g. in galaxies) determines the magnetic part of the $\vec{a}$-tensor (e.g. Pouquet et al. 1976 and Kleeorin & Rogachevskii 1999):

 \begin{displaymath} \vec{a}= \vec{a}^{V} - \vec{a}^{B}, \end{displaymath} (15)

where V and B denote the kinematic and magnetic parts of $\vec{a}$, respectively. We calculate the nonuniform form of these tensors following Moffat (1978) and Otmianowska-Mazur (2003):
 
$\displaystyle {a}^{V}_{ij}(\vec{x}, t) = \epsilon_{\rm ilm} \int_0^\infty{\over... ...}, t) \frac{\partial v_{\rm m}(\vec{x}, \tau)}{\partial x_{\rm j}}} \rm d\tau},$     (16)

where $\vec{v}$ is the small-scale velocity of gas and $\vec{x}$ is the position vector. Averaging here has the same meaning as in Sect. 2. We introduce a similar form to describe the magnetic part $\vec{a}^{B}$:
 
$\displaystyle {a}^{B}_{ij}(\vec{x}, t) = \epsilon_{\rm ilm} \int_0^\infty{\over... ... b_{\rm m}(\vec{x}, \tau)}{\partial x_{\rm j}}}{\rm d}\tau / \overline{\rho}} ,$     (17)

where $\vec{b}$ is the small-scale magnetic field and $\overline{\rho}$ is the averaged density of gas. We integrate back in time all 9 components of the magnetic and kinematic tensors and compute the full a-tensor. The dynamo coefficients obtained by the kinematic method are denoted as "aij'' in the rest of paper.

   
4 Proposals for modifications of the standard methods

In order to derive the dynamo coefficients from the electromotive forces obtained in our numerical simulations of the cosmic-ray driven dynamo (Hanasz et al. 2004) we use the statistical fitting methods, which we present below. The direct motivation for the modification of standard methods comes from the fact that none of the methods described above enables a satisfactory reconstruction of the electromotive forces and the mean magnetic field on the base of the determined dynamo coefficients. This will be demonstrated in Sect. 6.

The modifications of standard methods include adding a new constant term to the formula for electromotive forces of BS02 (Eq. (7)), implementation of the multi-linear regression method (MR), with and without the constant term, and other minor changes such as implementation of a running window in the computation of the components of dynamo tensors. The major aim of these modifications is to search for an accurate parameterization of the dynamo cosmic-ray driven dynamo process. As we shall see, one of improvements in the EMF reconstruction comes from the implementation of the multi-linear regression method and another from the addition of the constant term. The latter modification attempts to include, with a single number, all the truncated higher order terms in the Taylor expansion of electromotive forces in Eq. (4). Although we cannot provide a direct physical interpretation of the new term, we attempt to implement a less restrictive ansatz for the electromotive forces and therefore to search for a possible direction for better estimations of the dynamo coefficients in future.

Among the minor changes, we examine different frames in the EMF fitting procedures: a running fitting window on the time axis, with 100 Myr as the characteristic time in our model (the correlation time for magnetic field and velocity components is of the order of 25-50 Myr), and single fitting in time intervals longer than 100 Myr. The running frame is used to obtain the best fit to the time evolution of EMF. We apply the long interval fitting frames in order to get only three mean values of the dynamo coefficients in the whole time period, which fit sufficiently well the EMF calculated on the basis of simulation data. Such values subsequently can be used in large-scale dynamo simulations. We chose the following three long time-intervals: 0-1370 Myr, 1370-1820 Myr and 1820-2300 Myr. In our calculations of the dynamo coefficients, we use two types of averaging: plane averaging and cube averaging. In the case of plane averaging, the magnetic and velocity fields are averaged over constant-z planes. As a result, all components of a and b tensors containing x- and y-derivatives vanish and the related components of the $\alpha $ and $\beta $ tensors also vanish. The second type of averaging is performed in rectangular boxes of 200 pc in each direction, which is comparable to or larger than the correlation length of cosmic-ray driven turbulence that was computed in our simulations. In this case, all components of $\alpha $ and $\beta $ tensors can be determined, although some components vanish because the mean vertical magnetic field component is conserved and equal to zero in the initial state.

4.1 Modified local method

In the modified Brandenburg-Sokoloff method, in addition to the standard BS02 method, we also take into account constant terms in Eq. (7),

 
$\displaystyle \emf_x = a_{xx} \overline{B}_{x} + a_{xy} \overline{B}_{y} + b_{xxz} \overline{B}'_{x} + b_{xyz} \overline{B}'_{y} + \emf_x^0 ,$      
$\displaystyle \emf_y = a_{yx} \overline{B}_{x} + a_{yy} \overline{B}_{y} + b_{yxz} \overline{B}'_{x} + b_{yyz} \overline{B}'_{y} + \emf_y^0 .$     (18)

This implies the following modification of the $\vec{M}$ and $\vec{C}$ matrices
$\displaystyle \vec{M} = \left( \begin{array}{ccccc} \langle {\overline{B}_x ~ \... ... \rangle & \langle { \overline{B}_y'} \rangle & 1 \nonumber \end{array}\right),$      
$\displaystyle \vec{E}^i = \left( \begin{array}{c} \langle {{\cal E}_i \overline... ...c} a_{ix} \\ a_{iy} \\ b_{ixz} \\ b_{iyz} \\ \emf_{i}^0 \end{array}\right).$     (19)

As we see, this extension is rather simple and can be solved in the same way as the standard BS02 method.

4.2 Multidimensional regression (MR)

We assume the linear relation of the quantity $\tens{Y}$ with k independent variables $\tens{X}_{k}$:

 
$\displaystyle \tens{Y} = \tens{a}_1 \tens{X}_1 + ... + \tens{a}_{k} \tens{X}_{k} + \tens{a}_{k+1} .$     (20)

In the present case $\tens{Y}$ represents $\emf_i$ and $\tens{X}_{k}$ contains all $\overline{\vec{B}}_{j}$ and $\partial_{x_k}\overline{\vec{B}}_{j}$ expressions. Note that we also include a constant term $\tens{a}_{k+1}$, which is independent of any $\vec{X}_i$. For n different and independent samples at each time step we have:
    $\displaystyle \tens{Y} = \left[ \begin{array}{c} \tens{y}_1 \\ \tens{y}_2 \\ ... ...{x}_{1n} & \tens{x}_{2n} & \ldots & \tens{x}_{kn} & 1 \\ \end{array} \right] ,$ (21)
    $\displaystyle \tens{A} = \left[ \begin{array}{c} \tens{a}_1 \\ \tens{a}_2 \\ \ldots \\ \tens{a}_{k+1} \end{array} \right] .$  

Our system of equations can be rewritten in matrix form:
 
$\displaystyle \tens{Y} = \tens{X} \tens{A} .$     (22)

The fitting is made with the least square method. As a result, we obtain the following matrix solution:
 
$\displaystyle \tens{A} = \left( \tens{X}^{\rm T} \tens{X} \right)^{-1} \tens{X}^{\rm T} \tens{Y}.$     (23)

Since we use simulation data, the standard deviation of measurements are unknown. In this case we use a derived equation to compute the uncertainties of fitted coefficients

 \begin{displaymath}\sigma_{\tens{A}}^2 = {\rm diag} \left( \tens{X}^{\rm T} \tens{X} \right)^{-1} \chi^2 / (n - k) , \end{displaymath} (24)

where the quantity $\chi {^2}$ measures the quality of fitting and is defined in the next subsection by Eq. (25) and (n-k) is the number of degrees of freedom. For a detailed discussion of the computation $\sigma^2$-coefficients, see Numerical Recipes (Press et al. 1997, Sect. 15). For reliable comparison when fitting to the same function but with a different number of terms or points, we divide $\chi {^2}$ by the number of points n, which results in a normalized quality of fitting value.

4.3 Levenberg-Marquardt method (LM)

A linear relationship can be fit by minimizing the $\chi {^2}$ residuals, which are defined as follows:

 \begin{displaymath}\chi^2 \left( \tens{X}_{\rm i}, \tens{a} \right) = \sum_{n}{ ... ...[ \tens{Y}_{\rm i} - \tens{a} ~ \tens{X}_{\rm i} \right] }^2 , \end{displaymath} (25)

Levenberg proposed an efficient method to find the minimum of the $\chi {^2}$ function using a gradient method far from the minimum and an expansion method in its vicinity. The method searches for the coefficients iteratively using formula

\begin{displaymath}\tens{a}_{\rm k+1} = \tens{a}_{\rm k} - \left[ \frac{1}{2} \nabla \chi^2 ( \tens{a}_{\rm k} ) \right] \tens{\kappa}^{\rm -1} , \end{displaymath} (26)

where elements of the $\tens{\kappa}$ tensor are defined as $\kappa_{kl}=\frac{1}{2} \frac{\partial^2 \chi^2}{\partial a_k \partial a_l}$. Minimization is performed using the MPFIT package for IDL (Markwardt IDL library). Both methods give us values of $\sigma$, $\chi {^2}$ and the other statistical parameters.

   
5 Computation of the dynamo coefficients with the aid of statistical fitting methods

5.1 The local formulation


  \begin{figure} \par\includegraphics[width=8.2cm,clip]{3582f1a.eps}\hspace*{3mm} \includegraphics[width=8.1cm,clip]{3582f1b.eps} \end{figure} Figure 1: The time evolution of $\alpha _{xx}$ ( left) and $\alpha _{yy}$ ( right) obtained with the running frame MR method with CP (the solid line) and without CP (the dash-dotted line) at two chosen cube heights.
Open with DEXTER


  \begin{figure} \par\includegraphics[width=8.2cm,clip]{3582f2a.eps}\hspace*{3mm} \includegraphics[width=8.4cm,clip]{3582f2b.eps} \end{figure} Figure 2: The time evolution of $\alpha _{zz}$ ( left) and $\beta _{xx}$ ( right) obtained with the running frame MR method with CP (the solid line) and without CP (the dash-dotted line) at two chosen cube heights.
Open with DEXTER


  \begin{figure} \par\includegraphics[width=8.4cm,clip]{3582f3a.eps}\hspace*{3mm} \includegraphics[width=8.4cm,clip]{3582f3b.eps} \end{figure} Figure 3: The time evolution of $\beta _{yy}$ ( left) and $\beta _{zz}$ ( right) obtained with the running frame MR method with CP (the solid line) and without CP (the dash-dotted line) at two chosen cube heights.
Open with DEXTER


  \begin{figure} \par\includegraphics[width=7.7cm,clip]{3582f4.eps} \end{figure} Figure 4: The comparison of the $\emf_y$ with the fitted curves obtained with MR methods with the running frame of 100 Myr calculated at the height of 400 pc as a function of time.
Open with DEXTER

In this subsection we present the results obtained with the volume averaging MR method. The volume averaging procedure is applied prior to all other steps of calculation in order to obtain local values of EMF in the given cube regions. The global LM method gives identical results to the MR fitting, so here we only show output from the latter. The time evolution of $\alpha $ and $\beta $ coefficients was calculated with the aid of the MR method with (solid line) and without (dash-dotted line) the constant part of the electromotive force $\vec{{\cal E}}$, at two chosen heights of the cube (see Figs. 1-3). We use the abbreviation CP for the "constant part'' term in the next sections. The graphs present the curves obtained by the running fitting procedure (with a frame of 100 Myr), which is applied to the volume-averaged quantities. This method gives the best quality of fit (see Fig. 4) and provides the dynamo coefficients as a function of time. The second method of estimation, mentioned in the beginning of Sect. 4, is the MR fit applied to the volume averaged EMF, in three long time periods: 0-1370 Myr, 1370-1820 Myr and 1820-2300 Myr. This calculation gives a much smaller number of dynamo coefficients, which are presented in Tables 1-3 (with CP) and Table 4 (without CP). Next, we compare their values with the mean values of the same coefficients obtained from the running frame (100 Myr) MR fitting method with CP (Table 5) averaged in the same long period. We will demonstrate that such mean values of the dynamo tensors do not properly fit the EMF (Fig. 5). In the tables we also present the statistical parameter $\chi {^2}$ and it is possible to compare its values when we fit the same EMF at the same cube height level.

Table 1: Fitting of the whole range with the constant term: coefficients for T = 0-1370 Myr. Units of z, $\alpha $ and $\beta $ coefficients are [pc], [km s-1] and [1026 cm2/s], respectively.

Table 2: Fitting of the whole range with the constant term: coefficients for T = 1370-1820 Myr. Units of z, $\alpha $ and $\beta $ coefficients are [pc], [km s-1] and [1026 cm2/s], respectively.

Figure 1 presents the time evolution of $\alpha _{xx}$ (left) and $\alpha _{yy}$ (right) derived with the MR method with and without CP at a height of $\pm$200 pc. During the entire evolution time, $\alpha _{xx}$ changes between +70 km s-1 and -30 km s-1. The coefficient  $\alpha _{yy}$ takes values that are one order of magnitude smaller than $\alpha _{xx}$. Its value fluctuates between -3 km s-1 (at a height of -200 pc) and +4 km s-1 at both heights. The calculations with CP (the solid line) and without it (the dash-dotted line) give slightly different values of $\alpha _{xx}$ and $\alpha _{yy}$.

Table 3: Fitting of the whole range with the constant term: coefficients for T = 1820-2300 Myr. Units of z, $\alpha $ and $\beta $ coefficients are [pc], [km s-1] and [1026 cm2/s], respectively.

Table 4: Fitting of the whole range without the constant term: coefficients for T = 1370-1820 Myr. Units of z, $\alpha $ and $\beta $ coefficients are [pc], [km s-1] and [1026 cm2/s], respectively.

Table 5: The mean values of the coefficients obtained with the aid of running fitting MR with CP averaged in the time period T = 1370-1820 Myr. Units of z, $\alpha $ and $\beta $ coefficients are [pc], [km s-1] and [1026 cm2/s], respectively.

In the left panel of Fig. 2, we show the time variation of the coefficient $\alpha _{zz}$ at two chosen cube heights ($\pm$200 pc) for the MR with CP (the solid line) and without CP (the dash-dotted line), fluctuates between +11 km s-1 and -8 km s-1. The next graphs (right panel of Figs. 2 and 3) show that all $\beta $ coefficients oscillate between positive and negative values. The Rädler (1980) prescription of the dynamo tensor $\beta $ applied in Eq. (4) is constructed in such way that this symmetric part of the diffusion tensor should be positive. The problem of the resulting negative values will be discussed later in the Discussion (Sect. 7). The coefficient $\beta _{xx}$ changes its value from $0.8\times 10^{26}$ to $-1.2\times 10^{26}{\rm ~ cm}^2/{\rm s}$ at $z=-200{\rm ~ pc}$ level, being slightly lower at +200 pc height. Figure 3 shows the time evolution of the last two diffusion coefficients $\beta _{yy}$ (left panel) and $\beta _{zz}$ (right panel) calculated by means of the MR method without (the dash-dotted line) and with (the solid line) CP resulting in quite similar values for both coefficients. The value of $\beta _{yy}$ changes during the entire time period, reaching $\pm0.8\times 10^{26}{\rm ~ cm}^2/{\rm s}$ at both chosen heights. The changes of $\beta _{zz}$ over time (Fig. 3, right panel) give slightly higher values as $\beta _{yy}$, reaching an absolute maximum of about $1.6\times 10^{26}{\rm ~ cm}^2/{\rm s}$ at both heights.

  \begin{figure} \par\includegraphics[width=7.6cm,clip]{3582f5a.eps}\hspace*{3mm} \includegraphics[width=7.6cm,clip]{3582f5b.eps} \end{figure} Figure 5: The time evolution of $\emf_y$ at the height of 200 pc ( left) and 400 pc ( right) shown in three time intervals in which the MRC method is fitted (the black line). The dashed line shows the EMF reconstruction with the help of the averaged coefficients in the given time periods (see text for the explanation).
Open with DEXTER

The time evolution of $\emf_y$ (the gray line) calculated with the volume averaging procedure at a height of 400 pc above the disk plane is shown in Fig. 4. The curve is overlaid with the fitting curve (the solid line) obtained with the running frame MR method with CP (as shown in Figs. 1-3). We also calculated the quality of fit over the entire time to compare our data with other methods. The obtained value of $\chi^2_{\rm y}=0.277\times 10^{-6}$ shows that the quality of fit is high, due to the small value of $\chi^2_{\rm y}$ in comparison with the other methods of estimation (see Tables 1-4). In order to compare this approach with the long time frame procedure we also calculate $\chi^2_y$ for the case of the running frame method in three periods: 0-1370 Myr, 1370-1820 Myr, 1820-2300 Myr. Their values are respectively equal to $0.616\times 10^{-6}$, $59.8\times 10^{-6}$, and $74.8\times 10^{-6}$. The first value is one order of magnitude smaller, the second is 1.5 times lower and the third one is again one order of magnitude smaller than values of $\chi^2_y$ at the 400 pc height shown in Tables 1-3, respectively. This comparison shows that fitting the EMF with the running frame method has a higher quality than estimation with the long time periods.

As a next step, we present the results obtained with the long time frame fit of the MR method with and without CP. The values of $\alpha $ and $\beta $ fitted by the first procedure in the frame: 0-1370 Myr, 1370-1820 Myr and 1820-2300 Myr are shown in Tables 1-3 at the four chosen cube heights. The second one is shown in Table 4 for one time period 1370-1820 Myr. We can see that $\alpha _{xx}$ changes its value from -1.5 to 0.0 km s-1 in the first period (Table 1), -2.90 to 0.28 km s-1 in the second period (Table 2) and from -3.25 to 0.54 km s-1 in the last period (Table 3). The coefficient $\alpha _{yy}$ is one order of magnitude lower than the first coefficient with its absolute maximum around 0.19 km s-1 in all time frames. The value of $\alpha _{zz}$ ranges from 0.18 km s-1 in the first period (Table 1) to the maximum absolute value 1.61 km s-1 at z=-200 pc in the second period (Table 2), being slightly lower in the last period (Table 3). The values of all three $\alpha $ tensors obtained with the MR method without the CP, presented in Table 4 (period 1370-1820 Myr), are similar to the values resulting from MR with the CP procedure visible in Table 2.

The first diffusion coefficient $\beta _{xx}$ changes its value from -0.02 to $0.0\times 10^{26}{\rm ~ cm}^2/{\rm s}$ in the first time interval (Table 1). In the next chosen period it ranges from -4.9 to - $0.1\times 10^{26}{\rm ~ cm}^2/{\rm s}$ (Table 2), and so is one order of magnitude higher than before. In the last period this coefficient ranges between -0.05 and $-0.03\times 10^{26}{\rm ~ cm}^2/{\rm s}$ (Table 3). The coefficient $\beta _{yy}$ has its maximum absolute value of about $0.05\times 10^{26}{\rm ~ cm}^2/{\rm s}$ during first two periods, but it decreases to $0.01\times 10^{26}{\rm ~ cm}^2/{\rm s}$ in the last period (1820-2300 Myr). The last diffusion coefficient $\beta _{zz}$ has the highest values ranging from an absolute maximum of about $0.14\times 10^{25}{\rm ~ cm}^2/{\rm s}$ in the first period (Table 1) and $0.05\times 10^{26}{\rm ~ cm}^2/{\rm s}$ in the second one (Table 2), growing to about $0.33\times 10^{26}{\rm ~ cm}^2/{\rm s}$ in the last period (Table 3). The MR method without the CP gives slightly different values for all $\beta $ tensors but the order of magnitude of these quantities is similar (see Table 4). Again we obtain negative values for the diffusion coefficients in all long time intervals (see Discussion). The values of $\chi {^2}$ computed for the three electromotive forces along the X, Y and Z axe at certain levels can be compared with the similar quantities visible in Table 2. We can see that almost all values in the first table are only slightly higher than in the second one. This means that we can use neither method, i.e. MR with and without CP, to obtain similar quality. Although the $\chi {^2}$ values obtained from running frame estimations are smaller than quantities resulting from long period frame simulations, the estimated curves obtained with the help of the last method fit quite well to the original EMF. In Fig. 5 we present one of them, $\emf_y$ (gray line) at two levels 200 pc and 400 pc above the disk plane, as examples. We can see that the solid line describes the $\emf_y$ at both levels quite well.

The values of the dynamo tensors computed with the help of the running frame methods have the highest quality of fit, but they change continuously over time. In order to use the obtained dynamo coefficient in the large-scale dynamo calculations of a galactic disk we need fewer of them. Normally, this is done by calculating the mean over the entire time time period. We apply such averaging in three frames: 0-1370 Myr, 1370-1820 Myr, and 1820-2300 Myr. The values are presented in Table 5. The fitted curves of $\emf_y$ obtained from such averaged dynamo tensors are also presented in Fig. 5 as the dashed lines. It is evident that the curves do not fit the original $\emf_y$ line at all, which is also shown by the very high values of $\chi {^2}$ (see Table 5). This means that the simple averaging of the estimated dynamo tensors give values (Table 5) of too poor quality to be used.

   
6 Comparison of different methods of computation of the dynamo coefficients

6.1 The quality of fitting measured by $\chi {^2}$ for MR, LM and BS02

In order to compare the different methods of coefficient estimation from electromotive forces, we apply the methods used in BS02 and Ziegler et al. (1996) to our model of the cosmic-ray driven dynamo (Hanasz et al. 2004). Now, we use the averaging of the EMF over planes $z={\rm constant}$ prior the other steps of calculation (see Otmianowska-Mazur 2003). To show the role of the constant part, we apply the MR method with and without CP. We also include a constant part of $\vec{{\cal E}}$ in the BS02, resulting in the modified BS02 method. We use here the long frame fitting method in the same time periods as before: 0-1370 Myr, 1370-1820 Myr and 1820-2300 Myr. Due to different methods of averaging, the time evolution of all EMF is also different.

  \begin{figure} \par\includegraphics[width=7.7cm,clip]{3582f6.eps} \end{figure} Figure 6: The comparison between different methods of estimation at z=200 pc. Dotted vertical lines designate boundaries of the single fitting regions. Fitting is performed in 3 time periods: 0-1370, 1370-1820, 1820-2300 Myr.
Open with DEXTER

Figure 6 presents the time evolution of the electromotive force in the Y direction ($\emf_y$, the gray line) over the whole time period. Because the MR method of estimation gives a solution identical to the LM method we show only this curve (as the dash-dotted line) in the figure. The solid line marks the same method but with CP. We can see that this curve fits the best to the actual $\emf_y$. The dots mark the fitting curve obtained with the BS02 method. The dashed line shows the same procedure but with CP taken into account. The mean values of the dynamo coefficients (fitted in the given time period) and the quality of fitting of all methods are presented in Table 6 for the time interval, 1370-1820 Myr. The value of $\chi {^2}$ in the case of MR and MR plus CP is nearly half the value in the case of BS02. The modified BS02 possesses a slightly smaller $\chi {^2}$ value than the standard BS02, but still it is higher than quantities resulting from both MR methods. For comparison we show also the values of a and b computed at a cube height of 400 pc. We can see that both MR methods gives very similar values for all coefficients. The value of CP is also lower for the case of MR plus CP than for the case of BS02 plus CP.

Neither BS02 methods, the regular (BS) and modified one (BSC), include terms with ayz and byzz coefficients. Both MR methods with these terms included give relatively high absolute values (about -3 km s-1 in case of ayz and about -0.2 cm2/s for byyz). However, the z component of the magnetic field is at least two orders of magnitude lower than the y component. This gives a relatively small contribution to the total electromotive force.

Table 6: Coefficients and $\chi {^2}$ obtained using different methods for time period T = 1370-1820 Myr. Units of a, b and ${\cal E}^0$ are [km s-1], [1026 cm2/s] and [km s-1 $\mu $G], respectively.

6.2 Comparison of reconstructed electromotive forces

We check the quality of our estimation methods by integrating the induction equation (Eq. (3)) rewritten in a simplified form

 \begin{displaymath}\frac{\partial \overline{B}_{x}}{\partial t} = - \frac{\parti... ...\partial \emf_{\rm x}}{\partial z} - q \Omega \overline{B}_{z} \end{displaymath} (27)

using the fitted electromotive forces (see Fig. 7). In this integration we use plane averaging, so all x and y derivatives vanish. Following Brandenburg & Sokoloff (2002) we assume that, due to the periodic boundary conditions in the horizontal directions, the mean magnetic field component $\overline{B}_{z}$ does not change during evolution of the dynamical system. This assumption implies $\overline{B}_{z}=0$ during the whole simulation, if the z component of the magnetic field is zero initially. Therefore we focus only on the evolution of the x and y components of the mean magnetic field. The last term in the second equation in Eq. (27) corresponds to the local velocity shear. We obtain the magnetic field components $\overline{B}_{x}$ (Fig. 7, left panel) and $\overline{B}_{y}$ (Fig. 7, right panel). Both graphs in the figure compare the two method of estimation applied together with the plane averaging: BS02 and MR, both with a running frame of 100 Myr. With the solid line we mark the components obtained in simulations, with the dotted line we mark the BS02 method, and with the dashed line we mark the regression method (MR). The quality of fitting is confirmed by comparison to the third dash dot line (see Fig. 7), which represents direct integration of the original electromotive force obtained from the simulation data. We see that the fitting process does not cause any significant lost of information and the character of ${\cal E}$ is preserved. However, comparison to the mean magnetic field components $\overline{B}_{x}$ and $\overline{B}_{y}$ obtained directly from the simulation data shows that the set of Eqs. (27) is too simple to describe the complicated processes included in the fast dynamo model. Neither do they include magnetic diffusion, which is presented in our simulations. These reasons could explain the systematic shifts between integrated and directly computed magnetic field components beyond a time of about 1000 Myr.

6.3 Non-local formulation according to the BS02 method

We apply the Fourier transform to obtain the non-local values of $\hat{a}$ and $\hat{b}$ prescribed in the form of Eq. (14) in our experiment. We use plane averaging, because in cube averaging we would have few points to make the Fourier analysis reasonable. In cube averaging we separate individual cubes in order to preserve the independence of the averages. The correlation length is estimated as about 200 pc, so separation was computed with these averages. It results in a spacial resolution reduced by the number of points in each direction included in an individual cube.

In Fig. 8 we show graphs presenting all coefficients of $\hat{a}$ (upper row) and $\hat{b}$ (lower row) as functions of the harmonic number n. The $\hat{a}_{xx}$ coefficient starts from a value of about -20 km s-1 for n=0. This means that the mean global value of $\hat{a}_{\rm xx}$ is of the same order. This value is much higher than the values of $\alpha _{xx}$ resulting from the MR methods with the long time frame. For higher n it oscillates around zero with a decaying amplitude of the order of a few km s-1 for the first harmonic numbers. The next graph presents the $\hat{a}_{xy}$ coefficient, which changes from -1.5 km s-1 for n=0 to smaller amplitudes for higher n. It decays slowly with the value of n. The $\hat{a}_{yx}$ coefficient has only positive values but a slightly smaller amplitudes than $\hat{a}_{xx}$. Its value decreases from 10 km s-1 for n=2 to 5 km s-1 for n=5. It then decreases to 0.0 for higher n. The coefficient $\hat{a}_{yy}$ starts from its maximum of nearly 2 km s-1 for n=0 and slowly decreases to 0.0 km s-1 for higher harmonic numbers. Again these values are one order of magnitude higher than quantities obtained with the long time frame fitting methods.

Figure 8 in the lower row shows the $\hat{b}$ coefficients as functions of n. We can see that all tensor components are positive for almost all values of n. The coefficient $\hat{b}_{xx}$ has a very high value of $24.8\times 10^{26}{\rm ~ cm}^2/{\rm s}$, which is a few times higher than the classical value of diffusion in the ISM and also two orders of magnitude higher than our values obtained from our fitting methods with the long period frames. Then $\hat{b}_{xx}$ decreases significantly reaching $0.5\times 10^{26}{\rm ~ cm}^2/{\rm s}$ for n=5. For higher n this coefficient drops almost to zero. The changes of $\hat{b}_{yx}$ start from the value of $7\times 10^{26}{\rm ~ cm}^2/{\rm s}$ for n=0 and drop quickly to $-0.5\times 10^{26}{\rm ~ cm}^2/{\rm s}$. Starting from n=5 its value is very close to zero. The last two diffusion coefficients $\hat{b}_{xy}$ and $\hat{b}_{yy}$ possess similar maxima of about $0.5\times 10^{26}{\rm ~ cm}^2/{\rm s}$ (which is close to the fitted values), then their values drop to almost zero for n>15. Although the local fitting gives us negative values of the diffusion coefficients we can see that for the non-local formulations we get mostly positive values (only one value is negative, see Discussion, Sect. 7).

  \begin{figure} \par\includegraphics[width=8cm,clip]{3582f7a.eps}\hspace*{3mm} \includegraphics[width=7.7cm,clip]{3582f7b.eps} \end{figure} Figure 7: Integration of the magnetic field components Bx and By with the help of the reconstructed electromotive force for the different methods of estimation. The solid line represents mean magnetic field components computed from the simulation data, the dotted line represents integration of the BS02 fitted ${\cal E}$, and the dashed line the MR fitted ${\cal E}$. Additionally we include integration of the actual ${\cal E}$ obtained from the simulation data (gray dash dot line).
Open with DEXTER

   
6.4 The kinematic methods

In this paper we present only the main coefficients: axx, ayy, azz to compare them with the results of computations with the MR methods of fitting. Figure 9 shows the kinematic (left) and magnetic (right) parts of the axx, tensors as a function of time. The averaging is done over the XY-planes. It is apparent that the kinematic part of axx is much larger then the magnetic component presented in the figure (right), reaching a maximum of 40 km s-1 at a height of 400 pc. It is also apparent that this component changes sign below the galactic disk, growing to similar high values. The obtained maximum of the absolute value of the magnetic part are much lower (4 km s-1) and the time evolution has chaotic behavior. The time evolution of the sum of the coefficients discussed above (Eq. (15)) is presented in Fig. 12, left, at the chosen height of 400 pc. In the figure we show three mean values of axx taken after 1370 Myr, then in the period 1370-1820 Myr, and in the last interval of 1820-2300 Myr. We see that the mean value of this coefficient grows from 4 km s-1 to 24 km s-1. Changes of the kinematic and magnetic parts of the ayy coefficient over time are visible in Fig. 10 left and right, respectively. Their values oscillate between -0.3 and +0.3 km s-1 and the mean values of their difference (Eq. (17)), presented in Fig. 12 (in the middle), give a maximum of only 0.017 km s-1. The third coefficient azz is presented in Fig. 11 with its kinematic part on the left side and with its magnetic part on the right side, again at two chosen heights in the disk. We can see that the values of both parts above the disk plane grow in time up to 15 and 12 km s-1 and to -8 and -5 km s-1 below the plane, similar to axx. Their difference, shown in Fig. 12 (right) at a height of 400 pc, gives the averaged values, which rise from 1.7 to 12 km s-1. The antisymmetric parts ayx and axy are rather small and change their value from +0.5 km s-1 to -0.5 km s-1 at random. We do not present their graphs.

The value of $\alpha _{xx}$ obtained with the kinematic estimation is one order of magnitude too high in comparison with the value obtained from the statistical methods of estimation, which fit well to the EMF. However, their qualitative behavior remains similar. The ayy coefficient is the smallest one, while axx and azz have higher mean values.


  \begin{figure} \par\includegraphics[width=5.7cm,clip]{3582f8a.eps}\hspace*{3mm} ... ...ps}\hspace*{3mm} \includegraphics[width=5.7cm,clip]{3582f8h.eps} \end{figure} Figure 8: The individual coefficients of $\hat{a}$ and $\hat{b}$ tensors as a function of harmonic number n obtained with the non-local formulation. The harmonic number n determines the value of k by Eq. (13).
Open with DEXTER


  \begin{figure} \par\includegraphics[width=6.8cm,clip]{3582f9a.eps}\hspace*{3mm} \includegraphics[width=6.8cm,clip]{3582f9b.eps} \end{figure} Figure 9: The time evolution of the kinematic ( left) and magnetic ( right) part of axx at two chosen heights of the cube.
Open with DEXTER


  \begin{figure} \par\includegraphics[width=6.8cm,clip]{3582faa.eps}\hspace*{3mm} \includegraphics[width=6.8cm,clip]{3582fab.eps}\ \end{figure} Figure 10: The time evolution of the kinematic ( left) and magnetic ( right) part of ayy at two chosen heights of the cube.
Open with DEXTER


  \begin{figure} \par\includegraphics[width=6.8cm,clip]{3582fba.eps}\hspace*{3mm} \includegraphics[width=6.8cm,clip]{3582fbb.eps} \end{figure} Figure 11: The time evolution of the kinematic ( left) and magnetic ( right) part of azz at two chosen heights of the cube.
Open with DEXTER


  \begin{figure} \par\includegraphics[width=6.5cm,clip]{3582fca.eps}\hspace*{3mm} ... ...m}\includegraphics[width=6.5cm,clip]{3582fcc.eps}\hspace*{3.5cm}} \end{figure} Figure 12: The time evolution of the coefficients axx ( left), ayy (in the middle) azz ( right) for the model at the height of 400 pc.
Open with DEXTER

   
7 Discussion

Although we do not solve all physical uncertainties connected with computing the dynamo tensors from the EMF, the presented methods give us the possibility of estimating them with the desired statistical quality. Our MR methods show that we can apply a running fitting frame, as well a long time frame (with lower quality). But it is not possible to use the mean dynamo tensor values obtained by averaging in the long time periods. The dynamo coefficients continually fluctuate over time which results from the running fitting method because the obtained curves do not properly fit the EMF.

The dynamo coefficient $\alpha _{xx}$ calculated in the long time frame fitting is between -3.3 and 0.8 km s-1 in the third time period (see Table 3), changing its sign locally. The second component $\alpha _{yy}$ is smaller, reaching mainly plus or minus 0.1 km s-1 (see Tables 1-3), and the absolute maximum value of $\alpha _{zz}$ is about 1.6 km s-1 (in the second time interval, Table 2). Although $\alpha _{yy}$ is small, there is an increase in the total magnetic energy in our cosmic-ray driven dynamo experiment. These values are much higher than in the calculations of Ziegler et al. (1996) which gave 2-6 m/s for all components of the $\alpha $ tensor. In our experiments, the maxima reached by three of the $\alpha $ components are of the order of a few or tenths of km s-1, which is only slightly lower than the values obtained by Ferrière (1993a,b, 1995, 1996, 1998) in her models without a magnetic back reaction of the gas motion.

The absolute values of $\beta $ in our model range from 0.0 to $0.33\times 10^{26}$ cm2/s. This maximum value is of the same order as the classical values of the diffusion tensor in ISM. The values of $\beta $ have negative and positive values. It is possible that the negative values of the diffusion coefficients $\beta $ obtained in our local formulation of the physically complicated model could mean that it is not always possible to work with the linear approach given in Eq. (4). The negative values of the diffusion coefficients were also present in BS02, however they finally reached positive mean and non-local values. Our analysis showed that the conditions for the linear approximation may be not strictly fulfilled and the electromotive force approximated by Eq. (4) should be extended by additional components involving small scale velocity and magnetic fields multiplied by the large scale ones ( $\overline{\vec{V}} \times \vec{b}$, $ \vec{v} \times \overline{\vec{B}}$, $\vec{v} \times \vec{b}$, see e.g. Rädler (1980) for a detailed analysis). This means that EMF could be understood as a function of $\overline{\vec{V}}$, $\vec{v}$ and $\vec{b}$, and not only $\overline{\vec{B}}$. Another issue supporting this claim is the magnitude of the constant term in fits that included CP. It can be less than or comparable to the $\alpha $ and $\beta $ terms, what mean, that higher order terms in the ${\cal E}$ expansion (see Eq. (5)) are also important. Further investigation of the assumption of a linear relation between ${\cal E}$ and components of the mean magnetic field and the physical sense of the constant term will be a subject of further study.

The non-local formulation gives much better results concerning the sign of the diffusion coefficient, in which it is mainly positive (see Fig. 8). The maximum value among $\hat{b}$ tensor coefficients has a $\hat{b}_{xx}$ of $25\times 10^{26}{\rm ~ cm}^2/{\rm s}$ for k=1. Such a high value is one order of magnitude higher than the classical diffusion coefficient in the ISM. The rest of the diffusion components are one order of magnitude smaller for the same k. The next dynamo tensor $\hat{a}$ shows that the highest values are present again for $\hat{a}_{xx}$ which gives -15 km s-1 for k=0. Similarly, the absolute value of $\hat{a}_{yx}$ is 10 km s-1, while the components $\hat{a}_{yy}$ and $\hat{a}_{xy}$ are one order of magnitude smaller. The resulting values of both coefficients obtained in the non-local formulation are one order of magnitude higher than values derived in the local approximation. This means that they could be not fit to the EMF.

Our simulations of the kinematic approximation in calculation of the dynamo tensors (Sect. 6.4) shows that the obtained values are too high in comparison with the values of a and b from the method of EMF fitting. This approximation gives similar values to those obtained from the non-local formulation.

We apply two possible approaches in the case of cube averaging: MR with and without CP, which give very similar results in reconstructing electromotive force $\vec{{\cal E}}$ from the fitting procedure (see Table 6 and Fig. 6). The values of $\chi {^2}$ are similar, so it is possible to use either of them and obtain similar quality. The BS02 procedure with CP also gives a tolerable fit to the EMF, but the quality of BS02 without CP is low (see Fig. 6).

   
8 Conclusions

Our model for estimation of the dynamo coefficients from the electromotive forces shows, for the first time, the possibility of using statistical fitting methods to solve this problem. Our study also gives a quantitative comparison between different methods of estimation.

The main results of our study can be summarized as follows:

1.
The fitted values of $\alpha $ in the long fitting frame and averaged in small cubes range from a few tenths to a few km s-1.

2.
The fitted values of $\beta $ in the long fitting frame has an amplitude of tenths of $10^{26}~{\rm cm^2/s}$, which is of the order of the classical dynamo diffusion coefficient value in the ISM.

3.
Although we did not solve all uncertainties of the theory, the presented statistical methods, LM and regression, give reliable results in estimations of the dynamo coefficients from the EMF.

4.
Methods that use the long time period averaging of the dynamo tensors, which fluctuate over time, do not give a proper fit to the EMF.

5.
The methods applied, MR with and without the constant part and LM, give us statistical measures and reliable estimation of the uncertainties.

6.
The statistical analysis shows that the new fitting methods presented in this paper give us the best approximation of the computed EMF.

7.
The kinematic and non-local approximations result in dynamo coefficients with values that were too high in comparison with those obtained by the fitting methods, but their character is similar.

Acknowledgements
This work was partly supported by the Polish Committee for Scientific Research (KBN) from the grants 1 P03D 004 26 and 1 P03D 002 28.

References

 
Copyright ESO 2006