# Interstellar Neutral Hydrogen in the Heliosphere: New Horizons Observations in the Context of Models

P Swaczyna, M Bzowski, K Dialynas, L Dyke, F Fraternale, A Galli, J Heerikhuisen, M Z Kornbleuth, Dimitra Koutroumpa, I Kowalska-Leszczyńska, et al.

► To cite this version:

P Swaczyna, M Bzowski, K Dialynas, L Dyke, F Fraternale, et al. Interstellar Neutral Hydrogen in the Heliosphere: New Horizons Observations in the Context of Models. The Astrophysical Journal Letters, 2024, 969, pp.L20. 10.3847/2041-8213/ad5832. insu-04633667

HAL Id: insu-04633667

https://insu.hal.science/insu-04633667v1

Submitted on 3 Jul 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Interstellar Neutral Hydrogen in the Heliosphere: New Horizons Observations in the Context of Models

P. Swaczyna<fnref n="1" />, M. Bzowski<fnref n="1" />, K. Dialynas<fnref n="2" />, L. Dyke<fnref n="3" />, F. Fraternale<fnref n="4" />, A. Galli<fnref n="5" />, J. Heerikhuisen<fnref n="6" />,
M. Z. Kornbleuth<fnref n="7" />, D. Koutroumpa<fnref n="8" />, I. Kowalska-Leszczyńska<fnref n="1" />, M. A. Kubiak<fnref n="1" />, A. T. Michael<fnref n="9" />, H.-R. Müller<fnref n="3" />,
M. Opher<fnref n="7" />, and F. Rahmanifard<fnref n="10" />

Received 2024 April 15; revised 2024 June 5; accepted 2024 June 12; published 2024 July 1

## Abstract

Interstellar neutral (ISN) hydrogen is the most abundant species in the outer heliosheath and the very local interstellar medium (VLISM). Charge-exchange collisions in the outer heliosheath result in filtration, reducing the ISN hydrogen density inside the heliosphere. Additionally, these atoms are intensively ionized close to the Sun, resulting in a substantial reduction of their density within a few astronomical units from the Sun. The products of this ionization—pickup ions (PUIs)—are detected by charged particle detectors. The Solar Wind Around Pluto instrument on New Horizons provides, for the first time, PUI observations from the distant heliosphere. We analyze the observations collected between 22 and 52 au from the Sun to find the ISN hydrogen density profile and compare the results with predictions from global heliosphere models. We conclude that the density profile derived from the observations is inconsistent with steady-state model predictions. This discrepancy is not explained by time variations close to the Sun and thus may be related to the temporal evolution of the outer boundaries or VLISM conditions. Furthermore, we show that the cold and hot models of ISN hydrogen distribution are not a good approximation closer to the termination shock. Therefore, we recommend a new fiduciary point based on the available New Horizons observations at 40 au from the Sun, at ecliptic direction (285°, 62, 1.°, 94), where the ISN hydrogen density is 0.11 cm⁻<fnref n="3" />. The continued operation of New Horizons should give better insight into the source of the discussed discrepancy.

*Unified Astronomy Thesaurus concepts: Interstellar atomic gas (833); Interstellar medium (847); Heliosphere (711); Interstellar medium wind (848); Interplanetary physics (827); Pickup ions (1239); Solar wind termination (1535); Collision physics (2065); Charge transfer (2218)*

## 1. Introduction

Hydrogen atoms are the most abundant species in the very local interstellar medium (VLISM) around the heliosphere (Frisch et al. 2011). Because the Sun moves relative to the VLISM, interstellar neutral (ISN) atoms flow into the heliosphere (Axford et al. 1963). The ISN atoms are coupled to the magnetized plasma through charge-exchange collisions (Ruciński et al. 1996). The mean free path between charge-exchange collisions in the VLISM is of the order of ~100 au, i.e., comparable with the size of the heliosphere. Consequently, while the collisions are not frequent enough to equilibrate the conditions between the plasma and neutrals, they impact the density distribution of ISN hydrogen around and in the heliosphere. Moreover, the ISN hydrogen population in the heliosphere is slower and warmer than the pristine VLISM conditions (Ripken & Fahr 1983; Lallement & Bertaux 1990).

The plasma density is significantly lower in the heliosphere than in the VLISM and outer heliosheath, making charge exchange less frequent (Gurnett & Kurth 2019; Richardson et al. 2019). The ionized ISN atoms have a distinctive energy in the plasma frame and are called pickup ions (PUIs; Möbius et al. 1985). Furthermore, neutralized energetic ions, such as PUIs, form energetic neutral atoms (ENAs). Even though ENAs contribute to the neutral hydrogen density in the heliosphere, their density is orders of magnitude lower than that of the ISN hydrogen (Bzowski & Galli 2019). PUIs are also created through photoionization and electron impact ionization (Bzowski et al. 2013). The region, inside which the ISN density is reduced by these ionization processes more than 1/e, is called the ionization cavity. The cavity evolves over the solar cycle, extending to ~4–5 au from the Sun in the upwind direction and more than 20 au on the downwind side (Ruciński & Bzowski 1995; Sokól et al. 2019b).

The filtration in the outer heliosphere and filtration close to the Sun require different modeling approaches. The first one is typically addressed by global modeling of the heliosphere using a kinetic model of the neutral component, self-consistently coupled with a magnetohydrodynamic (MHD)

<figure><img src="image_3.png" /></figure>

