Modeling the Time Dependent Interaction between the Solar Wind and the Local Interstellar Medium with Kinetic Neutral Hydrogen and Helium Atoms

We present recent advancements in our 3D modeling of the interaction between the solar wind and the local interstellar medium (LISM). The latest model results (Fraternale et al., ApJ, 2023) have raised a question about the electron density of the LISM near the heliopause. We have shown that the presence of helium ions leads to a significant underestimation of this parameter compared to the past simulations and Voyager 1 PWS observations. The latter observations, with over 12 years’ worth of LISM data, offers a robust constraint on our models. Here we present additional simulations in support of the idea that the LISM proton density may need to be revised from approximately 0.054 cm–3 to values around 0.07 cm–3 or higher. Additionally, we have developed and successfully tested a new version of the kinetic code suitable for simulating time-dependent solutions.


Introduction
The interaction between the solar wind (SW) and the local interstellar medium (LISM) creates the heliosphere, serving as a natural laboratory for the investigation of various physical processes that occur under plasma conditions that are impossible to replicate in a terrestrial setting.Only a handful of missions provide data from the remote solar wind (SW) and the local interstellar medium (LISM), such as the still operational Voyager 1 and 2 (V1, V2), the Interstellar Boundary EXplorer (IBEX), and the New Horizons (NH) missions.
The LISM is influenced by the presence of the heliosphere in several ways.It extends over made by NH, have been focused on reevaluating the estimates of LISM parameters, e.g., see [6][7][8].A possibility also exists that the LISM in the solar neighborhood is composed of a mixture of the Local Interstellar Cloud (LIC) and the G-Cloud [9].
We have recently developed a fully 3D MHD-plasma/kinetic-neutrals global model, which, for the first time includes both neutral hydrogen (H) and helium (He) atoms self-consistently [4].The most recent version also includes electrons treated as a separate, charge-neutralizing fluid [10,11].
We have shown that the inclusion of He in a SW-LISM interaction model requires substantial reassessments of the H and H + distributions in the LISM with important consequences for the global heliosphere.A suggestion was made to revise also the LISM proton density to values higher than 0.054 cm −3 [10].In addition, the heliospheric community accumulated enough data (22 years of Ulysses and IBEX observations combined, as well as 11 years of V1 sampling the VLISM in situ) to investigate the solar cycle effects on the neutral atom transport.
In this paper, we further elaborate on the discussion presented in [10] and provide updated versions of the model also suitable for time-dependent simulations.

