Constraining the Evolution of the Proton Distribution Function in the Heliotail

Following the work of Zirnstein et al. (2018b), we solve the stationary Parker transport equation with charge-exchange source terms, including (1) the advection of PUIs with the bulk plasma, (2) the adiabatic heating of the plasma derived from our 3D MHD/kinetic model of the heliosphere, and (3) a velocity diffusion term that is left as a free parameter. The transport equation is given as

$$\begin{eqnarray}\begin{array}{rcl}{u}_{p}(s)\displaystyle \frac{\partial {f}_{p}(s,v)}{\partial s} & = & \displaystyle \frac{1}{{v}^{2}}\displaystyle \frac{\partial }{\partial v}\left[D(v){v}^{2}\displaystyle \frac{\partial {f}_{p}(s,v)}{\partial v}\right]\\ & & +\,\displaystyle \frac{v}{3}({\rm{\nabla }}\cdot {{\boldsymbol{u}}}_{p})\displaystyle \frac{\partial {f}_{p}(s,v)}{\partial v}+P(s,v)-L(s,v),\end{array}\end{eqnarray} \tag{ 1 }$$

where v is the proton speed in the plasma frame, u_p is the bulk flow speed along the streamline distance coordinate s, f_p is the isotropic proton distribution function, D(v) is the velocity diffusion coefficient, which we assume is a free parameter, ${\rm{\nabla }}\cdot {{\boldsymbol{u}}}_{p}$ is the bulk plasma flow divergence, P is the production (source) term due to charge exchange (photoionization is negligible, and electron-impact ionization is ignored), and L is the loss (sink) term due to charge exchange with neutral hydrogen atoms. As in our previous work, the plasma flow speed, flow divergence, and neutral hydrogen distribution used to calculate the charge-exchange source terms are taken from our MHD-plasma/kinetic-neutral simulation of the heliosphere (Zirnstein et al. 2016). The velocity diffusion coefficient is written in the form D(v) = D₀ (v/v₀)^γ, where D₀ is the diffusion rate amplitude with units km² s⁻³, γ is the velocity spectral index, v is the particle speed, and v₀ = 1 km s⁻¹ is a normalization speed. Following prior works (e.g., Chalov et al. 2003; Fahr & Scherer 2004; Chalov et al. 2016), we assume that the spatial diffusion time of ∼keV PUIs is much larger than their advection time with the bulk plasma flow in the IHS.

The charge-exchange source terms are given by (see Zirnstein et al. 2018b for details)

$$\begin{eqnarray}&&P(s,v)=\displaystyle \frac{1}{4\pi }\iint \eta (s,{\boldsymbol{v}}){f}_{H}(s,{\boldsymbol{v}}+{{\boldsymbol{u}}}_{p})\sin \theta d\theta d\phi ,\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}\begin{array}{rcl}L(s,v) & \cong & \displaystyle \frac{{f}_{p}(s,v)}{4\pi }\displaystyle \sum _{i}{n}_{H,i}(s){\sigma }_{{ex}}({v}_{{rel},i})\\ & & \times \,\iint \sqrt{\displaystyle \frac{4}{\pi }{v}_{\mathrm{th},H,i}^{2}+{| {\boldsymbol{v}}+{{\boldsymbol{u}}}_{p}-{{\boldsymbol{u}}}_{H,i}| }^{2}}\sin \theta d\theta d\phi .\end{array}\end{eqnarray} \tag{ 3 }$$

Equations (2) and (3) are derived by assuming that the proton distribution function is isotropic in the plasma frame; therefore, we integrate over solid angle in velocity space (see also Malama et al. 2006). Equation (3) is further simplified by the fact that the majority of neutral hydrogen atoms in the heliotail are (1) interstellar neutral atoms from the VLISM and (2) neutral atoms that are produced by charge exchange in the slowed and compressed outer heliosheath (OHS), where both populations are approximately Maxwell–Boltzmann (e.g., Izmodenov 2001; Heerikhuisen et al. 2006). It is important to include the second population since these neutral atoms continue to enter the heliotail from outside the heliopause for long distances from the Sun. Thus, by assuming these two neutral populations are Maxwell–Boltzmann, the charge-exchange loss term can be written in the form of Equation (3), where the double integral has an analytic solution (see Equation (7) of Zirnstein et al. 2018b). The method for solving Equation (1) is briefly summarized in the Appendix, and more details can be found in Zirnstein et al. (2018b).

