Effect of Upstream ULF Waves on the Energetic Ion Diffusion at the Earth's Foreshock. I. Theory and Simulation

A collisionless shock is thought to be an efficient particle accelerator in space. One of the most plausible acceleration models of galactic cosmic rays is the diffusive shock acceleration (DSA) model (Bell 1978; Blandford & Ostriker 1978). In that model, turbulent electromagnetic waves that are present both upstream and downstream of the shock scatter some charged particles back and force across the shock. Because the turbulent waves are convected by background flows with different flow speeds upstream and downstream, the particles gain energy through the scattering process. The process occurs stochastically, leading not only to diffusion in momentum space, but also to diffusion in real space. The spatial diffusion has an especially important role in the DSA model. A smaller spatial diffusion coefficient leads to a larger maximum energy attained by the model. This is because that the particles are well-confined near the shock, due to the smaller diffusion, and have a better chance of crossing the shock many times.

Quasi-linear theory (QLT) can provide the spatial diffusion coefficient, or alternatively, mean free path as a function of the wave turbulence spectrum. The spatial diffusion takes place as a result of pitch-angle scattering across 90° in pitch angle. Mathematically, the spatial diffusion coefficient is calculated by integrating the pitch-angle diffusion coefficient. Jokipii (1966) first developed the QLT in a so-called magnetostatic slab turbulence model. In the slab model non-propagating Alfvén wave, i.e., time stationary and purely transverse fluctuating field with wave vector along the ambient field, is considered. That theory combined with the slab model is often referred to as the Standard QLT. After Jokipii (1966), many authors discussed the validity of the QLT, and developed several theories in various turbulence models, such as a purely two-dimensional model and a composite model of 1D slab and 2D geometries (see Shalchi 2009 for details). The QLT has often assumed a simple power-law spectrum of the wave turbulence, and has been successfully applied to problems of diffusion and acceleration of cosmic rays in interplanetary and interstellar media. Near the Earth's foreshock, however, wave spectra have an intense peak rather than a simple power-law (Fairfield 1969; Childers & Russell 1972; Trattner et al. 1994; Kis et al. 2007, 2017).

Quasi-monochromatic waves categorized as terrestrial foreshock ULF waves with a period of about 30 sec are often observed in the Earth's foreshock. A number of studies have discussed the characteristics of the foreshock waves and their associated ion populations, such as field-aligned beams (FABs) and diffuse ions. There are a number of reviews on the physics of foreshock (e.g., Greenstadt et al. 1995; Eastwood et al. 2005; Burgess et al. 2012; Burgess & Scholer 2015). The foreshock ULF waves often lead to a peak in wave spectrum. The quasi-monochromatic 30 s waves are often observed near the foreshock boundary. They are nearly transverse fluctuations, propagating sunward away from the shock, and have intrinsically right-handed polarity in the solar wind frame. Hence, these ULF waves are identified as waves generated by FAB ions backstreaming away from the shock through resonant mode instability (Gary et al. 1981; Winske & Leroy 1984). This wave generation mechanism has been confirmed by in situ observations of low-frequency waves and ion distributions detected by Cluster spacecraft (Mazelle et al. 2003; Meziane et al. 2003; Kis et al. 2007). This fact indicates that the FAB ions generate the monochromatic ULF waves. As one moves deeper into the foreshock, the quasi-monochromatic waves steepen into shocklets with whistler-mode wave trains (Hoppe et al. 1981). Also, during convection, these waves grow in amplitude into pulsations, termed short large-amplitude magnetic structures (SLAMS) (Lucek et al. 2008). According to the changes in waveform, the peak of the ULF power spectrum becomes broader when the spacecraft further approaches the shock (Greenstadt et al. 1995). These ULF waves will efficiently scatter foreshock energetic ions whose parallel velocities are comparable to that of FAB ions.

The diffusion of upstream energetic ions is described by the so-called diffusion-convection equation (e.g., Skilling 1975). In steady state, the diffusion-convection equation tells us that density of energetic particles increases exponentially when approaching the shock from the upstream side; density reaches its maximum at the shock, then keeps a constant value downstream (e.g., Jokipii 1971). Here, the e-folding distance is determined by the balance between upstream spatial diffusion and solar wind convection. This feature was confirmed by a number of in situ observations of Earth's bow shock. Trattner et al. (1994) showed the result of statistical analysis of about 300 events where the AMPTE/IRM satellite observed foreshock energetic ions. The analyzed energy range was 10–67 keV e⁻¹. Later, Kis et al. (2004) and Kronberg et al. (2009) improved a specific event by using multi-spacecraft data. They used data obtained by Cluster satellites in the same period when two of the four satellites separate almost along the upstream magnetic field line. Using the Cluster data, one can estimate the parallel spatial gradient of energetic ion density for a single event. Kis et al. (2004) analyzed the data of lower-energy protons from 10 to 32 keV, whereas Kronberg et al. (2009) analyzed those of higher energy ions from 42 to 123 keV protons and helium.

All these in situ observations showed an exponential profile of upstream energetic ion density, as expected from the classical diffusion-convection theory. Kronberg et al. (2009) demonstrated that the values of e-folding distances obtained by Kis et al. (2004) and Kronberg et al. (2009) lie approximately on the line as a function of ion energy, confirming that these two observations are actually from the same single event. Further, the e-folding distances estimated by Kis et al. (2004) are ∼0.5–2.8 Re for ∼11–27 keV ions—apparently smaller than the values obtained by Trattner et al. (1994). This fact reveals that the efficient scattering of ions takes place in the event analyzed by Kis et al. (2004). The events analyzed by Kis et al. (2004) and Kronberg et al. (2009) contained an intense ULF peak in the spectra. In the event, the peak frequency was 0.03 Hz in the spacecraft rest frame, and the velocity of FAB ions was about 1600 km s⁻¹ in the plasma rest frame, corresponding to a FAB energy of ∼5 keV in the spacecraft rest frame (which is described in Section 5). Hence, diffuse ions whose energies are higher than 5 keV will be affected by the ULF waves generated by the FAB.

In this paper, we theoretically and numerically investigate the influence of a ULF wave peak in the spectrum of magnetic fluctuations to spatial diffusion of energetic ions. In Section 2, we model the spectrum of magnetic fluctuations containing features of ULF waves as a superposition of multi-power-law spectra that have different power-law indices. For the modeled spectrum, the diffusion coefficient of energetic particles is estimated by using quasilinear theory in Section 3 to discuss the qualitative effect of the ULF peak. Test particle simulation is then performed in Section 4 to quantitatively estimate the diffusion coefficient for various parameters. In Section 5, we apply our test particle scheme to the Cluster observation reported in the companion paper (Kis et al. 2017; hereafter referred to as Paper II) for the 2003 February 18th event, to show that the effect of FAB in the scattering of diffuse ions can be significant. We also discuss the relation between the energetic ion diffusion and wave spectral evolution toward the shock. Finally, a summary and our conclusions are given in Section 6.

