Impact of Supernovae on the Interstellar Medium and the Heliosphere

It is well known that interstellar clouds have line widths which are consistent with supersonic motions. Since turbulence decays very fast, in fact in less than a free fall time (e.g., Mac Low et al.  1998; Stone et al.  1998), irrespective whether it is hydro (HD)- or magnetohydrodynamical (MHD), energy has to be constantly supplied. This energy is drained from the average fluid motion via Reynolds stresses, giving rise to a stress tensor, whose diagonal elements can be interpreted as turbulent pressure , whereas the off-diagonal elements may be considered as turbulent viscosity , which can formally be included into the averaged momentum equation by adding the viscous stress tensor (last term in Eq.  1) to the normal and shear stresses \(\mu \frac{\partial \langle u_{i}\rangle } {\partial x_{j}}\) for an incompressible Newtonian fluid, where μ is the dynamic viscosity and \(\langle \ldots \rangle\) denotes the statistical ensemble average. In principle the velocity field could be calculated from the continuity, momentum, and energy equations, were it not for the existence of the Reynolds stresses . The straightforward procedure would be to calculate them by taking higher-order moments of the equations, which, as is well known, runs into a closure problem. Although the Navier-Stokes equations are deterministic, in a statistical approach the system is underdetermined and cannot be solved without additional ad hoc information. Such a model is, for example, given by the concept of eddy viscosity, where the term eddy loosely describes the swirling motion of the velocity field \(\mathbf{u}(\mathbf{x})\), whose size is given by the correlation length \(r = \vert \mathbf{r}\vert\) (for isotropic turbulence), derived from the autocorrelation function \(A(\mathbf{r}) = 1/V \int \mathbf{u}(\mathbf{x})\mathbf{u}(\mathbf{x} +\mathbf{ r})dV \equiv \langle \mathbf{ u}(\mathbf{x})\mathbf{u}(\mathbf{x} +\mathbf{ r})\rangle\), averaged over some volume V. As is shown in Eq. ( 1), viscous stresses “look” similar to shear stresses, although the former is a property of the turbulence and the latter of the fluid.

Stretching of fluid elements due to internal stresses increases their vorticity , \(\boldsymbol{\omega }= \nabla \times \mathbf{ u}\), since due to the conservation of angular momentum, the rotation must increase, if the moment of inertia is decreased by thinning out the vortex blob into a tube. This is the essence of turbulent mixing in the ISM, which is usually so much faster than molecular diffusion. In fact, the turbulent viscosity is related to the kinematic viscosity ν by \(\nu _{t} =\mathrm{ Re}\nu = uL\), where the Reynolds number \(\mathrm{Re} = uL/\nu\) measures the ratio of inertial to viscous forces and is generally very large in the ISM, i.e., Re ∼ 10 ⁶ − 10 ⁷ (Elmegreen and Scalo  2004).

Due to the strain field, which describes the distortion of a fluid element, its kinetic energy and vorticity increase as it becomes thinner. Since energy has to be conserved during this process, this must occur at the expense of the kinetic energy on the larger scale, thus corresponding to a transfer of energy from large to small scales. If such a process is self-similar, energy cascades through an inertial range (Re ≫ 1) from large scales, and it is deposited on the smallest scale, where it is dissipated (Re ≤ 1). The latter must occur, because Re becomes progressively smaller as the scale decreases, until a scale η with turbulent velocity v is reached, where \(\mathrm{Re} = v\eta /\nu \sim 1\), and viscosity takes over. Thus essentially two parameters enter the problem: the viscosity ν and the constant kinetic energy dissipation rate \(\epsilon \sim u^{3}/L \Rightarrow u \sim L^{1/3} \propto k^{-1/3}\). Dimensional analysis shows that the spectral energy density \(E(k) \propto k^{-5/3}\), the famous Kolmogorov ( 1941) law , since \(\int E(k)dk \propto \langle (\varDelta u)^{2}\rangle \sim k^{-2/3}\), where 〈( Δ u( r)) ²〉 is the second order structure function , representing the mean kinetic energy in a turbulent eddy of size r. A direct conclusion, which can be drawn from these scaling laws, is the coupling of scales over a wide range. Since \(\epsilon \sim u^{3}/L \sim \nu v^{2}/\eta ^{2} \Rightarrow \eta \sim l\mathrm{Re}^{-3/4}\), resulting in \(l/\eta \sim 10^{4} - 10^{5}\) for interstellar conditions. Therefore turbulence in the ISM is more than just adding an extra pressure term to the Navier-Stokes equations. The vorticity field is self-advecting. It is generated by a local velocity field, which in turn depends via Biot-Savart’s law on the global distribution of vorticity due to \(\mathbf{u}(\mathbf{x}) = \frac{1} {4\pi }\int \frac{\boldsymbol{\omega }({\boldsymbol x^{{\prime}}})\times \mathbf{r}} {r^{3}} d{\boldsymbol x^{{\prime}}}\), \(\mathbf{r} =\mathbf{ x} -{\boldsymbol x^{{\prime}}}\); here \(r = \vert \mathbf{r}\vert\) and \(\mathbf{x}^{{\prime}}\) is the point in space, at which the vorticity distribution is given and which is integrated over all these points to obtain the velocity field at point \(\mathbf{x}\). It has been shown that in view of the large decay rate, the most efficient driver of turbulence in the ISM is SNe (Mac Low and Klessen  2004). Eddies on the so-called integral scale L evolve by purely inertial forces and break up due to dynamical instabilities into smaller-size eddies ( L ≫  l ≫  η) on a time scale \(\tau \sim l/\varDelta u \sim l/u\), where Δ u is the velocity fluctuation within an eddy of size L and τ is the turnover time. Using the Kolmogorov scaling law for incompressible turbulence, τ ∝  l ^2∕3, implying that within a turnover time on the largest scale, dynamical instability has generated a swarm of smaller eddies within, like in a matryoshka doll, thus promoting an energy cascade down to the viscous scale.

