Astrophysikalisches Institut Potsdam, An der Sternwarte 16, 14482 Potsdam, Germany

Abstract
Context. Supernovae are known to be the dominant energy source for driving turbulence in the interstellar medium. Yet, their effect on magnetic field amplification in spiral galaxies is still poorly understood. Previous analytical models, based on the evolution of isolated, non-interacting supernova remnants, predicted a dominant vertical pumping that would render dynamo action improbable.
Aims. In the present work, we address the issue of vertical transport, which is thought to be the key process that inhibits dynamo action in the galactic context. We aim to demonstrate that supernova driving is a powerful mechanism to amplify galactic magnetic fields.
Methods. We conduct direct numerical simulations in the framework of resistive magnetohydrodynamics. Our local box model of the interstellar medium comprises optically-thin radiative cooling, an external gravitational potential, and background shear. Dynamo coefficients for mean-field models are measured by means of passive test fields.
Results. Our simulations show that supernova-driven turbulence in conjunction with shear leads to an exponential amplification of the mean magnetic field. We found turbulent pumping to be directed inward and approximately balanced by a galactic wind.

Key words: turbulence - magnetohydrodynamics (MHD) - ISM: supernova remnants - ISM: magnetic fields

1 Introduction

Large scale magnetic fields, as observed in disk galaxies, are thought to be the result of a turbulent dynamo. Mean-field models based on the turbulent $\alpha$ -effect can indeed explain most of the observed features in many nearby disk galaxies (cf. Beck et al. 1996). However, uncertainties remain in the estimation of the turbulent transport coefficients and the nonlinear saturation process. Additionally, the helicity conservation in ideal magnetohydrodynamics (MHD) could lead to the so-called catastrophic quenching (Brandenburg & Subramanian 2005), a saturation process that scales with the magnetic Reynolds number. Helicity transport by a wind, as considered by Sur et al. (2007), would be one solution to this scenario.

Recently, corrobative evidence has been found that the degree of regularity of magnetic fields decreases with increasing star formation activity (Chyzy 2007). This has renewed the theoretical interest in dynamo models based on the interaction of supernova-driven turbulence with galactic differential rotation.

Despite the success of mean-field models, until now there has been no clear numerical verification of the dynamo process in the galactic context, and semi-analytical models, which are based on doubtful assumptions, only have limited predictive power (see Gressel et al. 2008, for a recent account).

Balsara et al. (2004) found growing small scale fields for a non-rotating isothermal gas. Similarly, de Avillez & Breit-schwerdt (2005) considered the most realistic model of the interstellar medium (ISM) to date, but did not include the rotation and shear necessary for a mean-field galactic dynamo to operate. Except for their neglect of thermal instability (TI), the simulations of Korpi et al. (1999) are very similar to ours. However, due to the limited computational resources available at the time, their setup suffered from a too small box size, which prohibited a long-term evolution into developed turbulence.

Parker (1992) proposed the magnetic buoyancy instability supported by cosmic rays as a driving mechanism for the galactic dynamo. Hanasz et al. (2006,2004) performed simulations of this process, employing the injection of cosmic ray energy into an adiabatic gas as the driving mechanism. In these simulations, cosmic ray transport is treated in the (anisotropic) diffusion approximation. The authors found an amplification of the mean magnetic field with an e-folding time on the order of $100~\rm Myr$ . Contrary to these cosmic ray models we consider the injection of thermal energy of clustered supernovae and apply a radiative cooling function to reflect the multi-phase nature of the ISM. Although currently one cannot yet resolve all the complex effects of the interstellar plasma, going beyond simple artificial forcing (relying on ad-hoc assumptions), these new approaches mark the first step towards realistic dynamo models based on fundamental physics.

2 Methods

We explore the evolution of the galactic magnetic field within the stratified, turbulent ISM by means of combined 3D, resistive MHD simulations employing the NIRVANA code (Ziegler 2004) and 1D mean-field models based on turbulence parameters.

2.1 Direct simulations

The adopted domain covers a vertically centred box of $0.8\times 0.8\times 4.267~\rm kpc^3$ , representing a local patch of the galactic disk. The physical model includes an external gravitational potential (Kuijken & Gilmore 1989), differential rotation in the shearing box approximation (Gressel & Ziegler 2007), and optically-thin radiative cooling. For the latter, we apply a simple piecewise power-law description (Sánchez-Salcedo et al. 2002; Slyz et al. 2005). In the following, we will briefly introduce the key parameters of our simulations. For a more detailed description of the model, we refer the reader to Gressel et al. (2008).

