Astrophys Space Sci (2021) 366:94
https://doi.org/10.1007/s10509-021-03993-9
REVIEW ARTICLE
Non-equilibrium ionisation plasmas in the interstellar medium
Dieter Breitschwerdt1·Miguel A. de Avillez1,2
Received: 14 September 2020 / Accepted: 26 August 2021 / Published online: 1 October 2021
© The Author(s) 2021
Abstract Obtaining astrophysical information from diffuse
cool, warm and hot plasmas in interstellar and intergalactic
media by electromagnetic radiation is based on highly non-
linear heating and cooling processes, which are largely de-
termined by atomic physical time scales and reaction rates.
To calculate spectra is further complicated by gas dynami-
cal interactions and processes, such as e.g. shock waves, fast
adiabatic expansion and catastrophic cooling. In essence this
leads to a non-linear coupling between atomic physics and
hydro- or magnetohydrodynamics, which renders radiative
cooling to become time- and space-dependent, contrary to
the often conveniently used assumption of collisional ion-
isation equilibrium for optically thin plasmas. Computing
power and new algorithms for high performance comput-
ing have made it possible to trace the dynamical and ther-
mal evolution of a sufficiently large section of interstellar
space over an appreciable time scale to derive character-
istic quantities like temperature and density distribution as
well as spectra, which can be compared to X-ray, UV and
optical observations. In this review we describe diffuse in-
terstellar plasma simulations, the physical processes which
drive the temporal and spatial evolution, and present high
resolution numerical simulations, including time-dependent
This article belongs to the Topical Collection: Plasma, Particles, and
Photons: ISM Physics Revisited. Guest Editors: Manami Sasaki,
Ralf-Jürgen Dettmar, Julia Tjus.
D. Breitschwerdt
M.A. de Avillez
mavillez@galaxy.lca.uevora.pt;
1Department of Astronomy and Astrophysics, Berlin Institute of
Technology, Hardenbergstr. 36, 10623 Berlin, Germany
2Department of Mathematics, University of Évora, R. Romão
Ramalho 59, 7000 Évora, Portugal
cooling, which further our understanding of the state and
evolution of interstellar (magnetised) plasmas. We also dis-
cuss briefly the rôle of cosmic rays and their interaction with
the plasma.
Keywords Interstellar medium: hydrodynamics ·
Magnetohydrodynamics ·Supernovae ·Turbulence
1 Introduction
The interstellar medium (ISM) in galaxies is a complex
multi-component system composed of a highly compress-
ible plasma, high energy particles (cosmic rays), magnetic
fields, photons and dust, just to name the most important
ingredients. This complexity makes a complete theoretical
description virtually impossible, let alone the fact that the
interaction between individual components is mostly non-
linear. Although numerical simulations on parallel comput-
ers have boosted the subject and our understanding during
the last two decades, we are, despite all the efforts, still far
away from a complete picture. Arguably, it might be tempt-
ing to throw all the components and their respective interac-
tions on a supercomputer and predict the temporal and spa-
tial evolution from the formation of stars until the explosion
of many generations of supernovae a billion years later for a
whole galactic disk and its halo. But the validity and the pre-
dictive power of such an endeavour may be disappointingly
small, because one has to sacrifice spatial resolution and re-
sort to “recipes” rather than treating all processes and inter-
actions of components by their appropriate physical equa-
tions. To illustrate this, consider a parcel of fluid which has
been overrun by a shock wave - an event that happens in the
ISM approximately about every 10 million years. It will be
94 Page 2 of 21 D. Breitschwerdt, M.A. de Avillez
compressed and heated up (and if it is magnetised, magneto-
hydrodynamic (MHD) waves will be driven into the ambient
medium), presumably be ionised and driven out of thermal
and ionisation equilibrium. All this happens on microscopic
scales, and introduces various length and time scales, such
as electron and ion gyroradius, collisional length and time
scales, ionisation and recombination time scales etc., which
numerical codes have to tackle – and often enough fail to
do so. Moreover, in interstellar plasmas, there exist quite a
number of hydrodynamic and plasma instabilities, which de-
rive from a large amount of free energy available in the sys-
tem that can be tapped by a plethora of physical processes.
On top of this, there is also a coupling between micro- and
macroscales by the ubiquitous turbulence, which is often su-
personic.1Therefore we believe that progress at this stage of
research is mainly achieved by analysing the main compo-
nents and their governing physical processes – this again is
debatable – on a reasonably sized computational grid. This
computational box should be large enough, especially per-
pendicular to the disk, to prevent the material from leaving
the computational domain during the simulation time, be-
fore some sort of dynamical equilibrium might have been
established. Conversely, it should be small enough, such that
spatial resolution does not become a severe issue in pro-
viding reliable results. For example density inhomogeneities
and clumpiness of the ISM should be resolved, because the
average of the number density, n, squared does not equal the
square of the average density, i.e. n2=n2. This is par-
ticularly important in the case of radiative cooling, which
depends on n2locally. Too low spatial resolution therefore
underestimates energy losses by escaping photons in an op-
tically thin plasma.
Before devising an elaborate model, it is important to sort
out what we believe are the most important physical pro-
cesses that should be included. First and foremost, all major
energy sources and sinks should be scrutinized. Since the
70’s supernovae (SNe) have been identified as the major en-
ergy input (see Cox and Smith 1974; McKee and Ostriker
1977, henceforth CS74 and MO77, respectively). Soon it
was realized that core collapse SNe mostly occur in their
parent clusters and are thus correlated in space and time,
and that many HI-shells were much larger (see Heiles 1979,
1984) than expected from a single explosion. The theoret-
ical description of superbubbles (e.g. McCray and Kafatos
1987; Mac Low and McCray 1988) gave a plausible expla-
nation. For some time it was discussed that these superbub-
bles (SBs) could not break out of the disk in the presence
of a disk-parallel magnetic field, because magnetic tension
forces would stall the SB expansion, as some numerical sim-
ulations suggested. However, these were two-dimensional,
1This leads to effective dissipation by (weak) shock waves. In any case
turbulence decays fast enough that it needs to be driven, thus requiring
a continuous and substantial energy input rate.
so that the shock heated bubble had to work against ever
increasing magnetic tension forces, giving rather an upper
limit to the effects of the magnetic field than a plausible sce-
nario for the SB expansion. As an aside, it is well-known
how difficult plasma confinement is in terrestrial fusion de-
vices such as tokamaks or stellarators. One reason is current
driven plasma instabilities. In the ISM there are also many
instabilities due to the amount of free energy available, as
has been mentioned before. In addition turbulent reconnec-
tion (Lazarian and Vishniac 1999; Lazarian et al. 2020) can
also lead to plasma escape. In essence, we believe that the
disk is rather porous for SN heated plasma flows.
Later, the observations of vertical dust lanes (Sofue 1987)
have led to models in which the energy transport perpen-
dicular to the disk into the halo plays a vital rôle (Norman
and Ikeuchi 1989), venting a lot of energy out of the disk,
and thereby decreasing volume filling factors in accordance
with observations. Some of the material blown out into the
halo might rain down again in form of a Galactic fountain
(Shapiro and Field 1976; Bregman 1980; Kahn 1981;de
Avillez 2000). Here the rising gas cools down as the work
done against the gravitational pull must be delivered from its
internal energy. In addition radiative cooling exacerbates the
energy crisis. Therefore the escape speed vesc =√−20,
where 0equals the gravitational potential at the base of
its journey, is only a lower limit for the necessary initial
speed in case of an unbound flow (wind). In fact additional
energy has to be dumped into the halo, e.g. by SNe Type
Ia, which have a more extended vertical distribution than
core collapse SNe. While the plane-parallel calculations of
Kahn (1981) point towards a round-trip taking 200 Myr in
the Milky Way near the solar circle, numerical simulations
(de Avillez 2000) show that the cycle is established over a
period of 150-200 Myrs.
Driving a galactic wind needs more than the average
energy input from SN heated gas at the base of the halo,
and will therefore occur only locally. However, it has been
shown that a low mass wind can be driven by the assistance
of cosmic rays (Ipavich 1975; Breitschwerdt et al. 1991;
Everett et al. 2008; Dorfi and Breitschwerdt 2012; Uhlig
et al. 2012; Wiener et al. 2017; Farber et al. 2018)bya
non-linear interaction between cosmic rays (CRs), their self-
excited waves and the thermal plasma. It is noteworthy that
the presence of CRs is inevitable, and hence their dynamical
effects cannot be ignored (as it has been done in most mod-
els in the past), because particles are accelerated in SN gen-
erated shock waves predominantly by the first order Fermi
mechanism (see Sect. 5). Streaming along field lines during
their escape from the Galaxy, and being strongly scattered
by field irregularities, e.g. magnetohydrodynamical (MHD)
waves, leads to a so-called streaming instability, because the
particle distribution function in the waves’ frame becomes
increasingly anisotropic thus leading to a resonant gener-
ation of waves that tend to remove this anisotropy. Waves
Non-equilibrium ionisation plasmas in the interstellar medium Page 3 of 21 94
with wavelengths of the order of the particle gyroradius are
excited resonantly, which in turn leads to a strong scattering
of the particles themselves. These scatterings carry momen-
tum from the CRs to the plasma via the waves as a media-
tor, and hence assist in driving a wind as the particles leave
the Galaxy. Since these processes are a very important topic
in Galaxy evolution and formation, as cosmological simula-
tions show, and because CRs form an important component
of the ISM, we will discuss them in more detail in Sect. 5.
In the 90’s, with increasing computing power, global sim-
ulations of the ISM were carried out, taking SNe and stel-
lar winds as its main driving sources. In a seminal paper
Rosen and Bregman (1995) were one of the first who set
up a two-dimensional model, in which star formation was
parametrized, and the interaction between stars and gas was
taken into account in a simplified way. The grid had a 2 kpc
extension in the disk and ±15 kpc perpendicular to it. They
could reproduce a multi-phase filamentary medium as ob-
served. Since turbulence is an integral part of such simu-
lations, three-dimensional models are necessary in order to
follow the realistic evolution of vorticity. Korpi et al. (1999)
presented such a simulation of a SN driven ISM including
a magnetic field on a grid with a horizontal extension of
0.5 kpc and a vertical height zof ±1 kpc. They find again
a multiphase gas with turbulent cell sizes up to 60 pc. Since
a box size of 1 kpc in either z-direction leads to the loss
of rising gas, of which most of it in reality would come
back in a fountain, it is necessary to extend the height to at
least ±10 kpc. Using a correspondingly larger box size, de
Avillez and Breitschwerdt (2004,2005) have carried out hy-
drodynamical and MHD simulations with 1 ×1 kpc in the
disk (x−yplane) and ±10 kpc in the vertical z-direction
(in later simulations ±15 kpc, cf. de Avillez and Breitschw-
erdt (2012a)) with a spatial resolution of 0.5 and 1 pc, re-
spectively. Although this is still above the critical wave-
length of thermal instability (so-called Field length, Field
(1965)), the probability density functions (pdfs) show that
the amount of gas at small scales does not dominate, and
hence the global error in the cooling function does not be-
come too large. Simulations at decreasing scales (increas-
ing resolution) have shown that quantities like volume fill-
ing factors seem to converge at scales of about 1 pc and
below. Apart from being sensitive to the clumpiness of the
ISM, the cooling function (T ,Z), where Tis the temper-
ature and Zthe metallicity, cannot always be assumed for a
plasma to be in collisional ionisation equilibrium (CIE) as it
is mostly done for convenience. In CIE, the rate of ionisa-
tions is exactly balanced by the number of recombinations at
any time, and therefore the cooling rate L=ρ2(T ), where
ρis the mass density, depends only on the temperature for
a given metallicity of the medium.2This allows observers
2In reality, this is an approximation, in which one assumes that the
plasma is optically thin, and that the major source of ionisation is due
to deduce the temperature of an interstellar plasma from the
emission measure in UV or X-rays, when the density is low
enough for it to be optically thin. Hence there is no feed-
back on the ionisation structure by reabsorption. In addition,
for densities below 104cm−3, collisional de-excitation does
not play a significant rôle, so that the line emission is de-
termined by collisional excitation (coronal approximation).
For example the detection of the ubiquitous OVI absorption
line is commonly identified with a plasma at ∼3×105K.
However, there is a caveat. The assumption of CIE is only
approximately valid at temperatures above 106K (Shapiro
and Moore 1976). For lower temperatures, and this concerns
most of the ISM in mass and volume, ionisation and recom-
bination rates are different, simply due to different temper-
ature dependent cross-sections. This has been realized long
ago by Shapiro and Moore (1976), who studied the cooling
of a 106K shock heated gas, showing the progressive devia-
tion from CIE. Shapiro and Moore (1977), also analyzed the
X-ray emission in solar flare coronal non-equilibrium ion-
isation (NEI) plasmas, impulsively heated to 108K. In the
ISM, adiabatic compression and expansion of the gas play
an important rôle, because shock and rarefaction waves are
ubiquitous due to SN explosions. This leads to fast colli-
sional heating or cooling of the electrons, which then subse-
quently share energy and momentum with the ions. There-
fore the kinetic energy of the electrons changes on a short
timescale, whereas ionisation and recombination lag behind.
There exist two extreme possibilities in the plasma state, viz.
under- and overionisation. If the plasma gets shocked, the
kinetic temperature of the electrons rises quickly, whereas
the ionisation state still resembles the upstream state. Con-
versely, fast adiabatic expansion of a hot plasma bears the
ionisation signature of the hot state, although the gas can be
cold (see Breitschwerdt and Schmutzler 1994,1999). As we
will show, this can result in the emission of soft X-rays at
temperatures of only 104K (see Sect. 4). The fundamen-
tal difference to CIE is now, that the ionisation state varies
in space and time, and so does the cooling function in non-
equilibrium ionisation (NEI) plasmas, i.e. radiative cooling
is time-dependent (de Avillez and Breitschwerdt 2012a).
Other groups have taken different routes and emphases
in global ISM simulations and have also used different tech-
niques. Because it is easy to use and quite versatile, even
for complicated geometries, it is Galilean-invariant and it
is has very good conservation properties, smoothed parti-
cle hydrodynamics (SPH) has become very fashionable. Al-
ready developed in the late 70’s by Gingold and Monaghan
(1977), Lucy (1977) and others, it has been used in many
fields of astrophysics, also for the ISM. It is a Lagrangian
to collisions without any significant external photon field being present.
Hence every free electron is likely to recombine again with an ion, and
the photons it generates can escape.
94 Page 4 of 21 D. Breitschwerdt, M.A. de Avillez
method, and its main feature is a discretization in mass in-
stead of space, it is meshfree and can therefore be applied
to free flows, and it can be easily parallelized. On the down-
side is its resolution dependence on mass, i.e. in the ISM,
the spatial resolution of the hot tenuous medium is usually
poor; also its ability of treating dynamical instabilities is
problematic. For a comparison to grid based methods, see
Agertz et al. (2007). Because of these shortcomings, hybrid
schemes have been recently developed, which try to com-
bine the advantages of both Eulerian grid and Lagrangian
particle methods, like e.g. the AREPO code (Springel 2010),
which was developed originally for Galilean invariant cos-
mological simulations. In short, a set of particles, which can
in principle move freely, is distributed as an unstructured
mesh by Voronoi tessellation. The hyperbolic equations are
solved by a Godunov scheme with Riemann solver, thus re-
taining the shock representation of grid schemes. As limiting
cases, stationary particles mimic a Eulerian grid, and freely
moving particle the SPH method.
Dobbs et al. (2011), using an SPH method, have simu-
lated a whole galactic disk in order to analyse the distri-
bution of the gas in general, and the formation of molecu-
lar clouds in a starforming disk with feedback, in particular.
They are able to reproduce the cloud spectrum and find that
more massive clouds require a spiral potential. In a series of
papers Walch et al. (2015), Girichidis et al. (2016), Peters
et al. (2017), Girichidis et al. (2018)takeanevenmoream-
bitious road (the so-called SILCC Project), studying the ISM
from molecular clouds scale including chemistry, effects of
radiation fields and SN feedback on the ISM to the launch
of outflows, on box sizes ranging from 500 pc in the disk to
±500 pc perpendicular to it. On time scales large enough to
complete the galactic fountain cyle, however, the box size
remains a critical issue (see above).
The structure of this review is as follows: Sect. 2de-
scribes the numerical setup and modelling of ISM plas-
mas on a mesoscale, as well as the results of HD and
MHD simulations with respect to the structure of the ISM.
In Sect. 3, a brief historical overview and the necessity
for non-equilibrium ISM simulations is discussed, while
in Sect. 4, one of the fundamental implications, i.e. time-
dependent cooling, is outlined in detail, and observational
consequences are presented. Section 5deals with the im-
portance of cosmic rays in the ISM and their interaction
with magnetic fields and the gas, with galactic winds be-
ing just one important example. In Sect. 6the Local Bub-
ble, in which our solar system is embedded, is modelled and
treated as a test laboratory for non-equilibrium ionisation
plasmas, showing some more recent results on the distribu-
tion of highly ionised species and electrons there. Section 7
concludes this review and summarizes the results.
2 Numerical modelling of the ISM
In order to follow the evolution of the ISM plasma over
a reasonable amount of time (about 400 - 500 Myr in the
simulations discussed here), it is necessary to carry out the
computations on parallel processors, e.g. using Open MPI.
In order to obtain sufficient spatial resolution, we use the
technique of adaptive mesh refinement (AMR), which has
been pioneered by Berger and Oliger (1984) and Berger
and Colella (1989). In a nutshell, AMR simulations start
with a coarse grid spread over the whole computational vol-
ume, which, according to a user defined refinement crite-
rion, e.g. density or pressure gradient thresholds, is over-
laid by a locally finer grid and evolved in time. During the
procedure, fluxes of conserved quantities have to be strictly
balanced across different refinement levels. This dynamical
grid adaptation procedure is repeated as often as necessary
to meet the criterion, and once it proves negative, the finer
grid can be removed. In essence, many astrophysical prob-
lems, which typically cover a wide range of spatial scales,
are only tractable with AMR. Specifically, in the ISM, we
have quite often interfaces (e.g. propagating shocks, contact
discontinuities), which demand a locally high spatial reso-
lution, but cover only a comparatively small volume of the
grid.
2.1 Hydrodynamic simulations
The physical and numerical setup of high resolution 3D
mesoscopic hydrodynamic simulations with AMR is as fol-
lows (for details see de Avillez and Breitschwerdt (2004),
with refinements later on, like local self-gravity and heat
conduction, see e.g. de Avillez and Breitschwerdt (2012a):
(i) background heating due to starlight depending on zis in-
cluded according to Wolfire et al. (1995), being constant par-
allel to the Galactic disk, initially balancing radiative cool-
ing at 9000 K; (ii) the cooling function has been tabulated
and taken from Dalgarno and McCray (1972) in earlier ver-
sions, with the ionisation fraction being 0.1 at temperatures
below 104K with a cut-off in temperature at 10 K; (iii) SNe
have been modelled according to Freeman (1987) with SNe
type Ia having a scale height of 325 pc; (iv) we use den-
sity and temperature thresholds of 10cm−3and 100 K, re-
spectively, for OB star formation, with their number, masses
and main sequence life times, τMS, being determined from
the initial mass function (IMF). About 60% of them explode
within clusters, while the remaining go off in the field with
a scale height of 90 pc, having been given a random ve-
locity upon star formation. The OB associations have been
assumed to form in dense regions with a scale height of 46
pc. Stars explode within a time scale given by τMS, so that
the explosion intervals are determined by their respective SN
rate. In more recent publications, we have also included self-
gravity of the plasma and thermal conduction, although we
Non-equilibrium ionisation plasmas in the interstellar medium Page 5 of 21 94
still believe that turbulent mixing is the more efficient trans-
port process. The boundary conditions are outflow at the top
and bottom of the box, and periodic at the sides (perpendic-
ular to the disk), thus avoiding a depletion of matter within
the long simulation times. This of course introduces an arte-
fact, e.g. in the disk flows due to the periodic boundary con-
ditions, but, since we can expect that the gas in the galactic
potential does not change significantly over scales of a few
kiloparsecs in the disk, the results should still be reasonably
close to reality. The initial setup includes the distribution of
cold, cool, warm, ionised and hot gas (see Ferrière 1998)in
a gravitational potential provided by the stars in the Milky
Way.
2.1.1 General results
An important finding of HD ISM simulations is the estab-
lishment of the Galactic fountain cycle, which operates on
a time scale of about τgf ∼150 Myr, implying that sim-
ulations of smaller time spans cannot be a sensible repre-
sentation of a realistic ISM. This value can be estimated
from the distance and typical velocity it takes to reach the
Galactic fountain return point, which from simulations is at
zc∼15 kpc, roughly corresponding to the critical point in a
stationary flow, which would hence be subsonic for z<z
c.
At a typical lower halo temperature of T∼106K, as seen
in X-rays, the speed of sound is cs∼130 km/s, and thus
τgf ∼zc/cs≈113 Myr, which is also the sound crossing
time. One could argue that since the amount of mass being
circulated through the halo is only of the order of 15%, the
dynamics of the disk gas is only slightly affected. However,
this would ignore the large amount of thermal and kinetic
energy that is vented into the halo, which therefore has a
substantial influence, e.g. on the compression of disk gas
and volume filling factors. In addition, one has to bear in
mind that the material blown away is predominantly hot and
chemically enriched. On the other hand, simulation times of
about 400 Myr, as we have used here, imply a substantial
shearing of the computational box due to differential Galac-
tic rotation, which was omitted for convenience. We have
performed test simulations in the past, employing shearing
boxes, which have shown that the effect on the distribution
and dynamics of ISM gas is comparatively small, although
the influence, like in most other simulations, of a weak spiral
density shock wave has been neglected. In molecular clouds
or very dense regions, however, shear may work against
large gravitational collapse (e.g. Kim and Ostriker 2017).
The basic results of a dynamic ISM evolution, according
to hydrodynamic (HD) numerical simulations can be seen in
Fig. 1, which represents the density distribution at a maxi-
mum resolution of 0.25 pc in the Galactic midplane of an
HD run after 300 Myr of evolution time, for a Galactic SN
rate. The colour coding of the logarithmic density distri-
Fig. 1 Gas density distribution for a Galactic supernova rate after 300
Myr of evolution in the midplane of the Galaxy (logarithmic represen-
tation, colour coding with black being the highest and white the lowest
threshold densities), with the coordinate origin being in the Sun. The
highest spatial resolution of the simulation has been 0.25 pc
bution shows highest densities in black (>100 cm−3) and
lowest in white (<10−3.5cm−3), ranging from 10−4.4to
600 cm−3.
One remarkable feature is the overwhelming amount
of small scale frothy filamentary features, which are a
conspicuous sign of a high level of turbulence, coupling
a vast amount of scales. In incompressible turbulence
(see Sect. 2.1.2) this range is directly proportional to the
Reynolds number Re ∼uL/νm,i.e.η=LRe3/4, where η
is the viscous length scale at which molecular dissipation
with kinematic viscosity νmoccurs, with uand Lbeing
the bulk velocity and length scale of the flow, respectively.
The Reynolds number in the ISM has been estimated to be
Re ∼106−107(Elmegreen and Scalo 2004).
2.1.2 Scaling laws for turbulence
Since the ISM flows have quite often large Mach num-
bers, turbulence is compressible, and the well-known Kol-
mogorov (1941) scaling law, where the spectral energy den-
sity is given by E(k) =C2/3k−5/3, has to be modified; Cis
supposed to be a universal constant in Kolmogorov’s theory,
and is the energy dissipation rate. Here it was assumed that
turbulence on scales l(inertial range, ηlL) is statis-
tically homogeneous and isotropic, and that the dissipation
rate is constant, as one would assume for a scale invariant
process. These scaling laws work surprisingly well in Na-
ture, although it has been known for a long time that the pro-
cess as such cannot be universal and statistically self-similar
Loading more pages...