We consider a wave spectrum composed of N segments with different power-law indices as

$$\begin{eqnarray}&&p(k)=\displaystyle \sum _{j=1}^{N}{p}_{j},{p}_{j}=p({k}_{j-1}\leqslant k\leqslant {k}_{j})={\alpha }_{j}| k{| }^{-{\gamma }_{j}},\end{eqnarray} \tag{ 1 }$$

where p(k) is a dimensionless power spectral density (PSD) of the magnetic field fluctuating in transverse to the background field, and j indicates the index of segment in wavenumber space. The sum (or integration) of p(k) gives total wave energy, ${\sigma }_{\perp }^{2}$, normalized to the background magnetic field, where a symmetry of the spectrum with respect to k = 0 is assumed. Figure 1 shows examples of modeled spectra with N = 5. This type of wave spectrum is indeed observed by Cluster satellite, as will be shown later in Figure 7. The bottom axis is the normalized wavenumber, ${{kv}}_{{\rm{A}}}/{\rm{\Omega }}$, using the Alfvén speed v_A and the ion gyrofrequency Ω. The top axis denotes the corresponding particle parallel velocity, v_r/v_A, which is calculated from the linear resonance condition $-{{kv}}_{r}={\rm{\Omega }}$, where the wave frequency is omitted by assuming v_r ≫ v_A. The minimum and maximum wavenumbers are ${k}_{0}={\rm{\Delta }}k=2\pi {\rm{\Omega }}/6620{v}_{{\rm{A}}}$ and ${k}_{5}=\pi {\rm{\Omega }}/{{v}}_{A}$, respectively, and the wavenumber values at the four break points of p(k) are (k₁, k₂, k₃, k₄) = (31, 55, 85, 118) Δk. In the resonant velocity, the break points are represented correspondingly as (v₀, v₁, v₂, v₃, v₄, v₅) = (1054, 34, 19, 12.4, 8.9, π⁻¹)v_A. The triangular portion in k₁ < k < k₃ corresponds to the ULF peak.

Zoom In Zoom Out Reset image size

Parameters for the modeled PSD include the power-law index of each segment, γ_j, total transverse wave energy, ${\sigma }_{\perp }^{2}(=\sum p(k){\rm{\Delta }}k)$, and the rate of the ULF peak energy to the total one, $\epsilon (=({\sigma }_{\perp }^{2}-{\sigma }_{0}^{2})/{\sigma }_{\perp }^{2})$, where σ₀² is the wave energy except for the triangular portion in k₁ < k < k₃. In the following, a number of wave models are taken into account. Their parameters are summarized in Table 1. The PSDs in Figures 1(a) and (b) correspond to runs C and D, respectively. The different triangles in the peak are for different values of . In the quasilinear theory and test particle simulations, the PSD is used to give the turbulent magnetic field, with the further assumption that each wave has a random phase. In addition, we assume that the wave electric field is negligibly small, leading to a static magnetic field turbulence. This is justified when the wave phase velocity, ∼v_A, is sufficiently smaller than the velocities of energetic ions to be considered. This turbulence model is the so-called magnetostatic slab model, which is often used as a solar wind turbulence model (e.g., Schlickeiser 2002; Shalchi 2009).

Table 1.  Parameters for the Wave Spectrum Used in the Test Particle Simulation

Run	${\gamma }_{1}$	${\gamma }_{2}$	${\gamma }_{3}$	${\gamma }_{4}$	${\gamma }_{5}$	${\sigma }_{\perp }^{2}$ a(${W}_{\perp }$ b)	$\,\epsilon$ %c
A	0.2	−6.8	9.4	0.2	1.5	0.01, 0.1, 1	70
B	0.2	−5.7	8.0	0.2	3.0	0.01, 0.1, 1	70
C1	0.2	0.2	0.2	0.2	1.5	0.2	0
C2	0.2	−3.4	4.9	0.2	1.5	0.29	30
C3	0.2	−5.0	7.1	0.2	1.5	0.40	50
C4	0.2	−6.8	9.4	0.2	1.5	0.67	70
D1	0.2	0.2	0.2	0.2	3.0	0.2	0
D2	0.2	−2.5	3.7	0.2	3.0	0.29	30
D3	0.2	−4.0	5.8	0.2	3.0	0.40	50
D4	0.2	−5.7	8.0	0.2	3.0	0.67	70
E1	0.48	−2.76	4.81	2.45	3.48	0.37 (6.98)	71
E2	−0.24	−6.32	6.39	1.90	4.74	0.50 (9.49)	37
E3	0.73	−0.84	2.09	2.08	2.08	0.57 (13.8)	22

Notes.

^aTransverse wave intensity normalized by the background field. ^bTransverse wave intensity in units of ${10}^{-11}$ erg cc⁻¹. ^cRate of the ULF peak intensity.

Download table as:  ASCIITypeset image

3.1. Spatial Parallel Diffusion Coefficient

We derive a spatial parallel diffusion coefficient with the multi-power-law shape of the wave turbulence spectrum, using the standard quasilinear theory (Jokipii 1971; Schlickeiser 2002; Shalchi 2009; Burgess & Scholer 2015). In that theory, a pitch angle diffusion coefficient, ${D}_{\mu \mu }^{\mathrm{QL}}$, is given by

$$\begin{eqnarray}&&{D}_{\mu \mu }^{{\rm{QL}}}=\displaystyle \frac{\pi {\rm{\Omega }}}{4}(1-{\mu }^{2})\displaystyle \frac{| {k}_{r}| P({k}_{r})}{{B}_{0}^{2}},\end{eqnarray} \tag{ 2 }$$

where μ is the pitch angle cosine, B₀ the ambient magnetic field, ${k}_{r}=-{\rm{\Omega }}/\mu v$ the resonant wavenumber with particle velocity v, and P(k) the transverse wave power spectrum, respectively. Here, P(k) normalized by using B₀ corresponds to (1). The parallel spatial diffusion coefficient, ${\kappa }_{\parallel }^{\mathrm{QL}}$, is obtained by integrating ${D}_{\mu \mu }^{\mathrm{QL}}$ as

$$\begin{eqnarray}&&{\kappa }_{\parallel }^{\mathrm{QL}}=\displaystyle \frac{{v}^{2}}{8}{\int }_{-1}^{1}\displaystyle \frac{{(1-{\mu }^{2})}^{2}}{{D}_{\mu \mu }^{\mathrm{QL}}}d\mu .\end{eqnarray} \tag{ 3 }$$

Because the segment j of the wave spectrum (Figure 1) contributes to ${\kappa }_{\parallel }^{\mathrm{QL}}$ via corresponding μ interval, the integration of μ is divided into N segments as

$$\begin{eqnarray}&&{\int }_{0}^{1}d\mu \simeq \displaystyle \sum _{j=1}^{N}\,{\int }_{{\mu }_{j}}^{{\mu }_{j-1}}d\mu \end{eqnarray} \tag{ 4 }$$

