text:
# Thermal X-ray emission from shocked ejecta in Type Ia Supernova Remnants. Prospects for explosion mechanism identification

## 1 Introduction

Type Ia supernovae (SNIa) play an important role in the chemical evolution of galaxies, originating most of the Fe group elements in the interstellar medium, and they have become a key tool in our understanding of the universe by providing evidence for its accelerated expansion (Perlmutter et al., 1999). Still, our knowledge of the physical processes involved in the actual explosions is far from being complete. Fundamental issues such as the mass of the white dwarf (WD) at the moment of the explosion, the location of the ignition and the propagation mode of the burning front have not been established yet, and a number of different models or paradigms are capable of reproducing with reasonable accuracy the optical spectra and light curves of observed SNIa and their fundamental physical properties (Branch and Khokhlov, 1995; Hillebrandt and Niemeyer, 2000, for reviews).

Given this situation, it is important to explore all the potential sources of information that can shed some light on the detailed workings of thermonuclear supernovae. The thermal X-ray emission from the shocked ejecta in young supernova remnants (SNRs) originated by these explosions is, in principle, well suited for this purpose: as the reverse shock advances into the ejecta, it compresses and heats them to X-ray emitting temperatures, producing spectra that depend on the composition, density, temperature and ionization state of the material. The calculation of the expected thermal X-ray emission from the shocked ejecta synthesized in the different explosion paradigms, however, is not straightforward. The propagation of the shocks and the density of the shocked gas is governed by the dynamic evolution of the SNR, which at early stages depends strongly on the ejecta density profile and is not well represented by similarity or unified solutions (Truelove and McKee, 1999). The densities involved are so low that the ensuing shocks are collisionless, and their physics, and particularly how the postshock internal energy is distributed among ions and electrons, is not well understood. Moreover, the ages of young SNRs (up to a few thousand years) are short compared with the time scales for thermal equilibration between ions and electrons and for the onset of collisional ionization equilibrium, and therefore the plasma is in a transient ionizing state that is difficult to calculate.

Our goal is to provide observers with the means to make meaningful analysis of the thermal X-ray spectra of young, ejecta-dominated type Ia SNRs and eventually to constrain the nature of the event that originated them. Thus, here we present the first detailed models for these remnants, evolved from a set of theoretical calculations of thermonuclear supernova explosions. We have focused our attention on the general properties of each explosion model, and their consequences on the emitted X-ray spectra. A detailed analysis of particular objects is a delicate task that requires fitting many observables and an exhaustive exploration of the model parameter space in each case. This task is left for forthcoming publications. In section 2 we give the relevant characteristics of the explosion models in our grid, which includes all the paradigms currently under discussion. The SNR simulation scheme is described in section 3: hydrodynamics, plasma ionization and heating, and synthetic spectra. Section 4 contains a discussion of our results and the relevant observable quantities that can be derived from the models, as well as a preliminary comparison with observations for Tycho SNR. The conclusions of this study are presented in section 5.

## 2 Supernova explosion models

In order to be able to discriminate between different explosion mechanisms it is important that all the models be calculated consistently, i.e. with the same physics included in the same hydro and nucleosynthetic codes. With this requirement in mind, we have computed a grid of thermonuclear supernova explosion models which is representative of the whole diversity of theoretical models currently under debate. All the calculations have been performed in one dimension, assuming spherical symmetry. For each model, we have followed the explosion with a supernova hydrocode until $`10^6`$ s after the ignition, when the expansion has reached a nearly homologous ($`rv`$ everywhere) stage. The deposition of the energy of radioactive decay of 56Ni on the ejecta is taken into account in this supernova code. After $`10^6`$ s, however, the 56Ni disintegration becomes dynamically irrelevant because: 1) most of it has already decayed to 56Co, and 2) an increasing fraction of the energy of the photons escapes the supernova due to the drop in opacity caused by the expansion. In addition to the detailed nucleosynthesis, we have computed the light curves in order to be able to compare with historical supernovae. The codes used in the explosion calculations are the same as in Bravo et al. (1996). More details are given in the appendix.

The different categories of SNIa explosion models that we have included in our grid are:

* DET: Pure detonation model. In this class of models, the flame starts close to the center of the WD, and propagates supersonically nearly through the whole white dwarf, incinerating most of it to Fe-group elements.
* SCH: Sub-Chandrasekhar mass model. A He detonation is started at the edge of a He envelope, which feeds a converging shock wave into the C-O core. Close to the center, the converging shock wave transforms into a C-O detonation which propagates outwards and processes the rest of the core. A sandwiched structure is produced, rich in Fe-group elements both in the inner and in the outer parts of the ejecta (in our model, below a Lagrangian mass of $`0.4\mathrm{M}_{\mathrm{}}`$ and above $`0.8\mathrm{M}_{\mathrm{}}`$, in what was the He envelope, see Fig. 1), and rich in intermediate-mass elements plus C-O in between.
* DEF: Pure deflagration models. In these models the deflagration propagates at the laminar flame velocity (a small fraction of the sound velocity) close to the center, until the Rayleigh-Taylor instability develops, deforms the flame surface, and accelerates the combustion. The flame velocity remains subsonic all the way, but when the material has reached velocities that are of the same order as that of the flame the expansion quenches the flame. A large mass of unburnt C-O is ejected in the outer layers.111Other deflagration models can be found in the literature in which the mass of unburnt C-O is small while the mass of intermediate-mass elements is large. This is the case, in particular, of the popular W7 model (Nomoto, Thielemann, and Yokoi, 1984).
* DDT: Delayed detonation models. In these models the flame propagates initially as a slow deflagration, but a transition to a detonation is induced at a prescribed flame density during the expansion of the white dwarf. As a result, the otherwise unburnt C-O is processed partly into Fe-group and partly into intermediate-mass elements.
* PDD: Pulsating delayed detonation models. They differ conceptually from delayed detonation models in that the transition to detonation is induced only after the white dwarf has pulsed. The pulsation is due to the inefficient burning produced by a slow deflagration, which is unable to rise the whole energy of the white dwarf (gravitational \+ internal + kinetic) above zero.

Our model grid incorporates examples of all of these paradigms, and it aims to cover the whole parameter space. For this paper we have selected eight models from this grid as a representative sample, including extreme cases of DEF, DDT, and PDD. Their properties are given in table 1. There, $`E_\mathrm{k}`$ is the kinetic energy of the ejecta, $`_{\mathrm{m}ax}`$ and $`\mathrm{\Delta }_{15}`$ are, respectively, the bolometric magnitude of the supernova at light curve maximum and the change in bolometric magnitude between maximum and 15 days later, and the rest of the columns are clearly defined by the header. The meaning of the parameters in the second column is given in the appendix.

In Figure 1 we show the chemical structure of each model after the short lifetime radioactive isotopes have completed their decays (all the Fe that appears on the plots, for instance, was synthesized as 56Ni in the explosion). Their density structures are shown in Figure 2. In the delayed detonation models, the transition from deflagration to detonation happened at a Lagrangian mass of $`0.2\mathrm{M}_{\mathrm{}}`$, where its imprint on the density profile can be seen. The chemical structure is dominated by Fe-group elements up to $`1.0\mathrm{M}_{\mathrm{}}`$ for DDTa, at which point the flame density was too low to incinerate matter to nuclear statistical equilibrium (NSE), leaving a buffer of intermediate-mass elements. For DDTe this buffer is much larger, because of the lower densities achieved by the detonation. The outermost zone dominated by C-O is also larger than in model DDTa. The chemical structure of the pulsating delayed detonation model PDDa is very similar to that of DDTa, with the transition to detonation happening at $`0.3\mathrm{M}_{\mathrm{}}`$. The different hydrodynamical histories of both models, however, are reflected in their density profiles (see Fig. 2). At the time of the transition, the density of the external layers of PDDa was on average about two orders of magnitude lower than in DDTa. As a result, PDDa displays a tail of very low density but larger radii than DDTa. In model PDDe the transition to detonation took place at $`0.3\mathrm{M}_{\mathrm{}}`$. Its chemical profile is similar to that of DDTe, while its density profile is similar to that of PDDa. In the deflagration model DEFa, the flame quenched due to the expansion of the WD at a lagrangian mass of $`0.7\mathrm{M}_{\mathrm{}}`$, after which a narrow ($`0.1\mathrm{M}_{\mathrm{}}`$ wide) region rich in intermediate-mass elements was formed. The density structure is very different from the delayed detonation models, due to the sudden termination of nuclear energy generation which results in the formation of a bump of unburnt material just above the quenching flame front. In model DEFf the results were similar to DEFa, but the flame was quenched at $`0.9\mathrm{M}_{\mathrm{}}`$.