The purpose of this study is to test different initial proton distribution functions downstream of the TS and analyze their effects on the evolution of the distribution in the heliotail and on ENA fluxes at 1 au. We test three analytic forms for the initial proton distribution function (f_p,0) and one distribution extracted from a particle-in-cell (PIC) simulation of the TS. The three analytic functions are (1) a single kappa function representing the entire proton distribution, while varying the kappa index, (2) a superposition of a Maxwell–Boltzmann distribution representing the core SW ions and a filled-shell distribution (Vasyliunas & Siscoe 1976) for PUIs, and (3) a superposition of three Maxwell–Boltzmann distributions representing core SW ions, PUIs that are transmitted across the TS, and PUIs that are singularly reflected at the TS, then transmitted (Zank et al. 2010). The analytic distributions and their properties are described in Table 1. In Section 4.3, we describe the PIC simulation results.

First, we use a kappa distribution of the form (e.g., Livadiotis & McComas 2013)

$$\begin{eqnarray}\begin{array}{rcl}{f}_{\kappa }(v) & = & \displaystyle \frac{{n}_{p,0}}{{\pi }^{\tfrac{3}{2}}{{\rm{\Theta }}}_{p,0}^{3}}\displaystyle \frac{1}{{\left[{\kappa }_{p,0}-\tfrac{3}{2}\right]}^{\tfrac{3}{2}}}\displaystyle \frac{{\rm{\Gamma }}({\kappa }_{p,0}+1)}{{\rm{\Gamma }}\left({\kappa }_{p,0}-\tfrac{1}{2}\right)}\\ & & \times {\left[1+\displaystyle \frac{1}{\left[{\kappa }_{p,0}-\tfrac{3}{2}\right]}\displaystyle \frac{{v}^{2}}{{{\rm{\Theta }}}_{p,0}^{2}}\right]}^{-{\kappa }_{p,0}-1},\end{array}\end{eqnarray} \tag{ 4 }$$

with initial density n_p,0, temperature T_p,0, kappa index κ_p,0, and characteristic speed ${{\rm{\Theta }}}_{p,0}=\sqrt{2{k}_{B}{T}_{p,0}/{m}_{p}}$. The initial density and temperature are assumed to be 0.0017 cm⁻³ and 1.9 × 10⁶ K, which are taken from our MHD simulation (i.e., downstream of our simulated TS in the VLISM downwind direction; Zirnstein et al. 2016). We vary the kappa index in the range 1.51 ≤ κ_p,0 ≤ 10. We note that the variable κ_p,0 here is equivalent to the 3D kappa index from Livadiotis & McComas (2013), not their invariant kappa index.

The second distribution that we test is a combination of a Maxwell–Boltzmann distribution for the core SW ions and a filled-shell distribution for PUIs (i.e., Vasyliunas & Siscoe 1976) that is compressed at the TS. Upstream of the TS (subscript "u"), the proton distribution is

$$\begin{eqnarray}\begin{array}{rcl}{f}_{p,u}(v) & = & \displaystyle \frac{{n}_{p,u}(1-\alpha )}{{(\pi {v}_{\mathrm{th},p,u}^{2})}^{3/2}}\exp \left[-{\left(\displaystyle \frac{v}{{v}_{\mathrm{th},p,u}}\right)}^{2}\right]\\ & & +\,\displaystyle \frac{3\alpha {n}_{p,u}{\rm{\Theta }}(v-{u}_{p,u})}{8\pi {u}_{p,u}^{3/2}{v}^{3/2}},\end{array}\end{eqnarray} \tag{ 5 }$$

where n_p,1 is the total upstream proton density, α is the PUI density fraction (assumed to be 25%; e.g., see Figure 2 in Zirnstein et al. 2017), u_p,u = 350 km s⁻¹ is the upstream SW speed, and v_th,p,1 is the upstream thermal velocity of the core SW (assumed to have a temperature of 10⁴ K, similar to that observed by Voyager 2 upstream of the TS; Richardson et al. 2008). Following Zank et al. (2010), we derive the downstream proton distribution by transporting the upstream distribution across the TS under the assumption that the particles and bulk flow are decelerated by the cross-shock potential, yielding the downstream particle speed (in the flow frame) to be v_d = rv_u (see Zank et al. 2010 for more details), where r = 3 is the shock compression ratio. Using Liouville's theorem, we derive the downstream proton distribution f_p,d in terms of the upstream distribution f_p,u as

