PICKUP ION MEDIATED PLASMAS. I. BASIC MODEL AND LINEAR WAVES IN THE SOLAR WIND AND LOCAL INTERSTELLAR MEDIUM

Most of the heliosphere by volume is partially ionized with the number density of neutral hydrogen (H) greater beyond ∼10 AU than that of the solar wind plasma (Zank 1999; Zank et al. 2013; Heerikhuisen et al. 2014). The collisional charge exchange mean free path in the outer heliosphere ensures that we do not have to directly couple neutral H and the solar wind plasma. However, neutral H in the outer heliosphere is responsible for the creation of pickup ions (PUIs) that eventually form a thermodynamically dominant component of the plasma despite being relatively tenuous. When PUIs are present, they typically play a critical role in almost all facets of the underlying plasma physics.

Several efforts have been made to incorporate PUIs in a self-consistent manner into models that describe the solar wind. A number of models, summarized in Zank (1999), from the 1970s (see especially Holzer 1972; Wallis 1971 for example) assumed a single component background plasma into which pickup protons were added via source terms in the continuity (assuming photoionization), momentum, and energy equations. A more formal derivation of the one-fluid model for the supersonic solar wind mediated by PUIs was provided by Khabibrakhmanov et al. (1996) based on a guiding center kinetic equation derivation. An important extension to the one-fluid models of the supersonic solar wind was introduced by Isenberg (1986). The one-fluid models assume essentially that wave–particle interactions (and particle collisions) proceed sufficiently quickly that PUIs are almost immediately assimilated into the thermal solar wind, becoming indistinguishable from thermal solar wind protons. The one-fluid model of the solar wind predicts a substantial increase in the plasma temperature or pressure (see, e.g., Figure 4.2 of Zank 1999) with increasing heliocentric distance. Such an increase in the temperature of solar wind protons is not observed in the outer heliosphere (Gazis et al. 1994; Richardson et al. 1995) since the Voyager plasma instrument cannot measure PUIs, but a modest increase in the thermal proton plasma temperature is found beyond about 20 AU. The modest temperature increase is due to the dissipation of PUI-generated turbulence (Williams et al. 1995; Matthaeus et al. 1999; Smith et al. 2001; Isenberg et al. 2003; Isenberg 2005; Smith et al. 2006, see also the three-dimensional (3D) models by Kryukov et al. 2012; Usmanov et al. 2014). Instead, as noted by Vasyliunas & Siscoe (1976) and Isenberg (1986), PUIs are unlikely to be assimilated into the supersonic solar wind. Pre-existing and PUI-excited waves will isotropize and stabilize the new-born PUI distribution to produce two approximately co-moving proton populations.

Coulomb collisions are necessary to equilibrate a background thermal plasma, such as the solar wind, and the PUI protons. Let us consider the general case of a background Maxwellian plasma comprised of thermal protons (denoted by the subscript s) and electrons (denoted by the subscript e). PUIs (denoted by the subscript p) then satisfy the ordering

$$\begin{equation*} v_{ts} \ll v_{p} < v_{te}, \end{equation*}$$

where v_ts/e denotes the background proton/electron thermal speed respectively and v_p the PUI speed. Denoting quantities derived from the background Maxwellian proton or electron distribution with the subscript "b = s/e," the collisional frequency for PUIs is given by (e.g., Zank 2013)

$$\begin{equation} \nu _s^{pb} = 2 \bar{\nu }_{pb} \frac{T_p}{T_b} \left(1 + \frac{m_b}{m_p} \right) \frac{G(x_b)}{x_p}, \quad \bar{\nu }_{pb} \equiv \frac{n_b e^4 \ln \Lambda }{4\pi \varepsilon _0^2 m_p^2 v_{tp}^3 }, \end{equation} \tag{ 1 }$$

where x_b ≡ v/v_tb, x_p ≡ v/v_tp, and $v_{ta} = \sqrt{2kT_a /m_a}$ is the thermal speed, T_p/b is the PUI/background electron or proton temperature, n_b the background electron or proton number density, ε₀ the permittivity of free space, e the electron charge, ln Λ the Coulomb logarithm (∼20), and m_p/b the PUI/proton or electron mass. The function G(x) is the Chandrasekhar function for which

$$\begin{equation*} G(x) \sim \left\lbrace \begin{array}{@{}cc}\sim \frac{2x}{3 \sqrt{\pi } } & {{\rm for} x \ll 1} \\ \sim \frac{1}{2x^2} & {{\rm for} x \gg 1} \end{array} \right.. \end{equation*}$$

For PUIs scattering collisionally off a Maxwellian distribution of background protons, the collision frequency becomes (using the large argument approximation of the Chandrasekhar function)

$$\begin{equation} \nu _s^{ps} = \frac{n_s e^4 \ln \Lambda }{2\pi \varepsilon _0^2 m_p^2 v^3} s^{-1}, \end{equation} \tag{ 2 }$$

illustrating the well-known v⁻³ dependence with PUI speed. By contrast, PUI collisional scattering off a Maxwellian electron background yields a larger collision frequency, which is given by

$$\begin{equation} \nu _s^{pe} = \frac{n_e e^4 \ln \Lambda m_e^{1/2} }{2 (2\pi)^{3/2} \varepsilon _0^2 (k T_e)^{3/2} m_p } \end{equation} \tag{ 3 }$$

after using the small argument approximation to the Chandrasekhar function. If the collisional timescale exceeds the characteristic flow time of the plasma region of interest, τ_f ≃ L/U, where L is the size of the region and U the characteristic velocity, then the PUI distribution will not equilibrate with the background plasma. Expressions (2) and (3) are necessary to determine whether one needs to introduce a model that distinguishes the PUIs from the background plasma protons. In the case of the supersonic solar wind, Isenberg (1986) used related arguments to justify the introduction of a form of multi-fluid model to describe a coupled solar wind–PUI plasma. Isenberg (1986) assumed that PUIs and the solar wind plasma are perfectly co-moving, i.e., U_p = U_s, where U_p is the PUI bulk velocity and U_s the thermal solar wind proton bulk velocity. This assumption corresponds to assuming instantaneous isotropization of the PUI distribution at the time of creation. While a reasonable approximation when considering the effects of PUIs on a large-scale flow such as the supersonic solar wind, to understand the impact of PUIs on the basic plasma physical system requires that the PUI-scattering timescale not be neglected. As we discuss below, in the distant outer heliosphere, PUI scattering yields a dissipative collisionless heat flux and viscosity in a collisionless coupled thermal plasma–PUI system. In the inner heliosheath (IHS) regions, the electron heat flux is expected to dominate.

Here, we consider three specific environments in which PUIs mediate the plasma properties: the supersonic solar wind in the outer heliosphere beyond ∼10 AU, the subsonic solar wind, i.e., the plasma in the IHS, and the plasma in the very local interstellar medium (VLISM). Each of these regions is mediated by PUIs, but the origin of the PUI population in each region is different in important ways. We discuss each in turn.

Supersonic solar wind beyond 10 AU. Interstellar neutral gas flows relatively unimpeded into the heliosphere, possibly experiencing some "filtration" at the heliospheric boundaries. Neutral interstellar hydrogen is especially susceptible to the effects of filtration, being decelerated and heated in passing from the VLISM into the heliosphere. The velocity of the neutral H drifting into the supersonic solar wind across the heliospheric termination shock is ∼23 km s⁻¹ or less after filtration in the heliospheric boundaries (Zank et al. 1996b, 2013). The interstellar neutral gas flowing into the supersonic solar wind can be ionized by either solar photons (photoionization) or solar particles (charge exchange, electron-impact ionization, etc.) and the new ions are accelerated almost instantaneously by the motional electric field of the solar wind, E = −U × B, where U is the solar wind flow velocity and B the ambient interplanetary magnetic field (IMF). In a Cartesian frame co-moving with the solar wind, the velocity of a pickup ion is simply v(t) = (− U_⊥cos Ω_it, U_⊥sin Ω_it, U_∥), where the IMF is oriented along $\hat{\bf z}$, U_∥ is parallel to $\hat{\bf z}$, U_⊥ is perpendicular to $\hat{\bf z}$, and Ω_i ≡ qB/m is the local pickup ion gyrofrequency (q denotes charge and m ion mass). The pickup ions therefore form a ring-beam distribution on the timescale $\Omega _i^{-1}$ that streams sunward along the magnetic field. Newly created pickup ions can therefore drive a host of plasma instabilities (e.g., Lee & Ip 1987; see Zank 1999 for a summary). Since the variation in both U and B in the outer solar wind can be substantial, the ring beam should be broad, although this is not assumed typically when investigating related instabilities. PUIs experience scattering and gradual isotropization by either ambient or self-generated low-frequency electromagnetic fluctuations in the solar wind plasma. Since the newly born ions are eventually isotropized, their bulk velocity is now essentially that of the solar wind, i.e., they advect with the solar wind flow, and are then said to be "picked up" by the supersonic solar wind. The isotropized pickup ions form a distinct suprathermal population of energetic ions (∼1 keV energies) in the solar wind whose origin is the interstellar medium (ISM; e.g., Holzer 1972; Lee & Ip 1987; Williams & Zank 1994; see Zank 1999 for an extensive review).

At 10 AU, the number density of PUIs is about 5% of the background thermal solar wind ions, with the fraction increasing to perhaps as much as 20% in the vicinity of the termination shock at ∼90 AU. For the supersonic solar wind at 10 AU, we take the background proton and electron number density to be approximately n_s/e ∼ 0.08 cm⁻¹ and the PUI speed to be ∼400 km s⁻¹. This yields a PUI–proton collision time $\tau _s^{ps} = (\nu _s^{ps})^{-1} \sim 8.34 \times 10^{10}$ s. Although the proton number density decreases as r⁻² with increasing heliocentric radius r, we may nevertheless compare this collisional time with the convection time of a solar wind parcel from 10 AU to 90 AU, which is τ_f ∼ 3.38 × 10⁷ s. Clearly, in the supersonic solar wind PUI–proton collisions cannot equilibrate the two proton populations. For PUI–electron collisions at 10 AU, we assume an electron and proton temperature, T_s ∼ T_e ∼ 20,000 K. Consequently, $\tau _s^{pe} \sim 8.9 \times 10^8$, which is ∼30 times greater than the characteristic supersonic flow time τ_f. These estimates neglect the decreasing solar wind number density with increasing heliocentric distance and the possibility that electrons are hotter than solar wind protons, both of which will increase $\tau _s^{pe}$. Thus, neither proton nor electron collisions can equilibrate a PUI-thermal solar wind plasma in the supersonic solar wind before the heliospheric termination shock.

The inner heliosheath—the subsonic solar wind. The situation in the IHS, i.e., the region between the heliospheric termination shock (HTS) and the heliopause (HP), is more complicated. The supersonic solar wind is decelerated on crossing the quasi-perpendicular HTS. The flow velocity is directed away from the radial direction and is ∼100 km s⁻¹. The interplanetary magnetic field remains approximately perpendicular to the plasma flow. Voyager 2 measured the downstream solar wind temperature to be in the range of ∼120,000–180,000 K (Richardson 2008; Richardson et al. 2008), which was much less than predicted by simple MHD models. Instead, the thermal energy in the IHS is dominated by pickup ions. There are two primary sources of PUIs in the IHS. The first is interstellar neutrals that drift across the HP and charge exchange with hot solar wind plasma. These newly created ions experience pickup in the IHS plasma in the same way that ions are picked up in the supersonic solar wind. The characteristic energy for PUIs created in this manner is ∼52 eV or ∼6.03 × 10⁵ K, which is about five times hotter than the IHS solar wind protons. The second primary source is the PUIs created in the supersonic solar wind and then convected across the HTS into the IHS. The PUIs convected into the HTS are either transmitted immediately across the HTS or are reflected before transmission (Zank et al. 1996a). PUI reflection was identified by Zank et al. (1996a) as the primary dissipation mechanism at the quasi-perpendicular HTS, with the thermal solar wind protons experiencing comparatively little heating across the HTS. This was confirmed by the Voyager 2 crossing of the HTS (Richardson 2008; Richardson et al. 2008). The transmitted PUIs downstream of the HTS have temperatures ∼9.75 × 10⁶ K (∼0.84 keV) and the reflected protons have a temperature of ∼7.7 × 10⁷ K (∼6.6 keV) (Zank et al. 2010). Clearly, PUIs, those transmitted, reflected, and injected, dominate the thermal energy of the heliosheath, just as PUIs do upstream of the HTS despite being only some 20% of the thermal subsonic solar wind number density. The IHS proton distribution function can be approximated by a three- (Zank et al. 2010; Burrows et al. 2010) or four-component distribution function (Zirnstein et al. 2014), with a relatively cool thermal solar wind Maxwellian distribution and two or three superimposed PUI distributions. In two important papers, Zirnstein et al. (2014) and Desai et al. (2014) have exploited this decomposition of the IHS proton distribution function to great effect in modeling energetic neutral atom (ENA) spectra observed by the IBEX spacecraft at 1 AU (see also Desai et al. 2012, who used the slightly simpler three-distribution model of Zank et al. 2010) to identify multiple proton distribution functions in the IHS and the VLISM. These multiple proton populations are identified as the various PUI populations described above and the thermal solar wind proton population.

For the subsonic solar wind plasma in the IHS, we take the background proton and electron number density to be approximately n_s/e ∼ 0.005 cm⁻¹ and the PUI speed to range from ∼1000 (reflected and transmitted PUIs) to 400 (transmitted-only PUIs) to ∼100 km s⁻¹ (PUIs created in the IHS). This yields a lower limit on the PUI–proton collisional time of $\tau _s^{ps} = (\nu _s^{ps})^{-1} \sim 2.1 \times 10^{10}$ s. For an IHS length scale of 100 AU, the characteristic flow time is τ_f ∼ 1.5 × 10⁸ s. Clearly, in the IHS PUI–proton collisions cannot equilibrate the two proton populations in the upwind region. For PUI–electron collisions, we assume an electron temperature, T_e ∼ 200,000 K. Consequently, $\tau _s^{pe} \sim 2.8 \times 10^{10}$, which far exceeds the characteristic flow time τ_f on a scale of 100 AU. Thus, neither proton nor electron collisions can equilibrate a PUI-thermal solar wind plasma in the subsonic solar wind or IHS on scales smaller than at least 10,000 AU. For the IHS, a multi-component plasma description that discriminates between PUIs and the subsonic solar wind plasma is therefore necessary. First attempts to write down such a model were presented by Avinash & Zank (2007) and Avinash et al. (2009). Fisk & Gloeckler (2014) have argued similarly that both PUIs and anomalous cosmic rays should be regarded as partially decoupled from the background thermal plasma.

VLISM. the interstellar plasma upwind of the heliopause is also mediated by energetic PUIs. It was noted already by Zank et al. (1996b) that energetic neutral H created via charge exchange in the IHS and fast solar wind could "splash" back into the VLISM where they would experience a secondary charge exchange. The secondary charge exchange of hot and/or fast neutral H with cold (∼6300 K; McComas et al. 2012) VLISM protons leads to the creation of a hot or suprathermal PUI population in the VLISM. The heating of the VLISM has been discussed in detail by Zank et al. (2013), who pointed out that the heating of the VLISM plasma would result in an increase of the sound speed with a concomitant weakening or even elimination of the bow shock (yielding instead a bow wave). Zirnstein et al. (2014) have extended the Zank et al. (2010) model by including the multiple PUI populations that contribute to the heating of the VLISM (Zank et al. 2013). For the purposes of this work, we are concerned primarily with neutrals created in both the IHS and supersonic solar wind, i.e., with typical speeds of ∼100 km s⁻¹ or ∼400 km s⁻¹, experiencing secondary charge exchange in the VLISM. The fast/hot neutrals determine the pickup speed in the VLISM. The PUIs therefore form a tenuous (n_p ≃ 5 × 10⁻⁵ cm⁻³; Zirnstein et al. 2014) suprathermal component in the VLISM.

For the VLISM plasma, we take the background proton and electron number density to be approximately n_s/e ∼ 0.08 cm⁻¹ and the PUI speed to be ∼100 km s⁻¹ (i.e., H created in the IHS). This yields a lower limit on the PUI–proton collision time of $\tau _s^{ps} = (\nu _s^{ps})^{-1} \sim 1.3 \times 10^{9}$ s. For a length scale of 75 AU, and typical flow speed of 15 km s⁻¹, the characteristic flow time is τ_f ∼ 7.5 × 10⁸ s. For PUI–electron collisions, we assume an electron temperature, T_e ∼ 10,000 K. Consequently, $\tau _s^{pe} \sim 3.2 \times 10^8$, which is comparable to the characteristic flow time τ_f on a scale of 75 AU. Thus, neither proton nor electron collisions can equilibrate a PUI–thermal plasma in the VLISM on scales smaller than at least 75 AU.

Besides the outer heliosphere and the boundaries with the ISM, the interaction of PUIs with the solar wind is of interest in several planetary environments, both unmagnetized (Venus, Mars, and possibly Pluto), and magnetized (Mercury, Earth (with respect to the AMPTE experiment), Saturn, Jupiter), and at comets (see the review by Coates 2012). For example, PUIs at Venus are created by charge exchange or photoionization from cold neutral ionospheric atoms (e.g., Luhmann 1986). Lu et al. (2013) used Venus Express (VEX) observations to show that the strength of the Venusian bow shock is weaker when PUI production is high. In particular, investigating the interaction of comets with the solar wind has been extraordinarily fruitful in furthering our understanding of PUIs in the solar wind—this includes the missions to Comet Halley and Giacobini–Zinner, and later missions to Grigg–Skjellerup and Borrelly (e.g., Coates 2012). With impending encounters of New Horizons with the Pluto–Charon system (2015) and Rosetta with comet Churumov–Gerasimenko (2014–15), active Venus (VEX) and Mars (Mars Express) missions, and of course Voyager 1 and 2 and IBEX, understanding in more detail the nature of PUI mediated plasmas is of far reaching importance.

The paper is organized as follows. In Section 2, we derive the basic equations describing a PUI mediated plasma. At the first-order, we find that the scattering of PUIs in a turbulent magnetofluid introduces a term analogous to that of heat conduction. The second-order correct set of equations describing PUIs in a multi-fluid context leads to the introduction of the viscous terms that define the PUI stress tensor. We systematically derive the system of multi-fluid equations that describe a background Maxwellian proton and electron plasma plus a non-Maxwellian PUI population. Because we assume Maxwellian distributions for the background protons and electrons, the background plasma contributes no heat flux or stress tensor terms. For completeness, we derive a "single-fluid" description analogous to the equations of magnetohydrodynamics (MHD) that describes a PUI mediated plasma. The "single-fluid" model possesses collisionless heat flux and viscous stress terms, unlike the MHD equations. In Section 3, we derive the dispersion relation for linear waves in a PUI mediated plasma and discuss the general properties of waves in a multi-fluid PUI mediated plasma. In Section 4, we present numerical solutions to the dispersion relation for the supersonic solar wind, the IHS plasma, and the plasma in the VLISM. The multi-fluid waves are also compared to the more familiar two-fluid plasma modes. In Section 5, we present an analysis and numerical solutions of the linearized single-fluid model, presenting results for the VLISM only. In the final section, we discuss and summarize our results.

2.1. First-order Correct Multi-fluid Model: Heat Conduction

In deriving a multi-fluid model that includes PUIs self-consistently, we shall assume that the distribution function for the background protons and electrons are each Maxwellian, which ensures the absence of any heat flux or stress tensor terms for the background plasma. The exact form of the continuity, momentum, and energy equations governing the thermal electrons and protons are therefore given by

$$\begin{eqnarray} &&\frac{\partial n_e}{\partial t} + \nabla \cdot \left(n_e {\bf u}_e \right) = 0; \end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray} &&m_e n_e \left(\frac{\partial {\bf u}_e }{\partial t} + {\bf u}_e \cdot \nabla {\bf u}_e \right) = -\nabla P_e - en_e \left({\bf E} + {\bf u}_e \times {\bf B} \right); \end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray} &&\frac{\partial P_e}{\partial t} + {\bf u}_e \cdot \nabla P_e + \gamma _e P_e \nabla \cdot {\bf u}_e = 0, \end{eqnarray} \tag{ 6 }$$

for the electrons, and

$$\begin{eqnarray} &&\frac{\partial n_s}{\partial t} + \nabla \cdot \left(n_s {\bf u}_s \right) = 0; \end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray} &&m_p n_s \left(\frac{\partial {\bf u}_s }{\partial t} + {\bf u}_s \cdot \nabla {\bf u}_s \right) = -\nabla P_s + en_s \left({\bf E} + {\bf u}_s \times {\bf B} \right); \end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray} &&\frac{\partial P_s}{\partial t} + {\bf u}_s \cdot \nabla P_s + \gamma _s P_s \nabla \cdot {\bf u}_s = 0, \end{eqnarray} \tag{ 9 }$$

