STABILITY OF A PICKUP ION RING-BEAM POPULATION IN THE OUTER HELIOSHEATH: IMPLICATIONS FOR THE IBEX RIBBON

The Interplanetary Boundary Explorer (IBEX) satellite, launched in 2008 October, is currently observing the structure of the outer heliosphere in energetic neutral atoms (ENAs). The first year of operation brought forth a number of results that were not predicted by theoretical models (McComas et al. 2009; Fuselier et al. 2009; Funsten et al. 2009; Schwadron et al. 2009). The most striking discovery was a three-fold ENA enhancement measured at ∼1 keV above the surrounding emission (smaller at higher and lower energies) from a region in the sky in the shape of a nearly circular "ribbon." The largest intensity enhancement was measured at ∼1 keV (Fuselier et al. 2009).

Recently, Heerikhuisen et al. (2010) argued that the enhancement could be due to secondary ENAs produced in the outer heliosheath (OHS) from a parent population of heliosheath and solar wind pickup ions (PUIs). In their model, ENA emission comes primarily from those regions where the local interstellar (LISM) magnetic field is nearly perpendicular to the line of sight, i.e., ${\mathbf B}\,{\bm{\cdot}}\, \hat{r}\simeq 0$ (see also Schwadron et al. 2009). The model-predicted ribbon is consistent with both the ENA observations and the magnetic field orientation based on studies of the field draping around the heliopause (Pogorelov et al. 2009). The peak in intensity enhancement at 1 keV is explained by the fact that OHS PUIs produced from solar wind ENAs would have approximately that energy (1 keV corresponds to a proton speed of 438 km s⁻¹).

One of the principal assumptions of the OHS PUI model of the ribbon is that the pitch-angle evolution of the ring ion distribution is slow enough to allow these particles to charge exchange (thus becoming secondary ENAs) before scattering completely onto a spherical shell. Heerikhuisen et al. (2010) used a "partial shell" distribution that is in between a narrow ring and an isotropic shell. The width of the ring was dependent on the angle between the parent ion velocity vector and LISM B implying, essentially, that scattering does not extend the velocity distribution to smaller pitch angles.

The time scale for charge exchange, τ_ce = (n_Hσv)⁻¹, where n_H ≃ 0.2 cm⁻³ is the density of neutral H in the OHS, σ ≃ 2 × 10⁻¹⁵ cm² is the charge-exchange cross section at 1 keV (Lindsay & Stebbings 2005), and v ∼ 400 km s⁻¹ is the velocity of fast PUIs, is on the order of two years. This time is long compared with the electromagnetic time scales in the system that govern wave generation and pitch-angle scattering. The latter is typically some function of a gyro-period τ_g = 2π/Ω_i (Ω_i is the proton gyrofrequency) which is on the order of 4 minutes in the OHS. Thus, there is a factor of ∼10⁵ disparity between temporal macroscales and microscales of the problem.

If waves are present in a plasma, a ring distribution will undergo pitch-angle scattering to become a shell (Williams & Zank 1994). Scattering will occur regardless of the origin of the waves, which could be part of ambient interstellar turbulence, generated near the heliopause, e.g., as a result of a large-scale instability (Axford et al. 1963; Fahr et al. 1986; Baranov et al. 1992), or excited by the PUI ring itself, which is known to be unstable to wave generation (Wu & Davidson 1972; Sharma & Patel 1986; Gary & Madland 1988). This paper deals primarily with PUI-excited waves. Below we investigate the process of wave excitation by a gyrating ion ring distribution in a warm plasma of the OHS using a linear Vlasov theory of electromagnetic fluctuations (Section 3) as well as nonlinear hybrid simulations of parallel (also called slab) modes in a plasma with thin and broadened ring-beams (Section 4). The latter approach allows us to examine self-consistently the scattering of PUIs on self-excited waves which leads, eventually, to their isotropization. We discuss the kinds of wave modes generated, the resulting spectrum of the waves, and the timescales for pitch-angle scattering. Finally, we compare the wave intensity generated as a result of the instability with possible ambient turbulence levels in the LISM (Section 5).

Using the published IBEX ENA intensities from the direction of the ribbon (Funsten et al. 2009; Heerikhuisen et al. 2010) we can estimate the parent population of PUIs produced by charge exchange of neutral solar wind atoms, which are presumably responsible for ENA emission at ∼1 keV. To do this, we note that in the limit of collisionless transport the ENA differential intensity at IBEX is an integral over the line of sight l,