$$\begin{eqnarray}\begin{array}{rcl}4\pi {u}_{p,d}{f}_{p,d}({v}_{d}){v}_{d}^{2}{{dv}}_{d} & = & 4\pi {u}_{p,u}{f}_{p,u}({v}_{u}){v}_{u}^{2}{{dv}}_{u}\\ {f}_{p,d}({v}_{d}) & = & \displaystyle \frac{{u}_{p,u}}{{u}_{p,d}}{f}_{p,u}({v}_{u})\displaystyle \frac{{v}_{u}^{2}}{{v}_{d}^{2}}\displaystyle \frac{{{dv}}_{u}}{{{dv}}_{d}}=\displaystyle \frac{1}{{r}^{2}}{f}_{p,u}\left(\displaystyle \frac{{v}_{d}}{r}\right).\end{array}\end{eqnarray} \tag{ 6 }$$

This yields the downstream proton distribution:

$$\begin{eqnarray}\begin{array}{rcl}{f}_{\mathrm{MS}}(v) & = & {f}_{p,d}(v)=\displaystyle \frac{1}{{r}^{2}}\displaystyle \frac{{n}_{p,u}(1-\alpha )}{{(\pi {v}_{\mathrm{th},p,u}^{2})}^{3/2}}\exp \left[-{\left(\displaystyle \frac{v/r}{{v}_{\mathrm{th},p,u}}\right)}^{2}\right]\\ & & +\,\displaystyle \frac{1}{\sqrt{r}}\displaystyle \frac{3\alpha {n}_{p,u}{\rm{\Theta }}(v/r-{u}_{p,u})}{8\pi {u}_{p,u}^{3/2}{v}^{3/2}}.\end{array}\end{eqnarray} \tag{ 7 }$$

Since n_p,0 is the total proton density downstream of the TS, then r × n_p,u = n_p,d = n_p,0. Thus, Equation (7) can be rewritten as

$$\begin{eqnarray}\begin{array}{rcl}{f}_{\mathrm{MS}}(v) & = & \displaystyle \frac{1}{{r}^{3}}\displaystyle \frac{{n}_{p,0}(1-\alpha )}{{(\pi {v}_{\mathrm{th},p,u}^{2})}^{3/2}}\exp \left[-{\left(\displaystyle \frac{v/r}{{v}_{\mathrm{th},p,u}}\right)}^{2}\right]\\ & & +\,\displaystyle \frac{1}{{r}^{3/2}}\displaystyle \frac{3\alpha {n}_{p,0}{\rm{\Theta }}(v/r-{u}_{p,u})}{8\pi {u}_{p,u}^{3/2}{v}^{3/2}}.\end{array}\end{eqnarray} \tag{ 8 }$$

The initial proton distribution parameters for Equation (8) are shown in Table 1, which are taken from our 3D MHD simulation of the heliosphere (Zirnstein et al. 2016).

The third distribution that we test is a superposition of three Maxwell–Boltzmann functions, representing (1) the core SW ions, (2) PUIs transmitted across the TS, and (3) PUIs reflected and then transmitted across the TS. This composite distribution is commonly used to model multicomponent ion species in the heliosheath (e.g., Zank et al. 2010; Zirnstein et al. 2014, 2017; Swaczyna & Bzowski 2017a; Swaczyna et al. 2017b). Unlike the case where we process the filled-shell distribution across the TS, here we simply assume that the downstream distribution is represented by three Maxwell–Boltzmann distributions with total density n_p,0 and effective temperature T_p,0. Following Zank et al. (2010), we directly write the downstream distribution f_3M as

