TURBULENCE AND MAGNETIC FIELD AMPLIFICATION IN SUPERNOVA REMNANTS: INTERACTIONS BETWEEN A STRONG SHOCK WAVE AND MULTIPHASE INTERSTELLAR MEDIUM

Recent discovery of the year-scale variability in the synchrotron X-ray emission of supernova remnants (SNRs) suggests that the magnetic field should be amplified in the SNR up to the level of milligauss (Uchiyama et al. 2007; Uchiyama & Aharonian 2008, see also Bamba et al. 2003, 2005a, 2005b; Vink & Laming 2003). Since the typical magnetic field strength in the interstellar medium (ISM) is on the order of microgauss, amplification beyond the simple shock compression is necessary to achieve a milligauss level of a magnetic field in a supernova (SN) shell.

The amplification of the magnetic field around the shock wave has been one of the main interests of plasma physics. Recently, Giacalone & Jokipii (2007) have demonstrated by using magnetohydrodynamic (MHD) simulation that density fluctuations in the preshock medium cause turbulence and magnetic field amplification in the postshock medium. The density fluctuations in their preshock state are given by a lognormal probability distribution with a Kolmogorov-like power spectrum, and the maximum magnetic field strength achieved in the postshock medium is larger than a hundred times the preshock field strength. The density fluctuations with the Kolmogorov power spectrum would be accomplished in the preshock medium owing to the transonic nature of turbulence expected in the diffuse ISM (e.g., Hennebelle & Audit 2007). However, applicability of the lognormal distribution of density fluctuations for the ISM is unknown, although their pioneering study is significant. Balsara et al. (2001) have studied the case of an expanding SN blast wave in the turbulent ISM, but have not found the magnetic field amplification beyond 50 μG. Thus, the turbulent field alone does not seem to be able to amplify the magnetic field to the level of milligauss.

It has been known that the ISM is not a uniform but a highly inhomogeneous medium composed of clumpy H i clouds embedded in a diffuse intercloud medium. Observations through the 21 cm absorption lines have shown that H i clouds with column densities N ∼ 10¹⁹ cm⁻², whose spatial scale corresponds to l ∼ 0.1 pc, are ubiquitous in the ISM (see, e.g., Heiles & Troland 2003). The coexistence of the H i clouds and intercloud gas can be understood as follows: in the ISM, the balance between the line-emission coolings and the heating due to ultraviolet background radiation determines two thermally stable equilibrium states with different temperatures that can coexist in pressure balance (Field et al. 1969; Wolfire et al. 1995). One of the equilibrium states corresponds to the diffuse intercloud medium, the so-called warm neutral medium (WNM; n ∼ 0.5 cm⁻³, T ∼ 8000 K), and the other one corresponds to the H i clouds, the so-called cold neutral medium (CNM; n ∼ 50 cm⁻³, T ∼ 100 K). Thus, the coexistence of H i clouds and surrounding diffuse intercloud gas is naturally achieved by the thermally bistable nature of the ISM. The high-resolution MHD simulation performed by Inoue & Inutsuka (2008) has shown that when a weak shock wave (v_shock ∼ 20 km s⁻¹) propagates into the diffuse WNM, it generates thermally unstable gas that evolves into the small-scale cloudlets of CNM embedded in the WNM via the thermal instability (Field 1965). The scale of the CNM cloudlets generated by the thermal instability is typically 0.1 pc that resembles the H i clouds observed in 21 cm absorption lines. Since the ISM is frequently swept up by weak shock waves, for example due to old SN blast waves at a rate on the order of once per millions of years (Cox & Smith 1974; McKee & Ostriker 1977), the typical ISM should be treated as a "two-phase medium" composed of the WNM and CNM rather than a single-phase polytropic gas.

In this paper, we examine the propagation of strong shock waves into the two-phase medium by using two-dimensional MHD simulation, in which the shock strength corresponds to the Sedov–Taylor phase of the SNR. In Section 2, we provide basic assumptions and numerical methods to generate the two-phase medium and to induce shock waves. The results of the simulations are shown in Section 3. Finally, in Section 4, we summarize our results and discuss the implications.

We solve the MHD equations for Cartesian geometry in a conservative fashion:

$$\begin{equation} \frac{\partial \rho }{\partial t}+\vec{\nabla }\cdot (\rho \,\vec{v})=0, \end{equation} \tag{ 1 }$$

