INTERSTELLAR PICKUP ION PRODUCTION IN THE GLOBAL HELIOSPHERE AND HELIOSHEATH

A neutral atom from the local interstellar medium (LISM) can be ionized while drifting into the heliosphere and turn into an interstellar pickup ion (PUI). The neutral atom can become ionized through charge exchange with an ion, photoionization by sunlight or impact ionization by solar wind (SW) electrons. Charge exchange is the dominant mechanism of the three. In the SW frame, a newly created ion is "picked up" and starts to gyrate about the direction of the magnetic field. Subsequently, they may be scattered into an isotropic distribution by either ambient preexisting or self-excited waves (Vasyliunas & Siscoe 1976; Isenberg 1987; Zank 1999). The resulting PUI velocity distribution function (VDF) in the reference frame of the SW is a spherical shell, at the lowest order. These PUIs are also convected outward by the expanding SW. As the PUIs travel outward through the heliosphere, the shell becomes filled by adiabatic cooling. Generally, one can easily identify PUIs by their VDFs, which are distinctly different from that of SW ions.

In view of the paucity of measurements of keV ions in the outer heliosphere, the transport of PUIs is not yet fully understood. Especially, the production of the suprathermal (high-energy non-Maxwellian) tails on the VDFs is a subject of much debate (Chalov et al. 2003; Fahr & Fichtner 2011). It was suggested that during their propagation through the outer heliosphere, PUIs experience pitch-angle scattering and stochastic acceleration by interactions with different kinds of SW turbulences (Chalov et al. 1997).

PUI kinetic transport theory describes the evolving PUI distribution in phase space. The transport equation, an equation of "motion" for the distribution function, is the basis for almost all work on PUI transport. An early paper by Vasyliunas & Siscoe (1976) calculated an isotropic PUI VDF without energy diffusion but including the effects of adiabatic deceleration. The VDF appears as a thick shell centered around the SW beam. Isenberg (1987) investigated how PUI distribution varies with radial distance and phase space speed. He has shown that the effects of energy diffusion can be important and can lead to substantial particle acceleration. Chalov et al. (1995) studied the influence of different representations for the energy diffusion coefficient, which would result from different radial variations of relative fluctuation amplitudes of MHD turbulence in the outer heliosphere. These authors have calculated the energy spectra of PUIs upstream of the heliospheric termination shock (TS) on the basis of realistic PUI production rates. They showed that second-order Fermi acceleration by means of Alfvenic turbulence produces the suprathermal tail. Fahr & Fichtner (2011) replaced the adiabatic cooling with the "magnetic cooling" process resulting from the conservation of the first and second adiabatic invariants. They have demonstrated that small second-order Fermi acceleration and magnetic cooling generates a PUI VDF with f ∝ v⁻⁵. Intriligator et al. (2012) first compared simulated PUI densities and temperatures with the SWICS measurements on board Ulysses. The model results matched well with the observations at about 5.2 au during both quiet periods and the disturbed period during the Halloween 2003 storm. Assuming that the PUI distribution functions are κ distributions, Fahr et al. (2014) proceeded from the phase space transport equation to a pressure equation for the κ parameter κ = κ(r) as function of heliocentric distance r. They obtained the range of possible radial variations of κ from relatively high values in the inner heliosphere to values between 1.5 and about 2 further out, depending on the diffusion coefficient.

The purpose of the present paper is to place PUI transport in the context of the global picture of the SW-LISM interaction. The treatment is based on the more conventional view, where wave–particle interactions ensure rapid isotropization of the distribution. A parallel can be drawn with the approach of Gamayunov et al. (2012), who used a grid-based model for the isotropic PUI distribution function and for the waves generated by the PUI ring anisotropization process, but neglected PUI momentum diffusion. Their model was applied to the supersonic SW upstream of the TS. Here we investigate the detailed spatial distribution of PUIs in both the supersonic SW and the heliosheath. We take into account the effects of convection with the SW, adiabatic cooling, second-order Fermi process and ionization. Similar to Gamayunov et al. (2012), the model combines the MHD treatment of the background SW and neutral atoms with a kinetic treatment of PUIs in the isotropic approximation.

