The Solar Wind with Hydrogen Ion Exchange and Large-scale Dynamics (SHIELD) Code: A Self-consistent Kinetic–Magnetohydrodynamic Model of the Outer Heliosphere

Neutral hydrogen has been shown to greatly impact the plasma flow in the heliosphere and the location of the heliospheric boundaries. We present the results of the Solar Wind with Hydrogen Ion Exchange and Large-scale Dynamics (SHIELD) model, a new, self-consistent, kinetic–MHD model of the outer heliosphere within the Space Weather Modeling Framework. The charge exchange mean free path is on the order of the size of the heliosphere; therefore, the neutral atoms cannot be described as a fluid. The numerical code SHIELD couples the MHD solution for a single plasma fluid to the kinetic solution for neutral hydrogen atoms streaming through the system. The kinetic code is based on the Adaptive Mesh Particle Simulator, a Monte Carlo method for solving the Boltzmann equation. The numerical code SHIELD accurately predicts the increased filtration of interstellar neutrals into the heliosphere. In order to verify the correct implementation within the model, we compare the results of the numerical code SHIELD to those of other, well-established kinetic–MHD models. The numerical code SHIELD matches the neutral hydrogen solution of these studies as well as the shift in all heliospheric boundaries closer to the Sun in comparison with the multi-fluid treatment of neutral hydrogen atoms. Overall the numerical code SHIELD shows excellent agreement with these models and is a significant improvement to the fluid treatment of interstellar hydrogen.


