Interstellar Neutrals, Pickup Ions, and Energetic Neutral Atoms Throughout the Heliosphere: Present Theory and Modeling Overview

Interstellar neutrals (ISNs), pick-up ions (PUIs), and energetic neutral atoms (ENAs) are fundamental constituents of the heliosphere and its interaction with the neighboring interstellar medium. Here, we focus on selected aspects of present-day theory and modeling of these particles. In the last decades, progress in the understanding of the role of PUIs and ENAs for the global heliosphere and its interaction with very local interstellar medium is impressive and still growing. The increasing number of measurements allows for verification and continuing development of the theories and model attempts. We present an overview of various model descriptions of the heliosphere and the processes throughout it including the kinetic, fluid, and hybrid solutions. We also discuss topics in which interplay between theory, models, and interpretation of measurements reveals the complexity of the heliosphere and its understanding. They include model-based interpretation of the ISN, PUI, and ENA measurements conducted from the Earth’s vicinity. In addition, we describe selected processes beyond the Earth’s orbit up to the heliosphere boundary regions, where PUIs significantly contribute to the complex system of the global heliosphere and its interaction with the VLISM.


Introduction
Interstellar neutrals (ISNs), pick-up ions (PUIs), and energetic neutral atoms (ENAs) are fundamental constituents of the heliosphere and its interaction with the neighboring interstellar medium (ISM). Here, we focus on selected aspects of present-day theory and modeling of these particles. However, the accompanying parts of this book greatly complement the material presented here, starting from the historical overview by Zank et al., discussion of ISNs and ENAs observations by Galli et al. and Dialynas et al., and overview of PUI measurements by Zirnstein et al., to list just a few. Atoms from the local interstellar medium (LISM) enter the heliosphere and interact with the solar environment inside. They are ionized by the solar wind (SW) and the ultraviolet radiation and next are picked up by the ambient magnetic field and propagate outward as PUIs. PUIs are born in the supersonic SW inside the termination shock (TS) and in the heliosheath between the TS and heliopause (HP) due to the charge exchange of the interstellar atoms and SW ions. The PUIs re-neutralized by charge exchange become ENAs, which constitute the so-called globally distributed flux (GDF). The GDF ENAs are just one population of the ENAs present in the heliosphere. Unfortunately, due to the volume limitation for the paper, other ENA populations, like, e.g., planetary ENAs or ribbon ENAs are outside our scope. All three populations of particles, ISNs, PUIs, and ENAs are present throughout the heliosphere. However, their role in the physics of the heliosphere varies with distance to the Sun (Zank 2016). Figure 1 shows the SW-LISM interaction pattern with the TS, HP, and bow shock (BS) according to the MS-FLUKSS model ; see Sect. 3), which is one of several heliosphere concepts, which are extensively described in the accompanying chapter by Kleimann et al.
A region of the LISM surrounding the Sun that is modified or mediated by heliospheric processes or material is commonly called the very local interstellar medium, VLISM (Zank 2015). Thus, VLISM is the LISM affected by the presence of the heliosphere. Fraternale and Pogorelov (2021) further extended this definition by referring to the LISM modification by the mere presence of the heliosphere, regardless of what physical processes  (Zhang and Pogorelov 2016) are involved or which LISM properties are affected. This is especially important for the polar and downwind directions. The part of the heliosphere in close vicinity to the Earth is an important region for the study of the heliosphere because the majority of the measurements have been done from Earth's vicinity. Also, the solar ionizing environment affects the interstellar particles significantly at these distances and cannot be neglected. At distances greater than a few astronomical units (au) from the Sun, the solar influence decreases, PUIs start to contribute substantially to the pressure balance, heat the SW, and slow it down. In the transition region, the SW mixes with the interstellar medium, and the population of PUIs can be subdivided into those transmitted or reflected at the TS and those being created via charge exchange of ISNs inside the inner heliosheath . They contribute to the ENAs observed from inside the heliosphere next to the ENAs from other sources, like, e.g., those creating the ribbon ).
The ENAs observed are hydrogen ENAs which are the result of charge exchange between SW protons and interstellar H (see chapters by Galli et al. and Dialynas et al.). However, there are also possibilities for the creation of ENAs by interaction with the heavier mass elements, like, e.g., He, Ne, O (see, e.g., Swaczyna and Bzowski (2017) and extensive references on this topic therein). As concluded by Swaczyna and Bzowski (2017) the so-called heavy ENAs could bring information on the global structure of the heliosphere based on the estimated spatial distribution of their fluxes and the energization processes at the TS. However, to successfully observe ENAs of masses greater than hydrogen, a detector with advanced mass resolution capabilities for low intensities is needed.
In the last decades, progress in the understanding of the role of PUIs and ENAs for the global heliosphere and its interaction with VLISM is impressive and still growing. The increasing number of measurements allows for verification and continuing development of the theories and models. We aim to provide a balanced context for these topics and thus we organize the content into the following parts: we start with a brief presentation of the current state of the kinetic model description (Sect. 2). Next, we present the role of PUIs in the MS-FLUKSS model description (Sect. 3). This part is followed with discussions of a fluid model of PUIs at TS (Sect. 4), shock crossing by non-Maxwellian plasma (Sect. 5), and the role of electron impact ionization downstream of the SW TS (Sect. 6). Next, we present heliosphere model solutions proposed by the SHIELD DRIVE Science Center (Sect. 7). In the next part of the chapter, we discuss research fields in which interplay between theory, models, and interpretation of measurements reveals the complexity of the heliosphere and its understanding. In Sect. 8, we discuss topics necessary to interpret ISN, PUI, and ENA measurements conducted from the Earth's vicinity, where the solar environment significantly modifies fluxes of these interstellar particles. They include ionization rates (Sect. 8.1), ISN and PUI distributions throughout the heliosphere (Sect. 8.2), the secondary population of ISNs (Sect. 8.3), PUIs close to the solar corona (Sect. 8.4), terrestrial ENAs (Sect. 8.5), and inner source PUIs (Sect. 8.6). While in Sect. 9 focuses on processes beyond Earth's orbit up to the inner heliosheath (IHS), where PUIs significantly contribute to the complex system of the global heliosphere and its interaction with the VLISM. We focus on SW mediation by PUIs (Sect. 9.1), waves due to newborn interstellar PUIs (Sect. 9.2), acceleration of PUIs throughout the heliosphere (Sect. 9.3), and shocks in the IHS (Sect. 9.4). We briefly summarize the chapter in Sect. 10.

Kinetic Model Description
The mean free path (mfp) for charge exchange of ISNs is comparable with the size of the heliosphere (Izmodenov et al. 2000) and the ISNs penetrate through the heliopause (HP) reaching the vicinity of the Sun. The ISNs coming from the undisturbed VLISM are called the "primary" component of the ISNs. Their distribution inside the heliosphere is calculated with the use of the Boltzmann equation for the atom distribution function f (r, v, t) dependent on position (r), velocity (v), and time (t), where m is the mass of an atom: ∂f (r, v, t) ∂t The loss term on the right hand side of Equation (1) is due to ionization along the particle trajectory, where β is the ionization rates, t time, r distance from the Sun, and b heliolatitute (more in Sect. 8.1). For H the main contributions to the ionization rate are charge exchange with SW protons and photoionization (see Sect. 8.1). The force acting on the particle is the Sun's gravity and, in the case of H, radiation pressure (e.g., Rucinski and Bzowski 1995;Kowalska-Leszczynska et al. 2018): where G is the gravitational constant, M s is the mass of the Sun, and μ is the ratio of repulsion from radiation pressure over gravitational attraction (μ = |F rad |/|F g |). In general, the parameter μ depends on time (t ), heliolatitude (b), and the radial component of the particle velocity (v r ). Due to the low mass of H and the resonant interaction of H with the strong solar Lyman-α line, the precise determination of μ is crucial for the distribution of the ISN H in the vicinity of the Sun (e.g., Rahmanifard et al. 2019;Katushkina et al. 2021b,a;Kowalska-Leszczynska et al. 2022). However, the parameter μ is approximately zero for He and O due to their heavier mass compared to H and different line properties (Ruciński 1985). Under the assumption that the thermal speed of the atom is much less than its bulk flow speed relative to the Sun, the distribution function for the ISNs in VLISM is described by a cold gas approximation (e.g., Fahr 1968a; Blum and Fahr 1970;Holzer 1972;Axford 1972). However, the thermal and bulk flow speeds of atoms in the VLISM are comparable and the gas temperature in the VLISM is an important quantity for the interpretation of the heliospheric observations. Thus, the cold model was extended by Fahr (1971), Thomas (1978), Wu and Judge (1979) to the so-called hot model, in which the source distribution function is assumed as a Maxwell-Boltzmann distribution. The particle trajectories are inferred from solving Kepler's equation and the moments of the distribution function give number density, velocity, and temperature. The hot model commonly assumes that radiation pressure and losses due to ionization with the solar environment are constant in time. Fahr et al. (1987) and Rucinski and Bzowski (1995) further developed the hot model by including the variations in time for the radiation pressure and ionization losses. Present-day models (e.g., Katushkina et al. 2015;Sokół et al. 2015a) also take into account spatial variations of the ionization losses inside the heliosphere.
The kinetic Equation (1) is a linear partial differential equation that is solved by the method of characteristics. The solution of this equation is where f ∞ (r ∞ , v ∞ ) is the velocity distribution function of an ISN atom at the outer boundary. The quantities r ∞ , v ∞ and t ∞ are the position, velocity, and time when the atom crosses the outer boundary. The integration in the exponential is performed along the atom's trajectory, which is the characteristic of Equation (1). The velocity distribution function of ISN atoms at the outer boundary can be obtained from the results of the global model of the heliospheric boundary (i.e., Alexashov 2015, 2020).

Kinetic Modeling of the Pickup Proton Distribution
The distribution of H PUIs in the heliosphere can be calculated using the plasma and H atom distributions obtained in the frame of the global model of the heliosphere, such as the 3D time-dependent kinetic-MHD model by . The method was developed by Malama et al. (2006) and employed in Chalov et al. (2016) and Baliukin et al. (2020). The velocity distribution of pickup protons in the SW rest frame is assumed isotropic following theoretical estimates (Vasyliunas and Siscoe 1976) and observations . The velocity distribution function f pui (t, r, v) is then given by: where v = V(r, t) + w, v and V are the velocity of pickup protons and bulk the velocity of the SW plasma in the heliocentric coordinate system, w is the velocity of the pickup proton in the SW rest frame, and (w, θ , ϕ) are coordinates of w in the spherical coordinate system. The kinetic equation for f * pui (t, r, w) in the general form is (e.g., Isenberg 1987;Chalov et al. 2003): where D(t, r, w) is the velocity diffusion coefficient, but with spatial diffusion neglected. In Equation (5), the source (S + ) and sink (S − ) terms are: Here, f H (t, r, v H ) is the hydrogen velocity distribution function and β is the total ionization rate. In this approach, H atoms experience charge exchange both with SW protons and PUIs, where the effective charge exchange cross section depends on the relative atom-proton velocity (e.g., Lindsay and Stebbings 2005). The production rate of ENAs in Equation (8), ν H , due to the charge exchange of protons with H atoms, is defined by For relatively high energies, such as |v − v H | ≈ v, the production rate of ENAs can be approximated as ν H (t, r, v) ≈ n H (t, r)vσ ex (v). The process of velocity (energy) diffusion, which leads to stochastic acceleration of PUIs, is ignored in the approach proposed by Baliukin et al. (2020), i.e. D = 0 is adopted. This corresponds to a quiet SW when the magnetic field fluctuation level is low (Chalov et al. 2003). Nevertheless, the process of velocity (energy) diffusion is connected with effective pitch-angle diffusion, which leads to isotropization of the velocity distribution function, and therefore should be taken into account.
In the case without diffusion, Equation (5) becomes a first-order linear differential equation that can be solved by the method of characteristics. The characteristic is the SW particle trajectory along which the PUI velocity changes: The solution of Equation (5), f * pui (t, r, w), includes the loss of PUIs due to charge exchange along trajectory. The trajectory of PUIs is reconstructed backward in time from (t, r(t), w(t)) to point (t 0 , r(t 0 ), w(t 0 )), which is close to the Sun, where it can be assumed f * pui (t 0 , r(t 0 ), w(t 0 )) = 0. The kinetic moments of the velocity distribution function of PUIs, the number density, and the temperature, at point (t, r) ∈ R 4 are the following: where m p is the mass of proton, and k B is the Boltzmann constant. The pressure of PUIs is calculated as p pui = n pui k B T pui .

