THREE-DIMENSIONAL KINETIC-MHD MODEL OF THE GLOBAL HELIOSPHERE WITH THE HELIOPAUSE-SURFACE FITTING

Global modeling of the solar wind (SW) interaction with the local interstellar medium (LISM) is essential: (1) in general, to understand how the Sun acts as a star and exchanges matter and energy with LISM, and (2) in particular, to explain relevant space experiment data (e.g., obtained on board both Voyager spacecraft, the Interstellar Boundary Explorer (IBEX), the Solar and Heliospheric Observatory (SOHO), the Hubble Space Telescope (HST), etc.).

Following the pioneering paper by Baranov & Malama (1993), self-consistent kinetic-gas dynamics models of the SW/LISM interaction were developed by our group in Moscow (e.g., see Izmodenov & Baranov 2006; Izmodenov 2009 for reviews). There are two features which make our model approach different from similar models developed by other groups. These features are as follow: (1) a rigorous kinetic treatment of the interstellar H atom component using a Monte-Carlo method with splitting of the trajectories, and (2) a Godunov-type numerical method used to solve the ideal gas dynamical or MHD equations. The method has the potential to fit some major discontinuities of the flow and to capture others.

The kinetic treatment of H atoms is essential because their mean free path is comparable to the size of the SW/LISM interaction region (see, e.g., Izmodenov 2000). The difference between kinetic and alternative multi-component models has been explored in detail by Alexashov & Izmodenov (2005) and Müller et al. (2008). At heliospheric distances smaller than 10–20 AU, solar gravitation and radiation forces make kinetic effects even more important. The velocity distribution function of H atoms is essentially non-Maxwellian (Izmodenov 2001; Izmodenov et al. 2013) in this region. This is very important for analyses of any measurements (direct or indirect) of H atoms.

To obtain an effective tool to calculate the H atom space and velocity distribution inside the heliosphere, Katushkina & Izmodenov (2010) suggested a two-step procedure. In the first step, the H atom distribution is calculated in the frame of a global self-consistent model of the SW/LISM interaction. As a particular result of this model, we obtain the velocity distribution of H atoms at a certain heliocentric distance (e.g., 90 AU) where local solar effects could be neglected. This distribution function is used as a boundary condition for the second step in which the velocity distribution inside the heliosphere is calculated with the resolution (in time, space, and velocity space) required for the data analyses. The kinetic equation at this step is solved using the method of characteristics in 3D time-dependent case (Izmodenov et al. 2013).

The two-step procedure described above has been used in the companion paper by Katushkina et al. (2015) in order to analyze the IBEX-Lo data for the interstellar hydrogen fluxes. The aim of this paper is to describe the global model of the SW/LISM interaction that has been employed in the companion paper during the first step of the two-step procedure. The model presented is a currently state-of-the-art model which we developed. It has not yet been described in the literature.

The original Baranov & Malama (1993) model is a stationary axisymmetric kinetic-gas dynamical model. This model has been advanced by considering additional physical components—galactic and anomalous cosmic rays (Myasnikov et al. 2000; Alexashov et al. 2004), interstellar helium ions and SW alpha particles (Izmodenov et al. 2003), and interstellar oxygen and nitrogen (Izmodenov et al. 2004). This model has also been extended in the tail direction to distances of ∼20000–30000 AU (Izmodenov & Alexashov 2003). In 2005, the model became time-dependent (Izmodenov et al. 2005b), and then realistic SW data were used in the model (Izmodenov et al. 2008). In 2005, the stationary version of this model was expanded to three dimensions and the interstellar magnetic field (IsMF) was taken into account (Izmodenov et al. 2005a). This model allowed us to explain quantitatively the difference in the directions of the interstellar hydrogen and helium flows inside the heliosphere which had been discovered by Lallement et al. (2005, 2010). The model provided an asymmetry of the distances to the heliospheric termination shock (TS), in agreement with Voyager-1 and Voyager-2 data (Izmodenov 2009). Vincent et al. (2014) used this model to analyze the backscattered Lyα spectra obtained by HST.

The Izmodenov et al. (2005a) model has been criticized because it did not take into account the heliospheric magnetic field (HMF) and variations of the SW parameters with the heliolatitude. The importance of HMF on the plasma flow in the inner heliosheath was realized very early by Axford et al. (1963) and Nerney & Suess (1975), and later in the kinematic approach in the frame of the Baranov and Malama model by Barsky (1999). The dynamical effects of HMF have been explored by other groups in the models of Washimi (1993), Linde et al (1998), Pogorelov (2004), Opher et al. (2004), and Ratkiewicz et al. (2006). The heliolatitudinal variations of the SW on the global structure of the SW/LISM interaction region were previously studied by Tanaka & Washimi (1999), Grygorczuk et al. (2010), Kim et al. (2012), Pogorelov et al. (2013), and Provornikova et al. (2014). More recent work on the SW/LISM modeling includes the studies of Heerikhuisen et al. (2014), Zirnstein et al. (2014), Zank et al. (2013), Ben-Jaffel & Ratkiewicz (2012), Florinski et al. (2013), Opher et al. (2015), Borovikov & Pogorelov (2014).

In this paper, our 3D kinetic-MHD model is extended by taking into account the dynamic effects of HMF in the inner heliosheath region, as well as the effect of the latitudinal variations of the SW parameters. Note that the model presented here has already been used in a number of studies. For example, Wood et al. (2014) employed the model to simulate Lyα absorption spectra in the tail of the heliosphere. Katushkina et al. (2015) used the model results to analyze the deflection of the direction of the H atom velocity vector in the heliosphere from its direction in LISM. Finally, in the companion paper (Katushkina et al. 2015), this model has been used for the analysis of IBEX-Lo data.

Despite the fact that the model has already been employed for data analyses, a detailed description of the model had not been presented in the literature, although it had been reported at several conferences. The aim of this paper is to fill this gap and present a detailed formulation of the model (Section 2) and its results (Section 3). The effects of the HMF and latitudinal variations of the SW on the global structure of the SW/LISM interaction region will be explored explicitly by comparisons with models where the effects are not taken into account. To our knowledge, this has not yet been done. The relevance/irrelevance of the model results to observational data is discussed in Section 4.

In the model, partially ionized interstellar plasma is considered as a two-component gas. The two components are (1) a neutral component consisting of atomic hydrogen and (2) a charged (or plasma) component consisting of protons, electrons, and helium ions. In the SW, the plasma component consists of protons, alpha particles, and electrons.

The neutral component is described kinetically by means of a velocity distribution function, ${f}_{{\rm{H}}}\left({\boldsymbol{r}},{{\boldsymbol{w}}}_{{\rm{H}}}\right),$ that is determined from the solution of the following kinetic equation:

