A Global MHD Simulation of Outer Heliosphere Including Anomalous Cosmic-Rays

A global MHD–neutrals–cosmic-rays simulation is conducted to investigate the effects of anomalous cosmic-rays (ACRs) on the large-scale structure of the outer heliosphere. In the model, the cosmic-rays are treated as a massless fluid that only contribute their pressure to the dynamics of the system. The diffusion of cosmic-rays in the interstellar medium is assumed to be much faster than inside the heliosphere, where it depends on the strength of the magnetic field. The results show that the influence of the cosmic-rays on the structure of the outer heliosphere depends on momentum and energy transfer from the solar wind plasma to the ACRs, accomplished through diffusive shock acceleration at the termination shock, and the subsequent loss of ACRs across the heliopause and their rapid escape into the interstellar medium. Under favorable conditions characterized by a large fraction of energy conversion and a high enough diffusion coefficient in the solar wind, the ACRs were found to reduce the width of the heliosheath by up to ∼18 au. Consequently, these results indicate that the effect of cosmic-rays is a potential key factor for the global structure of the outer heliosphere in a computer model that could partially explain the timing of the heliopause encounters of the two Voyager probes.


Introduction
The heliosphere is formed by the interaction between the expanding supersonic solar wind and the Sun's motion through the local interstellar medium (LISM), conveniently viewed as the relative motion of the LISM past the stationary solar system. Both the solar wind and the LISM are multicomponent in nature, consisting of relatively low temperature plasma, hot pickup ions (PUIs), neutrals atoms, and very energetic cosmicrays. The interaction between the flows produces the heliospheric interface consisting of the termination shock (TS), the inner heliosheath, the heliopause (HP), and the outer heliosheath. The Voyager 1 and 2 are the only two space probes that had crossed both the TS and HP, and are currently providing in situ measurements of charged particles and magnetic fields in the outer heliosheath.
Theoretical modeling of the outer heliosphere has been performed for more than five decades, following the pioneering work by Parker (1961) and Baranov et al. (1971). Since the early influential work by Baranov & Malama (1993), computer models have greatly improved, helped by the development of theoretical knowledge and accumulation of in situ and remote observations (Zank 1999;Izmodenov & Baranov 2006). Benefitting from the progress in computers, numerical magnetohydrodynamic multicomponent models of the outer heliosphere were developed that treated the neutral atoms as a separate component coupled with the plasma via charge exchange (Linde et al. 1998;Pogorelov et al. 2004;Izmodenov et al. 2005;Opher et al. 2006;Washimi et al. 2007).
Cosmic-rays are thought to have an effect on the flow of plasma inside the heliosphere (Ko & Webb 1988;le Roux & Ptuskin 1995;Florinski & Jokipii 1999;Fahr et al. 2000). Myasnikov et al. (2000b) and Myasnikov et al. (2000a) showed that galactic cosmic-rays had an insignificant effect on heliospheric plasma flows compared with H atoms. The anomalous cosmic-rays (ACRs), produced in the heliosphere, appear to have a stronger influence on the structure of the outer heliosphere. Their pressure leads to a formation of a smooth precursor upstream of the TS and a change in the radial distance to the gasdynamic subshock (Fahr et al. 2000;Florinski et al. 2004). However, no apparent changes were seen in the position of the HP for the case of uniform diffusion coefficients in both solar wind and interstellar regions (Alexashov et al. 2004).
Recent Voyager 1 observations showed that some ACRs can cross the HP and escape into the LISM Florinski et al. 2015), where the interstellar magnetic field is nearly laminar below the scale of several au (Burlaga et al. 2015). It can be inferred from these observations that cosmicrays diffuse much faster in the LISM than in the solar wind. Therefore, it is not appropriate to impose a uniform diffusion coefficient throughout the heliosphere and LISM as was done in Alexashov et al. (2004). In order to overcome this deficiency, we conducted a numerical simulation including ACRs, in which the plasma and the ACRs are decoupled in the interstellar medium (Guo et al. 2018). However, the diffusion coefficient is equivalent to zero at HP, and might have negative impact on the ACR transport across the boundary.
Similar to the previous work (Guo et al. 2018), we have performed numerical simulations to investigate the effects of ACRs on the heliosphere using an MHD-neutrals-cosmic-ray coupled model, where the cosmic-rays are diffusive and are governed by a transport equation. The main difference between the model of Alexashov et al. (2004) and this model is the inclusion of both the interplanetary and interstellar magnetic fields and the use of a variable diffusion coefficient with rapid changes across the TS and HP. Compared with Guo et al. (2018), we introduced a new method to model the strongly diffusive environment of the LISM. We also compared the model derived flow parameters with Voyager 1 observations in the inner and outer heliosheaths, and discuss their similarities.