## 3 Simulation scheme

### 3.1 Overview

The density profiles from the explosion models described in the previous section are used as an input for a standard hydrodynamics code (see section 3.2 for details) that follows the interaction of the ejecta with a uniform ambient medium (AM). Hydrodynamical simulations give the evolution of the radius, velocity, density and internal energy per unit mass for each fluid element as a function of time, i.e., its dynamic history. This dynamic history, together with the chemical composition of the fluid element obtained from the explosion models is then used to calculate the evolution of the plasma heating and the nonequilibrium ionization (NEI) in a self-consistent way (section 3.3), resulting in electron temperatures and ion fractions for all the ions of the relevant chemical elements at each time. The density, electron temperature and ionization structure within a certain region of the SNR at a certain time then lead to a synthetic or predicted thermal X-ray spectrum, which can be readily calculated with a spectral code (section 3.4) and convolved with any particular instrumental response in order to compare the models with observations. In section 3.5 we will review the approximations made at each stage of our simulation method and ascertain its validity for the study of young, ejecta-dominated Type Ia SNRs.

### 3.2 Hydrodynamics

We have built a 1D Lagrangian hydrodynamics code such as described by Truelove and McKee (1999), with an ideal gas equation of state ($`\gamma =5/3`$), nonlinear pseudoviscosity and no external energy sources or sinks, so the SNR model remains adiabatic throughout its evolution. We have followed the interaction of the ejecta from various explosion models with a constant density ambient medium (AM), with $`\rho _{\mathrm{A}M}=10^{24},\mathrm{\hspace{0.17em}5}\times 10^{25}`$ and $`10^{25}`$ g cm-3 as a sample representative of the interstellar medium (ISM) conditions in most of the Galaxy and the Magellanic Clouds.

### 3.3 Plasma model and ionization calculations

The theory for modeling transient heavy element plasmas in young SNRs was developed by Hamilton and Sarazin (1984, henceforth HS84). These authors describe the plasma with three Maxwellian populations: one for the ions (including neutrals), and two for the electrons, a ’cold’ component for the electrons produced in the ionization process and a ’hot’ component for the electrons already present in the preshock gas which were heated at the collisionless shock front. The relative importance of the hot electron population thus depends on the amount of collisionless heating at the shock, on the preshock ionization state of the ejecta and on the composition of each layer, and in general it will decrease with time as more ’cold’ electrons are produced by postshock ionization, especially if collisionless electron heating is not efficient. The efficiency of this heating remains a controversial issue: while theoretical calculations of the effect of plasma instabilities on the electron population suggest a high degree of collisionless heating at the shock (Cargill and Papadopoulos, 1988), this is difficult to reconcile with observations of high velocity shocks, and this prediction is being revised (see Laming 2000, for a discussion of collisionless electron heating). We allow only for a small amount of collisionless heating at the reverse shock, so we have represented the electrons by a single Maxwellian population $`T_\mathrm{e}`$. In this two-fluid model, the postshock temperature ratio between ions and electrons, $`T_{\mathrm{e},s}/T_{\mathrm{i},s}`$, is a measure of the efficiency of collisionless heating, varying between 0 (no heating) and 1 (total temperature equilibration). For the models presented in this paper, we have assumed $`T_{\mathrm{e},s}/T_{\mathrm{i},s}=0`$, but we will explore the effect of different assumptions in future work.

All the elements in the unshocked ejecta are assumed to be in the singly ionized state. A low preshock ionization state agrees with the observational evidence provided by Wu et al. (1983) for SN 1006, one of the prototype galactic Type Ia SNRs, but we have found that varying the charge of the unshocked ions has little impact in the postshock ionization history, at least for charge states below ten. For a more in-depth discussion of preshock ionization, see HS84.

In the shocked plasma, ions and electrons interact through Coulomb collisions which gradually equilibrate the Maxwellian populations. At the same time, the ongoing ionization of the ejecta increases the electron population and therefore modifies the collision rate. Let $`\rho `$ and $`\epsilon `$ be the time-dependent mass density and specific internal energy per unit mass of a given fluid element as calculated with the hydrodynamic code, and $`n_\mathrm{e},n_\mathrm{i}`$ the number densities of electrons and ions, respectively. At any given time, the total internal energy is distributed among ions and electrons so that $`\epsilon _\mathrm{i}+\epsilon _\mathrm{e}=\epsilon `$, and the respective temperatures are related to the internal energies per unit mass in each population by $`\epsilon _{\mathrm{e},i}=(3kT_{\mathrm{e},i}n_{\mathrm{e},i})/2\rho `$. We can represent the chemical composition and ionization state of the fluid element with the normalized number abundances for each element X $`f_X=n_X/n_\mathrm{i}`$ and each ion q $`f_{X^q}=n_{X^q}/n_X`$ with $`q=0`$ for neutral and $`q=Z_X`$ for bare ions, so that $`_Xf_X=1,_qf_{X^q}=1`$. The plasma then evolves according to the following equations (adapted from HS84):

$$\frac{d\epsilon _\mathrm{i}/\epsilon }{dt}=\frac{d\epsilon _\mathrm{e}/\epsilon }{dt}=\frac{1}{\rho \epsilon }\frac{2^{5/2}\pi ^{1/2}e^4n_\mathrm{e}n_\mathrm{i}\overline{Z}^2\mathrm{ln}\mathrm{\Lambda }}{m_\mathrm{e}\overline{A}m_\mathrm{u}k^{1/2}}\frac{(T_\mathrm{i}T_\mathrm{e})}{\left(\frac{T_\mathrm{i}}{\overline{A}m_\mathrm{u}}+\frac{T_\mathrm{e}}{m_\mathrm{e}}\right)^{3/2}},$$
(1)
$$\frac{df_{X^q}}{dt}=\frac{\overline{Z}\rho }{\overline{A}m_u}\left[I_{X^{q1}}f_{X^{q1}}+R_{X^{q+1}}f_{X^{q+1}}\left(I_{X^q}+R_{X^q}\right)f_{X^q}\right],$$
(2)
where $`\overline{A}=_Xf_XA_X`$ is the average mass number, $`\overline{Z}=_Xf_X_qqf_{X^q}`$ the average ion charge in the fluid element, $`m_\mathrm{u}`$ the atomic mass unit, $`\mathrm{ln}\mathrm{\Lambda }`$ the Coulomb logarithm, all fundamental constants take their usual values and all time derivatives are Lagrangian. $`I_{X^q}`$ and $`R_{X^q}`$ represent the ionization and recombination rates *from* ion q of element X, respectively, and are taken from the recent compilation by Mazzotta et al. (1998). Both equations are obviously coupled in transient plasmas dominated by heavy elements (where $`d\overline{Z}/dt0`$), and they have to be integrated simultaneously. At the same time, the processes are controlled by the precalculated variation of the hydrodynamic quantities $`\rho `$ and $`\epsilon `$. The resulting stiff differential equation system is solved with an implicit scheme, resulting in a time series of temperatures and ion fractions for each fluid element. The values of $`T_\mathrm{e}`$ and $`T_\mathrm{i}`$ are recomputed from the updated values of $`\epsilon _{\mathrm{e},i}`$ and $`n_{\mathrm{e},i}`$ after each new iteration of equations (1) and (2). All the chemical elements with $`f_X10^3`$ have been included in the calculations.