The rest of this paper proceeds as follows. Section 2 explains how we simulate the SW-LISM interactions. The PUI transport model is introduced in Section 3. Section 4 presents the simulated PUI phase space density and energy spectra. Finally, Section 5 mentions the weaknesses of the presented model and proposes several improvements to the current methods.

The heliosphere is a low-density bubble embedded in the local interstellar cloud. The SW slows down when interacting with the LISM. At the TS, a standing shock wave, the SW speed falls below the effective speed of sound (that includes a contribution from PUIs) and becomes subsonic. The TS causes compression, heating, and a change in magnetic field. The SW continues to slow down as it passes through the heliosheath, a transitional zone bounded by the TS and the heliopause, where the pressures from the LISM and SW are balanced. The LISM pressure causes the heliosphere to develop into a comet-like structure. The nose of the heliopause is defined as the direction opposite to the Sun's motion through the Local Interstellar Cloud.

We use a three-dimensional MHD model with four neutral hydrogen atom populations to obtain the plasma background. The neutral hydrogen atoms are separated into four species according to the region in which they were created. Studies comparing multiple neutral fluid and the kinetic Monte Carlo approaches to model neutrals have been done and the four-fluid model shows excellent agreement for the plasma properties (Pauls et al. 1995; Zank et al. 1996; Heerikhuisen et al. 2006). The four neutral fluids have distinct properties. The Neutral 1 population is made up of the primary interstellar atoms. The Neutral 2 population consists of atoms born in the outer heliosheath. Neutral gas produced in the inner heliosheath constitutes the Neutral 3 population, while Neutral 4 is created in the supersonic SW. The computer model is based on a hexagonal spherical geodesic grid that ensures a uniform partitioning of the surface of a sphere, combined with a concentric nonuniform radial grid (Florinski et al. 2013). We solve a modified set of MHD equations using a finite volume method on this grid. We use 40,962 Voronoi polygons on the sphere and 528 radial shells. Because most PUIs are produced at small heliocentric distances and convected outward, we placed the inner boundary at 1.5 au and the outer boundary at 800 au. The radial cells are smaller near the origin; their width increasing monotonically with radial distance. The z-axis points northward of the solar equator, and the x-axis is in the plane defined by the interstellar helium flow direction (Lallement et al. 2005) and the z-axis. The x, y, z-axes constitute a right-handed orthogonal system.

The model heliosphere corresponds approximately to the solar-minimum conditions. We assumed the slow wind latitudinal extent angle of 36°. At 1 au, in the slow SW, the bulk velocity u = 430 km s⁻¹ and density n = 5.0 cm⁻³; in the fast SW, u = 725 km s⁻¹ and n = 1.4 cm⁻³ (Ebert et al. 2009; Jian et al. 2011). At 1 au, the radial magnetic field component B_r = 24.5 μG. The azimuthal magnetic field component B_ϕ ∝ 1/u, where u is the SW speed. Since the heliospheric current sheet is too thin to affect the dynamics of the plasma on large scales, we do not include it in our simulation and use a unipolar magnetic field model (Czechowski et al. 2010; Borovikov et al. 2011; Izmodenov & Alexashov 2015; Opher et al. 2015).

The interstellar neutral atom density was 0.11 cm⁻³ and its velocity vector was (−25.9, 0, −2.196) km s⁻¹ in our coordinates (McComas et al. 2012). We assumed the interstellar magnetic field direction toward the center of the Interstellar Boundary Explorer (IBEX) ribbon, with Cartesian components (1.834, −1.866, 1.422) μG (Funsten et al. 2013). The interstellar plasma and neutral temperatures were both equal to 6300 K. The plasma interacts with the hydrogen atoms through charge exchange. These atoms are modeled using a gas-dynamic model in a fashion similar to our MHD system except that the magnetic field is set to be zero. We use the charge exchange terms in Pauls et al. (1995). The MHD code is run until a steady state is reached, which provides the plasma-neutral simulation background for the next step.

The PUIs are non-thermal, and their preferred treatment is kinetic. We assume that PUIs are scattered into isotropic distribution by either ambient preexisting or self-excited waves, so the VDF depends only on the absolute value of velocity (in the plasma frame).