$$\begin{eqnarray}&&{{\boldsymbol{w}}}_{{\rm{H}}}\cdot \displaystyle \frac{\partial {f}_{{\rm{H}}}}{\partial {\boldsymbol{r}}}+\displaystyle \frac{{{\boldsymbol{F}}}_{{\rm{r}}}+{{\boldsymbol{F}}}_{{\rm{g}}}}{{m}_{{\rm{H}}}}\cdot \displaystyle \frac{\partial {f}_{{\rm{H}}}}{\partial {{\boldsymbol{w}}}_{{\rm{H}}}}={-\nu }_{\mathrm{ph}}{f}_{{\rm{H}}}\left({\boldsymbol{r}},{{\boldsymbol{w}}}_{{\rm{H}}}\right)\\ &&-{f}_{{\rm{H}}}\cdot \displaystyle \int | {{\boldsymbol{w}}}_{{\rm{H}}}-{{\boldsymbol{w}}}_{{\rm{p}}}| {\sigma }_{\mathrm{ex}}^{\mathrm{HP}}(| {{\boldsymbol{w}}}_{{\rm{H}}}-{{\boldsymbol{w}}}_{{\rm{p}}}| ){f}_{{\rm{p}}}\left({\boldsymbol{r}},{{\boldsymbol{w}}}_{{\rm{p}}}\right)d{{\boldsymbol{w}}}_{{\rm{p}}}\\ &&+{f}_{{\rm{p}}}\left({\boldsymbol{r}},{{\boldsymbol{w}}}_{{\rm{H}}}\right)\displaystyle \int | {{\boldsymbol{w}}}_{{\rm{H}}}^{*}-{{\boldsymbol{w}}}_{{\rm{H}}}| {\sigma }_{\mathrm{ex}}^{\mathrm{HP}}(| {{\boldsymbol{w}}}_{{\rm{H}}}^{*}-{{\boldsymbol{w}}}_{{\rm{H}}}| ){f}_{{\rm{H}}}\left({\boldsymbol{r}},{{\boldsymbol{w}}}_{{\rm{H}}}^{*}\right)d{{\boldsymbol{w}}}_{{\rm{H}}}^{*}.\end{eqnarray} \tag{ 1 }$$

Here, ${{\boldsymbol{F}}}_{{\rm{r}}}$ and ${{\boldsymbol{F}}}_{{\rm{g}}}$ are the forces of solar radiation pressure and solar gravitation, respectively. It is assumed in the present calculations that $| {{\boldsymbol{F}}}_{{\rm{r}}}| /| {{\boldsymbol{F}}}_{{\rm{g}}}| =1.258.$ Despite the fact that this parameter (and its temporal variations) is very important for the distribution of H atoms in the heliosphere inside a sphere of 10–20 AU (e.g., Katushkina et al. 2015), it does not influence the global structure of the SW/LISM region discussed in this paper. ${f}_{{\rm{p}}}\left({\boldsymbol{r}},{{\boldsymbol{w}}}_{{\rm{H}}}\right)$ is the local Maxwellian distribution function of protons with gas dynamic values of $\rho \left({\boldsymbol{r}}\right),$ ${\boldsymbol{V}}\left({\boldsymbol{r}}\right)$, and $T\left({\boldsymbol{r}}\right).$ ${\sigma }_{\mathrm{ex}}^{\mathrm{HP}}(u)$ is the effective charge-exchange cross section and ${\sigma }_{\mathrm{ex}}^{\mathrm{HP}}(u)={\left(2.2835\cdot {10}^{-7}-1.062\cdot {10}^{-8}\mathrm{ln}(u)\right)}^{2}$ cm² (Lindsay & Stebbings 2005), where u is the relative atom-proton velocity in cm s⁻¹. ${\nu }_{\mathrm{ph}}=1.67\cdot {10}^{-7}{({R}_{{\rm{E}}}/R)}^{2}\;{{\rm{s}}}^{-1}$, (${R}_{{\rm{E}}}=1\;\mathrm{AU}$) is the photoionization rate.

The charged (plasma) component is described in the context of an ideal MHD approach:

$$\begin{eqnarray}&&{\rm{\nabla }}\cdot (\rho {\boldsymbol{V}})={q}_{1},\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&{\rm{\nabla }}\cdot \left[\rho {\boldsymbol{VV}}+\left(p+\displaystyle \frac{{B}^{2}}{8\pi }\right){\boldsymbol{I}}-\displaystyle \frac{{\boldsymbol{BB}}}{4\pi }\right]={{\boldsymbol{q}}}_{2},\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&{\rm{\nabla }}\cdot ({\boldsymbol{VB}}-{\boldsymbol{BV}})=0,\;\;{\rm{\nabla }}\cdot {\boldsymbol{B}}=0,\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}&&{\rm{\nabla }}\cdot \left[\left(E+p+\displaystyle \frac{{B}^{2}}{8\pi }\right){\boldsymbol{V}}-\displaystyle \frac{({\boldsymbol{V}}\cdot {\boldsymbol{B}})}{4\pi }{\boldsymbol{B}}\right]={q}_{3},\end{eqnarray} \tag{ 5 }$$

where ${\boldsymbol{B}}$ is the magnetic field induction vector, ${\boldsymbol{a}}$ ${\boldsymbol{b}}$ is the tensor production of two vectors ${\boldsymbol{a}}$ and ${\boldsymbol{b}}$, ${\boldsymbol{I}}$ is the unity tensor, "·" is the scalar product, $E=\displaystyle \frac{\rho {V}^{2}}{2}+\frac{p}{\gamma -1}+\frac{{B}^{2}}{8\pi },$ and γ = 5/3. The total plasma density is determined as $\rho ={m}_{{\rm{p}}}{n}_{{\rm{p}}}+{m}_{\mathrm{He}}{n}_{\mathrm{He}},$ where n_He denotes the He⁺ number density in the interstellar medium and the He⁺⁺ number density in the SW. To determine the number densities, we solve the continuity equations for the helium ion component in the LISM and for alpha particles in the SW. Then, the proton number density can be obtained as ${n}_{{\rm{p}}}=(\rho -{m}_{\mathrm{He}}{n}_{\mathrm{He}})/{m}_{{\rm{p}}}.$ p is the total plasma pressure, defined as $p=(2{n}_{{\rm{p}}}+3{n}_{{\mathrm{He}}^{++}}){k}_{{\rm{B}}}{T}_{{\rm{p}}}$ in the SW and as $p=2({n}_{{\rm{p}}}+{n}_{{\mathrm{He}}^{+}}){k}_{{\rm{B}}}{T}_{{\rm{p}}}$ in the LISM. T_p is the plasma temperature and k_B is the Boltzmann constant.

The influence of charge exchange with the interstellar H atoms was taken into account on the right-hand side of the MHD equations by the source terms q₁, ${{\boldsymbol{q}}}_{2}$, and q₃. The sources are calculated as the following integrals of the H-atom velocity distribution function:

