A Comparison of Particle-in-cell and Hybrid Simulations of the Heliospheric Termination Shock

We compare hybrid (kinetic proton, fluid electron) and particle-in-cell (kinetic proton, kinetic electron) simulations of the solar wind termination shock with parameters similar to those observed by Voyager 2 during its crossing. The steady-state results show excellent agreement between the downstream variations in the density, plasma velocity, and magnetic field. The quasi-perpendicular shock accelerates interstellar pickup ions to a maximum energy limited by the size of the computational domain, with somewhat higher fluxes and maximal energies observed in the particle-in-cell simulation, likely due to differences in the cross-shock electric field arising from electron kinetic-scale effects. The higher fluxes may help address recent discrepancies noted between observations and large-scale hybrid simulations.


Introduction
The heliospheric termination shock, located ≈85-95 au from the Sun (as measured by Voyager 2 in 2007 and Voyager 1 in 2004, respectively), abruptly changes the properties of the solar wind (Decker et al. 2005;Richardson et al. 2008).Both spacecraft made multiple crossings of the shock in a short period (Burlaga et al. 2008), likely due to its nonstationarity and, in particular, self-reformation driven by the accumulation of reflected ions in the foot region ahead of the shock ramp (Lembège & Savoini 1992).
To a reasonable approximation, the termination shock is a fast-mode MHD shock, across which the magnetic field embedded in the solar wind increases in strength, while the wind itself is compressed, slowed, and heated.Interestingly, the thermal component downstream of the shock contains only ≈20% of the upstream ram energy, with much of the remaining energy likely deposited in nonthermal (energies 1 keV) pickup protons-the nuclei of hydrogen atoms of interstellar origin that are ionized within the heliosphere and then coupled to the magnetic field embedded in the solar wind-that cannot be directly detected by the Voyager instruments (Zank et al. 1996).Multifluid simulations of the Voyager 2 crossing support this view by demonstrating that the shock can transfer a significant fraction of the upstream ram energy to the pickup ions (PUIs; see Figures 5(b) and 7(b) in Zieger et al. 2015).
Multiple works have considered the processes controlling the acceleration of charged particles at the termination shock, which is, on large scales, quasi-perpendicular (i.e., the upstream magnetic field makes an angle >45°with the shock normal).Of particular concern is the injection problem, the difficulty of accelerating low-energy particles in the context of the theory of diffusive shock acceleration.1D hybrid (particle ion, fluid electron) simulations showed that injection does not pose a critical barrier for quasi-parallel shocks in the presence of PUIs (Kucharek & Scholer 1995;Lipatov & Zank 1999).While the termination shock is quasi-perpendicular, the same papers argued that the presence of turbulence can trigger the development of local quasi-parallel regions that could account for the observed fluxes at high energies.A Monte Carlo simulation of PUI acceleration with an empirically based scattering rate also suggested that the injection problem can be circumvented (Ellison et al. 1999).Later, Giacalone (2005) included the effects of large-scale turbulence explicitly, demonstrating that even thermal solar wind ions can be accelerated to nonthermal energies.In further hybrid simulations, Giacalone & Decker (2010) concluded that observations from Voyager 2 reported by Decker et al. (2008) are consistent with the acceleration of PUIs by the termination shock via a combination of shock drift acceleration and particle scattering in meandering magnetic fields.
Full-particle simulations (particle ion, particle electron), which incorporate a more complete treatment of the underlying physics at the expense of additional computational cost, have also considered the problem.(Although only hybrid and fullparticle methods are discussed in detail here, other approaches have also been investigated, e.g., the direct solution of the Vlasov equation by Ariad & Gedalin 2013.)Simulations of a strictly perpendicular 1D shock that included a population of PUIs showed the development of an extended foot and crossshock electric field that resembled those found in previous hybrid simulations (Scholer et al. 2003;Chapman et al. 2005), but did not find efficient PUI acceleration via shock surfing.Later work demonstrated that the concentration of PUIs has a significant effect on the transfer of bulk dynamic energy to the thermal energy of solar wind protons (Matsukiyo & Scholer 2014;Lembège & Yang 2016).Companion 2D full-particle simulations recovered similar results regarding the structure and dynamics of the shock front and the dynamics of the partition of energy (Yang et al. 2015).Larger particle-in-cell (PIC) simulations concluded that PUIs upstream of the shock are energized by adiabatic compression of the solar wind as well, due to an enhanced level of turbulence in a broad foreshock region (Kumar et al. 2018).
If the hybrid and full-particle approaches were known to yield similar results, the former would generally be preferred, since larger computational domains can be studied for the same expense.However, assessing the agreement between the two approaches in previously published works is complicated by the lack of consistency in the parameters controlling the various runs.In this study, we examine the similarities and differences between hybrid and PIC simulations of the solar wind, magnetic field, and particle distributions near the termination shock at the time of the Voyager 2 crossing.