### 3.4 Synthetic spectra

To produce the synthetic spectra, we have used an updated version of the Hamilton & Sarazin (HS) code incorporated in the XSPEC software package (Arnaud, 1996), which is suitable for modeling plasmas that depart significantly from the collisional ionization equilibrium (CIE). Given the density, ionization state and temperature of a fluid element at a certain time, the HS code produces a predicted spectrum that can be convolved with the response matrix of any appropriate instrument. It is important to mention that this code relies on modern atomic calculations for the Fe L shell emission from Li-like to Ne-like Fe ions (Liedahl, Osterheld, and Goldestein, 1995). This code includes no atomic data for Ar, and has other limitations which restrict its applicability, but problems of this kind are common to most of the currently available NEI codes. The limitations of the spectral code should be considered when comparing the predicted spectra with observations, but the quality of our calculated spectra is satisfactory for our present purposes.

We will only discuss the integrated ejecta spectra: detailed models for the shocked AM spectra and spatially resolved spectroscopy of the ejecta are beyond the scope of this paper, even though they can certainly be addressed using the simulation scheme we have just described. We have used the response matrix for the CCD MOS-1 EPIC camera onboard the XMM-Newton observatory, which provides an excellent spectral resolution and is a reasonable choice for observing extended objects like SNRs. The Doppler broadening of lines due to thermal and bulk motions have not been included because their effects are small at the spectral resolution of XMM EPIC-MOS.

### 3.5 Simulation scheme discussion

It is important to emphasize that our simulations assume spherical symmetry for both the SNIa ejecta and the interaction with the AM. Our aim is to show that young SNRs can be a valuable tool for studying SNIa explosions, and the assumption of one-dimensional dynamics is a first step towards that goal. In fact, it can be a satisfactory approach in many individual cases, especially since Type Ia SNRs tend to be, on average, more symmetric than those originated by core collapse supernovae. In contrast to core collapse SNRs, they seem to retain the original ejecta stratification (compare the results of Decourchelle et. al. 2001 on Tycho and Lewis et. al. 2003 on N103B with those of Hughes et. al. 2000 on Cas A). Whenever deviations from spherical symmetry affect the shape or the spectrum of a given object, a comparison with our models might be useful as a description of the overall ionization state and some bulk properties of the remnant, but it should be done with care. Moreover, the effect that the Raleigh-Taylor instabilities acting on the contact discontinuity (CD) might have on the ionization and heating of the outer ejecta layers is difficult to ascertain. As explained in Wang and Chevalier (2001), the averaged density profile is smoothed by the instabilities, but the fingers of shocked ejecta retain a higher density than the surrounding shocked AM, so the results of our one dimensional simulations discussed in the following sections might not deviate much from a more realistic case.

Even if spherical symmetry is assumed, it is worthwhile to note that, while SNIa progenitor systems are not supposed to substantially modify the surrounding ISM, the constant density hypothesis does not agree with current presupernova evolution models. These models imply a substantial modification of the AM, with profound effects in the dynamic evolution of the remnant, as discussed in Badenes and Bravo (2001). This possibility has not been taken into account here, but will be explored in future work.

Another issue that is of special concern is the validity of the adiabatic hypothesis in the hydrodynamic calculations: whereas radiative losses are generally not important for solar composition plasmas within time scales of a few thousand years, heavy element plasmas will radiate at a much faster rate, and the losses might have noticeable effects on the dynamics much earlier in the evolution of the SNR. The problem is complicated by the fact that losses are composition and ionization state dependent, so they can only be evaluated *a posteriori*, after the ionization calculations described in section 3.3 are completed, and it is not possible to include their effect in the hydrodynamics with our simulation scheme due to the fact that both calculations are decoupled. This *a posteriori* monitoring of the radiative losses in a heavy element NEI plasma can be done using the atomic data from Summers and McWhirter (1978) as described by Laming (2001). For evaluation purposes, we define the time scale for the onset of radiative losses $`t_{\mathrm{r}ad}`$ , as the time when the calculated *a posteriori* losses exceed 10% of the internal energy in a number of layers that amount to 5% of the total ejecta mass. This time scale is below 5000 years only for the deflagration models: $`t_{\mathrm{r}ad,DEFa}=3.03\times 10^{10}s;t_{\mathrm{r}ad,DEFf}=2.37\times 10^{10}s`$ for an interaction with $`\rho _{\mathrm{A}M}=10^{24}`$ g cm-3, with $`t_{\mathrm{r}ad}`$ increasing for lower AM densities. Because radiative cooling is a runaway process, the validity of these models beyond $`t_{\mathrm{r}ad}`$ is difficult to determine, but it is worthwhile to notice that the losses are usually confined to a very small volume and the effect on the overall dynamics of the SNR should not be important (Hamilton, Sarazin, and Szymkowiak, 1986).

## 4 Results

### 4.1 Hydrodynamics

The interaction of Type Ia ejecta with a constant density AM, based on a grid of one dimensional thermonuclear supernova explosion models, was first explored in Dwarkadas and Chevalier (1998, hereafter DC98). In that paper, the detailed dynamics of the interaction of six ejecta density profiles from thermonuclear supernova models was examined and compared with three analytical density distributions: an exponential, a power law of index $`n=7`$ with a constant density core, and a constant density profile. Except for the exponential, these analytical functions have been widely used before to approximate SNIa ejecta density as a function of velocity. Approximate temperature profiles were also calculated assuming solar abundances for the supernova ejecta. DC98 conclude that in all cases the density of the shocked ejecta increases from the reverse shock towards the contact discontinuity, that this density rise is usually coupled with a drop in temperature in the area close to the CD and that the sharp structures in the ejecta profiles, especially in the He detonation models, give rise to secondary waves propagating in the interaction region which could affect the instantaneous X-ray emissivity of the remnants.

When analyzing the ejecta-AM interactions, it is important to note that, once the ejected mass $`M_{\mathrm{e}j}`$, kinetic energy $`E_\mathrm{k}`$ and density profile of an explosion model are fixed, the interaction with the AM follows a scaling law for the AM density $`\rho _{\mathrm{A}M}`$ (Gull, 1973). We have tested this scaling law by computing the interaction of the supernova models with the three different AM densities mentioned in section 3.2. Our results follow the scaling law up to the precision of the hydrodynamic calculations. Therefore, it is sufficient to present the hydrodynamic calculations for each model with a certain value of $`\rho _{\mathrm{A}M}`$ and then use the characteristic magnitudes defined by DC98 (eqs. , , and ) to scale the results to any other $`\rho _{\mathrm{A}M}`$ value that might be of interest.

$$R^{}=\left(\frac{M_{\mathrm{e}j}}{\frac{4\pi }{3}\rho _{\mathrm{A}M}}\right)^{1/3},$$
(3)
$$V^{}=\left(\frac{2E_\mathrm{k}}{M_{\mathrm{e}j}}\right)^{1/2},$$
(4)
$$t^{}=\frac{R^{}}{V^{}}=\frac{M_{\mathrm{e}j}^{5/6}}{\left(\frac{4\pi }{3}\rho _{\mathrm{A}M}\right)^{1/3}\left(2E_\mathrm{k}\right)^{1/2}}.$$
(5)
In Figure 3 we show our results for the ejecta-AM interaction with $`\rho _{\mathrm{A}M}=10^{24}`$g cm-3. We have plotted the time evolution of the forward and reverse shock radii ($`r_{\mathrm{f}wd},r_{\mathrm{r}ev}`$), the velocity of the forward shock $`u_{\mathrm{f}wd}`$, the velocity of the reverse shock in the rest frame of the expanding ejecta $`u_{\mathrm{r}ev}=(r_{\mathrm{r}ev}/t)(dr_{\mathrm{r}ev}/dt)`$ and the expansion parameters for both shocks, defined as $`\eta _{\mathrm{f}wd,rev}=d\mathrm{ln}(r_{\mathrm{f}wd,rev})/d\mathrm{ln}(t)`$. The time axis spans between 50 and 5000 years after the explosion, and the reverse shock parameters have been plotted only up to the time of rebound. In fact, our simulations of the interaction of the ejecta with the AM start at $`10^7`$ s after the explosion, but we present our results only from year 50 on because,generally speaking, younger remnants are not expected to emit appreciably in X-rays.