However, interstellar turbulence is in general not incompressible. Instead, high Mach number flows are ubiquitous, due to large pressure gradients generated by SNe, stellar winds, etc., as mentioned before, acting like a piston on the flow. Therefore in a sufficiently large region in the ISM, there will always be an ensemble of shock waves . Thus an energy cascade, in which the energy will be simply transferred from one eddy to the next, will not exist, because dissipation will occur on any scale of the size of a shock width, which is of the order of a collisional mean free path, λ _mfp. But owing to the low density in the ISM, λ _mfp is exceedingly large. Therefore shocks, which are observed to exist in the ISM, are collisionless , and the randomization of ordered gas motion occurs by turbulent electromagnetic fields, induced by, e.g., the Weibel instability (Weibel  1959) , due to an anisotropic plasma distribution in momentum space.

In general the velocity field in some region will look like a series of step functions, and taking the Fourier transform of it will result in a power spectrum E( k) ∝  k ⁻², as opposed to \(k^{-5/3}\) in the incompressible case. An ensemble of shocks in homogeneous isotropic turbulence will then have a similar slope. This may not seem significant by numbers, but bears a fundamental difference: whereas in Kolmogorov turbulence energy decays through a cascade, in supersonic turbulence all scales are coupled in a dissipative shock due to its k ⁻²-spectrum, i.e., its nonlocality in k-space. It has been mentioned already by von Weizsäcker ( 1951), and shown later by Fleck ( 1996), that in a density hierarchy of compressible fluid elements of two successive levels, i.e., \(\rho _{j}/\rho _{j-1} = \left (l_{j}/l_{j-1}\right )^{-3\alpha }\) where j is the level of hierarchy and α the compressibility of the fluid (0 <  α < 1), a Kolmogorov-type scaling law can be recovered, if we use the energy density per unit volume instead of per unit mass, i.e., \(\epsilon _{V } =\rho \epsilon =\rho u^{3}/l \Rightarrow w =\rho ^{1/3}u\), where w is the mass-weighted velocity. Consequently, \(w^{3}/l \sim const. \Rightarrow u \sim l^{1/3}/l^{-\alpha }\sim l^{\alpha +1/3}\), so that eventually \(E(k) \sim k^{-1}u^{2} \sim k^{-5/3-2\alpha }\), and the correction due to compressibility is k ^−2α . The compressibility parameter α is related to the fractal dimension by \(D = 3 - 3\alpha\), such that 3 α is the number of dimensions in which the element is compressed. Numerical simulations of isothermal compressible turbulence (Kritsuk et al.  2007) have shown that E( k) ∼  k ^−1. 97, which is in agreement with our previous estimate E( k) ∝  k ⁻², and the simple scaling model of Fleck ( 1996), if α ≈ 0. 15, and hence D = 2. 55. In the limit of incompressible flow ( α → 0), D → 3, corresponding to a space filling turbulence without intermediate dissipative structures. Discrepancies of these scaling models with respect to experiments and also numerical simulations are due to their neglect of intermittency , i.e., the sudden disappearance and outburst of turbulence, which shows up particularly in higher-order structure functions (e.g., She and Leveque  1994).