Description of Simulations
We compare results from two computational models.One, the hybrid model, includes kinetic protons, but treats the electrons as a massless, charge-neutralizing fluid.The second, the PIC model, treats both protons and electrons kinetically.The inclusion of kinetic electron physics leads to additional computational burdens.Both codes must resolve the key ion lengthscales-e.g., the proton inertial length d p = c m n e 4 p p 2 p , where n p and m p are the proton number density and mass, respectively-and timescales such as the proton cyclotron time p 1 W -, where Ω p = eB/m p c is the cyclotron frequency based on the local magnetic field.However, PIC simulations must resolve smaller lengthscales-the Debye length usually being the smallest-and timescales, such as the electron cyclotron period.As a consequence, while PIC simulations employ a more complete physical representation, they also incur a greater computational expense.The degree to which the two models agree is closely related to the importance of electron kinetic effects in the system being modeled.
The hybrid code used here is described in Giacalone & Decker (2010) and has seen wide use over nearly 30 yr in studies of space plasmas-shocks in particular (see Winske & Omidi 1996 and references therein).For the fully kinetic simulations, we use the PIC code p3d (Zeiler et al. 2002).In both codes, a reference magnetic field strength B 0 and density n 0 define the Alfvén speed v B m n 4 .To simplify notation, the "0" will henceforth be dropped from most subscripts.
The initial conditions resemble those in Giacalone & Decker (2010), although the computational domain in this work, which has dimensions L x × L y = 409.6dp × 25.6d p , ≈0.020 × 0.0012 au, is significantly smaller, because their 500 × 4000d p domain is impractical for a PIC simulation.While the system is invariant in the third dimension (∂/∂z = 0), particles can move in that direction, making these 2 × 3v or 2.5D simulations.In the short dimension (L y ), the boundaries are periodic, while the two ends in the long direction are capped by hard conductors that specularly reflect particles.At t = 0, all of the plasma is moving in the positive x-direction.Reflection from the x = L x boundary drives counterstreaming flows that form a shock.The disturbance quickly reaches a time-averaged steady state that moves upstream for the remainder of the simulation, which ends once perturbations from the x = 0 boundary reach the shock.
The initial magnetic field is uniform, of unit magnitude, and is inclined at an angle of θ Bn = 70°with respect to the long dimension (which parallels the shock normal), while the initial electric field is taken to be the usual motional one given by −v × B/c, so that, given the geometry, the only initial component is in the z-direction.The initial B lies completely in the x − y plane.Inclining the field into the out-of-plane direction could, because of the 2.5D geometry, affect the development of fluctuations in the system and, ultimately, the shock itself.However, Kumar et al. (2018) investigated the effects of magnetic field orientation and concluded, after comparison with fully 3D simulations, that the choice made here offers the best representation of the full system.The chosen value of θ Bn is somewhat more oblique, i.e., farther from θ Bn = 90°, than that used in previous simulations based on Voyager data (Giacalone & Decker 2010;Giacalone et al. 2021)-although it does agree with the value determined in Li et al. (2008).Regardless, as discussed in Giacalone et al. (2021), the blunt shape of the shock (i.e., its departure from spherical symmetry) and the ripples that develop at the shock front mean that deviations of θ Bn from a more perpendicular incidence can occur.Since we are chiefly concerned not with a direct comparison to Voyager data, but instead to the similarities and difference between the models, the crucial point is that the same θ Bn is used in both.
The initial proton distribution combines thermal and pickup components in a density ratio of 3:1, both uniformly distributed in space.The thermal component has a Maxwellian velocity distribution with a temperature of T m c 0.25 p A