We choose an initial midplane density of $n_0=1~\rm cm^{-3}$ at a mean molecular weight of $\bar{\mu}=0.6$ . This sets the equilibrium midplane pressure at a value of $p_0/k_{\rm B}=6000~\rm K~\rm cm^{-3}$ . The initial stratification is an exact hydrostatic equilibrium solution in the absence of the kinetic pressure caused by the SNe. Contrary to previous isothermal solutions, our approach applies an effective equation of state given by the additional constraint that the initial disk is in equilibrium with respect to the radiative cooling. The weak initial magnetic field has a radial and azimuthal component with ${B_{\phi}/B_R=-10}$ . The magnitude of the field as a function of height is scaled with (p(z)/p₀)^1/2 to yield a constant plasma parameter of $\beta_{\rm P}=2\!\times\!10^7$ throughout the disk. The rotation period is measured in units of ${\Omega_0=25~\rm km~s^{-1}~\rm kpc^{-1}}$ , while differential rotation is quantified via the shear parameter $q={\rm d}\ln\Omega/{\rm d}\ln R$ . We apply a constant kinematic viscosity of $\nu=5\!\times\!10^{24}~~\rm cm^{2}~s^{-1}$ along with a microscopic magnetic diffusivity of $\eta=2\!\times\!10^{24}~~\rm cm^{2}~s^{-1}$ , resulting in a magnetic Prandtl number of ${\rm Pm}=2.5$ . For reference, we also define the magnetic Reynolds number ${\rm Rm}=L_x^2\Omega_0/\eta \simeq 2500$ .

Turbulence is driven solely via the localised injection of thermal energy. While type-I SNe are exponentially distributed with a scale height of $325~\rm pc$ , for the more frequent type-II SNe, we adopt a prescription where the probability distribution is determined by the vertical gas density profile. This avoids an artificial fragmentation of the inner disk observed in conjunction with a static SN distribution. The reference explosion frequencies are $\sigma_0^{\rm I}=4~\rm Myr^{-1}~\rm kpc^{-2}$ and $\sigma_0^{\rm II}=30~\rm Myr^{-1}~\rm kpc^{-2}$ , respectively, and the associated energies are 10⁵¹ and $1.14\!\times\!10^{51}~\rm erg$ , where the higher value for the latter accounts for the contribution from a stellar wind of the massive progenitor (Ferrière 2001).

2.2 Mean-field model

We follow the common paradigm where the magnetic field $\vec{B}$ and the velocity $\vec{u}$ are separated into mean part $\overline{\vec{B}}$ and $\overline{\vec{u}}$ plus fluctuations denoted by $\vec{B}'$ and $\vec{u}'$ . In this framework, the amplification of the mean field is described via an averaged induction equation that comprises an additional term $\mathcal{E}$ subsuming the action of the turbulent flow, i.e.,

$\begin{displaymath}\partial_{t}\overline{\vec{B}}= \nabla \times \left[\ \left... ...mathcal{E}- \eta \nabla\!\times\!\overline{\vec{B}}\ \right] , \end{displaymath}$

(1)

From the direct simulations described above, we compute the turbulent mean electromotive force $\mathcal{E}= \overline{\vec{u}'\!\times\!\vec{B}'}$ by applying horizontal averages. From the obtained vertical profiles $\mathcal{E}(z,t)$ and $\overline{u}_z(z,t)$ we successfully reconstruct the mean magnetic field $\overline{\vec{B}}(z,t)$ via Eq. (1) demonstrating the applicability of the chosen mean-field approach.

As a simple closure for the mean induction equation, that via $\mathcal{E}$ still depends on the small-scale flow $\vec{u}'$ and magnetic field $\vec{B}'$ , we assume a standard parameterisation of the form

$\begin{displaymath}\mathcal{E}_i = \alpha_{ij} \bar{B}_j - \tilde{\eta}_{ij}\v... ...tial_k \bar{B}_l , \quad i,j \in \left\{R,\phi\right\}, k=z~, \end{displaymath}$

(2)

$\begin{figure} \includegraphics[width=7.6cm,clip]{0195fig1} \end{figure}$	Figure 1: Dynamo coefficients for model H4. Starting at $t=1.3~\rm Gyr$ the TF-method is applied to four consecutive time intervals of $25~\rm Myr$ each. Quantities indicated by the ordinate labels are plotted in dark ( $\alpha _{RR},\dots$ ) or light ( $\alpha _{\phi \phi },\dots$ ) colours, respectively. Shaded areas indicate $1\sigma$ -fluctuations.
Open with DEXTER

While the diagonal elements of $\vec{\alpha}$ produce the dynamo process itself, its antisymmetric off-diagonal part is responsible for the vertical pumping $\gamma_z=\frac{1}{2}(\alpha_{\phi R}-\alpha_{R\phi})$ . Similarly, the diagonal elements of $\tilde{\vec{\eta}}$ are interpreted as turbulent resistivity $\eta_{\rm t}=\frac{1}{2}(\tilde{\eta}_{RR}+\tilde{\eta}_{\phi\phi})$ , while its antisymmetric off-diagonal components can lead to $\vec{\Omega}\times \vec{J}$ -type dynamo effects (Rädler 1969), determined by $\delta_z=\frac{1}{2}(\tilde{\eta}_{R\phi}-\tilde{\eta}_{\phi R})$ .