$$\begin{eqnarray}\begin{array}{rcl}{f}_{3M}(v) & = & \displaystyle \frac{{n}_{p,0}(1-\alpha )}{{(\pi {v}_{\mathrm{th},p,\mathrm{SW}}^{2})}^{\tfrac{3}{2}}}\exp \left[-{\left(\displaystyle \frac{v}{{v}_{\mathrm{th},p,\mathrm{SW}}}\right)}^{2}\right]\\ & & +\,\displaystyle \frac{{n}_{p,0}\alpha (1-\beta )}{{(\pi {v}_{\mathrm{th},p,\mathrm{tr}}^{2})}^{\tfrac{3}{2}}}\exp \left[-{\left(\displaystyle \frac{v}{{v}_{\mathrm{th},p,\mathrm{tr}}}\right)}^{2}\right]\\ & & +\,\displaystyle \frac{{n}_{p,0}\alpha \beta }{{(\pi {v}_{\mathrm{th},p,\mathrm{ref}}^{2})}^{3/2}}\exp \left[-{\left(\displaystyle \frac{v}{{v}_{\mathrm{th},p,\mathrm{ref}}}\right)}^{2}\right],\end{array}\end{eqnarray} \tag{ 9 }$$

where α is the PUI density fraction (25%), β is the fraction of PUIs that are reflected at the TS (8%; Zank et al. 2010; Chalov & Fahr 1996), and v_th,p,SW, v_th,p,tr, and v_th,p,ref are the initial thermal speeds of core SW ions, transmitted PUIs, and reflected PUIs, respectively (${v}_{\mathrm{th},p,i}=\sqrt{2{k}_{B}{T}_{p,i}/{m}_{p}}$). We choose the thermal speeds (i.e., temperatures) for each component such that their total thermal pressure satisfies the thermal pressure given by our MHD simulation downstream of the TS, n_p,0k_BT_p,0,

$$\begin{eqnarray}&&{T}_{p,0}=\left(\displaystyle \frac{{n}_{p,\mathrm{SW}}}{{n}_{p,0}}{{\rm{\Gamma }}}_{p,\mathrm{SW}}+\displaystyle \frac{{n}_{p,\mathrm{tr}}}{{n}_{p,0}}{{\rm{\Gamma }}}_{p,\mathrm{tr}}+\displaystyle \frac{{n}_{p,\mathrm{ref}}}{{n}_{p,0}}{{\rm{\Gamma }}}_{p,\mathrm{ref}}\right){T}_{p,0},\end{eqnarray} \tag{ 10 }$$

where n_p,0 = 0.0017 cm⁻³, T_p,0 = 1.9 × 10⁶ K, and the temperature fractions Γ_p,SW = T_p,SW/T_p,0, Γ_p,tr = T_p,tr/T_p,0, and Γ_p,ref = T_p,ref/T_p,0. Knowing their thermal pressure fractions (Table 1; see also Zank et al. 2010; Zirnstein et al. 2017), the initial temperature of each component is T_p,SW = 1.0 × 10⁵ K, T_p,tr = 4.1 × 10⁶ K, and T_p,ref = 4.4 × 10⁷ K.

Examples of each distribution (Equations (4), (8), and (9)) are shown in Figure 2 for the same total density and temperature (for the kappa distribution, we show several cases for κ_p,0). Using these distributions as initial boundary conditions, we solve the Parker transport equation (Equation (1)) for a range of velocity diffusion coefficients. In the following sections, we compute ENA fluxes at 1 au from the model proton distributions and provide a quantitative comparison with IBEX data to determine their physical likelihood.

Zoom In Zoom Out Reset image size

In this section, we analyze how the diffusion coefficient in the heliotail and model ENA spectrum at 1 au depend on the initial proton distribution assuming it is a kappa function (the first of three cases, see Equation (4)). Figure 4 shows model ENA fluxes assuming that f_p,0 is a kappa distribution for a selection of kappa index values in the range 1.51 ≤ κ_p,0 ≤ 10. Figure 4(b) shows the case where the velocity diffusion coefficient is close to the best fit derived by Zirnstein et al. (2018b) for κ_p,0 = 1.63. For small κ_p,0 (1.51) the resulting ENA spectrum is lower than those with larger proton kappa indices. While one might expect to see more high energy ENAs from a proton distribution with a small kappa index, protons in a κ_p,0 ∼ 1.5 distribution are at much higher energies than what IBEX can observe. A distribution with 5 < κ_p,0 < 10 actually produces the most ENAs observable by IBEX in the ∼0.1–5 keV energy range.

Zoom In Zoom Out Reset image size