Kinetic Model Jump Condition at the TS
The Liouville's theorem, which states the phase space flow conservation, together with conservation of the first adiabatic invariant (magnetic moment), and assumption of the weak scattering leads to the following condition on the velocity distribution function at the shock Siewert 2011, 2013a): where In Equation (12), f * pui,u and f * pui,d are velocity distribution functions upstream and downstream the TS. The upstream shock-normal angle ψ(t, r) (between the magnetic field and normal to the shock surface) and the shock compression factor s(t, r) = n d /n u depend on the position r and moment t of the TS crossing. As observed by Voyager 1 (for TS-2 crossing; see Richardson et al. (2008)) the compression ratio is 2.4. As follows from the jump condition (12) the downstream/upstream ratios for the moments can be obtained: It should be noted that some PUIs can experience reflections at the shock front due to an abrupt change in the magnetic field. These ions gain energy from their drift motion along the TS in the direction of the induced electric field. The PUIs can also experience the pitchangle scattering upstream and downstream of the TS, which allows transmitted particles to return to the shock front, so multiple reflections can occur. The process of reflection leads to the anisotropy of the velocity distribution (in the SW rest frame) in the vicinity of the TS, and, therefore, the transport equation for the anisotropic velocity distribution function should be solved (see, e.g., Chalov et al. (2016)).

Additional Energetic Population of Pickup Protons
An additional energetic PUI population is seen in the form of the energetic "tails" in the velocity distribution . It is produced by the processes such as shock-drift acceleration or reflection from the cross-shock potential ("shock-surfing" mechanism, see also Sect. 9.3). Based on the kinetic model (Sect. 2.1), the velocity distribution function can be calculated everywhere downstream the TS. This distribution is the so-called filled shell f * sh (t, r, w), and its critical velocity where V is the SW bulk velocity, V H is H atom bulk velocity, t TS and r TS = r(t TS ) are the moment and position of a PUI crossing the TS. For w ≥ w c the power-law distribution f * tail (t TS , r TS , w) ∝ w −η is assumed. If the additional parameter ξ , which is the density fraction of PUIs of the "tail" distribution, is introduced, the velocity distribution function downstream the TS is where and f tail (t TS , r TS , w) = 0 for w < w c , n pui,d is the number density downstream the TS. Equation (15) implies that η > 3, otherwise the number density of the PUIs in the "tail" distribution is infinite.
In the IHS, the solution of the kinetic equation (Equation (5)), which can be derived using the method of characteristics, with the boundary condition downstream the TS (Equation (14)) is employed. The first and second terms of this condition can be associated with transmitted and additional energetic populations of PUIs, respectively. This partitioning of PUIs into sub-populations implies the conservation of the total number density, while the pressure balance is not satisfied. By introducing the energetic "tail", additional energy is injected into the system, so this approach is not completely self-consistent. This additional energy can be "pumped" by the interaction of PUIs with the TS or produced by fluctuations of the heliospheric magnetic field.
The parameters ξ and η, which were artificially introduced in the approach described above, depend on the local TS properties, since the efficiency of acceleration processes depends on the shock-normal angle ψ and shock compression factor s, in particular. For the specified dependencies of parameters ξ and η on the position at the TS the velocity distribution function of pickup protons in the heliosphere can be calculated.
Based on the numerical kinetic model and observations of the ENA fluxes from IHS by the IBEX-Hi instrument, Baliukin et al. (2022) characterized the pickup proton distribution and provided estimations on the properties of the energetic pickup proton population downstream of the TS. They have performed the parametric study by varying ξ and η in different directions (upwind, downwind, flanks, solar poles) and in a wide range of values to minimize the difference between the model results and data. The quantitative estimates for the parameters of the energetic population of pickup protons obtained in their work can be useful for testing and verifying hybrid models that simulate the acceleration of protons at the heliospheric termination shock. The results of the ENA flux simulations based on these estimations are consistent with IBEX-Hi observations and discussed in Sect. 2.5.

Modeling of the ENA Fluxes
ENAs, which are born in charge exchange between the SW protons and H atoms in the IHS, form the GDF. The directional differential ENA flux in the solar inertial reference frame is the line-of-sight integral: where is the survival probability, t obs is the moment of observation, r obs is the position of the observer, LOS is the unit vector in the line of sight direction, v = v · LOS is the velocity of an ENA, m H is the mass of H atom, and f = f sw + f pui is the velocity distribution of protons  Baliukin et al. (2022) in the inner heliosheath (the sum of core SW proton distribution, which is assumed to be Maxwellian, and PUI distribution). The integration is performed along the ENA trajectory in the inner heliosheath, with ds = vdt being the differential path length. The values s TS and s HP correspond to atom intersections of the TS and HP, respectively.

Kinetic Model Results for ENAs and Comparison with IBEX-Hi Data
For the comparison with modeling results, the data sets of GDF observed by the IBEX-Hi detector during the IBEX mission were used. The data are presented by  and available via the public Data Release 8. 1 To calculate the differential fluxes measured by IBEX-Hi at different energy channels, the energy transmission of IBEX-Hi electrostatic analyzers should be taken into account since it leads to the re-distribution of fluxes between the energy channels and spreading of the observed spectrum. The Compton-Getting and survival probability corrections were applied for the IBEX-Hi data, so in the calculations, both the relative motion of the IBEX-Hi spacecraft and the ionization losses of ENAs in the region of the supersonic SW are not taken into account. The model maps of ENA fluxes were obtained for the lines-of-sight when the sensor views the heliosphere in the ram direction (i.e., in the direction of spacecraft motion), as it was done by . Figure 2 presents the comparison of the IBEX-Hi data at the top five energy channels with the ENA flux maps obtained by Baliukin et al. (2022) based on the time-dependent heliospheric model by . The model fluxes were averaged over 2009 -2013. The first and third columns present the model results, the second and fourth columns show the measurements. From the comparison of the simulation results with the data the following conclusions can be made: 1. The model reproduces the geometry of the multi-lobe structure seen in the IBEX data.
There is a good qualitative agreement between the model results and the observed fluxes. 2. A single structure of enhanced fluxes from the tail region is seen in the IBEX data at low energy channels (∼0.71, 1.1, 1.74 keV). At high energies, this structure "splits" into two parts (to the north and south from the solar equatorial plane). The model results qualitatively reproduce such "splitting" behavior of the lobes. The lines-of-sight within these lobes of enhanced ENA fluxes intersect the regions of the inner heliosheath, where (a) the fast SW, initially emitted from the solar polar coronal holes and propagated to the heliospheric tail, is collimated, and the plasma velocity and temperature are high, and (b) the heliosheath thickness is large. 3. The presence of the low flux areas from the flanks of the heliosphere (so-called Starboard and Port lobes) observed in the IBEX data is also seen in the model results. These low fluxes lobes are located in the vicinity of the solar equatorial plane, where the slow SW dominates and the thickness of the heliosheath is small Zirnstein et al. 2016). 4. There are a few differences between the model calculations and the IBEX-Hi data: (a) for quantitative agreement with the data, the model fluxes should be significantly scaled (multiplied by 1.54), which may be an indicator of a lower H number density used in the heliospheric model than it is in reality; (b) the simulated fluxes at ∼0.71 keV energy channel are lower compared to the data. These differences may be induced by several assumptions made by Baliukin et al. (2022), such as (a) the neglect of the velocity (energy) diffusion, (b) the isotropic form of the distribution function throughout the heliosphere, and (c) the weak scattering at the TS.