Schwadron et al. (1996) obtained the mean field strength B₀ from one-day averages measured by Ulysses, and averaged the square of the variations over 5 days, ${\eta }^{2}=\langle {(B-{B}_{0})}^{2}/{B}_{0}^{2}\rangle$, during the 150 days from 1992 February 19 to July 18. They found that Ulysses crossed flux tubes with different values of η², in the range from a high value of ${\eta }_{H}^{2}\approx 0.12$ (highly turbulent) to a low value of ${\eta }_{L}^{2}\approx 0.005$ (almost laminar). As long as the observed distribution functions are averaged over more than about 10 days, they can be viewed as a superposition of these high and low values of η². Following this idea, we separate the PUI distribution into two components. The first component ("core") is found in laminar flux tubes, while the second component ("tail") resides in turbulent flux tubes, where they experience stochastic acceleration.

The distributions f_core and f_tail satisfy the following transport equations.

$$\begin{eqnarray}&&\displaystyle \frac{\partial {f}_{\mathrm{core}}}{\partial t}+{\boldsymbol{u}}\cdot {\rm{\nabla }}{f}_{\mathrm{core}}-\displaystyle \frac{{\boldsymbol{\nabla }}\cdot u}{3}\displaystyle \frac{\partial {f}_{\mathrm{core}}}{\partial \mathrm{ln}v}={\rm{\Gamma }}\ast S,\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial {f}_{\mathrm{tail}}}{\partial t}+{\boldsymbol{u}}\cdot {\rm{\nabla }}{f}_{\mathrm{tail}}-\displaystyle \frac{{\boldsymbol{\nabla }}\cdot u}{3}\displaystyle \frac{\partial {f}_{\mathrm{tail}}}{\partial \mathrm{ln}v}\\ &&\quad -\displaystyle \frac{1}{{v}^{2}}\displaystyle \frac{\partial }{\partial v}\left({v}^{2}D\displaystyle \frac{\partial {f}_{\mathrm{tail}}}{\partial v}\right)=(1-{\rm{\Gamma }})\ast S\end{eqnarray} \tag{ 2 }$$

with Γ = 0.98, where v is PUI velocity in the non-inertial frame moving with the plasma, ${\boldsymbol{u}}$ is the bulk velocity of the SW plasma, D is the velocity diffusion coefficient, and S is the source term, describing PUIs produced from charge exchange or some other ionization process. The left-hand sides of the transport equations describes the explicit time dependence of the PUI velocity distribution, convection, adiabatic cooling or heating, and velocity diffusion, respectively. PUIs typically have energies between 100 eV and 100 keV. We neglect spatial diffusion and drift motions since the typical spatial diffusive scale of PUI is much smaller than our simulation scale, and the drift speed is small in this low-energy range (Rucinski et al. 1993).

The velocity diffusion coefficient is $D={v}^{2}\delta u/(9{\xi }_{c})$, where δu is the rms fluctuating velocity of the plasma and ξ_c is the correlation length (Chalov et al. 1997). It includes a multitude of effects, such as transit-time damping (Fisk 1976), parallel electric fields in two-dimensional turbulence (le Roux et al. 2002), and large-scale compressive structures (Chalov et al. 2003). The values of δu and ξ_c were (29 km s⁻¹, 0.3 au) in the inner heliosheath and (24 km s⁻¹, 1 au) in the supersonic SW. These parameters were chosen to produce suprathermal tails in a way that would be consistent with the Voyager observations (see below).

In the source term, charge exchange and photoionization are taken into account, but electron impact ionization is neglected. The distribution of PUIs (variables bearing a subscript "i") created by the charge exchange process holds imprints from the parent atom distribution (index "n") and the background plasma (index "p"). We assume that both are Maxwellian with densities n, mean velocities ${\boldsymbol{u}}$, and thermal speeds ${v}_{T}={(2{kT}/m)}^{1/2}$,