for the protons. Here n_e/s, u_e/s, and P_e/s are the usual macroscopic fluid variables for the electron/proton number density, velocity, and pressure respectively, γ_e/s the electron/proton adiabatic index, E the electric field, B the magnetic field, and e the charge of an electron.

Pickup ions initially form an unstable distribution that excites Alfvénic fluctuations. The self-generated fluctuations and in situ turbulence serve to scatter PUIs in pitch angle. The Alfvén waves and magnetic field fluctuations both propagate and convect with the bulk velocity of the system U = U(u_e, u_s, u_p, n_e, n_s, n_p, m_e, m_p), where n_p and u_p refer to PUI variables. The PUIs are governed by the Fokker–Planck transport equation with a (for now unspecified) collisional term δf/δt|_c,

$$\begin{equation} \frac{\partial f}{\partial t} + {\bf v} \cdot \nabla f + \frac{e}{m_p} \left({\bf E} + {\bf v} \times {\bf B} \right) \cdot \nabla _v f = \left. \frac{\delta f}{\delta t} \right|_c, \end{equation} \tag{ 10 }$$

for average electric and magnetic fields E and B. We assume that the velocity v of PUIs is always non-relativistic. The transport Equation (10) has to be transformed into a frame that ensures there is no change in PUI momentum and energy due to scattering. For the present, assume that the cross-helicity σ is nonzero and let

$$\begin{equation} {\bf v} = {\bf c} + {\bf U} + \sigma {\bf V}_A \quad \Longleftrightarrow \quad {\bf c} = {\bf v} - {\bf U} - \sigma {\bf V}_A, \end{equation} \tag{ 11 }$$

where V_A is the Alfvén velocity and c is the random velocity. The transport equation is therefore

$$\begin{eqnarray} && \frac{\partial f}{\partial t} + \left(U_i + \sigma V_{Ai} + c_i \right) \frac{\partial f}{\partial x_i} + \left [ \frac{e}{m_p} \left({\bf E} + {\bf U} \times {\bf B} \right)_i\right. \nonumber \\ &&\left. \quad + \frac{e}{m_p} \left({\bf c} \times {\bf B} \right)_i - \frac{\partial U_i}{\partial t} {} - (U_j + \sigma V_{Aj} + c_j) \frac{\partial U_i}{\partial x_j}\right. \nonumber \\ &&\left. \quad - \sigma \left(\frac{\partial V_{Ai} }{\partial t} - \sigma (U_j + \sigma V_{Aj} + c_j) \frac{\partial V_{Ai} }{\partial x_j} \right) \right] \frac{\partial f}{\partial c_i} = \left. \frac{\delta f}{\delta t} \right|_c. \nonumber\\ \end{eqnarray} \tag{ 12 }$$

The velocity U is still unspecified, so we choose U such that E' ≡ E + U × B = 0. This assumption corresponds to choosing

$$\begin{equation} {\bf U}_{\perp } = {\bf U} - {\bf U}_{\parallel } = \frac{ {\bf E} \times {\bf B}}{B^2} \equiv {\bf U}, \end{equation} \tag{ 13 }$$

since we choose U_∥ = 0 (U_∥ is parallel to B and therefore arbitrary). This corresponds to expressing (10) in the guiding center frame. The transformation to the velocity U then yields

$$\begin{eqnarray} && \frac{\partial f}{\partial t} + \left(U_i + c_i \right) \frac{\partial f}{\partial x_i} + \left[ \frac{e}{m_p} \left({\bf c} \times {\bf B} \right)_i - \frac{\partial U_i}{\partial t}\right. \nonumber\\ &&\left. \quad - (U_j + c_j) \frac{\partial U_i}{\partial x_j} \right] \frac{\partial f}{\partial c_i} = \left. \frac{\delta f}{\delta t} \right|_c, \end{eqnarray} \tag{ 14 }$$

after setting the cross-helicity σ = 0. By taking moments of (14), we can derive the evolution equations for the macroscopic PUI variables, such as the number density n_p = ∫fd³c, velocity n_pu_pi = ∫c_ifd³c, and so on. Although unspecified for now, we shall assume that moments of the collisional term δf/δt|_c are zero. This can be checked against the particular scattering model that we use below. The zeroth moment of (14) yields the continuity equation for PUIs,

$$\begin{equation} \frac{\partial n_p}{\partial t} + \frac{\partial }{\partial x_i} (n_p (U_i + {u_p}_i)) = 0, \end{equation} \tag{ 15 }$$

where u_p is the PUI bulk velocity in the guiding center frame. For the first moment, we multiply (14) by c_j and integrate over velocity space. This yields, after a little algebra,

$$\begin{eqnarray} && \frac{\partial }{\partial t} (n_p (U_j + {u_p}_j)) + \nabla \cdot [ n_p {\bf U} (U_j + {u_p}_j) + n_p {\bf u}_p U_j ] \nonumber\\ && \qquad + \frac{\partial }{\partial x_i} \int c_i c_j f d^3c = \frac{e}{m_p} n_p \varepsilon _{jkl} {u_p}_k B_l, \end{eqnarray} \tag{ 16 }$$

where ε_ijk is the Levi–Civeta tensor.

To close Equation (16), we need to evaluate the PUI distribution function f, which requires that we solve (14). In solving (14), we assume (1) that the PUI distribution is gyrotropic, and (2) that scattering of PUIs is sufficiently rapid to ensure that the PUI distribution is nearly isotropic. We can therefore average (14) over gyrophase, obtaining the so-called focused transport equation for non-relativistic particles (Isenberg 1997). Details of the derivation can be found in Chapter 5 of Zank (2013), and the explicit expression is given in Appendix A. To solve the gyrophase-averaged transport equation requires that we specify the scattering or collisional operator. We make the simplest possible choice, which is the isotropic pitch-angle diffusion operator,

$$\begin{equation} \frac{\partial }{\partial \mu } \left(\nu _s(1 - \mu ^2) \frac{\partial f}{\partial \mu } \right), \end{equation} \tag{ 17 }$$

where μ = cos θ is the cosine of the particle pitch-angle θ and $\nu _s = \tau _s^{-1}$ is the scattering frequency. The form of the scattering operator (17) allows us to solve the focused transport equation (A1) using a Legendre polynomial expansion of the distribution function f. This is summarized in Appendix A and details can be found in Chapter 5 of Zank (2013). The first-order correct solution to the gyrophase-averaged form of Equation (14), i.e., (A1), is

$$\begin{eqnarray} &&f \simeq f_0 + \mu f_1; \end{eqnarray} \tag{ 18 }$$

$$\begin{eqnarray} &&f_0 = f_0 ({\bf x}, c, t); \end{eqnarray} \tag{ 19 }$$

$$\begin{eqnarray} &&f_1 = -\frac{c \tau _s}{3} b_i \frac{\partial f_0}{\partial x_i} + \frac{DU_i}{Dt} \frac{\tau _s}{3} b_i \frac{\partial f_0}{\partial c}, \end{eqnarray} \tag{ 20 }$$

where c = |c| is the particle random speed, b ≡ B/B is a directional unit vector defined by the magnetic field, and D/Dt ≡ ∂/∂t + U_i∂/∂x_i is the convective derivative. Both f₀ and f₁ are functions of position, time, and particle random speed c, i.e., independent of pitch-angle μ (and of course gyrophase ϕ). Of particular importance is the retention of the large-scale velocity U acceleration and shear terms. These terms are often neglected in the derivation of the transport equation describing f₀ (for relativistic particles, the transport equation is the familiar cosmic ray transport equation). Thus, the second term in (20) is typically neglected, although it is known as the relativistic heat inertia term in the relativistic transport theory of cosmic rays (Webb 1985, 1987, 1989). As will be seen below, retaining these terms is absolutely essential to derive the correct multi-fluid formulation for PUIs. By introducing

$$\begin{eqnarray*} \int c_i c_j f d^3c &=& \int (c_i - {u_p}_i) (c_j - {u_p}_j) f d^3c + n_p {u_p}_i {u_p}_j \nonumber\\ &\equiv & \int c_i^{\prime } c_j^{\prime } f d^3c + n_p {u_p}_i {u_p}_j \\ &\simeq & \int c_i^{\prime } c_j^{\prime } (f_0 + \mu f_1) d^3c + n_p {u_p}_i {u_p}_j, \end{eqnarray*}$$

we can show that

$$\begin{eqnarray} \frac{\partial }{\partial x_i} \int c_i^{\prime } c_j^{\prime } f_0 d^3c &=& \frac{1}{m_p} \frac{\partial }{\partial x_i} (\delta _{ij} P_p), \quad {\rm and} \quad \nonumber\\ \int c_i^{\prime } c_j^{\prime } \mu f_1 d^3c &=& 0, \end{eqnarray} \tag{ 21 }$$

where

$$\begin{equation} P_p \equiv m_p \frac{4\pi }{3} \int {c^{\prime } }^2 f_0 {c^{\prime } }^2 dc. \end{equation} \tag{ 22 }$$

Consequently, the PUI stress tensor is identically zero at first-order and there exists only an isotropic pressure tensor δ_ijP_p. We show in the following section that retaining the second-order terms in the Legendre polynomial expansion of the gyrophase-averaged equation (A2) does in fact yield a non-zero collisionless stress tensor. The PUI momentum equation to first-order can therefore be expressed as

$$\begin{eqnarray} && \frac{\partial }{\partial t} (n_p (U_j + {u_p}_j)) + \frac{\partial }{\partial x_i} \left[ n_p ({\bf U} + {\bf u}_p) (U_j + {u_p}_j) + \frac{1}{m_p} \delta _{ij} P_p \right] \nonumber\\ && \quad = \frac{e}{m_p} n_p \varepsilon _{jkl} {u_p}_k B_l. \end{eqnarray} \tag{ 23 }$$

To derive the transport equation for P_p, we multiply (14) by (1/2)c² and integrate over d³c. We then use (18)–(20) to evaluate the various integrals. Introducing c' ≡ c − u_p as before, we find

$$\begin{eqnarray*} &&\int \frac{1}{2} c^2 f_0 d^3c = \frac{3}{2} \frac{1}{m_p} P_p + \frac{1}{2} n_p u_p^2, \end{eqnarray*}$$

for example. Similarly, we find that the heat flux q(x, t) can be expressed as

$$\begin{eqnarray} q_i ({\bf x}, t) &\equiv & \int \frac{1}{2} {c^{\prime } }^2 c_i^{\prime } f d^3c^{\prime } = \frac{1}{2} \int c^2 c_i f d^3c \nonumber\\ && - \frac{5}{2} \frac{1}{m_p} {u_p}_i P_p - \frac{1}{2} n_p u_p^2 {u_p}_i. \end{eqnarray} \tag{ 24 }$$

It then follows that

$$\begin{eqnarray*} &&\int \frac{1}{2} {c^{\prime } }^2 c_i^{\prime } f_0 d^3c = \pi \int {c^{\prime } }^3 \mu b_i f_0 {c^{\prime } }^2 dc^{\prime } = 0, \end{eqnarray*}$$

and

$$\begin{eqnarray} \int \frac{1}{2} {c^{\prime } }^2 c_i^{\prime } \mu f_1 d^3c^{\prime } &=& {-}\frac{2 \pi }{3} \int {c^{\prime } }^2 \kappa _{ij} \frac{\partial f_0}{\partial x_j} {c^{\prime } }^2 dc^{\prime } \nonumber \\ &=& {-}\frac{1}{2} K_{ij} \frac{\partial P_p }{\partial x_j} = q_i ({\bf x},t). \end{eqnarray} \tag{ 25 }$$

In (25), we introduced the spatial diffusion coefficient

$$\begin{equation} \kappa _{ij} \equiv b_i \frac{c^2 \tau _s}{3} b_j, \end{equation} \tag{ 26 }$$

together with PUI speed-averaged form K_ij. The collisionless heat flux for PUIs is therefore described in terms of the PUI pressure gradient and consequently the averaged spatial diffusion introduces a PUI diffusion time and length scale into the multi-fluid system. The diffusion coefficient, i.e., the coefficient for the PUI heat flux, is proportional to the particle scattering time τ_s, and therefore a function of the background turbulent intensity. A separate calculation, possibly based on quasi-linear theory for the parallel diffusion coefficient or the nonlinear guiding center theory for the perpendicular diffusion coefficient, is necessary to obtain reasonable estimates of the scattering time (Matthaeus et al. 2003; Zank et al. 2004).

The remaining terms are straightforwardly evaluated. We find

$$\begin{eqnarray*} && \frac{e}{m_p} \varepsilon _{ijk} B_k \int \frac{1}{2} c^2 c_j \frac{\partial f}{\partial c_i} d^3c = - \frac{e}{m_p} \varepsilon _{ijk} n_p {u_p}_i B_k {u_p}_j; \\ && \quad -\frac{DU_i}{Dt} \int \frac{1}{2} c^2 \frac{\partial f}{\partial c_i} d^3c = n_p {u_p}_i \frac{DU_i}{Dt}; \\ && \quad \frac{\partial U_i}{\partial x_j} \int \frac{1}{2} c^2 c_j \frac{\partial f}{\partial c_i} d^3c = \frac{5}{2} P_p \frac{\partial U_i}{\partial x_i} + \frac{1}{2} n_p u_p^2 \frac{\partial U_i}{\partial x_i}\\ && \quad + n_p {u_p}_i {u_p}_j \frac{\partial U_i}{\partial x_j}. \end{eqnarray*}$$

On combining these results, we obtain, after some algebra, the transport equation for the PUI pressure

$$\begin{eqnarray} &&\frac{\partial P_p}{\partial t} + ({\bf u}_p + {\bf U}) \cdot \nabla P_p + \frac{5}{3} P_p \nabla \cdot ({\bf u}_p + {\bf U}) = \frac{1}{3} \nabla \cdot ({\bf K} \cdot \nabla P_p),\nonumber \\ \end{eqnarray} \tag{ 27 }$$

illustrating that the PUI heat flux yields a spatial diffusion term in the PUI equation of state. The PUI system of equations is properly closed and correct to the first-order. The second-order correct PUI equations, which includes the PUI stress tensor, is given in the following subsection. For completeness, the PUI total energy equation has the form

$$\begin{eqnarray} && \frac{\partial }{\partial t} \left(\frac{3}{2} P_p + \frac{1}{2} n_p (u_p + U)^2 \right) + \frac{\partial }{\partial x_i} \left[ \frac{1}{2} n_p (u_p + U)^2 ({u_p}_i + U_i)\right. \nonumber\\ &&\left. + \frac{5}{2} P_p ({u_p}_i + U_i) - \frac{1}{2} K_{ij} \frac{\partial P_p}{\partial x_j} \right] {} = \frac{e}{m_p} \varepsilon _{ijk} n_p {u_p}_j B_k ({u_p}_i + U_i). \nonumber\\ \end{eqnarray} \tag{ 28 }$$

The full system of PUI equations is given by (15), (23), and (27) or (28). It is not particularly illuminating to work in the guiding center frame, and we may simplify (15), (23), and (27), (28), by letting

$$\begin{eqnarray*} &&{\bf U}_p = {\bf u}_p + {\bf U}. \end{eqnarray*}$$

The right-hand side (RHS) of Equations (23) and (28) is proportional to u_p × B, which becomes

$$\begin{eqnarray*} &&\left({\bf U}_p - {\bf U} \right) \times {\bf B} = {\bf E} + {\bf U}_p \times {\bf B}, \end{eqnarray*}$$

since E was perpendicular to B by construction initially. Hence the PUI fluid equations can be written in the more familiar form

$$\begin{eqnarray} \frac{\partial n_p}{\partial t} + \nabla \cdot (n_p {\bf U}_p) &=& 0; \end{eqnarray} \tag{ 29 }$$

$$\begin{eqnarray} \frac{\partial }{\partial t} (n_p {\bf U}_p) + \nabla \cdot [ n_p {\bf U}_p {\bf U}_p + {\bf I} P_p ] &=& \frac{e}{m_p} n_p ({\bf E} + {\bf U}_p \times {\bf B}); \nonumber\\ \end{eqnarray} \tag{ 30 }$$

$$\begin{eqnarray} && \frac{\partial }{\partial t} \left(\frac{3}{2} P_p + \frac{1}{2} n_p U_p^2 \right) + \nabla \cdot \left[ \frac{1}{2} n_p U_p^2 {\bf U}_p + \frac{5}{2} P_p {\bf U}_p \right.\nonumber\\ &&\left. \quad\quad - \frac{1}{2} {\bf K}\cdot \nabla P_p \right] = \frac{e}{m_p} n_p {\bf U}_p \cdot ({\bf E} + {\bf U}_p \times {\bf B}), {} \end{eqnarray} \tag{ 31 }$$

which is the form we use below. Similarly, we have

$$\begin{equation} \frac{\partial P_p}{\partial t} + {\bf U}_p \cdot \nabla P_p + \frac{5}{3} P_p \nabla \cdot {\bf U}_p = \frac{1}{3} (\nabla \cdot {\bf K} \cdot \nabla P_p). \end{equation} \tag{ 32 }$$

The full thermal electron-thermal proton–PUI multi-fluid system is therefore given by Equations (4)–(9) and (29)–(31) or (32), together with Maxwell's equations

$$\begin{eqnarray} &&\frac{\partial {\bf B} }{\partial t} = - \nabla \times {\bf E}; \end{eqnarray} \tag{ 33 }$$

$$\begin{eqnarray} &&\nabla \times {\bf B} = \mu _0 {\bf J}; \end{eqnarray} \tag{ 34 }$$

$$\begin{eqnarray} &&\nabla \cdot {\bf B} = 0; \end{eqnarray} \tag{ 35 }$$

$$\begin{eqnarray} &&{\bf J} = e (n_s {\bf u}_s + n_p {\bf U}_p - n_e {\bf u}_e), \end{eqnarray} \tag{ 36 }$$

where J is the current and μ₀ the permeability of free space.

2.2. Second-order Correct Multi-fluid Model: the Stress Tensor

As shown above, the zeroth- and first-order solutions for the pressure tensor yields an isotropic scalar pressure P_ij = P_pδ_ij only. Consider now the second-order Legendre polynomial expansion of f,

$$\begin{equation} f \simeq f_0 + \mu f_1 + \frac{1}{2} (3 \mu ^2 - 1) f_2. \end{equation} \tag{ 37 }$$

As before, we need to evaluate

$$\begin{eqnarray*} \int c_i c_j f d^3c &=& \int c^{\prime }_i c^{\prime }_j \left(f_0 + \mu f_1 + \frac{1}{2} (3 \mu ^2 - 1) f_2 \right)\nonumber\\ && \times d^3 c^{\prime } + n_p {u_p}_i {u_p}_j, \end{eqnarray*}$$

from which we find

$$\begin{eqnarray*} \int c^{\prime }_i c^{\prime }_j f_0 d^3 c^{\prime } &=& \frac{1}{m_p} P_p \delta _{ij}; \\ \int c^{\prime }_i c^{\prime }_j \mu f_1 d^3 c^{\prime } &=& 0. \end{eqnarray*}$$

Although not discussed explicitly above, since the PUI pressure is defined in the frame of the bulk PUI velocity u_p, the distribution function over which the integral is taken needs to evaluated in this frame. Since the expression (A7) for f₂ is a function of the guiding center velocity U, we need to transform to the frame U' = U + u_p. On using the solution (A7) for f₂, we obtain

$$\begin{eqnarray} \int {c_x^{\prime }}^2 \frac{1}{2} (3\mu ^2 - 1) f_2 d^3 c^{\prime } &=& \int {c_y^{\prime }}^2 \frac{1}{2} (3\mu ^2 - 1) f_2 d^3 c^{\prime } \nonumber\\ &=& \frac{\eta }{15} \left(b_i b_j \frac{\partial U_j^{\prime } }{\partial x_i} - \frac{1}{3} \frac{\partial U_i^{\prime } }{\partial x_i} \right);\quad\quad \end{eqnarray} \tag{ 38 }$$

$$\begin{eqnarray} &&\int {c_z^{\prime }}^2 \frac{1}{2} (3\mu ^2 - 1) f_2 d^3 c^{\prime } = -\frac{2 \eta }{15} \left(b_i b_j \frac{\partial U_j^{\prime } }{\partial x_i} - \frac{1}{3} \frac{\partial U_i^{\prime } }{\partial x_i} \right);\quad\quad \end{eqnarray} \tag{ 39 }$$