=
and an associated β SW = 0.5 that is somewhat larger than the value of 0.139 used in Giacalone & Decker (2010).Meanwhile, the pickup component is distributed uniformly over a thin shell in velocity space with a radius given by the initial upstream solar wind speed in the shock frame.The addition of a small thermal spread to the PUIs, drawn from a Maxwellian with T = 0.025, somewhat smears their initial distribution.This choice resembles the distribution of a freshly ionized PUI population; more sophisticated PUI distributions have been suggested (e.g., Chen et al. 2013), but the rapid evolution at the shock seen in this work and others (Lembège & Yang 2016) suggests the details have minimal effects on the results.While Giacalone & Decker (2010) include a third, even higher-energy population, we do not, because the Larmor radii of such ions would equal or exceed L y .The electrons have a temperature of 0.25 and a density equal to the total ion density.In order to reduce the separation of scales, and hence the computational expense, the electron mass is taken to be m e /m i = 0.01, which implies the electron inertial length d e = 0.1d p .(Matsukiyo et al. 2007 investigated the impact of the mass ratio on quasi-perpendicular shocks in the relevant parameter regime, concluding that it has little effect on the overall structure, but does affect the details of the instabilities at the shock front.)For similar computational reasons, the ratio of the speed of light to the Alfvén speed (equivalently, the ratio of the proton plasma and cyclotron frequencies) is taken to be 20.Giacalone & Decker (2010;and, in a similar run, Giacalone et al. 2021) also included a spectrum of magnetic and velocity fluctuations at t = 0 as a representation of solar wind turbulence, with lengthscales approaching the system size of ≈10 4 d p .The smaller size of our domain does not permit the inclusion of such turbulence.Instead, the initial state is laminar, although small-scale fluctuations in all quantities develop almost immediately once the simulations begin.