$$\begin{equation} \frac{\partial \rho \,\vec{v}}{\partial t}+\vec{\nabla }\cdot \left(p+\frac{B^2}{8\,\pi }+\rho \,\vec{v}\otimes \vec{v}-\frac{1}{4\,\pi }\,\vec{B}\otimes \vec{B}\right)=0, \end{equation} \tag{ 2 }$$

$$\begin{equation} \frac{\partial e}{\partial t}+\vec{\nabla }\cdot \left\lbrace \left(e+p+\frac{B^2}{8\,\pi }\right)\,\vec{v}-\frac{\vec{B}\cdot \vec{v}}{4\,\pi }\vec{B}\right\rbrace =\vec{\nabla }\,\cdot\, \kappa \vec{\nabla } T-L(n,T), \end{equation} \tag{ 3 }$$

$$\begin{equation} \frac{\partial \vec{B}}{\partial t}=\vec{\nabla }\times (\vec{v}\times \vec{B}), \end{equation} \tag{ 4 }$$

$$\begin{equation} e=\frac{p}{\gamma -1}+\frac{\rho \,v^2}{2}+\frac{B^2}{8\,\pi }, \end{equation} \tag{ 5 }$$

where κ is the thermal conductivity and L(n, T) is the net cooling function whose details are noted in the following section. We impose the ideal gas equation of state to close the equations. We use the second-order Godunov scheme (van Leer 1979) for solving MHD equations, in which the hydrodynamic equations with the magnetic pressure term are solved based on the solution of the Riemann problem (Sano et al. 1999), the magnetic tension terms are solved using the method of characteristics (MOC) for Alfvén waves (Stone & Norman 1992), and the induction equations are solved using the consistent MOC with the constrained transport algorithm (Clarke 1996). We use the second-order explicit time integration for the cooling/heating and thermal conduction. The two-dimensional computational domain of x × y = 2² pc² (0 ⩽ x, y ⩽ 2 pc) is used with the uniform 4096² cells (Δx = Δy = 4.9 × 10⁻⁴ pc).

Computations are performed in the following two stages: (1) the generation stage of the two-phase medium by the thermal instability and (2) the propagation stage of the shock wave through the two-phase medium. In the following, we describe the detailed numerical settings of each stage.

2.1. Stage 1: Generation of a Preshock Two-phase Medium

In Stage 1, we take into account the effects of the cooling/heating and thermal conduction. We use the cooling/heating functions given by Koyama & Inutsuka (2002) that are obtained by fitting various line-emission coolings (C ii 158 μm, O i 63 μm, etc.) and photoelectric heating from dusts, which can adequately describe the effects of cooling/heating in the ISM in the range 10 K ≲T ≲ 10⁴ K. The thermal equilibrium state with this cooling/heating is shown in Figure 1. Since this stage deals with a weakly ionized medium, the isotropic thermal conductivity due to the neutral atomic collisions (κ = 2.5 × 10³ erg cm⁻¹ s⁻¹ K⁻¹; Parker 1953) is used.

Zoom In Zoom Out Reset image size

In order to generate the two-phase medium, we initially prepare a uniform gas in thermally unstable equilibrium with n = 2.0 cm⁻³ and p/k_B = 2900 K cm⁻³ (the initial point is plotted as a cross in Figure 1). Flat-spectrum density perturbations are added to the unstable gas to seed the thermal instability. The initial magnetic field strength is uniform with the orientation parallel to the x-axis and the strength of B₀ = 6.0 μG. Note that, as shown in Inoue & Inutsuka (2008), such a condition is naturally produced by the shock compression of the diffuse WNM and the subsequent isochoric radiative cooling. In Figure 1, we illustrate a schematic evolutionary track from the WNM to a thermally unstable medium based on the results of Inoue & Inutsuka (2008). We use the periodic boundary conditions. The calculation of the first stage is stopped at t = 4.0 Myr. Since the growth timescale of the thermal instability is approximately equal to the cooling timescale that is approximately a few Myr in the typical ISM, a typical two-phase medium is formed at t = 4 Myr as a consequence of the thermal instability.

2.2. Stage 2: Injection of a Shock Wave