$$\begin{eqnarray} &&\int c_i^{\prime } c_j^{\prime } \frac{1}{2} (3\mu ^2 - 1) f_2 d^3 c^{\prime } = 0, \quad (i \ne j), \end{eqnarray} \tag{ 40 }$$

where the coefficient of viscosity η is defined as

$$\begin{eqnarray} \eta &\equiv & \frac{4 \pi }{15} \int \frac{\partial }{\partial c^{\prime } } ({c^{\prime } }^4 c \tau _s) f_0 dc^{\prime } \end{eqnarray} \tag{ 41 }$$

$$\begin{eqnarray} &\simeq & \frac{4 \pi }{3} \int {c^{\prime } }^2 \tau _s f_0 {c^{\prime } }^2 dc^{\prime } \end{eqnarray} \tag{ 42 }$$

$$\begin{eqnarray} &\simeq & \frac{P_p \tau _s }{m_p}. \end{eqnarray} \tag{ 43 }$$

Equation (41) is the formal definition of the coefficient of viscosity for the PUI gas. If we assume (probably reasonably) that |c| ≫ |u_p|, then we obtain (42), which may be regarded as a PUI pressure moment weighted by the PUI scattering time. Finally, if we assume that τ_s is independent of c, we then obtain the "classical" form (43) of the viscosity coefficient. The pressure tensor may therefore be expressed as

$$\begin{equation} (P_{ij}) = P_p (\delta _{ij}) + \left(\begin{array}{ccc}1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & -2 \end{array} \right) \frac{\eta }{15} \left(b_k b_{\ell } \frac{\partial U_k^{\prime } }{\partial x_{\ell } } - \frac{1}{3} \frac{\partial U_m^{\prime } }{\partial x_m} \right). \end{equation} \tag{ 44 }$$

The pressure tensor may be written in a more revealing form if we introduce a "viscosity matrix,"

$$\begin{equation} \left(M_{k \ell } \right) \equiv \left(\eta _{k \ell } \right) = \left(\frac{\eta }{15} b_k b_{\ell } \right) \simeq \left(\frac{1}{15} \frac{P_p \tau _s b_k b_{\ell } }{m_p} \right), \end{equation} \tag{ 45 }$$

and note that η_ij = η_ji and η/15 = η₁₁ + η₂₂ + η₃₃ = η_ijδ_ij (since b² = 1). Then

$$\begin{eqnarray} && \frac{\eta }{15} \left(b_k b_{\ell } \frac{\partial U_k^{\prime } }{\partial x_{\ell } } - \frac{1}{3} \frac{\partial U_m^{\prime } }{\partial x_m} \right) = \frac{\eta _{k \ell } }{2} \left(\frac{\partial U_k^{\prime } }{\partial x_{\ell } } + \frac{\partial U_{\ell }^{\prime } }{\partial x_k } \right) \nonumber \\ && \quad - \frac{1}{3} \eta _{k \ell } \delta _{k \ell } \frac{\partial U_m^{\prime } }{\partial x_m } = \frac{\eta _{k \ell } }{2} \left(\frac{\partial U_k^{\prime } }{\partial x_{\ell } } + \frac{\partial U_{\ell }^{\prime } }{\partial x_k } - \frac{2}{3} \delta _{k \ell } \frac{\partial U_m^{\prime } }{\partial x_m } \right), \quad\quad \end{eqnarray} \tag{ 46 }$$

which yields the pressure tensor as the sum of an isotropic scalar pressure P_p and the stress tensor, i.e.,

$$\begin{eqnarray} (P_{ij}) &=& P_p (\delta _{ij}) + \left(\begin{array}{ccc}1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & -2 \end{array} \right) \frac{\eta _{k \ell } }{2} \nonumber\\ && \times \left(\frac{\partial U_k^{\prime } }{\partial x_{\ell } } + \frac{\partial U_{\ell }^{\prime } }{\partial x_k } - \frac{2}{3} \delta _{k \ell } \frac{\partial U_m^{\prime } }{\partial x_m } \right) \equiv P_p {\bf I} + \Pi _p. \quad\quad \end{eqnarray} \tag{ 47 }$$

The stress tensor is a generalization of the "classical" form in that several coefficients of viscosity are present, and of course the derivation here is for a collisionless charged gas of PUIs experiencing only pitch-angle scattering by turbulent magnetic fluctuations. Use of the pressure tensor (47) yields a "Navier–Stokes"-like modification of the PUI momentum equation,

$$\begin{eqnarray} && \frac{\partial }{\partial t} (n_p {\bf U}_p) + \nabla \cdot \left[ n_p {\bf U}_p {\bf U}_p + \frac{1}{m_p} {\bf I} P_p \right] = \frac{e}{m_p} n_p ({\bf E} + {\bf U}_p \times {\bf B}) \nonumber \\ && \quad {} - \frac{1}{m_p} \nabla \cdot \left(\begin{array}{ccc}1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & -2 \end{array} \right) \frac{\eta _{k \ell } }{2} \left(\frac{\partial {U_p}_k }{\partial x_{\ell } } + \frac{\partial {U_p}_{\ell } }{\partial x_k }\right. \nonumber\\ &&\left. \quad - \frac{2}{3} \delta _{k \ell } \frac{\partial {U_p}_m }{\partial x_m } \right) = \frac{e}{m_p} n_p ({\bf E} + {\bf U}_p \times {\bf B}) - \frac{1}{m_p} \nabla \cdot \Pi _p, \end{eqnarray} \tag{ 48 }$$

where we used U_p = u_p + U ≡ U' as before. The full momentum equation with the second-order stress tensor correction is included for completeness but in the linearized wave analysis below, we use only the first-order correct equations, i.e., only the heat conduction term is included.

2.3. Reduced "Single-fluid" Model

For some problems, such as the investigation of turbulence in the outer heliosphere, IHS, or VLISM, the full multi-fluid model is far too complicated to solve. By making the key assumption that U_p ≃ u_s, we can reduce the multi-fluid system above to an MHD-like set of model equations. The assumption that U_p ≃ u_s is quite reasonable since (1) the bulk flow velocity of a plasma is always dominated by the protons, and (2) the pick-up process itself forces newly created PUIs to essentially co-move with the background plasma flow. Accordingly, we let U_p ≃ u_s = U_i be the overall proton (i.e., thermal background protons and PUIs) velocity. The thermal proton and PUI continuity and momentum equations are therefore trivially combined as

$$\begin{eqnarray} &&\frac{\partial n_i}{\partial t} + \nabla \cdot \left(n_i {\bf U}_i \right) = 0; \end{eqnarray} \tag{ 49 }$$

$$\begin{eqnarray} m_p n_i \left(\frac{\partial {\bf U}_i}{\partial t} + {\bf U}_i \cdot \nabla {\bf U}_i \right) &=& - \nabla (P_s + P_p) + en_i \left({\bf E} + {\bf U}_i \times {\bf B} \right) \nonumber\\ && - \nabla \cdot \Pi _p, \quad \end{eqnarray} \tag{ 50 }$$

where n_i = n_s + n_p. Since the PUIs are not thermally equilibrated with the background plasma (T_s ≠ T_p), we need to deal separately with the P_s and P_p equations. These become

$$\begin{eqnarray} &&\frac{\partial P_s}{\partial t} + {\bf U}_i \cdot \nabla P_s + \gamma _s P_s \nabla \cdot {\bf U}_i = 0; \end{eqnarray} \tag{ 51 }$$

$$\begin{eqnarray} &&\frac{\partial P_p}{\partial t} + {\bf U}_i \cdot \nabla P_p + \gamma _p P_p \nabla \cdot {\bf U}_i = \frac{1}{3} \nabla \cdot ({\bf K} \cdot \nabla P_p). \end{eqnarray} \tag{ 52 }$$

By combining the proton Equations (49)–(52) with the electron Equations (4)–(6), we can obtain an MHD-like system of equations. On defining new macroscopic variables,

$$\begin{eqnarray} \rho &\equiv & m_e n_e + m_p n_i; \nonumber \\ q &\equiv & -e (n_e - n_i); \nonumber \\ \rho {\bf U} &\equiv & m_e n_e {\bf u}_e + m_p n_i {\bf U}_i; \nonumber \\ {\bf J} &\equiv & -e \left(n_e {\bf u}_e - n_i {\bf U}_i \right), \end{eqnarray} \tag{ 53 }$$

we can express

$$\begin{eqnarray} n_e &=& \frac{\rho - (m_p /e) q }{m_p(1 - \xi) } \simeq \rho /m_p; \nonumber \\ n_i &=& \frac{\rho + \xi (m_p /e) q }{m_p(1 + \xi) } \simeq \rho /m_p; \nonumber \\ {\bf u}_e &=& \frac{\rho {\bf U} - (m_p/e) {\bf J} }{\rho - (m_p/e) q } \simeq {\bf U} - \frac{m_p}{e} \frac{\bf J}{\rho }; \nonumber \\ {\bf u}_i &=& \frac{\rho {\bf U} + \xi (m_p/e) {\bf J} }{\rho + \xi (m_p/e) q } \simeq {\bf U}, \end{eqnarray} \tag{ 54 }$$

where the smallness of the mass ratio ξ ≡ m_e/m_p ≪ 1 has been exploited. Use of the approximations (54) allows us to combine the continuity and momentum equations in the usual way and to rewrite the thermal electron and proton pressure in terms of the single-fluid macroscopic variables. Thus,

$$\begin{eqnarray} &&\frac{\partial \rho }{\partial t} + \nabla \cdot \left(\rho {\bf U} \right) = 0; \end{eqnarray} \tag{ 55 }$$

$$\begin{eqnarray} \rho \left(\frac{\partial {\bf U} }{\partial t} + {\bf U} \cdot \nabla {\bf U} \right) &=& -\nabla (P_e + P_s + P_p) \nonumber\\ && + {\bf J} \times {\bf B} - \nabla \cdot \Pi; \end{eqnarray} \tag{ 56 }$$

$$\begin{eqnarray} &&\frac{\partial P_s }{\partial t} + {\bf U} \cdot \nabla P_s + \gamma _s P_s \nabla \cdot {\bf U} = 0; \end{eqnarray} \tag{ 57 }$$

$$\begin{eqnarray} && \frac{\partial P_e }{\partial t} + {\bf U} \cdot \nabla P_e + \gamma _e P_e \nabla \cdot {\bf U} = \frac{m_p}{e \rho } {\bf J} \cdot \nabla P_e \nonumber\\ && \quad \quad + \frac{\gamma _e m_p}{e} P_e \nabla \cdot \left(\frac{\bf J}{\rho } \right), \end{eqnarray} \tag{ 58 }$$

where

$$\begin{eqnarray*} \Pi _{k \ell } = \left(\begin{array}{ccc}1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & -2 \end{array} \right) \frac{\eta _{k \ell } }{2} \left(\frac{\partial U_k }{\partial x_{\ell } } + \frac{\partial U_{\ell } }{\partial x_k } - \frac{2}{3} \delta _{k \ell } \frac{\partial U_m }{\partial x_m } \right). \end{eqnarray*}$$

Since we may assume that the current density is much less than the momentum flux, i.e., |J| ≪ |ρU|, we can simplify (58) further by neglecting the RHS. By assuming that γ_e = γ_s = γ, we can combine the thermal proton and electron equations in a single thermal plasma pressure equation with P ≡ P_e + P_s,

$$\begin{equation} \frac{\partial P }{\partial t} + {\bf U} \cdot \nabla P + \gamma P \nabla \cdot {\bf U} = 0. \end{equation} \tag{ 59 }$$

Note that at this point, no assumptions about either the thermal electron or proton pressures (or temperatures) have been made.

Finally, we need an equation for the electric field E. To do so, we multiply the respective momentum equations by the electron or proton charge, sum, and use the approximations (54) to obtain

$$\begin{eqnarray*} && \xi \left(\frac{m_p}{e} \right)^2 \frac{1}{\rho } \left[ \frac{\partial {\bf J} }{\partial t} + \nabla \cdot ({\bf J U} + {\bf U J}) \right] = \frac{m_p}{e \rho } (\nabla P_e - {\bf J} \times {\bf B} \nonumber\\ && \quad - \xi \nabla (P_s + P_p) - \xi \nabla \cdot \Pi) + {\bf E} + {\bf U} \times {\bf B}. \end{eqnarray*}$$

The generalized Ohm's law is therefore

$$\begin{equation} {\bf E} = -{\bf U} \times {\bf B} - \frac{m_p}{e \rho } (\nabla P_e - {\bf J} \times {\bf B} - \xi \nabla P_p), \end{equation} \tag{ 60 }$$

where we have retained the PUI pressure since in principle it can be a high temperature component of the plasma system and ξP_p may be comparable to the P_e term. For our cases of interest, however, T_p ∼ 10⁶ K, T_e ∼ 10⁴ K in the supersonic solar wind, T_p ∼ 10⁷ − 10⁸ K, T_e ∼ 10⁵ K in the IHS, and T_p ∼ 10⁶ K, T_e ∼ 10⁴ K in the VLISM, and n_p < 0.2n_e. Thus, we find that the P_p term can be neglected in Ohm's law (60). Neglect of the electron pressure and Hall current term then yields the usual form of Ohm's law.

The reduced single-fluid model equations may therefore be summarized as

$$\begin{eqnarray} &&\frac{\partial \rho }{\partial t} + \nabla \cdot \left(\rho {\bf U} \right) = 0; \end{eqnarray} \tag{ 61 }$$

$$\begin{eqnarray} &&\rho \left(\frac{\partial {\bf U} }{\partial t} + {\bf U} \cdot \nabla {\bf U} \right) = -\nabla (P + P_p) + {\bf J} \times {\bf B} - \nabla \cdot \Pi; \quad \end{eqnarray} \tag{ 62 }$$

$$\begin{eqnarray} &&\frac{\partial P }{\partial t} + {\bf U} \cdot \nabla P + \gamma P \nabla \cdot {\bf U} = 0; \end{eqnarray} \tag{ 63 }$$

$$\begin{eqnarray} &&\frac{\partial P_p }{\partial t} + {\bf U} \cdot \nabla P_p + \gamma _p P_p \nabla \cdot {\bf U} = \frac{1}{3} \nabla \cdot ({\bf K} \cdot \nabla P_p); \end{eqnarray} \tag{ 64 }$$

$$\begin{eqnarray} &&{\bf E} = -{\bf U} \times {\bf B}; \end{eqnarray} \tag{ 65 }$$

$$\begin{eqnarray} &&\frac{\partial {\bf B} }{\partial t} = - \nabla \times {\bf E}; \end{eqnarray} \tag{ 66 }$$

$$\begin{eqnarray} &&\mu _0 {\bf J} = \nabla \times {\bf B}; \end{eqnarray} \tag{ 67 }$$

$$\begin{eqnarray} &&\nabla \cdot {\bf B} = 0. \end{eqnarray} \tag{ 68 }$$

The single-fluid description (61)–(68) differs from the standard MHD model in that a separate description for the PUI pressure (64) is required. The PUIs introduce both a collisionless heat conduction and viscosity into the system. Notice that the derived electric field (65) is consistent with that used in the transformation of the PUI Fokker–Planck equation (10) i.e., the bulk velocity U derived for the reduced system corresponds to that used in the velocity transformation to eliminate the electric field in (10).

The model Equations (61)–(68), despite being appropriate to non-relativistic PUIs, are interestingly related to the so-called two-fluid MHD system of equations used to describe cosmic ray mediated plasmas (Webb 1983). However, the derivation of the two models is substantially different in that the cosmic ray number density is explicitly neglected in the two-fluid cosmic ray model and a Chapman–Enskog derivation is not used in deriving the cosmic ray hydrodynamic equations. Nonetheless, the sets of equations that emerge are the same suggesting that a more careful analysis of the cosmic ray two-fluid equations would show that the system can in fact include the cosmic ray number density explicitly and therefore be capable of describing the full plasma system.

The PUI mediated plasma model derived above lends itself to a variety of applications. Here, we analyze linear waves in the three regions discussed in the Introduction, viz., the supersonic solar wind in the outer heliosphere, the subsonic solar wind in the IHS, and the VLISM. As we discuss below, the wave properties in these three regions are significantly different from those in the inner heliosphere, which is essentially unmediated by PUI physics. The closest correspondence is to multi-fluid models that include heavy ions as a minority component (Mann et al. 1997; Verscharen et al. 2013), but the heavy ions do not typically dominate the thermal energy budget of the plasma system.

The linearization of the PUI multi-fluid equations proceeds in the standard way and we assume

$$\begin{eqnarray*} n_{e/s/p} &=& n_{e0/s0/p0} + \delta n_{e/s/p}; \quad {\bf u}_{e/s} = \delta {\bf u}_{e/s}; \quad \nonumber\\ {\bf U}_p &=& \delta {\bf U}_p; \quad P_{e/s/p} = P_{e0/s0/p0} + \delta P_{e/s/p}; \\ {\bf B} &=& {\bf B}_0 + \delta {\bf B} = B_0 \hat{\bf z} + \delta {\bf B}; \quad {\bf E} = \delta {\bf E}. \end{eqnarray*}$$

We seek plane wave solutions ∝exp [i(ωt − k · x)], where ω is the frequency, and k = (k_x, 0, k_z) the wave number. Some algebra allows us to express the linearized system in terms of the background protons, electrons, and PUI matrices M_s, M_e, and M_p through the electric field fluctuations and respective velocity fluctuations,

$$\begin{eqnarray} && \delta {\bf E} = B_0 {\bf M}_s \delta {\bf u}_s; \qquad {\rm (protons)} \nonumber \\ && \delta {\bf E} = B_0 {\bf M}_e \delta {\bf u}_e; \qquad {\rm (electrons)} \\ && \delta {\bf E} = B_0 {\bf M}_p \delta {\bf U}_p, \qquad {\rm (PUIs)} \nonumber \end{eqnarray} \tag{ 69 }$$

where

$$\begin{eqnarray} {\bf M}_s = \left(\begin{array}{@{}ccc@{}}-i\frac{C_s^2 k_x^2}{\omega \Omega _p}+i\frac{\omega }{\Omega _p}, & -1, & -i\frac{C_s^2 k_x k_z}{\omega \Omega _p} \\ +1, & +i\frac{\omega }{\Omega _p}, & 0 \\ -i\frac{C_s^2 k_x k_z}{\omega \Omega _p}, & 0, & -i\frac{C_s^2 k_z^2}{\omega \Omega _p}+i\frac{\omega }{\Omega _p} \end{array} \right); \quad \end{eqnarray} \tag{ 70 }$$

$$\begin{eqnarray} {\bf M}_e = \left(\begin{array}{@{}ccc@{}}+i\frac{C_e^2 k_x^2}{\omega \Omega _p}-i\xi \frac{\omega }{\Omega _p}, & -1, & +i\frac{C_e^2 k_x k_z}{\omega \Omega _p} \\ +1, & -i\xi \frac{\omega }{\Omega _p}, & 0 \\ +i\frac{C_e^2 k_x k_z}{\omega \Omega _p}, & 0, & +i\frac{C_e^2 k_z^2}{\omega \Omega _p}-i\xi \frac{\omega }{\Omega _p} \end{array} \right); \quad\quad \end{eqnarray} \tag{ 71 }$$

$$\begin{eqnarray} {\bf M}_p = \left(\begin{array}{@{}ccc@{}}-i \dsty\frac{C_p^2}{\Omega _p} \dsty\frac{k_x^2}{\omega -\dsty\frac{i}{3}{\bf k} \cdot {\bf K} \cdot {\bf k} }+i\dsty\frac{\omega }{\Omega _p}, & -1, & -i\dsty\frac{C_p^2}{\Omega _p} \dsty\frac{k_x k_z}{\omega -\dsty\frac{i}{3}{\bf k} \cdot {\bf K} \cdot {\bf k} } \\ +1, & +i\dsty\frac{\omega }{\Omega _p}, & 0 \\ -i\dsty\frac{C_p^2}{\Omega _p}\dsty\frac{k_x k_z}{\omega -\dsty\frac{i}{3}{\bf k} \cdot {\bf K} \cdot {\bf k} }, & 0, & -i\dsty\frac{C_p^2}{\Omega _p}\dsty\frac{k_z^2}{\omega -\dsty\frac{i}{3}{\bf k} \cdot {\bf K} \cdot {\bf k} }+i\dsty\frac{\omega }{\Omega _p} \end{array} \right). \nonumber\\ \end{eqnarray} \tag{ 72 }$$

The final relation needed to complete the derivation of the dispersion relation follows from the definition of the current,

$$\begin{equation} \alpha \delta {\bf u}_s + (1-\alpha)\delta {\bf U}_p -\delta {\bf u}_e = \frac{1}{B_0}{\bf M}_E \delta {\bf E}, \end{equation} \tag{ 73 }$$