Physical model and boundary conditions
Our physical model includes an MHD description of plasma and a kinetic description for neutral H and He atoms.The most recent version is described in [10,11].The plasma description is based on the solution of a set of the conservation laws for the mixture of thermal protons, pickup protons , electrons, and helium ions, hereinafter referred to as the plasma mixture: Here ρ s , p s , T s , are the mass density, thermal pressure and temperature of each plasma component s, u and B are the bulk velocity and magnetic field vectors, respectively, and u and B are their magnitudes.The total energy density is E = ρ(ϵ + u 2 /2) + B 2 /8π, where ϵ is the specific internal energy (ρϵ = p/(γ − 1)) and the adiabatic index γ is equal to 5/3 for all species.The terms S ρ,m,E,p s represent the sources in the mass, momentum, energy, and pressure equations due to charge exchange and photoionization, for each species s.These terms are computed via two Monte Carlo kinetic solvers [12].In the electron pressure equation, the S p e source term includes photoionization only.The source terms Q C s represent the energy transfer rate to species s due to Coulomb collisions with the other species, as described in [10].
As regards helium ions, here we introduce a modified version of the model presented in [10].Equations ( 5)-( 6) are now applied to two difference species in two different environments: He 2+ ions (alpha particles) in the SW and He + ions in the LISM.In fact, in the SW at 1 AU the density of He 2+ ions is higher by approximately a factor of 1000 than that of He + ions.On the other hand, the charge-exchange cross sections σ(He + He 2+ ) at ∼ 5 keV are only ∼ 3.5 times smaller than σ(He + He + ).Therefore, here we neglect He + in the SW and include He 2+ ions instead.The latter carry a significant amount of ram pressure (8% to 20% of that of SW protons) which contributes to the pressure balance at the HP.The opposite is true in the VLISM, where He + ions carry about 40% of the mass of all ions, while He 2+ ions are comparatively negligible.Therefore, in Eqs. ( 5)-( 6), z = 2+ in the SW regions and z = + in the LISM.This modification is the first step toward the inclusion of both He 2+ and He + everywhere in our model.Here, we have taken into account the ionization rates of He due to photoionization and charge exchange.For the latter process, for alpha particles, only the reaction He 2+ + He = He + He 2+ is included.We fit the the helium cross-section data published by Barnett, et al. [13], obtaining the following expression σ cx (He + He 2+ ) = (3.142− 0.490 log E) 0.967 × 10 −16 (cm 2 ), (9) where E is the collision energy expressed in keV.The above expression fits the Barnett, et al. data with the mean relative error equal to 2.59% in the energy range E < 100 keV, and the maximum relative error less than 5.5%.For hydrogen, we have been using the cross sections from Lindsay & Stebbings (LS05) [14], However, due to significant uncertainties in the available data, especially at low energies, no consensus has been reached so far and various formulas have been proposed (see the discussion in [15]).In this study, see Section 3.1, we provide a comparison of the results obtained using Eq.(10) with those obtained using the updated cross-sections provided by [16,17]: The above expression is expected to provide more accurate cross-section values in the lowenergy range (0.01 eV -100 eV), while LS05 may be more accurate at high energies [17].In the new approach, Eq. ( 11) is employed at E < 1 keV while LS05 is used for E ≥ 1 keV.A different fit based on Barnett et al. data [13] was proposed in [15], where the effects on global modeling were also discussed.
Our fully 3D, hybrid-parallelized kinetic code builds on the previous hydrogen-only versions [18][19][20][21].The Boltzmann equations for the transport of neutral atoms are solved via a Monte Carlo method [12,19,22] using two independent modules for H and He atoms, respectively.The solution is obtained through an iterative process consisting of two stages: the fluid stage and the kinetic stage.In the fluid stage, a fixed distribution of source terms is utilized, which is obtained from the neutral modules.On the other hand, during the kinetic stage, the distribution of plasma is fixed, and kinetic particles are run.Each stage has its own time interval that depends on whether we seek a steady-state solution or conduct a time-dependent simulation.
For a steady solution on a Cartesian grid with the base grid size of 10 AU cubed, we typically run the kinetic module every τ p = 0.3 years (about 8 time steps of the base grid for plasma) to avoid spurious instabilities or other numerical artifacts [10,22].To gather sufficient statistics for source terms, the kinetic modules run for τ N = 700 − 900 years (using ∼ 0.4 × 10 9 macroparticles per species).Both time intervals (τ p , τ N ) must be carefully chosen to avoid numerical artifacts.These artifacts may include, among others, (i) the generation of spurious wave-like modes in the heliosheath when τ p is too large or (ii) the enhancement of instabilities at the HP and increased stiffness in the solution when the source terms become too noisy due to low statistics, which occurs when τ N is too small.
In a time-dependent scenario, both the time interval of the kinetic modules and the time interval of the plasma need to be significantly reduced to accurately account for the spatio-temporal scales present in both plasma and neutral atom distributions, as discussed in Section 3.3.
The kinetic modules require the VDFs of ions to be specified.In the previous studies with helium, we assumed a Maxwellian plasma [4,10].In Section 3.1, we discuss the effects of employing Lorentzian (kappa) VDFs for protons in the SW regions.
We use the level-set method [23,24] to track the position of the HP.Basically, we solve the advection equation to track the position of a surface (or 'level'), ∂ t ψ(x, t) + u∇ψ = 0.The SW and LISM regions are initialized at t = 0 with the values ψ = −1 and ψ = 1, respectively (also used as boundary conditions), so that the HP is represented by the surface ψ = 0 separating the two plasmas [25], see the discussion in Section 3.3.Note that this approach is not equivalent to that used, for instance in [26], where the continuity equation for the concentration of a passive scalar is used.The difference is substantial in compressible or nonuniform media.
The boundary conditions used for the simulations presented in this study are summarized in Table 1.The bottom rows of the table list the simulations, labeled from A to D. Simulation A corresponds to the best of models presented in [10] (S3 in that study), where the plasma is assumed to be Maxwellian everywhere for the calculation of source terms, the LS05 cross sections are employed, and the proton density in the LISM (n p,LISM ) is equal to 0.054 cm −3 .Simulations B1 to D use a higher density of 0.068 cm −3 , the motivation for which is discussed in Section 3.1.Simulations B1 and B2 differ only in the expressions for the H charge-exchange cross sections, i.e., in B1 we use LS05 at all energies and in B2 we use Eq.11 at low energies.With simulation C, we test the effect of using a kappa distributed protons in the SW regions.Finally, in simulation D, we test a new time-dependent kinetic model by simulating a nominal solar cycle (11-year periodic changes in the latitudinal extent of the fast and slow wind boundary), see details in Section 3.3.n p,LISM = 0.054 cm −3 , was obtained with the aid of steady-state global simulations assuming, as a constraint, the fixed nucleon LISM density of 0.089 cm −3 , which yielded reasonable agreement between the observed heliopause (HP) position in Voyager 1 (V1) data (121.6AU) and the one in global simulations [27].In simulations without helium, the assumption n p,LISM = 0.089 cm −3 resulted in the accurate electron density reproduction (because n p = n e ), see [34].However, if one uses n p,LISM = 0.054 cm −3 or lower, the V1 data are not reproduced, even when helium ions are included [4].
Recent V1 data shown in Fig. 1(a) confirm that a radial gradient exists in the electron density, with values remaining as high as 0.15 cm −3 in 2022 [33].Simulation A (black curve) underestimates these values by approximately 36%.To mitigate this discrepancy in global simulations, several factors can be adjusted, including n p,LISM , n He + ,LISM , or both, as well as modifying the density of neutrals.
Based on the simulations presented in the Appendix of [10], varying the LISM H density from 0.15 cm −3 to 0.22 cm −3 does not produce significant changes in the peak of n e , which remains around 0.1 cm −3 in those simulations.On the other hand, increasing the He + density alone would introduce excessive ram pressure due to the mass ratio.Consequently, we have first explored the effects of increasing the LISM density of protons to ∼ 0.07 cm −3 (degrees of ionization χ H = 0.26, χ He = 0.368, and (n He /n H ) LISM = 0.052).This values are also compatible with the estimates provided in [35].
Including a realistic amount of alpha particles in the SW is sufficient to maintain the heliopause (HP) at a distance similar to previous simulations (compatible with Voyager observations).As illustrated in panel (a), the underestimation of n e is reduced in the new simulations B1 and B2, although the values still align with the lower bound of observations, indicating the need for further improvement.Importantly, it is worth noting that the comparison with NH data remains robust in these simulations.In all new simulations, as far as the plasma is concerned, the most evident changes as compared to case A are seen in the inner heliosheath, where the temperatures of the plasma mixture and protons (thermal + PUIs) immediately behind the TS becomes higher by ∼ 12.5% and ∼ 17.5%, respectively.Using kappa-distributed protons in the SW only (case C) does not lead to significant changes to the heliosheath thickness (compare the blue and green curves).The reason for this lies in the energy dependence of charge exchange cross sections, as explained in [36][37][38], and recently discussed in [39].In fact, the energy transfer due to charge exchange, associated with the energetic tails of the Lorentzian distributions, is reduced due to smaller cross sections at high collision energy.However, notable changes in the H atoms distributions are obtained, as discussed in Section 3.2.
Using kappa-distributed protons in the SW (case C) does not produce considerable changes in the positions of the TS and HP (compare the blue and the green curves).The reason for this resides in the energy dependence of the charge-exchange cross sections [36][37][38].In fact, the energy transfer due to charge exchange associated with the energetic tails of the Lorentzian distributions, is mitigated due to the smaller cross sections at high energy.In contrast, the properties of neutral H atoms vary considerably, see for instance the density distributions in Fig. 2.