where ${\mu }_{j}={v}_{j}/v$ for v > v_j and μ_j = 1 for v < v_j, and v_j(=Ω/k_j) is the resonant velocity at the break points of p_k. Here, k_j > k_j−1; then v_j < v_j−1, equivalently μ_j < μ_j−1. The minimum pitch-angle cosine, μ_N, is finite because p_k is truncated at the highest wave number, k_N. Figures 2(a) and (b) shows the integration regions in velocity and wavenumber spaces, respectively, for the case of N = 5. Here, the regions corresponding to each other are filled in the same way. First, let us consider the particles with velocity v₁ < v < v₀. In this case, μ₀ = 1 because v < v₀, otherwise μ_j < μ_j−1 < 1 for j = 2–5. The particle with the pitch-angle cosine of μ₁ < μ < μ₀ resonates the wave with the wavenumber of Ω/v < k < k₁ (shown by the shaded region). Similarly, the particles with μ < μ₁ interact with the corresponding resonant waves as shown in Figure 2. Thus, in this case, the particles can resonate with the waves in every segments from 1 to 5.

Zoom In Zoom Out Reset image size

Next, let us consider the particles with the velocity of v₂ < v < v₁. In this case, μ₀ = μ₁ = 1 because v < v₁ < v₀, otherwise μ_j < μ_j−1 < 1 for j = 3–5. Then the integration of the first segment vanishes due to μ₀ = μ₁ = 1, signifying that the particles with v < v₁ cannot resonate with the waves in the first segment, i.e., k₀ < k < k₁.

Finally, the diffusion coefficient κ_∥^QL in (3) reads

$$\begin{eqnarray}&&{\kappa }_{\parallel }^{\mathrm{QL}}=\displaystyle \frac{{B}_{0}^{2}}{\pi }\displaystyle \sum _{j}^{N}\displaystyle \frac{{v}^{3-{\gamma }_{j}}}{{\alpha }_{j}{{\rm{\Omega }}}^{2-{\gamma }_{j}}}\displaystyle \frac{{\delta }_{\max ,j}-{\delta }_{\min ,j}}{(2-{\gamma }_{j})(4-{\gamma }_{j})},\end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}&&{\delta }_{\max ,j}=(4-{\gamma }_{j}){\mu }_{j-1}^{2-{\gamma }_{j}}-(2-{\gamma }_{j}){\mu }_{j-1}^{4-{\gamma }_{j}},\end{eqnarray} \tag{ 6a }$$

$$\begin{eqnarray}&&{\delta }_{\min ,j}=(4-{\gamma }_{j}){\mu }_{j}^{2-{\gamma }_{j}}-(2-{\gamma }_{j}){\mu }_{j}^{4-{\gamma }_{j}}.\end{eqnarray} \tag{ 6b }$$

When ${\mu }_{j-1}={\mu }_{j}=1$, δ_max,j = δ_min,j = 2, the contribution of j th segment vanishes. When γ_j < 2, (5) recovers the standard QLT result as shown in (7).

3.2. Asymptotic Solution

Here, we consider asymptotic solutions of (5) for N = 5. First, let us consider a high-energy particle with v₁ ≪ v ≪ v₀. Here, well-separated k₀ and k₁ are assumed. In this case, μ₀ = v₀/v = 1 otherwise μ_j ≪ 1. When all γ_j < 2, ${\mu }_{0}^{2-{\gamma }_{1}}=1$ and ${\mu }_{j}^{2-{\gamma }_{j}}\ll 1$. Then ${\delta }_{\max ,1}=2$, otherwise ${\delta }_{\max ,j}\to 0$ and ${\delta }_{\min ,j}\to 0$. Hence, the term of j = 1 only contributes ${\kappa }_{\parallel }^{\mathrm{QL}}$, and the asymptotic solution becomes

$$\begin{eqnarray}&&{\kappa }_{\parallel }^{{\rm{QL}}}\simeq \displaystyle \frac{2}{\pi }\displaystyle \frac{{B}_{0}^{2}}{{\alpha }_{1}}\displaystyle \frac{{{\rm{\Omega }}}^{{\gamma }_{1}-2}{v}^{3-{\gamma }_{1}}}{(2-{\gamma }_{1})(4-{\gamma }_{1})},\end{eqnarray} \tag{ 7 }$$

The above solution is consistent with the standard quasilinear theory (e.g., (6.64) of Burgess & Scholer 2015). In general a power-law wave spectrum has the steepest slope at dissipation range, i.e., at the highest wavenumber range. Thus, all γ_j < 2 is satisfied if γ₅ < 2. However, that is not the case with intense peak in the spectrum. In a case such as Figure 2(b), γ₃ > 2, and the asymptotic formula (7) fails.

Second, let us consider a low-energy particle with v₅ < v < v₄. In this case, all μ_j are unity except for μ₅ < 1. Thus, the terms of j = 1–4 vanish, and the asymptotic solution of ${\kappa }_{\parallel }^{\mathrm{QL}}$ becomes

$$\begin{eqnarray}&&{\kappa }_{\parallel }^{{\rm{QL}}}\simeq \displaystyle \frac{{B}_{0}^{2}}{\pi {\alpha }_{5}}\displaystyle \frac{{{\rm{\Omega }}}^{{\gamma }_{5}-2}{v}^{3-{\gamma }_{5}}}{(2-{\gamma }_{5})(4-{\gamma }_{5})}\left(2-(4-{\gamma }_{5}){\left(\displaystyle \frac{v}{{v}_{5}}\right)}^{{\gamma }_{5}-2}\right),\end{eqnarray} \tag{ 8 }$$

where δ_max,5 = 2 and the largest term of δ_min,5 is taken as ${\delta }_{\min ,5}=(4-{\gamma }_{5}){\mu }_{5}^{2-{\gamma }_{5}}$. In (8), if v₅ ≪ v, the last term of the right-hand side (RHS) vanishes for γ₅ < 2. Thus, (8) again recovers the standard QLT without the peak in the spectrum. However, in our case, μ₅ = v₅/v is not infinitesimally small because k₅ is finite; therefore, the last term is not vanished.

Finally, let us consider the case with γ₅ > 2. In this case, ${\mu }_{j}^{2-{\gamma }_{j}}$ is not less than unity anymore, and the term for the largest pitch-angle, i.e., ${\mu }_{5}^{2-{\gamma }_{5}}$ is the largest term. Thus, δ_max,j and δ_min,j have not vanished. The most dominant two terms are

$$\begin{eqnarray}{\kappa }_{\parallel }^{{\rm{QL}}} & \simeq & \displaystyle \frac{2}{\pi }\displaystyle \frac{{B}_{0}^{2}}{{\alpha }_{1}}\displaystyle \frac{{{\rm{\Omega }}}^{{\gamma }_{1}-2}{v}^{3-{\gamma }_{1}}}{(2-{\gamma }_{1})(4-{\gamma }_{1})}\\ & & +\,\displaystyle \frac{{B}_{0}^{2}}{\pi {\alpha }_{5}}\displaystyle \frac{{{\rm{\Omega }}}^{{\gamma }_{5}-2}{v}^{3-{\gamma }_{5}}}{({\gamma }_{5}-2)}{\left(\displaystyle \frac{v}{{v}_{5}}\right)}^{{\gamma }_{5}-2}.\end{eqnarray} \tag{ 9 }$$