where

$$\begin{equation} {\bf M}_E = \frac{i V_A^2}{\omega \Omega _p}\left(\begin{array}{ccc}k_z^2, & 0, & -k_x k_z\\ 0, & k^2, & 0 \\ -k_x k_z, & 0, & k_x^2 \end{array} \right). \end{equation} \tag{ 74 }$$

In Equations (69)–(74), we have introduced the proton gyrofrequency Ω_p = eB₀/m_p, the mass ratio ξ ≡ m_e/m_p, the electron number density n₀, the density ratio or concentration α ≡ n_s0/n₀ and thus n_p0 = n₀(1 − α). The square of the Alfvén speed is defined with the electron number density n₀ as

$$\begin{equation} V_A^2 = \frac{B_0^2}{\mu _0 n_0 m_p}, \end{equation} \tag{ 75 }$$

and the thermal proton, electron, and PUI sound speeds are defined as

$$\begin{eqnarray} && C_s^2 = \frac{\gamma _s P_{s0}}{\rho _{s0}} = \frac{\gamma _s P_{s0}}{n_{s0} m_p}; \nonumber \\ && C_e^2 = \frac{m_e}{m_p}\frac{\gamma _e P_{e0}}{\rho _{e0}} = \frac{\gamma _e P_{e0}}{n_0 m_p}; \\ && C_p^2 = \frac{\gamma _p P_{p0}}{\rho _{p0}} = \frac{\gamma _p P_{p0}}{n_{p0} m_p}, \nonumber \end{eqnarray} \tag{ 76 }$$

and γ_e/s/p refers to the appropriate adiabatic index. Note that the definition of the electron sound speed contains the ratio m_e/m_p. Finally, note that the scalar k · K · k is present in (72). The diffusion tensor is assumed to be of a simple diagonal form (i.e., we do not include the off-diagonal terms associated with drift and curvature—see the discussion in Zank 2013) and we specify

$$\begin{eqnarray} {\bf K} = \left(\begin{array}{ccc}\kappa _\perp & 0 & 0 \\ 0 & \kappa _\perp & 0 \\ 0 & 0 & \kappa _\parallel \end{array}\right); \end{eqnarray} \tag{ 77 }$$

$$\begin{eqnarray} &&{\bf k} \cdot {\bf K} \cdot {\bf k} = k^2 (\kappa _\perp \sin ^2\theta + \kappa _\parallel \cos ^2\theta). \end{eqnarray} \tag{ 78 }$$

Recall that the formal definition of the heat conduction spatial diffusion coefficient is (see the definitions (25) and (26))

$$\begin{eqnarray*} K_{ij} &\equiv & \int {c^{\prime }}^2 \kappa _{im} \frac{\partial f_0}{\partial x_m} {c^{\prime }}^2 d c^{\prime } / \frac{\partial }{\partial x_j} \int {c^{\prime }}^2 f_0 {c^{\prime }}^2 d c^{\prime } \nonumber\\ &\equiv & \frac{\langle c^2 \tau _s \rangle }{3} \simeq \frac{\langle c^2 \rangle \langle \tau _s \rangle }{3}, \end{eqnarray*}$$

where the ensemble average is taken over characteristic PUI speeds. Accordingly, if we crudely approximate the scattering time by the inverse gyrofrequency, we have κ_∥ ≡ 〈c²〉/3Ω_p, where 〈c²〉 is an ensemble averaged characteristic speed of the PUIs. We therefore introduce the definitions

$$\begin{equation} \kappa _\perp = \chi \frac{1}{3\Omega _p} \langle c^2 \rangle, \quad \kappa _\parallel = \frac{1}{3\Omega _p}\langle c^2 \rangle, \end{equation} \tag{ 79 }$$

and χ = 0.01. In estimating the diffusion coefficients (79), we need to choose a typical PUI speed for the system of interest weighted appropriately according to the definition of K, e.g., in the supersonic solar wind, one might choose 〈c²〉^1/2 to be some fraction of 400 km s⁻¹. Although we assume that the scattering time can be approximated by the PUI inverse gyrofrequency, the real scattering time may well be much slower. Our choice of χ = 0.01 is motivated by models of the perpendicular diffusion coefficient derived from nonlinear guiding center theory (Matthaeus et al. 2003), which, when evaluated for explicit turbulence models, yields a ratio of ∼0.01 between the perpendicular diffusion and parallel coefficients (Zank et al. 2004).

From the expression for the current (73), and relations (69) et seq., the determinant is

$$\begin{eqnarray*} &&\left| \alpha {\bf M}_s^{-1} + (1-\alpha) {\bf M}_p^{-1} - {\bf M}_e^{-1} - {\bf M}_E \right| = 0, \end{eqnarray*}$$

from which the dispersion relation can be expressed as

$$\begin{equation} \left| \alpha {\bf M}_p + (1-\alpha) {\bf M}_s - {\bf M}_p\cdot \left({\bf M}_e^{-1} + {\bf M}_E\right)\cdot {\bf M}_s \right| = 0. \end{equation} \tag{ 80 }$$

The dispersion relation (80) is formidably complicated, containing as a subset the full two-fluid model (Stringer 1963; Sahraoui et al. 2012 and references therein). Nonetheless, neglecting electron inertia yields a somewhat "tractable" 11th-order polynomial dispersion relation, whereas in the limit of small PUI abundance α → 1, electron inertia effects can be retained. These various limits are discussed below.

3.1. Dispersion Relation in the Absence of Electron Inertia

The neglect of electron inertia in (80) yields an 11th-order polynomial dispersion relation,

$$\begin{equation} \sum _{n=0}^{11} A_n \left(\frac{\omega }{\Omega _p}\right)^n=0, \end{equation} \tag{ 81 }$$

where the coefficients are given in Appendix B. The entropy mode ω = 0 is already eliminated from the dispersion relation (81). The dispersion relation is written in a form that facilitates expression in dimensionless units, e.g., $\widetilde{\omega }=\omega /\Omega _p$, $\widetilde{k}=k d_i = k V_A/\Omega _p$ and $\widetilde{{\bf k} \cdot {\bf K} \cdot {\bf k} }={\bf k} \cdot {\bf K} \cdot {\bf k} /\Omega _p$.

3.1.1. Parallel Propagation, θ = 0°

On setting θ = 0°, the dispersion relation (81) can be factored into three separate dispersion relations. The first corresponds to dispersive Alfvén waves, i.e., ion-cyclotron or whistler modes satisfying

$$\begin{eqnarray} && \omega ^4-(2k^2V_A^2+\frac{k^4V_A^4}{\Omega _p^2})\omega ^2+k^4V_A^4 = \left(\omega ^2+\frac{k^2V_A^2}{\Omega _p}\omega -k^2V_A^2 \right) \nonumber\\ && \quad \times \left(\omega ^2-\frac{k^2V_A^2}{\Omega _p}\omega -k^2V_A^2 \right) = 0, \end{eqnarray} \tag{ 82 }$$

with solutions

$$\begin{eqnarray} && \omega = \pm \frac{k^2 V_A^2}{2\Omega _p} + kV_A\sqrt{1+\left(\frac{k V_A}{2\Omega _p}\right)^2} \quad ({\rm yields} \quad \omega >0); \quad \end{eqnarray} \tag{ 83 }$$

$$\begin{eqnarray} && \omega = \pm \frac{k^2 V_A^2}{2\Omega _p} - kV_A\sqrt{1+\left(\frac{k V_A}{2\Omega _p}\right)^2} \quad ({\rm yields} \quad \omega <0). \quad \end{eqnarray} \tag{ 84 }$$

The ion-cyclotron mode satisfies ω → ±Ω_p in the limit that k → ∞, whereas the whistler mode satisfies ω → ±∞ in the same limit (since electron inertia is neglected).

The second factorization corresponds to a hybrid proton–PUI mode

$$\begin{equation} \omega = \pm \Omega _p. \end{equation} \tag{ 85 }$$

The remainder of the polynomial dispersion relation couples the proton and PUI sound modes, and is given by

$$\begin{eqnarray} && \omega ^5-\frac{i}{3}{\bf k} \cdot {\bf K} \cdot {\bf k} \omega ^4 - \left(C_e^2+C_s^2+C_p^2 \right)k^2\omega ^3 \nonumber\\ && \quad +\frac{i}{3}{\bf k} \cdot {\bf K} \cdot {\bf k} \left(C_e^2+C_s^2 \right)k^2\omega ^2 \nonumber \\ && \quad + \left((1-\alpha)C_e^2 C_s^2 + C_p^2 \left(C_s^2+\alpha C_e^2 \right) \right) k^4\omega \nonumber\\ && \quad -\frac{i}{3}{\bf k} \cdot {\bf K} \cdot {\bf k} (1-\alpha)k^4 C_e^2 C_s^2 = 0. \end{eqnarray} \tag{ 86 }$$

For the cases that we consider, α is quite close to 1, typically different by some 10%–20%. The dispersion relation is therefore usefully rewritten as

$$\begin{eqnarray} && \left(\omega ^4-\frac{i}{3}{\bf k} \cdot {\bf K} \cdot {\bf k} \omega ^3- \left(C_e^2+C_s^2+C_p^2 \right)k^2\omega ^2\right. \nonumber\\ &&\left. \quad\quad +\frac{i}{3}{\bf k} \cdot {\bf K} \cdot {\bf k} \left(C_e^2+C_s^2 \right)k^2\omega + C_p^2 \left(C_e^2+C_s^2 \right)k^4 \right) \omega \nonumber \\ && \quad {} = (1-\alpha)k^4C_e^2 C_s^2 \left(\frac{i}{3}{\bf k} \cdot {\bf K} \cdot {\bf k} - \omega \right). \end{eqnarray} \tag{ 87 }$$

It is easily seen that setting α = 1 yields as solutions to (87), the

$$\begin{eqnarray} &&{\rm entropy\, mode}: \quad \omega = 0; \end{eqnarray} \tag{ 88 }$$

$$\begin{eqnarray} &&{\rm two-fluid\, sound\, mode}: \quad \omega = \pm k \sqrt{C_s^2+C_e^2}; \end{eqnarray} \tag{ 89 }$$

$$\begin{eqnarray} {\rm PUI\, sound\, mode}: \quad \omega &=& \pm \sqrt{k^2 C_p^2 - \left(\frac{{\bf k} \cdot {\bf K} \cdot {\bf k} }{6}\right)^2} \nonumber\\ && + \frac{i}{6}{\bf k} \cdot {\bf K} \cdot {\bf k}. \end{eqnarray} \tag{ 90 }$$

In the absence of spatial diffusion, (90) corresponds to a simple PUI sound mode ω = ±C_pk. The second term in (90) is imaginary with positive sign and therefore describes how PUIs also damp the PUI sound mode. For α ≠ 1, the modes (88)–(90) are coupled.

Further insight into (86) is obtained by introducing the phase speed V_p ≡ ω/k. The diffusion term k · K · k can be expressed as κ_∥k² for parallel wave propagation. For a characteristic phase speed V₀, we may exploit the length scale introduced by the diffusion coefficient κ_∥ to consider long wavelength and short wavelength modes, using either

$$\begin{equation} \quad \frac{\kappa _{\parallel } k }{V_0} \ll 1 \quad {\rm or} \quad \frac{\kappa _{\parallel } k }{V_0} \gg 1, \end{equation} \tag{ 91 }$$

respectively. On using $\bar{V}_p \equiv V_p /V_0$, $\bar{C}_{s,e,p} \equiv C_{s,e,p} /V_0$, the normalized form of (86) becomes

$$\begin{eqnarray} && \bar{V}_p^5 - \frac{i}{3} \frac{\kappa _{\parallel } k }{V_0} \bar{V}_p^4 - \left(\bar{C}_e^2 + \bar{C}_s^2 + \bar{C}_p^2 \right) \bar{V}_p^3 + \frac{i}{3} \left(\bar{C}_e^2 + \bar{C}_s^2 \right) \frac{\kappa _{\parallel } k }{V_0} \bar{V}_p^2 \nonumber \\ && \quad\quad + \left((1 - \alpha) \bar{C}_e^2 \bar{C}_s^2 + \bar{C}_p^2 \left(\bar{C}_s^2 + \alpha C_{e}^2 \right) \right) \bar{V}_p \nonumber\\ && \quad\quad - \frac{i}{3} (1 - \alpha) \bar{C}_e^2 \bar{C}_s^2 \frac{\kappa _{\parallel } k }{V_0} = 0. \end{eqnarray} \tag{ 92 }$$

We can solve (92) in the long wavelength limit by setting ε ≡ O(κ_∥k/V₀) ≪ 1 and using $\bar{V}_p = \bar{V}_{p0} + \varepsilon \bar{V}_{p1}$. Similarly, we can solve (92) in the short wavelength limit using ε ≡ O((κ_∥k/V₀)⁻¹) ≪ 1 and using $\bar{V}_p = \bar{V}_{p0} + \varepsilon \bar{V}_{p1}$. For the long wavelength limit, we obtain five solutions, one of which is a purely damped mode

$$\begin{eqnarray} \omega &\simeq & \frac{i}{3} \kappa _{\parallel } k^2 \frac{ (1 - \alpha) C_e^2 C_s^2 }{(1 - \alpha) C_e^2 C_s^2 + C_p^2 (C_s^2 + \alpha C_e^2) } \nonumber \\ &\simeq & \frac{i}{3} \kappa _{\parallel } k^2 \frac{ (1 - \alpha) C_e^2 C_s^2 }{C_p^2 (C_s^2 + C_e^2) }, \end{eqnarray} \tag{ 93 }$$

with damping rate proportional to the PUI concentration. Another solution is the PUI acoustic mode that is coupled weakly to the thermal plasma and is diffusively damped,

$$\begin{eqnarray} \omega &\simeq & \pm C_p k \left(1+(1-\alpha)\frac{C_e^2(C_p^2-C_s^2)}{2C_p^2(C_p^2-C_s^2-C_e^2)} \right)\qquad \nonumber\\ && + \frac{i}{6}\kappa _\parallel k^2 \left(1-(1-\alpha)\frac{C_e^2(C_p^2-C_s^2)^2+C_s^2 C_e^4}{C_p^2(C_p^2-C_s^2-C_e^2)^2} \right).\qquad \end{eqnarray} \tag{ 94 }$$

The final long wavelength limit mode admitted by (92) is a two-fluid acoustic mode that is coupled weakly to the PUIs and is weakly diffusively damped, with

$$\begin{eqnarray} \omega &\simeq & \pm (C_s^2+C_e^2)^{1/2} k\left(1-\frac{(1-\alpha)C_e^2(C_p^2-C_s^2)}{2(C_s^2+C_e^2)(C_p^2-C_s^2-C_e^2)}\right) \nonumber\\ && + \frac{i}{6}\kappa _\parallel k^2 \frac{(1-\alpha)C_p^2 C_e^4}{(C_s^2+C_e^2)(C_p^2-C_s^2-C_e^2)^2}. \end{eqnarray} \tag{ 95 }$$

The converse short wavelength limit shows that there is no PUI mode. Only the two-fluid acoustic mode persists and it is very weakly coupled to the PUIs via the concentration (1 − α) and is damped by PUIs. Thus, we have

$$\begin{eqnarray} \omega &\simeq & \pm (C_s^2+C_e^2)^{1/2} k \left(1-(1-\alpha)\frac{C_s^2C_e^2}{2(C_s^2+C_e^2)^2} \right) \nonumber\\ && + \frac{3i}{2\kappa _\parallel }(1-\alpha)\frac{C_p^2C_e^4}{(C_s^2+C_e^2)^2}. \end{eqnarray} \tag{ 96 }$$

A second slowly propagating mixed thermal acoustic mode exists, which is also damped by PUIs,

$$\begin{equation} \omega \simeq \pm \sqrt{1-\alpha }\frac{C_s C_e}{\sqrt{C_s^2+C_e^2}}k + \frac{3i}{2\kappa _\parallel }C_p^2\left(\! 1-(1-\alpha)\frac{C_e^4}{(C_s^2+C_e^2)^2} \!\right). \end{equation} \tag{ 97 }$$

Finally, in the short wavelength limit, (92) admits a purely damped non-propagating mode

$$\begin{equation} \omega \simeq \frac{3i}{2\kappa _\parallel }C_p^2 \left(1-(1-\alpha)\frac{C_e^2(C_p^2-C_s^2)}{C_p^2(C_s^2+C_e^2)} \right). \end{equation} \tag{ 98 }$$

3.1.2. Perpendicular Propagation, θ = 90°

The reduced dispersion relation for perpendicularly propagating waves is a fifth-order polynomial,

$$\begin{eqnarray} && \omega ^5-\frac{i}{3}{\bf k} \cdot {\bf K} \cdot {\bf k} \omega ^4 -\omega ^3\left [k^2 \left(C_e^2+C_s^2+C_p^2 \right)\right. \nonumber\\ &&\left. \quad +\Omega _p^2\right ] +\frac{i}{3}{\bf k} \cdot {\bf K} \cdot {\bf k} \omega ^2 \left [ k^2 \left(C_e^2+C_s^2+V_A^2 \right)+\Omega _p^2\right ] \nonumber \\ && \quad +\omega \left \lbrace k^4\left [ (1-\alpha)C_s^2 \left(C_e^2+V_A^2 \right)+C_p^2 \left(C_s^2+\alpha C_e^2+\alpha V_A^2 \right)\right ]\right. \nonumber\\ &&\left. \quad +k^2\Omega _p^2\left [ (1-\alpha)C_p^2+\alpha C_s^2+C_e^2+V_A^2\right ] \right \rbrace \nonumber \\ && \quad -\frac{i}{3}{\bf k} \cdot {\bf K} \cdot {\bf k} \left [ k^4(1-\alpha)C_s^2 \left(C_e^2+V_A^2 \right)\right. \nonumber\\ &&\left. \quad +k^2\Omega _p^2 \left(\alpha C_s^2+C_e^2+V_A^2 \right)\right ] = 0. \end{eqnarray} \tag{ 99 }$$

It is instructive to consider the limit α = 1. In this case, we recover the two-fluid fast mode

$$\begin{equation} \omega = \pm \sqrt{V_A^2+C_s^2+C_e^2} k, \end{equation} \tag{ 100 }$$

together with a third-order polynomial describing the separate propagation of PUI modes,

$$\begin{equation} \omega ^3-\frac{i}{3}{\bf k} \cdot {\bf K} \cdot {\bf k} \omega ^2-\omega \left(C_p^2 k^2+\Omega _p^2\right)+\frac{i}{3}{\bf k} \cdot {\bf K} \cdot {\bf k} \Omega _p^2 =0. \end{equation} \tag{ 101 }$$

In this case, the term k · K · k = κ_⊥k². To analyze (101), we again introduce the phase velocity V_p ≡ ω/k and normalize with a characteristic speed V₀, and consider the κ_⊥k/V₀ ≪ 1 and κ_⊥k/V₀ ≫ 1 limits. The long wavelength limit yields a purely damped PUI mode with

$$\begin{equation} \omega \simeq \frac{i}{3} \kappa _{\perp } k^2 \frac{\Omega _p^2}{C_p^2 k^2 + \Omega _p^2}, \end{equation} \tag{ 102 }$$

and a damped PUI-cyclotron wave

$$\begin{equation} \omega \simeq \pm \sqrt{ C_p^2 k^2 + \Omega _p^2} + \frac{i}{6} \kappa _{\perp } k^4 \frac{C_p^2}{C_p^2 k^2 + \Omega _p^2}. \end{equation} \tag{ 103 }$$

The short wavelength limit possesses two purely damped modes, one of which is

$$\begin{equation} \omega \simeq \frac{i}{3} \kappa _{\perp } k^2 \frac{\Omega _p^2}{C_p^2 k^2 + \Omega _p^2}, \end{equation} \tag{ 104 }$$

together with

$$\begin{eqnarray} \omega &\simeq & \pm k \sqrt{C_p^2-\left(\frac{\kappa _\perp k}{6}\right)^2 } + \frac{i}{6}\kappa _\perp k^2 \nonumber\\ & =& \pm \frac{i}{6}\kappa _\perp k^2 \sqrt{1-\left(\frac{6C_p}{\kappa _\perp k}\right)^2}+\frac{i}{6}\kappa _\perp k^2. \end{eqnarray} \tag{ 105 }$$

In the short wavelength limit, we may expand the last expression to obtain the two limiting purely imaginary cases,

$$\begin{equation} \omega \simeq +\frac{i}{3}\kappa _\perp k^2; \quad \omega \simeq +\frac{3i}{\kappa _\perp }C_p^2. \end{equation} \tag{ 106 }$$