Distributions of neutral hydrogen density and the effect of updated cross sections for hydrogen charge exchange
Figure 2 shows the distribution of the total density of H atoms (panel a) and the density of population 2 H atoms (those produced in the IHS, panel b), along with a comparison with the density derived in [7] based on NH data.The comparison between simulations B1 and B2 reveals that using the updated cross sections does not result in substantial differences in the plasma distributions and standoff distances of the HP and TS.However, as shown in Fig. 2, the density of H atoms becomes higher when the updated cross sections are applied.This increase varies with distance, ranging from 9.5% at 2 AU to 5.5% at the TS.Consequently, the H densities in simulations B2 and C are close to that in simulation A. The density of population 2 H atoms that reach 2 AU in simulation C decreases by 57-67%, compared to case A, however these neutrals are hotter.

Time dependent simulations with kinetic neutral hydrogen and helium atoms
Here we present the a preliminary time-dependent simulation (case D) with both kinetic neutral H and He atoms.As mentioned in Section 2, we would like to reduce the iteration time of the kinetic neutrals step to about 1 year, while keeping the same statistics, in order to resolve solar-cycle variability and large-scale fluctuations that may arise both in the plasma and in the neutral atom distributions.Because it is computationally unfeasible to increase the number of macroparticles ∼ 10 3 times, a different strategy needs to be developed.We achieved the objective by (1) creating, for each particle (He and H atoms), many realizations of the same iteration interval, slightly randomizing the initial position and speed of the particles and eventually averaging the collected source terms.Note that the occurrence of a charge exchange event may vary in terms of both timing and spatial location, for the different realizations of one individual particle; (2) we tune the parameters of particle splitting in front of the HP to increase the actual number of macroparticles in the SW regions by a factor of ∼ 3.
As far as step ( 1) is concerned, a similar approach was proposed in [40] and used in [41].After collecting the source terms, only the set of newly born particles resulting from the first realization is retained, and all particles are advanced by one year.Careful consideration is given to the randomization of velocities and positions.Specifically, one must avoid (i) smoothing out gradients or inhomogeneities in the neutral atom distributions that may be present in the time-dependent scenario and (ii) excessively impacting the neutrals' VDF.To ensure a proper position randomization, we establish maximum displacements that increases linearly with the heliocentric distance, denoted as δ = min(|x N | 20 850 , δ max ) AU, with δ max equal to 20 AU in the LISM and 5 AU in the SW regions.For each stochastic realization j, we calculate the new particle position x N,j = x N,0 + a • δ (where a ∈ [−1, 1] is a random number drawn from a uniform distribution).Regarding velocities, randomization only modifies the velocity vectors by mere 1%, exclusively within the LISM.
The same randomization strategy is applied in step when particles undergo the new splitting (2).The new splitting process occurs in two steps.The first step closely resembles the original version described by [42].In this step, macroparticles are enabled to split when they reach predetermined specified locations in the physical domain.However, the actual splitting, i.e. the generation of a new particle with appropriately weighted mass, only takes place when a charge exchange event occurs.In the case of helium, where the mean free paths for charge exchange are large, this strategy may lead to fewer particles than desired.The second step involves a 'boost' that ensures the generation of a new physical particle, only for those particles for which splitting was enabled in the first step but in the absence of a charge exchange event.In this case, after splitting, the new position of a particle and its velocity are slightly randomized following the criteria mentioned above.
In simulation D, the domain is 1680 AU cubed, with the Sun centered at (x, y, z) = (1000, 850.850) AU.The plasma grid is Cartesian, has 5 levels of refinement, and the base grid size of ∆x = ∆y = ∆z = 10 AU.The best resolution is ∼ 0.3 AU in all directions, applied everywhere around the Sun at R < 20 AU, and in the upwind region (x > 0) extending into the LISM radially by ∼ 50 AU.In the distant tail, the resolution is 1.25 AU cubed everywhere.The source terms from the kinetic modules are collected on parallelepiped-shaped Cartesian kinetic grids, as described in [20], with the finest level size of 1.25 AU.
Figure 3 presents 2D distributions of electron density, plasma pressure, and speed in the B−V plane (defined by the velocity and magnetic field vectors in the unperturbed LISM, intersecting the Sun), and in the cross sections of the heliotail perpendicular to the x-axis at 900 AU from the Sun, from simulation D.
For the sake of simplicity in code testing, the heliospheric magnetic field (HMF) at the inner boundary is maintained unipolar and does not undergo polarity reversal every 11 years.A similar simulation was performed in [43] to demonstrate that even under the unipolar HMF assumption, the solar cycle influences drastically the structure of the tail because it affects the SW collimation within the two polar lobes, as compared with the steady state simulations with spherically-symmetric SW.Moreover, as discussed in [10], fluctuations and turbulence in the tail are also different from those obtained in simulations with steady and uniform b.c.'s.Even in this unrealistic scenario, the tail does not exhibit a two-fold split configuration, as shown by the level set contours in the right panels of Fig. 3.The HP position is formally represented by the contour ψ = 0.In this simulation, the level set is intentionally not reinitialized and is allowed to diffuse numerically to highlight the presence of a diffusion/mixing layer.In the right panels of Fig. 3, we display the contours at ψ = −0.9(in green) and ψ = +0.9.It is important to note that (i) the presence of the mixing layer results from a combination of both numerical diffusion and actual mixing due to large-scale motions or instabilities at the HP and (ii) the level set can provide only partial information about the mixing.Moreover, it is known that in the simulations with bipolar HMF, solutions frequently exhibit small HMF values near the HP, which also favours the growth of instabilities [44].We further clarify here that the concept of a comet-shaped heliotail does not exclude the presence of mixing of the SW and LISM plasma as discussed in [45,46].