With these basic ideas at hand, we should now be able to interpret the results of SN-driven turbulence in ISM models.

Observationally, the ISM in galaxies appears to be frothy and filamentary, exhibiting structures on all scales, like, e.g., in a recent combined infrared image of Herschel and Spitzer data (Galametz et al.  2016). This is most likely a result of turbulence , acting on and coupling a huge range of scales. To model this behavior, 2D and 3D simulations, either global or of a box, cutting out a representative patch of the ISM, have been performed by several groups (e.g., Rosen and Bregman  1995; Korpi et al.  1999; Wada and Norman  2001; de Avillez and Breitschwerdt  2004,  2005; Dobbs et al.  2011; Girichidis et al.  2016). Common to these models is that turbulence is SN driven, which, as has been argued before, is dominant in regions where the SFR is sufficiently large, like in the star-forming disk of the Milky Way. In the outer Galaxy, where the SFR drops and the disk begins to flare, turbulence may still be generated by the magnetorotational instability (Piontek and Ostriker  2005). In the following, we will describe such a model in more detail in order to see how the ISM will react on stellar feedback processes. In order to track discontinuities in the flow (shocks and contact surfaces), these simulations are performed with a parallelized hydrocode (e.g., de Avillez and Breitschwerdt  2004; de Avillez et al.  2012), based on a conservative finite-volume scheme with adaptive mesh refinement (AMR).

The numerical setup of a mesoscale ISM simulation is as follows: consider a typical patch of ISM centered on the solar system, with a computational box which extends 1 kpc within and ± 15 kpc perpendicular ( z-direction) to the star-forming disk on either side and with a highest AMR resolution of 0.25 pc. The large z-extension has been chosen in order to capture the Galactic fountain flow, which arises from an overpressured gas disk due to SN explosions. As an initial model, a vertical density distribution is taken, representing molecular, neutral, ionized, and hot gas phases (e.g., Ferrière  2001). Initially the gas, which is immersed into the galactic gravitational field, is in hydrostatic equilibrium, but collapses into a thin disk, until sufficient pressure by star formation processes is built up to puff up the disk. Star formation is implemented by using an initial mass function (IMF), within a mass range of 8–60 M _⊙, respectively, which controls the distribution of stars as a function of mass. Whenever in a certain region in the box the density exceeds 10 cm ⁻³ and the temperature falls below 100 K, corresponding to a Jeans mass of \(M_{J} \simeq \frac{\pi ^{5/2}} {6} (\frac{k_{B}T} {\bar{m}} )^{3/2}1/(G^{3/2}(n\bar{m}))^{1/2}) \approx 10^{5}\,\mathrm{M}_{ \odot }\,(T/100\,\mathrm{K})^{3/2}\,(n/10\,\mathrm{cm}^{-3})^{-1/2}\), gravitational collapse is initiated, and after fragmentation stars begin to form. This process is restricted such that on average an SFR corresponding to the Kennicutt-Schmidt law is obtained and that the stellar masses (both high and low) are distributed according to a Salpeter IMF (Salpeter  1955). The stars evolve according to their main-sequence lifetime (Fuchs et al.  2006), and spectral types O and B have a preceding stellar wind phase during which a kinetic energy of about 10 ⁵⁰ erg is injected into the ISM. All stars receive a random drift velocity of ∼ 5 km∕s (Blaauw  1964) when they are born. After their main-sequence evolution, in accordance with observations, 60 % of the stars explode in a cluster and the rest uncorrelated in the field as SNe types Ib, Ic, and II (Cappellaro et al.  1999) and as type Ia randomly in an exponential disk (Freeman  1987), each of them releasing 10 ⁵¹ erg. A diffuse interstellar radiation field, which provides some background heating, has been implemented (Wolfire et al  1995). The boundary conditions are periodic on the sides of the box perpendicular to the disk and outflow parallel to it (top and bottom sides). The simulations show that about 15 % of the initial mass in the box is lost after 400 Myr of simulation time or equivalently 25 % of the galactic SFR (Breitschwerdt et al.  2012) compared to 17 %–41 %, derived from Spitzer/IRAC GLIMPSE data (Robitaille and Whitney  2010). The simulations have to be run at least for 150 Myr, the average time it takes to close the galactic fountain cycle, during which material rises into the halo and partly escapes as a wind, but most of it returning back to the disk due to radiative energy losses and insufficient pressure support against gravity. After this time the gaseous disk relaxes into a dynamical equilibrium state, if the SFR is kept constant.