Although algebraically tedious, the fifth-order dispersion relation (99) is similarly amenable to an analysis in the long and short wavelength regimes. This has the virtue of identifying the basic couplings of the thermal plasma and PUI components, but is no more revealing of the basic modes than the above analysis. We do not therefore reproduce the analysis for (99).

3.1.3. Oblique Propagation, θ ≠ 0°, 90°

The only limit of the full dispersion equation that is analytically tractable is α → 1. In this case, in the absence of electron inertia, the dispersion relation (81) factors into the standard two-fluid dispersion relation without electron inertia (Stringer 1963; Sahraoui et al. 2012),

$$\begin{eqnarray} && \frac{\omega ^6}{\Omega _p^6}-\frac{\omega ^4}{\Omega _p^4}\left [\! k^4\frac{V_A^4}{\Omega _p^4}\cos ^2\theta +k^2\frac{V_A^2}{\Omega _p^2}\left (\! 1+\cos ^2\theta +\frac{C_e^2}{V_A^2}+\frac{C_s^2}{V_A^2} \!\right) \!\right ] \nonumber \\ && \quad +\frac{\omega ^2}{\Omega _p^2}\cos ^2\theta \left [\! k^6\frac{V_A^6}{\Omega _p^6}\left (\!\frac{C_e^2}{V_A^2}+\frac{C_s^2}{V_A^2} \!\right) +k^4\frac{V_A^4}{\Omega _p^4}\left (\! 1+2\frac{C_e^2}{V_A^2}+2\frac{C_s^2}{V_A^2} \!\right) \!\right ] \nonumber\\ && \quad -k^6\frac{V_A^6}{\Omega _p^6}\cos ^4\theta \left (\frac{C_e^2}{V_A^2}+\frac{C_s^2}{V_A^2}\right)= 0, \end{eqnarray} \tag{ 107 }$$

and a PUI-only dispersion relation,

$$\begin{eqnarray} && \omega ^4-\frac{i}{3}{\bf k} \cdot {\bf K} \cdot {\bf k} \omega ^3- \left(C_p^2 k^2+\Omega _p^2 \right)\omega ^2+\frac{i}{3}{\bf k} \cdot {\bf K} \cdot {\bf k} \Omega _p^2\omega \nonumber\\ && \quad\quad +C_p^2\Omega _p^2 k^2\cos ^2\theta = 0. \end{eqnarray} \tag{ 108 }$$

The solutions to (107) are the usual dispersive two-fluid modes, viz., slow, fast, and Alfvén modes with the standard properties. The PUI dispersion relation can be solved in the long and short wavelength limits. In the long wavelength regime (using κ ≡ κ_∥cos ²θ + κ_⊥sin ²θ), we find two solutions for the arbitrary angle θ that yield positive real frequency,

$$\begin{eqnarray} \omega \simeq && \sqrt{\frac{1}{2} \left(C_p^2 k^2 +\Omega _p^2 \right) \pm \frac{1}{2}\sqrt{\left(C_p^2 k^2 +\Omega _p^2 \right)^2 -4C_p^2\Omega _p^2 k^2\cos ^2\theta } } \nonumber \\ && +\frac{i}{12}\kappa k^2 \left(1\pm \frac{C_p^2k^2-\Omega _p^2}{\sqrt{ \left(C_p^2 k^2 +\Omega _p^2 \right)^2 -4C_p^2\Omega _p^2 k^2\cos ^2\theta }} \right). \nonumber\\ \end{eqnarray} \tag{ 109 }$$

Note that the two plus signs and the two minus signs act in concert. The two quasi-parallel propagating modes in expression (109) are expressed in a more revealing form as

$$\begin{eqnarray} \omega &\simeq & \pm C_p k \left(\frac{C_p^2 k^2 - \Omega _p^2 \cos ^2 \theta }{C_p^2 k^2 - \Omega _p^2 } \right)^{1/2} \nonumber\\ && + \frac{i}{6} \kappa k^2 \frac{C_p^4 k^4 - C_p^2 \Omega _p^2 k^2 (2 - \sin ^2 \theta) + \Omega _p^4 }{C_p^4 k^4 - 2C_p^2 \Omega _p^2 k^2 \cos ^2 \theta + \Omega _p^4 }; \nonumber \\ \omega &\simeq & \pm \Omega _p \left(\frac{C_p^2 k^2 \cos ^2 \theta - \Omega _p^2 }{C_p^2 k^2 - \Omega _p^2 } \right)^{1/2} \nonumber\\ && + \frac{i}{6} \kappa k^2 \frac{ C_p^2 \Omega _p^2 k^2 \sin ^2 \theta }{C_p^4 k^4 - 2C_p^2 \Omega _p^2 k^2 \cos ^2 \theta + \Omega _p^4 }, \end{eqnarray} \tag{ 110 }$$

and the two quasi-perpendicular propagating modes as

$$\begin{eqnarray} \omega &\simeq & \pm \left(C_p^2 k^2 + \Omega _p^2 \right)^{1/2} \left(1 - \frac{C_p^2 \Omega _p^2 k^2 \cos ^2 \theta }{(C_p^2 k^2 + \Omega _p^2)^2 } \right)^{1/2} \nonumber\\ &&+ \frac{i}{6} \kappa C_p^2 k^4 \frac{C_p^2 k^2 + \Omega _p^2 \sin ^2 \theta }{C_p^4 k^4 + 2C_p^2 \Omega _p^2 k^2 \sin ^2 \theta + \Omega _p^4}; \nonumber \\ \omega &\simeq & \pm \frac{C_p \Omega _p k \cos \theta }{(C_p^2 k^2 + \Omega _p^2)^{1/2} } + \frac{i}{6} \kappa k^2 \Omega _p^2 \nonumber\\ && \times \frac{C_p^2 k^2 \sin ^2 \theta + \Omega _p^2 }{C_p^4 k^4 + 2C_p^2 \Omega _p^2 k^2 \sin ^2 \theta + \Omega _p^4}. \end{eqnarray} \tag{ 111 }$$

Finally, the short wavelength limit yields four damped modes, with

$$\begin{eqnarray} \omega &\simeq & \pm k\sqrt{C_p^2-\left(\frac{\kappa k}{6}\right)^2} +\frac{i}{6}\kappa k^2 \nonumber\\ &=& \pm \frac{i}{6}\kappa k^2 \sqrt{1-\left(\frac{6C_p}{\kappa k} \right)^2} +\frac{i}{6}\kappa k^2; \end{eqnarray} \tag{ 112 }$$

$$\begin{eqnarray} &&\omega \simeq \frac{3i}{\kappa } C_p^2 \cos ^2 \theta. \end{eqnarray} \tag{ 113 }$$

At sufficiently small scales (112) is purely damped and becomes

$$\begin{equation} \omega \simeq +\frac{i}{3}\kappa k^2; \quad \omega \simeq +\frac{3i}{\kappa } C_p^2. \end{equation} \tag{ 114 }$$

3.1.4. Inclusion of Electron Inertia

The inclusion of electron inertial terms renders the dispersion relation unreasonably large to express as a polynomial, and so we directly solve the linearized system of equations numerically. However, it is illuminating to consider, as before, the limit α → 1. In this case, the entire dispersion relation separates into a two-fluid dispersion relation (Stringer 1963; Sahraoui et al. 2012) and a separate PUI dispersion relation. The two-fluid dispersion relation is the well-known sixth-order polynomial,⁴ given by

$$\begin{equation} A_6\frac{\omega ^6}{\Omega _p^6}-A_4\frac{\omega ^4}{\Omega _p^4}+A_2\frac{\omega ^2}{\Omega _p^2}-A_0=0, \end{equation} \tag{ 115 }$$

where the coefficients are

$$\begin{eqnarray} && A_6 = (1+\xi)\left(1+\xi +\xi \frac{k^2 V_A^2}{\Omega _p^2}\right)^2; \nonumber \\ && A_4 = \xi ^2\frac{k^6 V_A^4 V_S^2}{\Omega _p^6}+\frac{k^4 V_A^4}{\Omega _p^4} \nonumber\\ && \qquad\times \left[ \cos ^2\theta (1+\xi ^3)+\xi (1+\xi)\left(1+\frac{2V_S^2}{V_A^2}\right)\right] \nonumber \\ && \qquad +\frac{k^2 V_A^2}{\Omega _p^2}(1+\xi)^2\left[1+\cos ^2\theta +\frac{V_S^2}{V_A^2}\right]; \nonumber \\ && A_2 = \cos ^2\theta \left[ \frac{k^4 V_A^4}{\Omega _p^4}(1+\xi)\left(1+\frac{2V_S^2}{V_A^2}\right) +\frac{k^6V_A^4V_S^2}{\Omega _p^6}(1+\xi ^2)\right];\nonumber \\ && A_0 = \frac{k^6V_A^4V_S^2}{\Omega _p^6}\cos ^4\theta. \end{eqnarray} \tag{ 116 }$$

In (116), we have introduced the combined thermal plasma sound speed $V_S ^2 = C_s^2 + C_e^2$.

The PUI dispersion relation may be expressed as

$$\begin{eqnarray} && \omega ^4-\frac{i}{3}{\bf k} \cdot {\bf K} \cdot {\bf k} \omega ^3- \left(C_p^2 k^2+\Omega _p^2 \right)\omega ^2+\frac{i}{3}{\bf k} \cdot {\bf K} \cdot {\bf k} \Omega _p^2\omega \nonumber\\ && \quad +C_p^2\Omega _p^2 k^2\cos ^2\theta = 0, \end{eqnarray} \tag{ 117 }$$

which, not surprisingly, is identical to (108). In the limit that α = 1, electron inertia therefore modifies only the two-fluid modes. For a non-vanishing pickup ion abundance, all modes are of course coupled and electron inertia modifies all waves.

For parallel wave propagation, solutions of the two-fluid dispersion relation (115)–(116) can be obtained readily. The dispersive Alfvén wave solutions (the ion-cyclotron and whistler modes) are given by

$$\begin{eqnarray} && \omega = \frac{\pm \frac{k^2 V_A^2}{2\Omega _p}(1-\xi) +kV_A\sqrt{1+\xi +\left(\frac{k V_A}{2\Omega _p}\right)^2 (1+\xi)^2}}{1+\xi \left(1+\frac{k^2V_A^2}{\Omega _p^2}\right)}, \quad\nonumber\\ && \quad\quad ({\rm yields} \quad \omega >0) \end{eqnarray} \tag{ 118 }$$

$$\begin{eqnarray} && \omega = \frac{\pm \frac{k^2 V_A^2}{2\Omega _p}(1-\xi) -kV_A\sqrt{1+\xi +\left(\frac{k V_A}{2\Omega _p}\right)^2 (1+\xi)^2}}{1+\xi \left(1+\frac{k^2V_A^2}{\Omega _p^2}\right)} \quad\nonumber\\ && \quad\quad ({\rm yields} \quad \omega <0). \end{eqnarray} \tag{ 119 }$$

Now, because of the inclusion of electron inertia, the whistler mode frequency is bounded by the electron cyclotron frequency in the limit k → ∞ and ω → ±(Ω_p/ξ) = ±Ω_e. For the ion-cyclotron mode, the limit k → ∞ yields the same result when electron inertia terms are neglected, i.e., ω → ±Ω_p.

Similarly, the two-fluid sound mode is a solution of (115)–(116) with

$$\begin{equation} \omega =\pm \frac{k V_S}{\sqrt{1+\xi }}. \end{equation} \tag{ 120 }$$

As in the case with no electron inertia, the hybrid proton-pickup ion mode ω = ±Ω_p and the ion-cyclotron and whistler modes (118) and (119) are also solutions of the full linearized coupled system (i.e., without taking the limit α → 1). The limit α → 1 is useful primarily to distinguish the two-fluid sound mode (120) and the pickup ion sound mode (90), which are otherwise coupled.

For perpendicular wave propagation, the fast two-fluid magnetosonic mode is modified by electron inertial effects to read

$$\begin{equation} \omega =\pm k \sqrt{\frac{V_A^2+V_S^2+\frac{\xi }{1+\xi }\frac{k^2V_A^2V_S^2}{\Omega _p^2}}{1+\xi +\xi \frac{k^2 V_A^2}{\Omega _p^2}}}. \end{equation} \tag{ 121 }$$

In the limit k → ∞, the frequency is $\omega =\pm k V_S/\sqrt{1+\xi }$.

Finally, for obliquely propagating modes, an additional horizontal "resonant" frequency asymptote is introduced, given by ω = Ω_ecos θ, as is clearly seen in numerical solutions (it is harder to see analytically). It is interesting to note that for purely parallel propagation, it is the fast (whistler) mode that has this limit. However, for oblique propagation (in which case the modes never intersect), the fast mode frequency continues to increase and it is the Alfvén mode that has this limit. For normal propagation angles, this is relevant only at very small scales, but for highly oblique angles (see, for example, our later Figure 1 for θ = 88°), this effect becomes more apparent at larger scales.

Zoom In Zoom Out Reset image size

4.1. Identification of Modes in the Numerical Solutions

To properly investigate the dispersion relation in full generality appropriate to the three PUI mediated regions of interest (the outer heliosphere, the IHS, and the VLISM), we need to obtain numerical solutions. The previous estimated solutions provide useful guidance to identify modes at the largest and smallest scales and weak concentrations. Considering only positive frequencies at long wavelengths (conveniently defined as in the vicinity of kV_A/Ω_p = 0.01), we introduce the following colors to correspond to the following (approximate) designation of wave modes. (1) The slow (red) and fast (blue) two-fluid magnetosonic modes:

$$\begin{equation} \omega _{s,f} = k\sqrt{\frac{1}{2} \left(V_S^2+V_A^2 \right) \pm \frac{1}{2}\sqrt{ \left(V_S^2+V_A^2 \right)^2-4V_S^2 V_A^2 \cos ^2\theta }}; \end{equation} \tag{ 122 }$$

(2) the Alfvén mode (green):

$$\begin{equation} \omega _A= k V_A \cos \theta; \end{equation} \tag{ 123 }$$

(3) the slow pickup ion mode (cyan):

$$\begin{equation} \omega _{sp}= k C_p \cos \theta; \end{equation} \tag{ 124 }$$

and the fast pickup ion mode (black):

$$\begin{equation} \omega _{fp}=\Omega _p. \end{equation} \tag{ 125 }$$

For the cases of interest here, the PUI sound speed C_p is larger than both the Alfvén speed V_A and combined thermal plasma sound speed $V_S = \sqrt{C_s^2 + C_e^2}$. Thus, for small propagation angles, the order of the modes will always be ω_s < ω_A < ω_f < ω_sp < ω_fp, i.e., red, green, blue, cyan, and black. However, as θ increases, the frequency of the fast proton (blue) mode increases too with $\omega _f\rightarrow k\sqrt{\vphantom{A^A}\smash{{V_S^2+V_A^2}}}$ and the frequency of the slow pickup ion (cyan) mode decreases with ω_sp → 0. Therefore, beyond some critical angle θ_c, the fast proton mode will be faster than the slow pickup ion mode (i.e., the blue mode will find itself above the cyan mode). This will happen at the critical angle

$$\begin{equation} \cos ^2\theta _c = \frac{V_A^2}{C_p^2}\left(1+\frac{V_S^2}{V_A^2}-\frac{V_S^2}{C_p^2} \right). \end{equation} \tag{ 126 }$$

Furthermore, because ω_A < ω_sp for all angles θ beyond the critical angle θ_c the order of the modes will always be ω_s < ω_A < ω_sp < ω_f < ω_fp, i.e., red, green, cyan, blue, and black. Of course, this is only approximate since we are discussing the strict limits α → 1 and k · K · k → 0. Nevertheless, for our parameters at 10 AU, for example, the critical angle corresponds to θ_c = 76.71°, which is why it is difficult to correctly identify the cyan and blue modes at 70° in this case. The corresponding angles for the VLISM and the IHS are θ_c = 71.86° and θ_c = 79.14°, respectively.

In all three of the particular cases that we consider below, we solve numerically the full dispersion relation including the electron inertial terms with ξ = m_e/m_p = 1/1836.

The parallel and perpendicular diffusion coefficients are specified in (79) and are taken to be proportional to the square of the weighted characteristic PUI speed and inversely proportional to the gyrofrequency. The scattering time is almost certainly slower than this, but it is a simple initial estimate and the magnitude of the diffusion coefficients can be adjusted as necessary. The magnitude of the diffusion coefficients therefore depends on the PUI speed and the strength of the magnetic field. As we illustrate below, the PUI diffusion coefficient is smallest in the VLISM, somewhat larger in the IHS, and largest in the supersonic solar wind of the outer heliosphere. For the purpose of ordering the discussion, we begin with the VLISM and conclude with the outer heliosphere.

4.2. The PUI-mediated VLISM

As discussed in the Introduction, we assume a PUI density n_p = 5 × 10⁻⁵ cm⁻³ and an electron density n₀ = 0.08 cm⁻³, so that n_p/n₀ = 6.25 × 10⁻⁴. Thus, the parameter α = 1 − n_p/n₀ = 0.999375. We assume a VLISM magnetic field B₀ = 3 μG = 0.3 nT, which yields an Alfvén speed $V_A=B_0/\sqrt{\mu _0 n_0 m_p}=23.14$ km s⁻¹. For the reasons discussed in Zank et al. (2013), we take the proton temperature to be T_s = 15,000 K, the electron temperature T_e = T_s and the polytropic indices γ_s = γ_e = γ_p = 5/3. The proton and electron sound speeds are $C_s^2=\gamma _s k_B T_s /m_p$ and $C_e^2=\gamma _e k_B T_e /m_p$, which yields C_s = C_e = 14.37 km s⁻¹. To estimate the sound speed for the VLISM PUI gas, we assume that fast neutrals from the supersonic solar wind and hot neutrals from the IHS are the primary origin of PUIs in the VLISM. The speed of the fast and hot neutrals relative to the VLISM plasma ranges from some 700–350 km s⁻¹ or 200–100 km s⁻¹. The PUI distribution may be approximated crudely as a shell distribution (f_p = n_pδ(c − U₀)/4πc² in the plasma frame) with some form of weighted particle speed U₀, a nested shell with a set of different values of U₀, or possibly even a partially filled shell distribution ($f_p = 3n_p c^{-3/2} / 8\pi U_0^{3/2}$ in the plasma frame), with a weighted U₀, related to the radius of the shell and the maximum particle speed relative to the background flow. The pressure for either a shell or filled-shell (Vasyliunus–Siscoe) distribution is $P_p = \rho _p U_0^2 /3$ or $P_p = \rho _p U_0^2 /7$. The PUI sound speed may be estimated from $C_p=U_0 \sqrt{\gamma _p/a}$, where a varies from 3–7. Below, we use C_p = 97.59 km s⁻¹, which for a shell corresponds to U₀ ∼ 130 km s⁻¹ or for a filled-shell to U₀ ∼ 200 km s⁻¹. Precisely what the correct weighted value of U₀ should be is uncertain but this provides a reasonable estimate for the PUI sound speed that is clearly distinguished from the values of C_s and C_e. We therefore do not expect the basic dispersion curves to be very much different for other values of C_p.

To illustrate the sensitivity of the plasma wave properties to the assumed PUI spatial diffusion coefficient in the VLISM, we consider in detail two cases. (1) Assume that 〈c²〉^1/2 = 200 km s⁻¹, and (2) that 〈c²〉^1/2 = 50 km s⁻¹. Both values appear to be quite reasonable. Since Ω_p = eB₀/m_p = 2.87 × 10⁻² s⁻¹, we may estimate the diffusion tensor as κ_∥ = 〈c²〉/(3Ω_p) = 4.64 × 10¹¹ m²/s or κ_∥ = 2.9 × 10¹⁰ m²/s. The perpendicular diffusion is suppressed relative to the parallel diffusion coefficient, and is given by κ_⊥ = χκ_∥ with parameter χ = 0.01.

Numerical solutions of the dispersion relation, including electron inertia, for the parameters above with 〈c²〉^1/2 = 200 km s⁻¹ are shown in Figure 1. Three columns of panels are shown in Figure 1. The leftmost column shows numerical solutions of the standard two-fluid plasma model (Stringer 1963; Sahraoui et al. 2012), i.e., thermal electrons and protons only. The middle column shows plots of the real frequency as a function of wave number for the PUI-mediated three-fluid plasma. The rightmost column shows the corresponding imaginary part of the frequency as a function of wave number of the PUI-mediated VLISM. The rows correspond to specific wave propagation directions, from parallel (top row) to perpendicular (bottom), as listed in the figure caption.