When v₁ < v < v₀, the first term arises from the segment j = 1 due to μ₀ = 1, and this term increases with v as ${v}^{3-{\gamma }_{1}}$. On the other hand, the second term arises from the segment j = 5, and this term diverges to infinity if the wave spectrum extends to ${k}_{5}\to \infty$, i.e., ${v}_{5}={\rm{\Omega }}/{k}_{5}\to 0$. This is the so-called 90° scattering problem. In the standard QLT, the pitch-angle diffusion coefficient becomes ${D}_{\mu \mu }\to 0$ for $\mu \to 0$, representing a free streaming motion of a 90° pitch-angle particle. When γ_N < 2, this free streaming motion does not affect ${\kappa }_{\parallel }^{\mathrm{QL}}$ obtained by integration of κ_μμ in μ space. This is because ${\kappa }_{\parallel }^{\mathrm{QL}}$ is mainly determined by the particles with 0° pitch-angle, as seen in (7) and the first term of the RHS in (8). On the contrary, the spectrum becomes steepened as γ_N > 2, ${\kappa }_{\parallel }^{\mathrm{QL}}$ is enhanced due to the free streaming motion of 90° pitch-angle particle, as seen in the second term of the RHS in (9), which shows that the 90° scattering problem causes ${\kappa }_{\parallel }^{\mathrm{QL}}\propto v$.

3.3. Comparison between Theory and Simulation

Figure 3(a) shows the diffusion coefficients obtained by (5) of the QLT and (12) of a test particle simulation, indicated by the solid lines and symbols, respectively. The detail method of the simulation is described in Section 4. The parameters used are from Run A in Table 1 (γ₁ = 0.2 and γ₅ = 1.5). The three lines or data set of symbols correspond to the cases ${\sigma }_{\perp }^{2}=0.01,0.1,1$ from top to bottom, where ${\sigma }_{\perp }^{2}$ indicates the total transverse wave energy normalized by ${B}_{0}^{2}$. Both the theory and simulation demonstrate that the diffusion coefficients anti-correlate with ${\sigma }_{\perp }^{2}$ and that the depression of κ_∥ appears near v = v₁. Most of the particles that have v = v₁ can resonate with the waves corresponding to the ULF peak, as shown by the resonance region indicated by the gray region being the largest for v = v₁ in Figure 2(a). Thus, the presence of a ULF wave peak produces a depression that is, equivalently, the concave shape of ${\kappa }_{\parallel }$. When the wave energy is small enough (${\sigma }_{\perp }^{2}=0.01$), the simulation result is in good agreement with the QLT, as expected. With increasing ${\sigma }_{\perp }^{2}$, the result deviates from the theory. In particular, the deviation becomes large at v ∼ v₁. The test particle simulations give a smaller κ_∥ and larger concave shape than the QLT. In the figure, the two dashed lines denote the asymptotic QLT solutions (7) and (8) for the higher and lower velocities, respectively, for ${\sigma }_{\perp }^{2}=1$. These asymptotic solutions match the theoretical solution (5) except for v ∼ v₁. The slope of the asymptotic solution at the higher velocity is ≈ 3 − γ₁ = 2.8, whereas that at the lower velocity is 1.58, which is slightly larger than $3-{\gamma }_{5}=1.5$ due to the second term in the parenthesis of (8), i.e., due to the truncation of p_k at k₅.

Zoom In Zoom Out Reset image size

Figure 3(b) is the same as 3(a), but γ₅ = 3.0 (Run B). Here, the dashed line denotes the asymptotic solution of (9) for ${\sigma }_{\perp }^{2}=1$. The difference between the QLT and the simulation is more significant and larger than one order at most. In theory, ${\kappa }_{\parallel }^{\mathrm{QL}}$ is proportional to v at v < 60v_A, and the slope gradually increases for increased v. This gradual change of the slope results from the 90° scattering problem, because the problem causes ${\kappa }_{\parallel }^{\mathrm{QL}}\propto v$. On the other hand, in the simulation the slope in ${\kappa }_{\parallel }$ clearly changes regardless of ${\sigma }_{\perp }^{2}$ near v = v₁, and at this velocity the deviation from the QLT becomes the largest again. If the second term of the RHS in (9) is not considered, the slope in the theory becomes 3 − γ₁ = 2.8, indicated by chained line's slope. When ${\sigma }_{\perp }^{2}=0.01$, the slope in the simulation for the higher velocities approaches 2.8. This fact indicates that the discrepancy between the simulation and theory is due to the 90° scattering problem. In other words, the 90° scattering problem in the QLT hides the effect of ULF wave peak on ${\kappa }_{\parallel }^{\mathrm{QL}}$, because the concave shape of ${\kappa }_{\parallel }^{\mathrm{QL}}$ vanishes due to the problem.

In summary, we extended the standard QLT using the multi-power-law wave spectrum to evaluate the spatial parallel diffusion coefficient, ${\kappa }_{\parallel }$. The extended QLT and the test particle simulations are in good agreement when the power-law index at the highest wavenumber, γ, is less than 2 and the wave energy is low. Otherwise, the κ_∥ predicted by the QLT can be an order of magnitude larger than that obtained by the simulation. This difference is caused by the so-called 90° scattering problem. The standard QLT considers the linear cyclotron resonance, $\omega -k\mu v={\rm{\Omega }}$, where the wave frequency, ω, is omitted, and Ω is the proton gyrofrequency. Accordingly, no particles can satisfy the resonance condition if their pitch angles become $90^\circ$, leading to ${D}_{\mu \mu }\to 0$ and ${\kappa }_{\parallel }^{\mathrm{QL}}\to \infty$ at μ = 0. Of course, including the effect of finite ω in the resonance condition eliminates the singularity at μ = 0. However, such waves with ω ∼ Ω are usually heavily dissipated at high wavenumber range. Also, the effect becomes negligibly small for ion velocities much larger than the Alfvén velocity. Other avenues to solve the problem include nonlinear effects, such as mirroring and the resonance broadening due to finite amplitude waves (Hada 2003). Even if the wave amplitude is relatively low, the resonance region significantly broadens (see Figure 5 of Kuramitsu & Krasnoselskikh 2005). Shalchi (2009) has developed second-order QLT, in which the resonance condition, including the effects of finite amplitude waves, is taken into account. They demonstrated that the second-order QLT can explain a parallel mean free path obtained from test particle simulation, despite γ > 2 (see Figures 6.5 and 6.6 of Shalchi 2009). Therefore, it may be worthwhile to apply the second-order QLT to our model, although it is beyond the scope of this paper.