An important process for the evolution of the interstellar plasma, which has hitherto been largely ignored, is its nonequilibrium ionization (NEI) structure. For convenience, it is often assumed that the plasma is in collisional ionization equilibrium (CIE), in which the ionization rate locally and instantaneously balances the recombination rate. However, this is at best a useful approximation for interstellar plasmas with T ≥ 10 ⁶ K. Below this temperature, where most of the radiative line cooling occurs, this cannot be true (for an early discussion of the problem see Kafatos  1973; Shapiro and Moore  1976), since the time scales for collisional ionization and radiative recombination can be vastly different. It has been shown that fast expansion of a hot plasma initially in CIE can lead to strong adiabatic cooling, lowering drastically the kinetic temperature of the electrons, while the recombination of highly ionized species lags behind (delayed recombination; see Breitschwerdt and Schmutzler  1994). The emerging spectra can contain many X-ray recombination lines, although the electron temperature might be as low as 4 × 10 ⁴ K (Breitschwerdt and Schmutzler  1994; de Avillez and Breitschwerdt  2012b). More dramatically, the interstellar cooling function, which accounts for all radiative energy losses, can deviate from the standard CIE case by more than two orders of magnitude and will in general be a function of space and time. The results will be presented in the next section.

The filamentary and frothy structure of the ISM, displayed here at t = 290 Myr, well after the galactic fountain cycle has been established, can be clearly recognized in density (see Fig.  1) and temperature (s. Fig.  2) images, representing cuts through the galactic midplane. The low-density regions are interstellar bubbles (or relics thereof) and carved out by SNe. In fact, the solar system resides inside a so-called Local Bubble (see, e.g., Breitschwerdt and de Avillez  2006; de Avillez and Breitschwerdt  2012a; Galeazzi et al.  2014), which is the closest part of ISM to the heliosphere . The highest densities (red color) are mainly associated with shock-compressed layers, resulting from supernova remnant (SNR) or superbubble shock waves. These transient filaments or shells exist only a few million years before they are destroyed by strong shear flows. They are the sites of new star formation, which again generates new bubbles of low density and recombining NEI gas. It can be seen that the average size of these bubbles is between 75 and 100 pc, corresponding to the integral scale, at which turbulence is fed in de Avillez and Breitschwerdt ( 2007). The latter can be derived from the scale of flattening of the second-order velocity structure function \(\langle \left [\varDelta u\right ]^{2}\rangle =\langle \left [\mathbf{u}(\mathbf{x} + r\,\hat{\mathbf{e}}_{x}) -\mathbf{ u}(\mathbf{x})\right ]^{2}\rangle\), which is a measure for the kinetic energy per unit mass contained within eddies of size r (or less).

In classical textbook ISM models (e.g., McKee and Ostriker  1977), the gas is distributed into three different “phases” (cold, warm, and hot, with molecular clouds being decoupled, because they are gravitationally bound), i.e., within certain regions in an n- T-diagram, which are stable with respect to entropy perturbations at constant pressure and coexist in pressure equilibrium. In the unstable regions in between, no or only very little transiting gas should be found. This can be tested by looking at the density and pressure probability density functions (pdfs). One would expect that for overall pressure equilibrium, the pressure pdf would be very narrowly peaked around \(P/k_{B} \sim 10^{4}\,\mathrm{cm}^{-3}\,\mathrm{K}\), typically. However, as can be seen from the pressure pdf (s. Fig.  3), the distribution of ISM pressures is fairly broad, and, likewise, the density pdf spans more than five orders of magnitude from 10 ⁻⁴ to 10 ² cm ⁻³ (s. Fig.  4), implying that the exchange of material between different “phases” at constant pressure is a myth. Instead, turbulence couples a large range of scales, mixing plasma of different density and temperature very efficiently by turbulent diffusion, which supersedes molecular diffusion and heat conduction by orders of magnitude (see Sect.  2.1). Hence a lot of material is found in thermally unstable regions, in fact, almost 50 % of the H i gas, in agreement with observations (Heiles and Troland  2003).