Numerical Models
We use a multifluid approach to represent the heliosphere consisting of plasma, neutrals, and cosmic-rays. For the plasma component the ideal MHD equations are used, while the neutrals are modeled via gasdynamic equations. The timedependent equations are integrated until a steady state is achieved. Cosmic-rays are treated as a single massless "fluid" that is entirely described via its pressure P e . This is an efficient approach when the details of the energy distribution are not of interest (e.g., Fahr et al. 2000;Myasnikov et al. 2000b). The combined equations can be written as where ρ, u, P, and B represent the plasma density, velocity, thermal pressure, and magnetic field, respectively. The plasma energy density E=P/(γ−1)+ρu 2 /2+B 2 /2 and the total pressure, excluding the cosmic-ray component, P T =P+B 2 /2. Likewise, E c =P c /(γ c −1) is the energy density of cosmic-rays, and κ is the isotropic energy-averaged diffusion coefficient. The adiabatic index γ=γ c =5/3 assuming nonrelativistic particles. The source terms Q N , Q M , and Q E represent the charge exchange of number, momentum, and energy between the plasma and neutral atoms (Pauls et al. 1995; see the Appendix A for more details). This model features four neutral populations which originate from the undisturbed interstellar medium, the outer heliosheath, the inner heliosheath, and the supersonic solar wind, respectively. To model the neutrals, the values of B and P c are set to zero in the above equations. Here we do not distinguish PUIs from solar wind protons (e.g., Wang & Richardson 2001), and assume that the PUIs are thermally assimilated into the solar wind. The quantity α is the local injection rate of ACR particles from the lower energy PUIs, which is simply assumed to be proportional to the thermal pressure, α=α′P, where α′ is a free parameter measuring the injection efficiency (Zank et al. 1993;Fahr et al. 2000;Alexashov et al. 2004) from thermal plasma to ACRs by means of the shock acceleration mechanism near the TS. Although some researchers proposed that α′ along the TS surface are affected by a change in the angle between the magnetic field and the normal to the shock (Chalov 1993); here we simply assume it to be uniform on the TS surface. We denote α′ as α in the following content for simplicity. There are six fluids in the coupled system in total (plasma, cosmic-rays, and four neutral populations). The system was solved numerically on a level 5 three-dimensional spherical geodesic mesh (Florinski et al. 2013) with 2562 hexagonal cells on the sphere and 648 concentric layers in the radial direction. The radial resolution of the mesh near TS and HP was about 0.5 and 0.7 au, respectively.
The inner boundary was a 10 au sphere centering at the Sun, and the outer boundary was set at 800 au. The solar wind conditions were calculated from the Parker solution (Parker 1958) based on the typical solar wind conditions at 1 au, which are chosen as follows: number density of protons n p =5.5 cm −3 , flow speed u p =420 km s −1 , temperature T p =64,000 K, and radial component of interplanetary magnetic field B r =2.45 nT. For the unperturbed LISM, we used the proton number density km s −1 , and its direction is given by the neutral He flow with (λ, β) He =(79°,−4°.98), where λ is the ecliptic longitude and β is the ecliptic latitude. The LISM hydrogen flow is deflected from its original direction (which is same as the He flow) by the interaction with the ionized component near the heliosphere and is observed at (λ, β) H =(72°.5,−8°.9) (Lallement et al. 2010). These two vectors form the hydrogen deflection plane (HDP). We assumed that the interstellar magnetic field B LISM bisects the angle between the direction of the IBEX ribbon center, (λ, β) ribbon =(219°.2, 39°.9) (Funsten et al. 2013), and its projection on the HDP. Figure 1 shows the magnetic field comparison along V1 direction between the simulation (blue solid line) and the observation (red solid line), including magnitude B, azimuthal angle λ, and elevation angle δ. Here we do not consider the effects of cosmic-rays, and the simulated HP is located at a radial distance of ∼133 au, i.e., farther away than observed. Note that the HP is identified numerically by the interface of the passive tracers, which have interstellar medium and solar wind origins (Florinski et al. 2013). For better comparison, the simulated HP is shifted to coincide with the observed 121.6 au, which is indicated by the vertical dashed line. The observations showed that the magnitude B and the azimuthal angle λ decreased gradually beyond the HP, while the elevation angle δ had an increasing trend from the year 2012.65 to 2017.30 (roughly 121.6-140 au; Burlaga et al. 2018). These trends of the two angles differed from those reported earlier using a shorter observational interval , and featured several numerical investigations (Isenberg et al. 2015;Zirnstein et al. 2016). The present results show that the magnitude B and the elevation angle δ have similar trends to the observation. However, the simulated elevation angle λ increased gradually, which is opposite to what was observed. Due to the relatively coarse grid spacing, the HP appears to be a diffusive transition region. The HP moves inward once the cosmic-rays are included in the simulation, as demonstrated in the next section.
The cosmic-ray pressure P c includes only ACRs, and the lower energy limit is assumed to ∼5 keV, so that the diffusive high energy PUIs in the solar wind are treated as ACRs as well (Decker et al. 2015). A zero gradient boundary condition for ACRs is used at the inner boundary and their pressure is set to zero at the outer boundary. Unlike the previous treatment of Guo et al. (2018), where the coupling between plasma and ACRs was switched off beyond the HP, here we instead increase the diffusion coefficient in the LISM to a much higher value compared to that in the solar wind. For simplicity and comparison, we set three scenarios for the energy-averaged diffusion coefficients in the solar wind: κ 1,2 =5.0×10 20 cm 2 s −1 (S1), κ 1,2 =5.0×10 21 cm 2 s −1 (S2), and κ 1 =5.0×10 21 cm 2 s −1 , κ 2 =5.0×10 20 cm 2 s −1 (S3). The subscripts 1 and 2 denote the unshocked solar wind region and the inner heliosheath, respectively. Because the solar wind is compressed and the magnetic field strength increased significantly in the inner heliosheath after the TS crossing (e.g., Richardson et al. 2013), the diffusion coefficients for cosmic-rays are expected to be lower than those in the supersonic solar wind. Therefore, S3 seems to be more realistic than the other two scenarios. In the LISM, the diffusion coefficient is expected to be extremely large for ACRs because of the weak turbulence (Burlaga et al. 2015), e.g., of order of ∼10 25 cm 2 s −1 as discussed in the previous paper (Guo et al. 2018). The very fast diffusion makes the explicit simulation of the Equation (5) nearly impossible because the time step for the numerical evolution becomes extremely small. Here we simply fix it to be 5.0×10 22 cm 2 s −1 , which is much smaller than the theoretical expectation, but still provide a fast diffusion for ACRs in LISM. This treatment is expected to still provide the trends for ACR pressure decrease beyond the HP. Detailed treatment of the diffusion coefficients can be found in the Appendix B. Figure 2 show the distribution of cosmic-ray pressure P c in the meridional (top) and equatorial (bottom) planes. From left to right, the three panels correspond to the three scenarios (S1-S3), respectively. Here the PUI injection efficient α was set to 1.0, which meant that a large amount of PUIs were injected and accelerated into ACRs. Injection took effect mainly near the TS where the plasma is compressed and ∇·u is negative. The ACRs are assumed to be accelerated from the lower energy PUIs through the mechanism of diffusive shock acceleration at the TS. Here the lower energy PUIs (5 keV) are included with the bulk thermal plasma. We use the term PUIs for this lowenergy component to distinguish it from accelerated PUIs, which are assumed to be diffusive and is included with the ACRs. After being accelerated, the newly born ACRs are convected with the shocked solar wind and simultaneously diffuse off the solar wind streamlines.