$$\begin{eqnarray}&&{f}_{p}({{\boldsymbol{v}}}_{p})=\displaystyle \frac{{n}_{p}}{{\pi }^{3/2}{v}_{{T}_{p}}^{3}}{e}^{-{({{\boldsymbol{v}}}_{p}-{{\boldsymbol{u}}}_{p})}^{2}/{v}_{{T}_{p}}^{2}},\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&{f}_{n}({{\boldsymbol{v}}}_{n})=\displaystyle \frac{{n}_{n}}{{\pi }^{3/2}{v}_{{T}_{n}}^{3}}{e}^{-{({{\boldsymbol{v}}}_{n}-{{\boldsymbol{u}}}_{n})}^{2}/{v}_{{T}_{n}}^{2}}.\end{eqnarray} \tag{ 4 }$$

The PUI production term from charge exchange is

$$\begin{eqnarray}{\left(\displaystyle \frac{\delta {f}_{i}({v}_{i})}{\delta t}\right)}_{{\rm{ex}}} & = & \displaystyle \int \sigma (| {{\boldsymbol{v}}}_{p}-{{\boldsymbol{v}}}_{n}| ){f}_{p}({{\boldsymbol{v}}}_{p}){f}_{n}({{\boldsymbol{v}}}_{n})| {{\boldsymbol{v}}}_{p}-{{\boldsymbol{v}}}_{n}| \displaystyle \frac{{d}^{3}{v}_{p}{d}^{3}{v}_{n}}{{d}^{3}{v}_{i}}\\ & = & \sigma ({\rm{\Delta }}{v}_{p}({{\boldsymbol{v}}}_{i})){f}_{n}({{\boldsymbol{v}}}_{i})\displaystyle \int {f}_{p}({{\boldsymbol{v}}}_{p})| {{\boldsymbol{v}}}_{p}-{{\boldsymbol{v}}}_{i}| {d}^{3}{v}_{p}\\ & = & {n}_{p}\sigma ({\rm{\Delta }}{v}_{p}({{\boldsymbol{v}}}_{i})){\rm{\Delta }}{v}_{p}({{\boldsymbol{v}}}_{i}){f}_{n}({{\boldsymbol{v}}}_{i}),\end{eqnarray} \tag{ 5 }$$

where ${\rm{\Delta }}{v}_{p}({{\boldsymbol{v}}}_{i})$ is the average relative speed between a given PUI and the proton population. To obtain the above equation, the condition for charge exchange, ${{\boldsymbol{v}}}_{n}={{\boldsymbol{v}}}_{i}$, and the experimental result that the cross section σ(v) is only weakly dependent on the relative velocity between the particles, were used. Numerical values of σ are available from Lindsay & Stebbings (2005). The average relative speed is computed as (Fahr & Müller 1967)

$$\begin{eqnarray}{\rm{\Delta }}{v}_{p}({{\boldsymbol{v}}}_{i}) & = & \displaystyle \frac{1}{{\pi }^{3/2}{v}_{{T}_{p}}^{3}}\displaystyle \int | {{\boldsymbol{v}}}_{p}-{{\boldsymbol{v}}}_{i}| {e}^{-{({{\boldsymbol{v}}}_{p}-{{\boldsymbol{u}}}_{p})}^{2}/{v}_{{T}_{p}}^{2}}{d}^{3}{v}_{p}\\ & = & {v}_{{T}_{p}}\left[\displaystyle \frac{1}{\sqrt{\pi }}{e}^{-{x}^{2}}+\left(x+\displaystyle \frac{1}{2x}\right){\rm{Erf}}(x)\right],\end{eqnarray} \tag{ 6 }$$

where $x=| {{\boldsymbol{v}}}_{i}-{{\boldsymbol{u}}}_{p}| /{v}_{{T}_{p}}$.

The PUI transport equation is solved in the plasma frame, so we transform (5) into that frame and average over a sphere. We introduce the PUI velocity in the non-inertial frame moving with the plasma, ${\boldsymbol{v}}{^{\prime} }_{i}={{\boldsymbol{v}}}_{i}-{{\boldsymbol{u}}}_{p}={\boldsymbol{x}}{v}_{{T}_{p}}$. Notice that in this frame (6) is isotropic and requires no averaging. Then

