Submitted to Astrophysical Journal

Long-Term Evolution of Supernova Remnants in MagnetizedInterstellar Medium
Hidekazu Hanayama1

arXiv:astro-ph/0507421 v1 19 Jul 2005

Department of Astronomy, Graduate School of Science, University of Tokyo, Tokyo
113-0033, Japan
Kohji Tomisaka2
Theoretical Astrophysics, National Astronomical Observatory, Mitaka, Tokyo 181-8588,
Japan
ABSTRACTThe evolution of supernova remnants (SNRs) is studied, with particular attention to the effect of magnetic fields with axisymmetric two-dimensional magnetohydrodynamical simulations. The evolution of magnetic SNRs is the sameas non-magnetic ones in the adiabatic Sedov stage. After a thin shell is formed,
the shell is driven by the pressure of the hot interior gas (bubble). Evolution in
the pressure-driven snow-plow phase is much affected by the magnetic field. The
shell sweeping the magnetic field lines thickens owing to the magnetic pressure
force. After 5 × 105 yr - 2 × 106 yr, the inner boundary of the thick shell begins
to contract. This compresses the hot bubble radially and maintains its thermalpressure. Thus, the bubble forms a prolate spheroidal shape and becomes thinner and thinner, since it expands in a direction parallel to the magnetic field forB0 >∼ 3µG. Finally, the bubble contracts. The porosity of the hot low-density
gas in ISM is reduced, taking the effect of the magnetic field into account.
Subject headings: ISM: bubbles — ISM: magnetic fields — ISM: supernova remnants
12hanayama@th.nao.ac.jp

also School of Physical Sciences, The Graduate University for Advanced Studies (SOKENDAI);
tomisaka@th.nao.ac.jp

–2–
1.

Introduction

Structure of the interstellar medium (ISM) is a key issue in the understanding of our
Galaxy. A large fraction of the interstellar space is occupied with warm (T ∼ 104 K) and5
hot (T >∼ 3 × 10 K) media. The origin of the hot component is the energy released fromsupernovae and/or stellar winds. Shock waves of the supernova remnants (SNRs) hit theinterstellar gas to heat it up and to sweep it into a dense shell. The evolution of SNRs
has been studied for many years (for example, Cox 1972; Chevalier 1974; Cioffi et al. 1988).
After the shock expansion speed decelerates below ≃ 260n0.08 km s−1 , the gas passed throughthe shock front is quickly cooled and a dense shell forms. The expansion of the shell is
decelerated as the interstellar gas accumulates. Finally, the expansion speed of the shock
front reaches the random speed of turbulent motion in the ISM and the shock wave decays
to a sound wave.
The interstellar magnetic field is seen at a glance to run in a direction parallel tothe galactic disk. As for the strength, combining rotation and dispersion measures, Rand
& Kulkarni (1989) obtained a turbulent field strength of ≃ 5µG and its cell size of ≃55
pc. The energy density of the magnetic field B02 /8π ≃ 10−12 dyn cm−2 (B0 /5µG)2 attains
that of the ISM turbulent motion vt with a density ρ0 as ρ0 vt2 /2 ≃ 10−12 dyn cm−2 (ρ0 /2 ×
10−24 g cm−3 )(vt /10km s−1 )2 and an interstellar thermal energy p = kT ρ/µH ≃ 3 × 10−12
dyn cm−2 (T /104 K)(ρ0 /2 × 10−24 g cm−3 ). These three energies are in approximate equipartition. This seems to indicate the magnetic field plays an important role in the ISM. At thesame time, the late evolutionary phase of SNRs must also be significantly affected by theinterstellar magnetic field.
Magnetic SNRs were first studied by Chevalier (1974) using spherically symmetric onedimensional hydro simulations including only the magnetic pressure force. Slavin & Cox(1992) pointed out the effects of the magnetic field as (1) an SNR forms a thick shell supported by the magnetic pressure force; (2) the hot gas heated by the shock wave in the
adiabatic phase occupies the interior of the SNR. This hot cavity is effectively re-compressedin the later phase by the effect of the magnetic pressure from their spherically symmetric calculations. Similar one-dimensional simulations have also been done by Smith & Cox (2001).
The tension force of the magnetic field has similar effects, broadening the shell and compressing the hot bubble, but is not included in these spherically symmetric simulations. At
least two-dimensional magnetohydrodynamical (MHD) simulations are necessary. Balsara
et al. (2001) have studied the evolution of SNRs expanding in the magnetic turbulent ISM
using three-dimensional MHD simulations and they have explained the evolution of SNRs
in the thermal X-ray and non-thermal radio continuum. However, their study is restricted
to the adiabatic stage. In the present paper, we explore the total evolution of SNRs in the

–3–
magnetized ISM using axisymmetric MHD simulations.
The late phase evolution of SNRs is essential to understanding of the structure of theISM (Koo & Kang 2004). For example, hot ionized gas is believed to be confined in thecavities of SNRs, especially in large-volume old SNRs. However, the volumes of the hotcavities must be strongly affected by the interstellar magnetic field strength. We devote
particular attention to the late-phase evolution of SNRs.
The plan of the paper is as follows: in §2, model and numerical method are described.
Section 3 shows numerical result. We pay special attention to the comparison between
models with strong and weak magnetic field. It is shown that the evolution of SNRs is also
affected by the ISM density and local heating rates. Expansion laws of the SNR shell and hot
cavity are also shown in §3. Section 4 is devoted to discussion on how the volume fractionof the hot ISM is affected by the magnetic field strength. The reason why the shell densityof SNR is determined is also discussed in §4.

2.

Model and Numerical Method

We assume that a uniform interstellar gas with density n0 (cm−3 ) and magnetic fieldstrength B0 (µG) is extending. The half-width at half-maximum of the thick Gaussian component and the scale-height of the exponential component of the HI disk are equal to 265
pc and 403 pc, respectively (Dickey & Lockman 1990). Since the size of the SNR is smaller
than the scale-height of the HI, HHI , in this paper, we ignore the effect of the density stratification in the z-direction, which becomes important after the size of the bubble exceeds
several times HHI (Tomisaka 1998). The scale-height of the magnetic field is believed to bemuch larger than that of the gas. We assume here a uniform density and a uniform magnetic
field strength.
The basic equations are the ideal MHD equations. Using the cylindrical coordinates(z, r, φ), assuming that we have no toroidal components of velocity and magnetic field vφ =Bφ = 0 and that variables are independent from the azimuth angle φ as ∂/∂φ = 0, the basic
equations are as follows:
∂ρ
∂
1 ∂
+
(ρvz ) +(rρvr ) = 0,
(1)
∂t ∂z
r ∂r


 

Br 2 − Bz 2
1 ∂
Br Bz
∂
∂ρvz
2
ρvz + p ++r ρvz vr −
=0
(2)
+
∂t
∂z
8π
r ∂r
4π
∂ρvr
∂t

∂
+
∂z



Br Bz
ρvr vz −
4π

–4–
 