$$\begin{eqnarray}&&{q}_{1}={m}_{{\rm{p}}}{n}_{{\rm{H}}}\cdot {\nu }_{\mathrm{ph}},{n}_{{\rm{H}}}\left({\boldsymbol{r}},t\right)=\int {f}_{{\rm{H}}}\left(r,{w}_{{\rm{H}}}\right){{dw}}_{{\rm{H}}},\end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray}&&{{\boldsymbol{q}}}_{2}=\int {m}_{{\rm{p}}}{\nu }_{\mathrm{ph}}{{\boldsymbol{w}}}_{{\rm{H}}}{f}_{{\rm{H}}}\left({\boldsymbol{r}},{{\boldsymbol{w}}}_{{\rm{H}}}\right)d{{\boldsymbol{w}}}_{{\rm{H}}}\;+\end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray*}&&\int \int {m}_{{\rm{p}}}{v}_{\mathrm{rel}}{\sigma }_{\mathrm{ex}}^{\mathrm{HP}}\left({v}_{\mathrm{rel}}\right)\left({{\boldsymbol{w}}}_{{\rm{H}}}-{\boldsymbol{w}}\right){f}_{{\rm{H}}}\left({\boldsymbol{r}},{{\boldsymbol{w}}}_{{\rm{H}}},t\right){f}_{{\rm{p}}}\left({\boldsymbol{r}},{\boldsymbol{w}}\right)d{{\boldsymbol{w}}}_{{\rm{H}}}d{\boldsymbol{w}},\end{eqnarray*}$$

$$\begin{eqnarray}&&{q}_{3}=\int {m}_{{\rm{p}}}{\nu }_{\mathrm{ph}}\displaystyle \frac{{w}_{{\rm{H}}}^{2}}{2}{f}_{{\rm{H}}}\left({\boldsymbol{r}},{{\boldsymbol{w}}}_{{\rm{H}}}\right)d{{\boldsymbol{w}}}_{{\rm{H}}}\;+\end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray*}&&\displaystyle \frac{1}{2}\int \int {m}_{{\rm{p}}}{v}_{\mathrm{rel}}{\sigma }_{\mathrm{ex}}^{\mathrm{HP}}\left({v}_{\mathrm{rel}}\right)\left({w}_{{\rm{H}}}^{2}-{w}^{2}\right){f}_{{\rm{H}}}\left({\boldsymbol{r}},{{\boldsymbol{w}}}_{{\rm{H}}}\right){f}_{{\rm{p}}}\left({\boldsymbol{r}},{\boldsymbol{w}}\right)d{{\boldsymbol{w}}}_{{\rm{H}}}d{\boldsymbol{w}}\end{eqnarray*}$$

$$\begin{eqnarray*}&&+\;{n}_{{\rm{H}}}{\nu }_{\mathrm{ph}}{E}_{\mathrm{ph}}.\end{eqnarray*}$$

Here, ${v}_{\mathrm{rel}}=| {{\boldsymbol{w}}}_{{\rm{H}}}-{\boldsymbol{w}}|$ is the relative velocity of an atom and a proton and E_ph is the mean photoionization energy (4.8 eV).

Source terms (6)–(8) have been calculated using the Monte Carlo method (Malama 1991) that was used to solve kinetic Equation (1). This method assumes that the velocity distribution function of the proton component is locally Maxwellian, although it allows a generalization to any isotropic distribution function as it was done by Malama et al. (2006).

Although our model and numerical code allow us to take into account the effects of electron impact ionization that might have an influence on the inner heliosheath plasma flow (e.g., Figure 3 in Baranov & Malama 1996) and on interstellar neutral filtration in the SW/LISM interaction region (Izmodenov et al. 1999), we do not take it into account here for the following two reasons. First, the temperature of the electrons is probably overestimated in the one fluid plasma approach employed here. A more rigorous multi-component plasma approach is needed (e.g., Malama et al. 2006; Izmodenov 2009). Second, in the current paper, we will explore the effects of HMF on the plasma heliosheath and, in particular, the effect of plasma density diminishing in the vicinity of the heliopause. The electron impact (if included) would make this effect less evident.

For the purpose of this paper, we performed model calculations for three different sets of boundary conditions. These sets are chosen in order to clearly explore the effects of the HMF and the dependence of the SW parameters on heliolatitude. Model 1 corresponds to the model where neither the HMF nor the heliolatitudinal dependence of the SW are taken into account. This is simply the model that was published previously (Izmodenov et al. 2005a, 2009) but with updated LISM boundary conditions. In model 2, we take into account the HMF, while the SW is still assumed to be spherically symmetric. Model 3 takes into account both the HMF and the heliolatitudinal variations of the SW.

To finish the formulation of the problem, the inner boundary conditions in the unperturbed (by LISM) SW and the outer boundary conditions in the pristine LISM should be formulated. It is important to note that the boundary conditions considered in this paper are stationary. Therefore, this paper focuses on the stationary solutions of the SW/LISM interaction problem. Nevertheless, our numerical method can be easily generalized for time-dependent (periodical) solutions, as was done for the axisymmetrical model by Izmodenov et al. (2005b, 2008).

2.1. Interstellar Parameters

The interstellar H atom and proton number densities are assumed to be n_H,LISM = 0.14 cm⁻³ and n_p,LISM = 0.04 cm⁻³, respectively. Both densities are somewhat smaller compared to the values of n_H,LISM = 0.18 cm⁻³ and n_p,LISM = 0.06 cm⁻³ which were accepted from previous studies (see, e.g., Izmodenov 2009). The reasons for such a reduction are discussed in Section 4.

The number density of the interstellar helium ions is chosen as ${n}_{\mathrm{He}+,\mathrm{LISM}}=0.003$ cm⁻³. This corresponds to the standard universal ratio of the total H to He, $({n}_{{\rm{p}},\mathrm{LISM}}+{n}_{{\rm{H}},\mathrm{LISM}})/({n}_{{\mathrm{He}}^{+},\mathrm{LISM}}+{n}_{\mathrm{He},\mathrm{LISM}})=10$ and n_He,LISM = 0.015 cm⁻³ (Gloeckler et al. 2004; Witte 2004). Interstellar helium is taken into account in all three of the models considered in this paper.

The interstellar bulk velocity and temperature are chosen as V_LISM = 26.4 km s⁻¹ and T_LISM = 6530 K. This velocity value is consistent with the ISM vector derived from analyses of Ulysses/GAS data for interstellar helium (Witte 2004; Bzowski et al. 2014; Katushkina et al. 2014; Wood et al. 2015). The direction ${\boldsymbol{V}}$_LISM is determined by the interstellar helium flow direction from Ulysses/GAS data. This value is (longitude = −102, latitude = −511) in the heliographic inertial coordinate system (HGI 2000). The velocity vector is in agreement with new results derived from IBEX data (see, e.g., McComas et al. 2015). The temperature of T_LISM = 6530 K is slightly smaller than the newly adoped temperature of 7500 K (McComas et al. 2015).

The magnitude of the IsMF is assumed to be B_LISM = 4.4 μ Gauss and its direction has α = 20° with the direction of the interstellar bulk flow. These values were obtained in the parametric study presented in Table 2 of Izmodenov (2009) to obtain asymmetry of the TS agreeable with the distances of the Voyager 1 and Voyager 2 TS crossings.