Original content from this work may be used under the terms of the [Creative Commons Attribution 4.0 licence](https://creativecommons.org/licenses/by/4.0/). Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.

<footnote n="1">¹Space Research Centre PAS (CBK PAN), Bartycka 18a, 00-716 Warsaw, Poland; pswaczyna@cbk.waw.pl</footnote>

<footnote n="2">²Center for Space Research & Technology, Academy of Athens, 10679 Athens, Greece</footnote>

<footnote n="3">³Department of Physics and Astronomy, Dartmouth College, Hanover, NH 03755, USA</footnote>

<footnote n="4">⁴Center for Space Plasma and Aeronomic Research, The University of Alabama in Huntsville, Huntsville, AL 35899, USA</footnote>

<footnote n="5">⁵Space Research and Planetary Sciences, Physics Institute, University of Bern, 3012 Bern, Switzerland</footnote>

<footnote n="6">⁶Department of Mathematics and Statistics, University of Waikato, Hamilton 3240, New Zealand</footnote>

<footnote n="7">⁷Astronomy Department, Boston University, Boston, MA 02215, USA</footnote>

<footnote n="8">⁸LATMOS-IPSL, CNRS, UVSQ Paris-Saclay, Sorbonne Université, 78280 Guyancourt, France</footnote>

<footnote n="9">⁹Johns Hopkins University Applied Physics Laboratory, Laurel, MD 20723, USA</footnote>

<footnote n="10">¹⁰Physics Department, Space Science Center, University of New Hampshire, Durham, NH 03824, USA</footnote>
description for plasma (e.g., Malama 1991; Baranov & Malama 1993; Heerikhuisen et al. 2006). The global approach to modeling of the ISN hydrogen typically does not provide good statistics close to the Sun. Moreover, many of these models do not consider details of the solar cycle evolution of ionization and radiation pressure due to computational limitations. To overcome this problem, the modeling in the proximity of the Sun uses time- and heliolatitude-dependent models evolved from the cold and hot models of the ISN gas in the heliosphere (Fahr 1968; Thomas 1978; Ruciński & Bzowski 1995; Tarnopolski & Bzowski 2009) based on boundary conditions at the termination shock.

The observations performed close to the Sun, i.e., within the ionization cavity, require detailed modeling of the ionization processes, including the temporal and latitudinal variations (Sokół et al. 2019a), and thus are prone to significant uncertainties. Therefore, observations from larger distances are crucial to derive the ISN hydrogen properties. For example, the density can be derived from observations of PUIs in the outer heliosphere (Zirnstein et al. 2022). The PUI observations from Ulysses from distances of ~5 au allowed Bzowski et al. (2008) to derive the ISN hydrogen density at the termination-shock density of ~0.09 cm⁻³. This value was consistent with estimations using other methods (Pryor et al. 2008; Richardson et al. 2008; Bzowski et al. 2009). However, Dialynas et al. (2019) utilized a unique combination of 28–3500 keV in situ ions from Voyager 2 and 5.2–55 keV ENAs from INCA on Cassini and predicted that the ISN hydrogen density should be ~30% higher (~0.12 cm⁻³) beyond the termination shock. Moreover, Swaczyna et al. (2020) found that the density is higher by ~40% based on the New Horizons observations performed up to ~38 au from the Sun, i.e., far beyond the ionization cavity. This Letter provides an update on the derivation of the ISN hydrogen density using New Horizons PUI observations up to ~52 au from the Sun (McComas et al. 2021, 2022) and a comparison of the results with the predictions of the global heliosphere models.

## 2. PUI Observations from New Horizons

NASA's New Horizons mission was primarily intended to study the Pluto system and the Kuiper Belt (Stern 2008). The spacecraft was placed on a trajectory escaping the solar system using Jupiter's gravity assist. Solar Wind Around Pluto (SWAP), a plasma instrument of the mission, was designed to measure the rarefied solar wind (SW) plasma at large distances from the Sun (McComas et al. 2008). SWAP collects ions entering the detector from a wide (276° × 10°) field of view in narrow (ΔE/E = 0.085 FWHM) energy-per-charge bins from 0.035 to 7.5 keV q⁻¹. This energy range covers SW protons, α particles, and PUIs created in the supersonic SW from ISN hydrogen atoms. In the early stage of the mission, the SWAP was on only during short intervals (McComas et al. 2010; Randol et al. 2012, 2013). Since 2012, SWAP operates most of the time, providing almost continuous observations starting at ~22 au from the Sun. The top panel of Figure 1 shows the New Horizons trajectory projected on the ecliptic plane. The spacecraft remains close to the ecliptic plane. Between 22 and 52 au (the range of radial distances used in our analysis), the spacecraft's ecliptic longitude changed slowly from 280.°03 to 287.°32, while the ecliptic latitude changed from 1.°81 to 1.°98.

The SWAP PUI observations were analyzed in a series of papers. McComas et al. (2017) analyzed the PUI observations up to ~38 au by fitting the Vasyliunas & Siscoe (1976) distribution function. The methodology was later refined, accounting for SW He⁺ ions near the PUI cutoff energy (Swaczyna et al. 2019a) and allowing for nonadiabatic cooling of PUIs according to the Chen et al. (2014) model (Swaczyna et al. 2020). McComas et al. (2021) found that fitting the PUI distribution cutoff velocity may further improve the goodness of the fit and analyzed the SWAP observations out to nearly 47 au. In 2019 February, the SWAP flight software was updated to provide a higher cadence of PUI data, improving the time resolution from 1 day to ~30 minutes (McComas et al. 2022).

This Letter updates the analysis by Swaczyna et al. (2020) with additional data (McComas et al. 2021, 2022). Because McComas et al. (2021) reprocessed all SWAP observations, the individual data points presented here may not overlap with those shown in Swaczyna et al. (2020). Nevertheless, the impact of these changes is negligible. Furthermore, we aggregate the data collected with the ~30 minute resolution into 1 day averages to match the earlier resolution. We select only those days with at least 42 half-hour data blocks. We apply the selection criteria discussed in Swaczyna et al. (2020).

The ISN hydrogen density profile may be derived from the PUI observations in the supersonic SW using two methods. The first method relates the PUI density with the column density of ISN hydrogen from the Sun to the distance at which the observations are collected. The ionization rate decreases with the distance r from the Sun as $β(r) = β_0(r_0/r)^2$, where $β_0$ is the ionization rate at $r_0 = 1$ au. The radially expanding SW advects the PUIs; thus, their accumulated flux in the SW is proportional to the column density:

$$u_{sw}n_{PUI}(r) = \left(\frac{r_0}{r}\right)^2 \beta_0 N_{H}(r), \quad (1)$$

where $u_{sw}$ is the SW bulk speed and $N_{H}(r) = \int_{0}^{r} n_{H}(r')dr'$ is the column density of the parent ISN hydrogen from the Sun to the distance r (see green line in Figure 1, top panel). As discussed in detail in Swaczyna et al. (2020), we include corrections for the slowdown of the SW and the change in the proton density due to the PUI production through photoionization. The most statistically significant SWAP data points are two or three points immediately below the energy cutoff, dominated by freshly ionized PUIs. Therefore, the fit parameters may be used to derive the local density of ISN hydrogen. We closely follow the previous implementation, and the densities obtained from each day of the SWAP observations are presented in bottom panels of Figure 1. The individual data points present a significant scatter in the estimated densities caused by the SW variations. Therefore, we average them into 1 au wide bins, each covering approximately 4 months of observations. The averaged densities are shown as black dots with error bars.

Swaczyna et al. (2020) found that the observed densities up to ~38 au are consistent with the cold model for the ISN hydrogen density of 0.127 cm⁻³ far from the Sun. The density of the ISN hydrogen in the cold model depends on the distance from the Sun r and angular distance from the upwind direction
<figure><img src="image_4.png" /></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 1.** Top panel: the trajectory of New Horizons (red line) projected on the ecliptic plane shown on top of the density modulation of ISN hydrogen calculated with the WTPM (Tarnopolski & Bzowski 2009; Sokół et al. 2015). The solid portion represents the almost continuous observations between 22 and 52 au used in our study. The orange and magenta ellipses show the Jupiter and Pluto trajectories, respectively. The black arrow indicates the ISN hydrogen flow direction. The green line represents the line along which the PUIs accumulated for observations at 45 au from the Sun. Bottom panels: column density (left panel) and local density (right panel) of ISN hydrogen in the heliosphere based on the New Horizons/SWAP observations of PUIs. The color of the data points corresponding to each day of SWAP observations represents the cooling index value (see Swaczyna et al. 2020). The averages over 1 au wide bins are shown using points with error bars (gray—excluded from the model fits; black—included in the model fits; see text). The cold model fit from Swaczyna et al. (2020) is shown with dashed gray lines. The new fit using a cold model and a model with a gradient (see text) are shown using purple and orange lines, respectively.</figcaption></figure>

$\theta$ (Thomas 1978):

$$n_H(r, \theta) = n_C \exp\left(-\frac{\lambda}{r} \frac{\theta}{\sin \theta}\right), \quad (2)$$

where $n_C$ is the density of ISN hydrogen far from the Sun and $\lambda$ is the size of the ionization cavity along the upwind direction. We adopt the direction of the ISN hydrogen inflow (252°2, 9°0) in ecliptic coordinates obtained from SWAN observations (Lallement et al. 2005). The adopted cavity size is consistent with the one used in the fitting of the PUI distribution function $\lambda = 4.0 \pm 0.5$ au (Sokół et al. 2019b; McComas et al. 2021). In fitting the models, we use only the observations collected from 22 au from the Sun, i.e., over the period of continuous SWAP observations. We removed the sporadic observations collected earlier because they represent shorter averaging times and

thus may be biased by the SW variability. The best-fit densities far from the Sun $n_C$ are $0.1278 \pm 0.0043_{(\lambda)} \pm 0.0011_{(\text{stat})} \text{ cm}^{-3}$ and $0.1217 \pm 0.0017_{(\lambda)} \pm 0.0010_{(\text{stat})} \text{ cm}^{-3}$ for the method using the column density and local density, respectively. As in Swaczyna et al. (2020), we adopt the mean of the results and the discrepancy between these two results as a measure of methodological uncertainty. Therefore, the result is $n_C = 0.1248 \pm 0.0030_{(\text{method})} \pm 0.0033_{(\lambda)} \pm 0.0010_{(\text{stat})} \text{ cm}^{-3}$. This result is well within the uncertainty range reported by Swaczyna et al. (2020), which should be expected considering that the data set used here includes the data set from their paper. Swaczyna et al. (2020) discussed that the subtraction of He$^+$ ion contribution adds a 1.2% relative uncertainty ($\pm 0.0015_{(\text{He}^+)} \text{ cm}^{-3}$), the adoption of the Lindsay & Stebbings (2005) charge-exchange cross section adds a 10% relative uncertainty ($\pm 0.013_{(\sigma)} \text{ cm}^{-3}$), and the
instrumental uncertainty is ~4% (±0.005(instr) cm⁻³). Therefore, the combined uncertainty is ±0.015 cm⁻³.

Due to the systematic nature of the dominant uncertainties, our analysis does not reduce the combined uncertainty compared to the previous study. Nevertheless, the larger span in the distances from the Sun provides a means to study departures of the measured density profile from the density profile predicted by the cold model. In Section 3, we show that, up to ~90 au from the Sun, the density profile can be described in all three considered models using the following ad hoc model with a gradient along the ISN hydrogen flow direction:

$$n_H(r, \theta) = n_G(1 + \gamma r \cos \theta) \exp \left( -\frac{\lambda}{r} \frac{\theta}{\sin \theta} \right), \quad (3)$$

where $\gamma$ is a parameter describing the gradient slope along the inflow direction. Because the best-fit parameters obtained from separate fits to the column and local densities are strongly correlated, we perform a combined fit using the densities from both approaches. The best-fit parameters are $n_G = 0.138 \pm 0.008(\lambda) \pm 0.007(\text{stat}) \text{ cm}^{-3}$ and $\gamma = -0.0036 \pm 0.0010(\lambda) \pm 0.0011(\text{stat}) \text{ au}^{-1}$. Both these uncertainty components are almost perfectly anticorrelated, with a correlation coefficient of -0.98. Nevertheless, the fits to the cold model and the model with a gradient give $\chi^2 = 96.6$ and 64.2, respectively, using 62 data points. The other uncertainty sources mentioned above also impact the parameter $n_G$ because they represent systematic effects that scale the densities. Combining all these sources of uncertainty yields a total uncertainty of ±0.019 cm⁻³.

### 3. ISN Hydrogen Density in Global Heliosphere Models

In this Letter, we use the ISN hydrogen densities obtained in three state-of-the-art hybrid kinetic-MHD global heliosphere models to put the New Horizons data in the context of global simulations. Model A has been developed over the years by the group at the University of Alabama in Huntsville (Pogorelov et al. 2008; Fraternale et al. 2021, 2023). Here, we use the most recent version, which incorporates hydrogen and helium atoms and ions treated self-consistently (see Fraternale et al. 2024). Model B is the product of the SHIELD Drive Science Center (Opher et al. 2015, 2017, 2020, 2023; Michael et al. 2021). Finally, Model C uses the same approach to the MHD solution as Model A but with different grid and boundary conditions (DeStefano & Heerikhuisen 2017; Heerikhuisen et al. 2019). Models A and C represent a traditional view of the heliosphere exhibiting a long tail in the downwind direction, while Model B challenges this view with a croissant-shaped heliosphere with the SW confined in two downwind jets. Table 1 summarizes the boundary conditions and selected assumptions employed in the models.

Figure 2 presents the ISN hydrogen density profiles along New Horizons trajectory in the global heliosphere models. All profiles have the highest densities in the outer heliosheath corresponding to the so-called hydrogen wall (Baranov & Malama 1993). The models, which use steady-state conditions, are not optimized to reproduce the density of ISN hydrogen close to the Sun. Therefore, we do not expect them to reproduce the solution close to the Sun (see Section 4).

We fit the cold model formula using the model predictions between 22 and 52 au, i.e., for the same range as in the fit to New Horizons data (Section 2), but we assume a fixed cavity size parameter $\lambda = 4$ au. The best-fit (asymptotic) ISN

hydrogen densities far from the Sun $n_C$ are 0.119, 0.132, and 0.128 cm⁻³ for models A, B, and C, respectively (purple lines in Figure 2). All models show a discrepancy in the radial gradient of the ISN hydrogen density. Therefore, we also fit the model given by Equation (3) to the global modeling results between 22 and 90 au. The additional term is the simplest possible extension of the cold model. The factor $\cos\theta$ is motivated by the global distribution of ISN hydrogen, suggesting a large-scale gradient along the inflow direction (see, e.g., Heerikhuisen et al. 2008, Figure 1). The fitting procedure resulted in the following best-fit parameters: $n_G = 0.112, 0.130$, and $0.121 \text{ cm}^{-3}$, and $\gamma = 0.0018, 0.0004$, and $0.0018 \text{ au}^{-1}$ for models A, B, and C, respectively.

The slope parameters found are positive, unlike the fit to the New Horizons observations. These fits provide a good description of the density profile inside the termination shock beyond the ionization cavity. However, the discrepancies between the numerical models and the cold model extrapolations near 90 au show that the extrapolation of the ISN hydrogen density to the termination shock using the cold model to “tune” the global models is incorrect. Such an approach was commonly used (e.g., Zirnstein et al. 2016; Bzowski & Heerikhuisen 2020).

Using the ISN hydrogen density profiles predicted by the global models, we find scaling factors for each model that best fit the observations using $\chi^2$ minimization with the observational data. The obtained values are 1.0184, 0.9277, and 0.9484 for models A, B, and C, respectively. Figure 3 compares the scaled density profiles with the observational data, the cold model fit, and the model with gradient fit found in Section 2. Note that the model with gradient fit to the data has an opposite gradient, and thus, the extrapolation beyond the current range of New Horizons data differs significantly from the fits shown in Figure 2. Applying these factors to scale the VLISM ISN hydrogen density assumed as the boundary conditions, the corrected densities are 0.202, 0.168, and $0.163 \text{ cm}^{-3}$ for models A, B, and C, respectively. The mean value is ~0.18 cm⁻³, i.e., about 10% less than the estimate of $0.195 \text{ cm}^{-3}$ from Swaczyna et al. (2020), which was based on the mean filtration factor provided by Müller et al. (2008). Considering nonlinear effects, e.g., through coupling with the plasma, the best-fit densities may differ slightly from these estimates.

The density in the VLISM cannot be constrained directly. Therefore, we recommend a fiducial point along the New Horizons trajectory at 40 au from the Sun, where the ISN hydrogen density is $0.11 \text{ cm}^{-3}$. The ecliptic coordinates are (285.°62, 1.°94). At this distance, the fits to the observed densities of all models considered in this study, including the analytic formulae, converge within less than 1%. This distance is far enough from the Sun that the consequences of the processes within the ionization cavity are negligible.

### 4. Discussion

Our analysis shows that the cold model reasonably approximates the ISN hydrogen density profile along the New Horizons trajectory. However, the model with gradient, given by Equation (3), informed by the density profiles in the global heliosphere models, provides a statistically better fit. Interestingly, the best-fit slope of this additional gradient found from the fit to the New Horizons data has the opposite sign and is 2 times higher than the slopes from the global models.
<table><caption>**Table 1**
Boundary Conditions and Assumption in the Global Models</caption><thead><tr><th></th><th>Model A</th><th>Model B</th><th>Model C</th></tr></thead><tbody><tr><td>Model Description (reference)</td><td>Fraternale et al. (2021, 2023, 2024)</td><td>Opher et al. (2015); Michael et al. (2021)</td><td>Heerikhuisen & Pogorelov (2010); Heerikhuisen et al. (2019)</td></tr><tr><td>Simulation Box</td><td>(1680 au)<sup>3</sup>, upwind boundary 680 au from Sun</td><td>(3000 au)<sup>3</sup>, upwind boundary 1500 au from Sun</td><td>Sphere with a 1000 au radius around Sun</td></tr><tr><td>Simulation Grid Resolution</td><td>Adaptive<sup>a</sup> (10 au)<sup>3</sup> → (0.3 au)<sup>3</sup></td><td>Adaptive<sup>a,b</sup> (6–12 au)<sup>3</sup> → (0.7 au)<sup>3</sup></td><td>Adaptive<sup>a</sup> r: 10 au → 0.4 au; Constant θ: 1.°<br>5, φ: 3°</td></tr><tr><td colspan="4">LISM Boundary Conditions Assumed at the Outer Boundary</td></tr><tr><td>LISM Flow Speed</td><td>25.4 km s<sup>-1</sup></td><td>26.4 km s<sup>-1</sup></td><td>25.4 km s<sup>-1</sup></td></tr><tr><td>LISM flow direction (in ecliptic coord.)</td><td>(75.°7, -5.°1)</td><td>(75.°4, -5.°2)</td><td>(75.°7, -5.°1)</td></tr><tr><td>LISM Temperature</td><td>7500 K</td><td>6519 K</td><td>7500 K</td></tr><tr><td>LISM H<sup>0</sup> Density</td><td>0.195 cm<sup>-3</sup></td><td>0.18 cm<sup>-3</sup></td><td>0.17 cm<sup>-3</sup></td></tr><tr><td>LISM He<sup>0</sup> Density</td><td>0.0154 cm<sup>-3</sup></td><td>...</td><td>0.0153 cm<sup>-3</sup></td></tr><tr><td>LISM H<sup>+</sup> Density</td><td>0.068 cm<sup>-3</sup></td><td>0.06 cm<sup>-3</sup></td><td>0.0494 cm<sup>-3</sup></td></tr><tr><td>LISM He<sup>+</sup> Density</td><td>0.00898 cm<sup>-3</sup></td><td>...</td><td>0.0089 cm<sup>-3</sup></td></tr><tr><td>LISM Magnetic Field</td><td>2.93 μG</td><td>4.4 μG</td><td>2.93 μG</td></tr><tr><td>LISM B-V Angle</td><td>39.°5</td><td>15.°1</td><td>39.°5</td></tr><tr><td>LISM B-V Plane Inclination to Ecliptic</td><td>52.°2</td><td>61.°4</td><td>52.°2</td></tr><tr><td colspan="4">SW Boundary Conditions</td></tr><tr><td>Location of Inner Boundary</td><td>10 au<sup>c</sup></td><td>30 au<sup>d</sup></td><td>10 au<sup>c</sup></td></tr><tr><td>SW Speed</td><td>420 km s<sup>-1</sup> at 1 au</td><td>417 km s<sup>-1</sup> at 30 au</td><td>411 km s<sup>-1</sup> at 1 au</td></tr><tr><td>SW Density</td><td>8 cm<sup>-3</sup> at 1 au</td><td>0.00874 cm<sup>-3</sup> at 30 au</td><td>6.74 cm<sup>-3</sup> at 1 au</td></tr><tr><td>SW Composition</td><td>H<sup>+</sup> + 3.5%<sup>e</sup> He<sup>2+</sup></td><td>Only H<sup>+</sup></td><td>H<sup>+</sup> + 4%<sup>e</sup> He<sup>2+</sup></td></tr><tr><td>SW Temperature</td><td>80,000 K at 1 au</td><td>108,680 K at 30 au<sup>f</sup></td><td>51,100 K at 1 au</td></tr><tr><td>SW Magnetic Field</td><td>Parker spiral B<sub>r</sub> = 37.5 μG at 1 au</td><td>Parker spiral B<sub>r</sub> = 65.5 μG at 1 au</td><td>Parker spiral B<sub>r</sub> = 37.5 μG at 1 au</td></tr><tr><td colspan="4">Model Assumptions</td></tr><tr><td>Plasma-neutral Coupling</td><td>Charge exchange + Photoionization</td><td>Charge exchange</td><td>Charge exchange + Photoionization</td></tr><tr><td>He Ions and Atoms</td><td>He<sup>2+</sup> in SW, He<sup>+</sup> in VLISM included self-consistently</td><td>...</td><td>He<sup>2+</sup> in SW, He<sup>+</sup> in VLISM comoving and thermalized with H<sup>+</sup></td></tr><tr><td>p-H<sup>0</sup> Charge Exchange</td><td>&lt;1 keV: Swaczyna et al. (2019b); &gt;1 keV: Lindsay & Stebbings (2005)</td><td>Lindsay & Stebbings (2005)</td><td>&lt;1 keV: Bzowski & Heerikhuisen (2020); &gt;1 keV: Lindsay & Stebbings (2005)</td></tr></tbody></table>

**Notes.**

<sup>a</sup> Adaptive resolution with larger cells near the outer boundary and smaller cells near the inner boundary.

<sup>b</sup> Grid for kinetic neutral component is coarser; see Michael et al. (2021) for details.

<sup>c</sup> Analytic plasma solution within 10 au from the Sun is employed; no inner boundary for neutral atoms.

<sup>d</sup> Conditions advectioned from 1 au using the Izmodenov et al. (2005) model.

<sup>e</sup> By number.

<sup>f</sup> Temperature of plasma with PUIs.

To check if the discrepancy in the slope is connected with the temporal evolution of the ionization and radiation pressure, we compare the results of two time-dependent hot models. The first one is the Warsaw Test Particle Model (WTPM; Tarnopolski & Bzowski 2009; Sokół et al. 2015) implementing state-of-the-art models of the temporal evolution of the ionization processes (Sokół et al. 2020; Porowski et al. 2022) and radiation pressure (Kowalska-Leszczynska et al. 2018a, 2018b, 2020) along the trajectories of modeled atoms. The other is the LATMOS hot model (Lallement et al. 2005; Quémerais et al. 2006), which considers stationary ionization conditions applied with a year-to-year variation of the total ionization rate and radiation pressure based on the SWAN latitude-dependent data (Koutroumpa et al. 2019). Figure 4 compares these models with the data. We proportionally scaled both hot models to match the densities at the termination shock

with the values obtained in our analysis. The figure shows a perfect agreement of the WTPM with the cold model. The LATMOS model predicts stronger variation as the annual ionization values are applied simultaneously along the atoms' trajectories. Therefore, we conclude that the simplifications of the cold model formula and the temporal variations of the ionization rates and radiation pressure are not the source of the discrepancy for the distance range covered in this study.

The global models predict positive slopes with the magnitude depending on the detail of the modeling, and the local models cannot explain the negative slope with time-depending ionization processes near the Sun. Consequently, this phenomenon requires a different explanation. The ISN hydrogen density may be related to temporal changes in the filtration at the outer heliospheric boundaries in response to the solar cycle and secular modulation of the SW (McComas et al.
<figure><img src="image_6.png" /><figcaption>**Figure 2.** ISN hydrogen density along New Horizons trajectory for the three models (in columns). The numerical results are presented using the green line. The cold model and the model with a gradient fitting the models over the limited range are shown using purple and orange lines, respectively. The bottom panels magnify the portion within 90 au from the Sun for easier comparison.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Figure 3.** ISN hydrogen density profiles (colored lines) from the fit analytic models and scaled numerical models used in this study compared to the observations (dots with error bars). The bottom panel shows the density profiles within the range covered by the New Horizons observations. The fiducial point distance and density are marked with dashed lines (see text).</figcaption></figure>

2019). These changes influence the position of the heliopause by a few astronomical units throughout the solar cycle (e.g., Washimi et al. 2017; Izmodenov & Alexashov 2020; Sokół et al. 2021) and, thus the filtration of ISN hydrogen in the heliosphere's outer boundaries. Another possibility is that the enhancement in the density is related to an inhomogeneity of the ISN hydrogen density in the VLISM.

## 5. Summary

New Horizons is the first spacecraft to provide PUI observations from the distant outer heliosphere, far beyond the ionization cavity. Consequently, the derived ISN hydrogen density profile is less affected by complex processes impacting the ISN hydrogen distribution within a few astronomical units from the Sun, such as ionization and radiation pressure. We find that the extended data set of New Horizons observations follows the predictions of the cold model with the density $n_C = 0.125 \text{ cm}^{-3}$ far from the Sun. However, the model with an additional density gradient along the inflow direction provides a statistically better fit. Nevertheless, we found that the gradient coefficient is negative, i.e., the density decreases as the spacecraft moves toward the VLISM inflow direction. This finding is inconsistent with the global heliosphere models that predict an opposite sign of such gradient within the region of the supersonic SW.

We employ three global heliosphere models to put the New Horizons observations in the perspective of the ISN hydrogen density profile through the heliosphere's outer boundaries. We find scaling factors for each model that best fit the observations. The scaled profiles and the fits of the analytic models give very similar ISN hydrogen densities, $0.11 \text{ cm}^{-3}$ (±1%) at 40 au from the Sun, i.e., approximately midway through the used New Horizons observations. We recommend this as a fiduciary point for future comparison with the predictions of the ISN
<figure><img src="image_8.png" /><figcaption>**Figure 4.** The ISN hydrogen density profiles predicted by the hot models and analytic formulas. The left and right panels show the column and local ISN hydrogen densities, respectively, as a function of distance from the Sun along the New Horizons trajectory. The hot models are scaled by a constant factor to match the data. The hot models do not explain the negative slope.</figcaption></figure>

hydrogen density models. The density at the termination shock is not a good fiduciary point because while the cold model prediction is close to the asymptotic value, the global models show strong gradients at these distances. The scaled ISN hydrogen density profiles from the global models using different boundary conditions estimate the ISN hydrogen density in the pristine VLISM between 0.163 and $0.202 \text{ cm}^{-3}$.

The apparent discrepancy between the global model predictions and the slope observed on New Horizons may result from transient effects, the filtration in the outer boundary regions of the heliosphere or a result of large-scale gradients in the pristine VLISM. Future observations of the SWAP instrument on New Horizons should allow further investigation of this problem.

## Acknowledgments

This research was supported by the International Space Science Institute (ISSI) in Bern through ISSI International Team project #541 (Distribution of Interstellar Neutral Hydrogen in the Sun's Neighborhood). P.S. is supported by a project cofinanced by the Polish National Agency for Academic Exchange within the Polish Returns Programme (BPN/PPO/2022/1/00017) with a research component funded by the National Science Centre, Poland (2023/02/1/ST9/00004). M.B. and M.A.K. are supported by the National Science Centre, Poland (2019/35/B/ST9/01241). K.D. acknowledges support at JHU/APL by NASA under contracts NAS597271, NNX07AAJ69G, and NNN06AA01C and by subcontract at the Center for Space Research and Technology (Academy of Athens) and useful discussions within the SHIELD DRIVE Science Center: https://shielddrivecenter.com/. L.D. and H.R.M. are supported by NASA grant 80NSSC21K1681. F.F. is supported by NASA grants 80NSSC24K0267 and 80NSSC22K0524. F.F. also acknowledges the Texas Advanced Computing Center (TACC) at The University of Texas at Austin for providing HPC resources on Frontera supported by NSF award CISE-OAC-2031611. M.Z. K. and M.O. are supported by NASA grant 18-DRIVE18_0029, Our Heliospheric Shield, 80NSSC22M0164, and also by NASA HGI grant 80NSSC22K0525. D.K. is supported by CNES through its Sun Heliosphere Magnetosphere (SHM) program. I.K.L. is cofunded from the statutory funds of CBK PAN (own research fund TT-420). F.R. is

funded by the IBEX mission as part of NASA's Explorer Program (80NSSC18K0237). The authors, except the corresponding author, are listed in alphabetical order. For the purpose of Open Access, the authors have applied a CC-BY public copyright license to any Author Accepted Manuscript (AAM) version arising from this submission.

## Data Availability

The SWAP derived data products used in this study were obtained by McComas et al. (2021, 2022) and can be downloaded from: https://spacephysics.princeton.edu/missions-instruments/swap/pui-data-2021 and https://spacephysics.princeton.edu/missions-instruments/swap/pui-data-2022.

## ORCID iDs

P. Swaczyna https://orcid.org/0000-0002-9033-0809

M. Bzowski https://orcid.org/0000-0003-3957-2359

K. Dialynas https://orcid.org/0000-0002-5231-7929

L. Dyke https://orcid.org/0009-0007-0108-2318

F. Fraternale https://orcid.org/0000-0002-4700-2762

A. Galli https://orcid.org/0000-0003-2425-3793

J. Heerikhuisen https://orcid.org/0000-0001-7867-3633

M. Z. Kornbleuth https://orcid.org/0000-0002-3479-1766

D. Koutroumpa https://orcid.org/0000-0002-5716-3412

I. Kowalska-Leszczyńska https://orcid.org/0000-0002-6569-3800

M. A. Kubiak https://orcid.org/0000-0002-5204-9645

A. T. Michael https://orcid.org/0000-0003-2227-1242

H.-R. Müller https://orcid.org/0000-0001-7364-5377

M. Opher https://orcid.org/0000-0002-8767-8273

F. Rahmanifard https://orcid.org/0000-0001-9316-0553

## References

Axford, W. I., Dessler, A. J., & Gottlieb, B. 1963, ApJ, 137, 1268

Baranov, V. B., & Malama, Y. G. 1993, JGR, 98, 15157

Bzowski, M., & Galli, A. 2019, ApJ, 870, 58

Bzowski, M., & Heerikhuisen, J. 2020, ApJ, 888, 24

Bzowski, M., Möbius, E., Tarnopolski, S., Izmodenov, V., & Gloeckler, G. 2008, A&A, 491, 7

Bzowski, M., Möbius, E., Tarnopolski, S., Izmodenov, V., & Gloeckler, G. 2009, SSRv, 143, 177

Bzowski, M., Sokół, J. M., Tokumaru, M., et al. 2013, in Cross-calibration of Far UV Spectra of Solar System Objects and the Heliosphere, ed.
E. Quémerais, M. Snow, & R.-M. Bonnet (New York: Springer Science + Business Media), 67

Chen, J. H., Bochsler, P., Möbius, E., & Gloeckler, G. 2014, JGRA, 119, 7142

DeStefano, A. M., & Heerikhuisen, J. 2017, JPhCS, 837, 012013

Dialynas, K., Krimigis, S. M., Decker, R. B., & Mitchell, D. G. 2019, GeoRL, 46, 7911

Fahr, H. J. 1968, Ap&SS, 2, 474

Fraternale, F., Pogorelov, N. V., & Bera, R. K. 2023, ApJ, 946, 97

Fraternale, F., Pogorelov, N. V., & Bera, R. K. 2024, JPhCS, 2742, 012011

Fraternale, F., Pogorelov, N. V., & Heerikhuisen, J. 2021, ApJL, 921, L24

Frisch, P. C., Redfield, S., & Slavin, J. D. 2011, ARA&A, 49, 237

Gurnett, D. A., & Kurth, W. S. 2019, NatAs, 3, 1024

Heerikhuisen, J., Florinski, V., & Zank, G. P. 2006, JGRA, 111, A06110

Heerikhuisen, J., & Pogorelov, N. V. 2010, in ASP Conf. Ser. 429, Numerical Modeling of Space Plasma Flows, Astronum-2009, ed. N. V. Pogorelov, E. Audit, & G. P. Zank (San Francisco, CA: ASP), 227

Heerikhuisen, J., Pogorelov, N. V., Florinski, V., Zank, G. P., & le Roux, J. A. 2008, ApJ, 682, 679

Heerikhuisen, J., Zirnstein, E. J., Pogorelov, N. V., Zank, G. P., & Desai, M. 2019, ApJ, 874, 76

Izmodenov, V., Malama, Y., & Ruderman, M. S. 2005, A&A, 429, 1069

Izmodenov, V. V., & Alexashov, D. B. 2020, A&A, 633, L12

Koutroumpa, D., Quémerais, E., Ferron, S., & Schmidt, W. 2019, GeoRL, 46, 4114

Kowalska-Leszczynska, I., Bzowski, M., Kubiak, M. A., & Sokół, J. M. 2020, ApJS, 247, 62

Kowalska-Leszczynska, I., Bzowski, M., Sokół, J. M., & Kubiak, M. A. 2018a, ApJ, 852, 115

Kowalska-Leszczynska, I., Bzowski, M., Sokół, J. M., & Kubiak, M. A. 2018b, ApJ, 868, 49

Lallement, R., & Bertaux, J. L. 1990, A&A, 231, L3

Lallement, R., Quémerais, E., Bertaux, J. L., et al. 2005, Sci, 307, 1447

Lindsay, B. G., & Stebbings, R. F. 2005, JGRA, 110, A12213

Malama, Yu. G. 1991, Ap&SS, 176, 21

McComas, D., Allegrini, F., Bagenal, F., et al. 2008, SSRv, 140, 261

McComas, D. J., Dayeh, M. A., Funsten, H. O., et al. 2019, ApJ, 872, 127

McComas, D. J., Elliott, H. A., & Schwadron, N. A. 2010, JGRA, 115, A03102

McComas, D. J., Shrestha, B. L., Swaczyna, P., et al. 2022, ApJ, 934, 147

McComas, D. J., Swaczyna, P., Szalay, J. R., et al. 2021, ApJS, 254, 19

McComas, D. J., Zirnstein, E. J., Bzowski, M., et al. 2017, ApJS, 233, 8

Michael, A. T., Opher, M., Tóth, G., Tenishev, V., & Drake, J. F. 2021, ApJ, 906, 37

Möbius, E., Hovestadt, D., Klecker, B., et al. 1985, Natur, 318, 426

Müller, H.-R., Florinski, V., Heerikhuisen, J., et al. 2008, A&A, 491, 43

Opher, M., Drake, J. F., Swisdak, M., Zieger, B., & Toth, G. 2017, ApJL, 839, L12

Opher, M., Drake, J. F., Zieger, B., & Gombosi, T. I. 2015, ApJL, 800, L28

Opher, M., Loeb, A., Drake, J., & Toth, G. 2020, NatAs, 4, 675

Opher, M., Richardson, J., Zank, G., et al. 2023, FrASS, 10, 1143909

Pogorelov, N. V., Heerikhuisen, J., & Zank, G. P. 2008, ApJ, 675, L41

Porowski, C., Bzowski, M., & Tokumaru, M. 2022, ApJS, 259, 2

Pryor, W., Gangopadhyay, P., Sandel, B., et al. 2008, A&A, 491, 21

Quémerais, E., Lallement, R., Bertaux, J.-L., et al. 2006, A&A, 455, 1135

Randol, B. M., Elliott, H. A., Gosling, J. T., McComas, D. J., & Schwadron, N. A. 2012, ApJ, 755, 75

Randol, B. M., McComas, D. J., & Schwadron, N. A. 2013, ApJ, 768, 120

Richardson, J. D., Belcher, J. W., Garcia-Galindo, P., & Burlaga, L. F. 2019, NatAs, 3, 1019

Richardson, J. D., Liu, Y., Wang, C., & McComas, D. J. 2008, A&A, 491, 1

Ripken, H. W., & Fahr, H. J. 1983, A&A, 122, 181

Ruciński, D., & Bzowski, M. 1995, A&A, 296, 248

Ruciński, D., Cummings, A. C., Gloeckler, G., et al. 1996, SSRv, 78, 73

Sokól, J. M., Bzowski, M., & Tokumaru, M. 2019a, ApJ, 872, 57

Sokół, J. M., Dayeh, M. A., Fuselier, S. A., et al. 2021, ApJ, 922, 250

Sokół, J. M., Kubiak, M. A., Bzowski, M., & Swaczyna, P. 2015, ApJS, 220, 27

Sokól, J. M., Kubiak, M. A., & Bzowski, M. 2019b, ApJ, 879, 24

Sokól, J. M., McComas, D. J., Bzowski, M., & Tokumaru, M. 2020, ApJ, 897, 179

Stern, S. A. 2008, SSRv, 140, 3

Swaczyna, P., McComas, D. J., Zirnstein, E. J., & Heerikhuisen, J. 2019b, ApJ, 887, 223

Swaczyna, P., McComas, D. J., & Zirnstein, E. J. 2019a, ApJ, 875, 36

Swaczyna, P., McComas, D. J., Zirnstein, E. J., et al. 2020, ApJ, 903, 48

Tarnopolski, S., & Bzowski, M. 2009, A&A, 493, 207

Thomas, G. E. 1978, AREPS, 6, 173

Vasyliunas, V. M., & Siscoe, G. L. 1976, IGR, 81, 1247

Washimi, H., Tanaka, T., & Zank, G. P. 2017, ApJL, 846, L9

Zirnstein, E. J., Heerikhuisen, J., Funsten, H. O., et al. 2016, ApJL, 818, L18

Zirnstein, E. J., Möbius, E., Zhang, M., et al. 2022, SSRv, 218, 28