Bz 2 − Br 2
1
Br 2 + Bz 2
1 ∂2
r ρvr + p +
=
p+
+
r ∂r
8π
r
8π


∂
Br (vz Br − vr Bz )B2∂
(ǫ + p)vz +
(ǫ +
) +
∂t
8π
∂z
4π
 

1 ∂
Bz (vz Br − vr Bz )+r (ǫ + p)vr −= Γ0 nH − Λ(T )nH 2
r ∂r
4π
1 ∂
∂Bz
=
{r(vz Br − vr Bz )},
∂t
r ∂r

(3)

(4)

(5)

∂Br
∂
= − (vz Br − vr Bz ),
(6)
∂t
∂z
where the variables have their ordinary meanings as ρ (density), v ≡ (vz , vr ) (velocity),p (pressure), T (temperature), B ≡ (Bz , Br ) (magnetic flux density), and ǫ ≡ p/(γ −1) + ρv 2 /2 (internal plus kinetic energy per unit volume). Equation (1) is the continuity
equation; equations (2) and (3) are the equations of motion; equation (4) represents theenergy conservation equation for internal plus kinetic energy. The induction equations forthe poloidal magnetic fields are equations (5) and (6). The terms Γ0 nH and Λ(T )nH 2 inequation (4) represent, respectively, the heating and radiative cooling rates, where nH means
the number density of the Hydrogen atoms. We use the radiative cooling rate estimated forthe gas in the collisional ionization equilibrium with the solar metallicity (Raymond, Cox,
& Smith 1976).

2.1.

Numerical Method

We employ here the modified Lax-Wendroff scheme (Rubin & Burstein 1967), an explicit
finite-difference method, to solve the basic equations (1)-(6). The number of grids and grid
spacings are chosen as (Nz , Nr ) = (2500, 2500) and ∆z = ∆r = 0.1 pc. Thus, our numerical
domain covers a region of 250 pc × 250 pc. Artificial viscosity is added to smooth the jumpat the shock front. For all the quantities U ≡ (ρ, ρv, ǫ + B 2 /8π, B), we added an extra
term to express artificial diffusion and solve equations like∂U
= ∇(c∇U),(7)∂twhere we take
"
2 
2 
2 
2 #1/2
∂vz
∂vz
∂vz
∂vr
c = Av
+
+
+
∆2 ,
(8)
∂z
∂r
∂z
∂r

–5–
with ∆ = ∆z = ∆r and Av = 4. The numerical scheme was tested by comparing known
solutions which have been analytically and numerically obtained. We checked the program
by comparing (1) the adiabatic SNR with the Sedov solution (Sedov 1959), (2) the sphericallysymmetric SNR with previously calculated (Slavin & Cox 1992) and (3) MHD shock tube
problems (Brio & Wu 1988). The analytical solution of (1) is reproduced within a 5% relative
error, for example, pressure and density distributions. The expansion law of the shock frontof hydrodynamical SNRs of (2) coincides with previous work (such as that by Slavin & Cox(1992)) and the difference is not more than 5% for t <∼ 1Myr.
We assume a uniform ISM with number density n0 and temperature T0 . A uniform
magnetic field is running in the z-direction as B 0 = (Bz , Br ) = (B0 , 0). As for the radiative
cooling rate, we use a spline function fit of cooling curves obtained by Raymond, Cox, &Smith (1976) (T > 104 K) and by Dalgarno & McCray (1972) (T < 104 K). We summarize
the parameters used in Table 1. The heating coefficient Γ0 is chosen as follows: in models
A-H (see Table 1), the heating rate Γ0 nH is balanced by the cooling rate Λ(T0 )nH 2 in theambient ISM. On the other hand, it is set to be 5×10−27 erg s−1 in models AW-HW (see Table1) as minimum heating rate given by cosmic-ray (Spitzer 1978). We begin the simulation
by adding thermal energy of E0 = 5 × 1050 erg within the sphere of 2pc in radius. We
assume a 10:1 H to He number ratio, which leads to a mean molecular weight for nuclei ofµ0 = 14/11 = 1.27 and that for fully-ionized ions and electrons of µt = 14/23 = 0.61, andthe number density of the Hydrogen atoms of nH = (10/11)n, where n means the number
density of the gas.

3.
3.1.

Results

Model with B0 = 5µG

In model A, we studied a magnetized ISM with B0 = 5µG, n0 = 0.2cm−3 , andT0 = 104 K, which mimics the warm ionized ISM. With a magnetic pressure of B02 /8π =10−12 dyn cm−2 (B0 /5µG)2 = 7.2 × 103 K cm−3 (B0 /5µG)2 and a thermal pressure of p =(µ0 /µt )n0 kT0 = 5.8×10−13 dyn cm−2 (n0 /0.2cm−3 ) (T0 /104 K) = 4.2×103 K cm−3 (n0 /0.2cm−3 )(T0 /104K), the plasma beta of the ISM is equal to β0 = 0.58(n0 /0.2cm−3 ) (T0 /104 K)(B0 /5µG)−2 . In Figure 1, we plotted cross-cut views along the z- (a and c) and the raxes (b and d). Upper panels (a and b) display density distributions, while the lower ones
(c and d) display temperature distributions. The shell formation time at which the SNRchanges its structure from the adiabatic Sedov phase (Sedov 1959) to the pressure-driven

–6–
snowplow phase is given (Cioffi et al. 1988) as3/14 E0n0 −4/7
4
tsf (E0 , n0 ) ≃ 7.8 × 10 yr.
5 × 1050 erg
0.2cm−3(9)

The first snapshot (t = 7.97 ×104yr; solid line) corresponds to the final stage of the adiabatic
phase or the transition stage to the pressure-driven snow-plow phase. Figure 1 shows thatthe differences in the z- and r-cross-cuts are slight at this stage. A dip in the temperature
distribution just after the shock front indicates that the radiative cooling begins to play apart at this stage.
The second snapshot at t = 1 × 105 yr shows clearly that a thin shell is formed [Fig.1(a)and (b)]. At this stage, Figure 1(a) and (b) indicates that the shell is thin and its widths
propagating in z- and r-directions are similar. The temperature distributions also look
similar to each other [Fig.1(c) and (d)]. Peak density in the shell reaches maximum at thisstage and declines gradually for the shell near the r-axis [Fig.1(b)], while it increases further
for the shell propagating near the z-axis [Fig.1(a)] until t = 2.52 × 105 yr.
At the age of 3.99 × 105 yr (dashed line), the SNR shows a shape far from sphericalsymmetry, that is, the shock front reaches R ≃ 62 pc while Z ≃ 57 pc. The shell width in
the r-direction ∆R becomes much wider than that in the z-direction ∆Z (short-dashed line)
as ∆R ≃ 14 pc while ∆Z ≃ 2 pc. The shape of the cavity which contains a hot low-density
5
gas indicates an elongated shape along the z-axis. That is, the bubble (T >∼ 10 K) half
radius is equal to Zbub ≃ 55pc while Rbub ≃ 47pc. This difference in z- and r-directions
comes from the effect of magnetic pressure which works to confine the gas in the r-direction.
The post-shock thermal pressure (not shown) reaches ≃ 4×104 K cm−3 on the z-axis, while itis equal to ≃ 2 × 104 K cm−3 for the shock propagating perpendicularly to the magnetic field(along the r-axis). In a case in which the magnetic field is parallel to the shock front, which
is applicable to the shock expanding in the r-direction, the post-shock to pre-shock pressure
ratio is given using the plasma beta of the pre-shock gas β0 and the Mach number of shockpropagation speed M [the Rankin-Hugoniot relation, see eq.(A8)]. Since part of the pressure
is attributed to the magnetic one, the post-shock thermal pressure of the magnetic shock
must be lower than that of the non-magnetic one (the ratio of post-shock thermal pressure of
the magnetic shock to that of the non-magnetic shock is plotted in Fig.11). With the plasma
β of the ISM β0 = 0.58 (model A) and the Mach number of the shock speed Vs ≃ 50km s−1of M ≡ Vs /cs = 3.3 at t = 3.99 × 105 yr (cs denotes the adiabatic thermal speed), equation(A8) anticipates that the thermal pressure of magnetized post-shock gas is 0.36 times smaller
than that of the non-magnetized one. This demonstrates that the post-shock magnetic fieldplays an essential role in determining the shock structure at this stage.
The shell width continues to increase with time in both z- and r-directions. Even in