In Stage 2, we examine the propagation of a shock wave through the two-phase medium generated as a result of Stage 1. We consider the propagation of parallel shock and perpendicular shock. In the case of the parallel shock, we set the hot plasma with n = 0.1 cm⁻³ at the boundary x = 0.0 by which the shock wave is induced. The free boundary condition is imposed at x = 2.0 pc, and we impose the periodic boundary conditions for the boundaries at y = 0.0 and 2.0 pc. In the case of the perpendicular shock, we set the hot plasma at the boundary y = 0.0. The free boundary condition is imposed at y = 2.0 pc, and the periodic boundary conditions are imposed for the boundaries at x = 0.0 and 2.0 pc. We measure the time since the shock is induced.

We study the effect of the shock strength by changing the thermal pressure of the hot plasma from p_th/k_B = 10⁸ to 10⁷ K cm⁻³. These models are summarized in Table 1. For convenience, we list the resulting average propagation speeds of the shocks v_shock in the fourth column of Table 1. The computations are stopped when the shock front reaches the opposite boundary. The effects of the thermal conduction, cooling, and heating that are considered in Stage 1 are omitted in this stage, since the dynamical timescale of this stage (∼1000 yr) is much shorter than the timescales of these effects (∼1 Myr). The hydrodynamic treatment would be reasonable on the scale considered in this paper, since the gyration radius for a thermal proton is estimated as

$$\begin{eqnarray} l_{\rm g} =&& 1.5\times 10^9\,\left(\frac{p/k_{\rm B}}{3\times 10^8\;\rm {K\;cm}^{-3}} \right)^{1/2}\,\nonumber\\ &&\times\left(\frac{n}{10\;\rm { cm}^{-3}} \right)^{-1/2}\,\left(\frac{B}{6\;\mu \rm {G}} \right)^{-1}\;\rm {cm}, \end{eqnarray} \tag{ 6 }$$

which is much smaller than the spatial resolution of our numerical simulations and justifies (magneto-)hydrodynamics approximation.

Table 1. Model Parameters

Model	pth/kB	Shock Type	Shock Speed
1	3.0 × 108 K cm−3	Perpendicular shock	1256 km s−1
2	3.0 × 108 K cm−3	Parallel shock	1289 km s−1
3	1.0 × 108 K cm−3	Parallel shock	726 km s−1
4	3.0 × 107 K cm−3	Parallel shock	397 km s−1

Download table as:  ASCIITypeset image

3.1. Results of Stage 1: Two-phase Medium

The resulting density structure of Stage 1 is shown in Figure 2. The cold and dense filamentary clumps (CNM) and the surrounding warm diffuse gas (WNM) are formed as a consequence of the thermal instability. The condensation driven by the thermal instability to form the CNM essentially arises along the magnetic field, since the motion perpendicular to the field line is easily stopped owing to the enhancement of the magnetic pressure. Thus, the thickness of the filaments can be described by the most unstable scale of the thermal instability (∼1 pc) times the compression ratio of the condensation that gives ∼0.1 pc, and the length of the filament is roughly given by the most unstable scale of the thermal instability.

Zoom In Zoom Out Reset image size

We stress that such a two-phase structure is quite natural and ubiquitously expected in the ISM not only from the theoretical but also from the observational point of view (Heiles & Troland 2003). A more detailed description of the evolution of the thermal instability can be found, for example, in Inoue et al. (2007), Hennebelle & Audit (2007), and Inoue & Inutsuka (2008).

3.2. Results of Stage 2: Effects of the Shock Wave

3.2.1. Turbulence in the Shell

In Figures 3 and 4, we show the results of Model 1 at t = 1425 yr and Model 2 at t = 1508 yr, respectively. The top and bottom panels, respectively, represent the structure of the number density and magnetic field strength. In the following, we use data at these snapshots for the analyses of Models 1 and 2.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In both Models 1 and 2, the Alfvén Mach number M_A and the sonic Mach number M_s in the WNM and CNM can be roughly estimated as

$$\begin{equation} M_{\rm A, WNM}\sim \frac{v_{\rm shock}}{B_0/\sqrt{4\,\pi \,\rho _{\rm WNM}}}\simeq 110, \end{equation} \tag{ 7 }$$

$$\begin{equation} M_{\rm s, WNM}\sim v_{\rm shock}/c_{\rm s, WNM}\simeq 170, \end{equation} \tag{ 8 }$$

$$\begin{equation} M_{\rm A, CNM}\sim \frac{v_{\rm shock}}{B_0/\sqrt{4\,\pi \,\rho _{\rm CNM}}}\simeq 760, \end{equation} \tag{ 9 }$$