The chosen values of the velocity and IsMF correspond to the supersonic and subalfvenic LISM. The sonic, alfvenic, and fast-magnetosonic Mach numbers are M = 2.17, A = 0.631, M_Z+,and _LISM = 0.628. M_Z+,LISM is calculated here along the LISM flow direction.

In the model, it is assumed that the plane which is determined by the vectors ${\boldsymbol{V}}$_LISM and ${\boldsymbol{B}}$_LISM (VB-plane) coincides with the plane determined by the velocity vectors of the interstellar helium, ${\boldsymbol{V}}$_He,LISM, and hydrogen, ${\boldsymbol{V}}$_H,LISM, inside the heliosphere—the so-called Hydrogen Deflection Plane (see, Lallement et al. 2005, 2010; for model justification see Izmodenov et al. 2005a, 2009). Let us define the positive normal to the plane as

$$\begin{eqnarray}&&{{\boldsymbol{e}}}_{y}=\displaystyle \frac{[{{\boldsymbol{V}}}_{\mathrm{LISM}}\times {{\boldsymbol{B}}}_{\mathrm{LISM}}]}{| [{{\boldsymbol{V}}}_{\mathrm{LISM}}\times {{\boldsymbol{B}}}_{\mathrm{LISM}}]| }=\displaystyle \frac{[{{\boldsymbol{V}}}_{\mathrm{He},\mathrm{LISM}}\times {{\boldsymbol{V}}}_{{\rm{H}},\mathrm{LISM}}]}{| [{{\boldsymbol{V}}}_{\mathrm{He},\mathrm{LISM}}\times {{\boldsymbol{V}}}_{{\rm{H}},\mathrm{LISM}}]| }.\end{eqnarray} \tag{ 9 }$$

Then, for the direction of the helium flow given above and for the hydrogen flow direction of (longitude = −466, latitude = −837) in HGI 2000 (Lallement et al. 2010), the angle between ${\boldsymbol{e}}$_y and the north solar pole direction is 1375. The plane and the angle α within this plane determine the direction of IsMF. In HGI 2000, this direction is (longitude = −1648, latitude = −1822).

2.2. SW Parameters at Earth's Orbit

The inner boundary for our model is located at 1 AU. For models 1 and 2, the SW is assumed to be spherically symmetric. Its number density and velocity are assumed to be n_p,E = 5.94 cm⁻³ and V_R,E = 432.4 km s⁻¹. These values are obtained by time-averaging the OMNI SW data. We also assume that the number density of the alpha particles He⁺⁺ is 3.5% of the proton number density.

The SW is hypersonic, and we assume that the Mach number is 6.44, which corresponds to a temperature of T_E = 188500 K. Our model results are not sensitive to this parameter.

In model 3, we take into account the heliolatitudinal variations of the SW density and speed at 1 AU obtained by time-averaging different experimental data sets. Namely, three sets of experimental data are used.

• 1.  
  In the ecliptic plane, we use data (SW density and speed) from the OMNI 2 database. The OMNI 2 data set contains hourly resolution SW magnetic field and plasma data from many spacecraft in geocentric orbit and in orbit about the L1 Lagrange point.
• 2.  
  Heliolatitudinal variations of the SW speed are taken from analysis of the interplanetary scintillation data. We use the results of Sokół et al. (2013) where one-year average latitudinal profiles of the SW speed with a resolution of 10° are obtained. Data are available from 1990 to 2011.
• 3.  
  Heliolatitudinal variations of the SW mass flux are derived from analysis of SOHO/SWAN full-sky maps of the backscattered Lyα intensities (Quémerais et al. 2006, Lallement et al. 2010; Katushkina et al. 2013). An inversion procedure (see Quémerais et al. 2006 for details) allows us to obtain the SW mass flux as a function of time and heliolatitude with a temporal resolution of approximately 1 day and an angular resolution of 10°. Data are available from 1996 to 2011.

From these data, we calculate the mass and momentum fluxes of the SW as functions of time (from 1996 to 2011) and heliolatitude. After that, we average these fluxes over time and then extract the heliolatitudinal variations of the SW density and speed from the averaged fluxes (see results in Figure 1). In the model, it is assumed that the thermal pressure of the SW at 1 AU is constant at all heliolatitudes. Therefore, corresponding heliolatitudinal variations of the SW temperature are applied (variations are about ±20%). Since the SW is hypersonic, the assumption on thermal pressure at 1 AU does not influence the results.

Zoom In Zoom Out Reset image size

For the HMF, the Parker spiral solution has been assumed at 1 AU:

$$\begin{eqnarray}&&{B}_{R}={\pm B}_{{\rm{E}}}{\left(\displaystyle \frac{{R}_{{\rm{E}}}}{R}\right)}^{2},\;{B}_{\theta }=0,\;{B}_{\varphi }={-\beta }_{s}{B}_{R}\left(\displaystyle \frac{R}{{R}_{{\rm{E}}}}\right)\mathrm{sin}\theta .\end{eqnarray} \tag{ 10 }$$

Here, (R, θ, φ) are the spherical coordinates connected with the solar equatorial plane: θ is the solar latitude, counted from north solar pole (0°) to south (180°). For R = R_E = 1 AU, we assume B_E = 37.5 μG, and for solar rotation parameter β_s = Ω R_E/V ≈ 1 for the SW velocity V = 432 km s⁻¹. Ω is the Sun angular velocity correponding to the solar rotation period of 25 days. φ is the solar longitude.

The sign of the B_R component should be chosen depending on the location of the heliospheric current sheet that is changing with time. However, this is not important for our study because the terms which are responsible for the influence of the magnetic field in the ideal MHD equations do not depend on the orientation of the magnetic field. Also, in the frame of the ideal MHD equation, we are not interested in the structure of the heliospheric current sheet assuming that it is a discontinuity where the magnetic field changes its orientation. This assumption has been confirmed so far by Voyager magnetic field measurements in the SW and in the inner heliosheath (e.g., Burlaga & Ness 2012) that demonstrate a quick change of the magnetic field orientation rather than a thick (comparing to the size of the heliosheath) transition region.

In addition, the reconnection of the HMF and IsMF at the heliopause is not possible in the frame of the ideal MHD approach. We consider the heliopause as a tangential discontinuity where B_n = 0 on both sides of the heliopause.

The solution of the problem considered does not depend on the orientations of the inner and outer magnetic fields along their magnetic field lines. Therefore, to perform numerical calculations, we can assume unipolar HMF at Earth's orbit ${B}_{R}={-B}_{{\rm{E}}}{\left(\frac{{R}_{{\rm{E}}}}{R}\right)}^{2}.$ The actual polarity of the HMF can be calculated after the self-consistent solution obtained in a similar manner as in Barsky (1999).

2.3. Coordinate System

2.4. Details of the Numerical Scheme and Computational Grid