PUIs in MS-FLUKSS
A number of one-plasma-fluid models (Baranov and Malama 1993;Pauls et al. 1995;Pogorelov et al. 2006;Heerikhuisen et al. 2006Heerikhuisen et al. , 2007Izmodenov et al. 2005Opher et al. 2006) take into account the effect of PUIs by assuming that they are in thermal equilibrium with the background plasma. Although this is not correct, the conservation laws of mass, momentum, and energy for the mixture of electrons, ions, and PUIs are satisfied approximately. According to , such approaches are correct only if the heat conduction tensor and the dissipation coefficient due to ion-PUI interaction are zero. This happens when PUIs are coupled with the thermal ions and the scattering time is infinitely small. Isenberg (1986) shows that the effect of PUIs can be quantified if they are treated as a separate fluid. The time-dependent simulations of the supersonic SW by Usmanov and Goldstein (2006), Usmanov et al. (2012) include PUIs as a separate fluid and also turbulence effects. However, these models still assume instantaneous isotropization, i.e., with PUI described as a filled shell, thus neglecting the nearly-isotropic character of the PUI distribution function when a finite scattering time is included .
In principle, any PUI fluid model is an approximation of their kinetic behavior (Malama et al. 2006), where it is assumed that PUIs are isotropic away from discontinuities. Recently, Gamayunov et al. (2012) investigate the evolution of the PUI distribution function in a pitchangle-averaged approximation together with the PUI-generated waves that heat the thermal SW ions. The extension of a 2-fluid SW model developed for the supersonic flows is not straightforward, since the pressure equation for PUIs used in such models is not valid, e.g., across the TS. Instead, some boundary conditions should be developed (Chalov and Fahr 1996;Zank et al. 1996b;Fahr et al. 2000b;Fahr and Chalov 2008;Fahr et al. 2012;Fahr and Siewert 2013b) since, as revealed when Voyager 1 crossed the TS (Richardson et al. 2008), most of the SW kinetic energy was absorbed by PUI.
Although the SW plasma consists of the thermal SW ions and PUIs, in Multi-Scale Fluid-Kinetic Simulation Suite (MS-FLUKSS) the MHD system is solved in the conservation-law form because this satisfies the shock boundary conditions efficiently. This requires solving the MHD system for the mixture of SW ions and electrons, and PUIs. To separate between the SW ions and PUIs, one can supplement the system with the continuity and pressure equations for PUIs . The distribution functions of protons, f p , and PUIs, f PUI to calculate the source terms in the MHD system are usually assumed to be Maxwellian. This is not true for PUIs. Heerikhuisen et al. (2019), DeStefano and Heerikhuisen (2017) studied the effect of the distribution function on the charge exchange source terms considering Lorentzian (kappa) distributions with different κ's for PUIs. The kappa distribution decreases the contribution of the most probable state while increasing the number of ions in the so-called "energetic tails". The conclusion derived from DeStefano and , DeStefano (2019) is that the dependence of charge-exchange cross-sections (σ ex ) on energy should be preserved when the integration is performed. This becomes crucial for κ approaching the minimum allowed value of 3/2. For κ → ∞, when the distribution function becomes Maxwellian again, the dependence of σ ex on energy becomes less critical for the integral evaluation. Otherwise, the integrals diverge at κ = 2 (see also, Scherer et al. (2017Scherer et al. ( , 2019, and Perri et al. this volume). MS-FLUKSS code avoids such divergences by performing the integration more accurately, avoiding exaggeration of the effect of PUIs on the heliosphere . That happens when PUIs are not treated separately from SW protons and a Maxwellian distribution function (with substantially higher temperature) is implied for the ion mixture.
There are essentially three approaches to include PUIs, which form the basis of the MS-FLUKSS models: (1) PUIs are not treated as a separate plasma component, but the conservation of mass, momentum, and energy for the mixture is preserved; (2) an isotropic distribution function for PUIs is chosen; PUIs are assumed to be co-moving with the SW ions, and the continuity and pressure equations are solved to describe the PUI motion; and (3) f PUI is assumed isotropic beyond shocks, at which it is intrinsically non-isotropic, and proper boundary conditions describing the PUI transition across the TS, or other shocks, treated kinetically, are applied.
The system describing the plasma mixture (thermal protons, PUIs, and electrons), neutral atom populations, and PUIs is presented below. Subscripts for charge exchange source terms are identified by the letter a, which is a fixed integer for the corresponding population of neutral atoms, p for thermal protons, and I for pickup protons: Neutrals Pick-up ions Here p is the pressure, ρ the mass density, u the bulk velocity vector, B the magnetic field vector, e the total energy density (includes the magnetic energy), p * the total pressure of the mixture (including magnetic energy), and I is the identity tensor. The source terms D, M, E, and P are written out in . The subscript refers to the mixture properties. It is reasonable to assume that there are no pickup electrons. Thus, introducing the number density n, we obtain p = n e kT e + n p kT p + n I kT I , n e = n p + n I = n, which reduces to Pogorelov et al. (2016) assume a multi-fluid model for H atoms. The system of MHD equations is solved using a Godunov-type method, where the fluxes of mass, momentum, energy and magnetic field are obtained as solutions of an MHD Riemann problem in the linear approximation. The scheme has second order of accuracy in space and time. Both MHD and Euler equations are solved in the framework of MS-FLUKSS (Pogorelov et al. 2008). Figure 3 shows the plasma density distributions (dimensionless units with the scale of VLISM plasma density) in the meridional plane (the plane formed by the Sun's rotation axis) for the single-ion-fluid model (left panel) and the case when the PUI and thermal ion fluids behave as co-moving, but distinguishable components (right panel). The HP moves inward towards the Sun and the TS moves outward away from the Sun with the inclusion of PUIs as a distinct component, i.e., the heliosheath width becomes quite significantly narrower when PUIs are treated distinctly from the thermal SW plasma. This result is consistent with that in Malama et al. (2006) and addresses the observation made by Voyager 1 that the HP location was closer than expected (Stone et al. 2013;Burlaga and Ness 2014). In addition, a comparison of MS-FLUKSS simulation based on the model with PUIs governed separately by the continuity and pressure equations and the turbulence model of Breech et al. (2008) with New Horizons (NH) observations is presented in Fig. 4.

Two-Fluid Model of PUIs at TS
Voyager 1 and 2 crossed the TS and entered the IHS at 94 and 84 au, respectively. The Voyager 1 crossing of the TS occurred during a data gap, and thus it was incapable of determining its detailed structure. The TS crossing by Voyager 2 returned many surprises and, most importantly, established the importance of non-thermal energetic particles. In contradiction to many predictions from theory and global MHD models, Voyager 2 showed that the thermal gas downstream of the TS is supersonic and much colder than expected. Richardson et al. (2008) suggested that the primary dissipation at the TS is not provided by proton thermal gas, but as was predicted by Zank et al. (1996b) and then studied by , the suprathermal PUIs which preferentially are reflected at the cross-shock electrostatic potential and are capable of getting almost all the heating across the shock.  ). The thermal gas temperature (Richardson et al. 2008) and PUI temperature ) upstream of the TS are about 0.2 × 10 6 K and 1.5 × 10 6 K, respectively, leading PUIs to account for the dominant internal energy. Mostafavi et al. (2017b) and  treated PUIs as a separate component and studied the structure of TS by including both collisionless heat flux and viscosity associated with energetic PUIs in their two-fluid model. Figure 5 represents the normalized PUI, thermal gas, and magnetic pressures to thermal gas pressure far upstream of the TS as a function of distance. It shows that the upstream pressure is almost all converted to PUI pressure, and thermal gas is only heated adiabatically. The thermal gas Mach number, and fast magnetosonic Mach number, are greater than 1 downstream of the TS (Fig. 5), meaning that they are insufficient to lead to a supersonic-subsonic transition across the TS. However, the combined thermal gas and PUI Mach number confirmed that the flow changes from supersonic to subsonic through the TS with the inclusion of PUIs. This illustrates that PUIs play a critical role in the dissipation process and smoothing of the TS (Oka et al. 2011;Yang et al. 2015;Mostafavi et al. 2017b,a;Lembège and Yang 2018).