The dynamics of the forward shock is affected by the differences in the density profiles during the first thousand years, then all the models converge towards the Sedov-Taylor solution ($`\eta _{\mathrm{f}wd}=0.4`$); this transition leads to the change of slope of the forward shock radii in the log-log plot (Fig. 3a). The shock trajectories of the deflagration models (DEFa, DEFf) lag behind the others because their $`E_\mathrm{k}`$ is lower, and they can also be easily distinguished by their high $`\eta _{\mathrm{f}wd}`$ values at early times, (Fig. 3e), about 50% higher than in the other models. Sudden increases in $`\eta _{\mathrm{f}wd}`$ can be seen in the PDD models around $`t=3\times 10^9`$ s, and in the SCH model at $`t=6\times 10^8`$ s (the rise is outside the range of Figure 3e, but its effect can be seen at the beginning of the plot) and $`t=6\times 10^9`$ s. In general, high $`\eta _{\mathrm{f}wd}`$ values are found in models that have high density material in the outermost ejecta, as the DEF models, or high density layers in the outer ejecta preceded by a low density tail, as PDDa at $`\text{r}=1.6\times 10^{16}`$ cm, PDDe at $`\text{r}=1.4\times 10^{16}`$ cm, and SCH at $`r=1.9\times 10^{16}`$ cm and $`r=1.3\times 10^{16}`$ cm (see Fig. 2). These high density layers in the outer ejecta transfer their momentum to the shocked material and to the forward shock, leading to the increased $`\eta _{\mathrm{f}wd}`$ that we have noted. The high density material, once shocked, stays close to the contact discontinuity, and the low density tail, if present, is also recompressed by the reflected shocks that ensue when the high density layers are overcome by the reverse shock. The result is that the density enhancement effect close to the contact discontinuity is stronger for the DEF, PDD and SCH models than for the DDT and DET models, where the momentum transfer is more gradual (for examples of density maps see Fig. 3 in Badenes and Bravo 2003b). The dynamics of some of these models are compared to those of the power law and exponential analytical profiles in Badenes and Bravo (2003a).

Summarizing, our conclusions agree almost completely with those of DC98: the density of the ejecta always peaks towards the contact discontinuity for all models at all times, and the rich internal structure of the explosion models produces a series of secondary shock waves that travel along the shocked ejecta and AM reheating and recompressing the material. Our results also agree with DC98 in that the mean temperature profile drops towards the contact discontinuity, but we have found that this behavior of the mean temperature cannot be extrapolated to the electron temperature, as is explained in the following section (see also Fig. 2 in Badenes and Bravo 2003b). We have also found that the models with dense layers in the outer ejecta, that is, the deflagration models and, to a lesser extent, the PDD models and the SCH model, tend to achieve higher densities in the region behind the contact discontinuity. This will have a profound impact in the ionization and temperature calculations.

### 4.2 Temperature, ionization, and emission measure

The spectral characterization of a young SNR is a complex issue: since the density, ionization state, temperature and composition of each fluid element are different, each region of the SNR will contribute differently to the total integrated spectrum, and so will each chemical element. A convenient way to measure these contributions is the emission measure (EM) for element X, defined as

$$EM_X=_{V_{\mathrm{s}h}}n_Xn_\mathrm{e}𝑑V,$$
(6)
where $`V_{\mathrm{s}h}`$ is the volume of shocked ejecta. For identical physical conditions and a common history, elements with equal emission measures contribute equally to the total integrated ejecta spectrum. But the electron temperature $`T_\mathrm{e}`$ and ionization time scale of the plasma, $`\tau =n_\mathrm{e}𝑑t`$, which play a key role in thermal NEI spectra, are different for each fluid element, resulting in different spectra produced by fluid elements with identical emission measures. Borkowski, Lyerly, and Reynolds (2001) approached this problem for the shocked AM in Sedov SNRs by introducing distribution functions, plots of $`T_\mathrm{e}`$ and $`\tau `$ versus EM, but the use of this approach for the shocked ejecta would call for an individual distribution function for each chemical element due to the nonuniform chemical composition. An incomplete, yet meaningful, description can be achieved by taking the first moment of the distribution functions and calculate, for each element X, an emission measure averaged electron temperature $`T_\mathrm{e}_X`$ and ionization timescale $`\tau _X`$. (Another quantity, an ionization timescale averaged electron temperature is generally necessary for a reliable modeling of X-ray spectra, but it is less important than $`T_\mathrm{e}_X`$ and $`\tau _X`$ ). We stress that we do not recommend single temperature, single ionization timescale spectral models to fit shocked ejecta in young SNRs, we are merely using the averaged quantities to describe average physical conditions of various layers of ejecta in our models. In all the spectra and calculations presented in the following sections we have used the full profile of electron temperatures and ionization timescales throughout the ejecta, not the averaged quantities.

The $`EM_X`$ evolution is shown in Figure 4 for $`\rho _{\mathrm{A}M}=10^{24}`$ g cm-3, between 20 and 5000 years after the explosion. The contributions from different chemical elements to the ejecta spectra depend strongly on both the composition profile of the models and their dynamic evolution. The density enhancement effect towards the contact discontinuity makes the chemical elements in the outer layers of the ejecta more prominent than those in the inner layers, so the importance of Fe in the model spectra is generally much less than might be expected on the basis of a Type Ia elemental composition alone. This relatively low prominence of Fe in spectra of Type Ia SNR candidates might have been noticed on a number of occasions, often accompanied by inordinately high apparent abundances of other elements (see Hendrick, Borkowski, and Reynolds 2003 for the SNRs 0548-70.4 and 0534-69.9 and Lewis et al. 2003 and van der Heyden et al. 2002 for N103B). The estimated ejected masses of various elements, which often indirectly rely on the assumption that chemical abundances are proportional to the fitted $`EM`$ for each element in the spectrum, are difficult to reconcile with the yields of theoretical explosion simulations. This discrepancy between $`EM`$ and ejected mass is most dramatic in the deflagration models (Fig. 4 e and f), whose spectrum is completely dominated by C and O, with emission measures of Fe about two orders of magnitude lower at all times, even though the ejected mass of Fe is higher than that of C or O. Note, however, that the peak value of the $`EM_{\mathrm{F}e}`$ in the deflagration models is about the same as in the other models, with the exception of PDDa, and that radiative cooling could reduce considerably the $`EM`$ of C and O in the deflagration models (see section 3.5). The prompt detonation model (Fig. 4g) is the only one whose spectrum is clearly dominated by Fe at all times, while in DDTa and PDDa (Fig. 4a and c) Fe takes over only after a few hundred years, with important contributions of Si and S (and O for DDTa) throughout the SNR evolution. The rise in the emission measures around $`t=10^{11}`$ s is due to the propagation of the reverse shock after it rebounds at the center, reheating and recompressing the ejecta; this rise happens earlier and is more gradual for the elements in the inner layers than for those in the outer layers. Model SCH is dominated by Fe only when the reverse shock is propagating through the He detonation layer, afterwards O takes over, because a second density enhancement region forms behind the interface between this layer and the rest of the exploded WD.