$$\begin{equation} J_{\rm ENA}=\int _{l_1}^{l_2}J_{\rm PUI}\sigma n_H dl\simeq J_{\rm PUI}\sigma n_H \Delta l, \end{equation} \tag{ 1 }$$

where Δl is the width of the emitting region. According to MHD models (e.g., Pogorelov et al. 2009) the width of the OHS is on the order of 200 AU, which can be used for Δl above. The PUI differential intensity is f_PUIp², which for a ring distribution of thickness Δp in momentum and pitch-angle width ΔΩ may be approximated as

$$\begin{equation} f_{\rm PUI}\sim \frac{n_b}{\Delta \Omega p^2\Delta p}, \end{equation} \tag{ 2 }$$

where n_b is the number density of the PUI ring beam. With ΔΩ ∼ 2π, Δp = ΔE/v, where ΔE ∼ 0.5 keV (Figure 4 of Heerikhuisen et al. 2010), and J_ENA = 200 (cm² s sr keV)⁻¹ at 1 keV we obtain

$$\begin{equation} n_b\simeq \frac{2\pi J_{\rm ENA}\Delta E}{\sigma n_H \Delta l v} \simeq 1.5\times 10^{-5}\;\mathrm{cm}^{-3}. \end{equation} \tag{ 3 }$$

Since the ambient plasma density in the OHS is about 0.08 cm⁻³ (Pogorelov et al. 2009), for all calculations in this paper we adopt a relative density of ring to background plasma ions of n_b/n_i = 1.8 × 10⁻⁴.

A warm Maxwellian electron–proton plasma supports three types of low-frequency waves propagating parallel to the LISM magnetic field B: an ion-acoustic longitudinal mode and two transverse electromagnetic waves (Alfvén-ion cyclotron and magnetosonic whistler). Relevant parameters are the plasma frequencies of the ions, assumed to be protons, and of the electrons, ω_p, ω_i, their cyclotron frequencies Ω_p, Ω_i, and their thermal speeds v_ti, v_te. For particles of species α these are given by

$$\begin{equation} \omega _{\alpha }=\sqrt{\frac{4\pi n_{\alpha }e_{\alpha }^2}{m_{\alpha }}},\;\;\; \Omega _{\alpha }=\frac{e_{\alpha }B}{m_{\alpha }c},\;\;\; v_{t\alpha }=\sqrt{\frac{2kT_{\alpha }}{m_{\alpha }}}, \end{equation} \tag{ 4 }$$

where m_α is the mass, e_α is the signed electric charge, T_α is the temperature, and n_α is the density of species α. Here we assume n_i ≃ n_e = 0.1 cm⁻³ and B = 2.5 μG, which are representative of the OHS region. With these parameters the Alfvén speed $V_a=B/\sqrt{4\pi m_i n_i}=17.2$ km s⁻¹ and the ion cyclotron frequency Ω_i = 0.024 s⁻¹.

The newly born PUIs represent an additional ring-beam component that is both drifting along B and gyrating in the perpendicular direction. We assume that the neutral solar wind hydrogen atoms are moving in a straight line radially from the Sun with a constant speed v = 400 km s⁻¹. We work in the standard plasma frame with the z-axis directed along B. Upon ionization in a charge transfer collision with an interstellar proton the newly born PUI's velocity in this frame is v_b = (v_⊥bcos Ω_bt, v_⊥bsin Ω_bt, v_∥b), where v_⊥b = vsin α, v_∥b = vcos α, with α being the pickup angle (the angle between the direction of travel of the parent hydrogen atom and the mean magnetic field). Here we designate parameters pertaining to the PUI ring beam with the letter "b." The distribution function of the PUIs in the plasma frame is

$$\begin{equation} f_b=\frac{n_b}{2\pi v_{\perp b}} \delta (v_{\parallel }-v_{\parallel b}) \delta (v_{\perp }-v_{\perp b}). \end{equation} \tag{ 5 }$$

The following analysis follows essentially that of Gary & Madland (1988), but using much lower ring-beam densities and background parameters appropriate for the OHS. We only consider parallel-propagating (slab) transverse electromagnetic waves because longitudinal electrostatic waves do not scatter ions in a pitch angle. For such modes the dispersion relation may be written as