Using the modeled magnetic field spectra shown in Section 2, we follow the motions of test particles. Here, we consider the motions of energetic ions in the solar wind rest frame. The so-called one-dimensional slab model is used as our wave turbulence model. In this system, particle energy is conserved because the electric field is negligible due to the time-stationary field. The wave vectors are parallel and anti-parallel to the x-axis, which is along the ambient magnetic field. The effect of convection by the solar wind will be discussed later.

The following equations are solved for a number of test particles.

$$\begin{eqnarray}&&\dot{{\boldsymbol{v}}}={\boldsymbol{v}}\times {\boldsymbol{b}}.\end{eqnarray} \tag{ 10a }$$

$$\begin{eqnarray}&&\dot{x}={v}_{x}\end{eqnarray} \tag{ 10b }$$

Here, time and space are normalized to the reciprocal of Ω and v_A/Ω, respectively. The magnetic field, ${\boldsymbol{b}}=(1,{b}_{y},{b}_{z})$, is normalized to the ambient field along the x-axis. Purely transverse fluctuating fields are given by

$$\begin{eqnarray}{b}_{y}+{{ib}}_{z} & = & \displaystyle \sum _{{k}_{\min }}^{{k}_{\max }}{b}_{k}\{\exp [i({kx}+{\phi }_{k,p})]\\ & & +\exp [i(-{kx}+{\phi }_{k,m})]\},\end{eqnarray} \tag{ 11 }$$

where ${b}_{k}=\sqrt{p(k){\rm{\Delta }}k/2}$, ${k}_{\min }={\rm{\Delta }}k=2\pi /L$, k_max = π, and system size L = 6620. The boundary condition is periodic, phases ϕ_k,p and ϕ_k,m are random between 0 and 2π, and p(k) is defined in (1).

Initially, mono-energetic ions with speed v and isotropic pitch angles are uniformly placed in space. Trajectories of ions are calculated by integrating (10). The spatial diffusion coefficient along the x axis (parallel diffusion) is evaluated as a function of elapsed timescale, τ:

$$\begin{eqnarray}&&{\kappa }_{\parallel }=\displaystyle \frac{\langle {\rm{\Delta }}{x}^{2}\rangle }{2\tau },\end{eqnarray} \tag{ 12 }$$

where the bracket denotes the ensemble average of 10,000 ions, and Δx is the ion displacement at the timescale, τ. For a sufficiently long timescale, the diffusion coefficient becomes constant, representing classical diffusion. We perform a number of runs with different v for the cases of different wave spectra represented in Table 1. We then obtain the classical diffusion coefficients depending on the velocity, or equivalently, on the energy. In each case, we have 10 runs with different v from $4\sqrt{2}{v}_{{\rm{A}}}$ to $128{v}_{{\rm{A}}}$ with a factor of $\sqrt{2}$ spacing.

4.1. Running Diffusion Coefficient

As an overview, the diffusion coefficients of ions that have $v=32\sqrt{2}{v}_{{\rm{A}}}$ are plotted in Figure 4(a) as a function of τ for Runs C1 ( = 0: dashed line) and C4 ( = 0.7: solid line). The evolutions of spatial displacement for some particles are shown for Run C1 (Figure 4(b)) and Run C4 (Figure 4(c)) up to Ωτ = 300. For both runs, the diffusion coefficients are increasing functions of τ when τ is small, although they show more or less constant values for large τ. Such a diffusive regime appears over a timescale longer than the pitch angle scattering time, τ_scatt, shown in Figure 4(a). At this timescale, the change in μ becomes Δμ ∼ 1, representing an isotropic scattering. Trajectories in Figures 4(b) and (c) indicate a Brownian motion along the ambient magnetic field through the isotropic scattering by the turbulent waves at τ > τ_scatt. The effect of the ULF peak manifests as a diminishing of κ_∥ at large τ and also as a shortening of τ_scatt. In the following, ${\kappa }_{\parallel }{\rm{s}}$ are evaluated on a sufficiently long timescale, Ωτ = 10⁴, for various (or ${\sigma }_{\perp }^{2}$) and γ₅.

Zoom In Zoom Out Reset image size

4.2. Effect of ULF Peak Wave

Let us see the dependence of κ_∥ on the ULF peak intensity by varying for relatively large ${\sigma }_{\perp }^{2}$ while keeping the spectral part of non-ULF peak fixed. Figure 5(a) shows κ_∥ obtained for Run C. Here, γ₅ = 1.5. The different markers correspond to different values of (=0, 0.3, 0.5, 0.7). Because ${\sigma }_{\perp }^{2}$ increases with , κ_∥ is smallest for Run C4 ( = 0.7). When the ULF peak exists ( > 0), κ_∥ denotes a concave shape similar to what was observed in Figure 3. The depression rate, R, is defined as

$$\begin{eqnarray}&&R=\displaystyle \frac{{\kappa }_{\parallel }^{\epsilon =0}-{\kappa }_{\parallel }^{\epsilon }}{{\kappa }_{\parallel }^{\epsilon =0}},\end{eqnarray} \tag{ 13 }$$

and plotted for  = 0.3, 0.5, and 0.7 in Figure 5(b). The maximum depression occurs at a velocity slightly larger than v₁ which is indicated by the vertical line at v = 34 v_A. Here, v₁ denotes the maximum parallel velocity with which a particle can resonate with the ULF peak waves. Note that this v₁ is apparently larger than the velocity of FAB, i.e., v_w = 19 v_A. The depression increases with , leading to the concave shape of κ_∥. However, the maximum depression occurs at the same $v\sim {v}_{1}$, regardless of . From Figure 2(a), particles that have parallel velocities in the region ${v}_{3}\lt {v}_{\parallel }\lt {v}_{1}$ (gray region) can resonate with the waves constituting the ULF peak. Among these, for the particles having v < v₁, the width of the resonant pitch angle (gray region) increases with increased v. On the other hand, for particles that have v > v₁ (stripe region), the width of resonant pitch angle decreases with increased v. This is why the maximum depression occurs at v ∼ v₁.

Zoom In Zoom Out Reset image size

Figure 6 is the same as Figure 5, but for Run D (γ₅ = 3.0). The basic feature is similar to that of Run C. However, we can confirm that κ_∥ is smaller (larger) for larger (smaller) v, compared to Run C. As seen in Figure 1, wave power for k < k₄ (k > k₄) in Run D is higher (lower) than that in Run C. Therefore, the particles with higher (lower) energy are easier (harder) to scatter.

Zoom In Zoom Out Reset image size