$$\begin{equation} M_{\rm s, CNM}\sim v_{\rm shock}/c_{\rm s, CNM}\simeq 1200, \end{equation} \tag{ 10 }$$

where we have used the conditions n_CNM = 50 cm⁻³, T_CNM = 100 K, n_WNM = 1 cm⁻³, and T_CNM = 5000 K in the evaluation.

The density structures in both models indicate that the shocked shells that are formed by piling up the two-phase medium are turbulent. As discussed in Giacalone & Jokipii (2007), the density bumps in the preshock region lead to the rippling of the shock front that consequently generates vortexes in the postshock flow. This type of vorticity generation is known as Richtmyer–Meshkov instability that is a kind of the Rayleigh–Taylor type instability (e.g., Brouillette 2002). Many mushroom-shaped structures that indicate the nonlinear growth of the Rayleigh–Taylor type instability are found in the shell (see the top panels of Figures 3 and 4).

In order to adequately calculate the rippling of the shock front that determines the vorticity of the postshock flow, we have to resolve the transition layers between the WNM and the CNM clumps at which the curvature of the rippling is the largest and the generated vorticity is the strongest. It is known that the scale of the transition layer is determined by the so-called "Filed length" $l_{\rm F}=\sqrt{\kappa \,T/L(n,T)}$, where L(n, T) is the cooling rate per unit volume (Field 1965; Begelman & McKee 1990). In the typical ISM, the filed length takes the value of l_F ∼ 0.01–0.1 pc (Inoue et al. 2006), which is much larger than the resolution of our simulations. Thus, our simulations can describe the postshock turbulence accurately.

The velocity dispersions in the shocked shells (δv = 500 km s⁻¹ for Model 1 and 450 km s⁻¹ for Model 2) are comparable to the sound speed in the shell:

$$\begin{equation} c_{\rm s} = 570\,\left(\frac{p/k_{\rm B}}{ 3\times 10^{8}\;{\rm K\; cm}^{-3}} \right)^{1/2}\left(\frac{n}{10\;{\rm cm}^{-3}} \right)^{-1/2}\;{\rm km\;s}^{-1}, \end{equation} \tag{ 11 }$$

where the velocity dispersions are estimated by the half-widths of the velocity distribution in the shell. Throughout this paper, we define the region in "the shell" such that p/k_B > 10⁵ K cm⁻³ and n > 1 cm⁻³ in the region. The former condition excludes the preshock region and the latter condition excludes the region filled by hot plasma with n ≃ 0.1 cm⁻³.

3.2.2. Magnetic Field Amplification

The magnetic field strength in the shell has huge amplitude fluctuations as depicted in the bottom panels of Figures 3 and 4. There exist many filamentary regions with the field strength far beyond the value expected only from the shock compression. The maximum field strengths achieved during the evolutions are 1280 μG in Model 1 and 860 μG in Model 2.

The amplifications of the magnetic fields would come from the turbulence. Since the preexisting magnetic field is much weaker than the postshock turbulence, the turbulent velocity field can easily stretch and deform the magnetic field lines, which creates the regions with the amplified magnetic field (Giacalone & Jokipii 2007). One can further understand the mechanism of the strong field amplification by the following equation that is obtained from the equation of continuity and the induction equation:

$$\begin{equation} \frac{d}{dt}\left(\frac{\vec{B}}{\rho } \right)=\frac{1}{\rho }\,(\vec{B}\cdot \vec{\nabla })\,\vec{v}. \end{equation} \tag{ 12 }$$

This equation shows that the magnetic field can be amplified, if the velocity has shear along the field line. Such a situation is realized in the postshock shell due to the turbulence. Especially at the transition layer between CNM and WNM, the velocity shear (vortex) is induced most strongly, since the propagation velocity of the shock wave is very different in WNM and CNM. This picture of the magnetic field amplification is reinforced by the following properties of the shocked shell. First, the regions where the magnetic field is strongly amplified roughly trace the transition layers between shocked CNM clumps and WNM. Figure 5 shows the closeup views of the magnetic field strength structure of Models 1 (top) and 2 (bottom). One can see that the thickness of the regions with the strong magnetic field is approximately 0.01 pc, which well agrees with the thickness of the transition layer (Inoue et al. 2006). Second, in the region where the strength of the magnetic field is maximum, the plasma β is on the order of unity (see Section 3.2.4). This reflects the fact that the back-reaction of the magnetic field stops amplification due to turbulent flows whose velocity dispersion is comparable to the sound speed.