Shock Crossing by Non-Maxwellian Plasma
In a collisionless shock, the directed flow energy is converted mainly into the thermal energy of plasma species. At the TS a substantial part of energy goes into the heating of PUIs, while heating of the SW protons is weaker than it would be without PUIs. The shock structure width is much less than the PUI diffusion scale, which makes the particle transport and acceleration more complicated than in the isotropic case. In addition, newly accelerated PUIs can generate waves, which in turn affect acceleration and transport of particles in the shock vicinity, typically a few diffusion lengths of energetic PUI (a few au for the majority of accelerated PUIs). One would have to solve an anisotropic particle transport equation based on pitch-angle diffusion. However, solving the focused transport equation (le Roux et al. 2007;le Roux and Webb 2009;Zuo et al. 2011) suitable for this purpose at the entire TS is computationally prohibitive.
While several different approaches have been developed to describe the boundary conditions for PUIs across the TS which describe the energy transfer in this region (Fahr and Chalov 2008;Fahr et al. 2012;Fahr and Siewert 2013b),  use a simplified, but not exactly justified, assumption that thermal protons experience adiabatic compression. However, more accurate approaches are implemented in Gedalin et al. (2020Gedalin et al. ( , 2021a. The extended MHD approach of Zank (2016) provides us with a physics-based width of the shock ramp, but it does not describe the change in the PUI spectrum at shocks. Since the shock transition itself is scatter-free, the fluid approach may miss the physics of relation between the kinetic and MHD scales. Accordingly, the obtained magnetic profiles become monotonic and cannot explain observed overshoots ). In the part of the heliosphere beyond ∼ 30 au the fraction of PUIs becomes substantial (Richardson and Stone 2009), so their contribution to the backstreaming ion beams may become dominant. On the other hand, test-particle simulations of Gedalin et al. (2021a) have shown, that backstreaming PUIs is negligible with the magnetic field angles θ Bn > 55 • if the overshoot is absent.
Previous attempts to take into account the ion kinetics at a shock front assumed magnetic moment conservation (Fahr and Chalov 2008;Fahr et al. 2012;Fahr and Siewert 2013b). However, it has been shown that magnetic moment is, in general, not conserved (Terasawa 1979;Gedalin et al. 2020), and such approximation may be satisfied only in the perpendicular regime. Since the distribution function becomes highly anisotropic in the vicinity of shocks, instead of the application of invalid equations, it may be beneficial to use kineticallyderived boundary conditions at the TS, which is of importance for the interpretation of ENA spectra. In the assumption of co-moving thermal and non-thermal ions, this requires a derivation of boundary conditions for the PUI density and pressure behind a shock, which should be accompanied by an appropriate turbulence model to explain SW ion heating by turbulence. Yet another approach may require us to relate the PUI distribution functions in front of and behind shocks, which may go further and treat the waves produced by unstable distributions of PUIs kinetically.
The proper establishment of the boundary conditions on MHD scales requires establishing a relation of the ultimately isotropic distributions to the collisionless ion dynamics within a shock front (Gedalin et al. 2020). The heating of both thermal species is nonadiabatic. The downstream pressure of the mixture is determined by conservation laws for the whole mixture. SW heating is sensitive to the details of the shock front structure, while heating of PUI is not. This made it possible to parametrize the behavior of PUIs across a quasi-perpendicular shock, such as the TS, using a test particle approach and validate those with full particle-in-cell simulations (Gedalin et al. 2021b). This parametrization requires only two parameters: the angle θ Bn between the magnetic field direction and the shock normal and magnetic compression (the ratio between the magnetic field magnitude downstream and upstream B d /B u ). The direct tracing of PUIs performed by Gedalin et al. (2021b) shows that the fraction of reflected PUIs may be 50% larger than that obtained with the cross-shock potential only. It was concluded that, in agreement with observational data, the cross-shock potential is not the only and probably not the main parameter that determines the intensity of PUI reflection, and gyration of ions in the magnetic field should be taken into account. As a result, the downstream temperature of reflected PUIs may be an order of magnitude lower than that predicted in the cross-shock potential estimates by . It was also shown that the overall heating is substantially stronger than the one predicted by magnetic moment conservation (Fahr and Chalov 2008).
Since the behavior of thermal protons across the TS depends on the details of the collisional shock structure, one needs to resort to pressure balance consideration to estimate the separation of energy between the SW protons and PUIs (Gedalin et al. 2021b). On neglecting the kinetic pressure of the SW and electrons in the upstream region and taking into account that the electron heating is much smaller than the SW proton heating (Schwartz et al. 1988), it has been shown that the relative SW heating is much larger than the relative PUI heating, whereas the situation is opposite for the absolute heating. The downstream SW temperature was shown to be strongly dependent on the fraction of PUIs, while the PUI temperature is practically independent of it. Meanwhile, the mean downstream temperature , varies weakly and remains of the order of 10 6 K.
Test particle simulations need validation with full particle models. Figure 6 shows typical distributions of quantities in front of and behind a collisionless shock (the TS) for different shock angles, obtained with the particle-in-cell code (Gedalin et al. 2021b). For lower angles, the PUI density increase immediately behind the shock ramp is slightly higher than that of the SW protons and that of the magnetic field, but eventually, it approaches the same value at a distance of the order of 100 ion inertial lengths, d i . The large increase in the upstream PUI temperature at θ = 60 • is due to backstreaming PUIs, which are not present at larger angles. The ratio of downstream to upstream T ⊥ increases slightly with shock angle and shock speed. The relative increase of SW proton temperature is much higher than that of the PUIs, although the downstream PUI temperature is still significantly higher than that of SW protons.

Electron Impact Kinetics Downstream of the Solar Wind Termination Shock
Different plasma fluids, like SW protons, pick-up protons, and electrons can be shockprocessed in a fluid-specific way at the SW TS (Fahr and Siewert 2007;Wu et al. 2009;Fahr and Siewert 2011). SW electrons when passing over the shock might experience extraordinary heating (Chalov and Fahr 2013;Fahr and Siewert 2013b;Fahr and Verscharen 2016). Fahr et al. (2015) have derived a theoretical basis for preferential heating of electrons with respect to ions due to conversion of their overshoot kinetic energies into thermal energies in the downstream bulk frame. The highly suprathermal electron distribution function contains a substantial portion of keV energetic electrons that can impact ionize neutral VLISM atoms, like H and He, during its passage through the heliosheath. During the last 50 years, it had been taken for granted that VLISM He atoms, due to the complete absence of ionization sources, cross the heliosheath without any interaction with the heliosheath plasma and get ionized only very close to the Sun (Weller and Meier 1981;Möbius et al. 2004;Chalov 2014). H atoms were thought to undergo charge exchange collisions with the heliosheath protons, but in such a negligible manner that the H density over the upwind heliosheath is practically unaffected.
However, the electrons leaving the TS in a downstream direction are strongly energized appearing in pronounced suprathermal distributions (Fahr et al. 2015). The evolution of the Fig. 6 Characteristic distributions of magnetic field, B, the SW and PUI densities, n SW and n PUI , and temperatures, T SW and T PUI , upstream and downstream a collisionless shock for different angles θ between the magnetic field direction upstream and the shock normal. Since the distribution functions of SW and PUIs are anisotropic, both parallel and perpendicular temperatures are shown. Figure from Gedalin et al. (2021b) SW electron distribution function downstream of the TS can be described under the influence of electron impact ionization of ISN H and He atoms that enter the heliosheath with the interstellar wind from the upwind hemisphere and continue onto the heliotail region. An adequate approach is starting from a kinetic phase-space transport equation in the bulk frame of the heliosheath plasma flow taking into account convective changes, cooling processes, whistler wave-induced energy diffusion, and electron injection into, and removal from velocity-space cells connected with electron impact ionization processes of ISNs.
The phase-space transport equation that describes the evolution of the isotropic electron distribution function f e along heliosheath flow lines after Fahr and Dutta-Roy (2019) takes into account convective changes (first term), the influence of the change in the magnetic field to the change in f e due to the conservation of the magnetic moment of the electrons (i.e., magnetic cooling; second term), velocity space diffusion due to non-linear electron-whistler wave interactions (third term), and electron impact ionization (fourth term) in the following form: Here U is the plasma bulk velocity. The temporal change of the distribution function f e due to electron energy losses is caused by electron impact ionization processes with ambient H or He atoms of the VLISM through the heliosheath. From the kinetic equation, one can ascend to an associated pressure moment of the electron distribution function and arrive at a pressure transport equation describing the evolution of the electron pressure in the bulk-velocity frame of the plasma flow. Assuming that the local electron distribution can be represented by a kappa function with a κ− parameter that varies with the streamline coordinate s, an ordinary differential equation for κ is obtained as a function of s. It turns out that the effect of H and He impact ionizations especially gives its signature to the electron distribution along the extended down-tail region of the heliosheath.
Results show that at the upwind hemisphere of the heliosheath the influences of electron impact ionization are marginal. Downstream it is assumed that the ISN densities are essentially constant. On the downwind side of the heliosheath and down the heliotail this assumption may be questionable. However, this assumption is not bad outside of a circumsolar region of the order of 15 au where the solar ionization reactions produce non-negligible gradients in the neutral atom densities (see Sect. 8.1;Jaeger and Fahr (1998)).

SHIELD Modeling of the Heliosphere
Classically, the heliosphere has long been thought to have a long comet-like tail extending for thousands of astronomical units along the direction of interstellar flow (Parker 1961;Baranov et al. 1971;Baranov and Malama 1993;Pauls et al. 1995;Zank et al. 1996c).  and  suggested this might not be the case and that the solar magnetic field in the heliosheath can collimate the SW plasma into two lobes, leading to a shortened heliotail. They showed that the tension of the twisted solar magnetic field can resist the stretching force of the flowing heliosheath plasma, driving the SW ions into two northward and southward columns. When the flow of the interstellar plasma is taken into consideration, these two columns are bent back into lobes. A few hundred astronomical units beyond the TS in the heliotail, reconnected magnetic field lines are present between the two lobes, and mixed VLISM and SW plasma sit on these field lines, giving rise to a shortened heliotail. While other models (Pogorelov et al. 2015;Alexashov 2015, 2018) have also noted this collimation when using a unipolar solar magnetic field configuration in MHD modeling, these works claim the heliotail to maintain the classical view of a long comet-like structure. Opher et al. (2017) investigated the effect of magnetic reconnection at the HP. The MHD model suppresses reconnection in the nose region while allowing it in the flanks. This is consistent with recent ideas about reconnection suppression from diamagnetic drifts. The stabilizing effect of diamagnetic drifts is an important kinetic effect missed by MHD models, which can develop at boundaries such as the Earth's magnetopause or the HP (Swisdak et al. 2003(Swisdak et al. , 2010. The stabilizing influence of these drifts has been extensively documented with SW and magnetospheric data (Phan et al. 2010(Phan et al. , 2013. The diamagnetic drift suppresses reconnection when the drift speed is larger than the Alfvén speed in the reconnecting magnetic field. The jump in plasma β (the ratio of the plasma pressure to the magnetic pressure) across the nose of the HP is much greater than in the flanks because the heliosheath β is greater at the nose than in the flanks (Opher et al. 2017). Large-scale reconnection is therefore suppressed in the nose but not at the flanks. The effect of this reconnection results in the twisting of the interstellar magnetic field to the azimuthal direction ahead of the HP. Moreover, the reconnection was used to explain the lack of rotation of the magnetic field observed by Voyager 1 (Burlaga and Ness 2014) and Voyager 2 (Burlaga et al. 2019). Kornbleuth et al. (2021b) investigated the effect of reconnection at the HP on the heliosphere and the plasma within. They compared the results of two different kinetic-MHD models from the Solar wind with Hydrogen Ion charge Exchange and Large-Scale Dynamics (SHIELD) model (herein referred to as the BU model) and the Moscow model . Both MHD models treated a single-ion plasma and used a Monte-Carlo approach for modeling ISN H streaming through the heliosphere. Additionally, identical boundary conditions were used for the two models to present a valid comparison. One  Kornbleuth et al. (2021b) key difference between the two models was the magnetic reconnection. The BU model allows for magnetic reconnection at the HP, as it was based upon the MHD model used in  and Opher et al. (2017), while the Moscow model does not. The results of the comparison show that the heliosphere is compressed at high latitudes in the BU model relative to the Moscow model due to the increased magnetic pressure from the twisting of the interstellar magnetic field. The results also indicate that while the BU and Moscow models have different tail configurations beyond 300 au from the Sun, both models displayed a split tail in the sense that the solar magnetic field confines the charged solar particles into distinct north and south columns that become lobes (Fig. 7). The different heliotail configurations arise because the BU model allows for reconnection and MHD instabilities (Opher et al. 2021), which are suppressed in the Moscow model, leading to ISM flowing between the two lobes as in .
Moreover,  separated the SW plasma into two ion fluids: thermal SW ions and PUIs created in the supersonic SW. After crossing the TS, the PUIs would be depleted from charge exchange with neutrals originating from the VLISM, which led to a cooling of the heliosheath plasma. The cooler heliosheath led to a shrinking of the heliosphere such that it became smaller and rounder in shape (Fig. 8).  presented IBEX ENA observations of the heliotail, and found two high latitude lobes of enhanced ENA flux at high energies (∼2 keV), with a deficit of flux in the low latitude tail.  suggested this profile seen at high energies, which was not seen at low energies, could reflect the SW profile from the solar minimum. With fast SW at high latitudes and slow SW at low latitudes in the heliotail, the observed high latitude lobes of ENA flux could reflect the fast SW flowing in the high latitude heliotail.  suggested the confinement of the SW by the solar  . Bottom: The HP as a 3D surface is shown by the yellow surface defined by the SW density of 0.005 cm −3 . The white lines represent the solar magnetic field. The red lines represent the interstellar magnetic field. See more in  magnetic field into two high latitude lobes could be a cause for the IBEX heliotail profile at high energies.  used a time-dependent simulation of the heliosphere with a long tail to model ENA fluxes at IBEX energies. In their ENA model, they partitioned the SW plasma into three distinct populations: thermal SW ions, PUIs transmitted across the TS without reflection, and PUIs reflected and energized at the TS. As per , these three populations were modeled with Maxwellian distributions to approximate the high energy tail of the SW particle distribution. They also included the effect of extinction, which restricts the distance out to which ENAs of particular energy can be observed due to the depletion of parent ions with the same particular energy due to charge exchange. Based on their work,  were able to qualitatively replicate many of the features observed by IBEX, including the high energy tail ENA profile noted by , and concluded that the profile of the SW was the primary driver for the observed ENA structures. Kornbleuth et al. (2018) employed a similar ENA model to  but used a short-tail heliospheric solution corresponding to  with uniform SW boundary conditions to investigate the effect of collimation of the SW plasma by the solar magnetic field on ENA maps. While they were unable to match the energy dependence of the ENA maps as observed by IBEX, Kornbleuth et al. (2018) qualitatively reproduce the IBEX heliotail observations at high energies with two high latitude lobes of enhanced ENA flux and a deficit of ENA flux in the low latitude tail. They concluded that the confinement of the SW plasma by the solar magnetic field leads to an enhancement of ENA flux. Kornbleuth et al. (2020) expanded on the work of Kornbleuth et al. (2018) by investigating four distinct heliospheric solutions with varying SW profiles and solar magnetic field intensities. Two of the solutions used a latitudinally-varying SW profile corresponding to  . Simulated sky maps are multiplied by a factor of 1.8. See more in Kornbleuth et al. (2020) yearly-averaged SW conditions from the year 2008 based on interplanetary scintillation observations (Sokół et al. 2015b), which reflected solar minimum conditions. One of these solutions used a solar magnetic field intensity at 1 au corresponding to that observed by OMNI during the same period, while the other solution neglected the solar magnetic field. When solar minimum conditions were used and the solar magnetic field neglected, Kornbleuth et al. (2020) found a significantly reduced ENA flux in the high latitude as compared with the same solution which included the solar magnetic field (Fig. 9). When the solar magnetic field is included, the modeled ENA maps showed good qualitative agreement with the IBEX observations. Moreover, they found that the ENA flux at high latitudes increased with increasing solar magnetic field intensity. Kornbleuth et al. (2021a) modeled ENA maps using the method of Kornbleuth et al. (2020) from the BU and Moscow MHD models (Kornbleuth et al. 2021b) at IBEX-Hi energies. They found that while the heliotails display differences beyond 300 au, these distances are beyond the cooling length of IBEX-Hi energies. Both MHD models produce a heliotail with heliosheath plasma organized by the solar magnetic field into two distinct high latitude columns, and the two models show qualitative and quantitative agreement with each other in modeled ENA maps. Both models differ quantitatively from IBEX GDF observations from  by approximately a factor of 2. To compare with IBEX GDF observations quantitatively, Kornbleuth et al. (2021a) extracted ratios of the ENA flux for the Voyager 1, port lobe, south pole, and south heliotail lobe directions relative to the downwind direction. The downwind direction was selected as a reference point as ENA flux The spectra show qualitative and quantitative agreement between the BU and Moscow models despite different heliotail shapes, but discrepancies remain with respect to the IBEX data. Adapted from Kornbleuth et al. (2021a) from the downwind direction is only limited by the cooling length and not by heliosheath thickness, for which there is a disagreement between observations and MHD models. The ratios allowed to compare ENA fluxes for different directions without the need to scale ENA fluxes. While there are differences in the ratios between the BU and Moscow models (Fig. 10), Kornbleuth et al. (2021a) found that the differences were largely insufficient to be observed by IBEX considering the estimated ∼ 20% systematic uncertainty of IBEX-Hi observations ). As such, the different heliotail configurations would not be distinguishable at IBEX-Hi energies, and would likely require higher energy observations to probe the shape of the heliotail.