$$\begin{equation} \left(\begin{array}{@{}cc@{}} C & iD \\ -iD & C \\ \end{array}\right) \left(\begin{array}{@{}c@{}} E_x \\ E_y \end{array}\right)=0, \end{equation} \tag{ 6 }$$

where E is the electric field and the coefficients C and D are given by

$$\begin{eqnarray} C &=& k^2c^2-\omega ^2 -\frac{\omega _i^2\omega }{2kv_{ti}} \left[Z\left(\frac{\omega -\Omega _i}{kv_{ti}}\right) +Z\left(\frac{\omega +\Omega _i}{kv_{ti}}\right)\right]\nonumber \\ &&-\frac{\omega _e^2\omega }{2kv_{te}} \left[Z\left(\frac{\omega -\Omega _e}{kv_{te}}\right) +Z\left(\frac{\omega +\Omega _e}{kv_{te}}\right)\right] \nonumber\\ &&+ \frac{\omega _b^2}{2} \left[\frac{\omega -kv_{\parallel b}}{\omega -kv_{\parallel b}-\Omega _b} +\frac{\omega -kv_{\parallel b}}{\omega -kv_{\parallel b}+\Omega _b}\right] \nonumber\\ &&+\frac{\omega _b^2}{4} \left[\frac{k^2v_{\perp b}^2}{(\omega -kv_{\parallel b}-\Omega _b)^2} +\frac{k^2v_{\perp b}^2}{(\omega -kv_{\parallel b}+\Omega _b)^2}\right], \end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray} D &=& -\!\frac{\omega _i^2\omega }{2kv_{ti}} \left[Z\left(\frac{\omega -\Omega _i}{kv_{ti}}\right) -Z\left(\frac{\omega +\Omega _i}{kv_{ti}}\right)\right]\nonumber\\ &&-\frac{\omega _e^2\omega }{2kv_{te}} \left[Z\left(\frac{\omega -\Omega _e}{kv_{te}}\right) -Z\left(\frac{\omega +\Omega _e}{kv_{te}}\right)\right] \nonumber \\ &&+\frac{\omega _b^2}{2} \left[\frac{\omega -kv_{\parallel b}}{\omega -kv_{\parallel b}-\Omega _b} -\frac{\omega -kv_{\parallel b}}{\omega -kv_{\parallel b}+\Omega _b}\right]\nonumber\\ &&+\frac{\omega _b^2}{4} \left[\frac{k^2v_{\perp b}^2}{(\omega -kv_{\parallel b}-\Omega _b)^2} -\frac{k^2v_{\perp b}^2}{(\omega -kv_{\parallel b}+\Omega _b)^2}\right], \end{eqnarray} \tag{ 8 }$$

where k is the parallel wavenumber, ω = ω_r + iγ is the complex frequency, and Z is the plasma dispersion function. The normal modes of Equation (6) are given by C ± D = 0. Numerically found unstable roots are plotted in Figures 1 and 2 for T_i = T_e = 0 and T_i = T_e = 20, 000 K, respectively, the latter value is appropriate for the OHS plasma. The cold plasma case is included because it conveniently illustrates the physics of the instabilities and their coupling to the background plasma modes. We use three quasi-perpendicular pickup angles α = 80°, 85°, and 90°. Note that k = Ω_i/V_a corresponds to a wavelength of ∼5 × 10⁸ cm.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

With very low ring-beam densities (n_b = 1.8 × 10⁻⁴ n_i) instabilities arise only near resonance frequencies given by the zeroes of the denominators of the last two terms in Equations (7) and (8), i.e.,

$$\begin{equation} \omega _r=kv_{\parallel b}\pm \Omega _b. \end{equation} \tag{ 9 }$$

There are several plasma beam instabilities known to operate in the low (near ion cyclotron) frequency regime (Gary et al. 1984; Winske et al. 1985). Of interest to us is the cyclotron resonant mode (Montgomery et al. 1976; Winske et al. 1984). The (mostly) right-hand part of this instability, labeled "R" in Figures 1(a) and 1(b), appears on the magnetosonic branch coupling it to the Alfvén wave, labeled "A" in that figure. This mode, corresponding to the minus sign in the resonance condition (9), can be either right-hand or left-hand polarized, with ω_r > 0 and ω_r < 0, respectively. The same dispersion branch features another instability at much higher (whistler) frequencies (Wu & Hartle 1974), which is of no interest here. The second instability corresponds to the plus in Equation (9) and is shown in the same figures as the "L" mode (Gary & Madland 1988). This left-hand instability generally has a smaller growth rate and is expected to be less important.