A consequence of SN-heated material punching holes through the disk and rising into the halo is a lower overall pressure in the disk, \(P/k_{B} \sim 3000\,\mathrm{cm}^{-3}\,\mathrm{K}\), as can be seen from the pressure pdf.

One might expect that an initially disk-parallel magnetic field could prevent breakout due to magnetic tension forces. These effects have been studied in ISM MHD simulations, similar to those discussed above (e.g., de Avillez and Breitschwerdt  2005; Korpi et al.  1999). In essence a magnetic field with a total strength (regular plus turbulent component) of 5 μG cannot inhibit the fountain flow significantly. Instead, the upward flow just shows more coherent and bubble-like structures, as compared to the unmagnetized case. Also the volume filling factors of the hot ( T > 10 ^5. 5 K) gas is for galactic SN rates about 20 %, both for the HD (de Avillez and Breitschwerdt  2004) and MHD (de Avillez and Breitschwerdt  2005) cases due to disk outflow, as compared to 70 %–80 % in the standard model (McKee and Ostriker  1977).

Apart from the magnetic field, whose energy density is comparable to the thermal and kinetic energies of the plasma, a similar amount goes into CRs, which can transfer some of their momentum and energy to the gas by scattering off magnetic irregularities (in general MHD waves), which are frozen into the plasma.

CRs are the high-energy component of the Universe and below 10 ¹⁸ eV of the ISM. The term is a misnomer, as the primary component consists of high-energy charged particles, with photons being only a by-product of their interaction with matter. CRs are predominantly protons, electrons (about 1 %), and heavier nuclei, close to ISM abundances, with notable differences in some elements, which are mainly the spallation products of heavier nuclei (e.g., iron) due to collisions with ISM gas. Observationally, CRs exhibit a power-law differential energy spectrum dN( E) ∝  E ^−γ   dE in energy from about 10 ⁹–10 ²¹ eV with two spectral breaks at 10 ¹⁵ (“knee”, γ ∼ 2. 7) and 10 ¹⁸ (“ankle”, γ ∼ 3. 0) eV. Particles with higher energies are of extragalactic origin, because their gyroradius, and hence any accelerator, exceeds the size of the Galaxy. This criterion, put forward by Hillas ( 1984), arises from the fact that the maximum energy a particle can attain can be derived from integrating over the Lorentz force, E _max ≃  ∫ ₀ ^L ( q∕ c)  v  B  dl, where q, c, L, B, and v are the electric charge, the speed of light, the size of the acceleration region, the magnetic field strength there, and the effective particle speed, respectively. Of course, it is not the magnetic fields which accelerate the particles, but the electric fields associated with it due to induction. CRs interact resonantly with MHD waves, arising from magnetic field perturbations, due to a gyro-resonance condition for wave numbers of the order of the gyroradius \(r_{L} = p/(qB)\), where p is the (relativistic) particle momentum. Loosely speaking, the particles pick out the waves from the Fourier spectrum, which fulfill the resonance condition, thus getting strongly scattered in pitch angle, essentially performing a random walk, i.e., diffusion, in the plasma. Therefore their effective drift velocity is of the order of the Alfvén speed , \(v_{A} = \sqrt{B/(4\pi \rho )}\), in the ISM typically less than 0.1 % of the speed of light. This also explains the near isotropy of CRs in arrival directions observed with satellite and terrestrial detectors.