Depending on the diffusion rate, ENA spectra resulting from proton distributions with different κ_p,0 can significantly deviate from each other at energies >0.5 keV, suggesting that different diffusion rates are required to best fit IBEX data. Figure 4(a) shows that if the diffusion rate amplitude D₀ was slightly lower (∼5 × 10⁻⁹ km² s⁻³), then it would appear that a larger kappa index would be needed to better match IBEX data; however, in this case none of the spectra can reproduce the entire IBEX spectrum. On the other hand, a larger diffusion coefficient amplitude (D₀ ∼ 2 × 10⁻⁸ km² s⁻³; Figure 4(c)) would produce similar ENA spectra regardless of κ_p,0, all of which overestimate IBEX measurements. The relationship between κ_p,0 and D₀ is illustrated in more detail in the left panel of Figure 5, which shows the reduced chi-square values between the model and IBEX data as a function of κ_p,0 and D₀, setting γ = 1.25. Figure 5 shows a broad region of reduced chi-square less than or equal to 1 when κ_p,0 ≲ 7 and 8 × 10⁻⁹ km² s⁻³ ≲ D₀ ≲ 1.3 × 10⁻⁸ km² s⁻³. As κ_p,0 approaches 1.5, a significantly larger diffusion rate is required to diffuse particles from low energies to the IBEX energy range. However, as κ_p,0 increases, the best diffusion rate reaches an asymptotic value of ∼9 × 10⁻⁹ km² s⁻³.

Zoom In Zoom Out Reset image size