The "R" mode has regions of stability that move to higher wavenumbers with increased α, while the left-hand part is unstable for the entire range of pickup angles considered here. The low wavenumber stability region becomes an oppositely propagating magnetosonic wave (not shown). In the singular case α = 90°, the left- and the right-hand unstable branches merge together. From Figure 1(c) it can be seen that this instability modifies the parallel Alfvén dispersion branch, while the magnetosonic wave is not at all affected by the presence of the gyrating ring (Sharma & Patel 1986).

Dispersion curves in a warm plasma with a PUI beam are similar (Figure 2), except that the wave speeds are modified somewhat from the cold plasma case. T_i = 20, 000 K corresponds to a plasma β (the ratio of ion thermal to magnetic pressures) of 1.1. The Alfvén wave is strongly damped beyond Ω_i/V_a. At pickup angles close to 90° the right-hand mode becomes uncoupled from the Alfvén branch as seen in Figure 2(c). Apart from that the ion ring-beam instabilities in a warm heliosheath plasma have nearly identical growth rates owing to the resonant nature of the interaction. Typical growth rates are on the order of 0.2Ω_i or 0.005 s⁻¹.

The linear analysis predicts that waves will be excited, but reveals nothing about their behavior past the brief exponential growth phase. We certainly do not expect this phase to persist for more than a few gyroperiods because scattering will cause particles to spread in the pitch angle, leaving fewer and fewer to resonantly interact with each particular mode. While initially scattering will occur, we do not know yet whether it will lead to a complete isotropization of the ring, and if so, on what timescale. To proceed, a nonlinear model taking into account ion scattering must be used.

To investigate the excitation of slab (k||B) modes by the instability and track its nonlinear development and the scattering of freshly ionized PUI rings we use a one-dimensional (1D) hybrid electromagnetic code that treats ions as individual particles and electrons as a charge-neutralizing adiabatic fluid (e.g., Winske et al. 1984; Terasawa et al. 1986). The magnetic field is in the x-direction and the plasma is at rest initially (note the difference with Section 3, where B was in the z-direction). Grid resolution is uniform and equal to half the ion inertial length (0.5 c/ω_i), which is sufficient to resolve the low-frequency electromagnetic modes. The simulation box length L_x is taken to be large enough to accommodate waves resonant with particles streaming along the field, i.e., L_x ≫ 2πv_∥b/Ω_i. A box size of ∼2000 c/ω_i is sufficient to include the longest waves in the problem, excited by the particles with a zero pitch angle. We use equal numbers of particles (typically 25,000 per cell) for the background protons and for the ring beam. The latter are given a weight of n_b/n_e = 1.8 × 10⁻⁴, whereas the weight of the background ions is adjusted to 1 − n_b/n_e to maintain charge neutrality. These simulations have excellent statistics (∼2 × 10⁸ particles) ensuring that the weak PUI signal is easily identified against the background noise. Relevant dimensionless parameters are β = 1.1 and ω_i/Ω_i = c/V_a = 1.74 × 10⁴. Periodic conditions are imposed at the spatial boundaries.

The first simulation starts with a narrow ring PUI distribution given by Equation (5) injected with v_b = 450 km s⁻¹ and α = 85°. Waves are excited from the thermal noise by the instability. These waves are initially counter-streaming with the beam and right-hand polarized, corresponding to the curve labeled "R" in Figure 2(b). Eventually, however, waves of both polarizations appear. Figure 3 shows several snapshots from the simulation at times t = 30/Ω⁻¹_i (1), t = 300/Ω⁻¹_i (2), and t = 4300/Ω⁻¹_i (3). The last panel corresponds to a real time of ∼2 days. The top panel for each time is the velocity distribution of the core (black) and the PUI ring (red), while the bottom panel shows the 1D power spectrum of the two transverse components of the magnetic field