Since the ionization and electron heating processes proceed faster at higher densities, the corresponding $`T_\mathrm{e}_X`$ and $`\tau _X`$ plots in figures 5 and 6 are also affected by the enhancement towards the contact discontinuity and the reverse shock rebound, albeit in different ways. Those models with stronger density enhancement close to the CD (DEFa, DEFf, PDDe, and, to a lesser extent, PDDa and SCH) tend to have higher ionization timescales for most elements, and in general, elements closer to the CD have higher $`\tau _X`$. In the deflagration models, for instance, C and O are always at a higher $`\tau `$ and initially hotter than other elements. The $`T_\mathrm{e}_{\mathrm{C},O}`$ drops at later times because in these models the density is so high that electron-ion temperature equilibration is achieved for most of the C and O in the ejecta before $`t=10^{10}`$ s and afterwards the electrons just cool due to adiabatic expansion of this region of the SNR. Iron is generally hotter and at a higher $`\tau `$ in the DDT and PDD models than in the DEF models. The anomalous behavior of the plots for some elements (for instance, Ni in DDTa and PDDa or Ca in SCH) is due to the averaging in $`EM`$ and can be understood by comparing the curves with the chemical composition profiles of Fig.1. As the reverse shock advances into regions with a much higher concentration of a given element, the newly shocked (and therefore cooler and less ionized) layers soon dominate the $`EM_X`$ and the averaged quantities shift their values accordingly. It is worth noting that the enhanced electron heating rate due to the higher densities towards the CD compensates for the lower specific internal energies found in that region, so we observe $`T_\mathrm{e}`$ profiles that always peak at the CD (i.e. electrons and ions are closer to thermal equilibrium at the CD, whereas $`T_\mathrm{e}T_\mathrm{i}`$ behind the reverse shock) in contrast to the behavior of the mean temperatures (DC98).

Our results for lower AM densities are similar, and we have found that the hydrodynamical scaling laws (section 4.1) also work reasonably well for the emission measures, except at extremely low AM densities. Scaled $`EM_X`$ agree to within a factor of 2, and so do the $`\tau _X`$. The temperature scaling seems more complex. The deviations are accounted for by the difference in average ionization state which is expected in SNRs that evolve in different $`\rho _{\mathrm{A}M}`$, because hydrodynamical scaling does not apply to ionization and electron heating processes. A detailed analysis of the ambient density effect on the spectral characteristics of our models will be the subject of future work.

### 4.3 Synthetic spectra

The integrated synthetic spectra from the ejecta in our SNR models are presented in figure 7 for $`\rho _{\mathrm{A}M}=10^{24}`$ g cm-3, 500, 1000, 2000, and 5000 years after the explosion (without interstellar absorption). As explained in the previous section, the contributions from the different elements to the spectra depend strongly on the details of the hydrodynamic evolution, and in particular on the density enhancement towards the CD. Note how the spectra of the deflagration models are always dominated by C and O, to the point that their continua ’veil’ the lines of the other elements at early times. Even though their composition profiles are similar, PDDa and PDDe have a richer line spectrum than DDTa and DDTe respectively, because their mean ionization state is more advanced due to the higher densities close to the CD (section 4.1). The strength of the Fe $`K\alpha `$ line varies from model to model, being important at all times in DDTa, PDDa, DET and SCH, which have more Fe in the outer layers of ejecta, while for DDTe and PDDe Fe K$`\alpha `$ is noticeable only after the first 1000 years, for DEFf after 2000 years and much later for DEFa. Oxygen can be easily identified in the models with high $`EM_\mathrm{O}`$ and low $`EM_{\mathrm{F}e}`$ (otherwise its presence is partially masked by the Fe L complex): DDTe and PDDe at early times, SCH up to 2000 years and the deflagration models at all times. These oxygen bright models with weak Fe lines could be easily mistaken for core-collapse SNRs. For the deflagration models in particular, their strong C and O continua could make spectral analysis difficult in presence of substantial interstellar absorption and in view of calibration problems and low spectral resolution of CCD detectors at low photon energies.

### 4.4 Comparison with observations: the Tycho SNR

A comparison with observations is essential for assessing validity of theoretical models as a spectral characterization of young type Ia SNRs and for learning about their progenitors. Our models are well suited for qualitative comparison with observed spectra, but a quantitative, detailed study of a given object is much more difficult. The excellent quality of present day observations from XMM-Newton and Chandra has surpassed our ability to model them accurately even with sophisticated theoretical models, and additional issues such as modeling of shocked AM spectrum, nonthermal spectrum subtraction, matching shock velocities, expansion parameters and apparent radii, and spatially resolved spectroscopy of the ejecta will have to be considered. Analysis of each particular SNR might call for some modifications in the calculations or for a more exhaustive exploration of the parameter space, and this task is deferred to future work. However, a preliminary comparison of our results with observations of a well studied SNR is very useful, as it allows us to assess the strengths and limitations of our models.

The best object for this comparison is Tycho, the remnant of SN 1572 (catalog ). Tycho is considered the remnant of a Type Ia supernova, and it has been extensively studied in radio, optical, and X-rays. The X-ray emission from the ejecta has been observed with high angular and spectral resolution by XMM-Newton and Chandra (Decourchelle et al., 2001; Hwang et al., 2002), and spatially integrated but with moderate spectral resolution by ASCA (Hwang, Hughes, and Petre, 1998; Hwang and Gotthelf, 1997, hereafter HG97). However, these extensive observational studies of Tycho resulted in just a few comparisons with SNIa explosion models. Itoh, Masai, and Nomoto (1988), and Brinkmann et al. (1989) calculated the predicted X-ray spectra at the age of Tycho using the deflagration model W7 (Nomoto, Thielemann, and Yokoi, 1984), and compared it to the Tenma and EXOSAT data to find that the observations could not be reproduced by the W7 model unless substantial modifications to its structure where introduced.

We have made no attempt to reproduce the apparent size of the remnant and the observed proper motion of Tycho; we just compared the X-ray spectra of some of our models with the spatially integrated XMM spectrum (Decourchelle et al., 2001, Figure 8). The contribution from the shocked AM has been modeled with the Sedov model in XSPEC (Borkowski, Lyerly, and Reynolds, 2001), with mean shock temperatures and shocked ISM emission measures derived from our hydrodynamical calculations for each case (see Tab. 2), and assuming no collisionless electron heating at the shock ($`kT_\mathrm{e}=0`$), in agreement with optical spectroscopy (Ghavamian et al., 2001). The ionization age, defined as the product of postshock electron density and the remnant’s age, is $`2.93\times 10^{10}`$ s cm-3. We note that the shock speeds we obtain (and therefore the mean shock temperatures $`kT_{\mathrm{s}h}`$ used in the Sedov models) are compatible with X-ray measurements (Hughes et al., 2000), but not with optical or radio observations (Hwang et al., 2002, and references therein). The distance to Tycho was assumed to be 2.3 kpc, and we have used a $`N_\mathrm{H}`$ of $`0.45\times 10^{22}`$ cm-2 as in HG97.

A simple visual comparison of Tycho’s spectrum with the model spectra at t=500 yr (Fig. 7) reveals that one can eliminate the DET model (absence of Si, S, and Ca lines) as well as DEFa, DEFf and DDTe (mainly due to the absence of the FeK$`\alpha `$ line). The contribution from the ISM (Sedov models) to the ejecta spectra varies a lot depending on the type of explosion. It is quite large, specially for the high energy continuum, in models DDTa, DDTe and DET, and quite small in the others, except PDDa which is an intermediate case. For the DEF models in particular, almost all the emission comes from the ejecta, while for DET and DDTe the contribution to the underpredicted lines from the ISM is small, and therefore these models remain unable to account for the spectrum of Tycho. The spectra of the most promising models have been plotted alongside the observed spectrum in Figure 8.

For a more quantitative comparison we have computed the main characteristics of the more prominent emission lines of our models, and compare them to the observational results of HG97 in Tables 2 and 3. The flux for the K$`\alpha `$ complexes of Si, S, Ca and Fe was obtained by adding the contributions from all lines due to atomic transitions from $`n=2`$ to $`n=1`$.

