M. A. de Avillez 1 - D. Breitschwerdt2
1 - Department of Mathematics, University of Évora,
R. Romão Ramalho 59, 7000 Évora, Portugal
2 -
Institut für Astronomie, Universität Wien,
Türkenschanzstr. 17, 1180 Wien, Austria
Received 8 October 2004 / Accepted 31 January 2005
Abstract
In star forming disk galaxies, matter circulation between
stars and the interstellar gas, and, in particular the energy input by
random and clustered supernova explosions, determine the dynamical
and chemical evolution of the ISM, and hence of the galaxy as a
whole. Using a 3D MHD code with adaptive mesh refinement developed
for this purpose, we have investigated the rôle of magnetized
matter circulation between the gaseous disk and the surrounding
galactic halo. Special emphasis has been put on the effect of the
magnetic field with respect to the volume and mass fractions of the
different ISM "phases'', the relative importance of ram, thermal
and magnetic pressures, and whether the field can prevent matter
transport from the disk into the halo. The simulations were
performed on a grid with an area of 1 kpc2, centered on the
solar circle, extending 10 kpc perpendicular to the galactic
disk with a resolution as high as 1.25 pc. The simulations were run
for a time scale of 400 Myr, sufficiently long to avoid memory
effects of the initial setup, and to allow for a global dynamical
equilibrium to be reached in case of a constant energy input rate.
The main results of our simulations are: (i) The
K gas is mainly concentrated in shock compressed layers,
exhibiting the presence of high density clouds with sizes of a few
parsecs and
K. These structures are formed in
regions where several large scale streams of convergent flow
(driven by SNe) occur. They have lifetimes of a few free-fall
times, are filamentary in structure, tend to be aligned with the
local field and are associated with the highest field strengths;
(ii) the magnetic field has a high variability and it is
largely uncorrelated with the density, suggesting that it is
driven by superalfvenic inertial motions; (iii) ram pressure controls the flow
for
K. For
K magnetic
pressure dominates, while the hot gas (
T> 105.5 K) in
contrast is controlled by the thermal pressure, since magnetic
field lines are swept towards the dense compressed walls; (iv) up
to
of the mass in the disk is concentrated in the classical
thermally unstable regime
K with
65% of the warm neutral medium (WNM) mass enclosed in the
K gas, consistent with recent
observations; (v) the volume filling
factors of the different temperature regimes depend sensitively on
the existence of the duty cycle between the disk and halo, acting
as a pressure release mechanism for the hot phase in the disk. We
find that in general gas transport into the halo in 3D is not
prevented by an initial disk parallel magnetic field, but only
delayed initially, for as long as it is needed to punch holes into
the thick magnetized gas disk. The mean volume filling factor of
the hot phase in the disk is similar in HD and MHD (the latter with
a total field strength of 4.4
G) runs, amounting to
17-21% for the Galactic supernova rate.
Key words: magnetohydrodynamics (MHD) - galaxies: ISM - ISM: evolution - ISM: bubbles - ISM: supernova remnants - ISM: structure
On the other hand, owing to the high sensitivity and large
throughput of the ROSAT XRT (Trümper 1983) and its PSPC instrument, a number of normal spiral galaxies with soft X-ray halos
were detected, e.g. NGC 891 (Bregman & Pildis 1994) and
NGC 4631 (Wang et al. 1995; Vogler & Pietsch 1996). In
some cases, even local correlations between H,
radio
continuum and soft X-rays were found (see Dettmar 1992), arguing for
local outflows as it had been suggested in the Galactic fountain
(Kahn 1981), chimney (Norman & Ikeuchi 1989) and the Galactic wind
model (Breitschwerdt et al. 1991; Breitschwerdt & Schmutzler 1994).
However, until a few years ago it was not possible to undertake an
extended 3D numerical study in order to follow the evolution of ISM
structures on small and large scales simultaneously, as was done
for the fountain by Avillez (1998). Recently, Korpi et al. (1999)
have numerically studied the evolution of superbubbles in a
magnetized disk and found that blow-out is much more likely than
previously thought, mainly because bubbles evolve in an
inhomogeneous background medium. Their 3D simulations are however
limited by the usage of a small grid (1 kpc perpendicular to the
disk), so that they were able to follow only the onset of a
supernova disturbed ISM for a short timescale of
100 Myr. As
a consequence, these results still bear the imprint of the initial
conditions, and the development of a disk-halo cycle could not be
followed.
It is a fortunate coincidence that exactly 100 years after the discovery of the ISM by stationary Ca II lines (Hartmann 1904), we are now for the first time in the position to give a detailed picture of the structure and evolution of the ISM on the scale of parsecs and below, that seems to converge (i.e., does not change substantially by merely increasing the numerical resolution). Our contribution here is to have carried out large scale high resolution 3D simulations of the ISM that include the Galactic magnetic field, background heating due to starlight and allow for the establishment of the duty-cycle between the disk and halo (commonly known as galactic fountain) by using a grid that extends up to 10 kpc on either side of the midplane. Our simulations capture both the largest structures (e.g., superbubbles) together with the smaller ones (e.g., filaments and eddies) down to 1.25 pc. We investigate, among other things, the variability of the magnetic field in the Galactic disk and its correlation with the density, the rôle of ram pressure in the dynamics of disk gas and the relative weight of the ram, thermal and magnetic pressures, the mass distribution and the volume filling factors of the different temperature regimes in the ISM. Other important issues like the variation of the volume filling factors of the ISM "phases'' with energy injection rate by SNe, the dynamics of the galactic fountain, the conditions for dynamical equilibrium and the importance of convergence of these results with increasing grid resolution have been treated in Avillez (2000) and Avillez & Breitschwerdt (2004, hereafter Paper I).
The outline of the present paper is as follows: in Sect. 2 the model and numerical setup used in the simulations is presented; Sect. 3 deals with the most important results obtained in these MHD simulations, their interpretation, and comparison with those in corresponding HD runs described in Paper I. In Sect. 4 a discussion of the results and comparison with other models is carried out. Section 5 closes the paper with a summary of the main results and some final remarks on future work.
We have run high resolution kpc-scale MHD simulations
of the ISM, driven by SNe (with a canonical explosion energy of
1051 erg) at the Galactic rate, on a Cartesian grid of
kpc size in the Galactic plane and extending
kpc into the halo. These runs use a modified version of the model
described in Avillez (2000) coupled to a three-dimensional numerical
AMR scheme that uses the MHD extension of the Piecewise Parabolic
Method (Dai & Woodward 1994, 1998) and the grid refinement procedure
discussed in Balsara (2001) assuring that B is divergence-free during the refinement process and simulation
time. The complete set of MHD equations that are solved numerically
can be written in the conservative form
The model includes supernovae types Ia (randomly distributed in the
space and having a scale height of 325 pc) and II (randomly
distributed or clustered in OB associations formed in regions with
density and temperature thresholds of cm-3 and
K, respectively). The number and masses of the OB stars are
determined by an appropriate IMF following Massey et al. (1995); each
star explodes after its main sequence lifetime has elapsed, using a
time estimate by Stothers (1972). The trajectories of all O and B
stars are followed kinematically by attributing to each a random
velocity at time of formation so that the location of each exploding star is known at any time. For further details we refer to
Paper I.
The initial conditions comprise the setup of (i) the fixed vertical
gravitational field
provided by the stars in the disk; of
(ii) the interstellar gas having a density stratification
distribution that includes the cold, cool, warm, ionized and hot
"phases'' in the Galaxy; and (iii) a magnetic field with uniform
and random components given by
and
,
respectively; here
G is the uniform
component field strength at the Galactic midplane, n(z) is the
number density of the gas as a function of distance from the
Galactic midplane, where it has a value of
.
Periodic boundary conditions are applied along the boundaries
perpendicular to the Galactic plane (i.e., the z-direction), while
outflow boundary conditions are imposed at the top (z=10 kpc) and
bottom (z=-10 kpc) surfaces (parallel to the galactic plane). The
resolution in the grid varies on the fly according to density and/or
pressure gradients thresholds with 1.25 pc being the highest
resolution of three levels of refinements, achieved in the region
kpc, while on the rest of the domain the
finest resolution is 2.5 pc (two levels of refinement). The
resolution of the coarse grid is 10 pc.
The initial evolution of the disk and halo gases is characterized by
the build-up of the random component of the magnetic field in addition
to the establishment of the duty disk-halo-disk cycle (see
Paper I). The random field component is generated during the first 15
Myr as a result of turbulent motions mainly induced by shear flows
due to SN explosions. For the remainder of the evolution time the
random field component and the mean field oscillate around their mean
values of 3.2 and G (Fig. 1), respectively,
and the average total field strength amounts to
G. The build-up of the disk pressure by
supernovae and the establishment of the duty cycle, is similar to that
seen in the HD runs discussed in Paper I (see also Avillez 2000), that
is, the initially stratified distribution does not hold for long as a
result of the lack of equilibrium between gravity and (thermal,
kinetic and turbulent) pressure during the "switch-on phase'' of SN
activity. As a consequence the gas in the upper and lower parts of the
grid collapses onto the midplane, leaving low density material behind.
However, in the MHD run it takes a longer time for the gas to descend
towards the disk. The complete collapse onto the midplane, like in the
hydrodynamical (HD) case, is prevented as a result of the opposing
magnetic pressure and tension forces. As soon as enough SNe have gone
off in the disk building up the required pressure support, transport
into the halo is not prevented, although the escape of the gas takes a
few tens of Myr to occur - somewhat longer than in the pure HD
case. The crucial point is that a huge thermal overpressure due to
combined SN explosions can sweep the magnetic field into dense
filaments and punch holes into the extended warm and ionized H I
layers as seen in Fig. 2, which shows the density and
magnetic field distribution in a plane perpendicular to the Galactic
midplane at time t=221.9 Myr (the warm neutral and ionized H I
gases are the light blue and greenish regions seen in the density
distribution). Once such pressure release valves have been set up,
there is no way from keeping the hot overpressured plasma to follow
the pressure gradient into the halo. As a consequence, transport into
the halo is not prevented and the duty disk-halo-disk cycle of the hot
gas is fully established, in which the competition of energy input and
losses into the ISM by SNe, diffuse heating, radiative cooling and
magnetic pressure leads the system to evolve into a dynamical
equilibrium state within a few hundred Myr (for an analytical estimate
see Paper I). This time scale is considerably longer than that quoted
in other papers (e.g., Korpi et al. 1999; Kim et al. 2001), because
these authors do not take the galactic fountain into account. The
vertical structure of the thick disk and halo is addressed in more
detail in a forthcoming paper.
![]() |
Figure 1: Time evolution of the random field component (red), the mean field (green) and total field strength (black). |
Open with DEXTER |
![]() |
Figure 2:
Two-dimensional slices through the three-dimensional data cube showing
the vertical (perpendicular to the midplane) distribution of the
density ( left panel) and magnetic field ( right panel) at time
t=221.9 Myr. The mean field in the midplane is ![]() |
Open with DEXTER |
![]() |
Figure 3: Two-dimensional slices through the three-dimensional data set at z=0 (Galactic midplane), showing density ( top left panel), temperature ( top right panel), pressure Phys. Rev. E/k ( bottom left panel), and magnetic field ( bottom right panel) after 358 Myr of evolution. Colour bars indicate the logarithmic scale of each quantity. Red refers to lowest density, highest temperature and pressure, while dark blue refers to the highest density and lowest temperature and pressure. In the magnetic field image red/blue correspond to the highest/lowest strength. Direct visual inspection shows a strong inverse spatial correlation between gas density and temperature (although the pressure diagram is still far from uniform), whereas the correlation between magnetic field strength and gas density is only present for the cold dense regions. For higher temperature regions this correlation becomes rather poor and vanishes completely for the hot gas (s. Fig. 4). |
Open with DEXTER |
Figure 3 shows slices of the three-dimensional data cube of the density, temperature, pressure Phys. Rev. E/k, and magnetic field distributions in the Galactic midplane at 358 Myr of evolution with a resolution of 1.25 pc. We can directly compare the morphology of the ISM between this run and the similar HD run discussed in Paper I. In both cases the highest density (and lowest temperature, if the gas had enough time to cool) gas tends to be confined into sheet-like structures (filaments in 2D) which constitute shock compressed layers (SCLs) as a result of the compression by shocks mainly driven by SNe in their neighbourhood. The SCL structures are thicker in the MHD run than in the HD run, because the compression of a shock wave at given Mach number is less as energy has now to be shared between gas and magnetic field. The T< 103 K gas in the MHD run tends to be aligned with the magnetic field, while in the HD run there is no preferred orientation. The highest densities up to 800 cm-3, observed in both runs occur in regions where several streams of convergent flow meet. The orientation of these streams is random. We note in passing, that if clouds are hit by SCLs from random directions, turbulence in the interior is generated. Tapping the turbulent ISM energy reservoir, which is huge, represents a neat way of sustaining supersonic turbulence inside molecular clouds, which as it is known, dissipates very efficiently. We will explore this idea in more detail in a forthcoming paper.
In the MHD run cold denser regions are dominated by the magnetic field and are more filamentary in structure than low density regions (which have low field strengths, although, as seen in Fig. 4, a wide range of B-field strengths is spanned) due to the anisotropy introduced by the strong field, which is much better correlated with the gas density for lower than for higher temperature gas.
During their initial expansion phase, young supernova remnants are
roughly spherical as long as the radius is below or at the order of
the magnetic field curvature scale, but as time progresses magnetic
tension forces increase as the radius of curvature of the field
lines decreases as
,
where
and
are the radius of field curvature and tension force,
respectively, and
is the uniform field component. It is
straightforward to estimate the average radius of a superbubble, or
respectively the time, when the shell stalls, because in the long
run magnetic tension forces beat the thermal pressure force in the
hot interior, since tension increases and pressure decreases with
time. This argument assumes some kind of ideal configuration, e.g. that the shock normal is perpendicular to the magnetic field
direction, which as can be seen in the simulations, does in general
not exist, as the field is tangled and inhomogeneous. It should also
be kept in mind that such an analytic analysis ignores the
possibility of MHD instabilities that would promote mixing with the
ambient medium and escape of hot bubble gas.
![]() |
Figure 4:
Scatter plot of B versus ![]() ![]() ![]() ![]() ![]() |
Open with DEXTER |
In the following we assume that the outer shell will have
cooled off and collapsed into a thin shell before slowing down
considerably, being subject to magnetic and thermal pressure. Let
be the shell radius,
be the polar angle between the
radial and the magnetic field direction
,
and
the
shell thickness, so that the bubble radius is given by
.
For
the outer shock is purely gasdynamic, as
compression along the B-field is unimpeded. However flux conservation
causes the field to wrap around the bubble and carry a tangential
component as
increases. The bending of the lines of force
near the polar region leads to a deflection of the incoming flow and
thus to a
-dependent momentum transport. As this reduces the
inertia of the shell near the polar region, the deceleration there is
more effective leading to a "dimple'' (s. Ferrière et al. 1991). Whereas the pressure is purely thermal here, the magnetic
pressure increases with
and will dominate completely for
.
Flux conservation at co-latitude
in the thin shell approximation requires
![]() |
= | ![]() |
|
= | ![]() |
(5) |
Numerical simulations show (e.g. Ferrière et al. 1991) that in
case of interstellar magnetic field strengths of the order 3-5
the evolution until the outer shock decays into a fast
magnetosonic wave can still be well described by a similarity
solution (assuming a homogeneous ambient medium) given by McCray &
Kafatos (1987). In the case of quasi-sphericity, the momentum
equation for self-similar flow, and using Eq. (6), reads
To solve Eq. (8) we need an expression for the relative
shell thickness, .
In the equatorial region, the pressure
in the shell at later times will be dominated by magnetic pressure
so that momentum conservation across the shock requires
As long as the ambient thermal pressure is small compared to the
pressure inside the bubble and the magnetic pressure in the shell,
and combining Eqs. (10)-(12), and inserting them into Eq. (8), we
obtain for the time
of maximum bubble expansion (for
,
where tension forces are strongest)
The radius of maximum expansion (perpendicular to the B-field
direction) is therefore given by
Perpendicular to the galactic plane the expansion of bubbles will be alleviated due to the density stratification of the extended disk and Rayleigh-Taylor instabilities can occur once the shell accelerates over a few density scale heights. The effect of the magnetic field then will be similar to surface tension in an ordinary fluid, limiting the growth of the most unstable modes to a maximum wave number. It is therefore expected that shell fragments or cloudlets should have a minimum size. A number of shells will stall before that happens, and these will appear in simulations and observations, like e.g. NGC 4631 (Wang et al. 2001) as extended loops.
During the evolution of the system, thermal and dynamical processes
broaden the distribution of the field strength with density so that
after the global dynamical equilibrium has been set up, the field
strength in the disk spans two orders of magnitude from 0.1 to 15 G (see Fig. 4). The spreading in the field strength
increases with temperature, being largest for the hot (
T>
105.5 K) and smallest (between 2 and 6
G) for the cold
(
K) gas. The lowest (0.1
G) and higher (larger
than 10
G) field strengths are associated with
T> 104.2gas and densities of 0.1-1 cm-3. The intermediate temperature
regimes have the bulk of their scatter points similarly distributed
between 2 and 6
G, although there is a clear distinction between
the different temperature regimes with respect to the density
coverage. It can be seen that the
K gas
covers a large fraction of the scatter points overlapping with the
coolest and partially with the hot regimes, and there, mainly with
those of low densities and high field strengths. It almost completely
covers points located in the
K both in
the field strength and in the density distribution. The thermally
unstable gas (
K) completely overlaps the
K scatter points, although the latter covers smaller
field strength intervals than the unstable gas. The
gas coverage with other (stable) regimes is an
indication that it is in a transient state trying to move towards
thermally stable regions at
K or to the
T>105.5 K range. Whether this gas becomes part of a
classically stable region, however, depends if it has enough time to
cool down (or to be heated up) with respect to dynamical time scales,
which change its location on the p-V-diagram. The higher field
strengths at the intermediate temperature regimes are related to the
increase in compression of the field lines in newly formed shells and
SCLs.
The large scatter in the field strength for the same density,
seen in Fig. 4, suggests that the field is being
driven by the inertial motions, rather than it being the agent
determining the motions. In the latter case the field would not be
strongly distorted, and it would direct the motions predominantly
along the field lines. In ideal MHD field diffusion is negligible, and
the coupling between matter and field should be perfect (we are
of course aware of numerical diffusion, which weakens the argument at
sufficiently small scales). Therefore gas
compression is correlated with field compression, except for strictly
parallel flow. The high field variability is also seen in the right
bottom panel of Fig. 3, which shows a highly turbulent field,
that seems to be uncorrelated with the density, and thus, the
classical scaling law
,
with
,
according
to the Chandrasekhar-Fermi (CF) model (1953) will not hold. It should
be kept in mind that in CF it was assumed that the field is distorted
by turbulent motions that were subalfvénic, whereas in our
simulations in addition both supersonic and superalfvénic motions
can occur, leading to strong MHD shocks. In other words, according to
the CF model, the turbulent velocity is directly proportional to the
Alfvén speed, which in a SN driven ISM need not be the case.
If ram pressure fluctuations in the ISM
are dominant as our simulations suggest, we would indeed not expect a
perfect correlation with
but a broad distribution of Bversus
.
Although in general
would be
expected, it should be noted that in reality heating and cooling
processes, and even magnetic reconnection could induce further
changes, making the correlation rather complex.
![]() |
Figure 5:
Comparison of the averages
![]() ![]() ![]() ![]() ![]() |
Open with DEXTER |
In order to gain insight into the driving forces in the simulated disk
(
pc) we turn now to the relative weight of
the ram, thermal and magnetic pressures as function of the temperature
averaged over bins of
K
(Fig. 5). The average distributions shown in the
figure have been calculated by using 51 snapshots with a time interval
of 1 Myr during the time evolution of 350 and 400 Myr.
The relative importance of thermal pressure over ram and magnetic
pressures increases with temperature. The gas with
K is mainly thermally driven, while at lower temperatures thermal
pressure has no influence on the dynamics of the flow. For
K,
indicating that the gas is
dominated by the Lorentz
forces, and
the magnetic field determines the motion of the fluid. At the
intermediate temperatures ram pressure dominates, and therefore, the
magnetic pressure does not act as a significant restoring force (see
Passot & Vázquez-Semadeni 2003) as it was already suggested by the
lack of correlation between the field strength and the density. It is
also noteworthy that in this temperature range (
K)
the weighted magnetic pressure is roughly constant, suggesting that
the magnetic and thermal pressures are largely independent, whereas
both thermal and ram pressures undergo large variations in this
temperature interval.
We therefore conclude that in galactic disk environments the magnetic field will have a strong influence on the cold gas, particularly by controlling cloud and star formation processes. On the other hand, the warm and hot phases of the ISM, which are dominating by volume but not by mass, are affected dynamically only to a minor extent by magnetic pressure forces. However, it should be kept in mind, that microscopic processes, like heat conduction, will be strongly modified by the presence of a magnetic field.
The averaged (over the period 350-400 Myr) volume weighted histograms
of the thermal, magnetic and total (
,
which is calculated at each grid cell, and then used in the
construction of the corresponding histogram) pressures of the
different temperature regimes in the magnetized ISM
(see Fig. 6) show that for volume fractions
the full distributions (shown by the dashed lines) have a pressure
coverage running from three (for
)
to more than five (for
)
orders of magnitude and having smooth (for
)
and steep (for
and
)
peaks located in the
dyne cm-2 interval. These
peaks are determined by the contributions of the
K and
K regimes followed by
the hot gas. The thermally stable regimes
K and
K, by this order, have the smallest contribution to the
peaks.
In the thermal and magnetic pressure histograms it is possible to
distinguish two regions (
and
dyne cm-2), which are determined by the gas with
and
T> 104.2 K (with the hot gas being the
main contributor for the high pressure region), respectively. In the
histogram (that has a conical structure covering the
dyne cm-2 region for
)
all regimes contribute to these two regions, although the
hot gas is the main contributor to the low pressure region and the
cold gas contributes mainly to the high pressure range.
A comparison between the top and middle panels of the figure indicates
that the highest values of the total pressure, i.e.,
dyne cm-2, are mainly of thermal origin. Contributions
come from the gas with
T>104.2 K and in particular from the
hot gas (whose histogram has a central peak located around
dyne cm-2 and two smaller peaks in the right wing of
the histogram at
dyne cm-2, resulting from
SN events). The gas with
has a very small
contribution and the cold gas practically none. For the remaining
pressure range (
dyne cm-2) both magnetic
and thermal pressures contribute, with the magnetic contributions
having a larger weight than the thermal ones.
![]() |
Figure 6: Averaged volume weighted histograms of the thermal ( top panel), magnetic ( middle panel) and total ( bottom panel) pressures for the different temperature regimes in the ISM. The total distribution of each pressure is shown by the dashed black line. The histograms have been calculated using 51 snapshots with a time interval of 1 Myrbetween 350 and 400 Myr. |
Open with DEXTER |
In order to gain insight into the distribution of the disk gas as a
function of density or temperature, it is useful to construct volume
weighted histograms. They give us information about the occupation
fraction of a certain "gas phase'' with respect to volume averaged
density and temperature. Therefore a peak in the histogram tells us,
at which density/temperature a particular phase is concentrated,
because here the occupation fraction is highest.
A comparison between the averaged (over the time 350 through 400
Myr) volume weighted histograms of the density and temperature of
the magnetized and unmagnetized disk gas (Fig. 7) shows
noticeable differences that include the (i) decrease/increase by an
almost order of magnitude in the histograms density/temperature
coverage; (ii) change in the relative weight of the dominant
temperature regimes (in the density histograms) and consequently;
(iii) changes in the pronounced bimodality of the total density and
temperature histograms.
![]() |
Figure 7: Averaged volume weighted histograms of the density for the different temperature regimes of the disk gas (the full distribution is shown by the dashed line) in the top panels and averaged temperature histograms in the bottom panels for the MHD run ( left column) and the HD run ( right column). The histograms have been calculated using 51 snapshots with a time interval of 1 Myr between 350 and 400 Myr. |
Open with DEXTER |
The density histograms in both HD and MHD runs are dominated by the
T> 105.5 K (less important),
K and
K regimes. The latter regime is the
most important in the HD run contributing to the pronounced peak
around
cm-3 in the total density distribution,
whereas in the magnetized ISM the latter two regimes interchange in
importance (there is a strong reduction in the peak of the
K gas). This leads to a shift of the peak in
the total density histogram to
cm-3, where it
is mainly dominated by the
K gas
contribution. The hot gas averaged PDFs have a similar shape but
different density coverage, being larger in the HD run by an order
of magnitude and peak location (being a factor of two higher in the
MHD run). This suggests that the areas underneath the histograms, that
is, the volume filling factors of the hot gas, are not much
different in the two runs. The smallest peaks seen in the density
histograms belong to the thermally stable regimes
K
(which has the smallest occupation fraction), followed by the
K regime. With the inclusion of the
magnetic field, there is an increase in the density coverage by the
cold gas histogram and thus of its occupation fraction in the disk.
The variations seen in the density weighted histograms are also seen
in temperature PDFs having different structures. In effect, while the
temperature PDF of the HD run has a bimodal structure (as it has two
peaks: one at 2000 K and another around 106 K), in the MHD run
the decrease/increase in importance of the
K/
K leads to the
reduction/increase of the occupation fraction of these regimes and
therefore to a change of the histogram structure appearing it to be
unimodal. This variation of the intermediate region appears to be an
effect of the presence of the magnetic field, with the smoothing
effect being less pronounced for lower field strengths in the disk.
![]() |
Figure 8: Time history of the volume occupation fractions of the different temperature regimes for the MHD ( top panel) and HD ( bottom panel) runs. |
Open with DEXTER |
![]() |
Figure 9:
Top and middle panels show the histories of the mass
fractions of the different thermal regimes for the magnetized ( top)
and unmagnetized ( middle) ISM. Bottom panel shows the history of the
fraction of mass of the WNM gas having
![]() |
Open with DEXTER |
Table 1: Summary of the average values of volume filling factors, mass fractions and root mean square velocities of the disk gas at the different thermal regimes for the HD and MHD runs.
![]() |
Figure 10: History of the root mean square velocity of the differenttemperature regimes of the disk gas in the MHD ( top panel) and HD( bottom panel) runs. |
Open with DEXTER |
During most of the history (t> 100 Myr) of ISM simulation the
occupation ()
and mass
(
)
fractions of the different thermal
regimes (Figs. 8 and 9 - top and middle panels)
have an almost constant distribution, varying around their mean
values (cf. Table 1). The thermally stable regimes with
K and
K have similar occupation
fractions of
5% and
10%, respectively, in both runs,
while the hot gas has an increase from
17% in the HD run to
20% in the MHD case. By far most of the disk volume is occupied by
gas in the thermally unstable regimes at
and
K with similar occupation
fractions
30% in the MHD run, while in the HD run these
regimes occupy
and
,
respectively, of the disk volume.
The major fraction of the disk mass (93% and
80% in the HD and MHD runs, respectively) is occupied by the gas
with
and transferred between the cold and the
thermally unstable regimes with the corresponding mass changes. The
cold (
K) regime during most of the simulation time
encloses more mass than in the MHD run until
320 Myr when
both runs show a similar amount of cold mass (
40%), while
the thermally unstable regime encloses the remaining mass. This
results from the fact that in the MHD run the gas compression is
smaller than in the HD case, because part of the shock energy goes
into work against magnetic pressure. The remaining ISM mass is
distributed between the other temperature regimes with the
K and
K regimes
enclosing a total of
7% and
16% of mass in the HD
and MHD runs, respectively, and the hot gas enclosing <1% of the
disk mass in both runs.
For most of the simulation time, in both runs, 60-70% of the warm
neutral mass (
K) is contained in the
K temperature range (bottom panel of
Fig. 9). This result in excellent agreement with the
observations of Fitzpatrick & Spitzer (1997) who found a value of
60% and Heiles & Troland (2003) who determined a lower
limit of
50% (see also Heiles 2001).
In both the MHD and HD runs the root mean square (rms) velocity,
,
which is a measure of the disordered motion of the gas,
increases with temperature as can be seen in Fig. 10,
which shows the
for the different thermal regimes in the MHD
(top panel) and HD (bottom panel) runs. The average rms velocity
(
)
in the last 100 Myr of evolution of both
media for the different thermal regimes is shown in Table 1. The
figure and table show that there are large fluctuations in the
of the different thermal regimes in the HD run, which are
reduced due to the presence of the magnetic field.
The rms velocities are in good agreement with those
estimated from observations of cold/cool neutral
(7-10 km s-1; e.g., Kulkarni & Fich 1985; Kulkarni & Heiles 1987), warm neutral (14 km s-1; Kulkarni & Fich 1985), warm ionized (
30 km s-1;
Reynolds 1985) and hot (
33-78 km s-1; e.g., Zsargó et al. 2003) gases.
Again, the near constancy of the rms velocity with time indicates the
presence of a dynamical equilibrium, with random motions, i.e. thermal and turbulent pressures adding to the total pressures,
provided that the energy injection rate remains constant on a global
scale.
In the previous sections we compared high resolution three-dimensional kpc-scale HD and MHD simulations of the ISM, driven by SNe at the Galactic rate, fully tracking the time-dependent evolution of the large scale Galactic fountain (up to 10 kpc on either side of the Galactic midplane) for a time sufficiently long (400 Myr) so that the memory of the initial conditions is completely erased, and a global dynamical equilibrium is established. These simulations show how important the establishment of the duty cycle is, and how it affects the global properties of the ISM, e.g., occupation fraction of the hot gas, just to name one. An essential feature of the hot gas is its escape into the halo establishing the fountain flow, thereby reducing the expansion volume in the underlying disk significantly. It is not possible to suppress such a flow even in the presence of an obstructing disk parallel magnetic field. All that is needed for break-out of the gas into the halo is a sufficient overpressure in the superbubbles with respect to the ambient medium. The effect of the magnetic field is just adding another pressure component. Its topology is almost irrelevant unless the galactic magnetic field is extremely high and/or the SN rate is very low. The field may inhibit the break-out of an individual remnant, but certainly not the high-pressure flow resulting from supernova explosions in concert within an OB association. Thus, the occupation fractions of the hot gas in the disk in the HD and MHD simulations are not too much different, although there is a slight increase in the MHD run as a result of the magnetic tension forces aiding to confine bubbles and Lorentz forces obstructing mixing with cooler gas. In Paper I we have shown that the volume filling factor of the hot regime strongly correlates with the supernova (and therefore star formation) rate. Although we have not carried out MHD simulations for higher supernova rates, it seems very plausible that this correlation will persist in the MHD case, since higher rates will produce more hot plasma which, as we have shown here, is not magnetically controlled.
We also emphasize the importance of carrying out three-dimensional simulations in order to describe the evolution of a magnetized ISM adequately. The idea that a disk parallel magnetic field could suppress break-out and outflow into the halo was mainly based on 2D simulations carried out in the last 15 years, owing to computing power limitations. It is obvious, that in 2D-MHD, the flow perpendicular to the magnetic field lines (and hence to the galactic plane) is subject to opposing and ever increasing magnetic tension and pressure forces. In 3D however, field lines can be pushed aside and holes and channels can be punched into the decreasing gas and field with z-height, allowing pressurized flow to circumvent increasing tension and pressure forces in z-direction. It is exactly this behaviour that we see in our simulations. This is less surprising if we think of similar problems of terrestrial plasma confinement, although here plasma instabilities prevail.
The lack of correlation between the field strength and density found
in the present simulations is supported observationally from
measurements of the magnetic field strength in the cold neutral medium
(Troland & Heiles 2001). A similar result is shown in the simulation
described in Kim et al. (2001), although these authors claim that the
magnetic field scales as
at densities n>1 cm-3.
Their Fig. 2 shows an almost order of magnitude variation of the
field with the density for the coolest gas. However, the distribution
of the scatter points is somewhat different than that seen in
Fig. 4 of the present paper. This discrepancy can be
explained by the (200 pc)3 box that Kim et al. have used,
centered in the Galactic midplane having no density stratification and
with periodic boundary conditions in all the box faces, driven by SNe
at a rate of 12 times the Galactic value and using an uniform field
strength of 5.8
G orientated along the x-direction; thus they
were unable to describe the disk-halo-disk circulation and did not
allow for a global dynamical equilibrium to be established (see also
Mac Low et al. 2004).
The distribution of the disk mass in the warm neutral medium (WNM),
in particular in the WNM fraction within the classical thermally
unstable regime seen in the simulations is strongly supported by
interferometric (Kalberla et al. 1985). Moreover optical/UV
absorption-line measurements (Spitzer & Fitzpatrick 1995;
Fitzpatrick & Spitzer 1997) indicate that a large fraction (63%) of the warm neutral medium (WNM) is in the unstable range
500<T<5000 K, whereas 21 cm line observations (Heiles 2001;
Heiles & Troland 2003) provide a lower limit of
for the WNM gas in this unstable regime. Direct numerical simulations of the
nonlinear development of the thermal instability under ISM
conditions with radiative cooling and background heating discussed
in Gazol et al. (2001) and Kritsuk & Norman (2002) show that about 60% of the system mass is in the thermally unstable regime.
However, it is unclear from their simulations what is the time
evolution of this mass fraction and what is explicitly the origin of
the unstable gas, although these authors suggest that ensuing
turbulence is capable of replenishing gas in the thermally unstable
regime by constantly stirring up the ISM. We have carried out
detailed numerical studies of the stability of the ISM gas phases
(Avillez & Breitschwerdt 2005), and could verify the hypothesis
that SN driven turbulence is capable of replenishing fast cooling
gas in classically unstable regimes.
Nonetheless, it is interesting to ask why such a large amount of gas can
exist in a thermally unstable regime. Obviously the simple Field
(1965) criterion, according to which instability is expected to set
in, viz.
,
where
is the heat loss function per unit mass,
and T and P are the temperature and thermal pressure of the fluid
element, respectively, is not adequate. It is the turbulence that
can have a stabilizing effect thereby inhibiting local condensation
modes. The situation is reminiscent to the existence of the solar
chromosphere, which consists of gas at around 105 K, again in the
thermally unstable regime. The reason is because heat conduction can
prevent thermal runaway on small scales. In other words, diffusion
processes have a stabilizing effect. In our case it is turbulent
diffusion that replaces the rôle of conduction, again, most
efficient for large wavenumbers. The turbulent viscosity
can be orders of magnitude above the
molecular viscosity, with Re being the Reynolds number of the
flow. What happens physically then, is that with increasing eddy
wavenumber
,
the eddy crossing time
(with
being the turbulent velocity
fluctuation amplitude) becomes shorter than the cooling time
,
where
is
the interstellar cooling function. Although not really applicable
here, it is instructive to see that in case of incompressible
turbulence following a Kolmogoroff scaling law, where the energy
dissipation rate is given by
,
we obtain a lower cut-off in wavelength,
where thermal instability becomes inhibited, if
A concern with the present simulations is the numerical resolution,
which is always a limitation, to which HD/MHD simulations are
subject. The crucial questions are (i) how do the results change as
the resolution is increased; and (ii) how is the amount of cold gas
affected by a resolution larger than the Field length,
pc in the WNM and 0.001 pc in the CNM (see e.g., Audit &
Hennebelle 2005), far below the limit one can handle in global
simulations. Note that both Koyama & Initsuka (2002) and Audit &
Hennebelle (without magnetic fields) used 2D boxes of 0.3 pc and 20 pc
in length, three and two orders of magnitude smaller than our box size
and cell resolutions of the order of 0.01 parsec. The first of these
questions has been discussed at length for the HD run in Paper I,
where it has been shown that the large scale structure and global
properties of the ISM of the different thermal regimes remain
essentially unchanged, once a minimum resolution of 1.25 pc is
achieved. With regard to (ii) the calculation of the cooling length
depends also on heat conduction, which may be suppressed in stochastic
magnetic fields as the ones seen in the present simulations (for a
detailed discussion see Paper I; see also Malyshkin & Kulsrud 2001). However,
as shown above in a dynamical ISM the Field criterion can be
superseded, as the cooling fluid element may be overturned by
turbulent fluctuations in the velocity and density field and
therefore, some of the gas that is heated up or cooling down to an
unstable regime may not participate in a transition to a stable phase
at lower temperature.
The present simulations still neglect an important component of the
ISM, i.e. high energy particles, which are known to be in rough
energy equipartition locally with the magnetic field, the
thermal and the turbulent gas in the ISM. The presence of CRs and
magnetic fields in galactic halos is well known and documented by
many observations of synchrotron radiation generated by the relativistic electron
component. There is increasing evidence that the fraction of these
cosmic rays that dominates their total energy is of Galactic origin
and can be generated in supernova remnants via the so-called
diffusive shock acceleration mechanism to energies up to 1015 eV (for original papers see Krymski et al. 1977; Axford et al. 1977;
Bell 1978; Blandford & Ostriker 1978, for a review Drury 1983, for
more recent calculations see Berezhko 1996). It is also common
wisdom that the propagation of these particles generates MHD waves
due to the so-called streaming instability (e.g. Kulsrud & Pearce
1969) and thereby enhances the turbulence in the ISM. In addition,
self-excited MHD waves will lead to a dynamical coupling between the
cosmic rays and the outflowing fountain gas, which will enable part
of it to leave the galaxy as a galactic wind (Breitschwerdt et al. 1991, 1993; Dorfi & Breitschwerdt 2005). Furthermore, as the cosmic
rays act as a weightless fluid, not subject to radiative cooling,
they can bulge out magnetic field lines through buoyancy forces.
Such an inflation of the field will inevitably lead to a Parker type
instability, and once it becomes nonlinear, it will break up the
field into a substantial component parallel to the flow (Kamaya et
al. 1996), thus facilitating gas outflow into the halo. It has been
suggested that eventually reconnection with an undisturbed halo
field component (Tanuma et al. 2003) and/or due to twisted
approaching flux tubes driven by Coriolis forces (Hanasz et al.
2002), will dissipate some of the field energy and change the field
topology. This may even provide the necessary ingredients for a
Galactic
-dynamo to operate, generating a large scale
poloidal field as was originally proposed by Parker (1992). It is
still a matter of debate, if Galactic halo conditions are suitable
for anomalous resistivity to occur in current sheets with a minute
thickness of the ion Larmor radius in order to promote fast
reconnection at about 10% of the Alfvén speed. In any case will
such a scenario alleviate the transport of hot plasma into the halo.
The effect of cosmic rays on large scale ISM evolution will also be
the subject of a forthcoming paper.
In the present paper we investigate the rôle of magnetic field and matter circulation between the disk and halo and its effect on the dynamical evolution of the disk gas. These results are compared to HD simulations of the same resolution that we have carried out as well. The main results of this paper can be summarized as follows:
Acknowledgements
M.A. and D.B. are partially supported by the ESO/FCT (Portuguese Science foundation) grant PESO/P/PRO/40149/2000. D.B. thanks G. Hensler for useful discussions on the subject. M.A. benifited from discussions with E. Vazquez-Semadeni and M.-M. Mac Low.