In order to obtain a self-consistent steady-state solution of the system (1)–(8) we use the global iteration method (Baranov & Malama 1993). Non-stationary versions of MHD Equations (2)–(5) are solved using a 3D modification of the finite-volume Godunov-type scheme that employs a Harten-Lax-van Leer Discontinuity MHD Riemann solver (see, for example, Miyoshi & Kusano 2005). A Chakravarthy–Osher TVD limiter with the possibility to choose compression parameters (Chakravarthy & Osher 1985) is used to increase the resolution properties of the first-order accuracy scheme.

The solver has been adopted to a 3D moving grid with the possibility to fit discontinuities. We use the soft fitting technique proposed by Godunov et al. (1979) and applied for two-dimensional (2D) MHD flows by Myasnikov (1997). The chosen numerical scheme allows us to fit all of the major discontinuities—the heliopause, the TS, and the bow shock (BS), when the latter exists. All of the other discontinuities (if any) can be captured by the scheme.

A steady-state solution of Equations (2)–(5) is obtained using the time-relaxation method. In other words, we solve non-stationary analogs of Equations (2)–(5) with stationary boundary conditions for long periods of time when the temporal derivaties of all of the parameters become negligiblely small and the flow pattern reaches steady state in the computational domain. It also has been checked that the obtained steady-state solution remains in a steady state over very long periods of time, and so possible instabilities do not destroy the numerical solution.

To satisfy the condition ∇ · ${\boldsymbol{B}}$ = 0 in our numerical solution of the stationary problem, we follow the well-known procedure suggested by Powell et al. (1999), which consists of adding to the right parts of the non-stationary versions of MHD Equations (2)–(5) those terms which are proportional to ∇ · ${\boldsymbol{B}}$. In the momentum equation, the term $[-{\boldsymbol{B}}({\rm{\nabla }}\cdot {\boldsymbol{B}})/4\pi ]$ has been added. In the energy equation, the term $[-({\boldsymbol{V}}\cdot {\boldsymbol{B}})({\rm{\nabla }}\cdot {\boldsymbol{B}})/4\pi ]$ has been added. In the magnetic field induction equation, we have added the term $[-{\boldsymbol{V}}({\rm{\nabla }}\cdot {\boldsymbol{B}})].$

In the calculations, we use a specific non-regular moving grid that allows us to perform exact fitting of the TS and heliopause. This grid allows us to decrease the sizes of the cells (e.g., increase the number of cells) in the vicinity of the discontinuities.

The example of a typical numerical grid in the (ZX) and (XY) planes is shown in Figure 2. The convergence of the numerical solution has been verified by changing the grid resolution up to two times in all directions and by changing cell sizes specifically in the vicinities of the discontinuities. Also, the numerical solution has been tested on the extended grids where the computational domain was increased up to two times. In particular, the tail region for model 3 varied in these calculations from 1000 AU up to 5500 AU.

Zoom In Zoom Out Reset image size

Kinetic Equation (1) for the H-atom component has been solved using a Monte Carlo method with spatial splitting of the particle trajectories (Malama 1991). Monte Carlo calculations have been performed using the same computational grid as described above for the plasma component. Since the TS and heliopause are the grid surfaces, we can define different populations of H atoms through their dependence on their place of origin (see, e.g., Malama 1991; Baranov et al. 1998). Optimization of the critical weights in the Monte Carlo method has been performed separately for each of the four populations (Malama 1990). This procedure increased the efficiency of the method.

First, we will explore the effects of the HMF. Figure 3 shows the results of numerical simulations for plasma in the ZX plane (see Section 2.3 for definition). Panels (A) and (C) present the results obtained in the frame of Model 1 (without HMF), while panels (C) and (D) represent Model 2 (with HMF). The panels on the top of the figure (panels (A) and (B)) show the plasma streamlines and the isolines of the plasma number density, while the panels at the bottom (panels (C) and (D)) show the magnetic field lines and isolines of HMF's magnitude.

Zoom In Zoom Out Reset image size

Visual comparison of the results shows that the HMF does not qualitatively change the global structure of the interaction region. The geometrical patterns and distributions of the plasma density and magnetic field look quite similar, especially outside of the heliopause.

However, it can clearly be seen that the heliospheric TS is about ∼10 AU closer to the Sun in Model 2 compared to Model 1. This effect is quantitatively seen in Figure 4, which shows one-dimensional (1D) distributions of the plasma and magnetic field in the direction of the Z-axis, i.e., toward the vector of the interstellar flow or in the upwind direction. In the upwind direction (and in the entire region close to the SW stagnation point that is shifted south of X = 0), the heliopause is located at the same distances for Models 1 and 2. However, the further one moves from the upwind region, the larger the difference between the HP distances in Models 1 and 2. In the crosswind (i.e., in the direction of positive X) the HP is about 100 AU further in Model 1 than in Model 2.

Zoom In Zoom Out Reset image size

This effect could be explained by a plasma depletion layer around the heliopause. The effect of plasma depletion in the vicinity of contact discontinuities is known for the magnetospheres (e.g., Zwan & Wolf 1976). The effect is connected with an increase in the perpendicular component of the magnetic field when plasma approaches the contact boundary (Cranfill–Axford effect for the heliosphere). The gradient of the magnetic field forces plasma to decelerate and flow around of the magnetic "wall" obstacle. This effect can clearly be seen in Figure 4 where the 1D distributions of the plasma and magnetic field parameters are shown for the upwind direction. All three components of the HMF increase toward the heliopause (plots (I), (J), and (K) in Figure 4). The plasma velocity components (plots (E), (F), and (G)) and density (plot (A)) decrease from the middle of the inner heliosheath to the heliopause. In the same region, both the Alfvenic number and plasma beta become less than 1, i.e., the flow is determined by the magnetic field. The effect of SW plasma density depletion around the heliopause can clearly be seen from the comparison of plasma density isolines in Models 1 and 2 (compare panels (A) and (B) in Figure 3).

Magnetic field pressure replaces the plasma pressure in such a way that the total pressure is conserved (plot (B) in Figure 4). This conclusion is valid for quite a wide region around the heliopause. However, because less plasma flows along the heliopause from the stagnation region to the flanks, some deficit of pressure appears at the flanks and the heliopause there moves toward the Sun. This explains why the heliopause in Model 2 has a slightly different (from Model 1) shape and is closer to the Sun in the crosswind and tail.

It is interesting to note that a sharp and strong jump of density at the heliopause exists in Model 2. This somewhat contradicts the Belov & Ruderman (2010) result for the smooth transition at the heliopause due to the charge-exchange effect. We explain this by the fact that the effect does not work in 3D geometry since the stagnation points are in different places inside and outside the heliopause.

We also note that we discuss plasma depletion from the inner side of the heliopause. A similar effect of plasma depletion should be pronounced from the outside of the heliopause, as discussed in Fuselier & Cairns (2013). We do not see the effect in the model because the interstellar H atoms strongly influence (due to charge exchange) the interstellar plasma flow around the heliopause. However, some small decrease of the plasma density that is observed in Figure 4(A) at the outer side of the heliopause could be associated with plasma depletion.