We apply our test particle scheme to understand the diffusion of upstream energetic ions observed by Cluster spacecraft reported in Paper II for the 2003 February 18th event. This event contains the high-intensity waves in the foreshock that are produced by the FAB. Figure 7 denotes the transverse wave power spectrum using data in the period of 18:30 UT–21:30 UT; this period is the main part of the data analyzed in Paper II. In the panels, from left to right, the spacecraft approaches the bow shock from the upstream side as time elapses: (a) 18:30 UT, (b) 21:06 UT, and (c) 21:30 UT, corresponding to distances of (a) 4.5 R_e, (b) 1.9 R_e, and (c) 1.3 R_e from the shock surface. Here, the distance of the spacecraft from the bow shock is described in Paper II. The average magnitudes of the magnetic fields are (a) B₀ = 6.5 nT, (b) B₀ = 6.5 nT and (c) B₀ = 7.8 nT, corresponding to ion gyrofrequencies f_i = 0.10 Hz and 0.12 Hz, respectively. Other solar wind parameters include plasma density ∼3.0 cm⁻³, solar wind speed v_sw ∼ 650 km s⁻¹, the Alfvén speed v_A = 74 km s⁻¹, the angle between the solar wind velocity and the magnetic field ψ ∼ 147°, the shock normal angle Θ_Bn = 20°, and the angle between the solar wind velocity and the shock normal Θ_sw,n = 231, respectively, as described in Paper II and Kis et al. (2004). The bottom axis indicates the wave frequency in the spacecraft rest frame, f_sc. The top axis denotes the wavenumber, k, calculated by using the Doppler shift relation,

$$\begin{eqnarray}&&\displaystyle \frac{{{kv}}_{{\rm{A}}}}{2\pi }=f=\displaystyle \frac{{f}_{\mathrm{sc}}}{1+{v}_{\mathrm{sw}}\cos \psi /{v}_{{\rm{A}}}},\end{eqnarray} \tag{ 14 }$$

where f is a wave frequency in the plasma rest frame. We assume a linear dispersion relation of the Alfvén wave propagating along the ambient magnetic field.

Zoom In Zoom Out Reset image size

It is clear that the wave spectrum is not time-stationary as the spacecraft approaches the shock, which equivalently represents the spatial variation of the wave spectrum. The spectrum clearly contains a peak at  f_sc ∼ 0.03 Hz, which corresponds to a wave frequency of f = −0.05 f_i in the solar wind rest frame. In a resonant mode instability, right-handed Alfvén waves, propagating in the same direction as the FAB are generated via cyclotron resonance. The resonance condition gives the resonant frequency in the plasma rest frame as

$$\begin{eqnarray}&&{f}_{\mathrm{res}}=\displaystyle \frac{{v}_{{\rm{A}}}{f}_{i}}{{v}_{{\rm{A}}}-{v}_{r}},\end{eqnarray} \tag{ 15 }$$

where v_r is the resonant parallel velocity. Using (14) and (15), the FAB velocity corresponding to f_sc = 0.03 Hz becomes v_r ∼ 1600 km s⁻¹ in the plasma rest frame, which roughly matches the FAB velocity of 900 km s⁻¹ in the spacecraft rest frame observed at ∼8 Re from the shock (Kis et al. 2007). Thus, the observed ULF waves are probably generated by the FAB ions. It can also be seen that the ULF peak is most remarkable at (a). As time elapses, the peak gradually decreases while the wave energy increases, with a frequency higher than the ULF.

The above three spectra are modeled by the multi-power law spectra, as shown by the colored solid lines in Figure 7. On each line, four break points are indicated by the open circles. The measured power-law index (γ) of each segment, the wave energy (${\sigma }_{\perp }^{2}$) normalized by the background field energy, and the rate of the ULF peak energy () are summarized in Table 1 as Run E. In that table, the wave energy (W_⊥) is displayed in units of erg cc⁻¹. From Run E1 to E3, ${\sigma }_{\perp }^{2}$ or W_⊥increase while decreases. Note the ratio of ${\sigma }_{\perp }^{2}$ in run E1 to that in run E3 is different from that of W_⊥, because the background field energy in Run E3 is 1.3 times larger than that in Run E1.

We numerically evaluate κ_∥ by performing additional test particle simulations for Runs E1 to E3. In evaluating κ_∥, here, care should be taken with the timescale τ. The simulation frame is the solar wind rest frame where the solar wind velocity is along the background field. Namely, we are looking at the particles scattered along the magnetic field line co-moving with the solar wind. Those particles may start being scattered when the solar wind enters the foreshock region, and they are observed by the spacecraft on a timescale, τ_conv. To realize the steady-state balance between the diffusion and convection, the particles should be scattered isotropically with the timescale, τ_scatt, before they are convected by the solar wind over timescale τ_conv. Namely, τ_scatt < τ_conv should be satisfied, and κ_∥ attains the classical regime at τ = τ_conv. Figure 8 shows τ_conv and τ_scatt as a function of the ion energy. From Figure 11 of Kis et al. (2007), the distance from the spacecraft to the foreshock boundary along the magnetic field is L_up ≈ 20 Re, leading to the convection timescale of ${\tau }_{\mathrm{conv}}={L}_{\mathrm{up}}/{v}_{\mathrm{sw},\parallel }\,\approx 200\,{\rm{s}}$ in the rough estimate. This corresponds to Ωτ_conv ≈ 120. The scattering timescale, τ_scatt, is obtained from numerical simulations as shown in Figure 4 (a). When E_sc < 42 keV, τ_scatt < τ_conv are almost always satisfied; thus, the diffusion concept is valid. On the other hand, τ_scatt for E_sc = 123 keV is slightly larger than τ_conv; this indicates that the particles leave the foreshock and enter the downstream region before being scattered sufficiently. We use classical values of κ_∥ to evaluate the e-folding distance, as assumed in the past studies (Trattner et al. 1994; Kis et al. 2004; Kronberg et al. 2009), although the diffusion approximation is doubtful at 123 keV ions. We discuss this uncertainty later.

Zoom In Zoom Out Reset image size

Figure 9 summarizes the numerical results of e-folding distances as compared to the Cluster observations (Paper II; Kronberg et al. 2009). The numerically obtained values for Runs E1–E3 are plotted, shown by the colored markers, versus the ion kinetic energy in the spacecraft rest frame. In the figure, observationally obtained values are superimposed by the solid squares. The four solid squares from the lowest energy and two solid squares from the highest energy indicate the values estimated by Paper II and Kronberg et al. (2009), respectively. Note that the data from Paper II shown here are slightly modified from those of Kis et al. (2004).

Zoom In Zoom Out Reset image size

In the diffusion-convection equation in the normal incident frame, the balance between the diffusion and convection terms yields L_n = κ_n/v_sw,n, where ${L}_{n}={L}_{\parallel }\cos {{\rm{\Theta }}}_{\mathrm{Bn}}$ the e-folding distance, v_sw,n = v_sw cosΘ_sw,n the solar wind velocity, and ${\kappa }_{n}={\kappa }_{\parallel }\cos {{\rm{\Theta }}}_{\mathrm{Bn}}^{2}$ the diffusion coefficient. The subscript "n" denotes the direction normal to the shock surface. Next, the e-folding distance along the background field, L_∥, is numerically estimated by

$$\begin{eqnarray}&&{L}_{\parallel }={\kappa }_{\parallel }/{v}_{\mathrm{sw},\parallel };\end{eqnarray} \tag{ 16 }$$