–7–
the z-direction, the shell width increases. This is understood as follows: In the case ofhydrodynamic shock, the shell is confined under the effect of the post-shock pressure andthe thermal pressure of the hot interior gas. Ram pressure working at the front of the shell,
which is proportional to the shock speed squared, decreases with time, since the expansionof the shock front is decelerated. At the same time, the thermal pressure in the hot interior
gas also decreases owing to the adiabatic expansion. These thicken the shell thick (Chevalier1974). In the r-direction, the magnetic force prevents the shell from being compressed andworks to support the shell. The width of the shell increases owing to the decrease in theconfining pressures and the difference between the shells propagating in z- and r-directions
comes from the magnetic force .
In Figures 2 and 3, we plotted density distributions, total (thermal + magnetic) pressure
distributions and magnetic field lines. For each snapshot, the cross-cut views are shown in
Figure 1. Figures 2(a) and 3(a) show the structure of SNRs at the final stage of the adiabatic
phase (t = 7.97 × 104yr). The interior pressure (mainly thermal) predominates over the ISM
pressure and the SNR holds the spherical symmetry. The shell thickness is approximately
equal to 10% of the radius, which decreases to form a cooled shell [Figs. 2(b) and 3(b);t = 1 × 105 yr]. The SNR looks spherically symmetric in the adiabatic phase, even with themagnetic field. Compared with panels (c) of t = 3.99×105 yr, it is shown that the shell width
increases with time in the pressure-driven snow-plow phase after the cool shell formation,
especially in the r-direction.
As is shown in Figure 1, the distance of the shock front from the explosion site reaches
R ≃ 62pc in the r-direction, while Z ≃ 57pc in the z-direction (panels c of t = 3.99 × 105 yr).
This means that the shock wave propagates faster in the direction perpendicular to themagnetic field than in the parallel direction. Since the magnetic field lines are bent reaching
the shock front after passing the shock front, the shock running outward is a fast shock. The
phase velocity of the fast wave increases in the direction perpendicular to the magnetic field.
It is natural that the fast shock propagates faster in the r-direction than in the z-direction.
Thermal pressure in the hot low-density cavity is twice as high as the total pressure p+B 2 /8πof the ISM at this phase. However, it becomes lower than the total pressure of the ISM aftert ∼> 5 × 105 yr. Around this equilibrium epoch, the hot cavity begins to contract in ther-direction, while in the z-direction, the inner boundary of the shell continues to expand.
Figures 2(d) and 3(d) show the structure at 2.52 × 106 yr, at which the snapshot is also
made in Figure 1 (dash-dotted line). As indicated by a magnetic field line running nearest
to the z-axis in these figures, the inner boundary of the shell contracts radially between
3.99 × 105yr (R ≃ 47 pc) and 2.52 × 106 yr (R ≃ 28 pc), while the outward-facing shock front
continues to expand from R ≃ 62pc at t = 3.99 × 105 yr to R ≃ 145pc at t = 2.52 × 106 yr.

–8–
The portion traced by the above nearest magnetic field line corresponds to the temperature
4
5
jump between a gas in the shell T <∼ 10 K and a hot gas in the cavity T >∼ 10 K. As for the
expansion in the z-direction, the temperature transition located near Z ≃ 55 pc (3.99×105 yr)
expands to Z ≃ 95 pc (2.52 × 106 yr) shown in Figure 1(c), while the cavity contracts in ther-direction as the radial size of the cavity decreases from R ≃ 47 pc at t = 3.99 × 105 yr to
R ≃ 28 pc at t = 2.52 × 106 yr shown in Figure 1(d). This indicates that the expansion of
the hot gas is prevented by the effect of the magnetic field in the direction perpendicular to
the magnetic field. As for the shock front, panels (d) show that the angle measured fromthe shock normal vector to the pre-shock magnetic field line is positive for r >∼ 40 pc. In
other words, the magnetic field component parallel to the shock front does not change its
sign passing through the shock front. This is characteristic of the fast shock, while in theregion near the z-axis (r <∼ 40 pc), the parallel components have different signs comparing
the pre- and post-shock magnetic fields. Thus, this part of the shock is the intermediateshock (Jeffrey & Taniuti 1964).
The structure at t = 5.64 × 106 yr is shown in Figures 1 (long dashed line), 2(e) and
3(e). The expansion of the front continues. In the thick diffuse shell propagating in the rdirection, the thermal pressure has an excess of only ∼ 10% over the ISM and the excess totalpressure in the shell is at most ∼ 20% of the ISM, even including the magnetic pressure. The
shock is weak. The expansion in the r-direction is steadily decelerated till t = 5.64 × 106 yr.
6Although the deceleration along the z-axis continues for t <∼ 10 yr, the expansion speed
6
becomes constant V ∼ 20km s−1 after t >∼ 10 yr. The shell propagating along the z-axis is6
not decelerated for t >∼ 10 yr. This may seem strange. However, this can be explained asfollows: The mechanical equilibrium in the r-direction requires that the thermal pressure in
the hot low-density cavity pcav , where the magnetic field is weak, is in pressure equilibriumwith the total pressure of the ISM as pcav ≃ p0 + B02 /8π in the late phase of SNRs. This leads
to a thermal pressure imbalance of pcav > p0 . This difference of B02 /8π pushes the shell near
the z-axis, where no magnetic force works. The equation of motion for the part expanding
along the z-axis is written
dσV
B2
= 0,
(10)
dt
8π
R
where σ represents the column density of σ = ρdz. Constant force per unit area drives the
mass-accumulating shell to move with a constant velocity as σ ∝ Z ∝ t.

–9–
3.2.

Anisotropy Driven by the Magnetic Field