3 Results

We here compare six models (see Table 1) with varying SN-rate $\sigma$ and rotation frequency $\Omega$ , keeping the shear parameter q and the ratio $\sigma^{\rm I}{:}\sigma^{\rm II}$ fixed. From the initial configuration turbulence builds up smoothly within about $50~\rm Myr$ and the disk reaches a quasi-stationary state. The vertical stratification is now determined by the additional kinetic pressure from the SNe.

$\begin{figure} \includegraphics[width=18cm]{0195fig2} \end{figure}$	Figure 2: Profiles of the regular radial and azimuthal field for model H4 at various times. Results of the simulation (grey lines) are compared with the fields computed from $\mathcal{E}(z,t)$ via Eq. (1) (black lines). At $t\approx 0.85~\rm Gyr$ , a field reversal with pronounced dipolar symmetry occurs.
Open with DEXTER

Table 1: Overview of conducted models. The letters ``Q'', ``H'', and ``F'' indicate quarter, half, and full SN-rate while numbers give the rotation rate. All models include shear with q=-1.

3.1 Dynamo parameters

The most distinctive finding persistent in all our runs is that despite the peaked distribution of the SNe the velocity dispersion increases with galactic height for the inner $2~\rm kpc$ of our box. This is also reflected in the lower two panels of Fig. 1 where we plot the $\vec{\alpha}$ and $\tilde{\vec{\eta}}$ parameters obtained via the test-field method. In excellent agreement with the prediction of second-order theory (SOCA, Rüdiger & Hollerbach 2004), the turbulent diffusivity closely follows the square of the turbulent velocity profile. From this comparison, we can infer a correlation time $\tau$ of the turbulence of $\approx$ $3~\rm Myr$ . We measure $\eta_{\rm t}$ to be around $2~\rm kpc~\rm km~s^{-1}$ , with only a weak dependence on the applied supernova rate. Improving over previous studies, we determine the off-diagonals of $\tilde{\vec{\eta}}$ and find a positive value of $\delta_z\approx 0.5~\rm kpc~\rm km~s^{-1}$ . According to the dispersion relation in Brandenburg (2005), a positive value for $\delta_z$ will, in principle, allow for a contribution of this term to the dynamo effect. Turbulent pumping is found to be directed inward and has an amplitude of $\vert\gamma_z\vert\approx 5~\rm km~s^{-1}$ . Most strikingly, in our simulations the inward pumping is balanced by a wind of approximately the same magnitude. Consequently, the effective vertical transport is drastically reduced, relieving the dynamo process of a major burden.

3.2 Field structure and amplification

In Fig. 2 we plot vertical profiles of the exponentially growing mean radial and azimuthal magnetic field for model H4. We find the field to be largely confined to the inner disk with a radial pitch angle of $-10\hbox{$^\circ$ }$ . Although this is somewhat smaller than the observed values of up to $-30\hbox{$^\circ$ }$ , this value is considerably higher than would be expected from classical $\alpha\Omega$ -dynamos. While the predominant symmetry with respect to the midplane is found to be quadrupolar, this mode is interrupted by field reversals showing dipolar symmetry. This distinct behaviour was successfully reproduced in a 1D toy model, where the frequency of the periodic field reversal depended critically on the interplay of the diamagnetic pumping and the mean vertical velocity profile.

To study the temporal evolution of the arising fields we introduce vertically integrated rms-values $\left<\right.\!\bar{B}_R\!\left.\right>$ and $\left<\right.\!\bar{B}_{\phi}\!\left.\right>$ . The exponential growth of the field is illustrated in Fig. 3. The e-folding time is on the order of $250~\rm Myr$ and varies with the reversals.

In our simulations, the turbulent component dominates over the regular by a factor of 2-3. In particular, we find values $\left<\right.\!\bar{B}\!\left.\right>{:}\left<\right.\!B\hbox{$^\prime$ }\!\left.\right>$ of $0.52~(\pm0.02)$ for model Q4, $0.40~(\pm0.03)$ for model H4, and $0.31~(\pm0.01)$ for model F4. This trend is consistent with observations of strong regular fields in the inter-arm regions of spiral galaxies (Beck 2007). Furthermore, from IR-based star formation rates, Chyzy (2007) observes a correlation $\log\ (B_{\rm reg}{:}B_{\rm tur}) \propto -0.32~(\pm 0.01) \log {\rm SFR}$ . From our values cited above we find a somewhat steeper slope of $-0.38~(\pm0.01)$ with respect to the SN-rate.