Ionization Rates
Sun drives the environment inside the heliosphere through SW outflow, solar EUV radiation, and resonant radiation pressure (e.g., Bzowski et al. 2013b). This ionizing environment influences particles of interstellar origin inside the heliosphere altering their fluxes (see, e.g., Blum and Fahr 1970;Thomas 1978;Ruciński and Fahr 1989;Sokół 2016). The primary ionization process that affects interstellar particles is charge exchange with SW protons; the next is photoionization by solar EUV flux, and the other is impact ionization by SW electrons. The intensity of the ionization processes varies in time in short time scales depending on active region distribution on the solar surface and in long time scales depending on the cycle of solar activity. Moreover, the ionization rates vary with heliographic latitude, especially for charge exchange, which depends on SW (see, e.g., Sokół et al. 2019a) speed and density which vary in latitude during the solar cycle (SC; see, e.g., McComas et al. (2000), ). These factors need to be accounted for in the model interpretation of measurements.

Fig. 11
Top: time series of SHOIR-calculated ionization rates due to various processes for H, O, Ne, and He in the ecliptic plane at 1 au with Carrington rotation averaged resolution in time for the last three SCs. A total of the rates is presented in red. Bottom: time series of the fraction of the ionization rates due to individual processes to the total ionization rates for a given species. The color code is presented below the bottom panel. See more in  Among the interstellar atoms observed from the vicinity of the Earth's orbit, H and O are the most prone to ionization, next is Ne, and the least affected is He. Figure 11 illustrates ionization rates for these four species in the ecliptic plane at 1 au for the last three SCs (22-24) calculated with the Sun-Heliosphere Observation-based Ionization Rates (SHOIR) model ). The SHOIR model calculates the ionization rates using available SW and solar EUV measurements, including SW variations in latitude, which allows reconstruction of the charge exchange rate variations in latitude after 1985. Figure 11 illustrates charge exchange with SW protons, photoionization, and electron impact ionization rates together with the sum of them all for each species.
Photoionization contributes the most to ionization rates for Ne and He, which follow SC-variations in time with maximum ionization rates during solar maximum and minimum ionization rates during solar minimum. For O, photoionization and charge exchange reactions can contribute equally to the total ionization rates (see Fig. 11). Thus for O, we observed at 1 au variations in time reflecting SC-variations and long-term changes characteristic for the SW temporal variations. For H, the total ionization rates reflect dominantly the SW in-ecliptic variations.
For all species discussed, the SHOIR-calculated total ionization rates show a decrease in time in the last decades, which is due to weakening of the SW flux (e.g., McComas et al. 2008b) and the solar activity, with the last SC (24) weaker than the previous two (22 and 23).
ENAs have been observed inside the heliosphere by several missions (e.g., IBEX, Cassini/INCA). However, before they reach detectors close to the Sun, their fluxes are attenuated by the solar environment, with more attenuation of low-energy ENAs. Thus, to retrieve the ENA flux magnitude in their source regions outside of the TS, a deconvolution of the ionization losses and radiation pressure attenuation is required (e.g., Bzowski 2008). In general, the survival probability of H ENA is calculated by integrating the ionization rate for the ENA exposure to the ionizing factor along its orbit. Calculation of the exposure time in the heliosphere requires solving the equation of motion which takes into account a joint action of solar gravity and resonant radiation pressure by solar Lyman-α photons (e.g. Bzowski et al. (2013b) with reference therein).

ISNs and Interstellar PUIs
The solar environment ionizes ISNs inside the heliosphere producing PUIs (Vasyliunas and Siscoe 1976).  studied the influence of the temporal and spatial variations of ionization rates for the density of ISNs and resulting production rates for PUIs at 1 au as a function of solar activity. The results show that the ISN density and PUI count rates are significantly modulated by variations of ionization rates along the Earth's orbit, which next affects the determination of the ISN flow direction based on the PUI cone and crescent locations. The bias is the greatest for O + PUIs (up to a few degrees depending on the phase of solar activity) and the smallest for He + PUIs. Additionally,  showed that heliolatitudinal anisotropy of ionization rates produces non-negligible asymmetries in the density distribution for O observed in the ecliptic plane in the downwind hemisphere. As further studied by Sokół et al. (2019b), these anisotropies for ISN and PUI densities along the ecliptic plane are present for all species with the greatest effects for O. It means, that the in-ecliptic observations are affected by the latitudinal structure of the ionization rates in the downwind hemisphere, and thus cannot be neglected in the data interpretation.
Distributions of ISN and PUI densities create large-scale structures throughout the heliosphere. Sokół et al. (2019b) studied those global distributions for H, He, Ne, and O for solar minimum and maximum conditions. A common large-scale feature among all species but H is the elongated enhancement of the density that forms behind the Sun in the direction opposite to the gas inflow, i.e. the focusing cone. According to the model results, the cone can reach as far as the TS distance. The maximum of the ISN density is located for various species at various distances from the Sun; the closest to the Sun for He, next for Ne and O, as illustrated in Fig. 12. Also, the location of the ISN density varies with the phase of solar activity being closer to the Sun during solar minimum and further away from the Sun during solar maximum. This variation is an effect of the solar activity variation of the ionization rates. During solar minimum, the ionization rates are smaller and fewer atoms are ionized, and more reach closer distances to the Sun, which results in greater ISN density at closer distances to the Sun. During solar maximum, the ionization rates are greater and more ISNs are ionized which results in smaller density at close distances to the Sun. This also results in an apparent shift of the maximum ISN density to greater distances from the Sun (see the top panel in Fig. 12).
In the upwind direction, the density of ISNs is depleted and drops below e −1 of its value assumed for TS at close distances to the Sun marking the so-called density depletion region (known also as an ionization cavity). The depletion due to ionization spans around the Sun at different distances (closer to the Sun in the upwind direction, and farther away from the Sun in the downwind direction). Also, the cavity location varies between species because of differences in the total ionization rates and a non-negligible radiation pressure for H. The ionization cavity reaches out to ∼ 4 au from the Sun in the upwind direction for H and O, ∼ 1.5 au for Ne, and < 0.4 au for He (Sokół et al. 2019b). In the downwind direction it can reach out as far as 20 au in the case of H (Rucinski and Bzowski 1995;Sokół et al. 2019b).
The production of PUI, which is a product of ISN density and ionization rate, is also modulated by the variations of ionization rates (Rucinski et al. 2003;Sokół et al. , 2019b. The maximum densities of interstellar PUIs are expected at distances from the Sun  Sokół et al. (2019b) greater than 1 au (see Fig. 12). Exceptions are He + PUIs and Ne + PUIs (the latter during solar minimum), for which the greatest densities are expected within 1 au. Similarly, as for the density of ISNs, the maximum of PUIs is at a closer distance to the Sun during solar minimum and at a greater distance during solar maximum.
In addition, the study of heavy atoms from the VLISM provides us with information on the evolution of the interstellar environment close to the Sun. The Milky Way is assumed to be rapidly enriched in Ne and O elements at the time of the formation of the Sun 4.6 Gyr ago. Thus, one expects to find a Ne/O ratio in the contemporary ISM close to the solar ratio. The Ne/O ratio has been also studied through PUI observations  and from X-ray spectroscopy of absorption edges in X-ray binaries (Juett et al. 2006). Moreover, ISN O and Ne have been measured at 1 au by IBEX-Lo instrument Park et al. 2014). However, derivation of the Ne/O ratio in ISM based on those measurements requires model support and interpretation . It is because Ne and O have very different ionization properties (Sect. 8.1) and their abundance ratio is highly variable spatially and temporally inside the heliosphere. While O undergoes charge exchange with protons in the HP, Ne travels practically unimpeded through the interface between heliosphere and VLISM. This difference leads to different kinetic properties of neutrals (and their respective PUIs) near 1 au. Therefore, comparing the ratio of these heavy ions is a tool to disentangle physical processes in the VLISM, in the boundaries of the heliosphere, or inside the heliosphere close to the Earth.

Secondary Population of ISNs
In the outer heliosheath (OHS), the interstellar plasma is heated and decelerated while flowing past the heliosphere. Charge exchange of the primary ISNs with the interstellar ions creates a new, secondary, population of the ISNs measurable in the Earth's orbit. These neutralized ions are decoupled from the very local interstellar magnetic field and then continue to flow away from their birth location. Since charge exchange reactions operate practically without momentum exchange between collision partners, the secondary ISNs inherit the kinematic properties of their parent ions. Thus, the secondary ISNs have slower bulk velocity and higher temperature compared to the primary ISNs. Since the spatial and velocity distribution of the secondary ISNs depends on the plasma properties in the OHS, it is important to measure the secondary ISNs separately from the primary component.
The creation mechanism of the secondary ISNs was anticipated theoretically by Baranov et al. (1981), Baranov and Malama (1993). There are three major species of secondary ISNs: H, He, and O. The existence of secondary ISN H in the heliosphere was confirmed by measurements of the backscattered Lyman-α radiation (Lallement et al. 2005(Lallement et al. , 2010. However, the secondary ISN H has not been unambiguously resolved from the primary population in in-situ measurements because the velocity distribution of the secondary ISNs is quite wide due to its higher temperature, and it overlaps the distribution of the primary ISNs Katushkina et al. 2015). Recently, Katushkina et al. (2021b,a) investigated the possible imprints of the secondary ISN H atoms in the full-sky maps of the H fluxes at the Earth's orbit. They showed that the maps of the H fluxes with energies 1-20 eV for the solar minimum conditions have two specific features: (1) the main peak of the fluxes is split into two blobs, and (2) there are two tails at high ecliptic latitude near the downwind direction. They claim that these features are pronounced due to a signal of the secondary ISN H.
The secondary ISN He was originally believed to be of negligible abundance compared to the primary population (Müller and Zank 2004). It was because the reaction assumed to be responsible for its creation, H e + + H → H e + H + , has a very low cross-section (Barnett et al. 1990). However,  pointed out that the cross-section of the charge exchange reaction between neutral He atoms and He + ions (H e + + H e → H e + H e + ) is comparable to the high cross-section of the charge exchange between H atoms and protons, and interpreted excess He fluxes in the ISN He observations by IBEX as a secondary component (see Fig. 13). Because of the relatively high abundance of He + ions in the VLISM near the heliosphere (Frisch and Slavin 2003), appreciable amounts of the secondary ISN He atoms should be produced in the OHS.  found that the He signal observed by the IBEX-Lo instrument on IBEX is explained as a superposition of the unperturbed ISN He gas and an additional flow component of He atoms, which they called the Warm Breeze (WB).
The existence of secondary ISN O was proposed in the two-shock heliospheric model by Izmodenov et al. (1997). For the secondary ISN O, the most important resonant reaction is the charge exchange between interstellar O + ion and ambient H atom (O + + H → O sec + H + ). The sky maps of the IBEX-Lo O fluxes reveal a tail-like emission expanding toward lower longitude and higher latitude from the observational peak of the ISN O and Ne gas flow Park et al. 2015). The extended tail-like emission could be explained as the imprints of the secondary ISN O population. Baliukin et al. (2017) presented  Figure 13 shows the ISN He and O flux maps from the IBEX-Lo observations and the results of the analytical model by . The emissions in the longitude of 150-210 • result from the secondary ISN populations, which appear as the extended tail-like emission.
One of the key points in the models of the secondary ISNs is the velocity distribution function at the outer boundary. The simplest case assumes that the velocity distribution functions of both the primary and secondary ISNs are non-interacting homogeneous Maxwellian distribution functions. This two-Maxwellian approximation has been used by several previous studies of the IBEX ISN He and O observations . Since the interstellar ions are decelerated, heated, and deflected in the OHS, this simple approximation does not represent the secondary ISN in the OHS. However, it is a reasonable approximation inside the TS because only a small subset of the whole secondary ISN population can penetrate inside the TS (e.g., ). On the other hand, the kinetic-MHD model of the SW interaction with the VLISM has been used to determine the distribution of ISN atoms at the outer boundary. Using this global model, Baliukin et al. (2017) calculated the distribution of the combined O atoms (primary and secondary ISN O) as a Maxwellian-like core with a wide plateau. The wide plateau consists of slow atoms (slower than the bulk speed of ISN gas flow, i.e., 26 km s −1 ) and it is a product of charge exchange in the OHS.