To see the departure from the spherically symmetric SNR more quantitatively, we define
the ellipticity of the SNR shell as
βsh ≡ Zsh /Rsh ,(11)
where Zsh and Rsh are the distances from the explosion site along the z- and r-axes and the
distance is measured for the point at which the density changes from that of the ISM in the
cross-cuts by more than 1% (Fig.1). As for the bubble, its ellipticity is calculated asβbub ≡ Rbub /Zbub ,(12)where Zbub and Rbub are the distances from the surface of the bubble to the explosion site.
The surface of the bubble is defined using the temperature isosurface of T = 105 K and theiso-density surface with the same density as the ISM n0 as
Rbub ≡ min(RT =105 K , Rn=n0 ),(13)

Zbub ≡ min(ZT =105 K , Zn=n0 ),(14)

where the former temperature isosurface indicates the surface in the pressure-driven phase5(t >∼ 10 yr) while the latter density isosurface does that in the adiabatic stage. We found theminimum of these two gives the bubble surface appropriately. It should be noted that these
two factors are chosen to be smaller than unity for ordinary magnetized SNRs (Zsh < Rsh
5
while Rbub < Zbub ). Figure 4 plots ellipticities βsh and βbub . After t >∼ 10 yr, βsh begins to
6decrease and reaches ≃ 0.8 after t >∼ 2 × 10 yr. This shows that the maximum difference
between Rsh and Zsh is 20% and the SNR is observed as elongated in the r-direction (anoblate spheroidal shape), if we take the shock front. Asymmetry in the bubble shape also
5appears t >∼ 10 yr and grows steadily. The bubble shrinks in the r-direction to form avery thin cigar-like structure βbub ≡ Rbub /Zbub ≪ 1. The hot cavity in the pressure-driven
expansion phase forms a prolate spheroid.
Figure 4 also plots the ratio of the bubble radius to the shell radius in respective directions as
αz = Zbub /Zsh,(15)andαr = Rbub /Rsh .

(16)

Since in the Sedov phase the shell size is approximately equal to one-twelfth of the shock
radius ∆R ≃ R/12 (Sakashita & Ikeuchi 1996), αz and αr keep constant ≃ 0.9 in the4adiabatic phase t <∼ 7 × 10 yr. As a cool shell forms, the difference between the shock radius

– 10 –and that of the hot cavity decreases. The quantities αz and αr reach ≃ 0.95 at the ageof t ≃ 105 yr. The shell width near the z-axis (∆Z = Zsh − Zbub ) continues to be thin for55t <∼ 4 × 10 yr and increases gradually after t >∼ 4 × 10 yr. On the other hand, the shell5
width in the r-direction (∆R = Rsh − Rbub ) increases steadily after t >∼ 10 yr (αr decreases
monotonically). Figure 4 shows us that asymmetry from the spherically symmetric expansion
begins just after the effective radiative cooling forms a thin shell. The diffuse shell formation
owing to the magnetic field is triggered just after this shell formation.

3.3.

Model with B0 = 1µG

In this section, we describe model C with a weak magnetic field of B0 = 1µG. Other
parameters are the same as model A. Since the magnetic pressure of the ISM is equal to
B02 /8π = 4 × 10−14 dyn cm−2 (B0 /1µG)2, the plasma β of the ISM is as large as β0 = 14.5.
Before and at the thin shell formation stage [Fig. 5(a) and (b)], the SNR is sphericallysymmetric and looks very similar to model A (first two panels of Figs. 2 and 3). However,
although the structure of the magnetized SNR in B0 = 5µG (model A) is far from being thespherically symmetric in the pressure-driven expansion stage [last three panels of Figs. 2
and 3], the SNR in a weak magnetic field (model C) remains its shape spherically symmetric[Fig. 5(c) and (d)]. For example, at t = 1.42 × 106 yr, the shell extends in the r-direction
from 65 to 89 pc and in the z-direction from 70 to 85 pc.
In Figure 4 (b), we plot the ratio of bubble-to-shock radii, and major-to-minor axis
ratio of the shock front and the bubble for model C. Compared with Figure 4 (a) of modelA, although the ratio of bubble-to-shock radii is as small as αr ≃ 0.5 at the age of 106 yr
for model A, it is αr ≃ 0.8 for model C. This means that the thickness of the shell, which
propagates in the direction perpendicular to the magnetic field (r-direction), is much reduced
in model C. The bubble begins to contract after 1.42 × 106 yr, which is much longer than thatof model A t ≃ 4 × 105yr. The difference comes from the effect of the magnetic force working
to compress the bubble. In both models A and C, αr decreases with time and reaches 0.1 int ≃ 4 × 106 yr (model A) and t ≃ 7 × 106 yr (model C).
As for the same ratio for the parallel (z-)direction, both models have similar values as
αz ≃ 0.93 for model A and αz ≃ 0.88 for model C at the age of 106 yr. However, although
6
model A remains at αz ≃ 0.85 for the final phase t >∼ 5 × 10 yr, αz continues to decreasein model C and αz ≃ 0.1 for t ∼ 7 × 106 yr. This means the shell propagating along thez-axis is confined to being thin in model A, while the shell thickness increases with timein model C. In the last paragraph of section 3.1, we already mentioned that the expansionof the shell propagating in the z-direction is driven by the thermal pressure in the bubble,

– 11 –
which is sustained by the radial compression due to the magnetic force. This seems to havethe effect of keeping the shell that is propagating in the z-direction thin.

3.4.

Effect of the Density

In this section, the models with high ISM density n0 = 1cm−3 are shown. Model E
has the same parameters as model A except for the interstellar density. The heating rate istaken to be 5 times larger than model A to keep T0 = 104 K.
Equation (9) gives the timescale to enter the pressure-driven snow-plow phase from
3/14
−4/7
the Sedov stage as tsf ≃ 3.1 × 104 yr (E0 /5 × 1050 erg)
(n0 /1cm−3 ), which is three−3times shorter than the models with low density n0 = 0.2cm . We plotted the snapshotof the adiabatic stage at t = 2.83 × 104 yr in Figure 6 (solid line). The next snapshot is at
t = 5.03 ×104 yr (dotted line), at which the density of the shell propagating in the r-direction
takes the maximum. This shows the thin shell phase of SNR.
At t = 5.03 × 105 yr (dashed line), Rbub takes the maximum Rbub = 29pc and the bubble
begins to contract in the r-direction. This timescale is almost the same as that of model A(t = 4 × 105 yr), even if the adiabatic to pressure-driven transition timescale tsf is three times
shorter than low-density model A. At this stage, the asymmetry between z- and r-directions
has developed. That is, the shell width of the r-direction ∆R = Rsh − Rbub ≃ 15pc and thatof the z-direction ∆Z = Zsh −Zbub ≃ 5pc differ. And as a consequence, the shell propagatingin the z-direction is thin, while that propagating in the r-direction is thick.
The shock wave propagates outwardly in both z- and r-directions. The next snapshot
is at t = 2.52 × 106 yr (dash-dotted line). The bubble has almost been crushed Rbub <∼ 10pc.
6In ≃ 3 × 10 yr the bubble is completely crushed. The contraction time-scale is much shorter
6
than that of the low-density models (>∼ 6 × 10 yr for model A).