$$\begin{eqnarray}{\left(\displaystyle \frac{\delta {f}_{i}({v}_{i}^{\prime })}{\delta t}\right)}_{{\rm{ex}}} & = & \displaystyle \frac{{n}_{p}{n}_{n}\sigma ({\rm{\Delta }}{v}_{p}({v}_{i}^{\prime })){\rm{\Delta }}{v}_{p}({v}_{i}^{\prime })}{4{\pi }^{5/2}{v}_{{T}_{n}}^{3}}\\ & & \times \,\displaystyle \int {e}^{-{({{\boldsymbol{v}}}_{i}^{\prime }-{{\boldsymbol{u}}}_{n}+{{\boldsymbol{u}}}_{p})}^{2}/{v}_{{T}_{n}}^{2}}d{{\rm{\Omega }}}_{i}^{\prime }\\ & = & \displaystyle \frac{{n}_{p}{n}_{n}\sigma ({\rm{\Delta }}{v}_{p}(x)){\rm{\Delta }}{v}_{p}(x)}{2{\pi }^{3/2}{v}_{{T}_{n}}^{3}}\\ & & \times \,{\displaystyle \int }_{-1}^{1}{e}^{-({\alpha }^{2}{x}^{2}-2\mu ^{\prime} \alpha {xh}+{h}^{2})}d\mu ^{\prime} \\ & = & \displaystyle \frac{{n}_{p}{n}_{n}\sigma ({\rm{\Delta }}{v}_{p}(x)){\rm{\Delta }}{v}_{p}(x)}{4{\pi }^{3/2}{v}_{{T}_{p}}{v}_{{T}_{n}}^{2}}\\ & & \times \,\displaystyle \frac{{e}^{-{(\alpha x-h)}^{2}}-{e}^{-{(\alpha x+h)}^{2}}}{{xh}},\end{eqnarray} \tag{ 7 }$$

where $x=v^{\prime} /{v}_{{T}_{p}}$, $h=| {{\boldsymbol{u}}}_{n}-{{\boldsymbol{u}}}_{p}| /{v}_{{T}_{n}}$, and $\alpha ={v}_{{T}_{p}}/{v}_{{T}_{n}}$, are the charge exchange source terms for an isotropic distribution of PUIs. Similarly, the photoionization source term is

$$\begin{eqnarray}&&{\left(\displaystyle \frac{\delta {f}_{i}({v}_{i}^{\prime })}{\delta t}\right)}_{{\rm{ph}}}=\displaystyle \frac{{n}_{n}\nu }{4{\pi }^{3/2}{v}_{{T}_{p}}{v}_{{T}_{n}}^{2}}\displaystyle \frac{{e}^{-{(\alpha x-h)}^{2}}-{e}^{-{(\alpha x+h)}^{2}}}{{xh}},\end{eqnarray} \tag{ 8 }$$

where $\nu ={\nu }_{\mathrm{ph}}{(1\mathrm{au}/r)}^{2}$ with ν_ph being the rate of photoionization at 1 au.

We solve the PUI transport Equations (1) and (2) by integrating them on a multi-CPU cluster simultaneously with the MHD system. A conservative form of these equations is used, namely,

$$\begin{eqnarray}&&\displaystyle \frac{\partial f}{\partial t}+{\boldsymbol{\nabla }}\cdot ({\boldsymbol{u}}f)-\displaystyle \frac{1}{{v}^{2}}\displaystyle \frac{\partial }{\partial v}\left[\displaystyle \frac{(v{\boldsymbol{\nabla }}\cdot u)}{3}{v}^{2}f+{{Dv}}^{2}\displaystyle \frac{\partial f}{\partial v}\right]=S,\end{eqnarray} \tag{ 9 }$$

and similar for Equation (2). The transport module is implemented on the same grid as the MHD simulation, using a finite volume method. Right and left interface values are computed using piecewise linear reconstruction with a WENO limiter (Jiang & Shu 1996; Friedrich 1998). The ionic components all have the same bulk speed available from the MHD solution. The velocity grid extends from 10 km s⁻¹ to 6000 km s⁻¹.

