Heliospheric Energetic Neutral Atoms: Non-stationary Modeling and Comparison with IBEX-Hi data
MMNRAS , 1–15 (2020) Preprint 30 September 2020 Compiled using MNRAS L A TEX style file v3.0
Heliospheric Energetic Neutral Atoms: Non-stationary Modelingand Comparison with
IBEX-Hi data
I. I. Baliukin, , , (cid:63) V. V. Izmodenov, , , and D. B. Alexashov Space Research Institute of Russian Academy of Sciences, Profsoyuznaya Str. 84/32, Moscow, 117335, Russia Moscow Center for Fundamental and Applied Mathematics, Lomonosov Moscow State University, GSP-1, Leninskie Gory, Moscow, 119991, Russia National Research University Higher School of Economics, Moscow, Russia Institute for Problems in Mechanics, Vernadskogo 101-1, Moscow, 119526, Russia
Accepted XXX. Received YYY; in original form ZZZ
ABSTRACT
The Interstellar Boundary Explorer (
IBEX ) has been measuring fluxes of the EnergeticNeutral Atoms (ENAs) using the
IBEX-Hi (0.3 – 6 keV) instrument since 2008. We havedeveloped a numerical time-depended code to calculate globally distributed flux (GDF) ofhydrogen ENAs employing both 1) 3D kinetic-MHD model of the global heliosphere and 2)reconstruction of atom trajectories from 1 au, where they are observed by
IBEX , to the pointof their origin in the inner heliosheath (IHS). The key factor in the simulation is a detailedkinetic consideration of the pickup ions (PUIs), the supra-thermal component of protons inthe heliosphere, which is “parental" to the ENAs and originates in the region of the supersonicsolar wind being picked by the heliospheric magnetic field. As a result of our study, we haveconcluded that (1) the developed model is able to reproduce the geometry of the multi-lobestructure seen in the
IBEX-Hi
GDF maps, (2) the GDF is extremely sensitive to the form of thevelocity distribution function of PUIs in the IHS, and the accounting for the existence of anadditional energetic population of PUIs is essential to explain the data, (3) despite a relativelygood agreement, there are some quantitative differences between the model calculations and
IBEX-Hi data. Possible reasons for these differences are discussed.
Key words:
ISM: atoms — ISM: magnetic fields — Sun: heliosphere
The solar system is surrounded by the local interstellar medium(LISM) and moving through it with a bulk velocity of ∼
26 km s − (e.g., Witte 2004; McComas et al. 2015). The interaction of thesolar wind (SW) with the ionized component of the LISM formsa complex structure that is called the heliospheric interface. Theheliopause (HP) is a tangential discontinuity that separates the SWand interstellar plasmas from each other. There are two shocks inthe heliospheric interface – the termination shock (TS), where theSW is slowed down from supersonic to subsonic speed, and a bowshock (BS), where the interstellar flow is slowed down, but theexistence of latter is under discussion (see, e.g., Izmodenov et al.2009; McComas et al. 2012; Zank et al. 2013). The region of thecompressed and heated plasma between the shocks is commonlycalled a heliosheath and the region between the TS and HP – theinner heliosheath (IHS).The neutral component of the LISM consists mainly of hydro-gen atoms. The interstellar H atoms have a large mean free pathfor charge exchange, which is comparable with the characteristic (cid:63) E-mail: [email protected] size of the heliosphere (Izmodenov 2001), and due to the relativemotion of the Sun and LISM, they can penetrate the heliosphere. Inthe heliosheath, H atoms may experience charge exchange with hotprotons, which results in the production of energetic neutral atoms(ENAs). Some of these ENAs propagate close to the Sun and Earth’sorbit, where they can be measured.The
Interstellar Boundary Explorer ( IBEX ) spacecraft waslaunched into a highly elliptical orbit around Earth in October 2008to obtain the first all-sky maps of the neutral gas/plasma interactionsat the heliospheric boundary and to directly sample the interstellargas flow through the inner heliosphere (McComas et al. 2009). Toachieve these goals
IBEX is a Sun-pointed spinning satellite, whichcarries two energetic atom sensors:
IBEX-Lo (0.01 – 2 keV) and
IBEX-Hi (0.3 – 6 keV). A detailed description of the
IBEX-Hi sen-sor may be found in Funsten et al. (2009). The
IBEX-Hi instrumentis measuring ENA fluxes and these data are one of the few sourcesof knowledge about the structure of the heliospheric boundary, im-posing significant limitations on the parameters of the heliosphericmodels. McComas et al. (2020) have examined IBEXâĂŹs globalENA observations over a full solar activity cycle (Solar Cycle 24),covering 11 years from 2009 through 2019.The observations of ENAs by
IBEX-Hi have revealed two pop- © 2020 The Authors a r X i v : . [ phy s i c s . s p ace - ph ] S e p I. I. Baliukin et al. ulations: one of them is emitted from a narrow circular part of thesky that is called “ribbon", and a globally distributed flux (GDF)that is controlled by processes in the heliosheath (McComas et al.2009). Schwadron et al. (2011, 2014) have performed the analy-sis of
IBEX-Hi data and developed the technique to separate ENAemissions in the ribbon from the GDF, which, as it turned out, hasa complex multi-lobe structure. It is commonly assumed that theGDF originates from the IHS, where the supersonic SW and pickupions (PUIs) are slowed and heated after crossing the TS. The PUIsare formed in a result of ionization of interstellar hydrogen atomsmainly through charge exchange in the heliosphere, where they arepicked by the heliospheric magnetic field. The ribbon, in turn, isformed by secondary charge exchange in the outer heliosheath (see,e.g., McComas et al. 2009; Chalov 2010; Heerikhuisen et al.2010).There are two populations of protons in the heliosphere – thecold thermal population of core SW protons and the hot supra-thermal population of PUIs. The PUIs in the inner heliosheath canbe divided into subpopulations, transmitted or reflected, dependingon its interaction with TS (Zank et al. 2010). In some works onthe modeling of ENA fluxes the authors made attempts to modelthe quite distinct populations âĂŞ thermal and pickup protons âĂŞusing one kappa-distribution (e.g., Heerikhuisen et al. 2008) or asuperposition of Maxwell distributions (Zank et al. 2010; Zirnsteinet al. 2017; Kornbleuth et al. 2018; Shrestha et al. 2020).In this paper, we investigate the GDF maps produced in theframe of the latest heliospheric model of the Moscow group (Iz-modenov & Alexashov 2015, 2020), and perform its comparisonwith
IBEX-Hi data (Schwadron et al. 2014). The key factor inthe simulation is a detailed kinetic consideration of the PUIs, thesupra-thermal component of protons in the heliosphere. The paper isorganized as follows. Section 2 describes the method of calculationof the ENA fluxes. In Section 3, the detailed description of the PUIdistribution model is provided. Section 4 describes the techniqueof partitioning charged particles into components from single-fluidplasma calculations. In Section 5, the results of the numerical cal-culations and their comparison with
IBEX-Hi data are presented. InSection 6, the qualitative effect of additional energetic PUI popula-tion on the ENA fluxes is discussed. Finally, Section 7 provides anoverall summary of our work.
The primary ENAs, which are the source of GDF, are born incharge exchange between the SW protons and H atoms in the innerheliosheath (IHS). The directional differential ENA flux (in the solarinertial reference frame) is a line of sight (LOS) integral: j ENA ( t obs , r obs , v , LOS ) = m H ∫ s HP s TS ν H ( t , r , v ) f ( t , r , v ) v S p , ENA ds , (1)where t obs is the moment of observation, r obs is the position ofthe observer, LOS is the unit vector in the line of sight direction, v = v · LOS is the velocity of an ENA, m H is the mass of H atom, f = f sw + f pui is the velocity distribution of protons in the IHS(the sum of core SW proton and PUI distribution functions). Theintegration is performed along the ENA trajectory in the IHS (in theregion between the TS and HP), with ds = v dt being the differentialpath length. The values of the variable s along the trajectory s TS and s HP correspond to atom intersections of the TS and HP, respectively.In principle, the velocity of an ENA is changing along the trajectory due to the influence of gravitational and radiation pressure forcesby the Sun, i.e. v = v ( s ) , but for the energies under consideration(in the IBEX-Hi energy range) the change of velocity is negligible,so we assume that the velocity is constant along the trajectory.In Equation (1) ν H is the production rate of ENAs due to thecharge exchange of protons with H atoms, which is defined as ν H ( t , r , v ) = ∫ ∫ ∫ f H ( t , r , v H )| v − v H | σ ex (| v − v H |) d v H , (2)where f H ( t , r , v H ) is the H velocity distribution function and σ ex is the effective charge exchange cross section that depends on therelative atom-proton velocity. In our calculations, the cross sectionfrom Lindsay & Stebbings (2005) was used. The extinction ofENAs is determined by survival probability S p , ENA = exp (cid:18) − ∫ t obs t ν ion ( τ, r ( τ ) , v ( τ )) d τ (cid:19) , where ν ion is the total ionization rate due to the ionization pro-cesses (charge exchange with protons, photoionization, and elec-tron impact), i.e. ν ion = ν ex + ν ph + ν imp . In our calculations weneglect electron impact ionization ( ν imp = ν ph = ν ph , E ( t , λ ) · ( r E / r ) , where λ is heliolatitude and r E = 1 au.The temporal and heliolatitudinal variations of the photoionizationrate ν ph , E ( t , λ ) at 1 au adopted in our model were obtained fromdifferent experimental data ( OMNI , SWAN/SOHO ) the same way asit was performed in Katushkina et al. (2015). For the stationarymodel calculations we assume constant photoionization rate at 1 au ν ph , E = . × − s − as it was taken in Izmodenov & Alexashov(2015). The charge exchange ionization rate ν ex is calculated as ν ex ( t , r , v ) = ∫ ∫ ∫ f ( t , r , v p )| v − v p | σ ex (| v − v p |) d v p , where v p is the velocity of proton.To calculate the differential fluxes measured by IBEX-Hi atdifferent energy channels, the energy transmission of
IBEX-Hi elec-trostatic analyzers were taken into account (for details, see AppendixA). Accounting for energy transmission leads to the re-distributionof fluxes between the energy channels and the spreading of the ob-served spectrum. Important to note, that it is possible to make somesimplifications in order to optimize the calculations. For energiesunder consideration the production rate of ENAs ν H can be safelyapproximated in Equation 2 as ν H ( t , r , v ) ≈ n H ( t , r ) v σ ex ( v ) , since | v − v H | ≈ v , so we use this simplification in our simulations.Thus, to calculate the fluxes of ENAs, we need to know thevelocity distribution function of protons in the heliosphere. The fol-lowing section will describe the method to calculate the distributionof PUIs using the plasma and interstellar H atom distributions inthe heliosphere obtained in the frame of the global kinetic-MHDmodel by Izmodenov & Alexashov (2020). The distribution function of PUIs is anisotropic in the case of weakscattering, i.e. when the SW turbulence level is low (the quiet solarwind). However, as follows from theoretical estimates and obser-vations (Gloeckler et al. 1994), the distribution function is almostisotropic in the disturbed SW, so the process of an effective pitch-angle diffusion must be operative. Therefore, we assume that thevelocity distribution of pickup protons in the SW rest frame is
MNRAS , 1–15 (2020) eliospheric ENAs: modeling and comp. with IBEX-Hi isotropic, and it is determined by the velocity distribution function f pui ( t , r , v ) in the heliocentric coordinate system by the expression: f ∗ pui ( t , r , w ) = π ∫ ∫ f pui ( t , r , v ) sin θ d θ d ϕ, where v = V ( r , t ) + w , v and V are velocity of pickup proton andbulk velocity of the plasma in the heliocentric coordinate system, w is the velocity of the pickup proton in the SW rest frame, and( w , θ , ϕ ) are coordinates of w in the spherical coordinate system.The kinetic equation for f ∗ pui ( t , r , w ) can be written in the followinggeneral form (see, e.g. Isenberg 1987; Chalov et al. 2003): ∂ f ∗ pui ∂ t + V · ∂ f ∗ pui ∂ r = w ∂∂ w (cid:32) w D ∂ f ∗ pui ∂ w (cid:33) + w ∂ f ∗ pui ∂ w div ( V ) + S ( t , r , w ) , (3)taking into account velocity diffusion (where D ( t , r , w ) is the ve-locity diffusion coefficient) but ignoring spatial diffusion. The es-timations of the spatial diffusion coefficient (e.g., Scherer et al.1998, its table 4) show that for the energies under considera-tion the spatial diffusion can be neglected (see, also, Rucinski etal. 1993; Chalov & Fahr 1997). The source term S ( t , r , w ) = S + ( t , r , w ) − f ∗ pui ( t , r , w ) S − ( t , r , w ) , where S + and S − are responsiblefor production and losses (extinction) of PUIs, respectively, and canbe calculated as S + ( t , r , w ) = π ∫ ∫ f H ( t , r , V + w ) ν ion ( t , r , V + w ) sin θ d θ d ϕ, (4) S − ( t , r , w ) = π ∫ ∫ ν H ( t , r , V + w ) sin θ d θ d ϕ. (5)In the IHS, neutrals can interact both with SW protons andPUIs through the charge exchange process, so ν ex = ν ex , sw + ν ex , pui and ν ex , sw ( t , r , V + w ) = ∫ ∫ ∫ f sw ( t , r , V + w sw )| w − w sw |·· σ ex (| w − w sw |) d w sw , (6) ν ex , pui ( t , r , w ) = ∫ f ∗ pui ( t , r , w pui )·· (cid:18)∫ ∫ | w − w pui | σ ex (| w − w pui |) w sin θ d θ d ϕ (cid:19) d w pui . (7) In this paper, we consider a simple model and adopt D = d r dt = V , (8)which for the stationary case is also a streamline, and PUI velocitychanges along it according to d w dt = − w ( V ) . (9) Equation (3) has the following solution f ∗ pui ( t , r ( t ) , w ( t )) = ∫ tt S + ( τ, r ( τ ) , w ( τ )) S p , pui ( τ, t ) d τ ++ f ∗ pui ( t , r ( t ) , w ( t )) S p , pui ( t , t ) , (10)where S p , pui ( τ, t ) describes the loss of PUIs due to neutralization ontheir way from point ( τ, r ( τ ) , w ( τ )) to ( t , r ( t ) , w ( t )) : S p , pui ( τ, t ) = exp (cid:18) − ∫ t τ S − ( ˆ τ, r ( ˆ τ ) , w ( ˆ τ )) d ˆ τ (cid:19) . (11)Using Equations (8) and (9) the trajectory of PUI is reconstructedbackward in time from the point of phase space ( t , r ( t ) , w ( t )) topoint ( t , r ( t ) , w ( t )) , where the characteristic is close to the Sun,and it can be safely assumed that f ∗ pui ( t , r ( t ) , w ( t )) =
0. In ourcalculations we use the following inner boundary: r ( t ) = R = . In the frame of the kinetic theory, the moments of the velocitydistribution function of PUIs at point ( t , r ) ∈ R are the followingvalues: • zero velocity distribution function moment – number density: n pui ( t , r ) = π ∫ f ∗ pui ( t , r , w ) w d w ; (12) • second velocity distribution function moment – kinetic tem-perature: T pui ( t , r ) = π m p n pui k B ∫ f ∗ pui ( t , r , w ) w d w , (13)where m p is the mass of proton, and k B is the Boltzmann constant.The pressure of PUIs can be calculated as p pui = n pui k B T pui . The usage of the Liouville’s theorem (phase space flow conservationover the shock), the conservation of the magnetic moment (firstadiabatic invariant), and assumption of the weak scattering leads tothe following jump condition at the shock (Fahr & Siewert 2011,2013): f ∗ pui , d ( t , r , w ) = sC / f ∗ pui , u (cid:18) t , r , w √ C (cid:19) , (14)where C ( s , ψ ) = ( A ( s , ψ ) + B ( s , ψ ))/ , A ( s , ψ ) = (cid:113) cos ψ + s sin ψ, B ( s , ψ ) = s / A . In Equation (14) f ∗ pui , u and f ∗ pui , d are the values of PUI dis-tribution function upstream and downstream the TS, ψ ( t , r ) and s ( t , r ) = n d / n u are the local upstream shock-normal angle (betweenthe magnetic field and normal to the shock surface) and the shockcompression factor that depend on the position r and moment t ofthe TS crossing. From the observations by Voyager 1 the compres-sion ratio is 2.4 (for TS-2 crossing, see Richardson et al. 2008a),and in the global model simulations by Izmodenov & Alexashov
MNRAS000
MNRAS000 , 1–15 (2020)
I. I. Baliukin et al. (2020) it is ≈ n pui , d / n pui , u = s , T pui , d / T pui , u = C , p pui , d / p pui , u = sC . (15)Important to note that some PUIs can be reflected at the TS, sothe reflection process leads to anisotropy (in the SW rest frame) ofthe velocity distribution of PUIs in some vicinity of the TS (Chalovet al. 2015). In principle, the condition (14) can be modified byintroducing at the TS the generation of distinct populations of PUIs(transmitted and reflected). To carry out calculations of the PUI distribution as described in theprevious sections, the global distributions of H atoms and plasmashould be known. We have performed global heliospheric simu-lations of SW/LISM interaction using kinetic-MHD model by Iz-modenov & Alexashov (2020) in the stationary and time-dependentcases. Hereafter we will refer to this model as
IA2020 . The inter-stellar parameters of the models are the following: the bulk velocityand temperature are V LISM = 26.4 km/s and T LISM = 6530 K; thedirection of V LISM is (longitude = 75.4 ◦ , latitude = -5.2 ◦ ) in eclip-tic (J2000) coordinate system; the H atom, proton, and helium ionnumber densities are n H , LISM = 0.14 cm − , n p , LISM = 0.04 cm − ,and n He + , LISM = 0.003 cm − , respectively. The SW parameters atEarthâĂŹs Orbit are the following: Mach number is 6.44 (corre-sponds to SW temperature T E = 188500 K); the number densityof the alpha particles He ++ is 3.5% of the proton number density.The heliolatitude and time variations of the SW were obtained fromdifferent experimental data ( OMNI 2 dataset, interplanetary scintil-lation data,
SWAN/SOHO full-sky Lyman- α maps), see AppendixA in Izmodenov & Alexashov (2020) for details. The distribu-tion of the solar wind proton number density and velocity at 1 AUas functions of time and heliolatitude are shown in Izmodenov &Alexashov (2020, Fig. A.1). For the heliospheric magnetic field,the Parker spiral solution has been assumed at 1 au with magneticfield magnitude B E = 37.5 µ G at 1 au. The configuration of the in-terstellar magnetic field is chosen as B LISM = 3.75 µ G in magnitudeand (longitude 125 ◦ , latitude = 37 ◦ ) in direction (HGI 2000). Thefluctuations of the heliospheric TS and the heliopause with time(in Voyager 1/2 directions) are shown in Izmodenov & Alexashov(2020, Fig. 2). The complete description of the model can be foundin Izmodenov & Alexashov (2015, 2020).For the charged particles, the models imply a single-fluid ap-proach, so “plasma" includes SW/LISM protons, pickup protons,electrons, α particles in SW and helium ions in LISM. In the simu-lations, the velocity distribution function moments of plasma (num-ber density n , the bulk velocity vector V , and kinetic temperature T ) have been calculated on specific non-regular moving grid thatallows to perform exact fitting of the TS and heliopause (for details,see Izmodenov & Alexashov 2015). Afterward, for the sake of sim-plicity, the kinetic moments were interpolated on the spherical grid,which is irregular by radius. For the points inside the inner bound-ary (1 au) we extrapolate the plasma solution with assumptions of(1) ∝ / r proportionality for the number density, (2) ∝ / r ( γ − ) proportionality for the temperature (adiabatic law, γ = / r = 70 au. Insidethe boundary sphere, the H distribution function was calculatedby solving the kinetic equation with the method of characteristics,which allows taking into account non-Maxwellian properties of theH distribution in the vicinity of the Sun; the heliolatitude and timevariations of the SW are considered as well. Outside the boundarysphere, the velocity distribution function of H atoms is assumed tobe the sum of anisotropic (3 component) Maxwellian distributionsof primary (population 4) and secondary (population 3) H atomsonly. By that, we ignore the production of PUIs due to the ionizationof minor populations (by its relative abundance; see, Malama etal. 2006; Izmodenov et al. 2009) of H atoms that originated inthe region of supersonic SW and the IHS (populations 1 and 2,respectively). The validity of such an assumption will be discussedin Section 5 with results.Therefore, starting from this point, we assume that the globaldistributions of H atoms and plasma are known. All the followingresults were obtained using the distributions of plasma and H atomscalculated in the frame of the IA2020 model in the stationary caseunless otherwise indicated.
In this section, we describe the results of calculations of PUI dis-tribution in the heliosphere. The simulations were performed usingthe charged particles partitioning technique described in Section 4.Figure 1 presents the evolution of the PUI distribution functionwith distance. Numerical calculations were performed using thestationary version of the
IA2020 model at different heliosphericdistances ( r = 40, 75, 90, 110 au) in the upwind (or Nose) direction(longitude λ = 255.4 ◦ and latitude β = 5.2 ◦ in ecliptic coordinatesJ2000). The velocity distribution function of PUIs is assumed to beisotropic and it depends on the velocity w in the plasma referenceframe.The black dotted line of Figure 1 presents the analytical solu-tion for PUI distribution function in the region of supersonic SWderived by Vasyliunas & Siscoe (1976), or the so-called filled shelldistribution (see, also, Zank et al. 2010): f ∗ pui ( r , θ, w ) = π N V V (cid:18) V sw w (cid:19) / λ r exp (cid:32) − λ r θ sin θ (cid:18) V sw w (cid:19) / (cid:33) , where w < V sw , r is the heliocentric distance, θ is the angle betweenthe upwind direction and radius vector r , w is the velocity of PUIin the plasma reference frame, N is the H atoms number density MNRAS , 1–15 (2020) eliospheric ENAs: modeling and comp. with IBEX-Hi −1 w/V SW,0 f * P U I , [ s / k m ] f *pui (r, w) evolution with distance in the upwind direction [V & S, 1976]: R = 40 auR = 40 auR = 75- au (before TS)R = 75+ au (after TS)R = 90 auR = 90 au (PUI originated in IHS)R = 110 au Figure 1.
The evolution of the PUI distribution function with distance in theupwind direction. The solid curves are the results of calculations at differentheliocentric distances ( r = 40, 75, 90, 110 au). The blue and red solid linesare the distribution function profiles just before and after the TS, respectively(TS is located at 75 au). The green dashed curve presents the distributionfunction at 90 au of PUIs that originated in IHS. The black dotted linepresents the analytical solution by Vasyliunas & Siscoe (1976) at r = 40 au. w is the velocity of PUI in the plasma reference frame, V sw , = 432 km s − . in the LISM, V is the velocity of H atoms relative to the Sun, V sw is the SW velocity, λ = r ν ion , E / V is ionization characteristicdistance, and ν ion , E is the ionization rate at r E = 1 au. To plot theblack dotted line of Figure 1 the following set of parameters wasused: r =
40 au, θ = N = 0.14 cm − , V = 23 km s − , V sw = 432 km s − , ν ion , E = . × − s − . From Figure 1 it can beconcluded that the calculations of the developed kinetic model ofPUI distribution (cyan solid line) reproduces the analytical solution(black dotted line) qualitatively well. The difference in quantities canbe explained by several simplifications made to derive the analyticalsolution, such as the assumption of a spherically expanding SW withconstant velocity, the neglect of a thermal spread in velocities ofSW protons and H atoms, a single neutral component assumption,the neglect of H atoms interaction with the heliospheric interface,etc. The blue solid line of Figure 1 is the distribution function justbefore (upstream) the TS (TS is located at 75 au). The differencebetween the maximal velocities of the cyan and blue curves, whichare ∼ V sw , , can be explained by the deceleration of the SW withdistance from the Sun. The red solid line of Figure 1 is the distribu-tion function profile right after (downstream) the TS. The transitionfrom blue to red curve represents the influence of the jump condition(14) on the PUI distribution function profile. As can be seen, afterthe TS crossing a PUI gets (cid:112) C ( s , ψ ) ≈ (cid:112) ( s + )/ ψ ≈ ◦ in the upwind region), so the fast (with w > V sw , ) PUIs exist.The green dashed curve presents the distribution function at r = 90 au of PUIs that originated only in the inner heliosheath(so-called “injected" PUIs; see, e.g., Zirnstein et al. 2014). Thecomparison of the green and orange curves, the latter of whichshows the distribution function of PUIs that originated both in theregion of the supersonic SW and IHS, demonstrates that the PUIsoriginated in the IHS have smaller velocity in the plasma referenceframe. In the IHS, the plasma flow is decelerated, and the relativevelocity between the plasma and H atom, which is parental to PUI,is smaller (compared to the region of the supersonic SW). With increasing radial distance from the Sun, the transitionfrom the red (75 au) to orange (90 au) and yellow (110 au) curvesis accompanied by a gradual decrease in the number of PUIs with w ∼ V sw , that originated in the supersonic SW and increase in thenumber of slow PUIs ( w ∼ . · V sw , ), “injected" in the IHS. Thedecrease can be explained by the extinction of PUIs (due to neu-tralization), especially in the IHS, and it is driven by the survivalprobability term (11). This effect is more pronounced for PUIs with w ∼ V sw , that originated in the supersonic SW since they travel thelonger time compared to the “injected" PUIs. The critical (or maxi-mal) velocity of PUIs in the inner heliosheath is w c ≈ (cid:112) C ( s , ψ ) V sw , (in the plasma reference frame), and it is different for the red, or-ange, and yellow curves since the corresponding streamlines inter-sect slightly different regions of the TS (due to the fact that thestructure of the heliosphere is perceptibly three-dimensional, so theupwind direction is not a streamline as it is in axisymmetric mod-els), while the shock compression factor s depend on the positionat the TS.Using the calculated velocity distribution function, the kineticmoments of PUIs can be obtained (see Equations 12 and 13), whichare presented in Figure 2. This figure shows the PUI moments up-stream the TS: the number density (plot D), the kinetic temperature(plot E), and the pressure (plot F), as functions of the sphericalangles θ , which is set in radial direction of the polar plots, and ϕ that is increasing in the counterclockwise direction. The angle θ iscounted from the upwind direction ( θ = ◦ direction is oppositeto V LISM ), and the angle ϕ is counted from the plane containingthe V LISM and B LISM vectors (so-called BV-plane; ϕ = ◦ , ◦ corresponds to BV-plane), such as projection of B LISM vector onthe ( θ = ◦ , ϕ = ◦ ) direction is negative.From Figure 2(D) we can see that the maximum of PUI numberdensity is in the Nose region, and that PUIs are concentrated in thesolar equatorial plane ( ϕ ≈ ◦ , ◦ ). It can be explained by thefact, that the number density of PUIs is proportional to (a) theH atom number density n H , which has its maximum in the noseand minimum in the tail directions, and (b) the charge exchangeionization rate ν ex , which has its minimum at solar poles (see, e.g.,Figure 2 from Katushkina et al. 2019). The temperature of PUIs(plot E) is highly correlated with the SW velocity, which is shown inthe plot C, and has its maxima in the solar pole directions since thefast SW originates from the coronal holes, which are mostly locatedat poles. The PUI pressure has its maximum in the upwind region(plot F).Figure 2 also presents the shape of the TS (plot A), which iselongated in the tail direction, and the upstream shock-normal angle ψ (in degrees) between the magnetic field direction and outwardnormal to the shock (plot B). As can be seen from Figure 2(B),the TS can be considered as perpendicular (i.e., ψ = ◦ ) at itsnose ( θ ≈ ◦ ) and tail ( θ ≈ ◦ ) parts only, and the TS has ablunt shape in the nose region since ψ > ◦ on the Starboard side( ϕ ≈ ◦ ) and ψ < ◦ on the Port side ( ϕ ≈ ◦ ). Therefore, itis important to note, that the jump condition at the TS (14) used inour calculations is not strictly correct everywhere at the TS, since itemploys the assumption of shock perpendicularity. The charged particles (plasma) distribution was obtained in theframe of the
IA2020 heliospheric model that considers the plasmacomponent in the context of an ideal MHD and imply a single-fluidapproach. Therefore, in the calculations, the plasma represents a
MNRAS , 1–15 (2020)
I. I. Baliukin et al. (A) R TS , [au] NSPSSPV1 V20 4590135 (B) , [ ]
NSPSSP0 4590135 (C) V sw , [km/s] (D) n pui , [10 cm ] NSPSSP0 4590135 (E) T pui , [10 K] NSPSSP0 4590135 (F) p pui , [10 Pa]
NSPSSP0 4590135
Figure 2. (A) Heliocentric termination shock distance R TS . (B) The upstream shock-normal angle ψ . (C) The SW velocity. (D) The PUI number density. (E)The PUI kinetic temperature. (F) The PUI pressure. All the parameters are presented upstream the TS as functions of the spherical angles θ , which is set inradial direction of the polar plots (in degrees), and ϕ (in degrees) that is increasing in the counterclockwise direction. The definitions of θ and ϕ angles arepresented in the text. The values were obtained in the frame of stationary version of Izmodenov & Alexashov (2020) model. NSP = north solar pole, SSP =south solar pole. mixture of SW/LISM protons, pickup protons, electrons, α parti-cles (He ++ ) in SW and helium ions (He + ) in LISM. For the heliumion component in the LISM and α particles in the SW, the continuityequations were solved separately. Then, the number density n rep-resenting a mixture of protons and electrons only was obtained as n = ( ρ − m He n He )/ m p , where ρ is the plasma density, n He denotesthe He + number density in the interstellar medium and the He ++ number density in the SW (for details, see Izmodenov & Alexashov2015).To calculate the PUI distribution function, the distribution ofSW protons should be known, according to Equation (6). We havedeveloped a technique to separate charged particles into compo-nents, or, to be more precise, to calculate the number density andkinetic temperature of both proton components (core SW protonsand PUIs) in the IHS. This method is based on the following as-sumptions: • The number density n from the global simulations representsa composition of SW protons, PUIs, and electrons. To be morespecific, n = n sw + n pui + n e m e / m p ≈ n sw + n pui , since m e / m p ≈ × − ( m e is the mass of electron). • The plasma is quasi-neutral, which leads to equation n e = n sw + n pui + n He ++ = n + n He ++ . • All the populations of charged particles are co-moving, i.e. V sw = V pui = V e = V He ++ = V , where V is the plasma bulkvelocity from the global calculations. • For all the components the distribution functions are isotropic.The total pressure of the ionized component is equal to the sum ofpartial pressures, i.e. p = p sw + p pui + p e + p He ++ . Therefore, the following equation can be derived ( n sw + n pui + n e + n He ++ ) T = ( n + n He ++ ) T == n sw T sw + n pui T pui + n e T e + n He ++ T He ++ , where T is the plasma temperature from the global modeling. • The number density of alpha particles is 3.5% of the protonnumber density, i.e. n He ++ = α ( n sw + n pui ) = α n , where α = .
035 (Izmodenov & Alexashov 2015). This assumption is notstrictly correct in the whole region of the SW because near theSun there are sources due to the ionization of the helium atoms.However, the difference rapidly decreases with distance from theSun and becomes insignificant in the IHS, in the region of thespecific interest. Also, we assume that T He ++ = T sw . • The temperature of electrons T e = β T sw , where β = β = . n sw = n − n pui , T sw = ( + α ) nT − n pui T pui ( + β + ( + β ) α ) n − n pui , (16)The velocity distribution function f sw of SW protons is as-sumed to be isotropic Maxwellian: f sw ( t , r , v ) = n sw ( c sw √ π ) exp (cid:32) − ( v − V sw ) c (cid:33) , where c sw = (cid:112) k B T sw / m p is the SW thermal velocity. MNRAS , 1–15 (2020) eliospheric ENAs: modeling and comp. with IBEX-Hi −4 −3 N m b e − d e n s i t y , [ c m − ] T e r m i n a t i o n S h o c k Charged particles moments in the upwind direction.
PlasmaPUIsSW protonsElectronsα particles10 T e m p e r a t u r e , [ K ] H e li o p a u s e Adiabatic law (T sw ∼ n γ−1sw )0 10 20 30 40 50 60 70 80 90 100 110 R, [a0] P − e ss − e , [ P a ] , ,0i + , s1 + , e + , He ++ Figure 3.
Charged particles moments in the upwind direction. The top panelpresents number densities, the middle panel – temperatures, the bottompanel – pressures. The plasma moments (from global simulations) havebeen plotted with black curves, PUIs – red curves and dots, SW protons– blue curves and crosses, electrons – cyan curves and dots, α particles –magenta dashed curve. For the sum of partial pressures (PUIs, SW protons,electrons, and α particles) the yellow curve was used. To calculate the moments of SW protons ( n sw , T sw ) and PUIs( n pui , T pui ) the following iterative algorithm was used.(i) Since the PUI distribution is unknown, initially assumed thatthere are no PUIs ( n pui = ν ex , pui =
0) and, according to Equation(16), n sw = n , T sw = + α + β + ( + β ) α T . As can be seen, T sw = T in case of β = n pui , T pui ) are calculated using Equations (12)and (13).(iii) Using the formulas (16) the parameters ( n sw , T sw ) can berecalculated.(iv) Steps ii and iii were repeated until the convergence is ob-served. During our numerical experiments, we have found that start-ing from the second iteration the number density and temperatureof PUIs change insignificantly. Therefore, we have concluded thatthe calculations using the initial assumptions (described in Step 1)approximate the genuine distribution well.Figure 3 presents the results of the calculations of the iterativealgorithm (described above) in the upwind direction. The calcula-tions were performed for the plasma and neutral component distri-butions obtained in the frame of IA2020 heliospheric model in the stationary case. The top panel of the figure presents the number den-sities, the middle panel – the temperatures, and the bottom panel –the pressures. The plasma distribution (from the global simulations)has been plotted with black curves, the red curves with dots wereused for the PUIs, blue curves with crosses – for the SW protons,cyan curves with dots – for the electrons, magenta dashed curves– for the α particles, and the yellow curve was used for the totalpressure of PUIs, SW protons, and electrons.As can be seen from Figure 3, in the region of the supersonicSW the number density of PUIs remains almost constant with dis-tance, which can be explained by the fact the production (due toionization) and loss (due to neutralization) rates nearly compensateeach other. The temperature of PUIs in this region is decreasingslowly with radial distance due to the decrease of SW velocity. Inthe IHS, the number density of PUIs is growing up with increas-ing radial distance since the production rate of the “injected" PUIsdominates over the extinction. The temperature of PUIs in the in-ner heliosheath is decreasing with increasing radial distance sincethe fast PUIs (with w ∼ V sw , ), which originated in the supersonicSW, extinct rapidly, making the distribution function more centeredaround small values of velocities (see Figure 1).The temperature of the SW protons (see blue curve in the mid-dle panel) decreases adiabatically up to 10 AU where it becomesvery low ( ∼ − K). According to Malama et al. (2006),in the region from 10 AU up to the TS the energy transferred tothe SW by photoelectrons becomes non-negligible, which results inthe formation of such plateau in the spatial distribution of the SWprotons temperature. Downstream the TS, the temperature of theSW protons T sw ≈ × K (consistent with
Voyager 2 observa-tions; Richardson et al. 2008b), which corresponds to SW thermalvelocity c sw = (cid:112) k B T sw / m p ≈
58 km/s ( ≈ T sw is almost constant with increasing radial distance. Since the SWthermal velocity c sw is smaller compared to the SW velocity in theinner heliosheath ( V sw ≈
100 – 150 km s − ), the charge exchangeof the SW protons appears to be an insignificant contributor to themeasured ENA fluxes at high energies ( (cid:38) β = . β = ≈ Voyager 2 observations.The sum of the partial pressures p pui + p sw + p e + p He ++ (theyellow curve in the bottom panel) equals the initial plasma pres-sure from the global simulations (black curve), which verifies thecorrectness and accuracy of our partitioning technique. For the comparison with our modeling results, we use the datasets of globally distributed flux observed by
IBEX-Hi during thefirst 5 years of the mission (from 2009 to 2013) and presentedby Schwadron et al. (2014), which are available on the webpageof the
IBEX public Data Release 8 ( http://ibex.swri.edu/ibexpublicdata/Data_Release_8/ ). We have calculated modelfull-sky maps of the ENA fluxes as described in Sections 2–4 andperformed the comparison with
IBEX-Hi data at the energy channels2–6 (with the central energies ∼ MNRAS000
IBEX-Hi data at the energy channels2–6 (with the central energies ∼ MNRAS000 , 1–15 (2020)
I. I. Baliukin et al. corrections were applied for the
IBEX-Hi data, so in our calculations,we did not take into account both the relative motion of the
IBEX spacecraft and the ionization losses of ENAs in the region of thesupersonic SW. Also, the model maps of ENA fluxes were obtainedfor the lines of sight when the sensor views the heliosphere in theram direction (i.e., in the direction of spacecraft motion), as it wasdone for the
IBEX-Hi data in Schwadron et al. (2014).In previous works on modeling the
IBEX-Hi observations, ascaling of the simulated values is applied. In this way, to perform acomparison with the
IBEX-Hi data, the scaling factors 2.5 and 1.8were used in Zirnstein et al. (2017) and Kornbleuth et al. (2020),respectively. We have estimated the scaling factor for our modelingresults (based on the χ minimization procedure; see Appendix B)and found it to be equal to 1.002 (for the time-dependent modelcalculations), which is close to 1.0, so we do not scale our modelresults, unless otherwise indicated.The results of comparison of the IBEX-Hi data with model arepresented in Figures 4 – 7. Figures 4 and 5 present the full sky mapsin ENA fluxes for
IBEX-Hi energy channels, while Figures 6 and 7show the ENA spectra in the specific directions of the sky.Figures 4 and 5 present the comparison of the
IBEX-Hi data atthe top five energy channels (2–6) with the ENA flux maps obtainedin the frame of the
IA2020 heliospheric model. The first column ofthe figures present the results of stationary calculations, the secondcolumn – the time-averaged (during 2009–2013) calculations usingthe time-dependent version of the model, the third column is the
IBEX-Hi data collected during the same time period. The mapsof Figure 4 are centered on the Nose longitude 255 . ◦ and zerolatitude, while the maps of Figure 5 – on the Tail longitude 75 . ◦ and zero latitude. The stationary model fluxes are multiplied bycorresponding best-fitting scaling factor 0.67. Let us note, that theENA fluxes are calculated for the specific IBEX-Hi energy channelswith the energy response functions of electrostatic analyzers (ESAs)taken into account according to Equation (A1). The modeled ENAsoriginate from the PUI population only. We do not present in themaps the fluxes of ENAs that originated from the SW protons,since this component provide negligible fluxes in the whole
IBEX-Hi energy range (2-3 orders of magnitude smaller as can be seenfrom Figure 6).In Figure 6 the ENA flux spectra as it was observed by
IBEX-Hi in the upwind direction are presented. The black solid line withcrosses shows the
IBEX-Hi data. The yellow solid line presents thecalculated spectrum of ENAs that originated from PUIs (using thetime-dependent model). The red line with dots is also the simu-lated ENA fluxes of the time-dependent model, but calculated forthe specific
IBEX-Hi energy channels with energy response of theinstrument taken into account according. The blue solid line is themodel fluxes produced by the ENAs that originated from the SWprotons (the stationary model was used). The green dashed curvepresents the stationary model fluxes of ENAs that originated fromPUIs, which were born in the IHS (“injected" PUIs). The cyansolid line is the ENA spectrum (from PUIs) calculated using thestationary version of the global heliospheric model.Figure 7 shows the ENA flux spectra as it was observed by
IBEX-Hi in the directions of the North/South heliotail regions withenhanced fluxes (plot A) and the Starboard/Port heliospheric flankswith low fluxes (plot B), or the so-called North/South and Star-board/Port lobes. In our study, these directions have taken as thecenter bin directions (for which the data is present) that are theclosest to the lobe directions estimated by Zirnstein et al. (2016,see its Table 1) at the fifth energy step ( ∼ Table 1.
The directions of the heliospheric lobes in ecliptic (J2000) cordi-nates. Lobe Ecliptic Eclipticlongitude [ ◦ ] latitude [ ◦ ]North 75 45South 87 -45Starboard 153 15Port 9 -15 solid lines with crosses present the IBEX-Hi data, while the red andblue solid curves shows the model spectra in the North/Starboardand South/Port directions, respectively. The orange and cyan dashedcurves with dots are the simulated ENA fluxes for the specific
IBEX-Hi energy channels (the energy response of ESAs was taken intoaccount).From the comparison of the simulation results with the
IBEX-Hi data we can make the following conclusions: • There is a good quantitative agreement between the time-dependent model results and the observed fluxes in the middlerange of energies (at energy steps 3 – 5), especially for the regionsof North/South heliotail lobes where the absolute values are wellreproduced by the model even at the energy channel 6. Let us ad-ditionally note that for the simulation results, obtained in the frameof our time-dependent model, a renormalization is not needed (thescaling factor is very close to 1). • The time-dependent version of the model explains the
IBEX-Hi data better than the stationary model, especially in the Tail region(see the first and second columns of Figures 4 and 5). From thecomparison of the yellow and cyan curves of Figure 6, we canalso see that accounting for time-dependence “flattens" the ENAspectrum by making the fluxes lower at low ( (cid:46) (cid:38)
IBEX-Hi data. This can be explained by the fact that thespectrum, calculated using the time-dependent version of the model,reflect the averaged plasma properties in the IHS (during 2009 –2013), when the SW speed was different, which creates the spreadin the ENA spectrum. • Both the stationary and time-dependent models are able toreproduce qualitatively the geometry of the multi-lobe structureseen in the
IBEX-Hi data.A single structure of enhanced fluxes in the Tail direction is seenin the
IBEX-Hi data at lower energy channels ( ∼ ∼ IBEX-Hi data is seen in the model results. The Ion and Neutral Camera (
INCA )on the
Cassini spacecraft (Krimigis et al. 2009), which provided
MNRAS , 1–15 (2020) eliospheric ENAs: modeling and comp. with IBEX-Hi ~ . k e V
60 S30 S030 N60 N30 90150210270330
Model ST x 0.67
NoseTail V1V2 60 S30 S030 N60 N30 90150210270330
Model TD
NoseTail V1V2 60 S30 S030 N60 N30 90150210270330
IBEX-Hi GDF (2009-2013)
NoseTail V1V2 ~ . k e V
60 S30 S030 N60 N30 90150210270330 NoseTail V1V2 60 S30 S030 N60 N30 90150210270330 NoseTail V1V2 60 S30 S030 N60 N30 90150210270330 NoseTail V1V2 ~ . k e V
60 S30 S030 N60 N30 90150210270330 NoseTail V1V2 60 S30 S030 N60 N30 90150210270330 NoseTail V1V2 60 S30 S030 N60 N30 90150210270330 NoseTail V1V2 ~ . k e V
60 S30 S030 N60 N30 90150210270330 NoseTail V1V2 60 S30 S030 N60 N30 90150210270330 NoseTail V1V2 60 S30 S030 N60 N30 90150210270330 NoseTail V1V2 ~ . k e V
60 S30 S030 N60 N30 90150210270330 NoseTail V1V2 60 S30 S030 N60 N30 90150210270330 NoseTail V1V2 60 S30 S030 N60 N30 90150210270330 NoseTail V1V2
Figure 4.
The Mollweide skymap projections (in ecliptic J2000 coordinates) of the ENA fluxes as it was observed by
IBEX-Hi at the energy channels 2–6(by rows). The modeled ENAs originate from the PUIs only. The first column presents results of calculations using the stationary model by Izmodenov &Alexashov (2020), the second column – the model results averaged over 2009–2013 using the time-dependent model, the third column is the
IBEX-Hi datacollected during the same period. The stationary model fluxes are multiplied by factor 0.67. The units of fluxes are ( cm sr s keV ) − . The maps are centered onthe Nose ecliptic longitude 255 . ◦ and zero latitude. measurements of ENAs at high energies (5.2 – 55 keV), has alsoobserved these areas that were called “basins" (see, e.g., Dialynas2013). These low fluxes lobes are located in the vicinity of the solarequatorial plane, where the slow SW dominates and the thicknessof the heliosheath is small. As it was suggested by McComas et al.(2013) and studied by Zirnstein et al. (2016), the side lobes areformed by the composition of the following effects: (a) the closerthe LOS to the upwind direction is, the thinner the IHS, producinglower flux; (b) in the Nose of the heliosphere, the enhancementof the ENA flux is forced by the compression and heating; (c) theemission of faster SW at high latitudes produce lobes at the flanksof the heliosphere. From Figure 7(B) we can see that the fluxes fromthe Port side are systematically lower than in the Starboard region,and the model reproduces such behavior, which is due to the smallerheliosheath thickness on the Port side of the heliosphere.As it is seen in Figure 7, the slope of the IBEX-Hi spectrum ismuch smaller in the North/South lobe directions (plot A) than in theStarboard/Port lobe directions (plot B). This behavior is reproduced by the model also and can be explained by the fact that in thedirections of the low latitude side lobes the lines of sight intersectthose regions of the IHS, where the SW is slower and colder (withrespect to the North/South lobe directions), which results in lowerfluxes at high energies (see, also, McComas et al. 2013). • While the model produces a comparable with the
IBEX-Hi data fluxes at energy channels ∼ ∼ IBEX-Hi datavalues (black crosses) at the energy channels 5 and 6 are ∼ ∼
10 times, respectively, higher than the model fluxes (red points).The possible explanation for this discrepancy is a lack of ENAs,which originated from the energetic PUI population produced bythe processes of shock-drift acceleration or reflection from the cross-shock potential. The presence of such energetic PUIs is not takeninto account in our modeling. • The model produces smaller fluxes (compared to the
IBEX-Hi data) from both the Nose and Tail regions of the sky at energy
MNRAS000
MNRAS000 , 1–15 (2020) I. I. Baliukin et al. ~ . k e V
60 S30 S030 N60 N3090150210 270330
Model ST x 0.67
Nose TailV1 V2 60 S30 S030 N60 N3090150210 270330
Model TD
Nose TailV1 V2 60 S30 S030 N60 N3090150210 270330
IBEX-Hi GDF (2009-2013)
Nose TailV1 V2 ~ . k e V
60 S30 S030 N60 N3090150210 270330Nose TailV1 V2 60 S30 S030 N60 N3090150210 270330Nose TailV1 V2 60 S30 S030 N60 N3090150210 270330Nose TailV1 V2 ~ . k e V
60 S30 S030 N60 N3090150210 270330Nose TailV1 V2 60 S30 S030 N60 N3090150210 270330Nose TailV1 V2 60 S30 S030 N60 N3090150210 270330Nose TailV1 V2 ~ . k e V
60 S30 S030 N60 N3090150210 270330Nose TailV1 V2 60 S30 S030 N60 N3090150210 270330Nose TailV1 V2 60 S30 S030 N60 N3090150210 270330Nose TailV1 V2 ~ . k e V
60 S30 S030 N60 N3090150210 270330Nose TailV1 V2 60 S30 S030 N60 N3090150210 270330Nose TailV1 V2 60 S30 S030 N60 N3090150210 270330Nose TailV1 V2
Figure 5.
The description is the same as for Figure 4, but the maps are centered on the Tail ecliptic longitude 75 . ◦ . channel 2 (central energy ∼ IBEX-Hi data at this energy channel since theatoms of these populations have quite low energy ( (cid:46) ∼ V sw (for details see, Izmodenovet al. 2009) and potentially they could be a source of PUIs that willbe parental to ENAs with energies (cid:46) IBEX-Hi energy range these scenarios produce ENA fluxes, which are 1-2orders smaller than the fluxes from the populations considered in ourmodeling, so they can be safely neglected. According to the resultsof Malama et al. (2006) there are no other H atom populations that can explain the lack of ENAs in the modeling at
IBEX-Hi energychannel 2 ( ∼ ∼ ∼ The main limitation of the described above model is the absence ofprocesses that produce energetic “tails” in PUI distribution, such asshock-drift acceleration or reflection from the cross-shock potential(“shock-surfing” mechanism). To see how the inclusion of additionalenergetic PUIs affects the modeled ENA fluxes, we consider a “toymodel" with a power-law “tail" in PUI distribution downstream theTS. To simulate the additional energetic PUI population, the fol-lowing approach will be employed. Using the method described in
MNRAS , 1–15 (2020) eliospheric ENAs: modeling and comp. with IBEX-Hi −1 Energy, (keV)10 −1 F l u x , ( c m s r s k e V ) − E - s t e p E - s t e p E - s t e p E - s t e p E - s t e p E - s t e p ENA spectra in the upwind direction
IBEX-Hi dataModel TD (2009-2013). ENA from PUIModel TD (2009-2013). ENA from PUI, T(E)Model ST. ENA from SW protonsModel ST. ENA from PUI (originated in IHS)Model ST. ENA from PUI
Figure 6.
The ENA flux spectra as it was observed by
IBEX-Hi in the upwind direction. The black solid line with crosses presents the
IBEX-Hi data. The yellowsolid line presents the calculated (using time-dependent model) spectrum of ENAs that originated from PUIs. The red line with dots is also the simulated fluxesof the time-dependent model, but for the specific
IBEX-Hi energy steps with energy response taken into account. The blue solid line is the fluxes produced bythe ENAs that originated from the SW protons (stationary results). The green dashed curve presents fluxes of ENAs that originated from PUIs that were bornin IHS (stationary model). The cyan solid line is the ENA spectrum (from PUIs) calculated using the stationary model. The calculations were performed usingthe Izmodenov & Alexashov (2020) model. TD = time-dependent, ST = stationary.
Section 3, the PUI velocity distribution function can be calculatedeverywhere downstream the TS. This distribution is the filled shell f ∗ sh ( t , r , w ) with critical velocity w c ≈ (cid:112) C ( s , ψ ) · | V sw − V H | , i.e. f ∗ sh ( t TS , r TS , w ) = w > w c , where s is the shock compressionfactor, V sw is the SW velocity vector, V H is H atom bulk velocity, t TS and r TS = r ( t TS ) are the moment and position of a PUI crossingthe TS. For velocities higher than w c we assume the power law dis-tribution f ∗ tail ( t TS , r TS , w ) ∼ w − η with index η that defines the slopeof the “tail". We also introduce the additional parameter ξ that isthe density fraction of the PUIs of the “tail" distribution. Therefore,the PUI velocity distribution function right after the TS is assumedas f ∗ pui , d = ( − ξ ) f ∗ sh + ξ f ∗ tail , (17)where f ∗ tail ( t TS , r TS , w ) = n pui , d ( t TS , r TS ) w − η π ∫ + ∞ w c w − η w d w = n pui , d ( t TS , r TS ) η − π w (cid:18) ww c (cid:19) − η , w ≥ w c , (18)and f tail ( t TS , r TS , w ) = w < w c , n pui , d is the PUI number density downstream the TS. Equation (18) implies that η > ξ and η , which were artificiallyintroduced in the approach described above, depend on the local TSproperties, such as shock-normal angle ψ and shock compressionfactor s , but for the sake of simplicity, we assume it to be constantin our study. In order to demonstrate the qualitative effect of theadditional energetic PUIs on the ENA fluxes, we have performedcalculations for the specific pair of parameters – ξ = . η = MNRAS000
Section 3, the PUI velocity distribution function can be calculatedeverywhere downstream the TS. This distribution is the filled shell f ∗ sh ( t , r , w ) with critical velocity w c ≈ (cid:112) C ( s , ψ ) · | V sw − V H | , i.e. f ∗ sh ( t TS , r TS , w ) = w > w c , where s is the shock compressionfactor, V sw is the SW velocity vector, V H is H atom bulk velocity, t TS and r TS = r ( t TS ) are the moment and position of a PUI crossingthe TS. For velocities higher than w c we assume the power law dis-tribution f ∗ tail ( t TS , r TS , w ) ∼ w − η with index η that defines the slopeof the “tail". We also introduce the additional parameter ξ that isthe density fraction of the PUIs of the “tail" distribution. Therefore,the PUI velocity distribution function right after the TS is assumedas f ∗ pui , d = ( − ξ ) f ∗ sh + ξ f ∗ tail , (17)where f ∗ tail ( t TS , r TS , w ) = n pui , d ( t TS , r TS ) w − η π ∫ + ∞ w c w − η w d w = n pui , d ( t TS , r TS ) η − π w (cid:18) ww c (cid:19) − η , w ≥ w c , (18)and f tail ( t TS , r TS , w ) = w < w c , n pui , d is the PUI number density downstream the TS. Equation (18) implies that η > ξ and η , which were artificiallyintroduced in the approach described above, depend on the local TSproperties, such as shock-normal angle ψ and shock compressionfactor s , but for the sake of simplicity, we assume it to be constantin our study. In order to demonstrate the qualitative effect of theadditional energetic PUIs on the ENA fluxes, we have performedcalculations for the specific pair of parameters – ξ = . η = MNRAS000 , 1–15 (2020) I. I. Baliukin et al.
Energy, [keV] −1 F l u x , ( c m s r s k e V ) − E - s t e p E - s t e p E - s t e p E - s t e p E - s t e p E - s t e p (A) North and South lobes IBEX-Hi data (North)IBEX-Hi data (South)Model (North)Model (North), with T(E)Model (South)Model (South), with T(E)
Energy, [keV] E - s t e p E - s t e p E - s t e p E - s t e p E - s t e p E - s t e p (B) Starboard and Port lobes IBEX-Hi data (Starboard)IBEX-Hi data (Port)Model (Starboard)Model (Starboard), with T(E)Model (Port)Model (Port), with T(E)
Figure 7.
The ENA flux spectra as it was observed by
IBEX-Hi in the directions of the North/South heliotail lobes (plot A) and the Starboard/Port lobes (plotB). The black and grey solid lines with crosses present the
IBEX-Hi data, the red and blue solid curves shows the model spectra in the North/Starboard andSouth/Port directions, respectively. The orange and cyan dashed curves with dots are the simulated ENA fluxes for the specific
IBEX-Hi energy channels (theenergy response of ESAs was taken into account). The calculations were performed using time-dependent version of Izmodenov & Alexashov (2020) modeland the results were averaged over 2009–2013. The chosen directions of the lobes are presented in Table 1. stream the TS (plot A) and the ENA flux spectra as it was observedby
IBEX-Hi in the upwind direction (plot B), calculated with differ-ent assumptions on PUI distribution right after the TS. The yellowcurves present the results of calculations using the filled shell dis-tribution downstream the TS, and the green curve – the sum ofthe filled shell distribution and power-law “tail" (with parameters ξ = . η = IBEX-Hi energy channels withthe energy response of ESAs taken into account. The black solidline with crosses presents the data. Figure 9 shows the comparisonof the
IBEX-Hi data at the top five energy channels (2–6) with thesimulated full-sky ENA maps in the frame of the time-dependentversion of
IA2020 heliospheric model with the “tail" in the PUIdistribution downstream the TS.As can be seen from Figures 8 and 9, the additional populationof energetic PUIs (simulated in our approach using the power-law“tail" in the PUI distribution) produce higher fluxes at the top en-ergy channels, which makes the model flux maps qualitatively andquantitatively more consistent with the
IBEX-Hi data. As can beconcluded, the GDF is extremely sensitive to the form of the veloc-ity distribution function of PUIs in the inner heliosheath, and theaccounting for the existence of additional energetic population ofPUIs is essential to explain the data. Therefore, a detailed paramet-ric study of the this population using the
IBEX-Hi data needs to beperformed, which is beyond the scope of this paper. It is planned tobe done in the future and will be published elsewhere.
In this work, we have calculated the ENA fluxes at EarthâĂŹs orbitand performed a detailed quantitative comparison with the
IBEX-Hi data. The main conclusions of these studies can be summarized asfollows.(i) In the model described in this paper, the PUI population isconsidered kinetically. Using the developed model, we were able tocalculate the full-sky ENA flux maps and reproduce the geometryof the multi-lobe structure seen in the
IBEX-Hi data. There is a goodquantitative agreement between the time-dependent model resultsand the observed fluxes in the middle range of energies (at energysteps 3 – 5), especially for the regions of North/South heliotail lobes,where the absolute values are well reproduced by the model evenat the energy channel 6. For the time-dependent model results ascaling of the fluxes is not needed.(ii) Despite a relatively good agreement, there are few quantita-tive differences between our model calculations and the
IBEX-Hi data: (a) a deficit of fluxes at highest energy channels 5 and 6 isobserved, especially from the Nose region; (b) the model producessmaller fluxes (compared to the
IBEX-Hi data) from both the Noseand Tail regions of the sky at energy channel 2 (central energy ∼ ∼ ∼ MNRAS , 1–15 (2020) eliospheric ENAs: modeling and comp. with IBEX-Hi −1 w / V SW,0 −1 f * P U I , ( s / k m ) (A) f *pui (w) downstream TS in the upwind Filled shell (ξ = 0)Fi))ed she)) + T i) (ξ = 0.3, η = 5) ∼ w E+e.gy, ((eV) F ) , ( c m s . s ( e V ) (B) ENA s−ec0. i+ 0he 1−2i+d di.ec0i,+ IBEX-Hi dataFilled shell TDFilled shell TD, with T(E)Filled shell + Tail (ξ = 0.3, η = 5), TDFilled shell + Tail (ξ = 0.3, η = 5), TD, with T(E)
Figure 8. (A) The velocity distribution function of PUIs downstream the TS and (B) the ENA flux spectra as it was observed by
IBEX-Hi in the upwinddirection, calculated with different assumptions on PUI distribution right after the TS: the yellow curve – filled shell distribution; the green curve – the sum ofthe filled shell distribution and power law “tail" with parameters ξ = . η = IBEX-Hi energy channels. The black solid line with crosses presents the
IBEX-Hi data. The simulationsof the PUI distribution functions (plot A) were performed in the stationary case, and the spectra (plot B) are calculated using the time-dependent model byIzmodenov & Alexashov (2020). V sw , = 432 km s − . TD = time-dependent.
60 S30 S030 N60 N30 90150210270330
IBEX-Hi GDF (2009-2013)
NoseTail V1V2 ~ . k e V
60 S30 S030 N60 N30 90150210270330
Model with tail TD x 1.04
NoseTail V1V2 60 S30 S030 N60 N30 90150210270330 NoseTail V1V2 ~ . k e V
60 S30 S030 N60 N30 90150210270330 NoseTail V1V2 60 S30 S030 N60 N30 90150210270330 NoseTail V1V2 ~ . k e V
60 S30 S030 N60 N30 90150210270330 NoseTail V1V2 60 S30 S030 N60 N30 90150210270330 NoseTail V1V2 ~ . k e V
60 S30 S030 N60 N30 90150210270330 NoseTail V1V2 60 S30 S030 N60 N30 90150210270330 NoseTail V1V2 ~ . k e V
60 S30 S030 N60 N30 90150210270330 NoseTail V1V2
60 S30 S030 N60 N3090150210 270330
IBEX-Hi GDF (2009-2013)
Nose TailV1 V2 ~ . k e V
60 S30 S030 N60 N3090150210 270330
Model with tail TD x 1.04
Nose TailV1 V2 60 S30 S030 N60 N3090150210 270330Nose TailV1 V2 ~ . k e V
60 S30 S030 N60 N3090150210 270330Nose TailV1 V2 60 S30 S030 N60 N3090150210 270330Nose TailV1 V2 ~ . k e V
60 S30 S030 N60 N3090150210 270330Nose TailV1 V2 60 S30 S030 N60 N3090150210 270330Nose TailV1 V2 ~ . k e V
60 S30 S030 N60 N3090150210 270330Nose TailV1 V2 60 S30 S030 N60 N3090150210 270330Nose TailV1 V2 ~ . k e V
60 S30 S030 N60 N3090150210 270330Nose TailV1 V2
Figure 9.
The Mollweide skymap projections (in ecliptic J2000 coordinates) of the ENA fluxes as it was observed by
IBEX-Hi at the energy channels 2–6(by rows). The modeled ENAs originate from the PUIs only. The first and third columns present the results of calculations using the time-dependent model byIzmodenov & Alexashov (2020) with the “tail" in the PUI distribution downstream the TS ( ξ = . , η = IBEX-Hi data collected during 2009 – 2013. The model fluxes were multiplied by the best-fitting scaling factor ˆ k = 1.04. The maps in the first and second columns arecentered on the Nose longitude 255 . ◦ , and the maps in third and fourth columns – on the Tail longitude 75 . ◦ .MNRAS000
IBEX-Hi at the energy channels 2–6(by rows). The modeled ENAs originate from the PUIs only. The first and third columns present the results of calculations using the time-dependent model byIzmodenov & Alexashov (2020) with the “tail" in the PUI distribution downstream the TS ( ξ = . , η = IBEX-Hi data collected during 2009 – 2013. The model fluxes were multiplied by the best-fitting scaling factor ˆ k = 1.04. The maps in the first and second columns arecentered on the Nose longitude 255 . ◦ , and the maps in third and fourth columns – on the Tail longitude 75 . ◦ .MNRAS000 , 1–15 (2020) I. I. Baliukin et al. the physical reasons for the lack of ENA fluxes at energy channel2 ( ∼ ACKNOWLEDGEMENTS
The authors would like to acknowledge Dr. Nathan Schwadron forproviding information on IBEX GDF uncertainties. The work ofI.I. Baliukin and V.V. Izmodenov was supported by the Founda-tion for the Advancement of Theoretical Physics and Mathematics“BASIS" 18-1-1-22-1. The numerical single-fluid MHD code of theSW/LISM global interaction used in this paper was developed byD.B. Alexashov in the frame of the Russian Science Foundationgrant 19-12-00383. The authors would like to thank for the discus-sions that appeared during online meetings in the frame of NASA18-DRIVE18_2-0029, Our Heliospheric Shield, 80NSSC20K0603project. The authors are no financial support collaborators of theDRIVE project.
DATA AVAILABILITY
The data underlying this article will be shared on reasonable requestto the corresponding author.
REFERENCES
Chalov, S. V., & Fahr, H. J. 1997, A&A, 326, 860-869Chalov, S. V., Fahr, H. J., & Izmodenov, V. V. 2003, J. Geophys. Res., 108,1266Chalov, S. V., Alexashov, D. B., McComas, D., et al. 2010, ApJL, 716, L99Chalov S.V., & Fahr H.J., 2013, MNRAS, 433, L40Chalov S.V., Malama, Y. G., Alexashov, D.B., Izmodenov V.V., 2015, MN-RAS, 455, 431-437Chalov, S.V. 2019, MNRAS, 485, 4, 5207-5209Dialynas, K., Krimigis, S. M., Mitchell, D. G., Roelof, E. C., & Decker, R.B. 2013, ApJ, 778, 40Fahr, H.J., & Chalov, S.V. 2008, A&A, 490, L35-L38Fahr, H. J., & Fichtner, H. 2011, A&A, 533, A92+Fahr, H. J., & Siewert, M. 2011, A&A, 527, A125Fahr, H. J., & Siewert, M. 2013, A&A, 552, A41Fahr, H. J., Sylla, A., Fichtner, H., and Scherer, K. 2016, J. Geophys. Res.Space Physics, 121, 8203-8214Fisk, L. A. & Gloeckler, G. 2007, Proceedings of the National Academy ofSciences, 104, 14, 5749-5754Funsten, H. O., Allegrini, F., Bochsler, P. et al., 2009, SSRv, 146, 75–103Gloeckler, G., Geiss, J., Roelof, E.C. et al. 1994, JGR: Space Physics, 99,A9, 17,637-17,643Heerikhuisen, J., Pogorelov, N. V., Florinski, V., Zank, G. P., & le Roux, J.A. 2008, ApJ, 682, 679Heerikhuisen, J., Pogorelov, N. V., Zank, G. P. et al. 2010, ApJ, 708, L126Isenberg, P. A. 1987, JGR, 92, 1067 Izmodenov, V. V., The Outer Heliosphere: The Next Frontiers, Edited by K.Scherer, Horst Fichtner, Hans Jörg Fahr, and Eckart Marsch COSPARColloquiua Series, 11. Amsterdam: Pergamon Press, 2001, 23Izmodenov, V.V., Malama , Y.G., Ruderman, M.S., et al., Kinetic-Gasdynamic Modeling of the Heliospheric Interface, SSRv, 2009, 146,329-351Izmodenov V.V., & Alexashov D.B., 2015, ApJS, 220, 32Izmodenov V.V., & Alexashov D.B., 2020, A&A, 633, L12Kallenbach, R., Hilchenbach, M., Chalov, S. V., Le Roux, J. A., and Bamert,K. 2005, A&A, 439, 1-22Katushkina, O. A., & Izmodenov, V. V. 2010, AstL, 36, 297Katushkina, O. A., Izmodenov, V. V., Alexashov, D. B., Schwadron, N. A.,and McComas, D. J. 2015, ApJS, 220, 33Katushkina, O. A., Izmodenov V. V., Koutroumpa, D., Quemerais, E., Jian,L. K., 2019, Solar Phys., 294, 17Kornbleuth M., Opher M., Michael A.T., Drake J.F. 2018, ApJ, 865, 84Kornbleuth M., Opher M., Michael A.T., et al. 2020, ApJL, 895, L26Kowalska-Leszczynska, I., Bzowski, M., Sokol, J. M., Kubiak, M. A. 2018,ApJ, 852, 2, 115Kowalska-Leszczynska, I., Bzowski, M., Kubiak, M. A., Sokol, J. M. 2020,ApJS, 247, 2, 62Krimigis, S. M., Mitchell, D. G., Roelof, E. C., Hsieh, K. C., & McComas,D. J. 2009, Science, 326, 971Lindsay, B. G., & Stebbings, R. F. 2005, JGRA, 110, A12213McComas, D. J., Allegrini, F., Bochsler, P., et al. 2009, Science, 326, 959McComas, D. J., Alexashov, D. B., Bzowski, M., et al. 2012, Science, 336,6086, 1291McComas, D. J., Dayeh, M. A., Funsten, H. O. et al. 2013, ApJ, 771, 77McComas, D. J., Bzowski, M., Frisch, P., et al. 2015, ApJ, 801, 28McComas, D. J., Bzowski, M., Dayeh, M. A., et al. 2020, ApJS, 248, 26Malama, Y. G., Izmodenov, V. V., Chalov, S. V. 2006, A&A, 445, 693-701Richardson, J. D., J. C. Kasper, C. Wang, J. W. Belcher, and A. J. Lazarus,2008, Nature, 464, 63-66.Richardson, J. D. 2008, GeoRL, 35, L23104Rucinski, D., Fahr, H. J., and Grzedzielski, S. 1993, Planet. Space Sci. 41,773Scherer, K., Fichtner, H., Fahr, H. J. 1998, J. Geophys. Res., 103, A2, 2105-2114Schwadron N. A., Allegrini F., Bzowski M. et al. 2011, ApJ, 731, 56Schwadron N.A., Moebius E., Fuselier S.A. et al., 2014, ApJS, 215, 13Shrestha, B. L., Zirnstein, E. J., Heerikhuisen, J. 2020, ApJ, 894, 102Vasyliunas, V.M., & Siscoe, G.L. 1976, J. Geophys. Res., 81, 7Witte, M. 2004, A&A, 426, 835 – 84Zank, G. P., Heerikhuisen, J., Pogorelov, N. V., et al. 2010, ApJ, 708, 1092Zank, G. P., Heerikhuisen, J., Wood, B., et al. 2013, ApJ, 763, 20Zirnstein, E. J., Heerikhuisen, J., Zank G. P. et al. 2014, ApJ, 783, 129Zirnstein, E. J., Funsten, H. O., Heerikhuisen, J. et al. 2016, ApJ, 826, 58Zirnstein, E. J., Heerikhuisen J., Zank G.P. et al. 2017, ApJ, 836, 238
APPENDIX A:
IBEX-HI
ENERGY TRANSMISSION
To calculate the differential flux J Mi in the LOS direction measuredby
IBEX-Hi at the energy channel J Mi ( LOS ) = ∫ E i , max E i . min j ENA ( E , LOS ) T i ( E ) dE , i = , ..., , (A1)where superscript “M" denotes the model, j ENA ( E , LOS ) is differ-ential spectra of ENA fluxes, E i , min and E i , max are the boundariesof ESA T i ( E ) is normalized energy responsefunction of the ESA ∫ E i , max E i , min T i ( E ) dE = IBEX-Hi sensor is not taken into ac-count in our modeling, so the fluxes are calculated for the centerdirections of the
IBEX skymap 6 ◦ × ◦ bins. MNRAS , 1–15 (2020) eliospheric ENAs: modeling and comp. with IBEX-Hi E n e r g y r e s p o n s e T ( E ) , [ e V ( ] IBEX-Hi energy response
ESA cen = 0.45 keVESA cen = 0.71 keVESA cen = 1.10 keVESA cen = 1.74 keVESA cen = 2.73 keVESA cen = 4.29 keV
Figure A1.
IBEX-Hi energy response as function of energy for all 6 energychannels. The calibration data files were taken from http://ibex.swri.edu/ibexpublicdata/CalData/Hi/ . APPENDIX B: MODEL SCALING FACTOR
The difference between the model and data can be described interms of the χ value (the weighted average of residuals): χ ( k ) = (cid:213) ESA i (cid:213) LOS j (cid:32) J Dij − k × J Mij σ ij (cid:33) , (B1)where k is the scaling factor, J Dij and J Mij are the
IBEX-Hi data andmodel ENA flux values (superscript “D" denotes the data), σ ij arethe uncertainties of observations, the summations are performed forthe top five IBEX-Hi energy channels ( i = , . . . ,
6) and all 60 × j = , . . . , k ,for which the χ ( k ) function takes its minimum, can be found asˆ k = (cid:205) ij J Mij J Dij / σ (cid:205) ij ( J Mij / σ ij ) , (B2)which is the solution of equation d χ / dk = χ , which is χ per degree offreedom, can be calculated as χ = χ / ν , where ν = N − M equalsthe number of observations N minus the number of fitted parameters M . In our study, the number of observations N = × = This paper has been typeset from a TEX/L A TEX file prepared by the author.MNRAS000