Numerical Results
From Figure 2 it is seen that the ACRs mainly reside in the inner heliosheath, and accumulate in the upwind direction. For S1, the rapidly expanding solar wind denies ACRs access to the inner heliosphere because of the small value of the diffusion coefficient. For the high diffusion scenarios S2 and S3, the socalled shock precursor develops in front of the TS, so that the shock has a broad transition and is not as sharply defined as in S1. Voyager 2 observations upstream of the TS appear to confirm the scenario where ACRs are efficiently back-scattered to the upstream of TS, where their effect was to slow down the flow of the solar wind (e.g., Florinski et al. 2009).
The same diffusion coefficients were set for S1 and S3 in the inner heliosheath, which is only one tenth of that for S2. This indicates that the ACRs for S2 will diffuse much faster than the other two scenarios. As we can see in Figures 2(B) and (E), the ACR pressure for S2 is lowest among all the three scenarios and the ACRs are not sufficiently accumulated in the inner heliosheath. This occurs mainly because of the less efficient acceleration at the shock and faster diffusion across the HP. The S1 and S3 distributions of ACR pressure look similar, except that in the former, TS is much sharper than in the latter. Note that the sudden drops of ACR pressure across HP for S1 and S3 were similar to the behavior of the measured intensities of ACR particles during the HP crossing by V1 Stone et al. 2013).
We now examine the distributions of plasma properties and cosmic-ray pressure in more detail. Figure 3 shows the profiles of the radial speed u r of plasma along the V1 direction for S1-S3. The profiles along V2 are not plotted here because they are similar to that of V1. The vertical black dotted and dashed lines indicate the positions of the TS and HP for the case with no ACRs (∼83 au and ∼133 au, respectively). The simulated TS crossing distance is closer to 85 au (Stone et al. 2008), than to 94 au for V1 (Decker et al. 2005). The HP location is within the range 130-157 au reported by other modeling groups (Linde et al. 1998;Pogorelov et al. 2004;Izmodenov et al. 2005;Opher et al. 2006). These simulated distances to the HP are significantly larger than the observed value of 121.6 au .
We next took the ACR effect into account by injecting the low-energy PUIs at the TS. The PUI injection efficient α is set to be 0.1, 0.5, and 1.0 in each scenario, which indicated that an increasing amount of PUIs were injected and accelerated into the ACRs. Different colors are used to represent the cases of the three α. The TS positions are not shown in Figure 3 for S1-S3, but can be inferred from the sharp jumps of the plasma speed. For the low diffusion scenario S1, the HP moves inward from ∼133 au to ∼128 au in response to increasing α. The TS does not move as much as the HP, although its effective width increases due to a precursor formation. The typical width of the shock precursor in these simulations is ∼1−5 au.
For S2, the subshock of the TS moves outward distinctly for α0.5, with a radial displacement of 5 au, as panel B shows. This response of TS is similar to the TS variation in the previous work (Alexashov et al. 2004), in which the κ increased from 3.75×10 20 cm 2 s −1 to 3.75×10 21 cm 2 s −1 , and a lower injection rate of ACRs α=0.1 is used. The largest displacement of the TS was ∼8 au for α=1.0, in which case the TS was located at ∼92 au. The precursor extended farther away from the subshock with a size of 24 au. Simultaneously, the HP moved inward significantly with a maximum displacement of ∼10 au, where it was located at ∼123 au. These distances are a close match to the distances to the TS and HP crossings by V1 (94 au and 122 au), respectively. The difference between the distances to the TS during its V1 and V2 crossings were attributed, in part, to different diffusion environments during the different phases of the solar cycle (Guo et al. 2018), but we do not assert that point of view here.
In the S3 case, the diffusion coefficient in the upstream solar wind is the same as in S2. The response of the TS to the different α is similar to S2, except the latter has a larger TS displacement with an increasing PUI injection, as seen in panels B and C of Figure 3. The narrowest heliosheath is produced for the case of α=1.0. However, the change in the distance to the HP is much less than for S2, and similar to S1. Because the diffusion coefficients in the inner heliosheath are the same for S1 and S3, this result indicates that the diffusion in the shocked solar wind plays an important role in controlling the distance to the HP. Table 1 presents the distances to the HP along the V1 direction for the three scenarios under different PUI injection conditions. We can see that unless a combination of large PUI injection and highly diffusive environment occurs simultaneously (e.g., α=1.0 for S2), the HP displacement is negligible compared with the typical scale of the outer heliosphere.
In order to demonstrate the behavior of the ACRs in the three scenarios we plot several pressure distributions in Figure 4. These consist of the "thermal" pressure, dominated by PUIs (P), the ACR pressure (P c ), the combined pressure (P+P c ), Figure 2. Cosmic-ray pressure P c at the meridional (top) and equatorial (bottom) planes for the three scenarios in the heliosphere: S1 (A, D), S2 (B, E), and S2 (C, F). The projections of V1 and V2 trajectories onto the meridional plane are shown with red lines. and the ratio P c /(P+P c ) along the V1 direction for S1. In the inner heliosheath, the maximum thermal pressure decreases from 1.4 to 0.7 pdyne cm −2 when the PUIs are injected at α=0.1 due to a transfer of energy to the diffusive component. Increasing α from 0.1 to 1.0 reduces the thermal pressure accordingly, eventually reaching a very low 0.3 pdyne cm −2 for the case of α=1.0. At the same time, the maximum ACR pressure increases from 0.7 to 1.1 pdyne cm −2 with increasing α. In each scenario the ACR pressure gradually increases along the radial direction after TS crossing until ∼102 au in the inner heliosheath, where it begins to drop finally reaching a very low value at the HP (∼128-133 au, as shown in Table 1), owing to fast diffusion beyond the HP. The trends are qualitatively consistent with the observed behavior of ACRs , except that the simulated ACR pressure drops earlier than in the observation.
On the interstellar side, the ACR pressure decreases to a very small amount but does not disappear completely, due to a relatively small diffusion coefficient used in the model. For the ACRs with the energy range 3.4-17.3 MeV, the observed fluxes disappeared completely only approximately a month after the HP crossing , which give an estimation of the penetration length of less than 0.5 au. Here we do not have similar signatures, as shown in the panel B, the pressure P c becomes extremely small and nearly constant right after the HP crossing. In panel C, the total pressure is nearly conserved in the inner heliosheath, except near the HP, where the ACRs diffuse very fast on the LISM side resulting in a depletion of the total pressure because of the energy and momentum loss by the plasma and ACRs. In order to justify the importance of ACR pressure, we plotted the ratio of P c /(P+P c ) as well shown in panel D. It seems that a higher injection efficiency leads to a larger ratio in the sheath, where the ratio gradually increases with radial distance and then decreases near the HP. The value of P c /(P+P c ) varied from 0.2 to 0.9 in this scenario, depending on the PUI injection rate.  The diffusion coefficient in S2 is larger than in S1 by a factor of ten. As seen in Figure 5(A), the thermal pressure P changes significantly in the inner heliosheath, and the distributions are both higher and narrower than in S1 because a higher diffusion corresponds to less ACR acceleration and a wider precursor followed by a subshock. Consequently, the TS expands outward and the HP moves inward. In panel B, the ACR pressure peaks at the TS and then decreases dramatically along the radial direction until reaching the HP, after which point the pressure continues to decrease but at a slower rate. This behavior is inconsistent with the observed intensities of the high energy ACRs which continued to increase well into the heliosheath after the TS crossing (Stone et al. 2008), but similar to the long-term variation of the low-energy ACRs (Decker et al. 2015). Compared with S1, the ACR pressure is much lower due to the faster diffusion of the energetic particles, because a fast diffusion indicates a less efficient acceleration from PUIs to cosmic-rays (Drury 1983), and a more efficient transport of accelerated particles to the other regions from the TS. In panel C, the total pressure P+P c occupies a narrower region with increasing PUI injections because of the reduction in the width of the inner heliosheath. The ratio of P c /(P+P c ) shown in panel D varies from ∼0.05 to 0.45 as more and more PUIs are injected and accelerated into ACRs.
From a physical perspective, S3 is more realistic than either S1 or S2 because the diffusion in the inner heliosheath is expected to be slower than that in the supersonic solar wind. Figure 6 presents the pressure distributions along V1 direction for S3. Compared to S1, the thermal pressure P is a little larger because of a less efficient acceleration from PUIs to ACRs in S3. On the contrary, the ACR pressure P c appears to be lower than in S1. The low PUI acceleration at the TS is simply related to the fact that the diffusion coefficient is ten times of that of S1 in the unshocked solar wind. In panel B, similar to that of S1, the ACR pressure increases after the TS crossing, reaching a low level near the HP. In panel C, the TS precursor is as pronounced as that of S2, and much wider than that of S1. The ratio of P c /(P+P c ) is similar to that of S1 (0.2-0.9, depending on the PUI injection rate).

