Interstellar Pickup Ion Production in the Global Heliosphere and Heliosheath
aa r X i v : . [ a s t r o - ph . S R ] S e p Draft version September 26, 2018
Preprint typeset using L A TEX style AASTeX6 v. 1.0
INTERSTELLAR PICKUP ION PRODUCTION IN THE GLOBAL HELIOSPHERE ANDHELIOSHEATH
Yihong Wu , Vladimir Florinski , and Xiaocheng Guo Center for Space Plasma and Aeronomic Research, University of Alabama in Huntsville, Huntsville, AL 35805, USA;[email protected] Department of Space Science, University of Alabama in Huntsville, Huntsville, AL 35899, USA
ABSTRACTInterstellar Pickup ions (PUIs) play a significant part in mediating the solar wind (SW)interaction with the interstellar medium. In this paper, we examine the details of spa-tial variation of the PUI velocity distribution function (VDF) in the SW by solvingthe PUI transport equation. We assume the PUI distribution is isotropic resultingfrom strong pitch-angle scattering by wave-particle interaction. A three-dimensionalmodel combining the MHD treatment of the background SW and neutrals with akinetic treatment of PUIs throughout the heliosphere and the surrounding local in-terstellar medium (LISM) has been developed. The model generates PUI power lawtails via second-order Fermi process. We analyze how PUIs transform across the he-liospheric termination shock (TS) and obtain the PUI phase space distribution in theinner heliosheath including continuing velocity diffusion. Our simulated PUI spectraare compared with observations made by
New Horizons , Ulysses , Voyager 1, 2 and
Cassini , and a satisfactory agreement is demonstrated. Some specific features in the
Wu et al. observations, for example, a cutoff of PUI VDF at v = V SW and a f ∝ v − tail in thereference frame of the SW, are well represented by the model. Keywords: acceleration of particles — plasmas — solar wind INTRODUCTIONA neutral atom from the local interstellar medium (LISM) can be ionized while drifting into theheliosphere and turn into an interstellar pickup ion (PUI). The neutral atom can become ionizedthrough charge exchange with an ion, photoionization by sun light or impact ionization by SolarWind (SW) electrons. Charge exchange is the dominant mechanism of the three. In the SW frame,a newly created ion is “picked up” and starts to gyrate about the direction of the magnetic field.Subsequently they may be scattered into an isotropic distribution by either ambient preexisting orself-excited waves (Vasyliunas & Siscoe 1976; Isenberg 1987; Zank 1999). The resulting PUI velocitydistribution function (VDF) in the reference frame of the SW is a spherical shell, at the lowest order.These PUIs are also convected outward by the expanding SW. As the PUIs travel outward throughthe heliosphere, the shell becomes filled by adiabatic cooling. Generally, one can easily identify PUIsby their VDFs which are distinctly different from that of SW ions.In view of the paucity of measurements of keV-ions in the outer heliosphere, the transport ofPUIs is not yet fully understood. Especially, the production of the suprathermal (high-energy non-Maxwellian) tails on the VDFs is a subject of much debate (Chalov et al. 2003; Fahr & Fichtner2011). It was suggested that during their propagation through the outer heliosphere, PUIs expe-rience pitch-angle scattering and stochastic acceleration by interactions with different kinds of SWturbulences (Chalov et al. 1997).PUI kinetic transport theory describes the evolving PUI distribution in phase space. The transportequation, an equation of “motion” for the distribution function, is the basis for almost all work
UIs in the Heliosphere and Heliosheath f ∝ v − . Intriligator et al. (2012) first comparedsimulated PUI densities and temperatures with the SWICS measurements on board Ulysses . Themodel results matched well with the observations at about 5 . κ distributions, Fahr et al. (2014) proceeded from the phase space transport equation to a pressureequation for the κ parameter κ = κ ( r ) as function of heliocentric distance r . They obtained therange of possible radial variations of κ from relatively high values in the inner heliosphere to valuesbetween 1.5 and about 2 further out, depending on diffusion coefficient.The purpose of the present paper is to place PUI transport in the context of the global picture ofthe SW-LISM interaction. The treatment is based on the more conventional view where wave-particleinteractions ensure rapid isotropization of the distribution. A parallel can be drawn with the approachof Gamayunov et al. (2012), who used a grid-based model for the isotropic PUI distribution function Wu et al. and for the waves generated by the PUI ring anisotropization process, but neglected PUI momentumdiffusion. Their model was applied to the supersonic SW upstream of the TS. Here we investigatethe detailed spatial distribution of PUIs in both the supersonic SW and the heliosheath. We takeinto account the effects of convection with the SW, adiabatic cooling, second-order Fermi processand ionization. Similar to Gamayunov et al. (2012), the model combines the MHD treatment of thebackground SW and neutral atoms with a kinetic treatment of PUIs in the isotropic approximation.The rest of this paper proceeds as follows. Section 2 explains how we simulate the SW-LISMinteractions. The PUI transport model is introduced in Section 3. Section 4 presents the simulatedPUI phase space density and energy spectra. Finally, Section 5 mentions the weaknesses of thepresented model and proposes several improvements to the current methods. THE SW-LISM INTERACTIONThe heliosphere is a low-density bubble embedded in the local interstellar cloud. The SW slowsdown when interacting with the LISM. At the TS, a standing shock wave, the SW speed falls belowthe effective speed of sound (that includes a contribution from PUIs) and becomes subsonic. TheTS causes compression, heating, and a change in magnetic field. The SW continues to slow down asit passes through the heliosheath, a transitional zone bounded by the TS and the heliopause, wherethe pressures from the LISM and SW are balanced. The LISM pressure causes the heliosphere todevelop into a comet-like structure. The nose of the heliopause is defined as the direction oppositeto the sun’s motion through the Local Interstellar Cloud.We use a three-dimensional MHD model with four neutral hydrogen atom populations to obtainthe plasma background. The neutral hydrogen atoms are separated into four species according tothe region in which they were created. Studies comparing multiple neutral fluid and the kineticMonte Carlo approaches to model neutrals have been done and the 4-fluid model shows excellent
UIs in the Heliosphere and Heliosheath . z -axis points northward of the solar equator, andthe x -axis is in the plane defined by the interstellar helium flow direction (Lallement et al. 2005) andthe z -axis. The x, y, z -axes constitute a right-handed orthogonal system.The model heliosphere corresponds approximately to the solar minimum conditions. We assumedthe slow wind latitudinal extent angle of 36 ◦ . At 1 AU, in the slow SW, the bulk velocity u = 430km s − and density n = 5 . − ; in the fast SW, u = 725 km s − and n = 1 . − (Ebert et al.2009; Jian et al. 2011). At 1 AU, the radial magnetic field component B r = 24 . µ G. The azimuthalmagnetic field component B φ ∝ /u , where u is the SW speed. Since the heliospheric currentsheet is too thin to affect the dynamics of the plasma on large scales, we do not include it in oursimulation and use a unipolar magnetic field model (Czechowski et al. 2010; Borovikov et al. 2011;Izmodenov & Alexashov 2015; Opher et al. 2015).The interstellar neutral atom density was 0 .
11 cm − and its velocity vector was ( − . , , − . − in our coordinates (McComas et al. 2012). We assumed the interstellar magnetic field di-rection toward the center of the Interstellar Boundary Explorer (IBEX) ribbon, with Cartesian
Wu et al. components (1 . , − . , . µ G (Funsten et al. 2013). The interstellar plasma and neutraltemperatures were both equal to 6300 K. The plasma interacts with the hydrogen atoms throughcharge exchange. These atoms are modeled using a gas-dynamic model in a fashion similar to ourMHD system except that the magnetic field is set to be zero. We use the charge exchange termsin Pauls et al. (1995). The MHD code is run until a steady state is reached, which provides theplasma-neutral simulation background for the next step. THE PUI TRANSPORT MODELThe PUIs are non-thermal, and their preferred treatment is kinetic. We assume that PUIs arescattered into isotropic distribution by either ambient preexisting or self-excited waves, so the VDFdepends only on the absolute value of velocity (in the plasma frame).Schwadron et al. (1996) obtained the mean field strength B from 1-day averages measured by Ulysses , and averaged the square of the variations over 5 days, η = h ( B − B ) /B i , during the 150days from February 19 to July 18, 1992. They found that Ulysses crossed flux tubes with differentvalues of η , in the range from a high value of η H ≈ .
12 (highly turbulent) to a low value of η L ≈ . η . Following thisidea, we separate the PUI distribution into two components. The first component (“core”) is foundin laminar flux tubes, while the second component (“tail”) resides in turbulent flux tubes, wherethey experience stochastic acceleration.The distributions f core and f tail satisfy the following transport equations: ∂f core ∂t + u · ∇ f core − ∇ · u ∂f core ∂ ln v = Γ ∗ S, (1) ∂f tail ∂t + u · ∇ f tail − ∇ · u ∂f tail ∂ ln v − v ∂∂v (cid:18) v D ∂f tail ∂v (cid:19) = (1 − Γ) ∗ S (2)with Γ = 0 .
98, where v is PUI velocity in the non-inertial frame moving with the plasma, u is UIs in the Heliosphere and Heliosheath D is the velocity diffusion coefficient, and S is the source term,describing PUIs produced from charge exchange or some other ionization process. The left-hand sidesof the transport equations describes the explicit time dependence of the PUI velocity distribution,convection, adiabatic cooling or heating and velocity diffusion, respectively. PUIs typically haveenergies between 100 eV and 100 keV. We neglect spatial diffusion and drift motions since the typicalspatial diffusive scale of PUI is much smaller than our simulation scale, and the drift speed is smallin this low energy range (Rucinski et al. 1993).The velocity diffusion coefficient is D = v δu/ (9 ξ c ), where δu is the RMS fluctuating velocity of theplasma and ξ c is the correlation length (Chalov et al. 1997). It includes a multitude of effects, such astransit-time damping (Fisk 1976), parallel electric fields in two-dimensional turbulence (le Roux et al.2002), and large-scale compressive structures (Chalov et al. 2003). The values of δu and ξ c were (29km s − , 0 . − , 1 AU) in the supersonic SW. Theseparameters were chosen to produce suprathermal tails in a way that would be consistent with the Voyager observations (see below).In the source term, charge exchange and photoionization are taken into account, but electronimpact ionization is neglected. The distribution of PUIs (variables bearing a subscript “i”) createdby the charge exchange process holds imprints from the parent atom distribution (index “n”) andthe background plasma (index “p”). We assume that both are Maxwellian with densities n , meanvelocities u and thermal speeds v T = (2 kT /m ) / , f p ( v p ) = n p π / v T p e − ( v p − u p ) /v Tp , (3) f n ( v n ) = n n π / v T n e − ( v n − u n ) /v Tn . (4)The PUI production term from charge exchange is (cid:18) δf i ( v i ) δt (cid:19) ex = Z σ ( | v p − v n | ) f p ( v p ) f n ( v n ) | v p − v n | d v p d v n d v i Wu et al. = σ (∆ v p ( v i )) f n ( v i ) Z f p ( v p ) | v p − v i | d v p (5)= n p σ (∆ v p ( v i ))∆ v p ( v i ) f n ( v i ) , where ∆ v p ( v i ) is the average relative speed between a given PUI and the proton population. Toobtain the above equation, the condition for charge exchange, v n = v i , and the experimental resultthat the cross section σ ( v ) is only weakly dependent on the relative velocity between the particles,were used. Numerical values of σ are available from Lindsay & Stebbings (2005). The average relativespeed is computed as (Fahr & M¨uller 1967)∆ v p ( v i ) = 1 π / v T p Z | v p − v i | e − ( v p − u p ) /v Tp d v p = v T p (cid:20) √ π e − x + (cid:18) x + 12 x (cid:19) Erf( x ) (cid:21) , (6)where x = | v i − u p | /v T p .The PUI transport equation is solved in the plasma frame, so we transform (5) into that frameand average over a sphere. We introduce the PUI velocity in the non-inertial frame moving with theplasma, v ′ i = v i − u p = x v T p . Notice that in this frame (6) is isotropic and requires no averaging.Then (cid:18) δf i ( v ′ i ) δt (cid:19) ex = n p n n σ (∆ v p ( v ′ i ))∆ v p ( v ′ i )4 π / v T n Z e − ( v ′ i − u n + u p ) /v Tn d Ω ′ i = n p n n σ (∆ v p ( x ))∆ v p ( x )2 π / v T n Z − e − ( α x − µ ′ αxh + h ) dµ ′ = n p n n σ (∆ v p ( x ))∆ v p ( x )4 π / v T p v T n e − ( αx − h ) − e − ( αx + h ) xh , (7)where x = v ′ /v T p , h = | u n − u p | /v T n and α = v T p /v T n , is the charge exchange source term for anisotropic distribution of PUIs. Similarly, the photoionization source term is (cid:18) δf i ( v ′ i ) δt (cid:19) ph = n n ν π / v T p v T n e − ( αx − h ) − e − ( αx + h ) xh , (8)where ν = ν ph (1 AU /r ) with ν ph being the rate of photoionization at 1 AU. UIs in the Heliosphere and Heliosheath ∂f∂t + ∇ · ( u f ) − v ∂∂v (cid:20) v ( ∇ · u )3 v f + Dv ∂f∂v (cid:21) = S, (9)and similar for equation (2). The transport module is implemented on the same grid as the MHD sim-ulation, using a finite volume method. Right and left interface values are computed using piecewiselinear reconstruction with a WENO limiter (Jiang & Shu 1996; Friedrich 1998). The ionic compo-nents all have the same bulk speed available from the MHD solution. The velocity grid extends from10 km s − to 6000 km s − . SIMULATION RESULTSThe TS and the heliopause are at 90 AU and 144 AU respectively along the
Voyager 1 spacecraftdirection in our simulations.
Voyager 1 actually crossed the TS and the heliopause at 94 AU and 122AU (Stone et al. 2005, 2013). These distances are appropriate for a model that is time independentand is based on solar-minimum conditions. In Figure 1, the left panel shows the simulated PUIVDF in the
Voyager 1 direction ( θ = 55 ◦ , φ = 0 ◦ , where θ denotes the co-latitude and φ denotes thelongitude), in the plasma frame. At r = 5 AU, the PUI distribution shows a rapid drop at around 430km s − , the bulk velocity of SW inside the TS, indicating that most particles are injected with thespeed of the SW V SW in the plasma frame. Some particles have filled in the shell at low energies dueto adiabatic cooling. Since particles with speeds v > V SW result from local acceleration in turbulentflux tubes, the mixture of core PUIs and tail PUIs yields this step like feature (Schwadron et al.1996). These distributions are to be interpreted in the time-averaged sense, over many flux tubespassing an observer.A power-law suprathermal tail develops at v > V SW . The tail extends to higher energies withincreasing radial distance. In view of the flux tube picture, an observer would alternately see distri-0 Wu et al. v (km/s)10 f ( s k m - ) -4 -2 r = 5 AUr = 30 AUr = 70 AUr = 100 AUr = 110 AU v (km/s)10 f ( s k m - ) -4 -2 r = 5 AUr = 30 AUr = 90 AUr = 120 AUr = 140 AU Figure 1 . PUI velocity distribution functions in the plasma frame along the
Voyager 1 direction ( θ =55 ◦ , φ = 0 ◦ ) (left panel) and along the north polar direction ( θ = 0 ◦ ) (right panel). butions with more or less developed tails; Figure 1 is their average over several days or even months(the model is not time dependent, therefore the averaging interval can be arbitrary long). A humpat around 100 km s − at 100 AU is clearly seen, which consists of the low-energy PUIs created in theinner heliosheath. The hump increases in height as the flow slows down toward the heliopause dueto the accumulation of low-energy PUIs produced in the heliosheath.The right panel of Figure 1 illustrates the variation of the PUI VDF along the polar direction( θ = 0 ◦ , φ = 0 ◦ ). The principle features are similar to the equatorial case. At r = 5 AU, the PUIdistribution shows a sharp change in slope around 725 km s − , the bulk velocity of the fast SW insidethe TS. A pronounced hump develops only beyond the TS (at about 90 AU). In both panels one cansee that the PUI gas gets compressed by the shock and an accelerated power law tail develops out ofthe PUI core. There is little change in the power-law tail in the heliosheath, which means stochasticacceleration cannot produce a spectrum that is any harder. The downstream value of the power lawindex is about − . UIs in the Heliosphere and Heliosheath Figure 2 . Plasma background density (cm − ) in the meridional plane. Distinct bands of fast and slow windin the high and low latitude regions are shown. The slow wind latitudinal extent angle is assumed to be36 ◦ . Fast SW prevailing at higher latitudes produces denser subsonic wind in the heliosheath. core and tail PUI spatial distributions. They show that from the inner boundary to the TS, thePUI density falls off slower than that of the SW core since PUIs are produced over the entire SW.One can clearly see that there is a great difference in PUI densities in the slow SW and in the fastSW. Figures 3(a) to 3(c) show that at low energies the PUI density is higher in the slow SW thanin the fast SW. The majority of the maps show an increase in PUI densities across the TS due tocompression. Naturally, PUIs created in a slow wind have energy lower than those produced in a fastwind. One can see that most low-energy PUIs are in the heliotail from Figure 3(a). This is becausethe PUIs produced in the direction of the heliotail have a lower energy on average because of greaterSW slowdown by mass loading before the TS.2 Wu et al. (a) (b) (c)(d) (e)
Figure 3 . Simulated plasma-frame PUI spatial distributions at (a) v = 30 km/s, (b) 100 km/s, (c) 200km/s, (d) v = 200 km/s, and (e) 800 km/s in the meridional plane. The model is tuned to fit the observations. For example, Gloeckler & Geiss (1998) reported anaveraged phase space density measured by SWICS in a 100-day interval in 1994 when
Ulysses wasat about 3 AU near the southern pole. To compare the observed velocity distribution with modelpredictions, we transformed our model distribution functions into the frame of SWICS, and thenintegrated over the SWICS field of view. The results are presented in Figure 4. The red line is thesimulated PUI phase space density after transformation. One can see that our model result matchesthe observed phase space density well.Recent measurements by Solar Wind Around Pluto (SWAP) on board
New Horizons were reported
UIs in the Heliosphere and Heliosheath −2 Ion Speed/Solar Wind Speed P h ase S p ace D e n s i t y , f ( s k m − ) Simulated PUI phase space densityProton (Solar Wind and PUIs) phase space densityfrom Ulysses/SWICS
Figure 4 . Phase space density of H + (including SW and PUIs) versus ion speed in the spacecraft frame.Individual data points are SWICS observation in a 100-day interval in 1994. Red line are the model resultat about 3 . in Randol et al. (2013). SWAP measured SW and PUI spectra from 11 AU to 22 AU. We trans-formed our model distribution function to the spacecraft frame and converted it into count rate. Thecomparison is given in Figure 5(a). The SW H + and He peaks can be seen at energy per charge(E/q) ≈
600 eV/e and 1100 eV/e, respectively As can be seen, the slope of the spectrum and thePUI step like feature are well reproduced. At all other energies, the PUIs are hidden beneath theSW background. The spectrum is below the data at v > V SW . This is partially due to contaminationfrom SW, but also due to local acceleration of the PUIs by variations in magnetic field or the SWbackground.Figure 5(b) compares the simulated PUI intensities j ( E ) with Voyager 1
Low Energy Charged4
Wu et al.
Particle (LECP) proton intensities averaged over one selected 78-day interval in the termination fore-shock (2004/167-2004/245) and one 78-day period immediately behind the TS (2004/349-2005/061)(Decker et al. 2008). We also compare simulated PUI intensity with
Voyager 2 low-energy protonintensity averaged over a selected 78-day interval in the inner heliosheath (2007/241-2007/319) inFigure 5(c). Note that the spectral shapes in Figure 5(b) and (c) are essentially the same, well repre-sented by j ∝ E − . corresponding to f ∝ v − , with the Voyager 2 ion spectrum being slightly harderthan that at
Voyager 1 behind the TS crossing. Each of the three spectra have the step feature. Thespectra shift upward once the TS is crossed.A partial overlap exists between the
Voyager 1, 2
LECP ion instrument (28 < E < Cassini
ENA Ion and Neutral Camera (INCA) sensor (5 . < E <
55 keV) (Krimigis et al. 2010).Krimigis et al. (2010) reported the INCA-inferred ion spectrum in the heliosheath and matched it tothe in situ measured
Voyager 2 spectrum. Figure 5(d) shows our simulated PUI spectrum at ∼ Voyager 2 direction compared with the
Cassini
ENA INCA-inferred ion spectrum. DISCUSSIONOur model introduced several improvements over the existing models. Chalov et al. (2004) haveinvestigated the spatial variation of PUI spectra, but only the upwind part of the heliosheath wasconsidered. Malama et al. (2006) have performed a multi-species simulation, but the magnetic fieldwas ignored and their model was two-dimensional. The model presented here computes PUI dis-tributions on a three dimensional grid. Usmanov & Goldstein (2006) considered the SW outside 1AU as a combination of three co-moving species, SW protons, electrons, and PUIs, but they onlycomputed the global structure of the SW from the coronal base to 100 AU without the TS. Ourmodel’s external boundary is well in the LISM covering supersonic SW region, the inner heliosheathand the outer heliosheath.
UIs in the Heliosphere and Heliosheath −2 −1 Energy−per−charge, E/q (eV/e) C oun t r a t e , C ( H z ) Simulated PUI spectrumProton spectrum measured by SWAP at about 22AU (a) −2 Energy, E (keV) I n t e n s i t y , j ( E ) ( p a r t i c l es / c m − s − s r − ke V ) Spectrum of simulated PUIs upstream of the TSV1 LECP 40−53 keV channel 04/167−04/245V1 LECP 53−85 keV channel 04/167−04/245Spectrum of simulated PUIs downstream of the TSV1 LECP 40−53 keV channel 04/349−05/061V1 LECP 53−85 keV channel 04/349−05/061Slope = −1.5 (b) −1 Energy, E (keV) I n t e n s i t y , j ( E ) ( p a r t i c l es / c m − s − s r − ke V ) Spectrum of simulated PUIs downstream of the TSV2 LECP 28−43 keV channel 07/241−07/319V2 LECP 43−80 keV channel 07/241−07/319Slope = −1.5 (c) −2 −1 −1 Energy, E (keV) I n t e n s i t y , j ( E ) ( p a r t i c l es / c m − s − s r − ke V ) Simulated PUI differential fluxCassini ENA INCA−inferred ion spectrum (d)
Figure 5 . (a) Coincidence count rate of SW and interstellar PUIs versus E/q spectra observed by SWAPat 11 .
44 AU. Red circles are the model result. (b) Energy spectra of simulated PUIs and 40-85 keV ionsat
Voyager 1 , averaged over selected 78-day periods before (pre-TS) and after (post-TS) TS crossings. (c)Energy spectra of simulated PUIs and 28-80 keV ions at
Voyager 2 , averaged over selected 78-day periods inthe inner heliosheath. (d) The simulated PUI differential flux at ∼
90 AU in
Voyager 2 direction comparedwith the
Cassini
ENA INCA-inferred ion spectrum (Krimigis et al. 2010). Wu et al.
Several weaknesses of the present model are now pointed out. Firstly, the two-population assump-tion is admittedly questionable. The Schwadron et al. (1996) paper reported a range of field strengthvariations, rather than a bimodal distribution. The observed distribution functions are averagedover longer times than the averaging interval in Schwadron et al. (1996), and therefore are a productof a superposition of high and low values of η . The simulated spectra shown in Figures 1 and 5should be compared with monthly or even yearly averages of the data. Secondly, our model assumesthe PUI VDF is isotropic, which is not always the case. For example, ion angular data from Voy-ager 1 observations during 2002.58 to 2003.10, 85 . . UIs in the Heliosphere and Heliosheath
IBEX distributed ENA sky maps. This will bring uscloser to explaining why the distributed ENA flux spectrum does not show a knee and why is it closeto a power law (Schwadron et al. 2011).This work was supported by NASA grant NNX12AH44G.REFERENCES