Zoom In Zoom Out Reset image size

In the top panel of Figure 9, we plot time evolutions of maximum (solid) and average (dotted) field strength in the shocked shell, where the average field strength is defined by $\langle |\vec{B}| \rangle _{\rm shell}$. The red and green lines, respectively, represent those of Models 1 (perpendicular shock) and 2 (parallel shock). The maximum field strengths grow fast and saturate at ∼1 mG in both Models 1 and 2. The average field strength also exceeds the value expected only from the shock compression. Does the growth of average field strength come only from the strong amplification at the narrow transition layers? Our answer is no, since the filling factor of the region where the magnetic field is amplified to the level of 1 mG is small. Furthermore, as seen in Figures 3 and 4, the magnetic fields are amplified to tens or hundreds of microgauss broadly in the downstream. This growth of the average field strength may be called "turbulent dynamo" (Batchelor 1950; Cho & Vishniac 2000; Cho & Lazarian 2003) that leads the magnetic field amplification through the stretching of the field line by turbulent eddies (vortices). Of course, the origin of the eddy is the rippling of the shock front that is due to the density inhomogeneity of the multiphase preshock gas. We discuss more details about the turbulent dynamo in the next section by taking spectra of the turbulent postshock medium.

The difference between the average field strengths between Models 1 and 2 is larger than the difference between the maximum field strengths. This is simply due to the difference in initial field amplification by the simple shock compression. The level of average field strengths is roughly an order of magnitude smaller than the maximum. Therefore, in some regions where the magnetic field is amplified to 1 mG, about a half of preshock kinetic energy is converted to magnetic energy, but on average the conversion factor is only 1% or less depending on the angle of the global shock normal and the magnetic field.

In Figure 6, we plot the volume fractions of the regions where magnetic field strength is larger than B in the shell at the times depicted in Figures 3 and 4. The volume fractions where |B| ⩾ 100 μG and 500 μG are, respectively, 27.6% and 0.5% in Model 1 and 9.6% and 0.03% in Model 2.

Zoom In Zoom Out Reset image size

3.2.3. Spatial Profiles and Power Spectra

In Figure 7, we show cross section profiles of Model 1 along the y-axis at x = 1.0 pc (solid) and of Model 2 along the x-axis at y = 1.0 pc (dotted). The number density, magnetic field strength, velocity component normal to the shock (v_y for Model 1 and v_x for Model 2), and velocity component tangential to the shock (v_x for Model 1 and v_y for Model 2) are plotted from the top to bottom panel. The shock front is located around the position y = 1.85 pc (Model 1) and x = 1.90 pc (Model 2). The left and right sides of the shock front are downstream and upstream, respectively. From the density cross section, we can see many clumps in the shell with well defined boundaries. Most of the clumps are the shocked CNM and some of the density jumps are generated by secondary shocks in the turbulence. The velocity cross sections also show large amplitude fluctuations, in which discontinuous jumps due to the secondary shocks can be seen.

Zoom In Zoom Out Reset image size

In order to understand the statistical properties of the turbulent fluctuations in the shocked shell, it is helpful to observe their power spectra. Figure 8 shows the power spectra of the velocity v²_k, magnetic field B²_k/4π, and $(\sqrt{\rho }v)_k^2$ obtained from the data of Model 1 in the region x ∈ [0.75, 1.25] and y ∈ [1.3, 1.8] (top), and Model 2 in the region x ∈ [1.3, 1.8] and y ∈ [0.75, 1.25] (bottom). Note that the shell integrals in the k space of B²_k/8π and $(\sqrt{\rho }v)_k^2/2$ give, respectively, the magnetic energy and the kinetic energy. Thus, the ratio of B²_k/4π and $(\sqrt{\rho }v)_k^2$ indicates the ratio of the magnetic to kinetic energies at scale k.

Zoom In Zoom Out Reset image size

The spectra of two models show qualitatively similar profiles. This is due to the fact that the turbulence in the shells is super-Alfvénic in both cases. In the case of the super-Alfvénic turbulence, as mentioned in the previous section, vortices (eddies) cause stretching and amplify magnetic field not only at the transition layers but also in a diffuse region (turbulent dynamo). In the following, we compare the spectra taken from our results and those of other simulations in which the driven turbulent dynamo is studied.