The effect of the magnetic increase (or magnetic wall) around the heliopause (plot (D) in Figure 3) and the plasma depletion leads to quite significant displacement of the heliospheric TS toward the Sun. This can be seen in Figure 4(B). In Model 2, the TS is ∼11 AU closer to the Sun compared with Model 1. Such a displacement of the TS is probably the major effect of the HMF on the global structure of the heliosphere.

Despite the fact that the geometrical pattern of the SW/LISM interaction region does not change dramatically, the HMF significantly influences the SW plasma flow. To illustrate this, we present the distribution of the plasma velocity and magnetic field along the heliopause (Figure 5). In this figure, the heliopause is projected onto a (projection) plane as follows. The origin of the coordinate system of this plane coincides with the intersection of the heliopause with the Z-axis. Polar coordinates (r, φ) in this plane are defined in such a way that the polar angle φ coincides with the polar angle in the (XY) plane. The radial distance r coincides with the length of the curve that is obtained at the intersection of the heliopause and the plane containing the Z-axis and determined by constant φ angle in the (XY) plane.

Zoom In Zoom Out Reset image size

Note that the (BV) plane (i.e., ZX plane) corresponds to the horizontal line Y = 0 in the projection plane. For Model 1, the plasma flow pattern is symmetric around this line because the (ZX) plane is the plane of symmetry in this model. This is not the case for Model 2. The main symmetry is observed around line φ = 424. This line is defined by the angle between the solar equatorial plane and the (BV) plane. Therefore, the SW plasma flow pattern in the heliosheath (and, specifically, at the heliopause) is determined by the direction of the solar rotation axis rather than by the directions of the interstellar flow and magnetic field vectors. It can be seen (panel (B) in Figure 5) that the SW velocity is smaller in the equatorial plane and the streamlines directed toward the poles. Some asymmetries with respect to the plane of φ = 424 are observed in the figure. This is due to the asymmetric shape of the heliopause distorted by the IsMF. However, these "interstellar" induced features are much less pronounced than the "Sun-connected" symmetry of the flow.

The same types of symmetries/asymmetries are also observed in the distribution of the HMF along the heliopause (panel (C) of Figure 5).

As opposed to the inner side of the heliopause, in the interstellar medium, the plasma and magnetic field distributions look very similar in Models 1 and 2 (panels (D)–(F) in Figure 5). Small differences in panels (D) and (E) are caused by differences in the heliopause shape in Models 1 and 2.

Before we move forward to discuss the results obtained in the context of Model 3, it is worthwhile to mention the very recent publication by Opher et al. (2015), where it was pointed out that the HMF tension-force in the inner heliosheath leads to the formation of two SW flux jets directed north and south. Furthermore, this tension force changes the topology of the heliopause, making it a more tube-like structure (see, also, axisymmetric model by Drake et al. 2015) as compared to our sheet-like structure. The model by Opher et al. (2015) assumes a unipolar magnetic field, as in our Model 2. Therefore, the effects of magnetic tension in the heliosheath can be seen in our Model 2. In order to explore these effects, we plot (Figure 6) the solar flux (ρV) in projection on a closed surface located in the inner heliosheath. The surface has been constructed by adding a fixed distance to the TS distance in any given direction. The distance has been chosen as 39 AU for Model 1 and 50 AU for Model 2. The difference in the distances chosen for Models 1 and 2 is connected with the different widths of the inner heliosheath in the models. Figure 6 presents the solar flux plotted for Model 1 and Model 2. Any point on the surface is characterized by two spherical angles θ and ϕ. The angle θ is counted from the Z-axis (i.e., from the upwind direction), while the angle ϕ is counted from the X-axis in the XY plane.

Zoom In Zoom Out Reset image size

Panel (1) in Figure 6 demonstrates the SW mass flux in the case of Model 1 when HMF is not taken into account. It is seen that the mass flux reaches its maximum for θ ∼ 100°–110° and ϕ ∼ 180^o. This nearly crosswind direction corresponds to the direction through which all of the solar mass flux from the upwind hemisphere passes into the tail. The maxima at ϕ ∼ 180^o corresponds to that part of the (BV) plane where the distance between the TS and HP is minimal. Conversely, panel (2) demonstrates that there are two mass-flux maxima (i.e., jets) in the north and south directions. These jets are clear results of the magnetic tension. Therefore, in this respect, our model results in qualitative agreement with the conclusion of Opher et al. (2015).

However, the conclusion on the tube-like shape of the heliopause made by Opher et al. (2015) is not supported by our calculations. The exact reasons for this will be explored in future studies. Here, we speculate on two possibilities. The first is connected with effects of charge exchange of solar protons and interstellar H atoms. In the axisymmetric model of Drake et al. (2015), where the spherically symmetric SW flows into the unmagnetized LISM resting with respect to the Sun, the axis of symmetry is the north–south pole axis. Due to the symmetry in any plane containing the axis of symmetry, there is a stagnation point at equator. Therefore, in 3D space, the stagnation line (i.e., the line of stagnation points) at the heliopause is a circle in the solar-equatorial plane. Since the tube-like structure remains in the case of Opher et al. (2015) when the relative Sun/LISM motion is taken into account, the stagnation line still exists, although it is deformed and not necessarily located in the equatorial plane. The interstellar H atoms penetrate into the stagnation line vicinity and begin charge exchange with the protons. As a result of the charge exchange, momentum is transferred from the H atom components into the plasma components. The momentum transfer results in a diminishing heliocentric distance to the heliopause in the upwind side of the heliosphere. The effect is very well known for the traditional sheet-like heliosphere (see, e.g., Figure 2 in Izmodenov 2000). In the tube-like heliosphere, the same momentum-transfer effect should result in the increase of the heliocentric distance to the heliopause in the downwind direction. The heliopause, in principle, could move further and further downwind, and so it can be treated as a sheet-like structure in the considered computational domain. However, while the effect of charge exchange is, in principle, taken into account in the multi-fluid approach, this effect should be verified by kinetic-MHD calculations. As has been shown by Alexashov & Izmodenov (2005), the multi-fluid and kinetic approaches may have qualitatively different results.

Another possible reason Model 2 does not have the tube-like shape of the heliopause could be due to the fact that our computational grid explicitly assumes the sheet-like structure of the heliopause. However, if the formation of the two streams would physically result in the tube-like shape geometry of the heliopause, then at some specific locations (where the SW mass flux becomes negligible) the different parts of our "sheet-like" heliopause should approach each other. Then, our computational grid should collapse and this would be an indication of the tube-like topology of the heliopause. However, grid collapse is not observed in our calculations. Therefore, we currently favor the first possibility, but the question of the topology of the heliopause is generally very interesting and should be explored in detail in the future.