We concentrate now on models DDTa, PDDa, PDDe, and SCH, which have not been eliminated previously. PDDa, PDDe, and SCH give too strong Si and S Ly$`\alpha `$ lines, with respect to the corresponding K$`\alpha `$ line complexes. These three models have in common the formation of an extended high-velocity envelope (corresponding to the accreted He envelope in the SCH model and to the external layers detached from the white dwarf in the pulse previous to the formation of the detonation in the PDDa and PDDe models), which is separated from the rest of the ejecta by an abrupt change in density and followed by dense layers rich in Si and S (see Fig. 1). These dense layers remain dense when shocked, speeding up the ionization processes inside them and producing more highly ionized Si and S and therefore higher Ly$`\alpha `$/K$`\alpha `$ ratios. Other model ratios compare reasonably well with the observed values, most of them within a factor two to three. The Fe K$`\alpha `$ line centroids show deviations larger than 50 eV only for models PDDa and SCH.

Given the limitations of the present work and the considerable uncertainties in the atomic data it is impossible to reach a definitive conclusion, but our preliminary results point to a delayed detonation explosion as the progenitor of Tycho, probably similar to the DDTa model, with a smooth density profile and a large amount of Fe in the outer ejecta. Pure deflagration and pure detonation models can be confidently eliminated. The only features that model DDTa does not reproduce well are the Ca K$`\alpha `$ / Si K$`\alpha `$ ratio and the S K$`\beta `$ / S K$`\alpha `$ ratio, suggesting that the ionization state of S in Tycho may be a little bit higher than in the model. The angular size of Tycho (8’) implies a radius of 2.8($`D`$/2.3 kpc) pc, while the radius of the DDTa model at the age of Tycho is about 3.2 pc (see Fig. 3). If the distance has been estimated correctly, then $`\rho _{\mathrm{A}M}`$ must be higher than $`10^{24}`$ g cm-3. In simulations with a moderately higher value of $`\rho _{\mathrm{A}M}`$, the ionization state of S in the DDTa model increases, bringing the S K$`\beta `$ / S K$`\alpha `$ ratio closer to the observed value while keeping the rest of the S and Si line ratios within reasonable limits. A more detailed study of Tycho including these and other issues will be the subject of future work.

## 5 Conclusions

We have performed self-consistent hydrodynamic and ionization state simulations of the interaction of a grid of SNIa explosion models with a uniform AM, producing a set of synthetic spectra for the integrated thermal X-ray emission from the shocked ejecta in young SNRs. The spectra show remarkable differences between the models: line strengths and strength ratios, overall luminosity, shape of the continuum and other parameters depend strongly on the type of supernova explosion that gave birth to the SNR. These differences stem from the hydrodynamic evolution, that is determined by the density profile of the ejecta synthesized in the explosion, and from the evolution of the electron heating and plasma ionization processes, which depend on both the dynamic history and the chemical composition of the ejecta. A close connection is thus established between supernova explosions and young SNRs, implying that SNRs have the potential to become a new tool to discriminate among various SNIa explosion models, a possibility that should be tested by comparing models with observations. It is also clear that any analysis of a thermal X-ray spectrum from the shocked ejecta in a young SNR, whatever its type, has to make realistic assumptions about the density and chemical composition profile for the ejecta, because using uniform densities or homogeneous compositions is unrealistic and will produce misleading results. The density enhancement effects caused by the hydrodynamic evolution can lead to emission measure ratios very different from elemental ratios (by number) in the material ejected by the explosion. This is because the density enhancement amplifies the contribution to the emitted spectrum by the elements in the outermost layers of ejecta, where the differences between the models are most pronounced. In contrast with the results of Dwarkadas and Chevalier (1998) for the mean gas temperature, we find that both density and electron temperature are highest at the contact discontinuity in all cases.

We have made a preliminary comparison between our theoretical spectra and the observed spectrum of the Tycho SNR. The best choice for Tycho is a delayed detonation model with a large kinetic energy and a high iron content. In particular, our delayed detonation model DDTa provides reasonable agreement with observed line strengths, line ratios, and Fe K$`\alpha `$ line centroid. Explosion models characterized by high density layers in the outer ejecta preceded by an extended low-density envelope, like pulsating delayed detonations and sub-Chandrasekhar models give rise to an excess of ionization of Si and S, which leads to high Ly$`\alpha `$/K$`\alpha `$ line ratios. Deflagration models characterized by a large buffer of C and O in the external layers would produce spectra with too weak lines and too strong continua. One important result is that the emission measure ratios (as well as the line ratios) cannot be safely used as a measure of the relative abundances of various elements in the supernova ejecta; they can in fact differ by as much as two or three orders of magnitude from the actual abundance ratios.

A potentially important issue is related to various hydrodynamical instabilities occurring during the explosion itself or shortly thereafter. Multi-dimensional simulations of Type Ia explosions show the presence of these instabilities (Khokhlov, 2000; Gamezo et al., 2003; Niemeyer et al., 2003; García-Senz and Bravo, 2003; Bravo and García-Senz, 2003), which may lead to inhomogeneities in freely expanding ejecta. There is indeed observational evidence for the presence of fast-moving ejecta clumps in Tycho SNR (Hwang and Gotthelf, 1997). But a quantitative understanding of the clumpiness of the freely expanding ejecta is lacking at this time, both from theoretical and observational perspective (see, for instance, Gamezo et al. 2003, and Thomas et al. 2002). A small degree of clumpiness is not likely to affect our results based on 1-D simulations, but the situation will be different for highly inhomogeneous ejecta. In extreme cases, dense clumps of shocked ejecta might significantly contribute to spatially-integrated X-ray spectra. The problem of clumpy supernova ejecta is clearly outside of the scope of this work, and should be addressed in the future through multidimensional simulations and through detailed observational studies of Type Ia SNR morphologies.

All the models, plots, and spectra shown here are available from the authors upon request.

An important part of this research was conducted during a stay of C.B. at North Carolina State University in the fall of 2002, funded by GENCAT and NASA grant NAG 5-7153. We are grateful to Jordi José for his calculations and Pasquale Mazzotta for providing us with his atomic data and routines. Una Hwang and Anne Decourchelle shared with us helpful discussions. We are very grateful to Anne for providing us with the XMM spectrum for Tycho, and to Martin Laming for his many helpful suggestions about radiative cooling and other issues, and his hospitality at NRL. C.B. and E.B. acknowledge funding from projects AYA2000-1785 and AYA2002-04094-C03-01. K.J.B. acknowledges support by NASA grant NAG 5-7153.

## Appendix A The Supernova Explosion Models

In this appendix we give more details about the SNIa explosion models. As has been mentioned before, the hydrocode, the nucleosynthesis, and the light curve codes, as well as the physics included (equation of state, nuclear reaction rates, etc) are the same as described by Bravo et al. (1996). The parameters of the models are given in Table 1.

The detonation model was obtained from a $`1.38\mathrm{M}_{\mathrm{}}`$ WD in hydrostatic equilibrium (composed of equal masses of 12C and 16O plus a 1% by mass of 22Ne) whose internal energy structure was adjusted to an adiabatic thermal gradient. The ignition was initiated by incinerating the mass in the central layer, and afterwards the detonation propagation was obtained consistently by solving the hydrodynamic and nuclear evolutionary equations. Details of a similar model can be found in Bravo et al. (1996).

The sub-Chandrasekhar mass model was obtained from a WD formed by a C-O core of $`0.8\mathrm{M}_{\mathrm{}}`$ surrounded by a He envelope of $`0.2\mathrm{M}_{\mathrm{}}`$. This envelope was the result of He accretion over the C-O core at a steady rate of $`3.5\times 10^8\mathrm{M}_{\mathrm{}}`$ yr-1. The hydrostatic evolution of the white dwarf subject to accretion was computed by J. José, who kindly provided us with the initial model for the supernova explosion calculation (private communication). In this initial model, the temperature at the base of the He envelope was high enough to induce a spontaneous He detonation. The evolution past this point was followed with the same hydrocode as above.

The rest of the explosion models started from the same initial configuration: an isothermal white dwarf in hydrostatic equilibrium, with the same chemical composition and central density as that used in DET model (the differences in the thermal structures of both configurations account for the small difference in total mass that can be seen in Table 1). Explosion models starting from different central densities $`\rho _c`$ (i.e. WD masses) do not produce substantially different energies or light curves (with the exception of a slight decrease of 56Ni yield with increasing central density due to a larger electron capture rate, Bravo et. al. 1993), so we did not consider variations in $`\rho _\mathrm{c}`$.