The two-fluid dispersion relation is included for reference, so that the similarities (and differences) in the PUI-mediated model are apparent. The two-fluid model for wave propagation angles θ < 90° shows that three waves exist in the system: the fast (blue) and slow (red) magnetosonic modes and the Alfvén mode (green). These are clearly distinguishable for 0 < θ < 90°. For large k, the parallel Alfvén mode tends to the asymptote ω = Ω_p. The parallel fast mode tends to Ω_e. For oblique angles, the Alfvén mode (green) eventually resonates at the electron frequency Ω_ecos θ, and the slow mode (red) resonates at the proton gyrofrequency Ω_pcos θ. The two-fluid perpendicular case is degenerate with only the fast mode remaining.

The PUI-mediated model possesses a richer set of wave modes than the two-fluid case. Shown in the middle panels, the PUI-mediated dispersion relation now exhibits two further modes, a fast PUI (black) and slow PUI (cyan) mode. The PUI sound speed is larger than the thermal plasma sound speeds and so these dispersion curves lie above their thermal counterparts. As illustrated in the middle column, the blue, red, and green curves are virtually unchanged from their corresponding two-fluid counterparts. Not surprisingly, since the abundance parameter α ≃ 1 for the VLISM, the thermal plasma modes are essentially decoupled from the PUI population, at least for this value of the diffusion coefficient. For parallel propagation, the dispersive Alfvén waves (blue and green) (83), sometimes called the ion-cyclotron and whistler modes, are unchanged from the two-fluid case. The fast PUI mode (ω = ±Ω_p) is now present (black curve). Two acoustic modes, the slow PUI mode (cyan; Equation (94)), and the acoustic thermal mode (red; Equation (95)) exist at long wavelengths (small k). As illustrated in the top panel of the right column, and supported analytically by expressions (94) and (95), the parallel PUI sound mode is damped more strongly than the parallel thermal acoustic mode. The parallel propagating PUI acoustic mode exists only for k ⩽ Ω_p/V_A, after which it is completely damped away and no longer propagates. By contrast, at short wavelengths, the slow thermal acoustic mode (Equation (96)) is coupled very weakly to the PUIs, the propagation speed being ${\sim } \sqrt{C_s^2 + C_e^2}$, but is damped by PUIs as illustrated by the red curve in the top panel of the right column (see also Equation (96)).

With increasing propagation angles θ, the five wave modes remain conceptually similar in terms of the classification but now, as illustrated in the panels of the rightmost column, all oblique modes are damped by PUIs, including the Alfvén mode. Also, despite not being damped when propagating parallel to B₀, the fast thermal acoustic mode is much more strongly damped than the slow thermal acoustic mode, despite the latter experiencing damping at θ = 0°. The fast PUI mode (black) deviates from the ω = Ω_p line with increasing k, but becomes evanescent shortly after peaking in frequency. The imaginary frequency of course becomes very large for this mode. The frequencies of the slow thermal and PUI acoustic modes (red and cyan) tend to the proton gyrofrequency asymptote ω = Ω_pcos θ. The fast PUI acoustic modes intersects the fast thermal acoustic mode at certain propagation angles. This is illustrated in the third panel down of Figure 1. For θ = 30°, the slow PUI acoustic mode (cyan) intersects both the thermal fast mode (blue) and the Alfvén mode (green), and converges to the same asymptote as the thermal slow mode. For greater obliquities, the PUI slow mode (cyan) intersects only the Alfvén mode (green). Finally, for perpendicular propagation (last panel down), only the fast PUI and fast thermal acoustic modes remain, both experiencing damping. Above a certain value of k, the short wavelength fast PUI mode becomes evanescent. The numerical results presented here are consistent with the simpler analytic derivations above.

Plotted in Figure 2 are the dispersion curves for the case of 〈c²〉^1/2 = 50 km s⁻¹. All other VLISM parameters identical, including the PUI sound speed and the PUI abundance. The figure format is identical to that of Figure 1 in terms of the columns and relationship of the color of the curves to the different wave modes. However, the dispersion relation curves are now fundamentally different from the 〈c²〉^1/2 = 200 km s⁻¹ example. The smaller diffusion coefficient ensures a coupling between the thermal VLISM plasma and the PUIs.

Zoom In Zoom Out Reset image size

Figure 2 shows the full complement of the dispersion curves but the level of detail makes it difficult to discern the differences with those curves shown in Figure 1. Illustrated in Figure 3 are the θ = 30° cases for the larger (top panel, 〈c²〉^1/2 = 200 km s⁻¹) and smaller (bottom panel, 〈c²〉^1/2 = 50 km s⁻¹) diffusion coefficients. The upper panel for the large diffusion example shows that the thermal fast (blue) and slow (red) magnetosonic modes and the Alfvén mode (green) are fundamentally unchanged from the two-fluid waves, and the PUI slow curve (cyan) intersects the thermal fast mode curve. Besides the obvious difference in the PUI fast mode (black curve), the bottom panel shows that (1) the PUI slow mode curve (cyan) no longer intersects the thermal fast mode curve but is "hyperbolic" to it; (2) the PUI slow mode curve completely changes character beyond k ∼ Ω_p/V_A, and now resembles the thermal fast mode curve in that ω_r (and k) increase asymptotically to the electron gyrofrequency Ω_ecos θ; (3) the thermal fast mode and the Alfvén wave curves approach one another hyperbolically and do not intersect; and (4) the thermal fast mode frequency increases with increasing k, tending asymptotically toward Ω_ecos θ, but the Alfvén wave (green) curve is now asymptotic to the proton gyrofrequency Ω_pcos θ. The behavior of the thermal and Alfvén modes is fundamentally different from the two-fluid model when the PUI heat conduction diffusion coefficient is small enough to ensure coupling of the thermal and PUI plasmas.

Zoom In Zoom Out Reset image size

The differences between the large and small heat conduction diffusion cases for different obliquities are now readily determined from Figure 2. Observe for example that the PUI slow mode is hyperbolic to the Alfvén mode when θ = 88°.

Consider the high-frequency short wavelength behavior of the two examples shown in Figures 1 and 2. Using the same parameters, we replot the θ = 30° and 88° cases in Figure 4, extending kV_A/Ω_p to 10⁴. The leftmost column shows the two-fluid dispersion curves, the middle column the three-fluid model where 〈c²〉^1/2 = 200 km s⁻¹ is used to estimate the diffusion coefficient (i.e., the Figure 1 case), and the rightmost column illustrates the three-fluid VLISM model using 〈c²〉^1/2 = 50 km s⁻¹ (i.e., the Figure 2 case). For the two-fluid plasma model, θ = 30°, notice that the Alfvén mode (green) resonates at the electron cyclotron frequency asymptote Ω_ecos θ as k → ∞ i.e., at small spatial scales. By contrast, the fast/whistler mode dispersion curve (blue) briefly tracks the electron gyrofrequency Ω_ecos θ asymptote for a range of wave numbers before eventually continuing to increase with increasing wave number. For highly oblique propagation angles (e.g., θ = 88°), the whistler mode no longer tracks the asymptote Ω_ecos θ at all and simply continues to increase monotonically.

Zoom In Zoom Out Reset image size

As discussed above at length, for the large diffusion coefficient case (〈c²〉^1/2 = 200 km s⁻¹), the middle panel shows again that the thermal modes in the three-fluid model are not influenced by the PUIs (except for damping) and the three-fluid model behaves as a separate two-fluid and PUI gas. By contrast, the three-fluid model with low PUI diffusion (right panels) exhibits a strong coupling of the thermal and PUI gas, leading to the distinctly different linear wave characteristics discussed in the context of Figures 2 and 3. The small wavelength dispersion curves are similarly different. For θ = 30°, the slow PUI mode curve (cyan) tracks the Ω_ecos θ line for a range of k values before increasing again with increasing k. The fast thermal mode curve (blue) increases more slowly than the slow PUI mode curve until it approaches the electron gyrofrequency Ω_ecos θ asymptotically for propagation angles less than the critical angle θ_c. By contrast, for highly oblique wave propagation angles greater than θ_c, the slow PUI sound mode (cyan) approaches the electron gyrofrequency Ω_ecos θ asymptotically and the fast thermal and fast PUI modes now continue to increase beyond the electron gyrofrequency line. Finally, the Alfvén wave dispersion curve, unlike the two-fluid and larger diffusion cases, now approaches the proton gyrofrequency line Ω_pcos θ asymptotically rather than the electron gyrofrequency.