In the large scales (k/2π ≲ 100 pc⁻¹), the velocity power spectra show the Kolmogorov spectrum v²_k ∝ k^{−(5/3)−(D−1)} (D = 3 in a three-dimensional case and D = 2 in a two-dimensional case), and the spectra of magnetic fields are flatter than those of velocities. These properties are also indicated from the simulations of the simple, one-phase, super-Alfvénic-driven turbulence (see, e.g., Cho & Lazarian 2003). The coincidence of the large-scale characteristics of v²_k and B²_k between our two-phase results and other one-phase results is not surprising, since the character of the two-phase medium arises only on small scales owing to the smallness of the CNM clumps. The only difference between the two models on large scales is that the power spectrum of the magnetic field in Model 1 has larger amplitude than that of Model 2. This can be also seen in the bottom panels of Figures 3 and 4, that is, the magnetic field in Model 1 has larger structures than in Model 2. In Model 1, the timescale of the interaction between the shock wave and CNM clumps is longer than that of Model 2 owing to the filamentary structure of the CNM clumps that would cause larger eddies and thus larger structures of the magnetic field. In contrast to the velocity, the power spectra of $\sqrt{\rho }v$ are slightly flatter than the Kolmogorov spectrum. This is due to the small-scale structures of the shocked CNM clumps, since the existence of the delta-function-like structures of the density makes the power spectra flatter than the Kolmogorov spectrum (e.g., Hennebelle & Audit 2007).

On small scales (k/2π ≳ 100 pc⁻¹), the velocity spectra decline more steeply than Kolmogorov. The level of the magnetic power B²_k and the kinetic power $(\sqrt{\rho }v)_k^2$ are closer than on the large scale. In the case of the transonic, super-Alfvénic turbulence (Cho & Lazarian 2003) and the incompressible, super-Alfvénic turbulence (Kida et al. 1991; Cho & Vishniac 2000), the power of the magnetic field is reported to be larger than that of the velocity field on small scales. Note that in the case of one-phase turbulence, the ratio of the velocity to magnetic field power spectra gives the ratio of the kinetic to magnetic energy spectra, since the density is normalized to unity. In our case, however, the magnetic field power spectrum does not dominate the kinetic power spectrum even on small scales, although the difference is small. A result similar to ours was also reported in Padoan & Nordlund (1999) in which highly supersonic, super-Alfvénic turbulence was studied. Thus, whether or not the magnetic power dominates the kinetic power on small scales seems to depend on the property of dynamics.

So far, we have not discussed time evolutions of the spectra. In the case of our simulation, it is difficult to analyze the detailed time dependence of the power spectra, since a shocked shell is too thin to take the spectra at an early stage. In the case of the results of Cho & Vishniac (2000) and Cho & Lazarian (2003), until the growth of an average magnetic field saturates, the power spectrum of the magnetic field shows an inverse cascade, that is, the scale of the magnetic field becomes larger as the field becomes stronger. In our case, the average field strength does not show clear saturation, though the growth rate is small compared to the earlier stage (see the top panel of Figure 9). Thus, the average field strength has some possibility to grow further and the scale at which the spectra change their properties (k ∼ 100 at the moment taking the spectra in Figure 8) possibly becomes larger. Thus, it remains imperative to do long-term calculations for studying the time evolution of the turbulence spectra.

Zoom In Zoom Out Reset image size

3.2.4. Dependence on Shock Strength

In Figure 9, we show the results of the parameter study. The top panel represents the evolutions of the maximum (thick lines) and average (thin lines) field strengths, and the bottom panel represents the local plasma βs at the point where the magnetic field strength is maximum. The red, green, blue, and black lines are the results of Models 1, 2, 3, and 4, respectively. The top panel shows that the magnetic field strengths decrease with increasing thermal pressure of the hot plasma or decreasing strength of the shock waves. As indicated in Section 3.2.1, in the case of Models 1 and 2, the velocity dispersion of the turbulence is on the order of sound speed of the postshock medium. In the case of Models 3 and 4, the velocity dispersions in the shocked shell are respectively 240 km s⁻¹ and 130 km s⁻¹, which are on the order of the sound speeds in each shocked shell (see Equation (11)). If one assumes that, as discussed in Section 3.2.2, the maximum magnetic field strength is determined by the strength of the turbulent flow that is comparable to the sound speed, the local plasma β where the magnetic field strength is maximum should be on the order of unity. The bottom panel clearly verifies this expectation.