The TS and the heliopause are at 90 au and 144 au, respectively, along the Voyager 1 spacecraft direction in our simulations. Voyager 1 actually crossed the TS and the heliopause at 94 au and 122 au (Stone et al. 2005, 2013). These distances are appropriate for a model that is time independent and is based on solar-minimum conditions. In Figure 1, the left panel shows the simulated PUI VDF in the Voyager 1 direction (θ = 55°, ϕ = 0°, where θ denotes the co-latitude and ϕ denotes the longitude), in the plasma frame. At r = 5 au, the PUI distribution shows a rapid drop at around 430 km s⁻¹, the bulk velocity of SW inside the TS, indicating that most particles are injected with the speed of the SW V_SW in the plasma frame. Some particles have filled in the shell at low energies due to adiabatic cooling. Since particles with speeds of v > V_SW result from local acceleration in turbulent flux tubes, the mixture of core PUIs and tail PUIs yields this step-like feature (Schwadron et al. 1996). These distributions are to be interpreted in the time-averaged sense, over many flux tubes passing an observer.

Zoom In Zoom Out Reset image size

A power-law suprathermal tail develops at v > V_SW. The tail extends to higher energies with increasing radial distance. In view of the flux tube picture, an observer would alternately see distributions with more or less developed tails; Figure 1 is their average over several days or even months (the model is not time dependent; therefore, the averaging interval can be arbitrarily long). A hump at around 100 km s⁻¹ at 100 au is clearly seen, which consists of the low-energy PUIs created in the inner heliosheath. The hump increases in height as the flow slows down toward the heliopause due to the accumulation of low-energy PUIs produced in the heliosheath.

The right panel of Figure 1 illustrates the variation of the PUI VDF along the polar direction (θ = 0°, ϕ = 0°). The principle features are similar to the equatorial case. At r = 5 au, the PUI distribution shows a sharp change in slope around 725 km s⁻¹, the bulk velocity of the fast SW inside the TS. A pronounced hump develops only beyond the TS (at about 90 au). In both panels, one can see that the PUI gas gets compressed by the shock and an accelerated power-law tail develops out of the PUI core. There is little change in the power-law tail in the heliosheath, which means stochastic acceleration cannot produce a spectrum that is any harder. The downstream value of the power-law index is about −5.2 in both directions.

The distribution of background plasma density is shown in Figure 2, while Figure 3 presents the core and tail PUI spatial distributions. They show that from the inner boundary to the TS, the PUI density falls off slower than that of the SW core since PUIs are produced over the entire SW. One can clearly see that there is a great difference in PUI densities in the slow SW and in the fast SW. Figures 3(a)–(c) show that at low energies the PUI density is higher in the slow SW than in the fast SW. The majority of the maps show an increase in PUI densities across the TS due to compression. Naturally, PUIs created in a slow wind have energy lower than those produced in a fast wind. One can see that most low-energy PUIs are in the heliotail from Figure 3(a). This is because the PUIs produced in the direction of the heliotail have a lower energy on average because of greater SW slowdown by mass loading before the TS.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The model is tuned to fit the observations. For example, Gloeckler & Geiss (1998) reported an averaged phase space density measured by SWICS in a 100-day interval in 1994 when Ulysses was at about 3 au from the southern pole. To compare the observed velocity distribution with model predictions, we transformed our model distribution functions into the frame of SWICS, and then integrated over the SWICS field of view. The results are presented in Figure 4. The red line is the simulated PUI phase space density after transformation. One can see that our model result matches the observed phase space density well.

Zoom In Zoom Out Reset image size

Recent measurements by Solar Wind Around Pluto (SWAP) on board New Horizons were reported in Randol et al. (2013). SWAP measured SW and PUI spectra from 11 au to 22 au. We transformed our model distribution function to the spacecraft frame and converted it into count rate. The comparison is given in Figure 5(a). The SW H⁺ and He²⁺ peaks can be seen at energy per charge (E/q) ≈ 600 eV/e and 1100 eV/e, respectively. As can be seen, the slope of the spectrum and the PUI-step-like feature are well reproduced. At all other energies, the PUIs are hidden beneath the SW background. The spectrum is below the data at v > V_SW. This is partially due to contamination from SW, but also due to local acceleration of the PUIs by variations in magnetic field or the SW background.