Results
Figure 1 compares the PIC and hybrid models at the ends of the simulations (t 40 p 1 = W -). Each panel shows a quantity averaged over the y-direction (i.e., the short dimension) for a portion of the x-direction (the long dimension) centered near the shock front.The PIC results have been shifted in x to align the curves.The models agree closely for the profiles of the total (thermal + pickup) proton density (panel (d)), the proton velocity along the shock normal (panel (c)), and |B| (panel (b)).The secondary peaks downstream of the shock front, most visible in |B| and n, have been noted in previous simulations (Lembège & Yang 2016;Kumar et al. 2018) and arise from the reflection and transmission of gyrating thermal ions.Perhaps the most notable differences occur in the profile and magnitude of E x , the cross-shock electric field, which is ≈two to three times larger and ≈two times sharper in the PIC simulation (note that this implies the cross-shock potentials are roughly equivalent).These differences are due to electron kinetic effects that are not part of the hybrid formulation.
Figure 2 shows the evolution of the proton density in the PIC simulation frame of reference.Moving upward, consecutive curves are separated in time by 2 p 1 W -, with the top curve the same as that in Figure 1(d).The basic features of the shock appear early in the simulation and exhibit only gradual evolution.(This contrasts with PIC simulations of shocks without PUIs, which typically show shock reformation and more dynamic stack plots.)A secondary peak exists ≈4d p downstream of the density maximum, with other lesser peaks appearing at offsets of ≈17d p and 30d p .These peaks, as will be seen in Figure 5, are primarily due to the thermal solar wind ions, although their characteristics-the heights, widths, and distances from the shock front-depend on the details of the upstream PUI distribution (Lembège & Yang 2016;Kumar et al. 2018).The velocity of the shock front, ≈ −2.6vA in the simulation frame, remains nearly constant.
2D images of the proton density from the PIC simulation, overplotted with magnetic field lines, are shown in Figure 3  Ripples in the shock front with a wavelength of ≈4d p develop early in the simulation (not shown) and persist throughout its further evolution.Similar structures have been noted in hybrid (Winske & Quest 1988;Lowe & Burgess 2003;Burgess 2006) and PIC (Yang et al. 2015) simulations of strong quasiperpendicular shocks and are also seen in the hybrid simulation discussed in this work, albeit at somewhat longer wavelength and lower intensity (e.g., peak densities of 5-6, rather than the 6-7 seen in Figure 3).Other PIC simulations (Lembège & Yang 2018) found that varying the concentration of PUIs can stabilize the rippling and shock reformation process.Ion reflection at the shock (see Figure 5) generates anisotropy in the proton distributions.The anisotropy excites an electromagnetic instability that grows to sizeable amplitude while rippling the shock surface, after which the perturbations convect downstream where the free energy goes into heat and electromagnetic waves.Analysis of magnetic turbulence by Winske & Quest (1988) shows that the Alfvén ion cyclotron instability, rather than the drift mirror instability, is most likely the responsible mode, although the large gradients at the shock complicate the identification.These ripples are sufficiently intense to modify the local θ Bn , in some locations making the shock quasi-parallel (θ Bn < 45°) and in others nearly fully perpendicular (θ Bn = 90°), with the former case facilitating the reflection of particles back upstream.While the box size is limited in the y-direction (L y = 25.6dp ≈ 10 −3 au for n = 10 −3 cm −3 ), similar variations in the shock geometry Interactions with the shock energize particles.The crossshock potential produces an electric field opposing the bulk motion of the incoming plasma.This field, with a scale on the order of the PUI gyroradius, redistributes the energy of the incoming flow., where U 1 , the plasma ram speed, is ≈6V A .In both panels, the black curves come from the hybrid simulation and the red curves from the PIC.The upstream distributions, shown in panel (a), include two components: protons from the thermal solar wind (dotted lines) and newly ionized PUIs (solid lines).The thermal protons have a Maxwellian distribution in the frame moving with the upstream bulk flow, while the shell distribution of the pickups has a sharp cutoff at twice the solar wind speed.Panel (b) shows the differential energy downstream of the shock at the end of the simulation, t 40 p 1 = W -. As expected, the PIC and hybrid simulations closely match in the top panel.The minor discrepancies (e.g., the features at low energies in the hybrid thermal population and the small difference in the pickup cutoff energy) have an insignificant effect on the subsequent evolution.The middle panel shows the  distributions downstream of the shock.While the low-energy components of both the thermal and pickup distributions are in good agreement, the PIC simulation achieves higher energies for both populations.The higher energies (a factor of ≈1.5) are presumably due to the shorter lengthscale associated with the larger cross-shock electric field (see Figure 1(a) and the discussion in Zank et al. 1996).The maximum energy that can be achieved is limited by the finite size of the computational domain (in particular L y = 25.6dp ) and the duration of the run.The Larmor radius of the highest-energy protons, those with energy E m U 20 2 , is ≈25d p .Larger simulations, run for longer times, would continue to fill the high-energy tail of the distribution, an effect clearly seen in PIC simulations with box sizes L x /2 × L y /2 and L x /4 × L y /2 (not shown).In significantly larger hybrid simulations (Giacalone & Decker 2010;Giacalone et al. 2021), the high-energy tail eventually merges into the low-energy end of the anomalous cosmic ray distribution, which has a different slope and cannot be captured here.The development of the distributions as a function of time in the PIC simulation is shown in panel (c) and demonstrates that it has reached a steady state.(Although not shown here, the hybrid simulation has as well.)The qualitative shape of the both the thermal and pickup components develops early (compare with panel (a)), with subsequent evolution primarily pushing the maximum energy to higher values.
The difference in the simulation fluxes at high energies may be relevant for a discrepancy noted in Gkioulidou et al. (2022) and Kornbleuth et al. (2023): termination shock ion fluxes between 3 and 20 keV inferred from IBEX and INCA energetic neutral atom observations are up to a factor of 5 larger than those seen in the hybrid simulation of Giacalone et al. (2021).This energy range, a few times the initial energy of the PUI distributions in Figure 4, is also where the fluxes from the PIC simulation exceed those from the hybrid.If the differences between the two models seen in Figure 4 persist in simulations with parameters similar to those examined in (Giacalone et al. 2021; e.g., much larger computational domains and, perhaps, the injection of upstream turbulence), the discrepancy with observations could be explained.
The interaction of the protons with the shock can be seen in Figure 5, which shows the x − v x phase space for the initial thermal (top) and pickup (bottom) populations at the end of the simulation.Velocities are shown in the simulation frame, with the shock velocity (negative since it is moving leftward) indicated by the horizontal dashed line.The vertical dashed line denotes the shock front, defined as the density peak in Figure 1.Reflected thermal ions extend ≈2d p upstream of the shock front and are the population responsible for exciting the Alfvén cyclotron instability that produces the ripples seen in Figure 3.The density peaks discussed above, at x ≈ 228, 232, 245, and 258, are most apparent in the thermal population, whose number density shows anticorrelations with those of the much less dense pickups.Heating at the shock, particularly of the pickup protons, is evident in the increase in vertical extent of the phase-space distributions to the right of the dashed line.The cross-shock potential enables the pickup particles, which are basically unmagnetized across the shock front, to be injected into the downstream with energy in excess of their upstream bulk motion, thus facilitating greater energy gain than the thermal protons.