Figure 7 presents the results obtained in the framework of Model 3 where the latitudinal variations of the SW parameters have been taken into account in addition to the HMF. It can be seen that both the global structure of the heliosphere and the plasma and magnetic field distributions are qualitatively the same as for Models 1 and 2. However, the comparison of panel (A) in Figure 7 with panel (B) in Figure 3 shows that the heliopause is more blunt in Model 3. The heliopause is closer (by about ∼10–15 AU) to the Sun in the upwind part of the heliosphere and further from the Sun in the crosswind. We associate this with the non-monotonic heliolatitudinal distribution of the SW dynamic pressure, with 30%–50% increase at the middle latitudes and further decrease toward the poles (see Figure 1). Such behavior of the dynamic pressure changes the plasma flow beyond the TS in such a way that more plasma evacuates to the flanks which makes the pressure smaller in the upwind part and larger further from the upwind part.

Zoom In Zoom Out Reset image size

Following the heliopause shape, the TS is also more elongated toward the solar poles. Since it is quite difficult to see this effect in the large scale of Figure 7(A), we present a table of the TS and heliopause positions in the upwind direction, and in the Voyager 1 and Voyager 2 directions. It can be seen that in the upwind direction, the TS is only slightly (∼1 AU) further away in Model 3 compared with Model 2. The effect is much larger in the directions of Voyager 1 and 2. The differences are 9 and 6 AU, respectively.

Panels (B), (C), (E), and (F) in Figure 7 show the projections of plasma velocities and magnetic fields along the heliopause. Panels (B) and (E) show the distribution inside the heliopause, panels (C) and (F) show it outside the heliopause. All of the main features are essentially the same as in Model 2. Inside the heliosphere, the flow is oriented in accordance with the solar equatorial plane. Although north–south asymmetry can clearly be seen in Figure 7(B). From the interstellar side of the heliopause, the flow and the magnetic field are oriented along the Y = 0 axis which is the projection of the (BV) plane. Some asymmetries are due to non-symmetric shape of the heliopause.

Since the shape of the heliopause has been changed in Model 3 and the heliopause can be considered as an obstacle to the interstellar flow, the interstellar plasma streamlines are somewhat different in Models 2 and 3. This is seen from the comparison of Figures 7(A) and 3(B). The interstellar atoms are coupled with the interstellar plasma by charge exchange. Therefore, the change in the plasma streamlines could influence the distribution of interstellar H atoms inside the heliosphere. Despite the fact that these effects are quite important, they are not visible in a global-scale picture. Qualitatively, the distributions of the interstellar H atoms in the (BV) plane for Models 2 and 3 look very similar to Model 1. The distributions can be found in Izmodenov & Baranov (2006, Figure 4.7(B)) or Izmodenov (2006, Figure 5(E)). To explore the difference in H atom distribution quantitatively, we present (Figure 8) the zero, first, and second moments of the velocity distribution function of H atoms at 90 AU as functions of the angle θ that is counted in the ZX plane from the positive direction of Z (i.e., from the upwind direction). It is worthwhile to note here that it is convenient to separate the interstellar atoms by primary and secondary populations. The primary population is the original interstellar atoms that penetrate into the heliosphere without charge exchange with protons. The secondary population is the atoms created in the disturbed LISM around the heliopause. Historically, in the case of the supersonic LISM, the atom is defined as a secondary atom if it is created by charge exchange in the region between the BS and the heliopause. Such a definition does not work for models with no BS, such as the models considered in this paper. In this case, we identify a newly originated H atom as the secondary atom if the interstellar plasma is sufficiently disturbed at the place of the atom's origin. More precisely, under sufficiently disturbed conditions, we assume that the temperature of the plasma is more than ∼10000 K.

Zoom In Zoom Out Reset image size

Although the total (primary + secondary) number density of the interstellar atoms remains the same for all of the considered models, the proportions of the primary and secondary populations are slightly different (Figures 8(A) and (G)). The changes in the interstellar plasma flow around the heliopause (due to different shape of the heliopause) lead to small (less than 1 km s⁻¹) changes in the velocities of the primary and secondary components (Figures 8(B) and (H)). The kinetic temperature components (which are defined as the second moments of the velocity distribution) for the secondary H atoms are changed within 1000–2000 K. This is a very noticeable effect. Another effect is in the deflection of the direction of the interstellar H-atom flow inside the heliosphere from the direction of the interstellar gas flow. The deflection can be measured by the deflection angle ${\alpha }_{\mathrm{defl}}=\mathrm{arccos}({V}_{z,{\rm{H}}}/{V}_{{\rm{H}}}),$ where V_z,H is the z component of the H atom velocity and V_H is the velocity. The deflection is associated with the change of the interstellar plasma streamlines around the heliopause due to the IsMF (Izmodenov et al. 2005a). Figures 7(C) and (I) show that the deflection angle between the models could be different by of the order of 1°, which is small but might be important for analyses of backscattered Lyα (see, e.g., Katushkina et al. 2015).

Before summarizing the results described in this paper, we have to give the reasons for our choice of the proton and hydrogen number densities in the LISM because this was not done in Section 2.1. Although we did not aim in this paper to discuss the relevance of the model to all available observational diagnostics of the SW/LISM interaction region, we would like to note that the presented model has only four free parameters—interstellar proton and hydrogen number densities, magnitude of the IsMF, and the direction of the magnetic field in the (BV) plane. Other model parameters are not free since they are based on available observational information discussed in detail in Sections 2.1 and 2.2.

In order to constrain the remaining free parameters, the following observational diagnostics can be used.

• 1.  
  The distances of the heliospheric TS in the directions of the Voyager 1 and Voyager 2 spacecraft. They were determined by the spacecraft crossings as 83.7 and 94.1 AU in 2004 August and in 2007 December, respectively.
• 2.  
  The number density of interstellar H atoms in the distant heliosphere (at 50–90 AU) of the order of 0.10 ± 0.01 cm⁻³. This value has been determined in different analyses of pick-up proton data (see, e.g., Geiss et al. 2006; Bzowski et al. 2008) and the SW deceleration at large heliospheric distances due to charge exchange with protons (e.g., Richardson et al. 2008).
• 3.  
  The difference of the interstellar H atom flow direction inside the heliosphere compared with the direction of the pristine circumsolar LISM (Lallement et al. 2005, 2010).
• 4.  
  The SW plasma density and velocity data in the heliosheath from the Plasma/Voyager 2 instrument.
• 5.  
  The distance to the heliopause in the direction of Voyager 1 of 122 AU in 2012 August–September.
• 6.  
  The magnetic field strength and direction observed by MAG/Voyager 1 after 2012 September (e.g., Burlaga & Ness 2012).
• 7.  
  The proton number density inferred from analyses of the plasma wave instrument by Gurnett et al. (2013).
• 8.  
  It is widely believed in the community (e.g., Funsten et al. 2013; Schwadron & McComas 2013; Heerikhuisen et al. 2014) that the center of the observed ribbon seen by IBEX is closely related to the direction of the IsMF and can serve as an additional constraint on the IsMF.

While this list of observational diagnostics is not complete (for example, energetic particle data should be added as well), it is quite evident that in the context of the current model paradigm it is impossible to find the model parameters to reconcile all of the available diagnostics.