$$\begin{eqnarray} P_{yy}(k)+P_{zz}(k) &=& 2\int _0^{\infty }[ \langle B_y(x) B_y(x+\xi)\rangle\nonumber\\ &&+\langle B_z(x) B_z(x+\xi)\rangle ]\cos (k\xi)d\xi, \end{eqnarray} \tag{ 10 }$$

where 〈 ⋅ ⋅ ⋅ 〉 denotes averaging with respect to x. In practice, a discrete Fourier transform of the two-point correlation function is taken (see Bieber et al. 1995). The spectrum above k = ω_i/c is not part of the PUI signal, but is essentially white noise of numerical origin.

Zoom In Zoom Out Reset image size

The narrow ring begins to broaden shortly after the start of the simulation. The waves initially form a narrowband spectrum centered at the resonant wavenumber k = (0.3–0.4) ω_i/c (Figure 3(a)). As the simulation progresses particles scatter to smaller pitch angles filling the forward part of the shell and exciting progressively longer waves in the process. The spectral maximum shifts to k = (0.1–0.2) ω_i/c (Figure 3(b)). Scattering across 90° is inhibited in this case in agreement with quasi-linear theory (e.g., Jokupii 1971). Some scattering is possible across the "resonance gap" due to a finite v/V_a ≃ 0.038 (Schlickeiser 1989). After one day the PUI population is a nearly isotropic shell (Figure 3(c)). The wave spectrum is similar to t = 300 Ω⁻¹_i and extends down to k ≃ 0.03 ω_i/c. A typical wave amplitude δB/B ∼ 0.05 near the end of the simulation. We note that the time scale for isotropization is almost 3 orders of magnitude larger than the characteristic timescale of the exponential wave growth phase (∼200 s).

Our second study case involves a broadened gyrotropic distribution of the form

$$\begin{equation} f_b=\frac{2\delta (v_{\parallel }-v_{\parallel b}) e^{-(v_{\perp }-v_{\perp b})^2/v_{\perp t}^2}}{v_{\perp t}^2\left\lbrace e^{-v_{\perp b}^2/v_{\perp t}^2} +\sqrt{\pi }v_{\perp b}/v_{\perp t}[1+\mathrm{erf}(v_{\perp b}/v_{\perp t})] \right\rbrace } \end{equation} \tag{ 11 }$$

(e.g., McClements & Dendy 1993), where v_⊥b is the velocity spread of the PUIs in the OHS. Broadening in perpendicular velocity occurs because some of the parent ENAs were born in a fast SW, while others were produced in a slow wind or in the heliosheath. The distribution (11) is therefore a more realistic representation of the actual OHS PUI population.

In this simulation, we have used the same v_b and α as in our narrow ring simulation, while the "thermal" spread of the distribution was half the mean speed v_⊥t = 225 km s⁻¹. Qualitatively, the results of this experiment, shown in Figure 4, are not very different from those obtained with a narrow ring. Both distributions have a similar number of wave-resonant particles (the resonance condition only depends on the parallel speed), which explains the similarity in the wave spectra in these two simulations. Note that in both cases the background ions, shown as black dots in Figures 3 and 4, were not affected by the low-density PUI ring.

Zoom In Zoom Out Reset image size

The third simulation (Figure 5) shows a simulation for a widened distribution with an initial spread in the pitch angle. This is essentially the partial shell model of Heerikhuisen et al. (2010) and is given by

$$\begin{equation} f_b=\frac{n_b}{4\pi \cos \theta _0 v_b^2}\delta (v-v_b),\;\;\; \theta _0<\theta <\pi -\theta _0, \end{equation} \tag{ 12 }$$

where θ is the pitch angle (cf. Equation (2)). A distribution of the form Equation (12) implies a range of pickup angles owing to (1) a difference in the direction of travel of parent heliosheath ENAs and (2) fluctuations in the direction of the magnetic field in the LISM. Here we used θ₀ = 60° (cf. Section 2). While the initial distribution (12) resembles a partial shell attained after a scattering time of 3.5 hr in our earlier simulation (Figure 3(b)), this simulation is quite different because the partial shell is introduced without the associated turbulence. This leads to a scarcity of short waves in the spectrum.

Zoom In Zoom Out Reset image size