In the Galaxy, the most likely accelerators are shock waves in the ISM, in particular those of SNRs. The process is called diffusive shock acceleration (DSA) or first-order Fermi acceleration , in which up- and downstream of the shock, MHD waves act as magnetic mirrors, reflecting the particles across the shock front like a ping-pong ball across the net. In the shock frame, the particles see the waves, which are frozen into the converging plasma, approaching the shock front on either side, thus converting plasma kinetic energy into particle energy. In the analogy, this would correspond to the ping-pong ball gaining energy by the rackets moving toward each other. Hence every collision is a head-on collision in contrast to the second-order Fermi process, in which the stochastic motions of the magnetic mirrors induce head-on and following collisions, so that the energy gain is of second order, as for a moving particle the number of head-on collisions dominates. The process is not very efficient per se, as the energy gain of a particle is only a small fraction of order V _s ∕ c, where V _s is the shock speed with respect to a fixed Eulerian frame. However, it is a probabilistic process, in which the probability that a particle gets reflected back and forth many times decreases with increasing energy (i.e., number of crossings). Let E =  β E ₀ be the energy after one collision, and let P be the probability of a particle remaining in the acceleration process, so that after one collision the number of particles is N =  PN ₀. Hence after k collisions, N =  N ₀ P ^k , and E =  E ₀ β ^k . Since then \(k =\log (E/E_{0})/\log \beta\), we have \(N/N_{0} = (E/E_{0})^{\log P/\log \beta }\), and the differential energy spectrum is a power law, \(dN \propto (E/E_{0})^{(\log P/\log \beta )-1}\), as observed. Let U be the downstream fluid velocity in the laboratory frame, which for a strong shock is given by \(U = 3/4\,V _{s}\). The energy gain per roundtrip is \(\varDelta E/E = 4/3\,U/c = V _{s}/c\) (averaging over pitch angles) and hence \(\beta = 1 +\varDelta E/E = 1 + V _{s}/c\), while the escape probability of the particles is related to the downstream and shock velocity, P _esc ∼  V _s ∕ c, due to advection away from the shock. The probability to stay in the acceleration process is \(P = 1 - P_{\mathrm{esc}} = 1 - V _{s}/c\), so that \(\log P \approx -V _{s}/c\), log β ≈  V _s ∕ c, and \(\log P/\log \beta \approx -1\) (see Bell  1978). Therefore the power-law index in the energy spectrum is \(\gamma \approx \log P/\log \beta - 1 \approx -2\), close to the observational value of γ ≈ −2. 7. The difference can be attributed to the energy-dependent diffusive transport of CRs, with a diffusion coefficient κ ∝  E ^−0. 7, because scattering becomes less efficient for particles with higher energies particles, as there exist fewer and fewer waves to scatter off.

The maximum energy reachable in DSA is E _max ∼ ( q∕ c)  V _s   B  L, which is about 3 × 10 ¹⁴ eV for protons in a SNR with a microgauss field and L ∼ 10 pc, making use of the fact that the highest energies are obtained in young SNRs, when the shock is still strong, and V _s  ∼ 10 ⁴ km∕s. This agrees remarkably well with the location of the knee in the spectrum.

As the CRs propagate by diffusion and advection from their sites of origin to the heliosphere, they suffer several energy losses. While protons and nuclei lose energy by spallation processes and π ⁰-production via collisions with ISM protons, synchrotron and inverse Compton losses (mainly with cosmic microwave background photons) affect mostly the electrons because of their lower mass. Ionization and bremsstrahlung losses are usually negligible in comparison for particles with energies above 100 MeV.

There are several phenomena that define the boundaries of the heliosphere, all of which are illustrated in Fig.  5 together with typical particle populations and the magnetic field structure (Opher  2016; Zank  2015). Toward the direction of the incoming interstellar medium, which has a density of 0. 2 cm ⁻³, the heliosphere is flattened. In the opposite direction – at the so-called heliotail – the structure is considerably more complex. So far, however, in situ measurements exist only from the nose direction, while the heliotail is uncharted territory.

In general, the extent of an astrosphere depends on the respective pressures of the stellar wind and that of the surrounding interstellar medium. Notably, the latter contains contributions from the ionized and neutral gas components, the (non-thermal) galactic cosmic ray (GCR) population, and the magnetic field. Since all these components are in approximate equipartition, neither may be neglected (Biermann and Davis  1960; Schlickeiser  2012).

The medium that streams toward the heliosphere is partially ionized and contains GCRs, i.e., a non-thermal particle population. Embedded into it is a weak magnetic field with a strength of 2–4 μG. The outermost heliospheric boundary is formed by the deceleration of the local interstellar medium, through which the heliosphere moves with approximately 26 km/s. That speed is close to the local sound speed, so that it is currently unknown precisely if the flow is super- or subsonic (McComas et al.  2012; Scherer and Fichtner  2014). Accordingly, either a bow shock or a bow wave is possible. The compressed and heated interstellar particles then act as partners for charge-exchange processes with inflowing neutral atoms, thereby forming the so-called hydrogen wall . Using Ly α absorption, this wall has indeed been detected. Unfortunately, however, the presence of the hydrogen wall does not prove the existence of a bow shock because it is also compatible with a bow wave. For the time being, therefore, that question remains unanswered, and more precise determinations of the sound speed will be required. In the region containing this so-called very local interstellar medium with a maximum distance of 0.01 pc from the Sun, the influence of heliospheric particles such as pickup ions is still palpable.