The right panel of Figure 5 shows three different ENA spectra when κ_p,0 = 1.51, 2, and 10, with different D₀ that give the best fit to IBEX data in each case. While the reduced chi-square value is slightly better for κ_p,0 ≲ 7, one can see that even the case where κ_p,0 = 10 matches IBEX data reasonably well. Our results suggest that there is not sufficient resolution in IBEX data to determine the best κ_p,0; however, this may be possible with higher resolution measurements from the future Interstellar Mapping and Acceleration Probe (IMAP) mission (https://imap.princeton.edu). Nonetheless, we are able to determine that the diffusion coefficient is limited to the range ∼8 × 10⁻⁹ km² s⁻³ ≲ D₀ ≲ 1.3 × 10⁻⁸ km² s⁻³.

So far, we have analyzed our model results under the assumption that the velocity diffusion rate in the heliotail is approximately proportional to v^1.3. While this may be true if κ_p,0 ∼ 1.63, this may not be true for smaller or larger κ_p,0. Figure 6 shows results from comparing model ENA fluxes with IBEX data where the diffusion coefficient D(v) scales between v^2/3 and v², which can arise from isotropic Kolmogorov (e.g., Isenberg 1987) and anisotropic Goldreich–Sridhar (e.g., Goldreich & Sridhar 1995; Eichler 2014) turbulence, respectively. We vary the diffusion amplitude to demonstrate how the model compares to IBEX data under different diffusion rates. As one can see in the left panels of Figure 6, the minimum reduced chi-square slightly depends on κ_p,0. For small κ_p,0 near ∼1.5, it appears that D(v) ∝ v^1.3 yields a slightly better fit to the data. This is clearer in the right panels of Figure 6 where we show ENA spectra for the best-fit cases when the diffusion spectral index γ = 0.67, 1.25, and 2. While the case for D(v) ∝ v^1.3 appears to be slightly better at lower energies, the model results for D(v) ∝ v^2/3 and v² are approximately within the 1σ measurement uncertainties. Similar results occur for larger κ_p,0 (κ_p,0 = 2, 10). Thus, while we cannot constrain the best κ_p,0 with significant certainty, we find that a diffusion coefficient that scales between ∼v^2/3 and ∼v² is most consistent with IBEX data for any κ_p,0. In order to balance the diffusion rate, the coefficient amplitude D₀ must be small (∼2–4 × 10⁻¹⁰ km² s⁻³) if D(v) ∝ v² and large (∼1.5–2 × 10⁻⁷ km² s⁻³) if D(v) ∝ v^2/3.

Zoom In Zoom Out Reset image size

In this section, we compare model ENA fluxes for different proton distribution functions: f_p,0 is represented by a kappa distribution (we assume κ_p,0 = 1.63 hereafter), a filled-shell distribution for PUIs, and a superposition of three Maxwell–Boltzmann distributions (see Table 1). We also vary the diffusion amplitude D₀ around ∼10⁻⁸ km² s⁻³ to illustrate how slight changes in D₀ may affect the results. The results are shown in Figure 7. Figure 7(a) shows that while each model produces similar results for energies <1 keV, at higher energies they deviate. The model that assumes PUIs are described by a filled-shell distribution downstream of the TS (f_MS) that produces too many ENAs at energies above ∼2 keV. This is also clear in the results shown by Zank et al. (2010) using a simplified transport model. In fact, this is independent of the diffusion coefficient amplitude (compare Figures 7(a)–(c)). If there is no velocity diffusion, this model would still produce too many ENAs above ∼2 keV and too few ENAs at lower energies. Thus, we can already see that this type of interaction at the TS is not realistic.

Zoom In Zoom Out Reset image size

For the case where the initial proton distribution is a superposition of three Maxwell–Boltzmann distributions representing core SW ions, PUIs transmitted across the TS, and PUIs reflected at the TS, the results are similar to the kappa distribution. By slightly varying D₀, Figure 7 shows that it is possible to choose a different D₀ to better fit either the f_κ or f_3M spectra to IBEX data, but in general they produce similar results. Moreover, the fact that both f_κ and f_3M can be used to reproduce IBEX data, but a filled-shell downstream distribution for PUIs (f_MS) cannot, suggests that if PUIs are in a filled-shell distribution upstream of the TS, they most likely experience diffusion in energy via interactions with electromagnetic fluctuations at the ion kinetic scale at the TS, producing a smoother energy distribution than that assumed by Equation (8). We test this hypothesis in the next section.

In the previous sections, we tested different analytic forms for the proton distribution function downstream of the TS, f_p,0. However, based on the results in Section 4.2, if the PUI distribution upstream of the TS is a filled-shell, just heating the distribution at the TS according to Equation (8) is not realistic. Thus, we are motivated to test whether a more realistic derivation of f_p,0 from an upstream filled-shell distribution can clarify how PUIs interact at the TS. To do this, we utilize the results of a fully kinetic PIC simulation of the TS from Kumar et al. (2018). This simulation utilizes data taken by Voyager 2 during its crossing of the TS in 2007 (Burlaga et al. 2008; Richardson et al. 2008) to constrain the upstream PUI distribution function. Kumar et al. (2018) showed that in order to reproduce Voyager 2 observations, including the compression ratio, energy partitioning, and magnetic field profile, the upstream PUI distribution can be represented by a filled-shell distribution, but with a cutoff approximately twice that of the SW flow speed just upstream of the TS (flow speed was ∼320 km s⁻¹), such that the PUI mean energy is higher. They suggest that this is due to the increase in magnetic field observed by Voyager 2 within ∼1 au of the TS, resulting in the adiabatic compression and heating of PUIs.

While the properties of the TS at the location of Voyager 2 are not necessarily the same as in the downwind direction, and time variations in the SW make exact comparisons difficult, testing the results of Kumar et al. (2018) may offer a better prediction for the processing of PUIs at the TS, compared to the analytic functions described previously. We extract the proton distribution function downstream of the TS simulated by Kumar et al. (2018), binning all protons in speed bins used in our transport model (see the Appendix). The speed distributions that we test as boundary conditions in our model are shown in Figure 8, normalized such that the total density is 0.0017 cm⁻³. Figure 8(a) shows the downstream proton distribution in the case where the upstream PUI distribution was a filled-shell with cutoff at 320 km s⁻¹ in the SW frame, i.e., the SW speed upstream of the TS at Voyager 2. Figure 8(b) is similar, except the upstream PUI filled-shell cutoff was 640 km s⁻¹, which Kumar et al. (2018) suggested as a more appropriate form for the PUI distribution at the time and location of Voyager 2's crossing of the TS due to its higher mean energy. In the speed range ∼100–800 km s⁻¹, corresponding to the IBEX energy range, the proton spectra are approximately power laws with speed slopes between −4 and −3.5. The distribution in Figure 8(b) has a shallower slope since more PUIs upstream of the TS exist at higher energies, and these PUIs are able to gain more energy at the TS (Kumar et al. 2018).

Zoom In Zoom Out Reset image size

We take the distributions shown in Figure 8 as f_p,0 in our proton transport model, solve Equation (1), and calculate ENA fluxes at 1 au. We assume that the diffusion coefficient is proportional to v^1.25, for consistency with our previous results. We also show the case when there is no diffusion in velocity. With no diffusion, the resulting ENA spectra for both initial distributions shown in Figure 8 are too low compared to IBEX data over most energies; interestingly, however, the 640 km s⁻¹ upstream cutoff case matches the IBEX data above ∼2 keV even without velocity diffusion (Figure 9(b)). We find that the 320 km s⁻¹ upstream cutoff case is consistent with IBEX data if D₀ ∼ 8 × 10⁻⁹ km² s⁻³ (Figure 9(a)), similar to the results presented in Section 4.1, but the other case cannot match IBEX data for any diffusion coefficient since it produces too many high energy ENAs. While we do not show the results for different diffusion spectral indices, it is clear that the PIC simulation results with an upstream filled-shell cutoff of 320 km s⁻¹ can produce results consistent with our prior conclusion that the diffusion coefficient in the heliotail scales between v^2/3 and v².

Zoom In Zoom Out Reset image size

The fact that the downstream proton distribution that best matches the data from Voyager 2's observations at the TS fails to produce ENAs observed by IBEX from the heliotail may seem counterintuitive. However, this may indicate that (1) the long-term average of the proton distribution upstream of the TS should be a filled-shell with a cutoff at the SW speed, and that Voyager 2 observations represent an example of temporal variations about this average, or (2) the upstream PUI distribution as assumed by Kumar et al. (2018) is not strictly a filled-shell with higher speed cutoff, but possibly a filled-shell with cutoff at the upstream SW speed and a suprathermal, but relatively steep, tail at higher energies produced by particle interactions with turbulence farther upstream of the TS. This last possibility, as noted by Kumar et al. (2018), would still need to have a mean energy similar to that of a filled-shell with cutoff at 640 km s⁻¹ in order to reproduce Voyager 2 observations.

In this study, we seek to understand how the functional form of the proton distribution processed by the TS in the downwind direction of the heliosphere (i.e., heliotail) affects ENA fluxes observed by IBEX at 1 au. By constraining the initial distribution function downstream of the TS with IBEX data, we may better understand how the TS processes PUIs, as well as understand the relative importance of energy-dependent processes occurring in the heliotail (charge exchange, adiabatic heating, and stochastic acceleration). We model IBEX ENA observations by solving the stationary Parker transport equation (Equation (1)) for the proton distribution in the heliotail assuming different functional forms for the proton distribution downstream of the TS. We include effects from charge exchange and adiabatic heating based on results from an MHD-plasma/kinetic-neutral simulation of the heliosphere, as well as a velocity diffusion term as a free parameter.

First, we test three analytic forms of the proton distribution function (f_κ, f_MS, f_3M, see Table 1). We find that, while differences in the distributions are apparent close to the TS, at the speeds coinciding with the IBEX sensor's energy range (∼100–1000 km s⁻¹) the distributions approach similar shapes after traveling a few hundred astronomical units from the TS due to the velocity diffusion term. As pointed out by Zirnstein et al. (2018b), a kappa distribution with index κ_p,0 = 1.63 does not significantly change with distance in the IBEX energy range since the best-fit diffusion rate approximately counteracts the loss of energetic protons by charge exchange. The other distributions, which are not initially smooth like the kappa distribution, change significantly over distance before smoothing in velocity space.

If the proton distribution downstream of the TS is a kappa distribution, we find that practically any initial kappa index in the range 1.51 ≤ κ_p,0 ≤ 10, where κ_p,0 = 10 is close to a Maxwell–Boltzmann distribution, is possible for a diffusion coefficient scaling close to ∼v^1.3, with amplitude D₀ ∼ 0.8–1.3 × 10⁻⁸ km² s⁻³. While it appears that D(v) scaling close to v^1.3 best reproduces IBEX data, a higher (D(v) ∝ v²) or lower (D(v) ∝ v^2/3) diffusion rate may be possible within our uncertainties. A change from the widely known isotropic Kolmogorov (D(v) ∝ v^2/3; e.g., Isenberg 1987) and anisotropic Goldreich–Sridhar (D(v) ∝ v²; e.g., Goldreich & Sridhar 1995; Eichler 2014) turbulence spectra suggests that the presence of PUIs, which constitute a significant fraction of the plasma density far from the TS, may be modifying the turbulence spectrum. We note also that particle interactions with compressive waves may also result in D(v) ∝ v² (see e.g., Toptyghin 1980; Chalov et al. 2003; le Roux et al. 2005; Fahr et al. 2012; Zhang & Lee 2013).

We are not able to constrain the kappa index of the proton distribution function downstream of the TS due to the degeneracy between the kappa index and diffusion rate. However, a kappa index of less than 2 appears to be most probable based on earlier analyses of IBEX ENA observations (Livadiotis et al. 2011), as well as from Voyager observations of 0.04–4 MeV protons (Decker et al. 2005), and models of kappa distribution transport (e.g., Fahr et al. 2014, 2016). If this is the case, the average velocity diffusion amplitude in the heliotail likely lies in the range D₀ ∼ 1–1.3 × 10⁻⁸ km² s⁻³ for a diffusion rate that scales as D(v) ∝ v^1.3.

A proton distribution constructed from a superposition of three Maxwell–Boltzmann distributions, commonly used to model multicomponent PUI distributions in the heliosheath (e.g., Zank et al. 2010), produces similar results to a kappa distribution. This is due to their similar initial shapes and processing by velocity diffusion. We also test a PUI filled-shell distribution that is compressed at the TS, producing a downstream filled-shell distribution with cutoff at large speed as proposed by Zank et al. (2010). We show that this distribution cannot reproduce IBEX observations, creating too many ENAs at ∼2–5 keV regardless of the diffusion rate in the heliotail. This is not a surprising result since more self-consistent models and simulations of particle interactions at quasi-perpendicular shocks show that PUIs experience diffusion in energy at the TS (e.g., Yang et al. 2015; Kumar et al. 2018), resulting in a smoother roll-over or tail rather than a sharp cutoff at higher energies. In fact, this was observed for the first time by the Solar Wind Around Pluto (SWAP) instrument on board New Horizons at an interplanetary shock ∼34 au from the Sun (Zirnstein et al. 2018c).

We test a more realistic processing of the PUI filled-shell distribution at the TS using the results from a fully kinetic PIC simulation from Kumar et al. (2018). First, we find that the PIC simulation produces a downstream distribution that, similar to the kappa and three Maxwell–Boltzmann distributions, cannot produce enough ENAs over most of the ENA spectrum without some form of velocity diffusion in the heliotail (see, however, Figure 9(b)). However, if the velocity diffusion coefficient in the heliotail does scale approximately as v^1.3, then the most likely diffusion amplitude is approximately 8 × 10⁻⁸ km² s⁻³ based on the PIC simulation results, which is similar to the results found earlier in this study (Sections 4.1–4.2) and by Zirnstein et al. (2018b). Interestingly, we find that a PUI filled-shell distribution with a large cutoff speed (here we tested a cutoff of 640 km s⁻¹), as inferred by Kumar et al. (2018) to be more realistic than a cutoff at the SW speed in order to reproduce Voyager 2 observations of the TS, is inconsistent with IBEX observations of ENAs from the heliotail. This suggests either (1) that the long-term average of the proton distribution upstream of the TS could be a filled-shell with cutoff at the local SW speed and that Voyager 2 observations represent an instance of a temporal variation about this average, or (2) the upstream PUI distribution, as assumed by Kumar et al. (2018), is not strictly a filled-shell with higher speed cutoff, but possibly a filled-shell with cutoff at the upstream SW speed plus a suprathermal, but relatively steep, tail at higher energies.

Thus, our results suggest that the PUI distribution upstream of the TS can be a filled-shell like that suggested by Vasyliunas & Siscoe (1976), in agreement with recent observations from New Horizons' Solar Wind Around Pluto (SWAP; McComas et al. 2008) instrument, which showed that even out to 40 au from the Sun the PUI distributions can be quantified using the filled-shell function (though not necessarily representing physically meaningful parameters; McComas et al. 2017b). Moreover, IBEX ENA measurements substantiate recent SWAP observations of the preferential heating of PUIs at an interplanetary shock and the production of a suprathermal PUI tail downstream of the shock (Zirnstein et al. 2018c).