The flame propagation velocity in the pure deflagration models was obtained as the maximum between the laminar flame velocity (as given by Timmes and Woosley 1992 and updated by Bravo and García-Senz 1999) and the turbulent velocity, $`v_{\mathrm{R}T}`$. The turbulent velocity was calculated as $`v_{\mathrm{R}T}=\kappa r_{\mathrm{f}l}/\tau _{\mathrm{R}T}`$, where $`r_{\mathrm{f}l}`$ is the flame radius, $`\tau _{\mathrm{R}T}`$ is the local Rayleigh-Taylor time scale at the flame location, and $`\kappa `$ is a parameter given in Table 1 (see Bravo et. al. 1996 for details). Furthermore, the mixing of matter and energy across the flame front was limited as proposed by Wheeler et al. (1987). In this formulation, the mixing is not allowed until the Rayleigh-Taylor front has propagated to a prescribed fraction, $`\theta `$, of the mass of the shell (in all our deflagration models we took $`\theta =0.5`$). Thus, if the mass of the shell ahead of the flame front is $`\mathrm{\Delta }M`$, the mixing was started when the condition $`_{t_0}^t\frac{dM_{\mathrm{R}T}}{dt}𝑑t=\theta \times \mathrm{\Delta }M`$ was held, with $`t_0`$ equal to the time of incineration of the last incinerated shell, and with $`\frac{dM_{\mathrm{R}T}}{dt}=4\pi r_{\mathrm{f}l}^2\times \rho _{\mathrm{f}l}\times v_{\mathrm{R}T}`$. Once this condition was fulfilled, the transfer of internal energy between the Rayleigh-Taylor unstable shells was allowed. Afterwards, the flame propagation was obtained consistently by the consequent increase in the nuclear energy generation rate and its feedback on temperature. At densities below a few times $`10^7`$ g cm-3 the energy generated goes predominantly to create and maintain an electron-positron pair gas rather than to increase the temperature, which is the ultimate cause of the quenching of the flame.

In the delayed detonation models the flame was propagated initially as a deflagration. In that phase, the flame velocity was taken as a constant fraction, $`\iota `$, of the local sound velocity (in our models we took $`\iota =0.03`$). As the flame traveled through lower and lower densities, the sound velocity and hence the deflagration speed decreased, and the flame was eventually quenched. Afterwards, the flame front was artificially accelerated up to a large fraction of the sound speed. Following this fast propagation through a few mass shells (typically 2-3 shells were enough) a detonation formed and propagated through the rest of the star. The location of the transition is univocally determined by the parameter $`\rho _{\mathrm{t}r}`$ (Tab. 1), that is the density ahead of the flame at which the sudden acceleration was imposed. In delayed detonation models as well as in the pure detonation and sub-Chandrasekhar models, Rayleigh-Taylor mixing was not allowed. The algorithm for the flame propagation in pulsating delayed detonation models was the same as in delayed detonation models, the only difference being that in PDD the detonation was not triggered until the WD had pulsated. Further details (in particular the nucleosynthetic output and more complete information about the light curves) are available to interested readers upon request.


format: text

text:
# Resistive MHD simulations of the Parker instability in galactic disks

## 1 Introduction