Introduction
Interstellar neutral hydrogen plays a major role in our understanding of the heliosphere. The heliosphere forms as the solar wind expands into the ionized portion of the local interstellar medium (ISM). The neutral atoms do not feel the boundaries separating the solar and interstellar plasma, and instead stream into the inner solar system. Hydrogen atoms are the dominant component of the ISM and are coupled to the plasma through resonant charge exchange (Wallis 1975). This removes up to 30% of momentum from the solar wind, decelerating the flow as it expands into the ISM. The influence of interstellar neutral hydrogen reduces the distance to the heliospheric boundaries by ∼150 au and significantly alters the structure of the termination shock (TS), the location where the solar wind becomes subsonic (Baranov & Malama 1993).
Neutral hydrogen also provides a mechanism to observe the global structure of the heliosphere. As a result of charge exchange, energetic neutral atoms (ENAs) are produced with the velocity of their parent ions. ENAs travel relatively unimpeded through the heliosphere and can be detected as far in as 1 au. ENA observations, such as those done by the Interstellar Boundary Explorer (IBEX; Schwadron et al. 2011Schwadron et al. , 2014 and the Ion and Neutral Camera instrument on Cassini (Krimigis et al. 2009), are able to indirectly probe the global structure of the heliosphere. Since neutrals play such a significant role in the system, including them self-consistently into models is crucial to understanding the global ENA maps and in situ measurements of the plasma at Voyager 1 and 2.
The solar wind-ISM interaction can be separated into four regions. The heliosphere consists of the supersonic solar wind in the inner heliosphere and the subsonic solar wind in the heliosheath (HS). If the Sun's motion through the ISM is supersonic, the pristine ISM will become subsonic at the bow shock (BS), before being diverted at the heliopause (HP). Baranov et al. (1981) was the first work to incorporate neutral hydrogen into computational models of the heliosphere. This model, as well as that of Pauls et al. (1995), treated the plasma and neutral hydrogen as two separate fluids, coupled by source terms that approximate the effect of charge exchange. Charge exchange in each region of the heliosphere will produce ENAs with distinct characteristics and charge exchange mean free paths. Each region has different characteristic parameters; therefore, a single-fluid description for neutrals is inadequate.
Subsequent models have employed a multi-fluid approximation for neutrals, defining a fluid population for each region of the heliosphere (Zank et al. 1996;Williams et al. 1997;Fahr et al. 2000;Pogorelov & Zank 2005;Opher et al. 2009). While the multi-fluid approach provides an incorrect solution for the neutral component, it is much more computationally efficient (McNutt et al. 1998). Furthermore,  were the first to show that the plasma solution differs only by a few percent in the upwind direction in comparison to a kinetic treatment of the neutrals. The fluid description assumes that there are sufficient H-H collisions to thermalize their distribution into a Maxwellian. However, H-H collisions are negligible throughout the heliosphere (Izmodenov et al. 2000), causing the distribution functions to be distinctly non-Maxwellian (Izmodenov et al. 2001). Additionally, since the charge exchange mean free path is on the order of the size of the heliosphere, the Knudsen number, Kn = l mfp /L, where l mfp is the mean free path and L is the characteristic length scale of the system-for neutral hydrogen it is around one. Therefore, a kinetic treatment is required to accurately model neutrals (Izmodenov et al. 2000).
The Monte Carlo method is one of the standard methods for numerically solving multidimensional kinetic equations. Bird (1994) has adapted it to model a rarefied gas whose dynamics are determined through a finite number of collisional processes such that the flow is unsteady. A main advantage of the Monte Carlo method is that it does not require the use of integrodifferential equations that describe the evolution of the distribution function. Instead, the Monte Carlo method follows the trajectories of many simulated macroparticles throughout the domain. Collisions are handled through particle pairs rather than by integrating the Boltzmann collision integral. This allows for a wider range of collisional processes to be included and makes the method a natural choice for outer heliosphere (OH) application.
The first Monte Carlo simulation of heliospheric neutrals was introduced by Malama (1991) and the first self-consistent kinetic treatment was implemented by Baranov & Malama (1993). In their model, Baranov & Malama (1993) solved the linear kinetic equation for the trajectory of the neutrals using a Monte Carlo approach. The Monte Carlo model was coupled to a single hydrodynamic plasma fluid, which was also done in later models by Müller et al. (2000) and Heerikhuisen et al. (2006). Subsequent models that use a kinetic solution for the neutrals have shown that it provides a more adequate description of the propagation of interstellar neutrals than the multi-fluid approach (Heerikhuisen et al. 2006;Müller et al. 2008;Alouani-Bibi et al. 2011). Charge exchange collisions are handled on a particle-byparticle basis while the source terms of the local plasma are the accumulation of events that occur in a particular location. The single plasma fluid has been extended to also include helium ions and alpha particles (Izmodenov et al. 2003) while anomalous cosmic rays have been treated as a massless, diffusive fluid by Alexashov et al. (2004). The Monte Carlo code can also be used to model additional interstellar atoms such as oxygen and nitrogen (Izmodenov et al. 2004), as well as interstellar dust (Alexashov et al. 2016).
The solar magnetic field has been shown to significantly impact the location of the heliospheric boundaries (Izmodenov & Alexashov 2015). Opher et al. (2015) also showed that the solar magnetic field can collimate the solar wind in the HS, possibly resulting in a two-lobe structure of the heliotail. An analytic solution for the collimation of the solar wind by the magnetic field was developed by Drake et al. (2015). They showed that the solar wind magnetic field strength affects the strength of the lobes as well as the thickness of the HS. Korolkov & Izmodenov (2021) have shown that stellar magnetic fields are also very important to the plasma flow structure in astrospheres around other stars. A self-consistent model of the heliosphere with kinetic neutrals coupled to a single magnetohydrodynamic (MHD) plasma fluid has been developed by Pogorelov et al. (2009) andAlexashov (2015). These kinetic-MHD models have been very successful at modeling the OH; however, they do have some limitations. Opher et al. (2017) showed that magnetic reconnection can occur in the eastern flank of the heliosphere, which organizes the magnetic field ahead of the HP. Izmodenov & Alexashov (2015) used a scheme that restricts any communication between the ISM and solar wind across the HP to limit nonideal-MHD effects and therefore cannot model this behavior. Additionally, Pogorelov et al. (2009) used a dipole description of the solar magnetic field that induces large numerical dissipation of the magnetic field and numerical reconnection across the HP. Since the scale of the current sheet is much smaller than the heliosphere, it is extremely difficult to resolve the structure. This results in the removal of magnetic field strength over large regions of the HS, making the results strongly affected by numerical artifacts (Michael et al. 2018).
In this work, we introduce the Solar Wind with Hydrogen Ion Exchange and Large-scale Dynamics (SHIELD) model, a global, self-consistent kinetic-MHD model of the heliosphere that utilizes the MHD solution of Opher et al. (2015). The first kinetic-MHD model was published by Izmodenov et al. (2005). The SHIELD code couples the OH and particle tracker (PT) components within the Space Weather Modeling Framework (SWMF; Tóth et al. 2005). We will first describe the codes and how they are coupled in a general sense before going into the specifics of the study used to verify their successful coupling and execution. The individual models are described in Section 2. The numerical scheme of the Monte Carlo code used in the PT component is presented in Section 2.2.1. Section 2.3 details how the Monte Carlo model is applied to the OH and coupled to the MHD solution. The study used to verify our model is presented in Section 2.5 and the results are presented in Section 3. Here, we adopt the same boundary conditions used in both  and Heerikhuisen et al. (2006) to compare our solution directly with their results. Finally, conclusions are presented in Section 4.

Model
The SHIELD code treats the plasma as a single fluid and couples the MHD solution with the kinetic solution of the neutrals to form a self-consistent model of the heliosphere. To this end, the numerical code SHIELD couples the OH and PT components within the SWMF because the SWMF has been designed for efficient code coupling. This greatly reduces the software development work as compared to an ad hoc coupling of the two models. While the SWMF has been developed for physics-based space weather modeling, it is a general and flexible software framework that incorporates physics models as components with minimal changes to their own optimal mathematical and numerical representations (Tóth et al. 2005). Components are easily coupled through the framework, allowing them to pass desired information back and forth. In addition to passing the physical quantities, the SWMF transforms between coordinate systems and interpolates among different numerical grids. This does not put any grid restrictions on the coupled components and allows each component to be run with an optimal mesh for that particular region. The SWMF compiles the desired models into a single executable and distributes the components onto a parallel computer, where they are executed simultaneously and coupled efficiently throughout the run.

The OH Component
The OH component is based on the Block-adaptive Tree Solar Wind Roe-type Upwind Scheme (BATS-R-US) solver, a highly parallel, 3D, block-adaptive, upwind finite-volume MHD code. BATS-R-US is the MHD code within the SWMF and is the core of the framework (Powell et al. 1999;Tóth et al. 2012). BATS-R-US can solve the ideal, resistive, and Hall MHD equations as well as describe semi-relativistic and multifluid plasma. BATS-R-US is a highly parallel, 3D code that can utilize adaptive mesh refinement (AMR) to produce a blockadaptive grid capable of running on Cartesian or generalized coordinate systems. The code is capable of explicit and implicit time stepping as well as solving the MHD equations with spatially first, second, and now fifth order accurate schemes (Chen et al. 2016). First-and second-order spatial schemes utilize a conservative, cell-centered, upwind, finite-volume method with a total variation diminishing (TVD) scheme to reduce numerical oscillations at shocks. The fifth-order spatially accurate scheme developed by Chen et al. (2016) is an MP5 limiter based, conservative, finite-difference scheme.
BATS-R-US was first adapted for the OH by Opher et al. (2003). Opher et al. (2003) instituted a 3D ideal-MHD, explicit, second-order solution that models the plasma component as a single fluid and neglects the interstellar magnetic field and neutrals. The OH model was further developed by Opher et al. (2004) to include both the solar and the interstellar magnetic field, making the model a full 3D ideal-MHD simulation. Opher et al. (2009) were the first to incorporate interstellar neutral hydrogen atoms into the model with the four-fluid approximation. When run in stand-alone mode, the OH component is a global 3D multi-fluid MHD simulation of the OH that describes the plasma and four neutral hydrogen species. The OH stand-alone model solves the ideal-MHD equations for the plasma and a separate set of Euler's equations for the different populations of neutral atoms. The neutral fluids are coupled to the plasma through source terms resulting from charge exchange as calculated by McNutt et al. (1998). Subsequent additions to the Opher et al. (2009) model include the sector structure of the heliospheric current sheet ) and time-and latitude-dependent solar cycle variations of the solar wind (Provornikova et al. 2014;Michael et al. 2015).
The OH component has been updated to include pickup ions (PUIs) as a separate fluid (Opher et al. 2020). Opher et al. (2020) solved the coupled, multi-fluid MHD equations for the thermal solar wind ions and PUIs as developed by Glocer et al. (2009) andTóth et al. (2012). The plasma components and neutral fluids were coupled through updated source terms depending on which fluids are undergoing charge exchange (Opher et al. 2020). In this paper, we chose to couple the kinetic model to the single plasma fluid approximation. This is in accordance with the kinetic-MHD models of the OH of Baranov & Malama (1993) and Heerikhuisen et al. (2006) and allows for a model verification before advancing to the multiion coupling. Future extensions of SHIELD will couple the kinetic code to the multi-ion model of Opher et al. (2020).

Multi-fluid Neutral Approximation
When run as a stand-alone code, the OH component is a global 3D multi-fluid MHD simulation of the OH that describes the plasma as a single fluid along with four neutral hydrogen species, one for each source region of the heliosphere.
The OH stand-alone model solves a set of Euler's equations for the different populations of neutral atoms and the ideal-MHD equations for the plasma. The conservative form of the equations is solved by an explicit, second-order spatiotemporal scheme, utilizing an HLLE-type Riemann solver developed by Linde (1998), for both the plasma and the neutral fluids. A TVD principle is employed to reduce oscillations at discontinuities. The OH component employs a monotized central limiter with a typical beta parameter set at 1.5 (Sokolov et al. 2006).
The impact of the charge exchange process on the plasma is approximated through source terms to the continuity, momentum, and energy equations of the ideal-MHD system: In the single-fluid approximation, where the newly formed PUI is immediately assimilated into the plasma, charge exchange does not alter the plasma's density. However, since the PUI has a different initial velocity from the original proton, the collision alters the momentum and energy of the plasma as the picked up ion is instantaneously accelerated to the plasma velocity. The change in momentum and energy due to charge exchange can be found by evaluating the moments of the collisional term (Equation (9)): The Maxwellian distribution of each fluid allows for analytic approximations for S ρv and S ò (Pauls et al. 1995;Williams et al. 1997;McNutt et al. 1998). This results in a fast, efficient calculation of the source terms, which are updated in every time step, and in quicker convergence of the solution for both steady-state and time-dependent problems. Analytic expressions for the source terms used in the OH stand-alone code were derived by McNutt et al. (1998), assuming both the neutral and ion distribution functions are Maxwellian. McNutt et al. (1998) Here i is the index for each neutral fluid species. The charge exchange frequency for each neutral fluid is given by , where σ ex is the charge exchange cross section taken from either Maher & Tinsley (1977) or Lindsay & Stebbings (2005). The bulk velocity and thermal speeds of the neutral and ion fluids undergoing charge exchange are U H(i) and U p and V th,H(i) and V th,p , respectively.
,H * ( ) are the relative velocities between the neutral and ion fluids used to evaluate the charge exchange cross section: To conserve energy and momentum, the source terms for the neutrals are equal and opposite to those in Equation (6).
The domain of the model is typically a 3000 au × 3000 au × 3000 au cube in a Sun-centered frame. The ISM enters into the simulation from the left-hand face (x = −1500 au). Our inner boundary is at 30 au for the ions, where the solar wind is taken to be either constant and uniform in both latitude and longitude or time-varying over the solar cycle (Provornikova et al. 2014;Michael et al. 2015). The time-varying solar wind conditions are extrapolated from 1 au to our inner boundary using a 1D model similar to that of Zubkov (2005), which takes into account the expansion of the solar wind and the charge exchange with interstellar hydrogen; however, we neglect the cosmic-ray component (Provornikova et al. 2014). For cases with uniform wind, we have derived our boundary conditions from similar models that take charge exchange into account, extracting values at 30 au to use them as solar wind conditions. We use a nonuniform Cartesian grid that can be adapted to ensure sufficient grid resolution anywhere in the heliosphere. The solar wind magnetic field is defined as a Parker spiral (Parker 1958) at the inner boundary. Outer boundary conditions are set such that the ISM enters into the domain from the x = −1500 au face, with outflow conditions, ensuring zero gradient of the fluid variables, at the others.
The external forces acting on the plasma and neutrals streaming into the heliosphere are the solar radiation pressure and gravity. At distances larger than where our inner boundary is applied, we assume that these forces are negligible and drop them from the governing equations. Additionally, the photoionization process is neglected within our calculation. Hot electrons in the HS could also contribute to the ionization of neutral atoms (Gruntman 2015); however, V2 observed that electrons have energies below the 10 eV instrument threshold (Richardson 2008). These electrons are not hot enough for electron impact ionization to have an influence on the mass loading of the solar wind in the HS; thus, this process has not been included within the OH component.

The Particle-tracking Component
When coupled in the SHIELD code, the PT component is applied to the OH and is used to kinetically solve for the neutral atom trajectories, including charge exchange collisions with the plasma. The PT component is based on the Adaptive Mesh Particle Simulator (AMPS) code, a global, kinetic, 3D particle code developed within the framework of the direct simulation Monte Carlo (DSMC) methods for the purpose of solving the Boltzmann equation for the motion and interaction of a dusty, partially ionized, multi-species gas in cometary comae (Tenishev et al. 2008). Among others, AMPS is capable of modeling collisions between species, photodissociation, energy exchange between internal vibrational and rotational states, and radiative cooling. AMPS employs a block-adaptive mesh that can utilize an AMR procedure to ensure that the cells are refined to resolve the local mean free path and capture the geometry of complex systems.
AMPS was developed to study the neutral components of cometary comae (Tenishev et al. 2008). The guiding principle that we followed when designing the code was to create a general-purpose tool for solving the Boltzmann equation for various space environments at conditions beyond those specified at the stage of designing the code. As a result, the range of the code's application spans from modeling the dynamics of neutrals, ions, and charged dust in planetary satellites and planetary exospheres on the lower-energy range (Valeille et al. 2009;Lee et al. 2014;Tenishev et al. 2014;Dong et al. 2018), up to modeling transport of solar energetic particles (Borovikov et al. 2017) and galactic cosmic rays on the other side of the energy range. A Monte Carlo method was also used by Heerikhuisen et al. (2006) and we use a similar method here to solve for the motion and distribution of interstellar neutral hydrogen.
AMPS is a fully functioning component of the SWMF, where it is coupled to the global MHD code BATS-R-US. In the work presented in this paper, AMPS is applied to model the neutral hydrogen population (both thermal atoms and ENAs) in the OH.

Numerical Scheme of AMPS
The numerical scheme in AMPS is described in detail by Tenishev et al. (2021); however, we provide a brief summary here. AMPS is a 3D DSMC model developed to solve the Boltzmann equation for a dusty, partially ionized, multi-species cometary coma. DSMC models are very effective at modeling collisional neutral fluids kinetically, which makes them a natural choice for modeling the neutrals streaming throughout the heliosphere. AMPS employs the DSMC method for solving the Boltzmann equation. The DSMC method was introduced by Bird (1994), and nowadays is one of the main numerical techniques used for solving the Boltzmann equation. The essence of the method is to separate the deterministic motion of the particles populating the simulated environment and affected by macroscopic forces from the stochastic interaction between those particles and the interaction of the particle with the ambient media. In the original application to the modeling of the dynamics of rarefied gas flows, the collision integral accounts only for binary collisions between molecules and atoms in the simulated gas flow. Mathematically, the collision integral is a term describing the source and loss of the distribution function, which in the application relevant to the work presented in the paper, accounts for the interaction of the neutral hydrogen population with the ambient plasma of solar wind via charge exchange.
The DSMC method employed by AMPS uses a series of Markov chains to solve for the evolution of the system (Tenishev et al. 2008(Tenishev et al. , 2021. AMPS accomplishes this by describing the neutral flow with a large but finite set of model particles. The number of interstellar neutrals traveling throughout the heliosphere is very large; therefore, a model particle within AMPS represents many real particles. The basic idea of the employed method is the separation of the translational motion of the neutrals through the heliosphere from their interaction with the plasma of the ambient solar wind. The latter is simulated probabilistically, calculating the time interval between successive charge exchange events in a particular computational cell according to the Poisson distribution. For modeling convection in the simulated environment, model particles are affected by the same physical laws as real neutrals in the heliosphere. That is equivalent to solving the kinetic equation with the method of characteristics. The model particles will follow the same trajectory as the real atoms would.
The velocities of the model particles injected into the computational domain are distributed according to the distribution function of real interstellar neutral atoms ahead of the heliosphere. The particles are assigned a weight that relates the number of real particles each model particle represents. Each neutral hydrogen particle is assigned the same weight. The weight is fixed, set when the model particles enter the domain, and it does not vary throughout the simulation. The weight is calculated to ensure that the number of injected model particles in each time step is equivalent to the number of real particles that enter the domain in that time period. This is derived from the neutral flux set by the boundary conditions. The velocity of each model particle is allocated according to a flux Maxwellian distribution, where the neutral hydrogen Maxwellian distribution is multiplied by the normal to the outer boundary. This ensures that the model particles reproduce the real distribution function of the atoms as they enter the simulation domain. The weight affects the statistical noise in the source terms the kinetic model supplies to the MHD simulation of the solar wind. The number of injected model particles can be varied to reach a level that provides a negligible level of statistical noise in the simulated neutral species' density profile.

Application to the OH
AMPS is applied to transport model hydrogen particles through the domain, tracking their trajectories through the heliosphere and determining the probability that a collision occurs for each model particle. If a collision does occur for a particle, the velocity is updated appropriately as well as the individual contribution to the plasma source terms due to the event. After the collisions, forces are applied to the particles and their locations are updated accordingly. The particles can be sampled at any time, producing the distribution function at any location. AMPS, generally, solves the time-dependent Boltzmann equation : for species s, with a distribution function, f s . Here F is an external macroscopic force and f t s c ¶ ¶ ( ) is the effect on the distribution function due to collisions as well as stochastic interactions with the ambient media.
The neutral composition of the ISM can be determined, as neutrals are unaffected by the heliospheric boundaries and penetrate deep into the heliosphere. These atoms, such as helium atoms that interact very weakly through charge exchange, can be measured directly by missions like Ulysses (Witte 2004) and IBEX (Möbius et al. 2012), or indirectly, as they are eventually ionized and picked up by the solar wind (Ruciński et al. 1996). Through these methods, H, He, N, O, Ne, and Ar have been detected. Neutral hydrogen constitutes 92% of the ISM by number. The rest of the ISM is comprised mostly of He atoms with N, O, Ne, and Ar comprising fractions of a percent (Gloeckler & Geiss 2001). Since helium atoms very rarely undergo charge exchange with ionized hydrogen, their impact on the structure of the heliosphere is minimal. Solar wind protons readily charge-exchange with oxygen atoms; however, the charge exchange cross section is an order of magnitude smaller than that of the H-H + process for energies typical of the solar wind (Lindsay & Stebbings 2005). Additionally, the density of oxygen is many orders of magnitude below the H density in the ISM, causing this reaction to be much more infrequent. We therefore model the neutral component of the ISM as a single hydrogen species gas and neglect all other species. Izmodenov et al. (2000) have shown that elastic hydrogenhydrogen or hydrogen-proton collisions are negligible over the spatial scale of the heliosphere when compared to the charge exchange mean free path. As a result, the only collisional process we model is charge exchange. The resulting collisional operator can be described in terms of the production and loss of neutral atoms: Here ΔE is the relative proton-neutral atom energy in keV calculated as m v 1 2 H rel 2 ( ) .
The external forces acting on the neutrals streaming into the heliosphere are the solar radiation pressure and gravity. At distances larger than where our inner boundary is applied, these forces are negligible. Consequently, there is no net external force acting on the atoms, allowing the last term on the lefthand side of Equation (8) to be dropped. This significantly reduces the complexity of the Boltzmann equation. When applied to the OH, AMPS therefore solves the resulting linear kinetic equation. The atoms travel on straight-line trajectories until a charge exchange event occurs and the velocity of the particles is altered.
The process of charge exchange between neutral species with the plasma of solar wind ions is simulated kinetically in AMPS considering the probability of a model particle participating in a charge exchange reaction according to the probability density of e 1 dt -t ( ) , where dt is the time step, and τ is the lifetime of the hydrogen atom at the location of the model particle. The lifetime is calculated individually for each model particle from the charge exchange loss rate of the particular particle, given by Since the plasma is described by a Maxwellian distribution function, Ripken & Fahr (1983) simplified the integral to where the average relative velocity over all protons, 〈v rel 〉, is given by A similar loss rate calculation was done by Malama (1991) followed by Heerikhuisen et al. (2006), with the lifetime given by τ = 1/ν ex . When a model particle representing the hydrogen atoms is selected for charge exchange, the velocity of the ion representing the solar wind plasma is generated. In general, taking the velocity distribution function of solar wind plasma at the location of the model particles, the velocity of the ion representing the solar wind is generated according to the frequency distribution function of a proton with velocity, v p , given by where v H is the initial hydrogen atom undergoing charge exchange, u p is the bulk plasma velocity, and v p th, = kT m 2 p H is the thermal speed of the plasma, as is defined in Malama (1991). A similar approach is used in Heerikhuisen et al. (2006).
Neutral atoms more frequently undergo charge exchange when the velocity difference between the proton and the neutral particle is maximized while charge exchange with solar wind particles whose velocity is near v H rarely occurs. That causes a depletion in the distribution around the neutral atom velocity, v H . An example 2D frequency distribution function, a 1D cut through the center of the distribution, and a comparison to the original Maxwellian distribution of the plasma are presented in Figure 1. In this example, the local bulk plasma velocity, u p , is less than the neutral hydrogen atom speed. Solar wind protons with velocities less than the neutral hydrogen atom's initial velocity are therefore more likely to undergo charge exchange than those moving faster than the neutral atom. The opposite is also true. If the local bulk plasma velocity is higher than the neutral hydrogen atom velocity, the probability of undergoing charge exchange is higher for protons with velocities higher than v H .
Rejection sampling, or the acceptance-rejection method, is applied to randomly select v p from Equation (14). The rejection method is a powerful, general technique that does not require taking the inverse of the cumulative distribution function (CDF), required for the inverse transform method. This avoids the issues with the discontinuity in Equation (14) when v p = v H and allows the selection to be done in 3D since the velocity components in Equation (14) are not independent. To sample a desired 1D function, f (x), the acceptance-rejection method generates a uniform distribution of points in both x and y. Points are accepted if they lie below the function such that y < f (x); otherwise they are rejected. To increase the acceptance rate of this method, another function, g(x), is introduced. The function g(x) is one that is easily sampled, such as a normal distribution, and, once multiplied by a constant, M, is always larger than the function, f (x), desired for sampling. The acceptance-rejection method uniformly samples from Mg(x) and accepts all the points that lie below f (x). This method is easily applied to 3D functions and has the benefit of only needing to evaluate functions at a given location instead of integrating to compute the CDF. Here we select g(x) to be a Maxwellian distribution and M is calculated to ensure that the proposed distribution is always larger than the frequency distribution function in Equation (14). The rejection method takes on the order of three iterations to return an acceptable value and can produce over 100,000 samplings in under a second.
Once the velocity of a solar wind proton participating in the charge exchange event is determined, a new neutral atom is generated and is assigned the selected proton velocity. The region of the heliosphere where the ENA is generated is also retained to track the propagation of each population of neutrals. The atoms can be sampled for each region to produce a solution for each population to determine the dominant population in each location as well as an easy comparison to the multi-fluid approach. One-dimensional and phase space distribution functions can also be produced for the total solution or for each population to determine how the populations evolve throughout the domain.
The DSMC model, however, allows for any form of the neutral distribution function. Following the work of Malama (1991), a statistical estimation of the source terms can be found by summing the changes in momentum and energy from individual charge exchange events in each cell over a given time interval. With this approach, the source terms at a given location, x, can be expressed as is the volume of cell k at location x, Δt is the time interval over which the charge exchange events occur, n ex is the number of charge exchange events that occur in the cell over that time span, v H,i is the velocity of the original hydrogen atom, and v p,i is the velocity of the initial proton selected according to the frequency distribution function in Equation (14). The parameter μ is the statistical weight, which is defined as the number of hydrogen atoms represented by each macroparticle that we use in the simulations. Following this definition, μ = (total flux of injected hydrogen atoms s 1 -[ ]) × (time step) / (the number of model particles that we want injected during a simulation time step). As all results of Monte Carlo modeling contain statistical noise, the latter was determined in our exploratory simulations such that the level of statistical noise in the numerical solution is small. In the example presented in Figure 1(a), the plasma bulk velocity is less than the neutral atom velocity. Protons with velocity near 0 km s −1 are much more likely to charge-exchange with this neutral. This results in an increase in energy for the plasma, which is characteristic in the disturbed ISM region of the OH. In the supersonic solar wind, the bulk plasma velocity is higher than the majority of parent neutral atom velocities. When this occurs, the most probable partner for the neutral is one having a higher relative speed, leading to the removal of momentum and energy from the plasma. The source terms need to be smooth from cell to cell in order for the resulting MHD solution to be stable. In order for this to occur, there need to be enough charge exchange events, n ex , within each cell to acquire accurate statistics. This can be achieved by increasing either the length of the time interval, Δt, before sampling or the total number of particles within the domain.

Solving the Self-consistent Problem
Due to the efficiency of analytic approximation for the source terms, a steady-state solution from the multi-fluid model is used as input to the SHIELD code. This brings the plasma closer to the final solution and reduces the amount of time needed to run the kinetic-MHD model. Figure 2 details the coupling between the Monte Carlo and MHD models over the course of one time step in the plasma fluid. The SHIELD code passes the MHD variables from BATS-R-US to AMPS. AMPS then pushes the neutral particles through the plasma solution. The neutrals are injected with a Maxwell-Boltzmann distribution from the west face of the domain. AMPS determines when and where charge exchange occurs for each particle, given the local plasma conditions. When an event occurs, an ENA is generated and the resulting change in momentum and energy of the neutral is added to the local plasma source terms within that cell. After the time interval, Δt, is reached, the source terms are then passed back to BATS-R-US and added to the momentum and energy equations. The plasma solution is then advanced a time step and the updated solution is then passed back to AMPS. This process continues until the solution relaxes into a new steady state.
Similar to those in Heerikhuisen et al. (2006), in the numerical code SHIELD, the source terms must be updated in every time step of the plasma as the solution progresses from the multi-fluid neutral solution, used to start the code, to a new steady state with kinetic treatment of the neutral hydrogen atoms. In this transition period, if the plasma solution is Figure 2. A schematic diagram summarizing the coupling between the OH and PT components within the SWMF. The parameters n, v p , B, and T p are the density, velocity, magnetic field, and temperature of the plasma. S ρ , S ρv , and S ò are the source terms to the continuity, momentum, and energy equations. advanced in time with constant source terms from a previous time step, then the SHIELD code develops a numerical instability in the flanks of the heliosphere. Waves form at the outer boundary, pulling the HS plasma out into the ISM. This instability does not occur in the Baranov & Malama (1993) model due to their treatment of the HP as a perfectly ideal, tangential discontinuity that does not allow the HS plasma to leak into the ISM.
The source terms are calculated over a time interval, Δt = Δn × dt, where dt is the time step of the neutral species determined from stability criteria and Δn is the number of steps or iterations over which the source terms are accumulated. For dynamic problems, the coupling frequency, and therefore Δt, is set by the minimum temporal scale resolved in the MHD solution. This ensures the source terms are updated frequently enough to resolve all desired phenomena. For short timescales, this requires a large number of simulated macroparticles to obtain the necessary statistics for smooth source terms. Steadystate problems are time-independent such that the terms in Equation (8) and the MHD equations do not vary with time. These runs can be handled differently from dynamic problems since we do not need to worry about resolving any temporal phenomena. The SHIELD code can therefore advance the Monte Carlo and MHD models at different rates, allowing AMPS to run for a longer period of time to develop smooth source terms before they are updated in the MHD equations and the plasma solution is advanced. A similar method was implemented by Heerikhuisen et al. (2006). This reduces the need for a very large number of macroparticles. The SWMF controls the frequency with which each component advances a time step. The SHIELD code cycles the components such that the MHD model can be restricted to advance a single time step over the Δt interval of the kinetic code. The interval Δt is set, along with the number of injected particles into the domain, by the user to ensure that the noise in the sampled moments of the neutral distribution function is on the order of 1%.
Additional techniques can also be applied to time-independent problems to increase efficiency and reduce numerical costs. Charge exchange in the heliosphere produces neutral populations with very different characteristic speeds. Setting a single time step for all particles can have challenges. In order to control the quality of the statistical sample, the time step is chosen such that particles cannot move through multiple cells within a single step. However, a time step that is too short could cause a large number of particles to build up when slower particles reside in the largest cells of nonuniform grids. This can cause the computational nodes to run out of memory. Therefore, we have allowed steady-state simulations to be able to separate the neutral hydrogen atoms into four energy bins. The four bins are separated in velocity space as follows: (1) |V| < 50 km s −1 , (2) 50 km s −1 < |V| < 150 km s −1 , (3) 150 km s −1 < |V| < 500 km s −1 , and (4) |V| > 500 km s −1 . Each neutral hydrogen species is allowed to have its own time step, dependent on the characteristic speed of that energy bin. This ensures that the aforementioned problems do not occur and allows the solution to relax to a steady state much sooner, significantly reducing the execution time of the model. The bins are selected to enhance the efficiency of the code and do not correspond to the populations used in the multi-fluid code. As mentioned in Section 2.3, the region of the heliosphere where each neutral atom is generated through charge exchange is also retained for direct comparisons between the models.

Model Comparison
There are no direct observations of neutral hydrogen atoms in the OH to validate the SHIELD code. Neither Voyager 1 nor Voyager 2 took direct measurements of neutral atoms. ENA maps also present a challenge for code validation. ENAs measured by IBEX and Cassini are produced primarily through neutral charge exchange with PUIs (Zank et al. 2010) and are line-of-sight measurements. While dynamic solar wind boundary conditions are included within our model (Michael et al. 2015), our MHD model does not include PUIs as a separate population. This makes direct comparison to ENA observations an ineffective way to verify the correct inclusion of the kinetic model.
In order to ensure the SHIELD code is coupled correctly, we compare our simulation to the results of the OH kinetic-MHD models of  and Heerikhuisen et al. (2006). Both of these models are 2D, axisymmetric, hydrodynamic simulations that treat the thermal solar wind ions, PUIs, and electrons as a single fluid. Each neglects photoionization and electron impact ionization; therefore, neutrals interact with the plasma only through charge exchange, as is done with the SHIELD code. In this work, we adopt the same boundary conditions for the solar wind and ISM as  and Heerikhuisen et al. (2006). We also neglect the magnetic field, as done in prior works. Our inner boundary is located at 30 au for the plasma unlike other models, which institute the boundary at 1 au. The solar wind conditions at 30 au are extracted from  using the WebPlotDigitizer software 4 and taken to be constant and uniform in both latitude and longitude with n SW = 0.008 cm −3 , V SW = 354.75 km s −1 , and T SW = 1.126 × 10 5 K. The ISM plasma enters the domain with the conditions n 0.06 p ISM = cm −3 , T ISM = 6527 K, and V ISM = 26.4 km s −1 parallel to the ecliptic plane. The interstellar neutrals are injected with a Maxwell-Boltzmann distribution with a density of n 0.18 H ISM = cm −3 and the same bulk velocity and temperature of the plasma.  and Heerikhuisen et al. (2006) both used the charge exchange cross section from Maher & Tinsley (1977) given by SHIELD is a 3D model. For this comparison the computational domain extends from −1000 au to 1000 au in all three directions. The ISM enters into the domain from the x = −1000 au face. In the nose of the heliosphere we use a grid resolution of 4 au inside a −400 to 400 au cube around the Sun to resolve the hydrogen wall and the disturbed ISM. The resolution increases to 12 au at the outer boundary. The nonuniform grid is static and not changed throughout the simulation. This grid is used for both the MHD and kinetic, and multi-fluid models to ensure consistency between the multi-fluid and kinetic approach. The mean free path of the neutral atoms is on the order of 100 au (Izmodenov & Alexashov 2015); therefore, a 4 au resolution is ample to resolve the neutral hydrogen solution within AMPS. Since our model is 3D instead of 2D, we will qualitatively compare the solution in the meridional cut of our 3D simulation to the works of  and Heerikhuisen et al. (2006). However, a 2D asymmetric model with the correct boundary conditions should yield the same solution as its 3D counterpart. This comparison is done for both the multi-fluid model, using the four neutral fluid approximation, and the SHIELD code.
The multi-fluid approximation is more computationally efficient than the K-MHD model; therefore, it is used to relax the plasma to a steady solution.  showed that the four neutral fluid approximation most closely matches the solution run with their K-MHD model. In an effort to obtain the plasma solution closest to the true solution to reduce the execution time of the SHIELD code, we also describe the neutrals with four separate neutral fluids. The neutral populations are separated as follows. Pristine ISM neutrals or Population IV neutrals enter the domain with values set by the boundary conditions. For these conditions, a BS forms within the solution. Population IV neutrals are defined as being supersonic with flow speed less than 140 km s −1 , distinguishing them from particles born in the supersonic solar wind. The region between the BS and the HP is the disturbed ISM, or Population I neutrals, defined to have a sonic Mach number less than 1 and a temperature below 7.5 × 10 4 K. Population II neutrals in the HS occur where the plasma is subsonic at temperatures higher than that of Population I neutral hydrogen. Finally, Population III neutrals are generated in the supersonic solar wind. The numbers of neutral populations are chosen to match previously published studies using this multi-fluid code. Unfortunately, this differs from the numeration used by ; therefore, we note that their results have been renumbered to match the definitions used here. A steady solution for the multi-fluid model is reached after the model is advanced for over 50,000 time steps.
The multi-fluid solution is used to start the fully coupled, self-consistent SHIELD code. The plasma solution from the multi-fluid model is passed to the SHIELD code. Neutrals are then propagated through the domain to generate a selfconsistent kinetic solution that can be used to update the plasma with the corresponding source terms. The AMPS domain is initially empty; therefore, neutrals are pushed through the fixed plasma solution of the multi-fluid model.
For a fixed plasma, the MHD model and Monte Carlo code are one-way coupled. The plasma solution is passed to the kinetic code allowing the kinetic solver to generate a neutral solution but not allowing the source terms from charge exchange to be added to the MHD equations and alter the plasma. AMPS is therefore able to simulate the transport of neutral particles through the plasma, allowing for the generation of an accurate neutral distribution, while holding the plasma solution constant. AMPS is run with the fixed plasma solution for 71,000 time steps to generate the neutral distribution. With the given time step, neutral particles can traverse the domain in 1200 time steps, allowing for many cycles of streaming particles. This generated neutral solution is used along with the plasma solution from the multi-fluid simulation to start the SHIELD model.
The problem is then treated as time-independent and the neutral hydrogen atoms are modeled as four separate ENA species for increased efficiency. One hundred thousand particles are injected at the x = −1000 au face in every time step, resulting in around 100 million particles residing within the domain at any given time. The source terms accrue for Δt = 10,000 time steps. This combination of injected particles and Δt is sufficient to reduce the statistical error in the Monte Carlo simulation to about 1% in the neutral density, velocity, and temperature. This ensures accurate source terms are obtained before the resulting source terms are passed back to the MHD solver and the plasma solution is updated. The OH component is cycled such that the MHD solution is advanced one time step for every Δt interval of time that AMPS is run to calculate the source terms within the PT component. This allows statistics to accumulate over a long period of time between each individual step in the plasma solution to reduce the number of model particles needed. When the simulation is run on 1000 processors, one Δt interval takes roughly 20 minutes of wall-clock time. The cycling procedure continues until the plasma relaxes to a new solution. This occurs after the MHD solution is advanced 890 time steps.

Results
As stated in the introduction, a consensus has not been reached among the most recent 3D kinetic-MHD models (Pogorelov et al. 2013;Izmodenov & Alexashov 2015) on the impact the different treatments of the solar magnetic field have on the global structure of the heliosphere. The varying plasma solutions between the models offer little in verifying that we have coupled the kinetic and MHD models correctly. We, therefore, rely on comparison to the simplified, hydrodynamic simulations of  for model verification, as done in Heerikhuisen et al. (2006). Future studies will explore a more detailed comparison of the 3D kinetic-MHD models.

Fixed Plasma Solution
It is important to ensure that the kinetic solution predicts the correct neutral solution before analyzing the fully coupled self-consistent SHIELD code. We, therefore, first look at how the neutral solution differs from the work of  using a fixed plasma solution using the method described at the end of Section 2.5. The plasma solution used for this comparison is taken from the steady state of the fully coupled SHIELD code. The neutral populations within the SHIELD code are selected with the same criteria as those for the multi-fluid model described previously. Figure 3 compares the neutral density, speed, and temperature for the kinetic solution of the SHIELD code to those for the kinetic solution of  for each population of neutrals. The numbers refer to the respective region of the heliosphere where they are created. In our scheme Populations IV, I, II, and III correspond to neutrals born in the pristine ISM, between the BS and HP, in the HS, and in the supersonic solar wind, respectively.
A comparison of our kinetic solution to that of  allows us to verify that our kinetic solver is determining the correct distribution for the interstellar neutrals. The data from  is extracted using the WebPlotDigitizer software (see footnote 4).  used a different numbering scheme for their neutral populations, referring to the pristine ISM, the region between the BS and HP, the HS, and the supersonic solar wind as Populations IV, III, II, and I, respectively. Their populations have been renamed to match the definitions used in this work.  used the fixed plasma solution from their 2D self-consistent K-MHD model. The kinetic solutions show very good qualitative and quantitative agreement.
Both kinetic models predict higher filtration of interstellar neutral populations, IV and I, into the heliosphere than their multi-fluid counterparts. The higher filtration of neutral particles into the heliosphere occurs within the kinetic model due to the large mean free path of hydrogen atoms. There are no collisions between hydrogen atoms; therefore, the large charge exchange mean free path (>100-200 au) allows kinetic neutral hydrogen atoms to travel into the upwind direction from the sides or flanks of the heliosphere. The collisions between neutrals dominate in the fluid description, causing the neutrals in the multi-fluid model to not propagate as far for particles that would be crossing flow streamline. This causes the multi-fluid model to more closely follow the streamlines. Less contribution Populations I, III, and IV match very well between the kinetic models throughout the entirety of the domain in both value and behavior. AMPS predicts a denser, cooler neutral population originating in the HS than . This population, however, is very sensitive to how the HP is defined when separating the populations. The difference in our solutions is most likely due to a different definition separating hydrogen wall neutrals from neutrals born in the HS. If some of the neutrals originating in the hydrogen wall are included in the HS population, the population will appear to be denser and cooler as in the SHIELD code. Overall the kinetic solution within SHIELD matches the results of  very well and we can be confident that AMPS provides an accurate representation of the interstellar neutrals as they propagate through the heliosphere.
Additionally, the kinetic treatment shows that the average speed of neutrals originating in the pristine ISM (Population IV) increases as they propagate into the heliosphere. Faster interstellar neutrals have larger charge exchange mean free paths as they travel through the hydrogen wall. As a result, slower particles undergo charge exchange more frequently. Charge exchange removes the slower neutrals from Population IV causing an increase in the average speed of the population. those of this work (f). The density, speed, and temperature are calculated by integrating the moments of the neutral atom velocity distribution function within the kinetic code. In order to calculate the total neutral solution from our multi-fluid approximation, the speed and temperature are taken to be weighted averages of the populations at each location while the density is calculated as the sum of the populations. As expected, the kinetic model produces a much different neutral solution in comparison to the multi-fluid model. The hydrogen wall in the multi-fluid model is denser, reaching a peak of n H = 0.3 cm −3 . The hydrogen wall within SHIELD has a peak value of n H = 0.26 cm −3 located closer to the HP. The neutrals within the hydrogen wall are faster in the kinetic model as well. The main feature seen in both K-MHD models is that the filtration of neutrals into the heliosphere is 20% higher. Both kinetic solutions predict a denser, slower population of neutrals propagating through the heliosphere in comparison to their multi-fluid counterparts.

Fully Coupled Self-consistent Kinetic Solution
Furthermore, SHIELD predicts the same penetration of neutrals into the heliosphere, which will act to increase charge exchange within the supersonic solar wind and HS. The main feature seen in both K-MHD models is that the filtration of neutrals into the heliosphere is 20% higher than that in their respective multi-fluid neutral models. A comparison of the plasma solution in the meridional plane between the kinetic neutral treatment in SHIELD and the multi-fluid model is presented in Figure 5. The kinetic treatment of the neutrals smooths out a few features in the tail that arise in the multi-fluid model, but overall, the solutions are qualitatively very similar. The bottom panels of Figure 5 provide 1D profiles along the nose of the heliosphere, similar to Figures 3 and 4. As seen in Figure 7. The plasma density, speed, and temperature in the upwind direction for the fully coupled, self-consistent kinetic model. The results of SHIELD are in black while the results of  are shown in blue. Figure 6. The neutral density, speed, and temperature in the upwind direction for the fully coupled, self-consistent kinetic model. The results of SHIELD are presented in black while the results of  are shown in blue.
the results of , the kinetic treatment of the neutrals causes the TS and HP to move toward the Sun by 1 au and 6 au, respectively, as a result of the additional charge exchange from the higher density of interstellar neutrals within the heliosphere in comparison to the multi-fluid model.  found that the kinetic treatment reduces the distance to all heliospheric boundaries; however, the BS moves outward by 6 au in the SHIELD model resulting in a wider plasma pileup region between the HP and the BS in comparison to the multi-fluid model. The outward motion of the BS in the kinetic model was also seen by Heerikhuisen et al. (2006). Figures 6 and 7 compare 1D profiles in the nose of the heliosphere from the fully coupled, self-consistent results from SHIELD for both the neutrals and the plasma, respectively. The kinetic solution of  is also shown. In comparing SHIELD to the well-established model results of , we see overall excellent agreement between the two approaches. Both models qualitatively predict similar structures within the neutral solution. The kinetic models have very good quantitative agreement at the peak of the hydrogen wall in the density, speed, and temperature of the neutrals as well as at the values within the heliosphere. The kinetic treatment of the neutrals is a significant improvement of our model to match the neutral solution of . Heerikhuisen et al. (2006) performed a similar run when developing the K-MHD model within the Multi-scale Fluidkinetic Simulation Suite. The TS, HP, and BS locations within all three kinetic simulations can be seen in Figure 8(b) and are presented in Table 1 along with the multi-fluid simulations performed in each study with their respective models. Here, the HP location in the SHIELD code is determined by the velocity streamlines. Both  and Heerikhuisen et al. (2006) observed kinetic neutrals cause the TS to move 5-6 au toward the Sun due to the increased hydrogen filtration into the heliosphere. The SHIELD code predicts this as well, although to a lesser extent, and has a final TS 2-3 au farther than the other kinetic models. Overall the SHIELD code matches the heliospheric locations of  very well. Figure 8 directly compares the model results from the SHIELD code to those of  and Heerikhuisen et al. (2006) along the nose of the heliosphere. The SHIELD code shows excellent agreement with both models. The SHIELD code closely reproduces the neutral density profile in . The BS in the SHIELD code is 5 au farther than that in , leading to a slightly larger hydrogen wall thickness. All three models predict similar neutral density filtration into the heliosphere. A comparison of the plasma temperature in the kinetic models is presented in Figure 8(b). The TS, HP, and BS locations can be seen within the models. The SHIELD code predicts a slightly cooler supersonic solar wind but shows very good agreement in the HS and in the shocked ISM, predicting a BS between those of the other two models. One of the main differences between the SHIELD code and the model of  is the Figure 8. Neutral density (a) and plasma temperature (b) profiles along the nose of the heliosphere. The plots compare the results of the kinetic models of blue) and Heerikhuisen et al. (2006;red) to those of the SHIELD code (black). Note. A&I denotes the work of  and HFZ that of Heerikhuisen et al. (2006). Values are shown for the kinetic and multi-fluid treatment of neutrals, represented by K and MF, respectively. discontinuity algorithm used in their model. This can be seen by the vertical lines in their solution at the HP. The solution within SHIELD shows a much more gradual decrease from the temperature within the HS to the disturbed ISM ahead of the HP. This gradient is also seen in the model of Heerikhuisen et al. (2006); however, it appears more evident in the SHIELD code due to a slightly coarser resolution in the area. The largest discrepancy within the models occurs between the BS and the HP. The SHIELD code predicts similar neutral hydrogen density and plasma temperature to those of . The model of Heerikhuisen et al. (2006) predicts a hotter region that pushes the BS farther from the Sun. This results in a wider hydrogen wall whose density peaks 30 au farther away from the HP than that of the SHIELD code. Overall, the SHIELD code has excellent agreement with the results of , validating that the PT component has been applied correctly to the OH and produces a correct solution that we can trust moving forward.

Summary
In this work, we present the SHIELD code, a 3D, global, self-consistent K-MHD model of the OH implemented within the SWMF. The SHIELD code couples AMPS, a Monte Carlo code, to our MHD solver to treat the neutrals streaming through the heliosphere kinetically. SHIELD codes the thermal protons, PUIs, and electrons as a single plasma fluid and includes hydrogen as the only neutral species. Here, H-H collisions, photoionization, and radiation pressure are neglected. Therefore, the only way the plasma and neutrals interact is through charge exchange.
We compare the results of SHIELD to the results of  to verify that the kinetic treatment of the neutrals is incorporated correctly. Overall, we find excellent agreement in the neutral distributions between the models. SHIELD can accurately reproduce the behavior of each population of neutrals as they propagate through the system and matches the higher filtration of neutrals into the heliosphere. The SHIELD code accurately reproduces the locations of the heliospheric boundaries predicted by the K-MHD results of , including the inward motion of the TS as well as the HP distance in the nose. The SHIELD code is a significant improvement over our multi-fluid model in matching the results of other models. The excellent agreement between the models makes us confident that SHIELD is coupled correctly and can be used for future work investigating the heliosphere.
The plasma is also strongly affected by the different modeling techniques of the neutrals. The plasma solution in the SHIELD code shows excellent agreement with the results of . Like other works, we conclude that the differences between the kinetic distributions and the multi-fluid approach are large and that the kinetic solution should be used when comparing model results to observations.