PUIs Close to the Solar Corona
Interstellar PUIs are omnipresent plasma phenomena connected with the presence of charge exchange reactions between plasma flows and neutral atoms. Fahr et al. ( , 2000a have shown that those PUIs also appear in regions close to the solar corona, and they can be seen there in their EUV resonance glow radiations at 304 Å as discussed by Paresce et al. (1981) for He + PUIs. ISN He atoms pass over the heliosphere and approach the Sun to within less than 10 solar radii before they get ionized. Next, the He + ions are picked up by the SW plasma flow and are transported out to larger solar distances. The background for this phenomenon was given in early papers by Fahr (1968a,b).
While ISN H atoms are strongly reduced in density by charge exchange processes with SW protons inside ∼3 au, ISN He atoms can penetrate much deeper towards the Sun and finally are extinguished and converted into He + ions (see Sect. 8.2). Initially, these ions have their original Keplerian velocities, but as ions, they are picked up by the ambient magnetic fields and are incorporated into the SW plasma flow. However, the pick-up process inside of 10 solar radii is very different from how it works at distances of several astronomical units. The main differences are because of the magnetic fields which are nearly unwound, but rather radial, and the magnitude of the B field dominates the kinetic energy of the ions.
These HeII ions are resonantly scattered by the solar EUV radiation at λ 0 = 303.8Å where the solar radiation has a very distinct and strong emission line; however, with a narrow spectral width of only λ 0 = (u c /c)λ 0 = 0.1Å. Thus, the strength of the resonantly scattered radiation by the He + PUI is strongly dependent on its Doppler-shifting ion velocity u pui . The resonantly scattered radiation at the solar He + resonance line of 303.8Å, taking it to be a Gaussian line emission, is given by the following expression: The calculated intensity of this resonance radiation seen from Earth along a line of sight perpendicular to the upwind-downwind line of the LISM at different solar impact distances in Rayleighs is given in Fig. 14.

Terrestrial ENAs
The shocked SW protons after the passage of the Earth's bowshock could act as a diffuse source of H ENAs. However, it requires the terrestrial H-exosphere to extend up to magne-tosheath distances of 8 to 10 Earth radii where the shocked SW plasma is present. Recently, Baliukin et al. (2019) derived from SOHO/SWAN Lyman-α data that the H geocorona extends well beyond the Earth's Moon. The existence of terrestrial H ENAs from the region ahead of the Earth's magnetopause has also been studied by measurements carried out by the IBEX satellite (Fuselier et al. 2010, their Fig. 1).  quantitatively calculate the H ENAs expected from the exosphere of the Earth using the Lyman-α studies from the TWINS/LAD detector on the H density of the terrestrial exosphere (see Zoennchen et al. 2011;Zoennchen et al. 2013). They modeled the magnetosheath plasma flow ahead of the solar part of the Earth's magnetopause and determined the source flow of H ENAs originating from the charge exchange collisions of terrestrial H atoms with magnetosheath protons. As a result, they obtained a quasi-isotropic source of H ENAs with keV-energies that partly bombard the exobase, elastically collide with the most abundant O atoms, and with a probability of 93% are reflected from the exobase with typical energies ofẼ H ENA (16/17) keV. With such super-escape energies, these reflected ENAs leave the exobase to outer space along straight trajectories, and their local densities fall off like n H ENA ∼ (r ex /r) 2 . Since the density of the normal exobase H atoms at the upper border of the exosphere falls off like n H (r) = n H 0 · (r/r 0 ) −2.7 , the relative abundance of H ENAs over H atoms permanently increases with height. This reveals that the (H ENAs/H atoms) cross over point is dependent on the location of the magnetopause r MP (it is further into the Earth for a low, and further out for a high magnetopause location; see e.g. Baumjohann and Treumann (1996)).

Inner Source PUIs
Inner source PUIs are a population of PUIs produced close to the Sun. Three main features differentiate them from interstellar PUIs. First, observations of C + PUIs -a species mainly ionized in the interstellar medium due to its high ionization rate . Second, observations of the radial profile of inner source C + decreasing with radial distance from the Sun . Whereas the radial profile of interstellar species such as O + has been shown to increase with radial distance (see, e.g., Fig. 12). Third, the velocity distribution function of inner source C + has a thermalized distribution that peaks near the SW speed (Schwadron and Gloeckler 2007). This indicates that they are picked up close to the Sun and cooled with the expanding SW. Compared to interstellar PUIs, inner source PUIs are still poorly understood.
The composition of inner source PUIs resembles the composition of the SW (Gloeckler et al. 2000;Allegrini et al. 2005;Schwadron and Gloeckler 2007). However, being single or double charged excludes them from being produced directly from the Sun. Instead, SW ions reduce their ionization through the interaction with interplanetary dust particles (IDPs). IDPs are dust grains from sources such as the asteroid belt, Edgeworth-Kuiper belt, comets, planets, and the VLISM. They are distributed throughout the heliosphere by a combination of fragmentation, electromagnetic forces, gravity, radiation pressure, and SW pressure.
The production of inner source PUIs is highly dependent on the composition, size, velocity, and density of the IDPs they interact with. The size distribution of the IDP population increases with smaller grain size due to the amount of dust-dust collisions (Ishimoto 2000). There is a large deficiency of IDPs of about 10 −15 g to 10 −12 g (i.e., beta-meteoroids) due to following hyperbolic orbits because of radiation pressure. However, grains smaller than beta-meteoroids (i.e., nanograins) have stable orbits close to the Sun (Mann et al. 2007).
There are five proposed mechanisms for inner source PUI production due to the interaction of SW ions with IDPs: SW neutralization (Wimmer-Schweingruber and Bochsler 2003), SW recycling (Schwadron et al. 1999), backscattering, sputtering, and sputteringinduced recycling (Quinn et al. 2018). In neutralization, SW ions gain electrons as they pass through nanograins and the edges of larger IDPs. In recycling, SW ions are implanted within larger IDPs and neutralized in the process. The surface of the IDP becomes saturated and the atoms eventually desorb from the grain, then become ionized and picked up by the interplanetary magnetic field (IMF). In backscattering, charge exchange occurs when SW ions reflect off the surface of IDPs. In sputtering, bombarding SW ions remove atoms that compose the IDP. Lastly, in sputtering-induced recycling, bombarding SW ions instead remove SW ions that were previously implanted into the IDP.
Which production mechanism is dominant is still uncertain. Observations of inner source C + and O + show that C + /O + > 1 and decreases with SW speed (Taut et al. 2015). Simulations of neutralization, backscattering, and sputtering-induced recycling tend to reflect these observations (Quinn et al. 2019). However, observations of the radial profile, latitudinal profile, size distribution, composition, and morphology of IDPs close to the Sun are necessary to gain more insight.  (Gloeckler et al. 1992) within the H ionization cavity, from about 1.4 to 5.4 au, identified numerous PUI species (H + , He + , N + , O + , Ne + ) (e.g., Gloeckler and Geiss (1998)). The distributions tended to be highly anisotropic , suggesting that the PUI scattering mfp was unusually long, although the precise reason for this is somewhat unclear. It is possible that the nature of the underlying turbulence responsible for the scattering of PUIs, being a superposition of a majority quasi-2D and a minority slab component, might be responsible for the long mfps.
The initial ring-beam PUI distribution is unstable to the excitation of Alfvénic modes (Lee and Ip 1987;Williams and Zank 1994;Zank 1999;Cannon et al. 2014b;Aggarwal et al. 2016;Hollick et al. 2018a), and therefore modifies the underlying low-frequency turbulence field responsible for scattering PUIs. The excitation of Alfvénic fluctuations by PUIs is key to understanding the observed radial evolution of magnetic turbulence in the distant SW (Zank et al. 1996a;Matthaeus et al. 1999;Breech et al. 2008;Oughton et al. 2011;Adhikari et al. 2014Adhikari et al. , 2015Wiengarten et al. 2016;. Furthermore the dissipation of low-frequency magnetic fluctuations heats the SW plasma and results in a non-adiabatic SW temperature profile and a slow temperature increase beyond ∼ 30 au . PUI generation of low-frequency turbulence in the outer heliosphere also affects cosmic ray mfps (Zank et al. 1998;Florinski et al. 2003;Engelbrecht and Burger 2013;Chhiber et al. 2017;Engelbrecht 2017;).
There is an intimate coupling of PUIs, thermal SW plasma, and the evolution of lowfrequency turbulence throughout the SW. The first PUI-mediated models of the SW assumed instantaneous assimilation and equilibration of newly created PUIs with the SW plasma (Wallis 1971;Holzer 1972;Holzer and Leer 1973;Isenberg et al. 1985) i.e., the combined thermal SW and PUI plasma was treated as a single Maxwellian (or at least isotropic) distribution. The early models predicted a dramatic SW temperature increase with the increasing heliocentric distance that was not observed. Isenberg (1986) pointed out that the SW plasma and PUIs could not form a single equilibrated distribution because the Coulomb collision time scale between suprathermal pickup protons and SW protons (and electrons) far exceeds the fluid advection time for a heliospheric supersonic flow of radial extent less than 100 au (i.e., the distance to the heliospheric TS). Instead the total distribution comprises a cold thermal SW core and a tenuous energetic pickup proton halo (Isenberg 1986;Zank et al. 1996b;Zank 1999;Zank 2015). As discussed by Isenberg (1986), , Zank (2016), a PUI-mediated model of the supersonic SW must include the thermal SW and the suprathermal pickup protons as distinct plasma components.  extended the models of Holzer (1972) and Isenberg (1986) by coupling a detailed model of turbulence transport and dissipation to a multi-fluid description of the SW plasma. This allowed them to examine the feedback between SW plasma heating, the modified large-scale SW velocity due to the creation of PUIs, and the driving of turbulence by SW and interstellar PUI sources.
The derivation of the equations for a PUI-mediated SW, including the structure of the TS, presented by , Zank (2016),  assumes that preexisting and self-generated Alfvén waves scatter PUIs and a form of collisionless Chapman-Enskog expansion is then used to derive model equations describing the evolution of the PUI density ρ p , velocity, and pressure P p . Since the scattered PUI distribution is nearly isotropic, the fluid description for PUIs contains both a heat conduction term ∝ K p · ∇P p , where K p is the collisionless PUI "thermal conduction" tensor and a stress tensor due to the collisionless viscosity η ij . The heat conduction term is the first-order anisotropy and the viscosity is the second-order anisotropy in the PUI distribution. Although the scale lengths introduced by both PUI heat conduction and viscosity are much smaller than the large-scale supersonic SW, where they can safely be neglected, they are important in determining the structure of the TS (Mostafavi et al. 2017b;. The PUI-mediated model of the supersonic SW is completed with the inclusion of turbulence, both pre-existing and excited by the creation of PUIs. The coupled model equations can be found in . An important distinction between earlier PUI-mediated SW models and the  model is in the 1D thermal plasma temperature transport equation, expressed as Adhikari et al. (2015) U dT s dr Here r is the radial distance, T s is SW temperature, γ s = 5/3 is the adiabatic index, n s is the SW density, k Boltzmann's constant, U is the large-scale radial flow velocity, and the turbulent heating term S t for the thermal plasma and the charge exchange-related term ν s c is given by Equations (31) and (12), respectively, in . The heating term S t expresses the dissipation of turbulent energy as measured in the quasi-2D and slab components ).  solve the PUI-mediated SW model numerically and restrict attention to a simple spherically symmetric steady-state model. This describes the evolution of the bulk SW and the associated turbulence from 1-75 au. The predicted steady-state solutions are compared to both Voyager 2 plasma and magnetometer data and NH plasma data. The largescale SW and PUI plasma number densities, pressures, and corresponding temperature, as well as the SW radial component of the flow velocity, are plotted in Fig. 15 as a function of heliocentric distance r. Each panel shows the theoretically predicted results together with Voyager 2 data from 1-75 au and NH SWAP data from 11.26-38 au. The two data sets are independent, being separated by three decades, but we show the spatial region of overlap. The same theoretical model curves are used in all plots. The Voyager 2 observations are taken from 1983-1992 and those of NH from 2008-2017 for the radial heliocentric distance interval 11-38 au, and the SC observed by NH was much weaker than that observed over this distance interval by Voyager 2 (e.g., Lockwood et al. 2011).
The SW density decays as usual according to n s ∝ r −2 , whereas n p increases to a maximum value at ∼ 12 au and then very slowly decreases as the SW expands radially outward. The model and Voyager 2 and NH SWAP observations exhibit a good correspondence, although some SC dependence is present in the thermal plasma density. This is more apparent in the Voyager 2 radial velocity plot. The PUI density observations exhibit shorter time scale variation, most likely associated with interplanetary shocks or streams.
The theoretical thermal and PUI pressures and temperatures plotted in Fig. 15 exhibit a change in the cooling rate at ∼ 20 au, which has been attributed to the additional dissipation of turbulence driven by PUIs. The effect of heating by turbulence is quite striking as can be seen from the dashed (no dissipative heating) and solid line (dissipative heating included) curves plotted in the pressure and temperature panels of Fig. 15. The agreement between the theoretical and observed Voyager 2 pressure and temperature is very good with the slow flattening/increase in the thermal pressure/temperature profile. A corresponding flattening/increase in the NH SWAP thermal plasma data appears to be consistent with predictions but more data is needed to confirm this. The PUI pressure and temperature profiles are consistent with NH SWAP observations, the PUI temperature being more than two orders of magnitude higher than that of the thermal pressure.
An interesting assumption about the temperature plots in the  model is that the dissipation of PUI-driven turbulence occurs only in the thermal plasma i.e., the dissipation of low-frequency turbulence does not heat the PUIs. However, the PUI temperature and pressure plots do indicate a turn-up in the PUI temperature near 38 au, and this is not consistent with the model presented by . McComas et al. (2021) confirm that the PUIs do appear to increase in temperature with increasing heliocentric distance. The reason for the heating is not entirely obvious. There are several possible explanations for this. The first is that we do not understand the nature of turbulent heating in a nonequilibrated plasma system, and it may be entirely possible that the hotter PUI distribution may experience some heating despite their already high temperature. A second, and perhaps more plausible, the explanation may be that the PUIs are further energized by some form of second-order Fermi acceleration related to counter-propagating Alfvén waves (Bogdan et al. 1991;Le Roux and Ptuskin 1998). A third possibility is that some form of further acceleration of PUIs related to the preferential reflection and subsequent energization of PUIs by interplanetary shocks (Zank et al. 1996b), as has been observed by the SWAP instrument (Zirnstein et al. 2018c), may account for the observed increase in the PUI temperature.
The SW radial velocity plot is highly variable compared to the steady-state solutions in Fig. 15. This is due to both SC variation in flow speeds, as well as the presence of shock waves and stream structure. The theoretical radial velocity plots exhibit a gradual deceleration of the SW due to charge exchange (Holzer 1972;Wallis 1971;Isenberg et al. 1985;Fig. 15 Comparison of a theoretical PUI-mediated SW model and observed bulk plasma observations for Voyager 2 data between 1 and 75 au (top plots of each panel) and NH SWAP data between 11.26 and 38 au. Top left: Solid lines show model predictions of the large-scale thermal plasma number density n s (black curve) and the PUI number density n p (red curve) together with Voyager 2 observations of n s (blue symbols, top plot) and NH SWAP observations of n s (green) and n p (red) (bottom plot). Top right: Same as previously, but for thermal plasma and PUI pressure for Voyager 2(top) and NH SWAP (bottom). The dotted curves show the pressure without turbulent heating. Bottom left: Thermal temperature (black line, blue symbols) and PUIs (red line). The dotted curves show the temperature without turbulent heating. Bottom right: Predicted and observed SW radial speed U . See more in   Khabibrakhmanov et al. 1996;Zank 1999), which is not easily observed, although see the analysis by Richardson et al. (2004).

Waves Due to Interstellar PUIs
ISNs which cross the heliosphere and get ionized represent a suprathermal energy source moving through the SW along the mean-field at the SW speed. They resonantly excite sunward propagating, right-hand polarized waves which become left-hand polarized in the spacecraft frame (Lee and Ip 1987). As the waves propagate at the Alfvén speed and the ions are scattered in pitch angle, the waves are seen at spacecraft-frame frequencies f sc ≥ f i,c where f i,c is the ion cyclotron frequency. While the process of scattering interstellar PUIs (top to bottom) The trace of the power spectral density matrix and spectrum of the magnitude time series (indicative of low compression) with the cyclotron frequency of He + and H + shown, degree of polarization and coherence, ellipticity showing the traditional (signed) polarization of the fluctuations, and the angle between the minimum variance direction and the mean magnetic field. The waves are seen at f sc > f i,c with a high degree of polarization and coherence, are left-hand polarized in the spacecraft frame, and have minimum variance directions parallel to the mean magnetic field. See more in Hollick et al. (2018a) excites magnetic wave energy that drives the turbulence, it is not possible to observe that wave energy when it is absorbed by the turbulent cascade. When the turbulent cascade is sufficiently weak to allow the wave energy to accumulate to observable levels over many hours, the waves can be observed in the magnetic field data.
Waves due to newborn interstellar PUIs have been observed by magnetometers on numerous spacecraft including ACE (Argall et al. 2015;, Ulysses (Murphy et al. 1995;Cannon et al. 2013Cannon et al. , 2014bAggarwal et al. 2016;, and Voyager (Joyce et al. 2010(Joyce et al. , 2012Hollick et al. 2018a,b,c) and are reviewed by . Figure 16 shows an example of Voyager 1 spectra exhibiting waves due to newborn interstellar H + and He + PUIs when the spacecraft was at 6.2 au. The waves due to H + offer the stronger signal, but the ellipticity at f sc > f H e,c does show the expected polarization signature. Both PUI species produce waves with E lip < 0 (left-hand polarized in the spacecraft frame) as expected from theory (Lee and Ip 1987). While waves due to interstellar He + PUIs can be seen at 1 au, interstellar neutral H does not penetrate inside 3 au in sufficient number to excite waves to an observable level. Waves due to newborn interstellar H + PUIs are seen in the Ulysses and Voyager data beyond 3 au.
The computed polarization in the spacecraft frame for wave intervals arising from H + PUIs in the Voyager 2 data (Hollick et al. 2018a) is shown in Fig. 17. The results are ob- Fig. 17 Spectral ellipticity averaged over the frequency range f p,c < f sc < 2f p,c as a function of the heliocentric distance. See more in Hollick et al. (2018a) tained by averaging the ellipticity spectrum over the frequency range f p,c < f sc < 2f p,c where f p,c is the proton cyclotron frequency. Approximately, half of the wave events are right-hand polarized beginning at 5 au which is contrary to the standard theory and most like attributable to incomplete or modified scattering of the PUI. Approximately 10% of the waves due to He + PUI seen by the ACE spacecraft and those excited by H + seen by the Ulysses spacecraft are found to be right-hand polarized in the spacecraft frame. The observation of right-hand polarized waves becomes more prominent beyond 5 au. This is thought to be due to incomplete or modified scattering of the PUI, but this has not been proven. The waves observed inside ∼ 3 au are unlikely to be due to interstellar PUIs and are thought to originate either from an inner source such as dust grains (see also Sect. 8.6) or ENAs. That study is ongoing.
The results of analysis of the ellipticity of waves due to H + PUI as seen by the Ulysses spacecraft (top) and Voyager 2 spacecraft (bottom) is combined in Fig. 18. The ellipticity is plotted as a function of the angle between the mean magnetic field and the radial direction, BR . Waves of both polarizations are consistently seen over the full range of field orientations, although the Ulysses results are a bit more limited in extent showing few right-handed waves at ellipticities < 25 • and > 75 • . It has been speculated that the waves are most often seen when the mean field is radial, owing in part to an approximation made in Lee and Ip (1987). Figure 18 shows that this is not the case, although right-hand polarized waves are not generally seen for nearly radial mean fields. The source of right-hand polarized waves is not yet determined.

Acceleration and Global Energy Distribution of PUIs
Although PUIs are born to have their thermal speeds nearly equal to the SW speed, observations by experiments on various spacecraft, e.g., ACE, WIND, Ulysses, Cassini, and NH, have shown that the speed of PUIs can go much higher even in periods of low SW disturbance. Quite often, an energetic tail of PUIs with a power-law distribution function close to ∝ p −5 can be observed (e.g., Zhang et al. 2019). PUIs are accelerated in the interplanetary space and the mechanisms for the acceleration are still under active discussion. Interplanetary shocks are one of the possible explanation (Giacalone et al. 1997;Summerlin and Baring 2006;Giacalone 2010). However, only a few interplanetary shocks were observed, and thus cannot be used to explain the more prevalent acceleration of PUIs in the quiet SW. The small-scale plasma waves or turbulence, which is ubiquitous in the SW, are another source to accelerate PUIs (Isenberg 1987;Schwadron et al. 1996;Le Roux and Ptuskin 1998;Gamayunov et al. 2015).
Fluctuations of electric fields in plasma turbulence appear to be random. If its amplitude is small enough, the effect of fluctuating fields leads to particle diffusion in momentum. It is called the second-order stochastic acceleration, similarly to the mechanism proposed Fig. 18 Computed wave ellipticity as a function of the angle between the mean field and the radial direction. (top) For waves due to interstellar H + PUIs seen by the Ulysses spacecraft (Cannon et al. 2014b). Only about 10% of the observations are right-hand polarized as is the case for waves due to He + seen by ACE . (bottom) Same as top but for Voyager 2 data (Hollick et al. 2018a) initially by Fermi (1949). The mechanism can accelerate particles to higher energies as well as cause particles to lose energy. Because we have more low-energy particles in typical space plasma or the gradient of particle distribution as a function of particle momentum is negative, the mechanism can cause a flow of particles towards higher energies, thus making the system gain energies statistically.
In the SW plasma, most of the plasma fluctuations are in the form of incompressible Alfvénic turbulence. The momentum diffusion coefficient D pp is related to particle pitch angle diffusion and parallel spatial diffusion coefficient κ through where V A is Alfvén speed, and the constant A is of the order 1/9.  and Gamayunov et al. (2013) considered plasma waves generation by PUIs and acceleration of PUIs self-consistently. They found that the acceleration of PUIs by Alfvén fluctuations is not fast enough to overcome the adiabatic cooling effect in the expanding SW, thus producing an insignificant number of particles above the SW speed. Zhang and Lee (2013) pointed out that the particle's parallel spatial diffusion coefficient κ is so large for PUIs that the momentum diffusion time p 2 /D pp is longer than adiabatic cooling time in the SW 3/ ·V, where V is the SW velocity. The mechanism, however, can heat up the core SW somewhat presumably due to a smaller particle diffusion coefficient at these low energies. Stochastic acceleration of particles can also occur when there are plasma motions that exhibit random compressions or expansions. If the particles can cause diffusion relative to the plasma, some particles may gain substantial energies by repeatedly going through plasma compression regions. The mechanism was investigated by Toptygin (1982, 1993) and Ptuskin (1988), who obtain a momentum diffusion equation using a quasilinear approximation. Ptuskin (1988) dismissed it as an effective mechanism to accelerate or reaccelerate cosmic rays in the interstellar medium based on an estimate of acceleration time using a typical cosmic ray spatial diffusion coefficient. Using an alternative term of "pumping mechanism," Gloeckler (2006, 2008) and Fisk et al. (2010) suggested that it can automatically produce the p −5 distribution of particles commonly observed in the heliosphere. But the quasilinear calculation by Jokipii and Lee (2010) disputes their claims. Zhang (2010Zhang ( , 2011 simulated particles going through trains of compressible plasma waves by solving the Parker transport equation. He found that particle acceleration by this mechanism can be fast and efficient even with amplitudes of plasma compression waves. It is particularly true for those particles with low energies. Significant particle acceleration occurs when the particle diffusion length is greater than the size of the individual plasma compression region but smaller than the entire plasma fluctuation region. In some of the simulations where shock wavelets are prescribed, it is found that particles need to go through as low as a few shocklets with small compression of a few percent to get accelerated a steady-state solution of a power-law distribution function. Compressible plasma fluctuations tend to propagate in the direction perpendicular to the magnetic field. The spatial diffusion coefficient should be much smaller than the parallel one. That is the main reason why acceleration by compressible plasma fluctuation is much faster than that by Alfvénic fluctuations. If the fluctuation is left to develop fully into many shocklets at compression fronts, the wave power spectrum becomes S(k) ∝ k 2 in a certain range of k between k min = 1/L 0 and k max = 1/L 1 .
The momentum diffusion rate D 0 is independent of particle momentum. With this momentum dependence of D pp , the momentum diffusion equation has a steady-state solution of the distribution function for accelerated particles of the form ∝ p −3 . Such a distribution, if left unmodified, contains an infinity amount of energy density and particle density. In reality, because particles suffer adiabatic energy loss and other possible loss mechanisms, the accelerated particles can never achieve this distribution over a wide range of particle momentum. As a result of acceleration, the pressure contained in accelerated particles can decrease the amplitude and wavenumber of compressible plasma waves, which in turn limits the momentum diffusion rate D 0 . Zhang (2010) and Zhang and Lee (2013) showed that the momentum diffusion and adiabatic cooling settle at some equilibrium rates so that the pressure of accelerated particles approaches a steady maximum value. At that point, a ∼ p −5 power-law distribution function is established without details of wave or particle parameters.
Based on the above mechanism of particle acceleration by compressible plasma fluctuations, Zhang and Schlickeiser (2012) built a bi-modal model for PUI acceleration throughout the heliosphere. The model considers that intense compressible plasma fluctuations only exist in a limited percentage of volume in the heliosphere. At the same time, the accelerated particles can populate other regions of the heliosphere without the presence of compressible plasma fluctuation through some kind of particle transport process or due to propagation of the regions. Most likely, regions of particle acceleration by compressible plasma fluctuations can form and disappear frequently. It is because compressible plasma fluctuations can rapidly convert energies to the internal energy of accelerated particles. That requires the generation of plasma compressible plasma fluctuations regularly through shear flows or conversion from Alfvenic turbulence. In Zhang and Schlickeiser (2012) model of PUIs acceleration, particles also experience various loss mechanisms.
With the condition that particles in particle acceleration regions can only produce a limit pressure in PUIs, Zhang and Schlickeiser (2012) established a maximum stochastic acceleration rate without an explosive increase of particle pressure. They showed under the maximum acceleration rate how particles could automatically establish a ∼ p −5 distribution function in suprathermal energy tails. Figure 19 shows the calculation results of the omnidirectional PUI distribution function at 1 au, 5 au, and near the TS, compared with observations by ACE and Ulysses. Figure 20 shows radial dependence of PUI density and pressure calculated from the model and comparison to those expected radial variations of the SW. The slow down of the SW due to the mass loading of PUIs is neglected. As can be deduced Fig. 19 The omnidirectional distribution function of PUI protons in the spacecraft reference frame in the SW at 1, 5, and 80 au and their comparison with measurements by ACE and Ulysses. Zone 2 is the region where active stochastic particle acceleration occurs in the presence of compressible plasma fluctuations. Zone 1 does not accelerate particles, but it receives accelerated particles through leakage from Zone. The calculation assumes that stochastic particle acceleration in Zone 2 is so intense that the back-reactions of particle pressure can moderate the amplitude of plasma fluctuation and hence particle acceleration rate. The velocity of the slow SW is 400 km s −1 , and the fast SW has 800 km s −1 . Figure adapted from Zhang and Schlickeiser (2012) from Fig. 20, the ratio of PUI to SW density and pressure increase with radial distances. By radial distances around 10 to 20 au, the pressure in PUIs dominates all magnetic field pressure, solar thermal pressure, and turbulence pressure all combined. It is only less than SW ram pressure. By the time it reaches the TS, the Mach number of the entire composite, which is equal to the square root of the ratio of SW ram pressure to PUI+SW internal pressure, is only around 3-4 at ∼ 80 au even with unaccelerated PUI pressure (Fig. 21). As a result, the TS is no longer a supercritical shock. Particle acceleration of PUIs during their propagation Fig. 20 Radial dependence of PUI protons, the SW, ISN H, and magnetic field parameters. The Parker model is used in the calculation of SW and magnetic field parameters up to the TS. The rate of PUI injection is calculated from the charge exchange between interstellar neutral and SW ions. PUI proton pressure is calculated without stochastic acceleration. The core SW (thermal) pressure is from adiabatic calculation assuming no interplanetary heating. The 2D turbulence pressure is assumed to follow the WKB approximation. The sound wave speed is calculated from the sum of core SW pressure and pickup proton pressure. Figure   ENAs' energy spectrum at 1 au in the Voyager 1 direction based on IBEX-Hi data and model results without and with shock presence by Mostafavi et al. (2019). The gray region represents an estimate of the ENA flux during solar minimum. See more in Mostafavi et al. (2019) through the interplanetary space partially contributes to PUI pressure buildup, thus further reducing the SW Mach number to below 2 in the quiet SW (Zone 1) at ∼ 80 au. In regions of active particle acceleration (Zone 2), the Mach number could dip below 1 outside 60 au. These numbers suggest that the TS is rather weak.

Shocks in the Inner Heliosheath
Nonthermal PUIs in the IHS introduce a large collisionless heat flux and viscosity by their pitch-angle scattering, and thus their contribution cannot be ignored . The theory estimated that the PUI number density is about 25% of the thermal gas density in the IHS , which leads to a significant PUI pressure (i.e., about 25 times than the thermal gas pressure; ). The comparison between the ENA spectra computed from  model with the IBEX observations showed that transmitted and reflected PUI downstream of the TS contribute to the ENA flux with 0.5-5 keV, and above 5 keV energy ranges, respectively (Desai et al. 2012). Later,  showed a significant increase in the ENA fluxes with energies below 0.5 keV by including the ENA production from the charge exchange between an ISN atom with a PUI outside the HP.  improved the multi-component treatment of the IHS plasma by accounting for the loss of PUIs through charge exchange as they propagate in the IHS. Their results suggested that the significant number of ENAs between 0.02 to 10 keV originate from outside the HP. However, there is a large discrepancy between these models and IBEX ENA observations, suggesting that they miss some PUI energization in the IHS. Zirnstein et al. (2018a,b) showed that including turbulence in IHS could increase the PUI energy and account for some part of the missing ENA flux.
Interplanetary shock waves generated by coronal mass ejections or corotating interaction regions (Pizzo 1978;Gosling and Pizzo 1999) propagate in the heliosphere, collide with the TS, and eventually produce one or more weakened transmitted shocks propagating in the IHS (Story and Zank 1997;Washimi et al. 2011;Provornikova et al. 2012). Any shock waves propagating in the IHS can be a significant source of PUI energization, and one should not ignore their presence in analyzing IBEX observations. Mostafavi et al. (2019) used a two-fluid (i.e., thermal gas and nonthermal PUIs) model presented by  to study the structure of IHS shock waves and their effect on the ENA flux. As at the TS, shocks in the IHS are mediated by energetic PUIs, and the PUI diffusion length scale determines their thickness. They assumed two shocks propagating in the IHS to calculate the ENA flux. The simulation results were compared with the first five-year average of IBEX observations  in the Voyager 1 direction. Figure 22 compares the shock model results and IBEX observations at 1 au. The presence of transmitted interplanetary shocks in the IHS increases the IHS PUI temperature primarily (Mostafavi et al. 2019). Thus, it leads to the more effective production of ENAs from the IHS. The extra production of the ENA flux due to the IHS shocks yields a better consistency with the ENA flux observed by IBEX at 1 au.

Summary
The Voyager 1 and Voyager 2 missions opened a window to characterize the outer solar system environment and search for the heliosphere boundary, the outer limits of the Sun's magnetic field, and the outward flow of the SW. IBEX was the next step to extensively explore the interstellar space by remote observations. And NH added significantly to the knowledge about the PUI mediation of the SW inside the heliosphere. They all greatly contributed to our understanding of the heliosphere and the development of the physics of PUIs and ENAs. However, still many questions remain open, especially when theory and models meet the measurements, as discussed here. More observational data is needed to advance and guide the theoretical picture of the heliosphere, its interior, and the processes at its boundary offered by the present-day models. The last decades show again that for such a complex system like the heliosphere modeling and measurement must coexist to satisfy progress in both.
NASA is preparing to launch the Interstellar Mapping and Acceleration Probe (IMAP) in 2025 , which goal is to better understand the boundary of the heliosphere, the role of PUIs, and the physics of ENAs. The concept of an interstellar probe that would leave our solar system and could make in-situ measurements in interstellar space would be the next huge and needed step in exploring our stellar environment. Several missions have been considered, like, e.g., the Interstellar Express 2 (Wu et al. 2019;Song 2021), the Interstellar Probe 3 , and the Interstellar Heliopause Probe 4 (Lyngvi et al. 2005). The coming decades with the new and innovative studies will bring answers to many questions and, for sure, will reveal new unknowns.

Conflict of Interest
The authors declare that they have no conflict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. 1