Conclusions
The most recent MHD-plasma/Kinetic-neutrals model implemented in MS-FLUKSS incorporates both hydrogen and helium atoms, treated self-consistently, and features separate fluid electrons and helium ions.This represents a state-of-the-art model for investigating global heliospheric processes and enables us to investigate the transport of neutral atoms throughout the heliosphere.Such model is critical for enhancing theoretical models and interpreting spacecraft data, and will provide predictions to assist the interpretation of future observations by the Interstellar Mapping and Acceleration Probe (IMAP) mission [47].
In this study, we presented recent updates of the model published recently in [10], that include a simplified, preliminary treatment of alpha particles; updated hydrogen cross sections from [16]; and the effects of kappa-distributed plasma.The new models suggest the need for a revision of interstellar plasma parameters, particularly the proton density, to align with Voyager PWS observations of electron density in the OHS.Furthermore, we have developed and successfully tested a time-dependent version of the kinetic model.

Figure 1 .
Figure 1.Comparison of plasma distributions in steady-state simulations and spacecraft data.(a) Electron density in the VLISM and comparison with V1/PWS data; (b) Proton and helium ions density and comparison with NH data; (c) Protons thermal pressure and ram pressure; (d) Proton temperature and comparison with NH data (PUIs in green and all protons in black).