Discussions
Gloeckler & Fisk (2016) have estimated the partial pressures of every dynamically important particle population in the heliosheath and LISM using the energetic neutral hydrogen spectra from IBEX (Galli et al. 2016), Cassini (Dialynas et al. 2013), and SOHO (Czechowski et al. 2008), as well as Voyager charged particle and magnetic field data. The pressure in the inner heliosheath is dominated by that of PUIs Figure 4. Profiles of thermal pressure P, cosmic-ray pressure P c , and the total pressure P+P c along V1 direction in panels A, B, and C, respectively, for the case of S1. The unit is pdyne cm −2 . The panel D shows the corresponding ratio of P c /(P+P c ).
(∼3.0 pdyne cm −2 ), while the ACR pressure with energies > 40 keV measured by Voyagers is relatively small (∼0.23 pdyne cm −2 ), giving the ratio of P c /(P+P c )∼0.07. However, the distinction between accelerated PUIs and ACRs is somewhat arbitrary. If the suprathermal tails of PUIs (> 5 keV) are counted as a diffusive population (i.e., ACRs), the pressure ratio could be as high as ∼0.4−0.5. For this reason we argue that our simulation results with α=0.1-0.5 in S1 (or S3), and α=1.0 in S2, are consistent with the observations mentioned.
The pressure of the accelerated PUIs is dominated by the particles in the 5-40 keV range that are not measured by Voyagers. The partial pressures of 0.028-3.5 MeV (V2) and 0.04-4.0 MeV (V1) have a decreasing tendency with the radial distance that might be related to the solar cycle (Decker et al. 2015). On the contrary, the intensities of the high energy ACRs (12-22 MeV) detected by the two probes continues to increase after the TS crossing (Stone et al. 2008). The present model does not distinguish particle by energy, and it is not possible to determine the average energy in a meaningful way without the knowledge of the gradients at different energies. For this reason we believe that the simulation results from at least one of the scenarios may correspond to the real situation. From Table 1, the largest displacement of the HP (∼10 au) occurs for S2 in the case of a high PUI injection efficiency α=1.0. For the other two scenarios, the HP does not move by much and the ACR effect is not important.
The observed two HP crossings at 122 au by the Voyager 1 and 119 au by Voyager 2 occurred earlier than predicted by global MHD simulation . It has been proposed that the Rayleigh-Taylor type instabilities triggered by the solar wind variation make the LISM penetrate deep inside the solar wind, thus locally reducing the apparent width of the heliosheath (Borovikov & Pogorelov 2014). Some researchers suggested adding a physical effect such as thermal conduction in the MHD equation, will decrease the thermal pressure in the heliosheath, and thus lead to a thinning of the inner heliosheath thickness (Izmodenov et al. 2014). From the above simulation results, we argue that these observations could also be interpreted as the heliosheath that is narrow in the upwind direction because of the momentum and energy transfer from the thermal plasma to the diffusive ACRs. Two conditions must be satisfied for this loss of energy by the solar wind plasma to be significant: (a) the large intensity of the lowenergy ACRs at the TS, and (b) fast diffusion of ACRs (including high energy PUIs) through the inner heliosheath.
The TS precursor observed by V2 had a size of ∼0.7 au (Richardson et al. 2008), which is close to the spatial resolution of ∼0.5 au near TS in the simulation. A finer grid spacing is expected to be better for the capture of TS precursor. Unfortunately, in the current explicit numerical scheme, a finer spatial resolution (e.g., 0.1 au near TS) can not be achieved in the simulation, which is highly time-consuming because the time step becomes extremely small due to the high diffusion in the interstellar medium. It will be valuable to investigate the TS precursor under a finer spatial resolution in the future.