Before concluding this section, we further clarify briefly the nature of the interaction of several of the wave modes as revealed in the dispersion relation plots Figures 1, 2, and 3 for oblique propagation. In Figure 5, we consider again the 30° case. The top and bottom panels show the cases illustrated in Figures 1, 2, and 3 for the fast thermal and slow PUI modes. For the 〈c²〉^1/2 = 200 km s⁻¹ case, as discussed, the two modes preserve their distinct identities and intersect, i.e., the modes behave as an approximately non-interacting thermal two-fluid and a PUI gas. Reducing the diffusion coefficient by assuming 〈c²〉^1/2 = 50 km s⁻¹ (bottom panel) shows that the wave modes no longer preserve their distinct thermal or PUI character—at larger values of wave number k, the thermal fast mode wave behaves more like the PUI slow mode, and the PUI slow mode assumes characteristics of the thermal fast mode. The two modes no longer intersect, as discussed above, as they do in the case of the larger diffusion coefficient. The panels in between illustrate the effect of varying the abundance parameter α (second and third panels down). This variation is somewhat academic since e.g., a value of α = 0.97 requires a VLISM PUI number density of 2.4 × 10⁻³ cm⁻³, which is unreasonably large for the VLISM, but may be reasonable in other settings. The sensitivity of the wave solutions to the abundance α is revealed in panels 2 and 3 of Figure 5, which show that for a value α = 0.98, the thermal fast and PUI slow dispersion curves intersect, whereas for α = 0.97 they do not. The fourth panel down further illustrates the sensitivity of the dispersion curves to the assumed value of 〈c²〉^1/2—in this case 〈c²〉^1/2 = 100 km s⁻¹, and the curves resemble the 〈c²〉^1/2 = 200 km s⁻¹ example. Figure 5 clarifies that there are two prerequisites for the three-fluid PUI mediated model to behave as a separate thermal two-fluid and PUI charged gas: (1) the PUI number density must be sufficiently low (α must be sufficiently close to (1), and (2) the PUI heat conduction diffusive coefficient must be sufficiently large (large value of 〈c²〉^1/2).

Zoom In Zoom Out Reset image size

4.3. The PUI-mediated Subsonic Solar Wind or the Inner Heliosheath

As discussed in the Introduction, the IHS plasma is mediated in a complex way by PUIs whose origin is both the ISM directly (through charge exchange of interstellar neutrals flowing into the IHS from the local interstellar medium) and PUIs created in the supersonic solar wind and transmitted through the heliospheric termination shock into the IHS. The transmission process is complicated in that some PUIs are transmitted directly and others only after experiencing specular reflection at the shock. The PUI properties inferred theoretically for the IHS as a consequence of these processes (Zank et al. 1996a, 2010; Burrows et al. 2010) have been confirmed and further refined observationally (Richardson 2008; Richardson et al. 2008; Desai et al. 2012, 2014; Zirnstein et al. 2014). The thermal proton density is taken to be n_s = 5 × 10⁻³ cm⁻³ and the ratio of pickup ion to proton density n_p/n_s = 0.2. This yields the parameter α = 1/(1 + n_p/n_s) = 0.833 and the electron density n₀ = n_s(1 + n_p/n_s) = 6 × 10⁻³ cm⁻³. The magnetic field is taken to be B₀ = 0.1 nT, which yields an Alfvén speed V_A = 28.16 km s⁻¹. We assume a proton temperature T_s = 150,000 K, an electron temperature T_e = T_s and polytropic indices γ_s = γ_e = γ_p = 5/3. The proton and electron sound speeds are therefore C_s = C_e = 45.43 km s⁻¹. We assume a composite PUI temperature T_p ∼ 10⁷ K, which yields the PUI sound speed classically as $C_p^2=\gamma _p k_B T_p /m_p$, or C_p = 370.91 km s⁻¹. For the diffusion tensor κ_∥ = 〈c²〉/(3Ω_p), we assume a characteristic PUI speed of 〈c²〉^1/2 = 300 km s⁻¹ and for the perpendicular diffusion coefficient, κ_⊥ = χκ_∥, we assume χ = 0.01.

Despite the greater mobility of the PUIs in the IHS yielding a larger diffusion coefficient than in the VLISM, the smaller parameter α ensures a strong coupling of the thermal plasma and PUIs. The plots of Figure 6 depicting the real frequency as a function of k therefore resemble those of Figure 2 in many respects. However, for parallel propagation, the two-fluid curves and the thermal modes of the three-fluid model are almost identical. For parallel propagation, the thermal and PUI modes behave as though they are uncoupled (although PUIs damp the thermal modes). The PUI slow mode (cyan) is evanescent for a range of k values, existing for long wavelengths satisfying approximately k ⩽ Ω_p/V_A and short wavelengths k ⩾ 10Ω_p/V_A.

Zoom In Zoom Out Reset image size

The dispersion curves change quite dramatically for obliquely propagating modes, and the two-fluid and three-fluid dispersion curves are now quite different. Whereas the parallel PUI slow (cyan) and fast (black) modes intersected, the obliquely propagating modes no longer do so. Instead, as illustrated for θ = 30°, the intersection point becomes "hyperbolic" and the PUI slow and fast modes swap characteristics. Now it is the obliquely propagating PUI fast mode (black) that is evanescent in a range of k space. The PUI slow mode (cyan) becomes quite complicated for oblique propagation, assuming characteristics of not only the PUI fast mode but also the thermal fast magnetoacoustic mode at shorter wavelengths. The PUI slow mode no longer possesses a region that is evanescent.

As illustrated in the two-fluid column of dispersion curves, the behavior of the oblique thermal fast (blue), slow (red) magnetoacoustic and Alfvén (green) modes changes significantly in the presence of PUIs (middle column). For example, both the thermal slow mode and the Alfvén wave tend to the proton gyrofrequency at short wavelengths—the short wavelength behavior of the oblique Alfvén mode is quite different from that exhibited in the two-fluid model.

Finally, the only perpendicularly propagating modes are the PUI and thermal fast modes.

4.4. The PUI-mediated Supersonic Solar Wind Beyond ∼10 AU

The thermal supersonic solar wind beyond the ionization cavity (i.e., >∼10 AU) is dominated by the presence of a suprathermal PUI population comprising primarily protons. Besides protons (Gloeckler et al. 1993), heavier interstellar PUIs are present but in much smaller quantities, first observed in situ by Ulysses (Geiss et al. 1994; Gloeckler et al. 1994) and since by New Horizons (McComas et al. 2008; Randol et al. 2013). We consider solar wind parameters appropriate to 10 AU, assuming that the magnetic field B₀ = 0.5 nT, density n₀ = 0.08 cm⁻³, proton temperature T_s = 20,000 K and estimate that the pickup proton abundance is n_p/n_s = 0.05 (5%) so that α = 1/(1 + n_p/n_s) = 0.952. Heavy PUIs are neglected since their number density is extremely small. At ∼10 AU, the Alfvén speed $V_A=B_0/\sqrt{\mu _0 n_0 m_p}=38.56$ km s⁻¹. We further assume that the polytropic indices γ_s = γ_p = γ_e = 5/3 and that the electron temperature T_e = T_s. The proton and electron sound speeds can be estimated from $C_s^2=\gamma _s k_B T_s /m_p$ and $C_e^2=\gamma _e k_B T_e /m_p$, which yields C_s = C_e = 16.59 km s⁻¹. We assume a solar wind speed U₀ = 400 km s⁻¹, from which we estimate the pickup ion sound speed $C_p=U_0 \sqrt{\gamma _p /7}=195.18$ km s⁻¹ under the assumption that the PUI distribution function may be approximated as a filled shell (Vasyliunas & Siscoe 1976; Zank 1999; Zank et al. 2010). The diffusion coefficient is estimated using 〈c²〉^1/2 = 400 km s⁻¹.

Plotted in Figure 7, using the same format, as before are the dispersion curves for the two-fluid and three-fluid models for the supersonic solar wind. Consider the top panels first (parallel propagation). The dispersion relation for the three-fluid model is of the 11th-order in ω and generally (at large scales) admits five forward propagating modes, five backward propagating modes and one mode which is purely imaginary. Until now, we have plotted only the five forward modes. Here we plot two additional modes. Starting at large scales (small k), the cyan crosses identify the backward propagating PUI slow sound mode. The real frequency of this mode is negative (and therefore not plotted in the logarithmic scale), but its absolute value is equal to the real frequency of the forward propagating mode (cyan line). The damping rates for the forward (cyan) and backward (cyan crosses) modes are naturally the same, as is illustrated in Figure 7. At kV_A/Ω_p ∼ 0.7, both the forward and backward PUI slow sound modes become so damped that they lose their identity and become purely imaginary. After this point, the damping rate of these modes is not equal. (Since the modes lost their identity, the decision about which is cyan and which has the cyan crosses is technically ad hoc.) The magenta curve represents the mode that is purely imaginary at large scales. At approximately kV_A/Ω_p ∼ 18 the damping rates for the magenta mode and the cyan crosses mode become equal. At this point, two new modes with ω_r ≠ 0 are generated, one forward (cyan crosses) and one backward (magenta).

Zoom In Zoom Out Reset image size

As shown in the second set of panels (θ = 30°), at approximately kV_A/Ω_p ∼ 0.9, the fast PUI mode (black) becomes so damped that it loses its identity and becomes purely imaginary. For this reason, we now plot the backward propagating PUI fast mode, represented by black crosses at large scales. After both modes lose their identity at kV_A/Ω_p ∼ 0.9, their damping rates become unequal and identification of these modes is an ad hoc choice. At approximately kV_A/Ω_p ∼ 20, the backward propagating PUI (black-crosses) mode matches the frequency of the purely imaginary (magenta) mode. At this point, a new forward propagating PUI fast mode (black crosses) and a new backward propagating mode (magenta line) is generated. The PUI fast (black line) mode now becomes the only purely imaginary mode existing in the system.

Like the θ = 30° case, at 70° the forward (black line) and backward (black crosses) propagating PUI fast modes lose their identities and become purely imaginary at about kV_A/Ω_p ∼ 6. However, the backward propagating PUI (black crosses) mode never reaches the purely imaginary magenta mode for the k values considered here and no new modes with nonzero real frequency are generated.

For completeness, we discuss briefly the nature of wave modes admitted by the reduced one-fluid model (61)–(68). The correspondence between these modes and those obtained from the multi-fluid description is noted. The linearization proceeds in the usual way, where we use

$$\begin{eqnarray*} \rho &=& \rho _0 + \delta \rho; \quad {\bf U} = \delta {\bf u}; \quad P = P_0 + \delta P; \quad\nonumber\\ P_p &=& {P_p}_0 + \delta P_p; \quad {\bf B} = B_0 \hat{\bf z} + \delta {\bf B}, \end{eqnarray*}$$

and suppose that the fluctuations are ∝exp i[ωt − k · x] and k = (k_x, 0, k_z). We neglect the higher order PUI viscosity terms and retain only the effects of PUI heat conduction. We obtain the dispersion relation in the form

$$\begin{eqnarray} && \left(\omega ^2 - V_{A}^2 k^2 \cos ^2 \theta \right) \left[ \omega ^4 - \left(V_S^2 + \frac{\omega }{\omega - i {\bf k} \cdot {\bf K} \cdot {\bf k} /3} C_{p}^2 + V_{A}^2 \right)\right. \nonumber \\ &&\left. \quad \times k^2 \omega ^2 {} + \left(V_S^2 + \frac{\omega }{\omega - i {\bf k} \cdot {\bf K} \cdot {\bf k} /3} C_{p}^2 \right) V_{A}^2 k^4 \cos ^2 \theta \right] = 0, \nonumber\\ \end{eqnarray} \tag{ 127 }$$

using the definitions (75)–(76) and $V_S^2 \equiv \gamma P_0 /\rho _0$, and θ is the wave propagation angle. The form of (127) is used since it expresses the correspondence with the standard dispersion relation of MHD. By setting k · K · k = 0 (the strong coupling limit) and identifying $C_{m}^2 \equiv V_S^2 + C_{p}^2$ as the sound speed, (127) is identical to the MHD dispersion relation. Alfvén wave propagation is unaffected by PUIs. The expression inside the square bracket of (127) is the dispersion relation for PUI mediated fast and slow magnetosonic modes, by analogy with the multi-fluid analysis.

For wave mode solutions, Equation (127) is usefully expressed as

$$\begin{eqnarray} && \left(\omega - \frac{i}{3} {\bf k} \cdot {\bf K} \cdot {\bf k} \right) \left(\omega ^4 - \left(V_S^2 + V_A^2 \right) k^2 \omega ^2 + V_S^2 V_{A}^2 k^4 \cos ^2 \theta \right) \nonumber\\ && \quad = \omega \left(\omega ^2 - V_{A}^2 k^2 \cos ^2 \theta \right) C_{p}^2 k^2. \end{eqnarray} \tag{ 128 }$$

Before solving (127) numerically, it is useful to introduce κ ≡ κ_∥cos ²θ + κ_⊥sin ²θ and derive long and short wavelength approximate solutions. We provide the parallel and perpendicular cases first. In the long wavelength limit, we obtain for θ = 0 (parallel propagation)

$$\begin{eqnarray} \omega &=& \pm \sqrt{V_S^2 + C_{p}^2 } k+ \frac{i \kappa k^2}{6} \frac{C_{p}^2 }{V_S^2 + C_{p}^2 }; \end{eqnarray} \tag{ 129 }$$

$$\begin{eqnarray} \omega &=& \pm V_{A} k; \end{eqnarray} \tag{ 130 }$$

$$\begin{eqnarray} \omega &=& \frac{i}{3}\kappa _\parallel k^2 \frac{V_S^2}{V_S^2+C_p^2}. \end{eqnarray} \tag{ 131 }$$

Equation (129) corresponds to the parallel sound mode propagating at the combined sound speed (i.e., the sum of the square of the thermal plasma sound speed and the PUI sound speed), which is damped by PUIs. Note that in the one-fluid model, the thermal sound and PUI waves are no longer distinct as they were in the multi-fluid model (Equations (94) and (95)). The other parallel mode admitted by (130) is unaffected by PUIs, and corresponds to the Alfvén wave solution obtained from the first bracketed term of (127). Finally, Equation (131) describes a purely imaginary non-propagating mode. For perpendicular propagation, θ = π/2, the long wavelength limits become

$$\begin{eqnarray} \omega &=& \pm \sqrt{V_S^2 + C_{p}^2 + V_{A}^2 } k + \frac{i \kappa k^2}{6} \frac{C_{p}^2 }{V_S^2 + C_{p}^2 + V_{A}^2 }; \end{eqnarray} \tag{ 132 }$$

$$\begin{eqnarray} \omega &=& \frac{i \kappa k^2}{3} \frac{V_S^2 + V_{A}^2 }{V_S^2 + C_{p}^2 + V_{A}^2 }. \end{eqnarray} \tag{ 133 }$$

The fast mode speed now includes the Alfvén speed, in accordance with MHD, as well as the PUI sound speed, and is damped by PUIs. The slow mode is non-propagating and purely damped. The comments regarding the comparison to the multi-fluid modes above hold for perpendicular propagation as well.

The short wavelength limits of (128) for parallel propagation are

$$\begin{eqnarray} \omega &=& \pm V_S k + \frac{3i }{2 \kappa } C_{p}^2; \end{eqnarray} \tag{ 134 }$$

$$\begin{eqnarray} \omega &=& \pm V_{A} k, \end{eqnarray} \tag{ 135 }$$

and for θ = π/2,

$$\begin{eqnarray} \omega &=& \pm \sqrt{V_S^2 + V_{A}^2 } k + \frac{3i}{2 \kappa } C_{p}^2; \end{eqnarray} \tag{ 136 }$$

$$\begin{eqnarray} \omega &=& 0. \end{eqnarray} \tag{ 137 }$$

The PUIs no longer contribute to the fast mode sound speed, as they did in the long wavelength regime, but the PUIs do damp the thermal plasma fast mode. The slow mode in the short wavelength limit is completely unaffected by the presence of PUIs.

Consider now the general case of oblique wave propagation. For all angles of propagation, we always have the classical Alfvén wave solution

$$\begin{equation} \omega = \pm V_A k \cos \theta, \end{equation} \tag{ 138 }$$

which is completely unaffected by the presence of a suprathermal PUI population. The obliquely propagating long wavelength mode is somewhat complicated, with

$$\begin{eqnarray} \omega &\simeq & k \sqrt{\frac{1}{2} \left(V_A^2+V_S^2+C_p^2 \right)(1\pm {\mathcal {A}})} {} + \frac{i}{6}\kappa k^2 C_p^2 \nonumber \\ && \times \frac{ \left(V_A^2+V_S^2+C_p^2 \right)(1\pm {\mathcal {A}}) -2V_A^2\cos ^2\theta }{\left(V_A^2+V_S^2+C_p^2 \right)^2(1\pm {\mathcal {A}}) -4V_A^2 \left(V_S^2+C_p^2 \right)\cos ^2\theta }, \nonumber\\ \end{eqnarray} \tag{ 139 }$$

where

$$\begin{eqnarray*} &&\mathcal {A} = \sqrt{1-\frac{4V_A^2 \left(V_S^2+C_p^2 \right)\cos ^2\theta }{ \left(V_A^2+V_S^2+C_p^2 \right)^2}}. \end{eqnarray*}$$

The purely imaginary mode in the long wavelength limit is

$$\begin{equation} \omega \simeq \frac{i}{12}\kappa k^2 \frac{C_p^2}{V_S^2+C_p^2}. \end{equation} \tag{ 140 }$$

For short wavelength modes propagating obliquely, we have the general result

$$\begin{eqnarray} \omega \simeq && k \sqrt{\frac{1}{2} \left(V_A^2+V_S^2 \right) \pm \frac{1}{2} \sqrt{\left(V_A^2+V_S^2 \right)^2 - 4 V_A^2 V_S^2 \cos ^2 \theta } } \nonumber \\ && + \frac{3i}{4\kappa } C_p^2 \left(1 \pm \frac{V_S^2+(1-2\cos ^2\theta)V_A^2 }{\sqrt{\left(V_A^2+V_S^2 \right)^2 - 4V_A^2 V_S^2 \cos ^2\theta } } \right). \quad\quad \end{eqnarray} \tag{ 141 }$$

The real frequency is unaffected by PUIs and is identical to the familiar MHD fast and slow magnetosonic speeds. Damping of the short wavelength fast and slow modes is directly due to PUIs, however. The parallel and perpendicular limits are consistent with our previous results.

The purely imaginary non-propagating short wavelength mode is given by

$$\begin{equation} \omega \simeq \frac{3i}{2\kappa }C_p^2 \frac{V_A^2\cos ^2\theta }{V_A^2+V_S^2}. \end{equation} \tag{ 142 }$$

For parallel and perpendicular propagation, we have

$$\begin{eqnarray} &&\omega \simeq \frac{3i}{2\kappa _\parallel }\left(V_S^2+C_p^2 \right); \end{eqnarray} \tag{ 143 }$$

$$\begin{eqnarray} &&\omega \simeq \frac{3i}{2\kappa _\perp }\left(V_A^2+V_S^2+C_p^2 \right). \end{eqnarray} \tag{ 144 }$$

Numerical solutions of the dispersion relation (127) are plotted in Figure 8 for VLISM parameters. The heat conduction diffusion coefficient was estimated assuming 〈c²〉^1/2 = 200 km s⁻¹. The left column shows the real frequency plotted as a function of k and the right column shows the corresponding imaginary frequency. We do not plot the purely imaginary mode. Only three real modes exist—the fast (blue) and slow (red) magnetosonic modes modified by PUIs, and the Alfvén wave (green). For parallel propagation, the slow mode coincides with the Alfvén mode, and there is a region around k ≃ Ω_p/V_A for which the fast mode is evanescent. For quasi-parallel propagation (e.g., θ ⩽ 30°), the region of evanescence for the fast mode is always present, but for quasi-perpendicular propagation (e.g., third panel of Figure 8), the fast mode exists for all wave numbers k. For the VLISM parameters used here, the slow mode and Alfvén wave are virtually indistinguishable, although the former is damped by PUIs.

Zoom In Zoom Out Reset image size

Pickup ions typically dominate the internal energy and do not equilibrate collisionally in plasma systems such as the VLISM, the IHS, and the supersonic solar wind beyond some 10 AU. This requires that we treat PUIs as a distinct plasma component. PUIs experience pitch-angle scattering from pre-existing and self-excited magnetic field turbulence. We have developed a collisionless form of Chapman–Enskog expansion that is correct to second order to effect the closure of the moments of the PUI transport equation. The first-order expansion of the PUI distribution function introduces the collisionless PUI heat flux through a heat conduction diffusion tensor. We show that the collisionless PUI heat conduction term is proportional to the gradient in the PUI pressure. The second-order expansion of the PUI distribution function yields the collisionless stress tensor for PUIs, whose structure generalizes the familiar form from hydrodynamics (e.g., Zank 2013). By deriving expressions for the PUI heat flux and stress tensor in terms of lower order moments (PUI pressure and bulk velocity) we can derive a closed form of multi-fluid equations for a plasma system comprised of thermal protons and electrons and suprathermal PUIs. The system of equations differs from the typical multi-fluid system in that we have made no ad hoc assumptions about the PUI distribution being isotropic but instead derived the heat flux and stress tensor corrections subject to weak departures from isotropy (i.e., assuming that the scattering of PUIs by turbulent magnetic fluctuations proceeds rapidly).

Since many problems are too complicated to solve using an elaborate multi-fluid model, we develop an analogue to the MHD equations that resembles a single-fluid description but possesses a separate equation of state for the PUIs, together with a separate equation of state for the thermal background plasma. In so doing, the single-fluid description now includes both collisionless heat conduction and stress tensor (collisionless viscosity) tensor terms associated with the PUIs. The derived single-fluid model possesses an interesting structural similarity with the so-called cosmic ray MHD model (Webb 1983). The single-fluid PUI-mediated MHD equations presented here are almost identical to the cosmic ray (CR) MHD equations but the derivation is entirely different. The derivation of the CR MHD equations explicitly neglects all the flow acceleration and shear terms (which, as we have shown above, is essential in the derivation of the closed fluid system describing the pitch-angle scattered charged particles), does not derive the underlying multi-fluid description from which an MHD description is developed, explicitly assumes that the cosmic ray proton number density is zero, and assumes an ad hoc decomposition of the total isotropic scalar pressure into thermal and cosmic ray components. Of particular interest is that the PUI number density is included in the PUI-mediated single-fluid equations, suggesting that a more careful corresponding derivation of the CR MHD equations could also lead to the explicit inclusion of the CR number density. The importance of this is that it allows the "injection" problem for the acceleration of cosmic rays at collisionless shocks to be addressed within the context of the single-fluid CR MHD equations (e.g., Zank et al. 1993).

We have used our new PUI-mediated plasma model to investigate the nature of linear waves in the context of the VLISM, the IHS or the subsonic solar wind, and the supersonic solar wind beyond ∼10 AU using the full multi-fluid description. The general dispersion relation in the presence of PUIs and a background thermal plasma comprised of protons and electrons was derived—this includes the effects of the collisionless heat conduction (diffusive) term, but we neglect the effects of the (viscous) stress tensor terms. The retention of the electron inertia terms yields a formidably complicated dispersion relation that does not lend itself to analytic representation. Consequently, the fully general dispersion relation is solved numerically only. Prior to solving the full dispersion relation (i.e., with electrons, thermal protons, and PUIs retained), we derive an analytically "tractable" dispersion relation under the assumption of massless electrons. An 11th-order polynomial dispersion relation is obtained in the low-frequency (ω  ≪  Ω_e) limit, which we investigate analytically for a variety of cases. The class of admissible wave solutions is of course richer than in the standard two-fluid case, with as many as five distinct propagating (real frequency) wave modes are possible: forward and backward propagating Alfvén waves, thermal fast and slow magnetoacoustic waves, and PUI fast and slow magnetoacoustic waves, as well as an entropy mode. At least one purely damped, non-propagating mode also exists. The thermal wave modes are characterized by their real frequency being proportional to the thermal sound speed $\sqrt{C_s^2 + C_e^2}$ to leading order, whereas the PUI modes are proportional to the PUI sound speed C_p to leading order. However, all the waves are damped by PUIs, whether thermal or PUI modes. These results were obtained from an analytic expansion of the 11th-order dispersion relation in the short and long wavelength limits for parallel, perpendicular, and oblique wave propagation. The analysis shows that the damping of both the long and short wavelength waves is a direct consequence of the collisionless PUI heat conduction, i.e., diffusive damping. In the short wavelength limit, we find that the PUI modes are absent, with only the thermal modes and a mixed mode existing.

The analytic limits are interesting and provide useful insight and guidance in understanding and classifying the numerical solutions that are necessary to understand the full complexity of the dispersion relations when electron inertial terms are retained. In solving the linearized full multi-fluid system numerically, we also examine the reduced two-fluid solutions that emerge from the abundance limit α = 1. In the α = 1 limit, the full dispersion relation separates into separate two-fluid and PUI fluid dispersion relations. The two-fluid dispersion relation is the familiar sixth-order polynomial relation first investigated by Stringer (1963; we have corrected a number of simplifications and typos that were present in earlier derivations of the two-fluid dispersion relation). Solutions of the two-fluid dispersion relation are useful in clarifying the numerical solutions when α < 1, and in showing explicitly the role that PUIs play in mediating wave properties. We considered parameter regimes for three distinct plasma environments that are mediated by PUIs created in part by interstellar neutral H.

A major application of the linear wave theory was to the plasma of the VLISM. As discussed, the IBEX observations appear to confirm the notion that the VLISM is comprised of a non-equilibrated admixture of thermal plasma and energetic PUIs created by fast neutrals originating in the supersonic and hot subsonic (IHS) solar wind (Zirnstein et al. 2014; Desai et al. 2014). The VLISM parameters can be reasonably estimated, including the abundance parameter α (Zirnstein et al. 2014), but the choice of value for the collisionless heat conduction term is less well constrained in the context of the multi-fluid model. Two estimates of the diffusion coefficient K were considered, one corresponding to the weak coupling of PUIs to the thermal plasma (i.e., a large value of the diffusion coefficient for the collisionless PUI heat flux) and the other to a slightly stronger coupling. Both values of the diffusion coefficient appear to be reasonable but the results for the two cases differ quite significantly. Both cases can yield as many as five distinct wave modes (depending on the wave number k). We find the following.

• 1.  
  For κ_∥ large, the thermal and PUI modes are essentially decoupled from one another in that the dispersion curves (the real frequency as a function of k) of the thermal modes (fast, slow, and Alfvén) very closely resemble those obtained from the two-fluid model (i.e., in the absence of PUIs). The fast and slow PUI modes are present and the thermal and PUI dispersion curves can intersect but there is no corresponding modification to either mode in the regions of intersection in (ω, k) space.
• 2.  
  Even when the diffusion coefficient is large and the thermal and PUI fluid essentially decouple, the PUIs nevertheless act to damp both the thermal and PUI modes, and the Alfvén mode (illustrated in the plots of the imaginary frequency as a function of k).
• 3.  
  The interplay between the fast and slow PUI modes is interesting for the larger choice of diffusion coefficient. For parallel wave propagation, the slow PUI wave is evanescent above a certain wave number and the fast PUI mode is undamped and satisfies ω = Ω_p. However, for oblique wave propagation, the fast mode is evanescent at short wavelength scales and the slow PUI wave is asymptotic to Ω_pcos θ at short wavelengths.
• 4.  
  For κ_∥ smaller, the PUI and the thermal modes are fully coupled and the two-fluid results are no longer a useful guide to the oblique wave properties in either the long or short wavelength regimes, particularly for wave numbers approximately >Ω_p/V_A. Only for parallel propagation do the thermal wave mode dispersion curves resemble the two-fluid dispersion relations.
• 5.  
  For smaller κ_∥, the various dispersion curves no longer intersect as they do for a larger diffusion coefficient. Instead, the dispersion curves in the vicinity of an apparent intersection point are modified in that they approach one another in a "hyperbolic" sense but do not cross. Furthermore, the properties of the converging dispersion curves "exchange properties" in these regions. For example, the dispersion curve for the slow PUI mode resembles what was the parallel thermal fast mode curve at short wavelengths.
• 6.  
  Of possible particular importance is that the dispersion curve for Alfvén waves completely changes character from the large diffusion coefficient case to the smaller diffusion case. In the former example, the Alfvén wave dispersion curve is asymptotic to Ω_ecos θ and in the latter to Ω_pcos θ.

We do not summarize the detailed properties of the linear wave analysis applied to the IHS or subsonic solar wind. The wave properties resemble the smaller diffusion coefficient VLISM case despite the nominal diffusion coefficient being larger in the IHS than in the VLISM. This is because the abundance parameter α is smaller than it is in the VLISM, so ensuring a stronger coupling of thermal and PUI plasma. Similarly, our analysis of the linear wave modes in the supersonic solar wind beyond ∼10 AU revealed that the wave properties were similar to those of the small diffusion coefficient VLISM example. As part of the supersonic solar wind analysis, we considered the backward propagating PUI modes as well, showing how fast forward and backward propagating PUI modes could lose their identity and become purely evanescent or imaginary. Correspondingly, when a backward propagating PUI mode matches the frequency of an imaginary PUI mode, then new forward and backward propagating modes are generated.

For completeness, we considered linear wave solutions of the reduced one-fluid PUI-mediated plasma model. Because the reduction to a single plasma with equations of state for both constituents (thermal plasma and PUIs) eliminates the distinction between the plasma components, the distinction between PUI and thermal wave modes is lost.

• 1.  
  The single-fluid or one-fluid model only admits three propagating waves—the fast mode, slow mode, Alfvén mode, and a purely imaginary mode.
• 2.  
  Long wavelength modes propagate at a mixed fast or slow magnetoacoustic speed where the sound speed is now a combination of the thermal sound speed and the PUI sound speed, i.e., $\sqrt{\vphantom{A^A}\smash{{V_S^2 + C_p^2}}}$.
• 3.  
  Short wavelength modes propagate at the fast or slow magnetosonic speed where the sound speed is given by the thermal sound speed only, i.e., PUIs do not contribute to the phase velocity of short wavelength waves.
• 4.  
  Both fast and slow mode, obliquely propagating modes are damped by the heat flux term associated with PUIs. For parallel propagation in both the long and short wavelength limits, the slow mode wave is identical to the Alfvén wave, and is undamped. This does not hold for oblique wave propagation.
• 5.  
  Unlike the multi-fluid model, Alfvén waves in the one-fluid reduced model are unaffected by PUIs and the dispersion curve continues to increase with increasing k, as in the case of ideal MHD.

The compressibility and polarization of the various waves was not discussed here. Furthermore, we note that the inclusion of Landau damping can modify the multi-fluid description in important ways, as described in Hunana et al. (2013).

One important implication to emerge from our three-fluid analysis of linear waves is related to the contrasting behavior of the two-fluid waves, particularly the Alfvén mode in the case of strong PUI and thermal plasma coupling. Our results suggest that turbulent fluctuations in the inner heliosphere may be fundamentally different from those in regions such as the supersonic solar wind beyond ∼10 AU, the IHS, or the VLISM. Precisely how this translates into e.g., observed magnetic field fluctuation spectra is an open question, but we note that Burlaga et al. (2006); Burlaga & Ness (2009) have observed "turbulence" in the IHS that appears to possess characteristics (distribution, compressibility) somewhat different from those observed in the inner heliosphere (e.g., Bruno & Carbone 2005). Initial observations of the VLISM by the Voyager 1 magnetometer (Burlaga & Ness 2014) suggest a surprisingly low level of magnetic field fluctuations. In the absence of significant external driving, it is possible, based on the results presented here, that PUIs act to damp waves in the VLISM significantly, leaving the medium relatively quiescent.

A second important implication is related to heavy ions that are minor components in the solar wind of the inner heliosphere and in cometary environments. Mann et al. (1997) used a three-fluid system of fluid equations to derive the dispersion relation for a proton–electron-heavy ion plasma and investigated the wave properties based on this model. However, the model is closed by assuming Maxwellian distributions for all the species. In both a cometary environment and in the inner heliosphere, heavy ions experience pitch-angle scattering off ambient and excited turbulent magnetic fluctuations and are not equilibrated with the background thermal plasma. In view of the important differences that we have identified in the dispersion curves as a consequence of deriving and including the collisionless heat flux associated with PUIs (and in principle the collisionless stress tensor) into a multi-fluid model, we suspect that that similar differences might well arise from extending our approach to heavy ions.

In subsequent papers, we plan to examine further the damping of waves in PUI mediated plasmas, including heavy ions, explore nonlinear waves and structures admitted by this system, and investigate the structure of shock waves in the presence of PUIs.

Finally, we note that the PUI mediated plasma model presented here can find application to a wide variety of astrophysical situations, ranging from the interaction of solar-like stars with different interstellar conditions (e.g., Müller et al. 2006; Florinski et al. 2004) to the propagation of shock waves in a partially ionized ISM (e.g., Raymond et al. 2008).

This work was supported by IBEX NASA/SwRI award NNG05EC85C, subcontract A99132BT and NASA grant NNX10AC17G. G.P.Z. appreciates discussions at the team meeting "Heliosheath Processes and Structure of the Heliopause: Modeling Energetic Particles, Cosmic Rays, and Magnetic Fields" supported by the International Space Science Institute in Bern, Switzerland.

The gyrophase-averaged equation for non-relativistic particles is given by

$$\begin{eqnarray} &&\frac{\partial f}{\partial t} + \left(U_i + c \mu b_i \right) \frac{\partial f}{\partial x_i} + \frac{1 - \mu ^2}{2} \left[ c \nabla \cdot {\bf b} + \mu \nabla \cdot {\bf U} - 3 b_i b_j \frac{\partial U_j}{\partial x_i} - \frac{2 b_i}{c} \frac{D U_i}{Dt} \right] \frac{\partial f}{\partial \mu } \nonumber \\ && \quad\quad {} + \left[ \frac{1 - 3 \mu ^2 }{2} b_i b_j \frac{\partial U_j}{\partial x_i} - \frac{1 - \mu ^2}{2} \nabla \cdot {\bf U} - \frac{\mu b_i}{c} \frac{D U_i}{Dt} \right] c \frac{\partial f}{\partial c} = \frac{\partial }{\partial \mu } \left(\nu _s (1 - \mu ^2) \frac{\partial f}{\partial \mu } \right), \end{eqnarray} \tag{ A1 }$$

where the direction vector b ≡ B/B. We use a Legendre polynomial expansion to solve the gyrophase-averaged equation (A1), expanding the gyrophase-averaged particle distribution function f as

$$\begin{eqnarray*} &&f({\bf x}, t, c, \mu) = \sum _{n=0}^{\infty } \frac{1}{2} (2n + 1) P_n (\mu) f_n ({\bf x}, t, c), \: {\rm where} \: f_n ({\bf x}, t, c) = \int _{-1}^1 f P_n (\mu) d\mu. \end{eqnarray*}$$

Following the derivation in Zank (2013), we obtain an infinite set of partial differential equations in the coefficients f_n of the Legendre polynomials,

$$\begin{eqnarray} && \frac{\partial f_m}{\partial t} + U_i \frac{\partial f_m}{\partial x_i} + \frac{c b_i}{2m+1} \left[ (m+1) \frac{\partial f_{m+1}}{\partial x_i} + m \frac{\partial f_{m-1}}{\partial x_i} \right] + c b_i b_j \frac{\partial U_j}{\partial x_i} \frac{1}{2} \frac{\partial f_m}{\partial c} {} - c b_i b_j \frac{\partial U_j}{\partial x_i} \frac{3}{2} \frac{(m-1)m}{(2m+1)(2m-1)} \frac{\partial f_{m-2}}{\partial c} \nonumber \\ && \quad - c b_i b_j \frac{\partial U_j}{\partial x_i} \frac{3}{2} \left(\frac{(m+1)^2}{2m+3} + \frac{m^2}{2m-1} \right) {} \frac{1}{2m+1} \frac{\partial f_{m}}{\partial c} - c b_i b_j \frac{\partial U_j}{\partial x_i} \frac{3}{2} \frac{ (m+2)(m+1) }{(2m+3)(2m+1)} \frac{\partial f_{m+2}}{\partial c} - \frac{c}{2} \frac{\partial U_i}{\partial x_i} \frac{\partial f_m}{\partial c} \nonumber \\ && \quad {} + \frac{c}{2} \frac{\partial U_i}{\partial x_i} \frac{(m-1)m}{(2m+1)(2m-1)} \frac{\partial f_{m-2}}{\partial c} + \frac{c}{2} \frac{\partial U_i}{\partial x_i} \left(\frac{(m+1)^2}{2m+3} + \frac{m^2}{2m-1} \right) \frac{1}{2m+1} \frac{\partial f_{m}}{\partial c} {} + \frac{c}{2} \frac{\partial U_i}{\partial x_i} \frac{ (m+2)(m+1) }{(2m+3)(2m+1)} \frac{\partial f_{m+2}}{\partial c} \nonumber \\ && \quad - \frac{DU_i}{Dt} \frac{b_i}{2m+1} \left[ (m+1) \frac{\partial f_{m+1}}{\partial c} + m \frac{\partial f_{m-1}}{\partial c} \right] {} + \frac{c}{2} \frac{\partial b_i}{\partial x_i} \left[ \frac{(m+1)(m+2)}{2m+1} f_{m+1} - \frac{m(m-1)}{2m+1} f_{m-1} \right] \nonumber \\ && \quad {} + \frac{1}{2} \frac{\partial U_i}{\partial x_i} \frac{(m+1)(m+2)(m+3)}{(2m+1)(2m+3)} f_{m+2} + \frac{1}{2} \frac{\partial U_i}{\partial x_i} \frac{m(m+1)}{(2m-1)(2m+3)} f_m {} - \frac{1}{2} \frac{\partial U_i}{\partial x_i} \frac{m(m-1)(m-2)}{(2m+1)(2m-1)} f_{m-2} \nonumber \\ && \quad - \frac{3}{2} b_i b_j \frac{\partial U_j}{\partial x_i} \frac{(m+1)(m+2)(m+3)}{(2m+1)(2m+3)} f_{m+2} {} - \frac{3}{2} b_i b_j \frac{\partial U_j}{\partial x_i} \frac{m(m+1)}{(2m-1)(2m+3)} f_m + \frac{3}{2} b_i b_j \frac{\partial U_j}{\partial x_i} \frac{m(m-1)(m-2)}{(2m+1)(2m-1)} f_{m-2} \nonumber \\ && \quad {} - \frac{b_i}{c} \frac{DU_i}{Dt} \left[ \frac{(m+1)(m+2)}{2m+1} f_{m+1} - \frac{m(m-1)}{2m+1} f_{m-1} \right] = - \nu _s \frac{ m (m+1) }{2} f_m. \end{eqnarray} \tag{ A2 }$$

For the f₂ approximation (i.e., assume f_n = 0 ∀ n ⩾ 3), we have, on setting m = 0,

$$\begin{eqnarray} &&\frac{\partial f_0}{\partial t} + U_i \frac{\partial f_0}{\partial x_i} - \frac{c}{3} \frac{\partial U_i}{\partial x_i} \frac{\partial f_0}{\partial c} = - \frac{\partial c b_i f_{1}}{\partial x_i} + \frac{DU_i}{Dt} b_i \frac{1}{c^2} \frac{\partial c^2 f_{1}}{\partial c} - \frac{1}{3} \frac{\partial U_i}{\partial x_i} \frac{1}{c^2} \frac{\partial }{\partial c} (c^3 f_2) {} + b_i b_j \frac{\partial U_j}{\partial x_i} \frac{1}{c^2} \frac{\partial }{\partial c} (c^3 f_2), \end{eqnarray} \tag{ A3 }$$

and on setting m = 1 and neglecting all terms with indices having i ⩾ 3, we obtain

$$\begin{eqnarray} &&\frac{\partial f_1}{\partial t} + U_i \frac{\partial f_1}{\partial x_i} +\frac{2}{3} c b_i \frac{\partial f_2}{\partial x_i} - \frac{2}{5} c b_i b_j \frac{\partial U_j}{\partial x_i} \frac{\partial f_1}{\partial c} + \frac{2}{5} c \frac{\partial U_i}{\partial x_i} \frac{\partial f_1}{\partial c} -\frac{2}{3} b_i \frac{D U_i}{Dt} \frac{\partial f_2}{\partial c} + c \frac{\partial b_i}{\partial x_i} f_2 \nonumber \\ &&{} + \frac{1}{5} \frac{\partial U_i}{\partial x_i} f_1 - \frac{3}{5} b_i b_j \frac{\partial U_j}{\partial x_i} f_1 -\frac{2 b_i}{c} \frac{D U_i}{Dt} f_2 = - \nu _s f_1 - \frac{c b_i}{3} \frac{\partial f_0}{\partial x_i} + \frac{DU_i}{Dt} \frac{b_i}{3} \frac{\partial f_0}{\partial c}, \end{eqnarray} \tag{ A4 }$$

where the f₀ Legendre coefficients are expressed as source terms in the evaluation of the next order Legendre coefficients f₁. To solve Equation (A4) for f₁ in terms of the lower order Legendre coefficient f₀, we assume that the zeroth-order coefficient f₀ is almost isotropic, implying that f₁ ≪ f₀. We also assume that $\nu _s = \tau _s^{-1}$ is large, i.e., rapid scattering of the charged particles (which is consistent with the assumption that the particle distribution is nearly isotropic), so that the term ν_sf₁ ∼ O(f₀). Subject to these assumptions, Equation (A4) can then be solved, yielding

$$\begin{equation} \nu _s f_1 \simeq - \frac{c b_i}{3} \frac{\partial f_0}{\partial x_i} + \frac{DU_i}{Dt} \frac{b_i}{3} \frac{\partial f_0}{\partial c}. \end{equation} \tag{ A5 }$$

Using (A5) in (A3) to rewrite f₁ in terms of f₀ yields the familiar transport equation in the presence of large-scale gradients and accelerations (see Zank 2013) if the f₂ terms are neglected,

$$\begin{equation} \frac{\partial f_0}{\partial t} + {\bf U} \cdot f_0 - \frac{c}{3} \nabla \cdot {\bf U} \frac{\partial f_0}{\partial c} + \nabla \cdot \left({\bf K} \frac{D {\bf U} }{Dt} \frac{1}{c} \frac{\partial f_0}{\partial c} \right) = \nabla \cdot \left({\bf K} \cdot \nabla f_0 \right), \end{equation} \tag{ A6 }$$

where K is the diffusion tensor.

The m = 2 expression can be written as

$$\begin{eqnarray*} &&\frac{\partial f_2}{\partial } + U_i \frac{\partial f_2}{\partial x_i} + \frac{2c}{5} b_i \frac{\partial f_1}{\partial x_i} + \frac{c}{2} b_i b_j \frac{\partial U_j}{\partial x_i} \frac{\partial f_2}{\partial c} - \frac{c}{5} b_i b_j \frac{\partial U_j}{\partial x_i} \frac{\partial f_0}{\partial c} - \frac{3c}{10} b_i b_j \left(\frac{9}{7} + \frac{4}{3} \right) \frac{\partial f_2}{\partial c} \\ && \quad\quad{} - \frac{c}{2} \frac{\partial U_i}{\partial x_i} \frac{\partial f_2}{\partial c} + \frac{c}{15} \frac{\partial U_i}{\partial x_i} \frac{\partial f_2}{\partial c} + \frac{c}{10} \frac{\partial U_i}{\partial x_i} \left(\frac{9}{7} + \frac{4}{3} \right) \frac{\partial f_2}{\partial c} - \frac{2}{5} b_i \frac{D U_i}{Dt} \frac{\partial f_1}{\partial c} - \frac{c}{5} \frac{\partial b_i}{\partial x_i} f_1 \\ && \quad\quad{} + \frac{1}{7} \frac{\partial U_i}{\partial x_i} f_2 - \frac{3}{7} b_i b_j \frac{\partial U_j}{\partial x_i} f_2 + \frac{2}{5c} b_i \frac{D U_i}{Dt} f_1 = -3 \nu _s f_2, \end{eqnarray*}$$

from which we obtain the approximate solution for f₂,

$$\begin{equation} f_2 \simeq \frac{c \tau _s}{15} \left(b_i b_j \frac{\partial U_j}{\partial x_i} - \frac{1}{3} \frac{\partial U_i}{\partial x_i} \right) \frac{\partial f_0}{\partial c}. \end{equation} \tag{ A7 }$$

The coefficients of the dispersion relation (81) are given by

$$\begin{eqnarray} && A_{11} = 1, \end{eqnarray} \tag{ B1 }$$

$$\begin{eqnarray} && A_{10} = -\frac{i}{3}\frac{{\bf k} \cdot {\bf K} \cdot {\bf k} }{\Omega _p}, \end{eqnarray} \tag{ B2 }$$

$$\begin{eqnarray} && A_9 = -k^4\frac{V_A^4}{\Omega _p^4}\cos ^2\theta -k^2\frac{V_A^2}{\Omega _p^2}\left [ 1+\cos ^2\theta +\frac{C_e^2}{V_A^2}+\frac{C_s^2}{V_A^2}+\frac{C_p^2}{V_A^2}\right ] - 1, \end{eqnarray} \tag{ B3 }$$

$$\begin{eqnarray} && A_8 = +\frac{i}{3}\frac{{\bf k} \cdot {\bf K} \cdot {\bf k} }{\Omega _p}\left [ k^4\frac{V_A^4}{\Omega _p^4}\cos ^2\theta +k^2\frac{V_A^2}{\Omega _p^2}\left (1+\cos ^2\theta +\frac{C_e^2}{V_A^2}+\frac{C_s^2}{V_A^2}\right)+1 \right ], \end{eqnarray} \tag{ B4 }$$

$$\begin{eqnarray} && A_7 = +k^6\frac{V_A^6}{\Omega _p^6}\cos ^2\theta \left [ \frac{C_e^2}{V_A^2}+\frac{C_p^2}{V_A^2}+\frac{C_s^2}{V_A^2}\right ] +k^4\frac{V_A^4}{\Omega _p^4}\left [ (1-\alpha)\frac{C_s^2}{V_A^2}\left (\frac{C_e^2}{V_A^2}+1\right) +\frac{C_p^2}{V_A^2}\left (\alpha +\alpha \frac{C_e^2}{V_A^2}+\frac{C_s^2}{V_A^2}\right)\right. \nonumber \\ &&\left. \qquad +\cos ^2\theta \left (2+(2-\alpha)\frac{C_p^2}{V_A^2}+(1+\alpha)\frac{C_s^2}{V_A^2}+2\frac{C_e^2}{V_A^2}\right) \right ] \nonumber \\ && \qquad +k^2\frac{V_A^2}{\Omega _p^2}\left [ 1+(1-\alpha)\left (\frac{C_p^2}{V_A^2}+\frac{C_s^2}{V_A^2}\cos ^2\theta \right) +\frac{C_e^2}{V_A^2}+\alpha \frac{C_s^2}{V_A^2}+\cos ^2\theta \left (1+\alpha \frac{C_p^2}{V_A^2}\right) \right ] \end{eqnarray} \tag{ B5 }$$

$$\begin{eqnarray} && A_6 = -\frac{i}{3}\frac{{\bf k} \cdot {\bf K} \cdot {\bf k} }{\Omega _p}\left \lbrace + k^6\frac{V_A^6}{\Omega _p^6}\cos ^2\theta \left [ \frac{C_e^2}{V_A^2}+\frac{C_s^2}{V_A^2}\right ] + k^4\frac{V_A^4}{\Omega _p^4}\left [ (1-\alpha)\frac{C_s^2}{V_A^2}\left (1+\frac{C_e^2}{V_A^2}\right) +\cos ^2\theta \left (2+2\frac{C_e^2}{V_A^2}+(1+\alpha)\frac{C_s^2}{V_A^2}\right) \right ]\right. \nonumber \\ &&\left. \qquad +k^2\frac{V_A^2}{\Omega _p^2}\left [ 1+\cos ^2\theta +\frac{C_e^2}{V_A^2} +\frac{C_s^2}{V_A^2} (\alpha +(1-\alpha)\cos ^2\theta) \right ] \right \rbrace \end{eqnarray} \tag{ B6 }$$

$$\begin{eqnarray} && A_5 = -\cos ^2\theta \left \lbrace k^8\frac{V_A^8}{\Omega _p^8}\left [ (1-\alpha)\frac{C_e^2 C_s^2}{V_A^4} +\frac{C_p^2}{V_A^2}\left (\alpha \frac{C_e^2}{V_A^2}+\frac{C_s^2}{V_A^2}\right)\right ]\right.\nonumber \\ &&\left. \qquad +k^6\frac{V_A^6}{\Omega _p^6}\left [ (1+\cos ^2\theta)\left (\frac{C_e^2}{V_A^2}+\frac{C_p^2}{V_A^2}+\frac{C_s^2}{V_A^2}\right) +2(1-\alpha)\frac{C_e^2 C_s^2}{V_A^4}+2\frac{C_p^2}{V_A^2}\left (\frac{C_s^2}{V_A^2}+\alpha \frac{C_e^2}{V_A^2}\right) \right ]\right.\nonumber \\ &&\left. \qquad +k^4\frac{V_A^4}{\Omega _p^4}\left [ (1-\alpha)\frac{C_s^2}{V_A^2}\left (\frac{C_e^2}{V_A^2}+\cos ^2\theta \right) +\frac{C_p^2}{V_A^2}\left (2+\frac{C_s^2}{V_A^2}+\alpha \left (-1+\cos ^2\theta +\frac{C_e^2}{V_A^2}\right) \right) +(1+\alpha)\frac{C_s^2}{V_A^2}+2\frac{C_e^2}{V_A^2} +1 \right ] \right \rbrace \end{eqnarray} \tag{ B7 }$$

$$\begin{eqnarray} && A_4 = +\frac{i}{3}\frac{{\bf k} \cdot {\bf K} \cdot {\bf k} }{\Omega _p}\cos ^2\theta \left \lbrace k^8\frac{V_A^8}{\Omega _p^8}(1-\alpha)\frac{C_e^2 C_s^2}{V_A^4} +k^6\frac{V_A^6}{\Omega _p^6}\left [ 2(1-\alpha)\frac{C_e^2 C_s^2}{V_A^4} +(1+\cos ^2\theta)\left (\frac{C_e^2}{V_A^2}+\frac{C_s^2}{V_A^2}\right) \right ]\right. \nonumber \\ &&\left. \qquad +k^4\frac{V_A^4}{\Omega _p^4}\left [ (1-\alpha)\frac{C_s^2}{V_A^2}\left (\frac{C_e^2}{V_A^2}+\cos ^2\theta \right) +1+2\frac{C_e^2}{V_A^2}+(1+\alpha)\frac{C_s^2}{V_A^2} \right ] \right \rbrace \end{eqnarray} \tag{ B8 }$$

$$\begin{eqnarray} && A_3 = +\cos ^4\theta \left \lbrace k^8\frac{V_A^8}{\Omega _p^8}\left [ 2(1-\alpha)\frac{C_e^2 C_s^2}{V_A^4} +2\frac{C_p^2}{V_A^2} \left (\frac{C_s^2}{V_A^2}+\alpha \frac{C_e^2}{V_A^2}\right) \right ]\right. \nonumber \\ &&\left. \qquad +k^6\frac{V_A^6}{\Omega _p^6}\left [ 2(1-\alpha)\frac{C_e^2 C_s^2}{V_A^4}+2\frac{C_p^2}{V_A^2}\left (\frac{C_s^2}{V_A^2}+\alpha \frac{C_e^2}{V_A^2}\right) +\frac{C_e^2}{V_A^2}+\frac{C_p^2}{V_A^2}+\frac{C_s^2}{V_A^2} \right ] \right \rbrace \end{eqnarray} \tag{ B9 }$$

$$\begin{eqnarray} && A_2 = -\frac{i}{3}\frac{{\bf k} \cdot {\bf K} \cdot {\bf k} }{\Omega _p}\cos ^4\theta \left \lbrace k^8\frac{V_A^8}{\Omega _p^8} 2(1-\alpha)\frac{C_e^2 C_s^2}{V_A^4} +k^6\frac{V_A^6}{\Omega _p^6}\left [ 2(1-\alpha)\frac{C_e^2 C_s^2}{V_A^4}+\frac{C_e^2}{V_A^2}+\frac{C_s^2}{V_A^2} \right ] \right \rbrace \end{eqnarray} \tag{ B10 }$$

$$\begin{eqnarray} && A_1 = - \cos ^6\theta k^8\frac{V_A^8}{\Omega _p^8}\left [ (1-\alpha)\frac{C_e^2 C_s^2}{V_A^4} +\frac{C_p^2}{V_A^2}\left (\frac{C_s^2}{V_A^2}+\alpha \frac{C_e^2}{V_A^2}\right) \right ] \end{eqnarray} \tag{ B11 }$$

$$\begin{eqnarray} && A_0 = +\frac{i}{3}\frac{{\bf k} \cdot {\bf K} \cdot {\bf k} }{\Omega _p} \cos ^6\theta k^8\frac{V_A^8}{\Omega _p^8}(1-\alpha)\frac{C_e^2 C_s^2}{V_A^4}. \end{eqnarray} \tag{ B12 }$$