The separation of the local interstellar medium from the solar wind plasma has been predicted to be a tangential discontinuity. Its location depends on the orientation relative to the local interstellar flow but, in upwind direction, has a distance from the Sun of at least 130 au. In 2012, observations of low-energy particles led to the impression that Voyager 1 had crossed the heliopause and entered the (very) local interstellar medium. There was, however, some debate as to whether Voyager 1 has indeed encountered the heliopause, mainly because the magnetic field direction has remained unchanged, in contrast to expectations. Therefore, the term “heliocliff” was coined to indicate a transition layer before the actual heliopause . Despite additional measurements, however, it is still not unanimously accepted that Voyager 1 is in the outer heliosheath . While a variety of promising theoretical models – including magnetic reconnection, turbulence, and magnetic bubbles – have been proposed, it will require additional observational evidence to resolve this issue. As of 2016, it is expected that Voyager 2 will soon reach the heliopause.

The highly turbulent region inside the heliopause that contains basically subsonic plasma flows is called the (inner) heliosheath . Here, the solar wind ions are diverted to stream toward the heliotail , thus remaining inside the heliopause. Note that, analogously, the region outside the heliopause is sometimes called the outer heliosheath, even though, strictly speaking, this term would require the presence of a bow shock, so there is no general consensus. Voyager 1 has revealed a variety of surprises, the most prominent of which is the heliosheath being surprisingly thin. Almost all models had overestimated the thickness by nearly a factor of two. As of yet, no all-encompassing solution to this mystery has been found, even though it is believed that magnetic reconnection may play an important role (Opher  2016). Reconnecting magnetic fields and turbulence may extract energy from the heliosheath, thereby reducing the pressure and leading to a region that is smaller than predicted by traditional models. Further difficulties arise from the observations of possibly compressible, magnetized turbulence that has been observed both by Voyagers 1 and 2. The understanding of particle observations is thus made significantly more complicated. This is even more so as the turbulence is partly supersonic, which results in a series of shocks propagating from the termination shock into the heliosheath.

In December 2004, Voyager 1 found the termination shock to be located at a radial distance from the Sun of 94 au (Richardson  2013). The distance at which Voyager 2 encountered the termination shock in August 2007, however, was only 84 au. This asymmetry can be explained either by the oblique angle of the interstellar magnetic field (see below) or by oscillations induced by the solar wind. However, these are responsible only for variations of up to 2–3 au. In each case, the shock crossing was indicated by a sudden decrease of the flux of low-energy particles while, at the same time, the flux of high-energy particles increased. Since the magnetic field direction did not change, the shock, which is collisionless in nature and surprisingly thin, appears to be almost perpendicular (Balogh and Treumann  2013). It was believed for a long time that the termination shock marks the transition of the solar wind from the supersonic to the subsonic regime. However, there are in fact observations indicating that, instead, pickup ions (see below) are heated, while the solar wind ions continue to stream at a supersonic speed (Zank  2015). Despite the fact that the termination shock is the most prominent boundary of the solar system, it is relatively weak in terms of the shock compression ratio, which is between 1.4 for the particle density and 1.7 for the magnetic field (Burlaga et al.  2013).

There is a variety of processes by which the interstellar medium affects the heliosphere , including forming its very shape. Likewise, the heliosphere also has – albeit to a lesser degree – a certain influence on the very local interstellar medium, mainly through the magnetic field.

The termination shock cannot be directly observed with Earth-based instruments or from (near-Earth) space. Indirectly, the anomalous component of the cosmic ray particles (ACRs) – originally named after their peculiar feature in the energy spectrum – as well as energetic neutrals may represent an exception. In contrast to ionized particles, neutral interstellar atoms stream unhindered into the heliosphere, owing to the collisionless nature of the limiting shocks. These particles are therefore particularly suited to characterize the conditions in the outer the heliosphere and the local interstellar medium. Some fraction of these neutrals – mostly hydrogen – undergo multiple collisions with energetic protons of the solar wind, thereby exchanging one or more electrons (Balogh and Treumann  2013). Initially, these pickup ions are cold but are typically heated up to form a supra-thermal population with energies of 10–100 MeV per nucleon. The particle spectrum therefore is a power law with an index between 1.5 and 3.2, depending on the location close to the termination shock or deep inside the heliosheath. Originally, ACRs (see, e.g., Fichtner  2001) were thought to be accelerated at the termination shock. So far, however, evidence for this scenario has been found only for particles with moderately low energies, which have thus been named termination shock particles (TSP). In contrast, Voyager 1 found that the intensity of ACRs kept rising without an intensity peak at the termination shock. This fact represents a challenge for the widely accepted picture of diffusive shock acceleration and hints at a source inside the heliosheath.