In this scenario, short waves are not excited because the resonant particles are on a locally isotropic (stable) part of the shell. Only particles at the edge of the distribution (i.e., near θ₀) generate the waves which explain the narrow spectrum in Figure 5(b). After some 35 hr (Figure 5(c)) the distribution has nearly relaxed to a spherical shell. Relaxation is faster in this case because there is no need for particles to scatter across 90°, as both hemispheres were equally populated by the PUIs from the start.

Selected pitch-angle distributions (PADs) of 1 keV PUI protons derived from the broad shell simulation at t = 3.5 hr and t ≃ 2 days are shown in Figure 6. PADs are constructed as phase angle averaged distributions per interval in velocity per interval in solid angle, namely, −ΔN/(2πv²ΔvΔμ), where μ is the pitch-angle cosine. At t = 3.5 hr (dashed line) the distribution is a partially filled shell with most particles traveling in the forward direction. The state at t = 50 hr consists essentially of two half shells (forward and backward) with little pitch-angle variation in each half shell. Despite some difficulty in scattering across μ = 0, particles do eventually migrate to the μ < 0 space with a final distribution invariably being an isotropic shell.

Zoom In Zoom Out Reset image size

The interstellar medium is known to be in a turbulent state from remote sensing measurements of electron density. The spectrum of the turbulence appears to be a power law with the Kolmogorov index −11/3, corresponding to a reduced 1D spectrum ∼k^−5/3, spanning some 10 orders of magnitude in wavenumber from parsec scales down to 10⁻⁸ cm⁻¹ (Armstrong et al. 1995). The power-law inertial range of the turbulence transitions to the energy containing range at the outer scale, where the energy input from large-scale motion of the interstellar gas occurs. Recent estimates place the outer scale at a few parsec, i.e., a few times 10¹⁸ cm (Haverkorn et al. 2008).

It is of interest to compare the ambient intensity of turbulent fluctuations at PUI resonant wavenumbers with the in situ turbulence generated by the ring instability. To perform this comparison we assume that the reduced (slab) LISM magnetic power spectrum is composed of a flat energy containing range followed by a Kolmogorov power law, i.e.,

$$\begin{equation} P_{{\rm sl}}(k)=\frac{\pi }{5}\delta B_{\rm {sl}}^2 l_c\left\lbrace \begin{array}{@{}ll} 1, & kl_c<1, \\ \ms (kl_c)^{-5/3}, & kl_c>1, \end{array}\right. \end{equation} \tag{ 13 }$$

where l_c is the bendover length, which is related to the outer scale of the turbulence, and δB²_sl is the mean square fluctuating field in the slab component. The spectrum (13) is normalized such that (cf. Equation (10))

$$\begin{equation} \delta B_{\rm {sl}}^2=\langle B_y^2(x)\rangle +\langle B_z^2(x)\rangle =\frac{1}{\pi }\int _0^{\infty }P_{\rm {sl}}(k)dk. \end{equation} \tag{ 14 }$$

We can also make a simple, but plausible, assumption that the turbulent magnetic field is a fixed fraction of the mean, i.e., that δB²/B² ∼ 0.1–1. The slab component will be on average about a third of the total, so we may well take δB²_sl/B² ∼ 0.03–0.3. In addition to interstellar turbulence waves may be excited due to large-scale motions in the OHS. Some possible drivers include the solar wind large-scale transient structures and the solar-cycle dynamic pressure variations (Zank & Müller 2003; Scherer & Fahr 2003; Izmodenov et al. 2005), and the instability of the heliopause due to a difference in charge-exchange rates across the boundary (Liewer et al. 1996; Zank et al. 1996; Florinski et al. 2005; Borovikov et al. 2008). The scale of these motions is expected to be significantly shorter than the interstellar turbulent scales, perhaps on the order of 50 AU.

Figure 7 compares the LISM spectra in the slab component of the turbulence computed using B = 2.5 μG, l_c = 10¹⁸ cm (blue line), the possible locally driven large-scale turbulence with l_c = 7.5 × 10¹⁴ cm (green), and the wave intensity taken from the hybrid simulations in Figure 3(c) (red line). In deriving both ambient spectra we took δB²_sl/B² = 0.3, an extreme scenario which assumes LISM turbulence is very strong. Even so, the disparity in spatial scales combined with a steep power-law spectrum implies that the amount of fluctuations available to scatter PUIs is some 3–4 orders of magnitude below that of the self-excited waves. This means that ambient waves are not likely to contribute to the PUI scattering process. This conclusion remains in force for locally produced turbulence that is a result of a plasma motion on macroscopic scales, such as the Rayleigh–Taylor instability of the heliopause.