For example, it is impossible to get the TS at 94 AU and the HP at 122 AU in Voyager 1 in the frame of the global model with the same realistic boundary conditions. Solar cycle variations of the SW parameters can not resolve this problem as is clearly seen from time-dependent simulations (see, e.g., Izmodenov et al. 2005b, 2008), which show that the TS fluctuates with the solar cycle within ±10 AU, while fluctuations of the HP are ±3–4 AU around its mean value. In order to obtain agreement between the observed distances to the TS and HP, Izmodenov et al. (2014) suggested adding a physical effect such as thermal conduction, which would decrease the thermal pressure in the heliosheath. Another example of diagnostic inconsistency is the direction of the IsMF. The direction of the IBEX ribbon center does not belong to the (BV) plane derived from backscatterred Lyα analyses by Lallement et al. (2010). The detailed analyses of the inconsistencies is outside the scope of this paper and should be a subject of further study.

Since the model described in this paper has been used to study interstellar H atoms (e.g., Katushkina et al. 2015, this issue), we choose to fulfill the first three of the diagnostics listed above. The first is the distances to the TS, which provide the correct size of the heliosphere. The second and third are directly related to the interstellar atoms, and it makes sense to use this parameters for analyses of the interstellar H atom fluxes, backscattered solar Lyα, and heliospheric absorption in Lyα. Since complete parametric studies are very numerically expensive for the considered 3D kinetic-MHD model, below we describe how we obtained the interstellar proton and H atom densities used in this paper.

The interstellar proton and H-atom number densities of ∼0.06 and ∼0.18 cm⁻³ were established in the parametric study by Izmodenov et al. (2003). The axis-symmetric Baranov–Malama model was used in that study (see, also, Figure 1 in Izmodenov 2009). Izmodenov (2009) performed a new parametric study in the context of a 3D kinetic-MHD model (Model 1) and concluded that the model with the following set of model parameters gives results agreeable with the observational diagnostics listed above: n_p,LISM = 0.05 ± 0.015 cm⁻³, n_H,LISM = 0.18 ± 0.018 cm⁻³, B_LISM = 2.5–3.5 μG, α_LISM = 15°–30°. Shortly afterward, Izmodenov (2009) reported another possible solution with stronger IsMF: B_LISM = 4.4 μG, α_LISM = 20°.

In the present study, we started with interstellar parameters obtained by Izmodenov (2009). However, due to the presence of the HMF in the model, the TS is ∼10 AU closer to the Sun. In order to obtain a TS distance in the model comparable with Voyager's, we have to reduce the interstellar pressure. The reduction of the IsMF magnitude is not possible because a smaller IsMF would lead to (1) less distortion of the heliopause and, therefore, smaller difference in the distances of the TS in the direction of Voyager 1 and 2; and (2) smaller deflection angle α_defl. Therefore, only the interstellar plasma pressure could be reduced. LISM velocity is not a free parameter because it is known from the Ulysses/Gas and IBEX observations. Therefore, to decrease the interstellar pressure, we have to decrease the interstellar proton number density n_p,LISM. The decrease of n_p,LISM would result in smaller filtration of H atoms in the SW/LISM interaction region and, therefore, in a larger amount of neutrals inside the heliosphere. The latter would contradict the second of the observational constraints listed above. Therefore, n_H,LISM should be reduced as well. In this manner, we ended up with the values of n_p,LISM = 0.04 cm⁻³, n_H,LISM = 0.14 cm⁻³, B_LISM = 4.4 μG, α_LISM = 20°. With these parameters, our Model 3 satisfies (with some uncertainties) all of the observational constraints listed above. For the first and second diagnostics, see Table 1 and Figure 7, respectively. The thrid diagnostic—the deflection angle—was intensively studied in Katushkina et al. (2015). It is important to note that we did not aim to have exact distances in the Voyager 1 and Voyager 2 directions, since they would change with the solar cycle variations of the SW parameters (see, Izmodenov et al. 2008).

In this paper, we presented the results of our latest kinetic-(ideal)MHD model. The novel parts of the model are the HMF and the heliolatitudinal variations of the SW density and velocity. The main results reported in this paper are as follow.

1. The HMF results in the following:

• i.  
  closer heliocentric distance to the TS by ∼10 AU compared to models without the field;
• ii.  
  change of the heliopause shape when the heliocentric distances to the heliopause remain the same in the upwind part of the heliosphere and become smaller in the crosswind and tail parts;
• iii.  
  increase of the HMF in the heliosheath toward the heliopause leads to the formation of the plasma depletion layer around the heliopause, similar to those observed in the magnetospheres of planets;
• iv.  
  the SW plasma flow in the inner heliosheath (and, in particular, at the heliopause) has general symmetry with respect to the solar equatorial plane, while some small asymmetries connected with the distortion of the heliopause shape by the IsMF are observed;
• v.  
  conversely, the interstellar plasma flow around the heliopause is oriented by the interstellar velocity and magnetic field directions, although asymmetries due to the three dimensions of the heliopause exist as well.

2. The latitudinal variations of the SW dynamic pressure adopted in the model have maxima at middle heliolatitudes. This results in a more blunt heliopause when the heliopause becomes closer to the Sun in its upwind part and further away from the Sun on the flanks. The general orientations of the plasma flows inside and outside the heliopause remain the same—the SW flow is oriented around the solar equatorial plane, and the interstellar plasma flow is oriented around the plane determined by the interstellar velocity and magnetic field vectors. We note here that in Model 3 there are no strong discontinuities in the magnetic field components through the heliopause in the Voyager 1 direction. However, the detailed comparison of the model results with Voyager magnetic field data is not the subject of this paper and will be explored separately.

3. We have found that the model with the following parameters n_p,LISM = 0.04 cm⁻³, n_H,LISM = 0.14 cm⁻³, B_LISM = 4.4 μG, α_LISM = 20° provides reasonable agreement with the distances of the heliospheric TS in the directions of the Voyager 1 and Voyager 2 spacecraft, the number density of interstellar H atoms in the distant heliosphere (at 50–90 AU) of the order of 0.10 ± 0.01 cm⁻³, and the difference of the interstellar H atom flow direction inside the heliosphere as compared with the direction of the pristine circumsolar LISM.

The model presented in this paper is not directly relevant to analyses of IBEX-Hi data because it does not consider the pick-up protons as a separate component, as was done, for example, by Malama et al. (2006) for the axis-symmetric case or by Chalov et al. (2010) for the 3D case. Such a model is in progress and will be reported in the near future. At the same time, the model presented here is quite appropriate for analyzing IBEX-Lo hydrogen data and has been used in the study presented in a companion paper by Katushkina et al. (2015).

This paper is devoted to the memory of our senior colleague and Teacher—Y.G. Malama. We thank Olga Katushkina for her help with the preparation of this paper. Numerical calculations were performed using the supercomputers "Lomonosov" and "Chebyshev" of the Supercomputing Center of Lomonosov Moscow State University. The heliospheric part of the paper has been supported by RFBR grant 14-02-00746. The numerical modeling of the global heliosphere/astrosphere was performed under Russian Science Foundation grant 14-12-01096.