3.5.

Effect of the Heating Rate

In this section, we compare models with the same parameters except for the heatingrate. In models AW-DW, we take Γ0 = 5 × 10−27 erg s−1 , which is ≃ 1/10 smaller thanprevious models A-D. In Figure 7, we plot snapshots at t = 2.52 × 106 yr of the temperature
distribution of models A (a) and AW (b) with B0 = 5µG and models C (c) and CW (d) with
B0 = 1µG. First, it is shown that the structures are very similar for models A and AW,
even if Γ0 differs ≃ 10 times. Although the height of the bubble Zbub differs several %, theradius of the bubble Rbub differs only slightly. On the other hand, the sizes of the bubbles

– 12 –
of models C and CW are different, as the bubble size differs ≃ 20% for both Rbub and Zbub .
To closely see model A, as for the structure along the r-axis, the shell (30pc < r <140pc) in model A has a temperature of T ≃ 104 K, while that for model AW is as low
as T ≃ 103 K owing to inefficient heating. As a result, the thermal pressure of model A[pth ≃ 4 × 103 cm−3 K] is much higher than that of model AW, and decreases from the postshock value of pth ≃ 4 × 103 cm−3 K to pth ≃ 3 × 102 cm−3 K near the bubble boundary.
Although the thermal pressure differs one order of magnitude, the expansion of the bubble
is unchanged. This is owing to the fact that magnetic pressure is dominant in the shell ofthese models.
As for the structure along the z-axis, the temperature of the shell of model AW (T <∼10 K) is much lower than that of model A (T ≃ 6 × 103 K). However, the thermal pressuresin the shell look similar, as follows: pth ≃ 3 × 103 cm−3 K (model A) and pth ≃ 2 × 103 cm−3 K(model AW). Further, those in the bubble are equal to pth ≃ 8 × 103 cm−3 K (model A) andpth ≃ 6 × 103 cm−3 K (model AW). Since the thermal pressure in the bubble is driving theshell in the z-direction, this similarity explains the fact that the expansion is very similareven if Γ0 changes 10 times. As a result, it is concluded that the evolution does not dependon the heating rate for models with strong magnetic fields of B0 >∼ 3µG.
3

In Figure 8, we plot the thermal pressure (a)-(b) and temperature distributions (c)-(d)of models C and CW. Before and at the shell formation stage (solid and dotted lines) thetwo evolutions are the same, since the heating rate has no effect on evolution at this stage.
After the shell formation (dashed, dash-dotted, and long dashed lines), these two models
show different evolutions. At t = 2.52 × 106 yr, the bubble occupies R <∼ 60pc in model C
while R <∼ 70pc in model CW, measured from the radius with a steep temperature gradient.
At t = 5.64 × 106 yr, these are R <∼ 35pc in model C and R <∼ 60pc in model CW (see also
Fig.9). This shows that the bubble shrinks to a greater degree in models with strong heating
rates.
Structures in the shell are significantly different: the thermal pressure and the temperature in model CW are much lower than those in model C. The low heating rate leads theshell to be cool and under low pressure. High pressure in the shell hinders the expansion of
the bubble in model C and thus the bubble size is reduced. The fact that the pressure in
the shell controls the bubble radius means the contraction of the bubble is partly owing to
the gas pressure in the shell. As a result, we conclude that for models with weak magnetic
field B0 <∼ 1µG, the shell thickness increases and bubble size decreases as the heating rateincreases.
Comparing models E and EW with n0 = 1cm−3 , even in the models with B0 =5µG, the– 13 –
size of the bubble depends on the heating rate as Rbub =26.3 pc (model E), 29.9 pc (model
EW) at t = 7.97 × 105 yr. The radius of the bubble differs ≃ 14% for Rbub between models
E and EW, and the total pressure along the r-axis in model E is twice as high as that inmodel EW at t = 7.97 × 105 yr. The difference seems to come from a large difference in Γ0(Γ0 =2.5 × 10−25 erg s−1 for model E and 5 × 10−27 erg s−1 for model EW). In high-density
models, the effect of Γ0 exceeds that of the magnetic field, and the evolution of the bubble
is different in models E and EW, in contrast to models A and AW.

3.6.

Expansion Law

In Figure 9, we plot the expansion law of the shell Zsh and Rsh for respective directionsof z- and r-axes as well as that of the bubble Zbub and Rbub for various models. Figure 9(a)shows that the expansion of the shock front Rsh is affected by the magnetic strength, if n0
and E0 are chosen to be the same. That is, the expansion is promoted with the strength ofthe magnetic field. This is related to the fact that the phase speed of the magneto-sonic wave
propagating in the direction perpendicular to the magnetic field increases with the magnetic
strength. However, the expansion of Rsh is independent of the value of the heating rate.
As for the shell expanding along the z-axis [Fig.9(b)], the models with B0 >∼ 3µG show
expansions different from those of from the models with weaker magnetic fields B0 <∼ 3µG.
The expansion of B0 >3µG
ispromotedbythestrengthofthemagnetic
field,
while
the∼expansion speed in Zsh for B <∼ 3µG finally approaches the adiabatic thermal speed cs . It has
been pointed out that the late-phase expansion of Zsh with constant speed comes from theexcess thermal pressure inside the bubble, which pushes the shell in the z-direction (§3.1).
Since the excess thermal pressure is well approximated by the ambient magnetic pressure as
B02 /8π ≃ 10−12 dyn cm−12 (B0 /5µG)2, this excess becomes greater than the ambient thermalpressure p0 ∼ 5.8 × 10−13dyn cm−12 (n0 /0.2cm−3 )(T0 /104 K) only for B0 >∼ 3.8µG. This seems
to explain the numerical result that the difference in Zsh disappears for B <∼ 3µG.
Figure 9(c) and (d) shows the bubble sizes measured along the r- and z-axes, respectively. Figure 9(c) shows the expansion of Rbub . This shows that the magnetic field suppresses
the expansion of the bubble radius Rbub , especially in the late phase evolution. This is anatural consequence of the fact that the contraction of the bubble is driven by the magnetic
pressure and tension forces. Besides this, there is a clear difference in the expansions ofmodels with strong heating (models A-D) and weak heating (models AW-DW). Especiallyin models B(W) - D(W) with B0 <∼ 3µG, models with weak heating have larger radii Rbub
compared with those with strong heating. Since Rsh is not affected with Γ0 , this means thatthe shell becomes thinner in models with lower heating rates, for weak magnetic field models