Zoom In Zoom Out Reset image size

Figure 5(b) compares the simulated PUI intensities j(E) with Voyager 1 Low Energy Charged Particle (LECP) proton intensities averaged over one selected 78-day interval in the termination foreshock (2004/167-2004/245) and one 78-day period immediately behind the TS (2004/349-2005/061) (Decker et al. 2008). We also compare simulated PUI intensity with Voyager 2 low-energy proton intensity averaged over a selected 78-day interval in the inner heliosheath (2007/241-2007/319) in Figure 5(c). Note that the spectral shapes in Figures 5(b) and (c) are essentially the same, well represented by j ∝ E^−1.5 corresponding to f ∝ v⁻⁵, with the Voyager 2 ion spectrum being slightly harder than that at Voyager 1 behind the TS crossing. Each of the three spectra have the step feature. The spectra shift upward once the TS is crossed.

A partial overlap exists between the Voyager 1, 2 LECP ion instrument (28 < E < 4000 keV) and the Cassini ENA Ion and Neutral Camera (INCA) sensor (5.2 < E < 55 keV) (Krimigis et al. 2010). Krimigis et al. (2010) reported the INCA-inferred ion spectrum in the heliosheath and matched it to the in situ measured Voyager 2 spectrum. Figure 5(d) shows our simulated PUI spectrum at ∼90 au in Voyager 2 direction compared with the Cassini ENA INCA-inferred ion spectrum.

Our model introduced several improvements over the existing models. Chalov et al. (2004) have investigated the spatial variation of PUI spectra, but only the upwind part of the heliosheath was considered. Malama et al. (2006) have performed a multi-species simulation, but the magnetic field was ignored and their model was two-dimensional. The model presented here computes PUI distributions on a three-dimensional grid. Usmanov & Goldstein (2006) considered the SW outside 1 au as a combination of three co-moving species, SW protons, electrons, and PUIs, but they only computed the global structure of the SW from the coronal base to 100 au without the TS. Our model's external boundary is well inside the LISM covering supersonic SW region, the inner heliosheath, and the outer heliosheath.

Several weaknesses of the present model are now pointed out. First, the two-population assumption is admittedly questionable. The Schwadron et al. (1996) paper reported a range of field strength variations, rather than a bimodal distribution. The observed distribution functions are averaged over longer times than the averaging interval in Schwadron et al. (1996), and therefore are a product of a superposition of high and low values of η². The simulated spectra shown in Figures 1 and 5 should be compared with monthly or even yearly averages of the data. Second, our model assumes the PUI VDF is isotropic, which is not always the case. For example, ion angular data from Voyager 1 observations during 2002.58–2003.10, 85.3–87.3 au showed large beam-like anisotropies (Decker et al. 2005). Third, our model is time independent. Incorporating time-dependent SW boundary conditions may improve the results and produce better agreement with the observations. Finally, charge exchange of the PUIs on interstellar H atoms was ignored. This process replaces one PUI with another, drawn from a different velocity distribution. While, in principle, including the loss of PUIs would be straightforward, the production term requires numerical integration over the VDF of PUIs in each computational cell, which is a costly procedure. PUI charge exchange may be important in the inner heliosheath, causing an energy redistribution in their VDF.

In spite of these limitations, our model does provide insights into the interpretation of the PUI data and may be used to predict PUI distribution at all locations inside the heliosphere. These distributions show the details that are directly comparable with those seen in spacecraft data. We have obtained the rapid drops in the spectra that appear to be required to match the observations. The model also features power-law tails in the energy, which are commonly observed in space. A velocity diffusion origin of these tails appears to be a valid interpretation.

The compressed SW and PUIs behind the TS create energetic neutral atoms (ENAs) via charge exchange. ENAs with energies high enough to overcome the outward flow speed can be directed back at Earth. Future work will use these PUI results to calculate ENA fluxes at 1 au. We plan to compare the simulated ENA fluxes with the IBEX distributed ENA sky maps. This will bring us closer to explaining why the distributed ENA flux spectrum does not show a knee and why is it close to a power law (Schwadron et al. 2011).

This work was supported by NASA grant NNX12AH44G.