From the left panel of Fig. 4, where we plot the time evolution of the regular and fluctuating components of the various models, we see that even the absolute value of the mean field $\left<\right.\!\bar{B}\!\left.\right>$ increases with decreasing SN-rate (models F4, H4, and Q4). This is consistent with the trend in the turbulent diffusivity $\eta_{\rm t}$ , which has values of 2.3, 1.7, and $1.4~\rm kpc~\rm km~s^{-1}$ , respectively.

For the range of parameters studied, we do not observe a significant dependence of the growth rate on the supernova frequency $\sigma$ but only on the rotation rate $\Omega$ . For the models F1-F8, we find e-folding times of $\approx$ 500, 140, 102, and $54~\rm Myr$ for the amplification of $\langle B\hbox{$^\prime$ }\rangle$ . With exception of model F1, which directly corresponds to the parameters used in Korpi et al. (1999), we find exponentially growing regular fields $\langle \bar{B}\rangle$ with e-folding times of 147, 102, and $52~\rm Myr$ for models F2, F4, and F8, respectively. For model F1 the regular field decays at $\approx$ $500~\rm Myr$ . The listed values have been obtained for a time frame of about $100~\rm Myr$ after the turbulence reaches a steady state. Due to the different time base these values are not directly comparable to the long-term growth rate of model H4, where field reversals become important. The initial amplification timescale for the models H4 and Q4 is $90~\rm Myr$ .

4 Discussion and conclusions

We have performed direct simulations of a supernova-driven galactic dynamo and find the rotation frequency to be the critical parameter allowing the dynamo to operate. For our setup, it turns out that $\Omega > 25~\rm km~s^{-1}~\rm kpc^{-1}$ is necessary, which coincides with the prediction of Schultz et al. (1994). Nevertheless, this value may still be dependent on the assumed gas density and the gravitational potential.

While the e-folding time of the amplification mechanism scales with the rotation period, within the range of parameters studied, it only shows a minor dependence on the supernova rate. With $\tau_{\rm e} \simeq 250~\rm Myr$ , the overall growth time for our model H4 is comparable to the value obtained for the cosmic ray driven dynamo (Hanasz et al. 2006) and about four times larger than the expected growth time for the magneto-rotational instability.

For similar models without shear (Gressel et al. 2008), we found no amplification of the mean field. Further simulations will have to show whether the differential rotation is indeed essential for the SN-dynamo to work, or if it was only that the critical values for dynamo action were not achieved in the rigid rotating case given the assumed density profile and gravitational potential. An estimation based on dynamo numbers leads to slightly subcritical values for mean-field dynamo action in the case of rigid rotation.

Although the effects of numerical resolution have to be investigated by further simulation runs, we are confident that the main features are indeed robust with respect to the numerical modelling. The good agreement of the direct simulations with the 1D toy model, as well as the absence of field amplification in the case without Coriolis forces, provide further evidence in support of the classical picture of cyclonic turbulence. That is, the Coriolis force plays the central role in giving the expanding remnants a characteristic handedness. The mean vertical velocity, together with the turbulent pumping term, add a new timescale to the system, which is probably responsible for the fast growth. This was already observed in mean-field models with a prescribed wind in the mean velocity (Elstner et al. 1995).

$\begin{figure} \includegraphics[width=8.8cm]{0195fig3} \end{figure}$	Figure 3: Evolution of the turbulent and regular magnetic field for model H4. For $\langle \bar{B}_R\rangle$ and $\langle \bar{B}_{\phi}\rangle$ we show the results from the direct simulation (grey lines) together with the reconstruction from $\mathcal{E}(z,t)$ (black lines).
Open with DEXTER

$\begin{figure} \includegraphics[width=8.8cm]{0195fig4} \end{figure}$	Figure 4: Comparison of the time evolution of the regular and fluctuating magnetic field strength over kinetic energy for the various setups (cf. Table 1). Model F4 $^{\star }$ is identical to model F4 with Coriolis forces disabled and demonstrates that no field amplification is obtained in the case of Cartesian shear.
Open with DEXTER

The relative strength of the regular field compared to the turbulent field is determined by the SN rate. The resulting slope of $-0.38~(\pm0.01)$ is similar to the observed one derived from comparisons of polarised over total intensity of radio synchrotron emission with SFRs (Chyzy 2007). The strong radio-FIR correlation also supports a strong dependence of the total magnetic field strength with the SFR. Therefore, the saturation level of the dynamo with respect to the SN-rate should be compared in the future.

Since we do not observe evidence for catastrophic quenching, our results suggest that the combined action of the diamagnetic pumping and the wind might cause an outward helicity flux while retaining the mean field. This poses a natural solution to the quenching catastrophe and will be the subject of further investigations at higher magnetic Reynolds numbers.

Direct simulations of a supernova-driven galactic dynamo