– 14 –
B0 <∼ 3µG. The fact that the shell density is higher in the weakly heated models BW-DW
than in the strongly heated B-D is explained by the fact that higher density is necessary toequilibrate heating and cooling in the shell, which will be discussed in §4.2.
Figure 9(d) shows the expansion of Zbub . There are two distinctly different evolutions of
the bubble expansion in the z-direction. In models with B0 >∼ 3µG, Zbub [models A(W) andB(W)] continues to expand. On the other hand, in models with B0 <∼ 1µG, Zbub expands
6
first but contracts after (1.5 − 3) × 10 yr [models C(W) and D(W)]. Although all models
C, CW, D, and DW show contraction in the final stage, the models with ineffective heating(CW and DW) have larger bubble sizes in the z-direction than those with effective heating(models C and D), respectively. This indicates models C and D have thicker shells thanmodels CW and DW, respectively. This is similar to the shells expanding in the r-direction.
In models with strong magnetic field, Zbub is larger for stronger magnetic field. This may
seem strange, since the magnetic field has no force working in the z-direction near the z-axis.
However, as already seen, the bubble is contracting in the r-direction owing to the magnetic
6
tension force for t >∼ 10 yr for models A(W) and B(W) [Fig.9(c)]. The thermal pressurein the hot interior, which mainly drives the shell near the z-axis, increases with increasing
magnetic field strength. This explains the fact that the inner boundary of the shell Zbub
expands faster for stronger magnetic fields in the models with a non-negligible magnetic fieldof B0 >
∼ 3µG.

4.
4.1.

Discussion
Porosity of ISM

To understand the structure of ISM, the volume fraction of the hot gas fh is a key issue.
The fraction of volumes covered by the hot cavities of SNRs is estimated by a quantity called
porosity, which is defined asZ∞Q ≡ rSNV (t)dt,(17)0where rSN is the supernova rate per volume as ∼ 10−13 yr−1 pc−3 (McKee & Ostriker 1977)
and V (t) represents the volume covered by a hot cavity. We assume all the SNRs end theirlives when the cavities contract and V (t) = 0. After obtaining the porosity, fh is estimatedas fh = 1 − exp(−Q). In Figure 10, we plot the porosity and the volume fraction of the
hot gas against the magnetic field strength. This shows clearly that the porosity decreases
with the magnetic field strength. Comparing models A-D (diamonds), Q = 0.6r−13 for nonmagnetic ISM while Q = 0.22r−13 for B0 = 5µG. Here r−13 represents the SN rate pervolume normalized by the standard value, 10−13 yr−1 pc−3 . The anticipated volume filling– 15 –
factor of hot gas fh = 0.45 for non-magnetic ISM (model D) is reduced to fh = 0.2 for
B0 = 5µG (model A). In models with ineffective heating rates (models AW-DW), the bubble
is compressed more weakly by the thick shell than models A-D. The porosity of these models
is greater than that of models A-D. In these models, the effect of the magnetic field in
reducing porosity is more remarkable, that is, Q = 1.95r−13 for B0 = 0 (model DW) while
Q = 0.28r−13 for B0 = 5µG (model AW). This anticipates fh = 0.86 (model DW) andfh = 0.24 (model AW).
Even in models with ISM density of n0 = 1cm−3 , the effect of the magnetic field in
reducing porosity is evident. Increasing B0 from 1µG to 5µG, porosity decreases by a factorof 1.7 from Q = 0.026r−13 to Q = 0.015r−13. In the case where the average density of the
ISM is high, the volume fraction of hot ISM is restricted unless the supernova rate is muchhigher than that of our Galaxy. Consequently, it is shown that the volume fraction of the
hot gas is restricted as fh ∼ 0.2 by the effect of the magnetic field, even if the average ISM
density is as low as n0 ≃ 0.2cm−3 .
Observationally, in the solar neighborhood, the “not strongly absorbing” HI gas contains
large holes, perhaps ranging up to 400 pc diameter, and the large holes occupy 10-20% of thevolume, fh ≃ 0.1 − 0.2 (Heiles 1980). For M31, another Sb galaxy, the filling factor which isderived for the HI holes, which are found in the range from 7 to 16 kpc from the nucleus, isabout 1%, fh ≃ 0.01 (Brinks & Bajaja 1986). For M33, an Sc galaxy, the filling factor for HI
holes, allowing a scaling factor of 2 for non-detected holes, is less than 40%, fh <∼ 0.4 (Deul& den Hartog 1990). HI holes arise by OB associations with multiple supernova explosions
that produce supercavities, and are also produced by the independent, random individual
supernova explosions (Heiles 1987). In this section, we considered the latter case, and theporosity and the filling factor of models for B0 = 5µG presumably looks consistent with theobservational data.

4.2.

Thermal Equilibrium

To understand how the shell width is determined, we consider the gas in thermal equilibrium. Balance between the radiative cooling Λ(T )n2 and heating Γ0 n leads
Λ(T )
kΓ0
=
,
T
p

(18)

4
where k represents the Boltzmann constant. The cooling function below T <∼ 10 K is anincreasing function of temperature and we assume Λ = Λ0 (T /104 K)q , ignoring the dependence of the ionization fraction x. Equation (18) gives the equilibrium pressure proportional– 16 –top ∝ Γ0 T 1−q∝ Γ0 1/q n(q−1)/q .

(19)
(20)

If we assume the post-shock gas is cooled and is in a thermal equilibrium [eq.(18)] under
the post-shock pressure, it is shown that the density decreases with the deceleration of theshock velocity asq/(q−1) 2q/(q−1)
pq/(q−1)n0Vsn ∝ 1/(q−1) ∝,(21)
1/(q−1)Γ0Γ0where we used the relation between the shock velocity Vs and the post-shock gas pressurep. Since the power q = 3, the shell density n decreases as the shell expands with time, since
the shock velocity Vs decreases with time.
In this discussion we have ignored the magnetic field, which works to prevent compression. In the thick shell sweeping the magnetic field, dynamics is controlled by the magnetic
field as shown in model A. In the extreme case, the change in density is controlled underthe magnetic field. Although equation (18) does not hold near the shock front, the temperature and the thermal pressure reach the equilibrium values expected from equation (18)near the shell-bubble boundary. We calculate the equilibrium density under the pressure p
in the shell, ρeq (r) along the r-axis. The difference from this equilibrium density is given as∆ = (ρ(r) − ρeq (r))/ρeq (r). The volume average of ∆ decreases with time as 0.91,0.67,0.14,
and 0.07 at t = 1 × 105 , 3.99 × 105 , 2.52 × 106 , and 5.64 × 106 yr respectively in model A. Thisshows that, with time, the density of the shell approaches the state of equilibrium expectedfrom equation (18).
Numerical calculations were carried out by VPP5000 in the Astronomical Data Analysis
Center, National Astronomical Observatory of Japan. This work was partially supported by
Grants-in-Aid for Scientific Research from MEXT, Japan (14540233,16036206). The authors
would like to thank the CANS (Coordinated Astronomical Numerical Software) project who
provided us with the MHD simulation code.

– 17 –A.

Rankin-Hugoniot Relation