The velocity and magnetic field power spectra do not show any qualitative difference between Models 1 and 2. The probability distribution functions (PDFs) of the magnetic field strength in the shocked shell (Figure 10) indicate that there exists a range in which the PDFs show the power-law dependence of the field strength (P[B] ∝ B^−p) with the index p ∼ 1. This power-law range is extended to higher field strength with increasing shock speed and increasing angle of the preshock magnetic field direction to the normal of the plane of the global shock wave.

Zoom In Zoom Out Reset image size

We have examined the MHD simulations of the shocked two-phase medium in two-dimensional geometry. We have shown that the shocked shell becomes turbulent owing to the preshock density inhomogeneity and the magnetic field amplification taking place in the shell. This is consistent with the work of Giacalone & Jokipii (2007). In the cases of Models 1 and 2, in which the propagation speeds of the shocks correspond to the Sedov–Taylor stage of the SN blast wave, the magnetic field strength is amplified up to the order of 1 mG that is roughly determined by the condition β ∼ 1.

The preshock density inhomogeneity in our model originates in the thermal instability that is naturally and ubiquitously expected in the ISM (Inoue & Inutsuka 2008). Thus, the process of the magnetic field amplification investigated in this paper is also expected in real SNRs. Recent X-ray observations of SNRs by Uchiyama et al. (2007) and Uchiyama & Aharonian (2008) discovered synchrotron X-ray hot spots or filaments. The magnetic fields are amplified up to ∼1 mG there, if the flux-decreasing timescale of the X-ray hot spots, which is on the order of a few years, is determined by the synchrotron cooling timescale of accelerated electrons, $t_{\rm synch}\simeq 1.5\,(B/\rm {mG})^{-1.5}\,(\epsilon /\rm {keV})^{-0/5}$ year, where is the energy of an emitted photon. In the following, we compare our simulation results with the observed characteristics of the X-ray hot spots.

• 1.  
  Spatial scale. The typical scale length of the observed X-ray hot spots is ℓ ∼ 0.05 pc, which is comparable to the scale of the regions where the magnetic field strength is amplified to the order of 1 mG in our simulation (see Figure 5, where the regions depicted in red have field strengths larger than 0.6 mG). This scale length is determined by the width of the transition layers between the CNM clumps and surrounding WNM (∼0.1–0.01 pc) at which the strongest shear is induced and amplifies the magnetic field.
• 2.  
  Location of the X-ray hot spot. The X-ray hot spots observed seem to be located at more than 0.1 pc behind the shock front, although the projection effect should be taken into account. This fact can be naturally explained in our scenario (see the bottom panels of Figures 3, 4, and 7). However, it may be difficult for high-energy electrons, which accelerated at the shock front and emitting synchrotron X-rays, to move from the shock front to the emitting region (hot spot) because of the synchrotron cooling. Moreover, the crossing time of the hot spots for such electrons, ℓ/v where v is the flow velocity at the hot spot, is much longer than the observed variability timescale of a few years. These facts imply that the electrons responsible for the X-ray hot spots are not accelerated at the shock front, but accelerated in situ or at least very nearby compared to the shock front. A possible acceleration site is the secondary shocks arising from the turbulent flow (e.g., collisions of turbulent flows) in the shocked shell. As indicated in Section 3.2.3, our simulations show the generation of many secondary shocks in the shocked shell (see, Figure 7).
• 3.  
  Variability timescale and electron acceleration. The maximum magnetic field strength in our simulation is large enough to explain the flux decreasing timescale of a few years. However, it is still uncertain whether the secondary shock that sweeps the region with a strong magnetic field can accelerate electrons enough to make the hot spot. In order to confirm this, further studies are necessary taking into account the acceleration of electrons at the secondary shocks.

The filling factor of the regions with a high magnetic field (B ∼ 1 mG) that correspond to the hot spots is less than 1% in our results. However, we emphasize here that it depends on the number of the CNM clumps embedded in the diffuse WNM. If we consider the region of ISM including a much larger number of CNM clumps, a shock wave that piles up such ISM to form an SNR is expected to show more hot spots. In particular, huge amounts of H i CNM clumps are expected around molecular clouds (see, e.g., Blitz et al. 2007 for observation and Koyama & Inutsuka 2002; Inoue & Inutsuka 2008; Hennebelle et al. 2008; Heitsch & Hartmann 2008 for simulations of molecular cloud formation). For example, SNR RXJ1713 is suggested to be interacting with molecular clouds (Fukui et al. 2008) and, thus, may show many hot spots than our simulations.