While the heliopause generally hinders the interstellar plasma from streaming into the heliosphere, neutral particles cross the heliopause to become pickup ions . In these charge-exchange processes, energetic neutral atoms (ENA) are created which, in turn, may move inward and thus can be detected at Earth. Based on the theoretical understanding of their origin, ENAs can therefore be used to probe the plasma conditions at the outer heliosphere where they had been created. In particular, a bright narrow band of energetic neutral atoms with energies around 0.7–1.7 keV has been discovered by the Interstellar Boundary Explorer (IBEX), which was therefore called the “IBEX Ribbon” (McComas et al.  2014). There are indications that it is aligned with the interstellar magnetic field as well as with solar latitude, thus pointing at a direct connection with the solar wind. Based on these findings, a number of theoretical models have been brought forward, some of which are based on multiple charge-exchange processes. In particular, aside from the so-called primary or classic pickup ions, a second population of higher-energy ions (dubbed “ACR pickup ions”) is required that are transported back to the termination shock (e.g., McComas et al.  2014; Siewert et al.  2012). However, the precise mechanism by which the Ribbon is produced is currently still under active debate.

Even though Voyager 1 seems to have crossed the heliopause , the magnetic field still appears to be that of the heliosphere. Therefore, direct observations are mainly important to measure the strength and the fluctuations of the local magnetic field. But the direction of the local interstellar magnetic field can be inferred indirectly to account for the asymmetries observed in situ by the Voyagers. However, to reconstruct the large-scale structure of the surrounding interstellar magnetic field, mainly remote observations are being used. For instance, anisotropies in TeV cosmic rays from the ISM can be taken to infer the orientation of the local interstellar magnetic field (Schwadron et al.  2014), which is in agreement with results based on the polarization of starlight. In addition, the IBEX Ribbon appears to be oriented in the direction perpendicular to the ambient interstellar magnetic field. There is no reason why the interstellar magnetic field should be aligned with the flow direction of the ISM, and indeed it is not. Accordingly, the shape of the heliosphere is not symmetric which, in turn, accounts for some of the differences in the observations of Voyagers 1 and 2. It is these findings that connect observations of GCRs to the structure of the heliosphere (Desiati and Lazarian  2013), thereby emphasizing the importance of high-energy particles not only to probe the conditions of supernova remnants but also as a diagnostic tool for the heliosphere.

While it may sound like a trivial statement, the Sun is the origin of the heliosphere for mainly two reasons (Zurbuchen  2007). First, a tenuous plasma – dubbed the solar wind – streams radially outward with a supersonic speed of 250–2200 km/s. It consists of two distinct components, the slow and the fast solar wind, with typical speeds of 400 and 750 km/s, respectively. The mass loss rate of 3 × 10 ⁻¹⁴ M _⊙ per year is low compared to that of massive OB stars, which can be nine orders of magnitude higher and thus amount up to 10 ⁻⁵ M _⊙ per year. Still, the solar wind is the basis for all phenomena discussed above, and more. For instance, shocks in the solar wind – so-called co-rotating interaction regions – are known to affect the spectrum of GCRs measured on Earth (Kóta and Jokipii  1995). Second, the solar magnetic field extends outward in the form of an Archimedean spiral – also known as the famous Parker ( 1958) spiral . The polarity reversals, which follow the 11-year solar cycle , affect the particle spectra in a variety of ways (Potgieter  2013). For instance, the cosmic ray energy spectrum becomes softer for a positive solar polarity, even though it has to be noted that particles of opposite charge will react differently. In addition, also the dynamics of the warped heliospheric current sheet , which separates the two hemispheres of opposite polarity, has an influence on the particle populations. All of these effects are summarized under the term solar modulation .

In conclusion, particle measurements in the heliosphere – and in particular near Earth – require a thorough understanding of the solar wind. Charge-exchange processes creating pickup ions, acceleration processes creating ACRs, and the solar modulation need to be taken into account. The interpretation of what is recorded inside the heliosphere does not directly reflect what will be found in the interstellar medium. However, our understanding of the interaction regions grows so that neutral particles as well as cosmic rays are increasingly recognized as being messengers of interstellar conditions. Nevertheless, given the history of the Voyager observations, it is most likely that further surprises will follow as the Voyagers go where no space probe has gone before.