In the case that the shock front is proceeding in the direction perpendicular to themagnetic field (x-direction), the mass flux conservation advected with fluid motion is writtenasρ1 vx1 = ρ2 vx2 ,(A1)
where subscripts 1 and 2 denote the quantities of pre-shock and post-shock, respectively, andvelocities are measured from the shock front which is moving with Vs . Other momentum,
energy, and magnetic flux conservations are written as
B12
B2
2
2
+ ρ1 vx1= p2 + 2 + ρ2 vx2
,
8π
8π




γ
vx1 B12γ
vx2 B22112
2
ρ1 vx1 +p1 += vx2
ρ2 vx2 +
p2 +
,
vx1
2
γ−1
4π
2
γ−1
4π
p1 +

vx1 B1 = vx2 B2 .

(A2)
(A3)
(A4)

Using
ρ2
vx1
B2
=
=
,
ρ1
vx2
B1
p2
y ≡
,
p1

x ≡

(A5)
(A6)

the above equations (A1)-(A4) give following relations (Priest 1982; Sakashita & Ikeuchi
1996)


2(2 − γ) 222
x + γ (γ − 1) M ++ 2 x − γ (γ + 1) M2 = 0,
(A7)
β1
β1


x2 − 1
1
2
−
+ 1,
(A8)
y = γM 1 −xβ1where M = Vs /(γp1 /ρ1 )1/2 and β1 = 8πp1 /B12 represent the Mach number of the shock
speed against the pre-shock sound speed and the plasma beta for the pre-shock medium,
respectively.
We plot the ratio of the normalized post-shock pressure y to that of non-magnetizedshock y(β1 → ∞) as Π ≡ y(β1)/y(β1 → ∞) against the shock Mach number M in Figure
11. Even if M is fixed at 10, Π increases with β1 as Π= 0.35, 0.67, 0.88, 0.96, and 0.99 forβ1 = 0.1, 0.3, 1.0, 3.0, and 10, respectively. This indicates that small M and β1 (slow shock
speed and strong magnetic field) induce the ratio Π to decrease, or in other words, the effectof the magnetic field becomes evident.

– 18 –
Using the Sedov solution, the expansion of the shock front is given as


1/5



1/5

E0
Rs (t) = 1.15
ρ0
for the radius andE0
Vs (t) = 0.46
ρ0t2/5 ,(A9)

t−3/5 ,

(A10)

for the expansion speed. These give the following typical values for Rs and Vs as
Rs
Vs

1/5 
2/5

E0
n0 −1/5t= 38.5pc,
5 × 1050 erg
0.2cm−3
105 yr
1/5 
−3/5
E0
n0 −1/5t−1= 151km s.
5 × 1050 erg
0.2cm−3
105 yr


(A11)
(A12)

Using the Mach number of the shock expansion speed M = Vs /cs and the plasma beta of the
pre-shock ISM, the post-shock gas pressure is obtained in relation to the pre-shock pressure.

– 19 –REFERENCES
Balsara, D., Benjamin, R. A., & Cox, D. P. 2001, ApJ, 563, 800
Brinks, E., & Bajaja, E. 1986, A&A, 169, 14
Brio, M., & Wu, C. C. 1988, Journal of Computational Physics, 75, 400
Chevalier, R. A. 1974, ApJ, 188, 501
Cioffi, D. F., McKee, C. F., & Bertschinger, E. 1988, ApJ, 334, 252
Cox, D. P. 1972, ApJ, 178, 159
Dalgarno, A., & McCray, R. A. 1972, ARA&A,10, 375
Deul, E. R., & den Hartog, R. H. 1990, A&A, 229, 362
Dickey, J. M., & Lockman, F. J. 1990, ARA&A, 28 215
Heiles, C. 1980, ApJ, 235, 833
Heiles, C. 1987, ApJ, 315, 555
Jeffrey, A., & Taniuti, T. 1964, Nonlinear Wave Propagation, (New York: Academic Pr.)
Koo, B. C., & Kang, J. H. 2004, MNRAS, 349, 983
McKee, C. F., & Ostriker, J. P. 1977, ApJ, 218, 148
Priest, E. R. 1982, Solar Magnetohydrodynamics, (Dordrecht: Reidel) Chap.5
Rand, R. J., & Kulkarni, S. R. 1989, ApJ, 343, 760
Raymond, J. C., Cox, D. P., & Smith, B. W. 1976, ApJ, 204, 290
Rubin, E., & Burstein, S. Z. 1967, Journal of Computational Physics, 2, 178
Sakashita, S., & Ikeuchi, S. 1996, Uchuryutairikigaku (Hydrodynamics in Astronomy) inJapanese (Tokyo: Baifukan) §3-7,§7-2
Sedov, L. I. 1959, “Similarity and Dimensional Methods in Mechanics,” (New York: Academic Press)
Slavin, J. D., & Cox, D. P. 1992, ApJ, 392, 131
Smith, R. K., & Cox, D. P. 2002, ApJ, 134, 283

– 20 –
Spitzer, L. J. 1978, Physical Processes in the Interstellar Medium, (New York: Wiley)
Tomisaka, K. 1998, MNRAS, 298, 797

This preprint was prepared with the AAS LATEX macros v5.2.

– 21 –

Table 1: Model parameters.
Model
A. . . . . .
B. . . . . .
C. . . . . .
D. . . . . .
AW. . . . . .
BW. . . . . .
CW. . . . . .
DW. . . . . .
E. . . . . .
F. . . . . .
G. . . . . .
H. . . . . .
EW. . . . . .
FW. . . . . .
GW. . . . . .
HW. . . . . .

n0(cm−3 )
0.2
0.2
0.2
0.2
0.2
0.2
0.2
0.2
1
1
1
1
1
1
1
1

T0
B0
(K) (µG)Γ0(erg s−1 )

104
104
104
104
104
104
104
104
104
104
104
104
104
104
104
104

5 × 10−26
5 × 10−26
5 × 10−26
5 × 10−26
5 × 10−27
5 × 10−27
5 × 10−27
5 × 10−27
2.5 × 10−25
2.5 × 10−25
2.5 × 10−25
2.5 × 10−25
5 × 10−27
5 × 10−27
5 × 10−27
5 × 10−27

5
3
1
0
5
3
1
0
5
3
1
0
5
3
1
0

– 22 –

(b)

101

101

100

100

10-1

10-1

n (cm-3)n (cm-3)

(a)

10-2
10-3
10-4
10-5
0

10-2
10-3
10-4

50

100
150
Z (pc)

200

10-5
0

250

50

109

109

108

108

107

107

106

105

104

104
50

100
150
Z (pc)

250

200

250

106

105

103
0

200

(d)

T (K)

T (K)

(c)

100
150
R (pc)

200

250

103
0

50

100
150
R (pc)

Fig. 1.— Density distributions along the z- (upper-left:a) and r-axes (upper-right:b) of SNRs
(model A). Those of temperature along the z- and r-axes are also shown in c (lower-left) andd (lower-right), respectively. Five snapshots are taken at the ages of t = 7.97 × 104 yr (solidline), 1 × 105 yr (dotted line), 3.99 × 105 yr (dashed line), 2.52 × 106 yr (dash-dotted line), and5.64 × 106 yr (long dashed line).

– 23 –
(a)