The issue of stability of magnetized astrophysical disks is of primary importance for their dynamical evolution. Parker (parker1 (1966), parker2 (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 (parker3 (1992)) and developed subsequently in a series of papers by Hanasz & Lesch (hanasz1 (1993), hanasz2 (1998), hanasz3 (2000)) and Hanasz, Otmianowska-Mazur & Lesch (hanasz4 (2002), hereafter referred to as Paper I), and by Brandenburg and Schmitt (brandenburg (1998)) and Moss et al. (moss (1999)).

Another relevant MHD instability is the magneto-rotational instability, referred to as Balbus-Hawley instability (balbus1 (1991), balbus2 (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. fleming (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`$, $`\varphi `$ 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 Section 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., barge (1984), Barghouty and Jokipii barghouty (1996), Zimbardo zimbardo (1995)). Our method is different from the one applied by Fleming et al. (fleming (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 and Hoernes (beck (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 and Tremaine binney (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:

$$\frac{\rho }{t}+\left(\rho v\right)=0,$$
(1)
$$\frac{v}{t}+\left(v\right)v=\frac{1}{\rho }\left(p+\frac{B^2}{8\pi }\right)+\frac{\left(B\right)B}{4\pi \rho }2\mathrm{\Omega }\times v+g,$$
(2)
$$\frac{B}{t}=\times \left(v\times B\eta j\right),$$
(3)
where $`\mathrm{\Omega }`$ is the angular rotation velocity, $`j=\times B`$ is the current density, $`g=g\widehat{e_\mathrm{z}}`$, $`g=\mathrm{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

$$p=c_\mathrm{s}^2\rho ,$$
(4)
where $`c_\mathrm{s}=\mathrm{const}`$ is the isothermal sound speed.

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

$$\eta (j)=\eta _1+\eta _2\left(j^2j_{\mathrm{crit}}^2\right)^{1/2}\mathrm{\Theta }\left(j^2j_{\mathrm{crit}}^2\right),$$
(5)
where $`\eta _1`$ and $`\eta _2`$ are constant coefficients, $`j_{\mathrm{crit}}`$ denotes the critical current density above which the resistivity switches on and $`\mathrm{\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 $`\mathrm{\Omega }`$ which takes a typical galactic value defined in section 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 equations (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 $`(X_i,Z_i)_S`$ of the $`i`$-th line in the S plane are identical as the intersection points $`(X_i,Z_i)_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 $`(X_i,Z_i)_S=(X_i,Z_i)_E`$ at $`t=0`$ then $`(X_i,Z_i)_S=(X_i,Z_i)_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(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 $`(\mathrm{\Delta }X_i,\mathrm{\Delta }Z_i)=(X_i,Z_i)_E(X_i,Z_i)_S`$ (see Fig. 1).

As soon as $`(\mathrm{\Delta }X_i,\mathrm{\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.

$`\overline{\mathrm{\Delta }x}`$ $`=`$ $`{\displaystyle \frac{1}{N}}{\displaystyle \underset{i=1}{\overset{N}{}}}\mathrm{\Delta }x_i,`$ (6)
$`\overline{\mathrm{\Delta }z}`$ $`=`$ $`{\displaystyle \frac{1}{N}}{\displaystyle \underset{i=1}{\overset{N}{}}}\mathrm{\Delta }z_i,`$ (7)
$`\sigma _x`$ $`=`$ $`\left({\displaystyle \frac{1}{N1}}{\displaystyle \underset{i=1}{\overset{N}{}}}(\mathrm{\Delta }x_i\overline{\mathrm{\Delta }x})^2\right)^{1/2},`$ (8)
$`\sigma _z`$ $`=`$ $`\left({\displaystyle \frac{1}{N1}}{\displaystyle \underset{i=1}{\overset{N}{}}}(\mathrm{\Delta }z_i\overline{\mathrm{\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 $`4^{\mathrm{th}}`$ order. This procedure is applied for a large sample of staring points $`\mathrm{P}_\mathrm{S}`$ 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 Section 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 parker1 (1966), Paper I). The initial equilibrium is characterized by an uniform vertical gravity with the vertical gravitational acceleration $`g=210^9\mathrm{cm}\mathrm{s}^1`$, the midplane gas density $`n_0=1\mathrm{cm}^1`$ and the sound speed $`c_\mathrm{s}=7\mathrm{km}\mathrm{s}^1`$, constant across the disk. The fixed ratio of the magnetic pressure to the gas pressure is given by $`\alpha =p_{\mathrm{mag}}/p_{\mathrm{gas}}`$, so that the total pressure is given by $`p_{\mathrm{tot}}=(1+\alpha )p_{\mathrm{gas}}`$. We assume that $`\alpha =1.0`$, corresponding to the magnetic field strength of $`4.5\mu \mathrm{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

$$\frac{n_0(z)}{n_0(0)}=\frac{p_0(z)}{p_0(0)}=\frac{B_0^2(z)}{B_0^2(0)}=\mathrm{exp}\left(\frac{z}{H}\right),$$
(10)
where $`H=\left(1+\alpha \right)c_\mathrm{s}^2/|g|159\mathrm{pc}`$ the vertical scale-height. We apply the angular velocity of the galactic rotation $`\mathrm{\Omega }=0.025\mathrm{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

$$v_\mathrm{z}=\{\begin{array}{ccc}v_0\mathrm{cos}\left(\frac{2\pi n_\mathrm{x}x}{L_\mathrm{x}}\right)\mathrm{cos}\left(\frac{2\pi n_\mathrm{y}y}{L_\mathrm{y}}\right)\mathrm{sin}\left(\frac{\pi z}{H}\right)\hfill & \mathrm{for}|z|H\hfill & \\ 0\hfill & \mathrm{elsewhere}\hfill & \end{array}$$
(11)
where $`n_\mathrm{x}`$, $`n_\mathrm{y}`$ are the numbers of the harmonic components of the X, Y directions, respectively. $`L_\mathrm{x}`$, $`L_\mathrm{y}`$ are the physical sizes of the simulation domain and the initial velocity amplitude is $`v_0=1.0\mathrm{kms}^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_\mathrm{x}=600\mathrm{pc}`$, $`L_\mathrm{y}=1800\mathrm{pc}`$, $`L_\mathrm{z}=1200\mathrm{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_\mathrm{x}\times n_\mathrm{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 $`B_\mathrm{y}`$. 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\mathrm{p}\mathrm{c}z960\mathrm{p}\mathrm{c}`$ 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: 30x90x60 (low) and 60x180x120 (high) as listed in Table 1. We assume $`\mathrm{pc}^2\mathrm{Myr}^1`$ as the unit of $`\eta `$.

### 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:

1. the maximum of vertical velocity $`v_{\mathrm{z}\mathrm{max}}`$ representing the amplitude of unstable Parker modes,
2. the ratio of mean energy densities of total and uniform magnetic fields in the disk ($`zH`$), $`U_{\mathrm{Bt}}/U_{\mathrm{Bu}}=<B^2>/<B>^2`$. The averaging denoted by $`<\mathrm{}>`$ is performed over a volume between $`z=0`$ and $`z=H`$,
3. the maximum resistivity $`\eta _{\mathrm{max}}`$ in the computational volume,
4. the maximum of current density $`j_{\mathrm{max}}`$ in the computational volume,
5. the dispersion $`\sigma _x`$ of magnetic field lines in the X direction,
6. 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.

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 $`\mathrm{\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 $`\mathrm{\Delta }\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. 2(b) and 3(a) 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 Subsection 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.

### 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_{\mathrm{z}\mathrm{max}}`$ presented in Figs. 3(a) and 4(a) demonstrates a common existence of three phases: the initial phase, the exponential growth phase and the saturation one.

During the initial phase, $`v_{\mathrm{z}\mathrm{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_{\mathrm{z}\mathrm{max}}`$ by factor of about 100 (from 0.1 to $`10.0\mathrm{kms}^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÷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 $`v_\mathrm{z}`$ reaches a level of about $`10.0\mathrm{kms}^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_{\mathrm{z}\mathrm{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.

Figs. 3(b) and 4(b) show the evolution of the ratio of mean energy densities of the total and uniform magnetic fields in the disk $`U_{\mathrm{Bt}}/U_{\mathrm{Bu}}`$ 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_{\mathrm{Bt}}/U_{\mathrm{Bu}}`$ 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_{\mathrm{Bt}}/U_{\mathrm{Bu}}`$ 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_{\mathrm{Bt}}/U_{\mathrm{Bu}}`$. 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. 3(b). This means that really high physical resistivity results in that the numerical one is insignificant.

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

The evolution of $`j_{\mathrm{max}}`$ (Fig. 3(d), 4(d) shows similar character as the ratio $`U_{\mathrm{Bt}}/U_{\mathrm{Bu}}`$. 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_{\mathrm{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_{\mathrm{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. 3(d)) 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. 3(e), 3(f) and Figs. 4(e), 4(f) 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 _\mathrm{x}`$ (Figs. 3(e), 4(e)) and Z direction $`\sigma _\mathrm{z}`$ (Figs. 3(f), 4(f)). 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. 2(c). The vertical current sheets lead to the dominating vertical dispersion of reconnected magnetic field lines. The maximum values of $`\sigma _\mathrm{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 _\mathrm{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 _\mathrm{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`$

Figs. 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 _{\mathrm{max}}`$) than the influence of the uniform resistivity.

The time evolution of $`v_{\mathrm{z}\mathrm{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_{\mathrm{crit}}=0.1`$) in the middle of the exponential growth phase. In the saturation phase, the maximum of $`v_\mathrm{z}`$ falls down less than in the uniform resistivity case.

Comparing to the uniform resistivity (Figs. 3(b) and 4(b)) the evolution of $`U_{\mathrm{Bt}}/U_{\mathrm{Bu}}`$ (Figs. 5(b) and 6(b)) 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_{\mathrm{Bt}}/U_{\mathrm{Bu}}`$ 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_{\mathrm{max}}`$ grows to large values, especially for smaller values of $`\eta _1`$ (Figs. 5(d) and 6(d)). For example, for $`\eta _1=0.0`$ and $`\eta _2=0.0`$ $`j_{\mathrm{max}}`$ reaches the value of 0.2 and 0.5 for the low and high resolution, respectively. This difference we discussed in Subsection 3.1. For a large $`\eta _2=100.0`$ the evolution of $`j_{\mathrm{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 _\mathrm{x}`$ and $`\sigma _\mathrm{z}`$ (Figs. 5(e), 5(f) and 6(e), 6(f) 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 $`t900`$ 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 _\mathrm{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 _\mathrm{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.

### 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 Section 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 $`\overline{\sigma _\mathrm{z}}`$ of magnetic field lines defined in Section 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 Figure 7). The computed values of $`\overline{\sigma _\mathrm{z}}`$ are listed together with errors $`\mathrm{\Delta }\overline{\sigma _\mathrm{z}}`$ in Table 2.

For the uniform resistivity models p0, q0, and r0 (see Table 2, low resolution), $`\overline{\sigma _\mathrm{z}}`$ has the same value of about $`150.0\pm \mathrm{\Delta }\overline{\sigma _\mathrm{z}}`$ pc. Independently of the magnitude of explicit resistivity, the time evolution of $`\sigma _\mathrm{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 $`\overline{\sigma _\mathrm{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 $`\overline{\sigma _\mathrm{z}}=A+B\mathrm{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 $`\overline{\sigma _\mathrm{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 $`\overline{\sigma _\mathrm{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 $`\mathrm{\Delta }\overline{\sigma _\mathrm{z}}=16.2`$ in our estimations. By going through the same procedure as described above, but now for values $`\overline{\sigma _\mathrm{z}}+\mathrm{\Delta }\overline{\sigma _\mathrm{z}}=158.8`$ and $`\overline{\sigma _\mathrm{z}}\mathrm{\Delta }\overline{\sigma _\mathrm{z}}=126.4`$, we estimated a possible range for the value of the numerical resistivity, which is $`\eta _{1\mathrm{numerical}}(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 _{1\mathrm{numerical}}`$.

## 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 Subsect. 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 NGC6946 and M51 the magnetic field appears to be much more uniform in the inter-arm regions then in arms (see model of Rohde et al. rohde (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 Subsect. 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:

1. 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.
2. 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.
3. 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.
4. The method can be applied to compare the numerical magnetic diffusivity of different codes and studies of different physical processes e.g. turbulence.

###### Acknowledgements.

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.