Non-equilibrium ionisation plasmas in the interstellar medium

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 determined by atomic physical time scales and reaction rates. To calculate spectra is further complicated by gas dynamical 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 ionisation equilibrium for optically thin plasmas. Computing power and new algorithms for high performance computing have made it possible to trace the dynamical and thermal evolution of a sufficiently large section of interstellar space over an appreciable time scale to derive characteristic 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 interstellar plasma simulations, the physical processes which drive the temporal and spatial evolution, and present high resolution numerical simulations, including time-dependent cooling, which further our understanding of the state and evolution of interstellar (magnetised) plasmas. We also discuss briefly the rôle of cosmic rays and their interaction with the plasma.


Introduction
The interstellar medium (ISM) in galaxies is a complex multi-component system composed of a highly compressible 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 nonlinear. Although numerical simulations on parallel computers 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 tempting to throw all the components and their respective interactions on a supercomputer and predict the temporal and spatial 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 predictive power of such an endeavour may be disappointingly small, because one has to sacrifice spatial resolution and resort to "recipes" rather than treating all processes and interactions of components by their appropriate physical equations. 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 compressed and heated up (and if it is magnetised, magnetohydrodynamic (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 derive from a large amount of free energy available in the system 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 supersonic. 1 Therefore we believe that progress at this stage of research is mainly achieved by analysing the main components and their governing physical processes -this again is debatable -on a reasonably sized computational grid. This computational box should be large enough, especially perpendicular to the disk, to prevent the material from leaving the computational domain during the simulation time, before 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 providing 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. n 2 = n 2 . This is particularly important in the case of radiative cooling, which depends on n 2 locally. Too low spatial resolution therefore underestimates energy losses by escaping photons in an optically thin plasma.
Before devising an elaborate model, it is important to sort out what we believe are the most important physical processes 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 energy 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 1979Heiles , 1984 than expected from a single explosion. The theoretical description of superbubbles (e.g. McCray and Kafatos 1987;Mac Low and McCray 1988) gave a plausible explanation. For some time it was discussed that these superbubbles (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 simulations suggested. However, these were two-dimensional, 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 scenario for the SB expansion. As an aside, it is well-known how difficult plasma confinement is in terrestrial fusion devices 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 reconnection (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 perpendicular 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 v esc = √ −2 0 , where 0 equals 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) by a non-linear interaction between cosmic rays (CRs), their selfexcited 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 models in the past), because particles are accelerated in SN generated 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 generation of waves that tend to remove this anisotropy. Waves 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 momentum from the CRs to the plasma via the waves as a mediator, 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 simulations 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 simulations of the ISM were carried out, taking SNe and stellar 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 observed. Since turbulence is an integral part of such simulations, 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 z of ±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 (2004Breitschwerdt ( , 2005 have carried out hydrodynamical and MHD simulations with 1 × 1 kpc in the disk (x − y plane) and ±10 kpc in the vertical z-direction (in later simulations ±15 kpc, cf. de Avillez and Breitschwerdt (2012a)) with a spatial resolution of 0.5 and 1 pc, respectively. Although this is still above the critical wavelength 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 become too large. Simulations at decreasing scales (increasing resolution) have shown that quantities like volume filling 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 T is the temperature and Z the 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 ionisations 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. 2 This allows observers 2 In 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 feedback on the ionisation structure by reabsorption. In addition, for densities below 10 4 cm −3 , collisional de-excitation does not play a significant rôle, so that the line emission is determined by collisional excitation (coronal approximation). For example the detection of the ubiquitous OVI absorption line is commonly identified with a plasma at ∼ 3 × 10 5 K. However, there is a caveat. The assumption of CIE is only approximately valid at temperatures above 10 6 K (Shapiro and Moore 1976). For lower temperatures, and this concerns most of the ISM in mass and volume, ionisation and recombination rates are different, simply due to different temperature dependent cross-sections. This has been realized long ago by Shapiro and Moore (1976), who studied the cooling of a 10 6 K shock heated gas, showing the progressive deviation from CIE. Shapiro and Moore (1977), also analyzed the X-ray emission in solar flare coronal non-equilibrium ionisation (NEI) plasmas, impulsively heated to 10 8 K. 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 collisional heating or cooling of the electrons, which then subsequently share energy and momentum with the ions. Therefore 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. Conversely, fast adiabatic expansion of a hot plasma bears the ionisation signature of the hot state, although the gas can be cold (see Schmutzler 1994, 1999). As we will show, this can result in the emission of soft X-rays at temperatures of only 10 4 K (see Sect. 4). The fundamental difference to CIE is now, that the ionisation state varies in space and time, and so does the cooling function in nonequilibrium 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 techniques. 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 particle hydrodynamics (SPH) has become very fashionable. Already 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. method, and its main feature is a discretization in mass instead of space, it is meshfree and can therefore be applied to free flows, and it can be easily parallelized. On the downside 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 combine 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 cosmological 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 retaining 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 simulated a whole galactic disk in order to analyse the distribution of the gas in general, and the formation of molecular 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) take an even more ambitious 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. 2 describes the numerical setup and modelling of ISM plasmas 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. timedependent cooling, is outlined in detail, and observational consequences are presented. Section 5 deals with the importance of cosmic rays in the ISM and their interaction with magnetic fields and the gas, with galactic winds being just one important example. In Sect. 6 the Local Bubble, 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 distribution of highly ionised species and electrons there. Section 7 concludes this review and summarizes the results.

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 volume, which, according to a user defined refinement criterion, e.g. density or pressure gradient thresholds, is overlaid 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 problems, 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 resolution, but cover only a comparatively small volume of the grid.

Hydrodynamic simulations
The physical and numerical setup of high resolution 3D mesoscopic hydrodynamic simulations with AMR is as follows (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 z is included according to Wolfire et al. (1995), being constant parallel to the Galactic disk, initially balancing radiative cooling at 9000 K; (ii) the cooling function has been tabulated and taken from Dalgarno and McCray (1972) in earlier versions, with the ionisation fraction being 0.1 at temperatures below 10 4 K 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 density and temperature thresholds of 10 cm −3 and 100 K, respectively, 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 velocity 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 selfgravity of the plasma and thermal conduction, although we still believe that turbulent mixing is the more efficient transport process. The boundary conditions are outflow at the top and bottom of the box, and periodic at the sides (perpendicular to the disk), thus avoiding a depletion of matter within the long simulation times. This of course introduces an artefact, e.g. in the disk flows due to the periodic boundary conditions, 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.

General results
An important finding of HD ISM simulations is the establishment of the Galactic fountain cycle, which operates on a time scale of about τ gf ∼ 150 Myr, implying that simulations of smaller time spans cannot be a sensible representation 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 z c ∼ 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 ∼ 10 6 K, as seen in X-rays, the speed of sound is c s ∼ 130 km/s, and thus τ gf ∼ z c /c s ≈ 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 Galactic 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 maximum 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 representation, 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.5 cm −3 ), ranging from 10 −4.4 to 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. η = LRe 3/4 , where η is the viscous length scale at which molecular dissipation with kinematic viscosity ν m occurs, with u and L being the bulk velocity and length scale of the flow, respectively. The Reynolds number in the ISM has been estimated to be Re ∼ 10 6 − 10 7 (Elmegreen and Scalo 2004).

Scaling laws for turbulence
Since the ISM flows have quite often large Mach numbers, turbulence is compressible, and the well-known Kolmogorov (1941) scaling law, where the spectral energy density is given by E(k) = C 2/3 k −5/3 , has to be modified; C is 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, η l L) is statistically 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 Nature, although it has been known for a long time that the process as such cannot be universal and statistically self-similar in the inertial range. If it were, the scale invariance would imply that the velocity increment δ u( r) = u( x + r) − u( x) for a scaling of λ, i.e. r → λ r, δ u( λr) would follow the same statistical distribution as λ α δ u( r), where the exponent α is the same for all scales. Hence the statistical moments of the velocity increment should obey (δ u( r)) p = C p ( r) p/3 , where p is the order of the moment, ... is the statistical (ensemble) average and the expression itself is called the structure function of order p. Now the success of the Kolmogorov scaling rests on the fact that the spectral energy density E(k) is related to the 2nd order structure function (p = 2), where deviations from scale invariance have been found to be small. However, higher order structure functions, such as skewness or kurtosis, show significant deviations from self-similar scaling, and therefore question the validity, in particular the universality of the constant C, of Kolmogorov's assumptions. A fundamental characteristic of a turbulent system is the phenomenon of intermittency, which is a sudden change between outburst and breakdown of fluctuations in the fluid variables. Mathematically, this points to a departure from a Gaussian distribution, as higher order moments do not converge: turbulence on the microscales is not isotropic and space filling as envisaged by Kolmogorov, since the increase of higher order moments is a sign of spatial concentration. More recently, higher order moments have also been considered in supersonic turbulence in molecular clouds by Padoan et al. (2016).
In ISM, as has been already mentioned, the scaling must be different, because shocks can dissipate energy at any scale. Hence it is remarkable that a scaling law does still exist. Based on an old idea of von Weizsäcker (1951) and von Hoerner (1951), Fleck (1996) has pictured a compressible non-self-gravitating flow as a hierarchy of mass densities ρ ν according to for successive levels of ν and ν − 1, with f being the volume filling factor, and 0 ≤ β ≤ 1 is a measure for the compression (β = 0 recovers the incompressible case, and β = 1 stands for isotropic compression). It can then be shown that, if instead of the energy rate per unit mass that is constant during transfer to smaller scales, one uses the energy density rate (energy rate per unit volume) V = ρ = ρv 3 /l, the scaling law re-emerges for the density weighted velocity u ≡ ρ 1/3 v with u 3 /l = const.. This has been confirmed by 3D numerical simulations (Kritsuk et al. 2007) for an isothermal gas up to Mach numbers 6. They also have shown that the probability density function (pdf) is lognormal over several decades (see also Passot and Vázquez-Semadeni 1998, and Sect. 2.1.4).
We have demonstrated (de Avillez and Breitschwerdt 2007) that for 3D ISM turbulence, which is not artificially driven, but is due to SNe and stellar winds stirring up the ISM, and without restrictions of the equation of state, including magnetic fields and covering scales from subparsec to kiloparsec, the most dissipative structures have Hausdorff dimensions D of the most intermittent structures varying from D = 1.98 − 2.2. This corresponds to shock surfaces in the simplest case. We could also show that the scaling of the structure functions is consistent with the log-Poisson statistics, emphasising once again the rôle of intermittency. In addition, the flattening of the second order structure function for small wave numbers k, or large scales l, constrains the eddy correlation length for the energy injection on the largest scales, and results in an integral (or outer) scale for turbulence of ∼ 75 pc for a SN driven ISM. This is, loosely speaking, the scale on which the interstellar mixing spoon operates, and represents the average size of closed bubble structures before break-up due to density and pressure gradients in the ambient medium. Of course turbulence can be generated in the ISM also on larger or smaller scales, depending on the sources (see Elmegreen and Scalo 2004), but SN driving is the dominant process (Mac Low and Klessen 2004).

Volume filling factors
Not surprisingly, the volume filling factor (see Fig. 2) for the hot gas (defined as gas in the temperature range T > 10 5.5 K) is of the order of f V ,h ∼ 0.2 for a Galactic SN rate, much lower than in the MO77 model, which predicted f V ,h ∼ 0.7 − 0.8 for an average temperature of T = 10 5.7 (McKee and Ostriker 1977). Even increasing the SN rate by a factor of 4 only increases f V ,h to 0.28. This is due to the fact that hot SNR generated gas with its overpressure with respect to the ambient medium has a strong tendency to follow the density gradient and expand perpendicular to the Galactic disk into the halo. But in which temperature range can the bulk of the interstellar gas be found? About 64% of it resides in the cool and warm medium with temperatures of 10 3 − 10 4 K and 10 4 − 10 5.7 K, respectively, with a substantial fraction of it being in classically thermally unstable regimes. For example, in agreement with observations by Heiles (2001) and Heiles and Troland (2003), who find a fraction of 44% in the unstable temperature range 500 < T < 1500 K, our simulations give 55%, again in disagreement with the MO77 model. This is consistent with the fact that more than 50% of the gas is rather cold (see Ferrière 2001). The effect of additional heating processes like e.g. photo-electric heating by dust grains (Hill et al. 2018) also influences the volume filling factors and the stability of the warm neutral medium (WNM) and cold neutral medium (CNM) phases. Although the supernova locations and explosion times are not correlated with gas density, which is unphysical, the authors wanted to isolate the diffuse heating rate from the supernova history. They find that over a wide range of parameters the gas resides in the two-phase medium (CNM and WNM), and that the volume of the hot gas depends on the diffuse heating rate.
In our view, the ISM should not be viewed as a SN regulated medium which is in distinct phases, which are roughly in pressure equilibrium, but rather as a SN driven turbulent medium, in which the gas is cycled through different density and temperature regimes, with broad probability density functions (pdfs) in density and pressure.

Probability density functions (pdfs)
The difference with respect to classical ISM models, in which the gas resides in distinct stable phases, which exchange mass, being roughly in pressure equilibrium, e.g. by evaporation or condensation, is best seen in analysing the probability density functions (pdfs). They represent, as shown here, the occupation fractions of densities (Fig. 3), pressures (Fig. 4) and temperatures (Fig. 5) in their respective regimes. The most conspicuous feature of all the pdfs is their broad distributions, spanning six orders of magnitude or more. This contrasts with the general picture of a quiescent ISM with small fluctuations about average values in density, pressure or temperature. If this were the case, the pdfs would be much narrower and sharply peaked, instead of having, as is seen in the figures, broad, almost plateau-like maxima. This behaviour can be most easily attributed to turbulence, coupling a large range of scales, channelling kinetic Probability density function for the electron density distribution in the midplane of the Galaxy, with its flat maximum spanning more than three orders of magnitude for an evolution time scale of 300 Myr. This is in contrast to the classical distribution of the gas into distinct phases; here (and in the following) x, σ/σ Gal and E 51 represent the maximum spatial resolution, the ratio of SN rate to the standard Galatic value, and the SN energy input in units of 10 51 erg, respectively energy through a cascade by a hierarchy of eddies. This is done by extracting energy from the mean flow by steadily increasing vorticity. Mathematically, this can be seen directly by rewriting the Navier-Stokes equation for an incompressible and Newtonian fluid in terms of vorticity, ω = ∇ × u, which can be interpreted as a measure of the local angular rotation rate of a fluid element: with ν m being the kinematic viscosity of the flow. The change of vorticity is due to two different processes. The first one (first term on the right hand side) describes the change of the moment of inertia by stretching a fluid element, whereas the second term represents the torque due to viscous stresses, i.e. spinning up (increasing vorticity) or slowing down (decreasing vorticity) of a fluid element. The increase of vorticity due to advection and diffusion leads to an ever increasing coverage of space with thinning out flux tubes down to the viscous scale, where the kinetic energy is dissipated due to molecular viscosity. Hence such a situation can only be maintained over a substantial amount of time if there is a source of continuous energy input on the largest (integral) scale, which, in the case of the ISM, is due to SN explosions (Mac Low and Klessen 2004).

Magnetoydrodynamic simulations
HD simulations are our basic tool of studying physical processes in the ISM. However, the ISM is a magnetised fluid, and hence, including magnetic fields is an important step towards greater realism. The simplest way to achieve this is Probability density function for the thermal pressure distribution in the midplane of the Galaxy, with its flat maximum spanning more than one order of magnitude, thus contradicting the standard assumption of ISM evolution in pressure equilibrium

Fig. 5
Probability density function for the thermal plasma temperature distribution in the midplane of the Galaxy, spanning several orders of magnitude. Its broadness contradicts the classical assumption of a distribution into distinct phases. The peak at T ≈ 10 7.5 K is due to the impulsive SN shock heating at the time of the snapshot, taken at 300 Myr after the start of the simulation by ideal (sufficiently large conductivity, neglecting resistivity) magnetohydrodynamics (MHD), which neglects static electric fields and displacement currents (being of order v/c compared to magnetic fields, where c is the speed of light), leaving us with Ampère's law. Electric fields are induced by fluid motions in a magnetised plasma, since these fields only vanish in a co-moving frame, and are therefore given The equation of motion is supplemented with terms of magnetic pressure forces and stresses. And finally, the magnetic field has to be divergence free, ∇ · B = 0, which in numerical schemes poses some prob-lems; we have used the method suggested by Balsara (2001) for MHD simulations with AMR. The box has been set up as before, i.e. with periodic boundary conditions on the faces of the vertical direction, and outflow boundary conditions at the top and bottom. The size is 1 kpc ×1 kpc in the x − y-plane and ±10 kpc in the z-direction perpendicular to the disk. The magnetic field has a uniform and a random component given by B u = (B u,0 (n(z)/n 0 ) 1/2 , 0, 0) and δ B = 0 (initially, building up over time), respectively, where B u,0 = 3 μG is the field strength of the uniform component at the Galactic midplane near the solar circle, n(z) is the number density of the gas varying in the vertical direction with a midplane value of n 0 = 1 cm −3 . The highest grid resolution in the MHD runs has been 1.25 pc, triggered again on the fly by AMR threshold values of density and pressure gradients.

Results
It is probably not fortuitous that a random field component, initially only represented by random noise, builds up to a field strength of 3.2 μG, comparable to the average value of the uniform field component of 3.1 μG. Like in the HD simulations, the initial gas distribution is not stable due to radiative cooling, which sets in immediately, whereas the build-up of thermal pressure -ram and turbulent pressures are zero initially -in order to stabilise the system against gravitational collapse needs time, typically of the order of a sound crossing time, which, for a vertical distance of 1 kpc, takes about 10 Myr for a 10 6 K gas. In contrast to HD runs, the initial collapse is slowed down due to the magnetic pressure of an initially disk parallel field. By the same token, the outflow into the halo is slowed down due to magnetic tension forces, which, however, can only delay but not suppress the upward flow. This behaviour can be attributed to the formidable field tangling, as can be seen in Fig. 6. Note that there are also large scale filaments, spanning several hundred parsecs across. In regions with higher field strength, like e.g. in the Galactic Centre, there will be naturally less field tangling and ordered structures should be more common there. Moreover, dynamo simulations show that field amplification is mainly due to a small scale field dynamo (Gent et al. 2020) and less to field tangling. Although flux freezing is always a good approximation in the ISM, as a low degree of ionisation is sufficient for a reasonable coupling of plasma to the field, the chaotic nature of the magnetic field, i.e. the random direction, allows for gas to flow into the halo without much hindrance. One might also speculate if the strongly tangled field provides lots of close anti-parallel field line encounters and therefore promotes turbulent magnetic reconnection, as has been suggested some time ago (Lazarian and Vishniac 1999;Lazarian et al. 2020).
A total magnetic field strength of about 4.46 μG is not sufficient to drastically alter the volume filling factors, pdfs and the ISM structure as a whole when compared to the HD simulations. The most notable differences are in the volume filling factors for the thermally unstable regimes of 200 < 10 3.9 K and 10 4.2 < 10 5.5 K (although added up, they have similar values). This is due to the additional magnetic pressure that reduces compression on average and thereby also affects radiative cooling, resulting in values of f V ,HD ∼ 0.46 versus f V ,MHD ∼ 0.29, as well as f V ,HD ∼ 0.22 versus f V ,MHD ∼ 0.33, respectively. This is consistent with the fact that ram pressure dominates thermal and magnetic pressure, except for the hottest and coldest regions (de Avillez and Breitschwerdt 2005), implying that turbulence is the dominant physical process that governs the ISM dynamics. Nevertheless, the flow in the magnetised case looks somewhat smoother on small scales, as one would expect from the increase of magnetic tension forces for smaller curvature radii when bending the field lines. Likewise the density distribution perpendicular to the disk (left panel of Fig. 6) shows more coherent structures on larger scales, which can be attributed to the stretching of field lines, and a large scale vertical field component underlying the small scale field tangling.

Non-equilibrium simulations of the ISM
In numerical simulations, radiative cooling is usually treated in the coronal approximation for an optically thin gas, which holds as long as the photons, which are generated through collisional ionisation and subsequent radiative de-excitation, can escape from the plasma. Otherwise, additional equations for radiation transport, e.g. in the form of ray tracing, have to be implemented into the code. The effect of radiative cooling can then be summarized in the energy equation as a loss term, which depends on the plasma density, ρ, the temperature T , and the metallicity Z. This implicitly assumes that the ionisation structure of the plasma is decoupled of its dynamics -which in reality is not the case. Schmutzler and Tscharnuter (1993) have shown that an element of fluid cooling down in a fully ionised plasma at T = 10 9 K, both isochorically and isobarically, will progressively deviate from collisional ionisation equilibrium (CIE) -as a rule of thumb one can state that plasmas below 10 6 K will be out of ionisation equilibrium, because in general the cooling time scale will be shorter than the corresponding recombination time scale. This has been recognized early on in exploratory works by Kafatos (1973) and Shapiro and Moore (1976), and was numerically taken into account by Cox and Anderson (1982), who, in a simple dynamical model, followed the expansion and the ionisation structure of a SN blast wave. They tracked 30 parcels of fluid starting from CIE to trace their ionisation structure in order to eventually derive the X-ray emission of the SN remnant (SNR).
In a series of papers Innes et al. (1987a,b,c) calculated the ionisation structure and line emission behind a shock, including a photoionised and heated precursor. In an iterative process between dynamics and distribution of ions, they calculated a self-consistent solution. When a shock overruns a cool plasma, collisional heating of the electrons is very fast and the ionisation lags behind, generating an "underionised" plasma, in which the distribution of the ionisation stages still mimics a cold plasma. Breitschwerdt and Schmutzler (1994) calculated the opposite case of an "overionised" plasma, in which due to fast adiabatic expansion, recombination lags behind the fast cooling electrons. Such a plasma resembles observationally a hot X-ray emitting diffuse gas, although with distinct recombination edges, but currently not resolvable with X-ray spectrometers due to a too low photon flux. They were able to reproduce the X-ray emission, in particular the C-and M-band (0.3 − 0.6 keV) in the Wisconsin Survey (McCammon et al. 1983;McCammon and Sanders 1990) as a plasma characterized by delayed recombination due to fast adiabatic expansion of a SB breaking out of a molecular cloud. It is thus possible to generate X-ray emission at temperatures as low as 4 × 10 4 K. For a more recent interpretation of the origin of the Local Bubble see Sect. 6.1.
One of the fundamental differences between a CIE and NEI plasma is, that in NEI the thermal history and its thermodynamic path are crucial for the plasma emission and cooling. As an example, take two parcels of fluid, which start at the same temperature, density, metallicity etc., i.e. are thermodynamical twins with equal ionisation structure, but evolve differently, e.g. one being shocked, the other one expanded. Even if they returned after some time to the same end state, their photon emission and radiative cooling would have been different, and hence their emergent spectra. In between the two states many irreversible processes have occurred, which have in turn changed the dynamics as well. This may sound like a nightmare to observers -and theoreticians who want to build as simple as possible models as well -when interpreting spectra. But one should keep in mind, that interstellar NEI spectra therefore contain much more information on the history of the plasma, as opposed to CIE, which ignores the history completely. This has been shown (de Avillez and Breitschwerdt 2012a) by taking into account electron impact ionisation, inner-shell excitation autoionisation, radiative and dielectronic recombination, chargeexchange reactions, continuum (bremsstrahlung and twophoton), free-bound and line (permitted, semi-forbidden and forbidden) emission in the wavelength range 1 Å-610 μ (see Sect. 4.1 for details).
The resulting cooling function will be lower than in CIE, because much of the radiative energy is still stored in high ionisation stages. Clearly, these time-dependent processes should have noticeable effects on the energy budget of the ISM. Since the cooling function is now also a function of space and time (and metallicity Z), the time-dependent ionisation structure needs to be calculated simultaneously and self-consistently with the HD/MHD evolution of the plasma.

Time-dependent cooling
As has been emphasised in the last section, the interstellar cooling function is in general dependent on both space and time (and metallicity), i.e. = ( x, t; Z). If the cooling is dominated by delayed recombination, there is still an appreciable amount of energy stored in highly ionised species as a form of "potential" energy, which is only gradually released at a longer time scale. Cooling down from higher to lower temperatures, leads to a smaller energy loss than would correspond to the already lower electron temperature of a CIE plasma. Therefore the cooling function lies below the CIE standard, and can differ from it by up to two orders of magnitude (s. Fig. 9), depending on the thermodynamic path each parcel of fluid takes.
In the general ISM, which is characterized by many energetic and impulsive processes, such as SNe, stellar winds, bursts of radiation etc., the assumption of CIE is more than dubious, since the majority of the gas by mass is below 10 6 K. Strictly speaking, only plasma with very low abundances can be approximated by CIE, because here the timescale for radiative cooling may be larger than for recombination.
The lack of detailed balancing inevitably drives the plasma out of ionisation equilibrium, even if it was so initially, simply because collisional ionisation (which does not need to obey the quantum-mechanical selection rules) and radiative recombination are not inverse processes, with the photons eventually leaving the system. This can be seen when realizing that collisional ionisation is a cooling process for the plasma, because the energy of thermal electrons is used for ionisation of atoms and ions, while recombination is also a cooling process, because the resulting photon energy is lost in optically thin plasmas. To maintain ionisation equilibrium would mean to balance cooling by a specially designed heating function.
As has been argued above, fast compression (e.g. by shocks) and expansion (e.g. by rarefaction waves) are events which lead to "underionised" and "overionised" plasmas, respectively, implying significant deviations from the CIE cooling function. The only way to keep reliably track of the energy loss is by following the ionisation structure selfconsistently with the dynamics. Let us illustrate this by considering a parcel of fluid, which has been expanded suddenly, e.g. because a hot plasma breaks out of a bubble, something occurring regularly in ISM simulations, as these bubbles frequently encounter steep density gradients. The internal energy has to be redistributed over a larger volume, and the gas additionally performs work on the surrounding medium. Therefore the internal energy and thus the temperature of the electrons, and with some delay of the ions, drops significantly, while the collisional ionisation and recombination are affected indirectly due to the lower number density of electrons and ions. They consequently lag behind in time, simply because the atomic timescales are much larger than the dynamical timescale. This reduces the cooling rate, and affects the energy equation, thereby changing the dynamics of the flow, thus closing a feedback loop.

X-ray emission at low temperatures
Current fit packages for X-ray spectroscopy already contain routines for NEI plasmas, e.g. to fit lines from SNRs. They are often simple approximations using (for a given metallicity) the ionisation parameter and a constant temperature (e.g. in the routine XSPEC) of a plasma to improve on CIE spectral fitting. However, they cannot capture the full nonlinearity of the problem and ignore the dynamical state of the plasma, hence are not physical but fit models. Such a detailed treatment can only be achieved by performing a selfconsistent thermal and dynamical numerical simulation of the plasma.
It is worth emphasising that for a sufficiently high spectral resolution, X-ray spectra will also contain the information of the dynamical history of the plasma, which has hitherto been ignored both in models and in data analysis. Since this is an inverse problem, it is difficult, and to some extent also not unique, to reconstruct the plasma history. Nevertheless, it would be a step forward to start with simple NEI simulations (not fit models) to extract information with respect to the deviation from ionisation equilibrium. To summarise the physical problem again: a SN driven ISM plasma evolves in a highly non-linear way with a high level of turbulence. This occurs on timescales which are independent of atomic timescales, such as ionisation or recombination rates, which differ for species and their respective ionisation stages. Line cooling depends on the recombination rate and the local density for a given local temperature, which changes the internal energy content, and hence the density and pressure, and thereby the flow velocity, which in turn alters the recombination rate and so forth.
To illustrate this effect, we have performed NEI simulations of an ISM plasma as described in Sects. 2.1 and 2.2, but now including the time-dependent ionisation and recombination rate equations (de Avillez and Breitschwerdt 2012a). The different plasma histories for the patches denoted in Fig. 7 lead indeed to individually different cooling functions, which can differ by an order of magnitude in the density normalised cooling rate. This by itself should be a sufficient reason not to neglect the time-dependence of radiative cooling. But what is even more striking is how the different cooling histories have markedly differing spectral characteristics, as can be seen in Fig. 8. The left (right) Fig. 7 A temperature regime of 10 4.1 − 10 6.1 K extraction from a NEI HD ISM simulation. The regions labelled F -J represent an initial temperature of 10 5.5 K and the labels K -O an initial temperature of 10 6 K, from which we follow the subsequent thermal and dynamical evolution, as the plasma is cooling down. The highest spatial resolution of the simulation has been 0.5 pc panel shows the emission starting at an initial temperature of T = 10 5.5 K (T = 10 6 ), arriving at later times at temperatures of 10 4.8 K, 10 4.4 K and 10 4.0 K. The corresponding CIE emission is shown as a dashed line in the lower right corner for the highest temperature, while it disappears from the graphs completely at lower temperatures. As a visible comparison the CIE temperature of a 10 6.2 K CIE gas (orange line) and the NEI cooling of a fully ionised static plasma (black solid line) are shown instead.
In both panels the effect of NEI line emission due to delayed recombination is obvious, and is even more stressed by exhibiting X-ray emission at a temperature of 10 4 K! This should be a caveat for observers to identify X-ray emission from diffuse plasmas automatically with ∼ 10 6 K gas. In summary, we reiterate our contention that radiative cooling of an optically thin plasma is characterised both by spatial and temporal dependence. Thus the application of static cooling functions in simulations and spectral fitting depending only on temperature and metallicity is a gross oversimplification and therefore bears the potential of significant misinterpretation of the plasma state and evolution in general.

Distribution of free electrons in NEI plasmas
The distribution of free electrons in the Galactic disk as derived from the observation of the average volume densities of the diffuse gas and the pdfs (Berkhuijsen and Fletcher 2008) reflect the NEI character of the turbulent ISM plasma. These densities were obtained by studying the dispersion measures and HI column densities in the direction of known pulsars and stars. A snapshot at a simulation time of 350 Myr of a NEI electron density distribution in the Galactic midplane, covering a range of four orders of magnitude, can be seen in Fig. 10, which again shows a frothy and fractal structure down to the smallest resolution scale of 0.25 pc (cf. de Avillez et al. 2012).
We have simulated the pdf for the electron density distribution (Fig. 11), as well as the volume weighted pdf (Fig. 12) for simulation times between 420 and 450 Myr. It can be clearly seen that the distribution is rather broad, in fact lognormal (which points towards the central limit theorem, with a large number of randomly distributed variables) and in agreement with the results of Berkhuijsen and Fletcher (2008), having a distinct hump at high densities, which may be a sign of delayed recombination in the plasma. Measuring the electron column density along lines of sight of lengths up to 4000 pc (cf. Fig. 13) from a given vantage point can be useful to determine not only the accuracy of the simulation results with regard to observations, but can also provide insight into the mean electron density in the Galactic disk and, in particular, in the Galactic midplane.
The top panel of Fig. 13 displays the minimum, average and maximum of the time average over the period 400- The simulated dispersion measures fall within the observed values and show a constant dispersion with respect to the mean along the lines of sight, indicating that the electrons are distributed in clumps rather than in clouds. The mean electron density has a large variation up to 500 pc from the vantage point, and tends to 0.03 cm −3 , as more and more ISM is sampled, consistent with the observations.

Cosmic rays
Cosmic rays (CRs) fill the Galactic disk and the halo and represent the high energy component of the ISM. It has been known for a long time that locally there exists a rough energy density equipartition of the order of ∼ 1 eV/cm 3 between thermal plasma, magnetic fields, radiation fields, turbulence and CRs. This is not a coincidence since all these components of the ISM strongly interact with each other, thus driving the system to a minimum energy configuration. The energy range of CRs is enormous, from below 100 MeV to beyond 10 20 eV, so that it is not surprising Fig. 8 Free-bound emission (colour curves) by the cooling gas at sites having an initial temperature of 10 5.5 K (similar to sites F -J in Fig. 7; left-column panels here), and 10 6 K (similar to sites K -O; right-column panels here), respectively. Each panel shows a snapshot of the cooling plasma at the temperature indicated in the upper right corner. Solid and dashed black lines represent the free-bound emission from NEI and CIE static plasmas, respectively, at the temperatures shown in each panel. The orange curve represents the emission from a CIE plasma at 10 6.2 K for comparison. Recombining ionisation stages responsible for the edges are indicated. The free-bound emission due to CIE evolution is lower than 10 −33 erg cm 3 s −1 Å −1 for T ≤ 10 4.4 K and hence is not displayed (from de Avillez and Breitschwerdt 2012a, ©AAS. Reproduced with permission) that there exists a variety of possible sources. The energy spectrum, which is a power law over almost 12 decades, with a bend-over at low energies due to solar modulation, a break ("knee") at ∼ 10 15 eV, and another break ("ankle") at ∼ 10 18 eV, suggests a transition form solar to interstellar to extragalactic CRs. The latter, according to the standard argument, occurs because the proton gyroradius at 10 18 eV becomes comparable to the radius of the Galaxy. CRs are a misnomer, as they consist of charged particles, about 90% protons, 9% helium plus heavier nuclei, and 1% electrons. As such, they strongly interact with magnetic fields, and no conclusion can be drawn for the site of their origin from their arrival directions, (except for the ultra-high energy CRs (UHECRs) which are poorly scattered), because pitch angle scattering randomizes their paths. Thus CR propagation can be viewed as a Lévy or Wiener type stochastic process, in which a particle performs a large number of successive displacements, which are random and independent (and similar even over different time intervals). However, this so-called CR diffusion is superposed by advection, as the CR pressure gradient acts preferentially out of the disk, thus promoting CR escape. The CR residence timescale has been measured from 10 Be isotope clocks to be τ esc = 15 +7 −4 Myr (Simpson and Garcia-Munoz 1988). Another important measurement is the "grammage", i.e. the average matter traversed by a particle, accompanied by nuclear interactions and spallation of heavy nuclei, which for a 1 GeV proton has been observed to be x obs ∼ 5 − 10 g/cm 2 . This implies that for an average ISM number density of 1 particle per cm 3 in the Galactic disk, the time the particle would spend there would only be a few Myr, contrary to observations of primary to secondary ratios. Hence, the particle "sees" an average density of 0.24±0.07 atoms cm −3 (Simpson and Garcia-Munoz 1988), i.e. by a factor 3 to 4 lower, and therefore spends most of its time in the lower density Galactic halo.
In order to explain the power law differential energy spectrum of CRs, a natural process would be a stochastic one in which particles have an energy dependent loss. This would ensure that for higher energies, less particles stay in the accelerator, as it is the case in the first order Fermi or diffusive shock acceleration (DSA). Since particles leave the Galaxy at a rate τ −1 esc , CRs have to be constantly re- Fig. 9 Total emission (colour curves, each with a different cooling history) by the cooling gas at sites having an initial temperature of 10 5.5 K (similar to sites F -J in Fig. 7; left figure) and 10 6 K (similar to sites K -O in Fig. 7; right figure), respectively. Solid dotted and dashed black curves represent the free-bound emission from CIE and NEI static plasmas, respectively, at the temperatures shown in each figure. Note that the energy loss rate of a dynamical NEI plasma (colour curves) cooling down is lower than that of a static CIE and NEI plasma Fig. 10 Electron density distribution in the midplane of the Galaxy (logarithmic representation, colour coding with red being the highest and blue the lowest density) at 350 Myr of disk evolution, with the coordinate origin being in the Sun. The highest spatial resolution of the simulation has been 0.25 pc plenished. The energy loss rate is F CR ∼ CR V conf /τ esc ∼ 4.7 10 40 erg/s, where V conf ∼ 1.8 10 67 cm 3 is the CR confinement volume. 3 This requires a powerful source, and leaves essentially SNRs, which can provide a total energy of ∼ 6.3 10 41 erg/s, for a SN rate of 2 per century having an average energy of 10 51 erg per explosion. This leads to a reasonable value of about 10% of the available SN energy to be converted into CR energy. Hence, DSA at SNRs can provide the necessary acceleration. In brief, suprathermal particles are deflected and scattered back and forth across the shock by magnetohydrodynamic waves, which are excited by the particles themselves as they stream along magnetic field lines by the so-called streaming instability (Lerche 1967;Kulsrud and Pearce 1969;Skilling 1975). In essence, scattering off waves isotropizes the particle distribution in the wave frame, but particle streaming induces a small anisotropy, which is removed by the resonant generation of waves, preferentially exciting waves with wavelengths of the order of the particle gyroradius. As the particle moves across the shock it "sees" a converging fluid, and gains energy of the order of the difference between up-and downstream velocity, v/c (where c is the speed of light), which is of first order. Before being convected downstream, the particle has a probability to stay in the acceleration process by being back-scattered across the shock, thus going through the same process over and over again, moving up secularly in energy. However, as the energy is increasing, the number of scatterers (waves) decreases, and the particle diffusion coefficient increases, so that the probability to escape is proportional to its energy. This naturally leads to a power law distribution in energy or momentum. One of the big successes of DSA is that it could predict correctly the power law index α of E −α to be α 2 in energy for particles with energies E ≤ 10 15 eV, where the difference to α ≈ 2.7 in observation is due to energy dependent diffusion ∝ E −0.7 . Chandra X-ray observations of Tycho's supernova remnant (SN1572) seem to confirm the DSA theory for particle acceleration up to the knee, according to the inferred high emissivity spacing regions inside the remnant, pointing at gyroradii of protons with energies of 10 14 − 10 15 eV (Eriksen et al. 2011). The anisotropy in the particle distribution function can be viewed hydrodynamically as a CR pressure gradient, and by scattering off the waves, which are magnetic field fluctuations, frozen into the plasma, energy and momentum is transferred from the CRs to the plasma via the waves as a mediator. Since the CRs eventually leave the Galaxy, their pressure gradient assists the thermal gradient to drive an outflow or Galactic wind. The associated CR adiabatic index is γ c = 4/3, as expected for ultra-relativistic particles. Therefore the CR pressure gradient pointing away from the Galaxy is less steep than the thermal pressure gradient for a γ g = 5/3 monatomic plasma, enabling CRs to drive a wind to larger heights with a lower mass loss rate but higher terminal wind speed (Breitschwerdt et al. 1991). In addition, CR protons, which are the main drivers, do not suffer radiative losses like the plasma, and can therefore accelerate the Fig. 11 Probability density function for the electron density distribution in the midplane of the Galaxy, averaged over 45 lines of sight (LOS, length 1 kpc), taken at times between 420 and 500 Myr. The highest spatial resolution of the simulation has been 0.5 pc outflow as long as pressure gradients exist. Although a mass loss rate of about 0.3 − 0.4 M /yr for the Galaxy is not dramatic, it can nevertheless have an important influence over ISM simulation times of several hundred Myr. Although the theory of CR driven galactic winds took quite some time to catch on, as we know from own experience, there meanwhile

The local bubble as a test laboratory
The solar system is embedded in a region of comparatively low density, n ∼ 10 −3 cm −3 , emitting copiously in soft Xrays (0.1 − 0.7 keV) and being responsible for a significant amount of the soft foreground emission in the ROSAT All Sky Survey (RASS). This region, called the Local Bubble, has an extension of ∼ 200 pc in the disk, and ∼ 600 pc perpendicular to it. It is located in the local ISM (LISM) within a larger star forming region, which was historically called the Gould's Belt for its ringlike appearance, but has been recently found from accurate distance measurements to gas clouds by the Gaia satellite, to form a large coherent wavelike structure (Alves et al. 2020). The LISM surrounding the solar system has been rather active in the past. For example, radio synchrotron emission from relativistic electrons, which are accelerated at SN/SB shock waves have been observed long ago (Berkhuijsen et al. 1971) and are indicative of extinct SNRs/SBs, e.g. the so-called radio Loops 2 and 3 as well as the still active Loop I (LI) SB. Note that recently there have been claims to associate LI with cavities much further away (Puspitarini et al. 2014), even with the Fermi Bubbles in the Galactic Centre (Kataoka et al. 2013). But in contrast to the LB there is still an active stellar cluster in the centre of LI, i.e. the Scorpio-Centauraus (Sco-Cen) association, and it has also been argued that the younger Upper Scorpius group is the product of triggered star formation in the Sco Cen association (Preibisch and Zinnecker 2007).
Some time ago it was claimed (Koutroumpa et al. 2009) that a substantial fraction of the local soft X-ray emission is due to charge exchange reactions on frozen-in solar wind ions at high ionisation stages (SWCX), in other words, the X-ray emission would be very local, i.e. from the solar system(!). Indeed there is some contribution to the emission, but it could not explain soft X-ray emission at higher latitudes. Later the contribution of solar wind charge exchange reactions to the local 1/4 keV X-ray background was estimated to be about 40% (Galeazzi et al. 2014). Based on sounding rocket experiments, the DXL mission also supports the existence of the LB after subtracting the SWCX (Liu et al. 2017), although assuming the LB to be in CIE, which, as we shall argue, is probably not the case. Therefore there remains a considerable amount of emission that is very likely of SN origin. Another strong argument for the existence of the LB comes from the observation of a Local Cavity, which exhibits a significant depression of HI emission. This has been observed in a number of ways, first as absorption line measurements towards individual stars with the EUVE satellite, giving low HI column densities, N(HI) ∼ 10 17 − 10 18 cm −2 up to about 100 pc, followed by a noticeable jump to N(HI) ∼ 10 20 − 10 21 cm −2 (which probably marks the LB shell) for larger distances (Fruscione et al. 1994). Later Sfeir et al. (1999) used NaI as a tracer which is much more sensitive to lower column densities, and more recently, differential extinction maps out to several kpc have been calculated using new Gaia data (Puspitarini et al. 2014;Capitanio et al. 2017;Lallement et al. 2018). Here and in the following we adopt the view that the LB and LI are local SBs in a still active starforming region of the ISM.

SN origin of the local bubble
The missing stellar content in the LB has put doubts on the SB origin of the LB, and it was suggested that instead it was an interarm region with gas flowing over from the LI SB (Frisch 1995); for a debate at that time, see Breitschwerdt et al. (1996). We have therefore explored the possibility whether the LB hosted a stellar moving group in the past (Berghöfer and Breitschwerdt 2002) and found that the subgroup B1 from the Pleiades crossed the LB region at about 15 Myr ago, while its still surviving members are now part of the Sco-Cen (see also Maíz-Apellániz 2001;Benítez et al. 2002). Indeed it has been shown (Fernández et al. 2008) that quite a number of stellar groups have moved through the LISM in the past. From their list, our group found e.g. a high probability for a SN about 14 Myr ago in the moving group β Pic (Göpel 2018, Bachelor Thesis, TU Berlin).
In the following, we briefly describe the method we have developed in order to trace the motion of stellar groups and the time and location of possible explosions therein. In Fuchs et al. (2006) we performed a thorough search in a sphere centred at the solar system, with a diameter of 400 pc, which marks the limit of distance inaccuracies for parallaxes of > 5 milli-arcseconds in the Hipparcos data. These were supplemented by the Arivel catalogue for radial velocities, such that all locations and velocity vectors in 3D space are known. Looking for young main sequence stars we selected those with B − V < −0.05, discarding dwarfs, subdwarfs and otherwise peculiar stars, resulting in an unbiased sample of 762 stars in total, out of which 610 have radial velocities in the Arivel data base. The others have B − V > −0.1 and are therefore not relevant for our analysis. Next, cluster membership was tested by selecting those stars which are close to each other in physical and velocity space in order to exclude projection effects, resulting in 236 stars, which, after careful de-reddening, were put into an HR-diagram and compared with isochrones from stellar evolution calculations from Schaller et al. (1992) to determine the cluster age. This was derived to be between 20 − 30 Myr, which is on the high side, if compared to other authors, and just about overlapping with Sartori et al. (2003), who determined it to be 16 − 20 Myr. We have therefore used as cluster age, t cl ∼ 20 Myr. Comparing the final sample with the extensive compilation by de Zeeuw et al. (1999), we found that in the subgroups of Sco-Cen there were 79 stars in common, belonging to Upper Scorpius (US), Upper Centaurus Lupus (UCL) and Lower Centaurus Crux (LCC).
Having all space and velocity components available, one can calculate the trajectory of the stellar group backwards in time, and study when and where SNe exploded. Since 20 Myr is less than 10% of a Galactic revolution at the solar circle, it is still possible to use the epicycle approximation for the kinematics of the stars. The next step is to find the number and masses of already exploded stars. Since there is no direct information about the SN progenitors other than the existence of SNRs and SBs, respectively, and indirect evidence from cosmic rays and some radioactive isotopes, 4 we need a statistical approach. The distribution of the number of stars over their masses when star formation takes place has been investigated for more than 60 years and has been approximated by a so-called initial mass function (IMF, cf. Salpeter 1955;Miller and Scalo 1979;Kroupa 2001), which has originally been derived from the luminosity function, and which seems to be rather common for the Galaxy (Krumholz 2014) and may even be universal (Bastian et al. 2010) over cosmic time after a few generations of star formation. Also the parent molecular clouds follow a similar distribution, which raises the question if eventually structure forming processes like turbulence in the ISM may be responsible for the IMF. The number of stars dN in a mass interval dm is given by where N 0 is a normalisation constant, and α = + 1 = 2.3 (Kroupa 2001), for stars with m ≥ 0.5 M , which compares well with the value of Salpeter (1955) of 2.35, who used a single power law. But for the lower mass end the distribution switches to lognormal 5 with a corresponding steepening of the IMF. Since we are only interested in the high mass end of the IMF either value can be used. As an aside, the break in the spectrum (from power law to lognormal distribution) may point to some difference in the star formation process between high and low mass stars. For practical purposes the total mass M(m 1 , m 2 ) in a given interval, limited by m 1 and m 2 , respectively, is important, and can be derived from the first moment of Eq. (3), i.e.
which at the lowest level allows us to calculate the mass in an individual bin in order to derive the mass of the SN progenitor. Using the 42 and the 27 members of UCL and LCC, respectively, one can calibrate the IMF (see Fuchs et al. 2006). Accordingly, 16 SNe with progenitor masses between 9 and 20 M have been derived, the first having exploded presumably about 13 Myr ago (Breitschwerdt et al. 2016), thereby marking the most likely explosion centre of the LB. There have been arguments that the LB should be due to a single explosion and also that the LB has been reheated just by one SN. This may explain the measured 60 Fe, generated in late type star nucleosynthesis, peak at about 2.2 Myr before present (BP) from the ferromanganese ocean crusts (found in many different locations in the sea and on the moon, which is therefore a global sign). However, the broadness of the peak, as well as an earlier peak 7 -9 Myr BP found by Wallner et al. (2016) hint at several explosions. Statistically, the most probable distribution of masses in the IMF gives each mass bin the same statistical weight, meaning that we use a variable bin size, such that each bin contains exactly one star (see Maíz Apellániz and Úbeda 2005). This naturally results in a star with maximum mass M max ∼ 19.86 M in our simulations, which, in co-evally born stars in a cluster, should have been the first in a series to create the LB. There is a caveat however. The procedure we have described here is the best, if we do not have any specific information depending on the time sequence and locations of explosions. Since the LB, as any other bubble, is always a singular realisation of a statistical event, it may well deviate individually from the statistically most probable model, and according to the Bayesian concept of probability, we should use any information available in the future in order to improve on our prior. In fact, the existence of a second peak in the 60 Fe fluence may be an indication to deviate from that, as preliminary calculations, we have performed, seem to sustain.
Having determined the cluster age from comparison with stellar isochrones, and the main sequence lifetime of the SN progenitors according to the IMF binning, the explosion time is simply the difference between the two, i.e. t exp = t MS − t cl . The explosion sites can be found by calculating backwards the stellar trajectories from the present to the past, using the epicyclic approximation (see Fuchs et al. 2006). A major result, as has been shown by Breitschwerdt et al. (2016), is the explanation of the measured 60 Fe peak fluence at 2.2 Myr BP, which corresponds naturally to the time of closest approach of the cluster to the solar system. To model the 60 Fe transport, it has been assumed that 60 Fe was incorporated into dust grains, which were able to penetrate deeply into the solar wind and subsequently sediment on the Earth's surface and the ocean floor. Its quantity is determined by so-called uptake factors (responsible for atmospheric and ocean sedimentation) and stellar mass depen-dent 60 Fe yields, which are taken from stellar evolution calculations (e.g. Woosley and Weaver 1995). In summary, our SN model of the LB can very well reproduce the measured 60 Fe fluence and the peak 2.2 Myr BP, as measured by Knie et al. (2004) and, more recently, by Wallner et al. (2016).

NEI simulations of the local bubble
Having determined the stellar content in the LB and the space and time sequence of explosions, one can perform hydrodynamic simulations self-consistently with the timedependent ionisation structure as discussed in Sect. 4. In de Avillez and Breitschwerdt (2012b), we compared it to the measurements of highly ionised Li-like species, such as CIV, NV and OVI, measured by FUSE, IUE, GHRS-HST and Copernicus.
To set up the simulation correctly, we remember that the cluster, which is responsible for the SNe in the LB, still has SN progenitors below 8.8 M which now belong to the Sco Cen. This cluster, on its own, has according to a similar IMF analysis, stars exploding at about the same time as in the LB, thus generating the Loop I SB (LI). From the analysis of the ROSAT soft X-ray emission and its anti-correlation with HI data, Egger and Aschenbach (1995) have shown the existence of an HI ring in conjunction with a depression of the soft X-ray flux in the direction towards the Galactic Centre, where LI is located, with the North Polar Spur being presumably part of its shell expanding out of the disk (but see the discussion by Puspitarini et al. (2014)). The ring itself is interpreted as the collision shell between the LB and the LI, giving rise to a dense shock compressed layer. Deriving the stellar content of LI with the same method as for the LB, numerical simulations  have shown that the collision of the bubbles happened about 3 Myr ago, and that the ring will start fragmenting in about 3 − 5 Myr from now. The reason is that the interval between SN explosions in the LI is somewhat smaller, and therefore the pressure in LI is higher. This leads to an acceleration of the interaction shell towards the LB, pushed by a hot tenuous gas, which therefore is subject to Rayleigh-Taylor instabilities (RTIs). It has been calculated by Breitschwerdt et al. (2000) that the RTI induced fragmentation can be the source of the local clouds (i.e. small fluffy cloudlets, called LIC, G-cloud etc.) through which the solar system is moving. The joint dynamically and thermally self-consistent evolution of the LB and LI, allows to fix the age range of the LB by comparing specific column density line ratios of CIV/OVI and NV/OVI to the aforementioned observations. It turns out that the LB should be in an evolutionary state between 0.5 and 0.8 Myr after the last SN reheated the LB so that N (OVI) < 8 × 10 12 cm −2 , and log N (CIV)/N (OVI)]< −0.9 and log[N (N V)/N (O VI)]< −1 in the LB interior. We note that according to the simulations explaining the 60 Fe fluence, the last SN inside the LB occurred about 1.5 Myr ago, which seems to be inconsistent with the NEI results. It should be noted, however, that the precise volume of the LB and its true ambient environment leave some uncertainties in the LB temperature and pressure. The essence here is rather that NEI, given some input data, constrains the evolution time rather well in contrast to CIE. Moreover, the errors due to poor knowledge of uptake factors, 60 Fe yields of SN progenitors, measurement errors of the Hipparcos data, also leave some scope and flexibility in the interpretation of the explosion time of the last SN. Hopefully all these input data will become more precise in the future, e.g. with the new Gaia data and improving stellar evolution calculations.

Summary and conclusions
ISM research has gone a long way since the field was established by the discovery of Hartmann (1904) from the Astrophysical Observatory Potsdam, who found so-called stationary CaII lines in the spectrum of a binary star, which turned out to be of non-circumstellar origin. Until the late 60's, before the launch of the UV satellite Copernicus in 1972, the ISM was thought to consist of a cold and weakly ionised warm gas (so-called 2-phase model of the ISM; Field et al. (1969)), in which energy input by low energy cosmic rays was balanced by line cooling driving the plasma by thermal instability (Field 1965) into two stable phases. With the discovery of the diffuse soft X-ray background (Bowyer et al. 1968), and with the detection of a widespread OVI absorption line (Jenkins and Meloy 1974), which was subsequently interpreted as the result of a SN dominated tunnel network of a hot phase (Cox and Smith 1974), the importance of a SN regulated ISM, as was described in the 3-phase model (McKee and Ostriker 1977), became widely acknowledged.
However, both the factual and conceptual differences of the "modern" ISM studies by numerical simulations could hardly be larger. The work of many researchers, whose decisive contributions could only be mentioned in passing here (or not even mentioned at all for lack of space), has led to a change in paradigm. The idea of an ISM distributed into distinct phases was helpful for the interpretation of early observations, but is obsolete now. Instead the new appropriate concept, which we have outlined in the previous paragraphs is the following: • The ISM is highly filamentary and its structure is determined by turbulence. • Since turbulence is dissipative, there has to be a constant energy input, which is mainly due to SNe, and which is of the order of 10 41 − 10 42 erg/s in the Galaxy. As a consequence, the Reynolds number is very large, i.e. Re ∼ 10 6 − 10 7 , many orders of magnitude higher than achievable in any laboratory experiment.
• The integral scale at which turbulence is fed in is about 75 pc according to our simulations, corresponding to the average diameter of superbubbles breaking up in a very inhomogeneous ISM (de Avillez and Breitschwerdt 2007). • Thermal instability can be suppressed by turbulence on scales of 1 pc and below (de Avillez and Breitschwerdt 2005), and therefore the amount of gas in thermally unstable phases can be large (more than 50% for HI). • The flow is ram pressure dominated, except for the dense clouds, in which magnetic pressure is higher, and for the hot tenuous gas above 10 6 K, which is controlled by thermal pressure. • The volume filling factors are very different from the MO77 model, especially the hot "phase" does not fill more than 20% of the disk ISM. • The probability density functions (pdfs), which are a measure for the amount of gas to be found in given parameter regions, e.g. for density, temperature or pressure distributions, are rather broad and lognormal as a result of turbulence, with a sizeable fraction of the gas occupying thermally unstable regions. This is, for example, in very good agreement with the broad electron distribution found in HI gas as derived from pulsar dispersion measures (Berkhuijsen and Fletcher 2008). • Hot gas cannot be kept down in the disk by disk parallel magnetic fields; instead the field lines are strongly tangled, both in the disk and the halo. • A dynamical equilibrium can only be approximately established if the SN energy input rate is constant, and only after the Galactic fountain duty-cycle has been set up after 150 Myr in case of the Milky Way. • Interstellar clouds appear to be shock compressed layers, and thus transient features with lifetimes of the order of 10 6 − 10 7 yr (see Ballesteros-Paredes 2004). • Due to different intrinsic timescales, the thermal and the ionisation structure of the ISM plasma are coupled and lead to emission and absorption line spectra of gas which is not in ionisation equilibrium (NEI). • Spectral lines with sufficient resolution contain more information than usually applied fitting procedures based on the assumption of collisional ionisation equilibrium (CIE) can extract. • The cooling function is time-dependent and determined by the plasma's thermal and ionisation history. • NEI evolution can change spectra dramatically, leading, e.g., to X-ray emission at low electron temperatures of 10 4 K, if the gas was previously hot enough. • The NEI electron distribution in the ISM also depends on the dynamical and thermal history of the plasma, and is in good agreement with observations of pulsar dispersion measurements. • It can be shown that in the Local Bubble, our local ISM test laboratory, line ratios in NEI simulations can be different from the usually used CIE simulations (or spectral fit programs), and can in principle explain observations, as well as constrain the age of the Local Bubble. • Cosmic Rays are an essential component of the ISM, which interact strongly with magnetic fields and magnetohydrodynamic waves, and thereby transfer momentum to the plasma, thus being capable of driving a galactic wind. Their inclusion in numerical simulations is difficult, because transport coefficients like the diffusion coefficient are energy dependent, vary in space and time, and are to a large extent determined by microscopic processes. Future simulations will show how to treat them on an appropriate level, and what their consequences for the ISM will be.
Thermal and ionisation histories of plasmas can be complex and therefore numerical simulations including the full NEI structure are more demanding. In fact, the attractiveness of using CIE lies in its simplicity. Nevertheless, we believe that numerical ISM simulations are now in a stage that timedependent cooling should be implemented routinely in order to learn more about the ISM evolution.