(b)

(c)

(d)

(e)

Fig. 2.— Structures of the magnetized SNR (model A). Density distributions (false color)and magnetic field lines (white lines) are shown. Five snapshots are taken at the ages of
7.97 × 104 yr (a), 1 × 105 yr (b), 3.99 × 105 yr (c), 2.52 × 106 yr (d), and 5.64 × 106 yr (e). Panels(a)-(c) cover the region of 150pc × 150pc, while panels (d) and (e) cover 250pc × 250pc.

– 24 –
(a)

(b)

(c)

(d)

(e)

Fig. 3.— Total pressure p + B 2 /8π distribution and the magnetic field lines for model A.
Five snapshots are taken at the ages of 7.97 × 104 yr (a), 1 × 105 yr (b), 3.99 × 105 yr (c),2.52 × 106 yr (d), and 5.64 × 106 yr (e). Panels (a)-(c) cover the region of 150pc × 150pc,
while panels (d) and (e) cover 250pc × 250pc.

– 25 –

(b)
1.2

1.0

1.0

0.8

0.8

0.6
0.4
0.2
0.0
104

α,β

α,β

(a)
1.2

βbub
βsh
αz
αr

0.6
0.4
0.2

105

106
time (yr)

107

0.0
104

βbub
βsh
αz
αr
105

106

107

time (yr)Fig. 4.— The ratios of the bubble radii to the shell radii as αz ≡ Zbub /Zsh and αr ≡ Rbub /Rsh .
We also plot the time variations of asymmetric factors βsh ≡ Zsh /Rsh and βbub ≡ Rbub /Zbub .
Shell radii Zsh and Rsh are calculated as the distance of the shell from the explosion site
along z- and r-axes, respectively. Bubble radii Zbub and Rbub are the sizes of the hot interior
cavity measured along the z- and r-axes. Models A (a: B0 = 5µG) and C (b: B0 = 1µG)
are plotted.

– 26 –

(a)

(b)

(c)

(d)

Fig. 5.— Structures of the SNR expanding in the weak magnetic field model C. Density
distributions (false color) and magnetic field lines (white lines) are shown. Four snapshots
are taken at the ages of 7.97 × 104 yr (a), 1.59 × 105 yr (b), 1.42 × 106 yr (c), and 2.52 × 106 yr(d).

– 27 –

(b)

101

101

100

100

10-1

10-1

n (cm-3)n (cm-3)

(a)

10-2
10-3
10-4
10-5
0

10-2
10-3
10-4

50

100

10-5
0

150

50

Z (pc)

150

100

150

(d)

108

108

107

107

106

106

T (K)

T (K)

(c)

105
104
103
0

100
R (pc)

105
104

50

100
Z (pc)

150

103
0

50
R (pc)

Fig. 6.— The same as Fig.1 but for SNR in n0 = 1cm−3 (model E). Structures of the SNR
expanding in the dense ISM with n0 = 1cm−3 . Five snapshots are taken at the ages of
2.83 × 104 yr (solid line), 5.03 × 104 yr (dotted line), 5.03 × 105 yr (dashed line), 2.52 × 106 yr(dash-dotted line), and 5.64 × 106 yr (long dashed line).

– 28 –

(a)

(b)

(c)

(d)

Fig. 7.— Comparison between models with strong heating rates and inefficient heating rates.
Temperature distributions are plotted for models A (panel a), AW (panel b), C (panel c),
and CW (panel d) at the age of 2.52 × 106 yr.

– 29 –

(b)

106

106

105

105
Ptot (cm-3K)

Ptot (cm-3K)

(a)

104
103

104
103

102
0

50

100
150
R (pc)

200

102
0

250

50

(c)
108

107

107

106

106

T (K)

T (K)

109

8

10

5

104

103
102
0

103
102
0

100
150
R (pc)

250

105

104

50

200

(d)

109
10

100
150
R (pc)

200

250

50

100
150
R (pc)

200

250

Fig. 8.— Comparison between models C and CW. Thermal pressure distribution along raxis of model C (a) and that of model CW (b) are plotted. Temperature distribution alongr-axis of model C (c) and that of model CW (d) are plotted. Snapshots are taken at the agesof 7.97 × 104 yr (solid line), 1.59 × 105yr (dotted line), 1.42 × 106yr (dashed line), 2.52 × 106 yr(dash-dotted line), and 5.64 × 106 yr (long dashed line).

– 30 –

(a)
250

250

Model A
Model AW
Model B
Model BW
Model C
Model CW
Model D
Model DW

200

150

Zsh (pc)

Rsh (pc)

200

(b)

100
50
0
0

Model A
Model AW
Model B
Model BW
Model C
Model CW
Model D
Model DW

150
100
50

2•106

4•106 6•106t (yr)

8•106

0
0

1•107

2•106

(c)
100

200
150
Zbub (pc)

Rbub (pc)

60
40

1•107

8•106

1•107

Model A
Model AW
Model B
Model BW
Model C
Model CW
Model D
Model DW

100
50

20
0
0

8•106

(d)
Model A
Model AW
Model B
Model BW
Model C
Model CW
Model D
Model DW80

4•106 6•106t (yr)

2•106

4•106 6•106t (yr)

8•106

1•107

0
0

2•106

4•106 6•106t (yr)

Fig. 9.— Expansion laws of the shell size, Rsh (a) and Zsh (b) as well as the size of the
bubble, Rbub (c) and Zbub (d) are plotted. Models A-D and AW-DW are shown. The figure
shows that the expansion of the shell is affected by the strength of the magnetic field B0 ,especially in Rsh . A stronger magnetic field induces faster expansion, while, the bubble size
(especially Rbub ) is affected by both the strength of the magnetic field and the heating rate.

– 31 –

(a)

(b)

2.0

1.0
Models AW-DW

Models AW-DW

0.8

Models A-D

1.5

Models EW-HW

0.6

Models E-H

1.0

Models E-H

fhQ

Models A-D

Models EW-HW

0.4
0.5
0.0
0

0.2

1

2

3
B (µG)

4

5

0.0
0

1

2

3
B (µG)

4

5

Fig. 10.— In panel (a), porosity of the volume occupied with the hot gas is plotted against thestrength of the magnetic field. We have assumed the supernova rate of rSN = 10−13 yr−1 pc−3 .
In panel (b), an anticipated volume-filling factor of the hot gas fh is shown. Diamonds(models A-D), squares (models E-H), triangles (models AW-DW), and crosses (models EWHW) represent models with different n0 and Γ0 .

– 32 –

1.0
0.8

Π

0.6
0.4
0.2
0.0
0

β1=0.10β1=0.30
β1=1.00
β1=3.00
β1=10.0

2

4

6

8

10

M
Fig. 11.— Rankin-Hugoniot relation for the MHD shock propagating perpendicular to themagnetic field. The post-shock gas pressure normalized by that attained in the non-magnetic
shock (Π ≡ y(β1 )/y(β1 → ∞)) is plotted against the shock Mach number (M). Each curve
corresponds to different plasma β.