Summary
We used a PIC simulation to explore the acceleration of protons at the heliospheric termination shock and compared the results to a hybrid model (kinetic protons, fluid electrons).The steady-state profiles show excellent agreement, particularly the variations in the density, plasma velocity, and magnetic field.The proton distributions downstream of the shock broadly agree, although the PIUs in the PIC simulation reach somewhat higher (≈50%-100%) energies and fluxes, likely due to the sharper lengthscales associated with the cross-shock potential and electric field.This suggests that larger hybrid runs-such as those described in Giacalone & Decker (2010) and Giacalone et al. (2021)-capture most of the physics necessary to simulate the termination shock.However, these hybrid simulations predict ion fluxes at moderate energies (≈3-20 keV; roughly equivalent to the highest energies seen in the simulations in this work) that are lower than those inferred from observations (Gkioulidou et al. 2022;Kornbleuth et al. 2023).Proving that the differences between the PIC and hybrid models explored here can explain the discrepancy between hybrid simulations and observations requires additional, much more expensive, PIC simulations that approach the scale studied in Giacalone et al. (2021).
In this work, the maximum energy of the PUIs changed by a factor of ≈4-5 after passage through the shock (see Figure 4), while in the hybrid simulation described in Giacalone & Decker (2010)-in which the simulation domain was ≈1.2 times longer in the long direction and 150 times longer in the short direction; approximately 0.021 × 0.17 au-the energy increase was a factor of ≈30.Reaching the scales of the actual termination shock would require substantially larger domains, implying either the need for extrapolation from current simulations or the use of novel computational methods.For instance, the kglobal code, described in Arnold et al. (2021), consists of an MHD fluid backbone with guiding-center particles and produces nonthermal power-law energy distributions without having to resolve kinetic scales.Its applicability to studies of the termination shock remains an open question.
normalized to the proton inertial length d p0 = c/ω p0 (where ω p0 is the proton plasma frequency) and times to the ion cyclotron time p0 1 W -. Electric fields and temperatures are normalized to v A0 B 0 /c and m v p A0 2 , respectively.For typical values near the termination shock, n = 10 −3 cm −3 and B = 0.05 nT, d p0 ≈ 7 × 10 3 km,

Figure 1 .
Figure 1.Comparison of hybrid (black) and PIC (red) simulations in a portion of the domain centered near the shock front and averaged over the y-direction.The PIC results have been shifted in x to align the curves.Panels (a)-(d) show E x , |B|, the proton velocity along the shock normal direction, and the total (thermal + pickup) proton density, respectively.
Figure 4 shows dJ/dE, the omnidirectional differential intensity (sometimes called "flux") as a function of the kinetic energy E of the protons in the two simulations.The measurements are made in the shock rest frame, with an energy unit of m U 2 i 1 2

Figure 2 .
Figure 2. The proton density, averaged over y, in the simulation frame.Each cut is separated from the previous by 2 p 1 W -and offset upwards.The x-axis has been shifted to match Figure 1.

Figure 3 .
Figure 3. 2D images of the total proton density, normalized to the upstream value, at t 30 p 1 = W -(top) and t 40 p 1 = W -, and overplotted with magnetic field lines (white).The x-axis has shifted to match Figure 1.

Figure 4 .
Figure 4. Omnidirectional differential intensity (flux) vs. kinetic energy of protons in the shock rest frame for the hybrid (black) and PIC (red) simulations.The energy unit, m U 2 i 1 2 , corresponds to the ram speed of the incoming flow.Panel (a) Upstream distributions showing thermal solar wind protons (dashed lines) and freshly ionized pickup protons (solid lines).Panel (b) Downstream distributions at t 40 p 1 = W -. Panel (c) Development in time of the distributions in the PIC simulation.The

Figure 5 .
Figure 5. Position-velocity (x − v x ) phase-space density, on a logarithmic scale, in the PIC simulation for protons at t 40 p 1 = W -taken over a narrow range in y (≈1d p ). Panel (a) The distribution of protons that were part of the thermal component at t = 0. Panel (b) The distribution of protons that were in the pickup component at t = 0. Velocities are in the simulation frame with the dashed horizontal line showing the shock velocity.The vertical dashed line shows the position of the density maximum seen in Figure 1.