Figure 3 .
Figure 3. Two-dimensional distributions of electron density (top panels), plasma pressure (middle panels) and velocity magnitude (bottom panels) in a time-dependent simulation with kinetic H and He atoms.The left panels show quantities on the B ∞ − V ∞ plane, while the right panels show cross cuts of the heliotail at 900 AU from the Sun.In the right panels, two level set contours for ψ = −0.9 and ψ = +0.9 are also shown with green and blue lines, respectively.

Table 1 .
Boundary conditions and list of simulations presented in this study.LISM = 0.054 cm −3 , Maxw.plasma, LS05 cross sections (Eq.10),S3 in [10] B1 n p,LISM = 0.068 cm −3 , Maxw.plasma, LS05 cross sections (Eq.10)B2 n p,LISM = 0.068 cm −3 , Maxw.plasma, updated cross sections C n p,LISM = 0.068 cm −3 , kappa SW protons, updated cross sections D n p,LISM = 0.068 cm −3 , nominal solar-cycle with unipolar B3.Results and Discussion3.1.Electron density in the VLISM and kappa distributionsOne of the findings discussed in[10]is that the new models with helium underestimate the electron density observed by V1 in the VLISM.Figure1presents linear distributions of plasma [30][31][32][33]ned from steady-state simulations A -C.In Panel (a), the electron densities in the V1 data and from our models are compared.Voyager provides indirect measurements of electron density through PWS observations of plasma oscillations and quasi-thermal emission at the plasma frequency[30][31][32][33]. A commonly adopted value for the LISM proton density, 15th International Conference on Numerical Modeling of Space Plasma Flows Journal of Physics: Conference Series 2742 (2024) 012011