Summary
In this work, we have conducted a combined global MHDneutrals-cosmic-rays simulation of the outer heliosphere using the typical solar wind and interstellar medium conditions. The cosmic-ray component was treated as a diffusive and massless fluid, so that their pressure was the only variable pertaining to the energetic particles included in the MHD equations and evolved in the time-dependent manner. The extremely high diffusion environment of the outer heliosphere was treated in a simplified fashion by assuming that the diffusion coefficient is 5.0×10 22 cm 2 s −1 , which is a factor of 10 or 100 larger than in the solar wind.
Three scenarios were constructed by adjusting the diffusive properties of the cosmic-rays. The results show that the ACRs could have a significant effect on the geometry of the outer heliosphere for a sufficiently large PUI injection efficiency and diffusion coefficient (κ=5.0×10 21 cm 2 s −1 ) in the solar wind. Compared to the case without cosmic-rays, the width of the inner heliosheath in the upwind direction is dramatically reduced because of the expansion of the TS (∼8 au) and the shrinking of the HP (∼10 au). This result may be used to interpret, at least in part, the early HP encounter by both Voyagers. For the case of a smaller solar wind diffusion coefficient (κ=5.0×10 20 cm 2 s −1 ) in the inner heliosheath, ACR effects on the HP are negligible amounting to a change in the distance to the HP by only 1-3 au.
The work of X.G. and C.W. was supported by NNSFC grants 41874171, 41674146, 41731070, and 41574159, in part by the Specialized Research Fund for State Key Laboratories of China and NSSC research fund for key development directions, and by the Strategic Pioneer Program on Space Science of CAS grant XDA15052500, XDA17010301, QYZDJ-SSW-JSC028. The work of V.F. was partially supported through NASA grant NNX17AB85G and 80SSC18K1209.

Appendix A Source Terms of Charge-exchanges
The interstellar neutrals are divided into the four populations according to their origins, the undisturbed interstellar medium (region 1), the outer heliosheath (region 2), the inner heliosheath (region 3), and the supersonic solar wind (region 4). The four neutrals couple with the plasma independently,