where ${v}_{\mathrm{sw},\parallel }={v}_{\mathrm{sw}}\cos {{\rm{\Theta }}}_{\mathrm{sw},n}/\cos {{\rm{\Theta }}}_{\mathrm{Bn}}$ is solar wind velocity in the de Hoffmann–Teller frame where a motional electric field vanishes across a shock (e.g., Sections 2.5 and 6.3 of Burgess & Scholer 2015). The particle energy in the spacecraft rest frame is then written as ${E}_{\mathrm{sc}}={m}_{p}({v}^{2}+{v}_{\mathrm{sw},\parallel }^{2})/2$, where v is the particle speed in the simulation frame and m_p is the proton mass. Here, an isotropic velocity distribution of ions in the solar wind frame is assumed.

Runs E1, E2, and E3 correspond to the red circles, blue triangles, and green squares, respectively. These simulation values are close to each other when E_sc increases, and the observed value for E_sc = 123 keV matches the simulation result. For E_sc < 42 keV, the observed L_∥ with the lower energy matches the simulated L_∥ with the PSD closer to the shock. This indicates that the energy dependence of L_∥ is closely affected by the spatial variation of the wave spectrum.

Let us now discuss each of these runs in more detail. The red arrows indicate the energies of 0° pitch-angle particles resonant with the lowest and peaked frequency in the ULF peak, shown by the red solid circles in Figure 7(a). In the energy range between the two arrows, the values of Run E1 are depressed due to the presence of the intense ULF peak (i.e.,  = 71%). The observed L_∥ at 42 keV matches with the depressed L_∥. In addition, the resonant frequency range for 123 keV is indicated by longer, horizontal bars in Figure 7. The range is obtained from (14) and (15). The minimum resonance frequency corresponds to a 0° pitch-angle of the 123 keV ion, while the maximum one is simply fixed as f_res = −f_i, corresponding to a 88° pitch-angle of the 123 keV ion, due to the negligible wave energy for the higher-frequency part. It can be seen that the resonance range fully covers the intense ULF peak. Hence, the particles with 42 and 123 keV may be efficiently and directly scattered by the ULF peak waves.

On the other hand, the observed values for relatively low energies (<21 keV) are in good agreement with run E3 using the PSD (Figure 7(c)) in the vicinity of the shock. The resonant frequency range for 11.5 keV is indicated by shorter, horizontal bars in Figure 7, the same as in the case for 123 keV. The minimum resonance frequency corresponds to a 84° pitch-angle of the 11.5 keV ion. This frequency range corresponds to segment 4 in the spectrum. In this range, the resonant wave energy increases from (a) 1.85 × 10⁻¹¹ erg cc⁻¹ to (c) 7.5 × 10⁻¹¹ erg cc⁻¹, and the rate of the resonant wave energy against ${\sigma }_{\perp }^{2}$ increases from (a) 26% to (c) 54%. This fact is approximately equivalent to a hardening of the spectrum, i.e., decreased γ₄, toward the shock, while the ULF peak rate, , decreases. This implies that the intense ULF wave energy transfers to a frequency higher than the ULF peak, which may result in the smaller e-folding distances at E_sc < 21 keV.

In summary, we used three different ULF wave spectra observed at different positions from the shock, and demonstrated that the spatial dependence of the wave spectrum significantly affects the e-folding distance of the diffuse ions. Kronberg et al. (2009) found that the e-folding distance of high-energy ions increases approximately linearly from lower energies as energy per charge is increased, whereas our result show a deviation from this linear dependence on the energy. Our result suggests that the spatial diffusion of the foreshock energetic ions does not obey the simple picture that the standard QLT predicts, and that the spatial or temporal variation of the wave spectrum should be taken into account. In our calculation to evaluate the diffusion coefficient as a function of time, we assumed that the wave spectrum does not change along the magnetic field. To include the spectral variation in more detail, particles should experience the evolving foreshock ULF waves as they travel toward the shock during the convection.

The sequence of the wave spectra detected by the spacecraft shows not only a gradual change but also some small variations in time. We carefully chose the time intervals in which the observed spectra explain the observed e-folding distances well, and obtained good agreement with the simulation and observations (Figure 9). We emphasize here that the wave spectrum far upstream (>2 Re from the shock) never reproduces the small e-folding distance for the diffuse ions with their energies <21 keV. This implies that the small e-folding distances at the lower energy range are caused by the developed ULF power spectrum near the shock, although their association with the shocklet and SLAMS are not clear. Our modeled waveforms are different from these structures due to the random phase approximation.

In our model, we have assumed a balance between convection and diffusion in a steady state. This would be appropriate for the kind of infinite planar shock acceleration described by Bell (1978) and Blandford & Ostriker (1978). However, the Earth's bow shock has a large curvature, so the process might be modified from a simple DSA model assuming one spatial dimension. Early ISEE spacecraft observations demonstrated that the particles far upstream of the Earth's bow shock move essentially scatter-free from the shock (Scholer et al. 1980). This indicates that the wave activity is not enough to scatter the particles, and that the diffusion concept breaks down far upstream. Terasawa (1981) modifies the ideal DSA model by introducing the scatter-free effect in the terrestrial bow shock. Eichler (1981) also theoretically discussed an effect of particle escape from the shock acceleration region. We now discuss the applicability of the ideal DSA to our model. Terasawa (1981) introduced a dimensionless parameter, $\zeta \equiv {v}_{\mathrm{sw},\parallel }{L}_{\mathrm{up}}/{\kappa }_{\parallel }={L}_{\mathrm{up}}/{L}_{\parallel }$, to characterize the effect of escape far upstream along the magnetic field. Beyond the distance, L_up, from the shock in the upstream region, the particles escape freely due to the absence of turbulent magnetic field. When ζ < 1 (weak scattering), particle escape becomes significant, and the diffusion concept may break. In our case, ζ ∼ 16 for 11.5 keV, whereas ζ ∼ 1.3 for 123 keV. This implies that the scattering at 123 keV is not sufficient to allow us to neglect the escape effect, and that the acceleration must be terminated somewhere around this energy.

Typical timescale to escape far upstream is written as ${\tau }_{\mathrm{esc},\parallel }={L}_{\mathrm{up}}^{2}/{\kappa }_{\parallel }=\zeta {\tau }_{\mathrm{conv}}$, where ${\tau }_{\mathrm{conv}}={L}_{\mathrm{up}}/{v}_{\mathrm{sw},\parallel }$. Also, a typical timescale to escape across the magnetic field from the shock region is approximated as ${\tau }_{\mathrm{esc},\perp }={a}^{2}/{\kappa }_{\perp }$, where a(=6 × 10⁴ km) is the radius of curvature of the bow shock, and ${\kappa }_{\perp }$ is the cross-field diffusion coefficient. The diffusion coefficient is determined from the result obtained by the test particle simulation using a non-compressional two-dimensional turbulence model (Otsuka & Hada 2009). From the above, the diffusion coefficient is written as ${\kappa }_{\perp }\sim \rho {{vK}}^{2}$ for K < 1, where the dimensionless parameter $K={{bL}}_{\perp }/\rho$, and the wave amplitude $b={\sigma }_{\perp }$. Here, the wave amplitude, σ_⊥, from Table 1 is used, and a perpendicular scale length of the turbulence is roughly assumed to be L_⊥ ≃ 1000 km from the foreshock wave observations (Narita & Glassmeier 2010). In the present study, K < 1 throughout the energy range from 11.5 to 123 keV. Accordingly, ${\tau }_{\mathrm{esc},\perp }{\rm{\Omega }}\sim ({a}^{2}/{L}_{\perp }^{2})/{b}^{2}$.