The evidence of magnetic field amplification in an nSNR has also been obtained from the thickness of X-ray filaments that typically shows B ∼ 100 μG (Bamba et al. 2003, 2005a, 2005b; Vink & Laming 2003). Since the X-ray filaments have apparently coherent features on parsec scales, it seems to be refracting the "average" field strength in the SNR. The average field strength in our simulation shows tens to a hundred microgauss (see Figure 9) that possibly explains the X-ray filaments. Note that a thickness of ∼ 0.01 pc of the X-ray filaments would indicate the spatial distribution of accelerated electrons. At this moment, this is not described in our simulations omitting nonthermal effects. Obviously, the nonthermal effects (e.g., Lucek & Bell 2000) should also be important in understanding the origin of the X-ray filaments.

Our simulations predict the spatial power spectrum of magnetic field distribution in the SNR shell that is slightly shallower than the Kolmogorov spectrum on the scales k < 100 pc⁻¹ and slightly steeper than the Kolmogorov on the scales k > 100 pc ⁻¹ as is the case of the super-Alfvénic turbulence (see Figure 8). Such fluctuations of the magnetic field should be reflected as a spatial inhomogeneity of the synchrotron radiations as the synchrotron X-ray images of SNRs depicted in Uchiyama et al. (2007) and Uchiyama & Aharonian (2008). Although the inhomogeneity of the radiation reflects the inhomogeneous distribution of accelerated particles, it may not be easy to directly measure the spatial power spectrum of the magnetic field from observations.

Magnetic field strength on the order of milligauss is also indicated in the center of the Galaxy where spectacular filamentary structures are observed with polarized radio emission. Although the mechanism of creation of a large-scale filamentary structure with nonthermal emission remains to be found, our result provides a hint of the origin of the ∼milligauss field strength because the Galactic center is expected to be occupied by many SNRs and, hence, should be the site for the magnetic field amplification studied in this paper. Spreading out of the locally amplified magnetic field into a large-scale volume is also expected in starburst galaxies that inevitably experience concurrent SN explosions. Thus, it might be interesting to study the long-term evolution of the strongly magnetized turbulent ISM with the effects of multiple strong shock waves in such environments.

In this paper, we have performed simulations in two-dimensional geometry in order to keep the spatial resolution. In principle, it is imperative to perform similar high-resolution simulations in three-dimensional geometry to analyze the realistic three-dimensional phenomena. We expect, however, that the possible difference is not significant, since the statistical properties of the turbulence in our simulation can be understood consistently with the three-dimensional simulations of super-Alfvénic turbulence (see Section 3.2.3). We also neglect the effect of explicit magnetic diffusivity in the calculation of Stage 2, since the resistivity of the collisionless plasma in SNRs is highly unknown. It is known that the strength of magnetic field amplification by turbulent dynamo depends on the resistivity or magnetic Reynolds number R_m, which determine the efficiency of magnetic field dissipation mainly via reconnection, that is, larger resistivity (smaller Rm) leads weaker amplification of the magnetic field. However we expect that the saturation level of magnetic field strength does not change, in the cases where R_m is roughly lager than 100 (see Figure 2 of Cho & Vishniac 2000). In our simulation, the reconnection of the magnetic field takes place at the scale of a numerical grid. The magnetic Reynolds number R_m can be estimated to the order of the cell number in the numerical domain that is ∼1000 in our simulation. Therefore, we think that our results should not change unless the resistivity in SNRs is huge to leads R_m < 100.

T.I. thanks T. K. Suzuki for the helpful discussion of MHD turbulence and numerical scheme. Numerical computations were carried out on NEC SX-9 and XT4 at Center for Computational Astrophysics (CfCA) of National Astronomical Observatory of Japan, and NEC SX-8 at YITP in Kyoto University. This work is supported in part by a grant-in-aid from the Ministry of Education, Culture, Sports, Science, and Technology (MEXT) of Japan, No. 18740153 and No. 19047004 (R.Y), and No. 15740118, No. 16077202, and No. 18540238 (S.I.).