Zoom In Zoom Out Reset image size

Independently using linear Vlasov theory and numerical hybrid simulations we have shown that a ring of pickup protons introduced into the OHS plasma by charge exchange with neutral solar wind hydrogen is unstable to wave generation by a resonant cyclotron mechanism. Hybrid simulations revealed that the generated waves in turn resonantly scatter the PUIs onto a shell. The time scale for isotropization is on the order of several days. This time, while significantly in excess of a proton gyroscale, is in turn substantially shorter than the time scale for charge exchange in the OHS, which is about two years. This appears to contradict the assumption of Heerikhuisen et al. (2010) that the PUI distribution remains narrow for an extended period of time. An anisotropic PUI distribution in the ribbon (i.e., the regions where ${\mathbf B}\,{\bm{\cdot}}\, \hat{r}\simeq 0$) is required for their mechanism to work. While some isotropy will produce a ribbon, the flux levels will be considerably lower than what IBEX observes.

Waves generated by the PUI instability could contribute substantially to turbulence in the OHS region. The waves would persist even after the free energy in the ring has been exhausted because the damping rate of the long wave right-hand mode (i.e., the magnetosonic wave) is essentially zero in the OHS. The waves could eventually develop into a turbulent spectrum with implications for low-energy cosmic-ray scattering in the OHS (e.g., Florinski et al. 2003) and heating of the LISM plasma.

The simulations presented here were performed without taking into account spatial transport of PUIs. Particles with small pitch angles can potentially stream away from their place of origin along the interstellar magnetic field lines. Escape along field lines may alter the PUI distribution in the ENA emitting region (i.e., the ribbon). The emitting region cannot be very narrow, however, given that the width of the ribbon is about 20° (Fuselier et al. 2009), which corresponds to a width of some 100 AU in the OHS. It would take close to a year for 1 keV PUIs to travel that distance, which is again very long compared with the scale for the isotropization. We must conclude that spatial transport will not change our result appreciably except perhaps very close to the edge of the ribbon.

Our simulations were limited by the geometry of the wave field (slab only). An oblique multi-dimensional geometry could alter wave growth rates compared with the parallel case. This scenario is currently under investigation. The slab model is adequate for studying scattering in the pitch angle, but is not very suitable for spatial transport (especially, transport across the mean magnetic field) and turbulence development. Because wave growth is usually largest in the parallel direction, the slab model should provide a reasonable upper limit estimate on the wave generation process.

Our theoretical results are not the final verdict on the PUI ring stability in the OHS. On several occasions PUI rings injected at nearly perpendicular angles were observed to be stable, in seeming contradiction to the theoretical expectations. For example, Tsurutani et al. (1989) reported an absence of cyclotron waves near the comet Giacobini–Zinner when the pickup angle was between 70° and 90°. They also reported a lack of pitch-angle scattering and persistently high anisotropies of cometary PUIs. Gloeckler et al. (1995) showed that interstellar pickup H distributions measured at high heliographic latitudes (presumably for magnetic field directions close to radial) in the solar wind retained a relatively high anisotropy, which they interpreted as being due to a lack of instability-generated waves (see also Möbius et al. 1998). Finally, Oka et al. (2002) showed that He PUIs in some circumstances would have a torus-like distribution for a broad range of pickup angles. It appears that some physical effects that could delay isotropization are missing from existing theoretical models.

Turbulence cascading may be the missing piece of the puzzle. Our 1D simulations do not adequately describe this process, which is an essentially 3D effect, involving a transfer of power between oblique modes. The characteristic nonlinear time in the OHS is τ_nl ∼ 1/(kV_a) ∼ 500 s, which is shorter than the isotropization time based on the results of Section 4. It is therefore possible that spectral transfer would channel the fluctuation power away from the resonant regime leaving fewer waves to scatter the PUIs. We will address this topic in the future using multi-dimensional electromagnetic plasma simulations (Khazanov & Singh 2008).

This research was supported, in part, by NASA grants NNX10AE46G and NNX09AG63G and by NSF grant AGS-0955700.