In Figure 8, the escape timescales are plotted as a function of the particle energy, E_sc. First, ${\tau }_{\mathrm{esc},\perp }$ is at least 8000 seconds—much larger than τ_conv. Hence, the escape due to the cross-field diffusion is negligible, and the particles are almost tied to the field lines during the convection. Second, τ_esc,∥ decreases for increased E_sc and drops below τ_scatt at E_sc = 123 keV. This implies that the effect of escape far upstream is not negligible, and the isotropic velocity distribution is not guaranteed. On the contrary, our simulation assumed a balance between convection and diffusion, which leads to the isotropic velocity distribution. Nevertheless, our resultant e-folding distances at E_sc = 123 keV are in good agreement with those obtained by Kronberg et al. (2009). We roughly assumed the free escape boundary as the ion foreshock boundary with L_up = 20 Re from the shock. However, the boundary must be fluctuating in time. Also, Terasawa (1981) described L_up ∼ 100 Re, and Scholer et al. (1980) predicted L_up < 200 Re from the in-site observations. For the case of L_up > 20 Re, ${\tau }_{\mathrm{esc},\parallel }$ and τ_conv is larger than τ_scatt for all the energies, and the diffusion concept is still acceptable at E_sc = 123 keV. Further investigation to precisely determine the free escape boundary is needed, because it is closely related to the ability of the DSA in the Earth's bow shock.

For the lower energy range at E_sc < 42 keV, the ideal DSA process is reliable. The e-folding distances of this range estimated by Trattner et al. (1994) are a few times larger than our result or the results given by the analysis in Paper II. There are two possible reasons for that. First, transverse wave energies are apparently smaller in the events analyzed by Trattner et al. (1994), although the PSD clearly contains the ULF peak from their Figure 3. In most of the events, they treated the transverse wave intensity as W_⊥ < 1 × 10⁻¹¹ erg cc⁻¹ (see Figure 5 of Trattner et al. (1994)). For the event reported here, however, W_⊥ < 13.8 × 10⁻¹¹ erg cc⁻¹ (Run E). Second, the power-law spectrum for frequencies higher than the ULF peak is harder for the event treated here. In Run E3, the index is γ₄ ≈ 2.0, while the corresponding index in Trattner et al. (1994) is roughly estimated as ∼3 from their Figure 3.

The cause of the variation of wave spectra should also be discussed. In the event we focus upon here, the transverse wave energy increases as it approaches the shock. This is probably due to inhomogeneity of the background plasma. It is common that wave amplitude generated far upstream of the shock increases while convected by the solar wind toward the shock. We also found that the power-law wave spectrum in segment 4 becomes harder as it approaches the shock. Two possibilities are considered here. First, nonlinear evolution of the large amplitude FAB generated waves contributes to the hardening of the high-frequency part of the wave spectrum. When the amplitude of ULF waves becomes sufficiently large with respect to the ambient magnetic field, nonlinear effects (such as parametric instabilities) often lead to the cascade of wave spectrum (Bruno & Carbone 2013, and references therein). As shown in Figure 7, the wave spectrum near the foreshock boundary has a rather large ULF peak (Figure 7(a)). As time has passed, relative intensity of the ULF peak has decreased (Figures 7(b) and (c)). Instead, relative intensity of the higher-frequency waves corresponding to segment 4 have increased from Figures 7(a) to Figure 7(c), accompanied by a hardening of this part of the power-law spectrum. This scenario may not work when the amplitude of ULF waves is not large enough, as in the cases discussed by Trattner et al. (1994). Second, wave generation by the diffuse ions themselves may also amplify the wave energy, with a frequency higher than the ULF peak, as discussed in Paper II. When the diffuse ions are effectively scattered by the ULF peak waves, the DSA process occurs efficiently, resulting in an increase of the diffuse ion density. Because the diffuse ions have an isotropic velocity distribution in the plasma rest frame, they can generate the waves with frequency higher than the ULF peak. This is a conventional self-consistent model of particle acceleration and wave excitation (Lee 1982; Burgess & Scholer 2015, and references therein). Obliqueness of the wave propagation and SLAMS near the shock may also contribute to the variation of the spectrum. To verify these possibilities, one must perform more advanced analysis of wave evolution in the foreshock and subsequent scattering of the energetic particles, which is a future issue.

We have discussed field-aligned diffusion of energetic ions upstream of the Earth's bow shock. In particular, we have focused on the effect of the foreshock ULF waves with frequencies <0.1 Hz. The waves are modeled as a superposition of a number of non-propagating Alfvén waves that have a multi-power-law spectrum containing a ULF peak with random phases. We extended the standard QLT using the modeled multi-power-law wave spectrum to formulate a spatial parallel diffusion coefficient, ${\kappa }_{\parallel }^{\mathrm{QL}}$. When ${\kappa }_{\parallel }^{\mathrm{QL}}$ is plotted as a function of particle velocity or energy, a concave shape appears due to the presence of the ULF wave peak. However, the QLT does not give reliable results if transverse wave power with respect to the ambient field becomes large, or if the power-law index at the highest wavenumber is above a critical value, γ > 2. The latter is caused by the so-called 90^◦ scattering problem in the QLT. Test particle simulation is performed for a variety of parameters by using the modeled multi-power-law wave spectrum. Qualitative features predicted by the QLT are confirmed. In addition, it is revealed that the lowest frequency of the ULF peak determines the particle energy at which the maximum depression of κ_∥ occurs due to the peak. The observed e-folding distances of the diffuse ions reported in Paper II of this series and Kronberg et al. (2009) are reproduced by the test particle simulation scheme developed here. The evolution of the ULF power spectrum toward the shock affects the energy dependence of the e-folding distance. The e-folding distances of 42 and 123 keV ions are explained by considering the effect of direct scattering by the observed ULF peak waves. On the other hand, the efficient scattering for the ions of 10–32 keV is due to their high intensity, as well as the hardening of the power-law wave spectrum at a frequency higher than the ULF peak in the near vicinity of the shock.

This work was partially supported by a Grant-in-Aid for JSPS Research Fellow 15J40063 (F.O.), JSPS Bilateral Joint Research Project (S.M.), the NKFIH NN 116446 (A.K.), and the HAS-JSPS NKM-95/2016 bilateral project (A.K.). We also offer our thanks to the CSA (Cluster Science Archive) for providing the magnetic field data.