Proton Energization by Phase Steepening of Parallel-propagating Alfvénic Fluctuations

In situ spacecraft measurements show that the solar wind is permeated by Alfvénic fluctuations in which the velocity and magnetic field are correlated mainly in the sense of Alfvén waves propagating away from the Sun (Coleman 1967; Belcher & Davis 1971; D'Amicis & Bruno 2015). Those fluctuations exhibit large amplitudes (comparable to the magnitude of the background magnetic field) and are characterized by negligible density fluctuations and a nearly constant magnetic field magnitude. These properties correspond to spherical polarization (Goldstein et al. 1974; Bruno et al. 2001; Matteini et al. 2015). Alfvénic fluctuations in the solar wind thus display a high degree of coherence that manifests itself not just by the velocity-magnetic field correlation that defines linear shear Alfvén waves, but also via an intrinsic degree of phase coherence among the oscillating fields that is necessary in order to maintain the locally constant intensity constraint (from now on constant-B). Such Alfvénic fluctuations also display typical turbulent features including a well developed energy spectrum and the ubiquitous presence of intermittent structures, that in turn provide suitable places where energy dissipation and particle energization can occur (Marsch 2006; Osman et al. 2010; Tessein et al. 2013; Perrone et al. 2016).

Rotational discontinuities are found at the steepened edges of arc-polarized Alfvén waves, a particular case of constant-B fluctuations in a 1D geometry (Barnes & Hollweg 1974; Tsurutani et al. 2005; Erofeev 2019). It is believed that rotational discontinuities are necessarily generated by wave steepening and the occurrence of abrupt changes in the wave phase (Cohen & Kulsrud 1974; Malara & Elaoufir 1991; Medvedev et al. 1997; Vasquez & Hollweg 2001; Tsurutani et al. 2018; Valentini et al. 2019). Steepened edges, often in the form of rotational discontinuities, have also been observed at the boundaries of switchbacks (Larosa et al. 2020), extremely large-amplitude Alfvénic structures in the magnetic and velocity fields (Horbury et al. 2020). However, how compressible effects are quenched and the nearly constant-B condition is maintained during the dynamical phase wave steepening still remains elusive, especially when the plasma beta is smaller than unity and strong couplings with compressible modes are expected (e.g., Malara & Elaoufir 1991; Jayanti & Hollweg 1993; Roberts & Wiltberger 1995; Malara & Velli 1996). On the other hand, large-amplitude Alfvén waves can survive even in the presence of temperature anisotropy from relativistic Vlasov–Maxwell equations (Barnes & Suffolk 1971). In fact, Alfvénic fluctuations coexist with a plasma that is out of thermodynamic equilibrium (Abraham-Shrauner & Feldman 1977). Proton distribution functions indeed display temperature anisotropies and signatures of preferential perpendicular heating (nonadiabatic expansion), and a stable field-aligned proton beam streaming ahead of the proton core population at the local Alfvén speed (Marsch et al. 1982; Matteini et al. 2013; Verniero et al. 2020).

It is the goal of this Letter to investigate the connection between Alfvén wave steepening, plasma compressibility, and the development of nonthermal features in a low-beta plasma. It has been shown via hybrid simulations that a large-amplitude monochromatic Alfvén wave is subject to parametric instabilities (Araneda et al. 2008; Valentini & Veltri 2009; Nariyuki et al. 2009; Matteini et al. 2010; González et al. 2020). The decay of the wave leads nonlinearly to enhanced proton heating and to a field-aligned beam driven by the field-aligned electric field generated during the wave decay. Here, we consider the more realistic situation of a broadband Alfvénic fluctuation and adopt a hybrid framework complemented by test-particle simulations to investigate the connection between Alfvén wave dynamics and proton energization. In this case, due to the presence of multiple wavelengths and dispersion, the initial fluctuation undergoes phase steepening resulting in the rapid collapse of the wave (Spangler 1989; Buti et al. 2000). We show that phase steepening and wave collapse is accompanied by the formation of rotational discontinuities embedded in compressional structures at steepened fronts. Proton heating and acceleration occur locally at the steepened fronts, and ultimately contribute to the dissipation of compressible fluctuations.

In this work we adopt a hybrid model of the plasma that describes electrons as a massless and isothermal fluid and protons as particles via the (nonrelativistic, quasi-neutral) Vlasov–Maxwell equations:

$$\begin{eqnarray*}&&\begin{array}{l}\displaystyle \frac{\partial {f}_{i}}{\partial t}+{\boldsymbol{v}}\cdot \displaystyle \frac{\partial {f}_{i}}{\partial {\boldsymbol{r}}}+\displaystyle \frac{e}{{m}_{i}}\left({\boldsymbol{E}}+\displaystyle \frac{{{\boldsymbol{v}}}_{i}}{c}\times {\boldsymbol{B}}\right)\cdot \displaystyle \frac{\partial {f}_{i}}{\partial {\boldsymbol{v}}}=0\\ \ \ \displaystyle \frac{\partial {\boldsymbol{B}}}{\partial t}=-c\ {\rm{\nabla }}\times {\boldsymbol{E}},\ \ \ \ {\boldsymbol{J}}=\displaystyle \frac{c}{4\pi }{\rm{\nabla }}\times {\boldsymbol{B}}\\ \ \ {\boldsymbol{E}}+\displaystyle \frac{{{\boldsymbol{u}}}_{i}}{c}\times {\boldsymbol{B}}=-\displaystyle \frac{{k}_{B}{T}_{e}{\rm{\nabla }}n}{{en}}+\displaystyle \frac{{\boldsymbol{J}}\times {\boldsymbol{B}}}{{en}}+\eta {\boldsymbol{J}},\end{array}\end{eqnarray*}$$

with c being the speed of light, e being the electron charge, k_B being the Boltzmann constant, and T_e being the electron temperature. The proton number density n and the proton bulk velocity u _i are computed from the moments of the distribution function (n = ∫f( r , v , t)d v and n u _i = ∫ v f( r , v , t)d v , respectively). Since the model assumes quasi-neutrality, the electric field is obtained through the generalized Ohm's law where contributions from inductive, Hall, electron pressure, and resistive terms are retained while electron-inertia and higher-order terms are not considered.

We have performed 2.5 D numerical simulations by means of the CAMELIA hybrid-PIC code (Matthews 1994; Franci et al. 2018). We have used periodic boundary conditions that are imposed in all directions of the computational box. Lengths are normalized to the proton inertial length d_i = c/ω_p with ${\omega }_{p}={(4\pi {{ne}}^{2}/{m}_{i})}^{1/2}$ the proton plasma frequency. Time is normalized to the inverse of the proton gyrofrequency ${{\rm{\Omega }}}_{{ci}}^{-1}\,={({{eB}}_{0}/{m}_{i}c)}^{-1}$ and velocities to the Alfvén speed ${v}_{A}={B}_{0}/{(4\pi {{nm}}_{i})}^{1/2}$. The plasma beta for both ions and electrons is defined as ${\beta }_{p,e}=8\pi {{nk}}_{B}{T}_{p,e}/{B}_{0}^{2}$. We have included explicit resistivity to improve energy conservation by avoiding energy accumulation at the grid scale, and the corresponding dissipation length scale (l_d ) is related to the Reynolds number and box size (L) through ${R}_{e}\sim {(L/{l}_{d})}^{4/3}$, that is chosen to be greater than the grid scale. We initialize the system with an exact nonlinear solution of the MHD system corresponding to a large scale, non-monochromatic and constant-B Alfvénic fluctuation propagating along the mean magnetic field B ₀, taken along the x-axis. The magnetic field of the wave is given by $\delta {b}_{z}=\delta {b}_{0}\cos (\phi ({k}_{0},x))$ and $\delta {b}_{y}=-\delta {b}_{0}\sin (\phi ({k}_{0},x))$, with ∣δ b ∣ = δ b₀ the amplitude of the wave normalized to the mean magnetic field magnitude B₀. The phase $\phi ({k}_{0},x)={k}_{0}x+\epsilon {\sum }_{m={n}_{i},m\ne {n}_{0}}^{{n}_{f}}\tfrac{{k}_{0}}{{k}_{m}}\cos ({k}_{m}x+{\phi }_{m})$, where ϕ_m is a random phase between [0, 2π). The main wavevector is k₀ = 2π n₀/L and the initial wave satisfies the Walen relation in the dispersionless limit δ u = − (ω₀/k₀)δ b . The wave frequency ω₀ is determined from the normalized dispersion relation ${k}_{0}^{2}={\omega }_{0}^{2}/(1-{\omega }_{0})$ for left-handed circularly polarized waves in parallel propagation. This initial condition corresponds to an initially uniform and isotropic plasma with a broadband Alfvénic fluctuation comprised of modes ranging from k_i = 2π n_i /L to k_f = 2π n_f /L. The parameter controls the deviation from the monochromatic case, which is recovered for = 0 (Malara & Velli 1996; González et al. 2020).

We present results for a squared box simulation of side L = 128 d_i with 1024 particles per cell and 1024² grid points, corresponding to a grid size Δ = d_i /8. We set the time step to move the particles as ${\rm{\Delta }}t=0.01\,{{\rm{\Omega }}}_{{ci}}^{-1}$ and used 30 substeps for the fields. The main wavevector is k₀ d_i = 0.2 and we set = 0.5. This yields an initial wave-packet with a power-law energy spectrum between the interval kd_i = [0.05, 0.5]. The plasma beta is β_p,e = 0.5 and the wave amplitude δ b₀/B₀ = 0.5. Finally, the explicit resistivity is η = 2 × 10⁻³ (in units of $4\pi {v}_{A}{c}^{-1}{{\rm{\Omega }}}_{{ci}}^{-1}$) that corresponds to l_d ∼ 2Δ.

Throughout the text we will use a field-aligned coordinate system and will refer to parallel and perpendicular in terms of the direction of the total magnetic field $\hat{{\boldsymbol{b}}}={\boldsymbol{B}}/| {\boldsymbol{B}}|$. The temperatures are defined in terms of the decomposition of the pressure tensor according to the direction of the total magnetic field as ${p}_{\parallel }={\boldsymbol{p}}:\hat{{\boldsymbol{b}}}\hat{{\boldsymbol{b}}}$ and ${p}_{\perp }={\bf{p}}:({\boldsymbol{I}}-\hat{{\boldsymbol{b}}}\hat{{\boldsymbol{b}}})/2$ and the pressure tensor p = ∫( v − u _i )( v − u _i )f( r , v , t) d³ v is obtained from the particle velocity distribution.

Figure 1 summarizes the time evolution of the system and the top panel shows the rms of the magnetic and velocity field fluctuations, as well as the average parallel and perpendicular proton temperatures. The initial Alfvénic fluctuation collapses after a few tens of proton gyroperiods by releasing its energy to the plasma in the forms of thermal and kinetic energy, the latter via the formation of a field-aligned beam, as we discuss below, until a steady-state is achieved at around $250\ {{\rm{\Omega }}}_{{ci}}^{-1}$.

Zoom In Zoom Out Reset image size

The disruption of the wave and the resulting proton energization and beam formation are due to the phase steepening of the initial wave. Because of dispersive effects, we observe a rapid phase steepening of the wave due to the larger scales catching up with the smaller ones, ultimately leading to a strong modulation of the magnitude of ∣ B ∣ and hence to localized steepened wave fronts, a process analogous to the collapse of localized Alfvénic wave packets due to modulational instability (Spangler 1985; Buti et al. 2000; Matteini et al. 2010). Note that in the case of a soliton the timescale of the disruption scales nearly like $\delta {b}_{0}^{-4}$ (Velli et al. 1999). With a general wave packet, the timescale of the disruption is instead dictated by a combination of nonlinearity and dispersion.

Departures from the initial constant-B state drive compressible fluctuations and a field-aligned electric field at the steepened edges of the Alfvénic fluctuation. The compressive fluctuations are displayed in the bottom panel of Figure 1, where the standard deviation ($\delta (X)=\ \lt {\left(X-\left\langle X\right\rangle \right)}^{2}{\gt }^{1/2}$) of proton density and the field-aligned electric field are plotted. The ratio of the standard deviation of the magnetic field magnitude (σ_F = δ(∣ B ∣)) to the magnetic field fluctuations (${\sigma }_{c}={[{\delta }^{2}({b}_{x})+{\delta }^{2}({b}_{y})+{\delta }^{2}({b}_{z})]}^{1/2}$) is also shown. All these quantities display a maximum during the initial stage of wave collapse/disruption and then minimize at saturation, pointing to the fact that nonlinear wave-particle interactions at the steepened edges reduce plasma compressibility.

In Figure 2 we show a contour plot (spatio-temporal diagram x–t plane ) of the magnitude of the magnetic field and the mean field-aligned electric field at the plane y_p = 64 d_i , showing the characteristics of the discontinuities generated at the wave edges as the result of phase steepening. These structures propagate at nearly the Alfvén speed and essentially along the mean magnetic field. One can observe that the shape of these structures changes as they propagate, and that they largely fade away at saturation (at around $t=300\ {{\rm{\Omega }}}_{{ci}}^{-1}$). The main contribution to the field-aligned electric field at those locations comes from the Hall term ( J × B = B · ∇ B − ∇( B ²/2)), in particular from the field-aligned component of the second term on the right-hand side (Matteini et al. 2010; González et al. 2020). The interaction of protons with these propagating discontinuities and the parallel electric field therein provides a suitable mechanism for the acceleration and heating of protons.

Zoom In Zoom Out Reset image size

The proton velocity distribution function (VDF) in the (v_∥ − v_⊥) plane is shown in Figure 3. The VDF is displayed early in the evolution when the strongest heating is taking place. As can be seen, the largest contribution to the parallel temperature comes from the field-aligned proton beam at the Alfvén speed. The proton beam is generated early in the evolution after the formation of magnetic pressure fluctuations and remains stable throughout the final steady state. A qualitatively similar behavior was found in Nariyuki & Hada (2014) by means of a reduced one-dimensional MHD–Vlasov model in a radially expanding geometry. In their simulations a beam also forms at steepened wave fronts generated by MHD nonlinearities, rather than by dispersive effects.

Zoom In Zoom Out Reset image size

To illustrate the dynamical phase steepening process and its implications on the particles VDF, we computed single-point measurements of some field quantities and the reduced parallel and perpendicular VDF at two different fixed points in the simulation domain. The two probes are located at r₁ = (32 d_i , 64 d_i ) and r₂ = (64 d_i , 64 d_i ), respectively, and correspond to the yellow lines marked in Figure 2. The time series at both probes is presented in Figure 4, where we show the transverse magnetic field fluctuations, the magnitude of the magnetic field, the particle density, and the field-aligned electric field as a function of time. The vertical dashed lines in the top panels of Figure 4 show the location of a single structure that crosses the two points. We estimated that the structure propagates at a speed of V = 0.82 v_A , which is the group speed at the wavenumber k₀, although the characteristics in Figure 2 show that the speed is not constant and some of the steepened fronts break down.

Zoom In Zoom Out Reset image size

The resulting signatures on the proton VDF due to rotational discontinuities on the particle VDF are different for the parallel and perpendicular components. Plume-like structures can be identified in the perpendicular component with enhancement of particles with large v_⊥ at the location at the steepened fronts. On the other hand, the velocity-space structures in the parallel component show the presence of beams in front of the discontinuities. As the faster particles in the beam are farther from the steepened fronts, they arrive at the probe earlier, hence the tilted structure that can be seen in the time series of the parallel VDF. Besides, phase-space holes and slower parallel particle velocity right in the location of the steepened fronts are evidenced.

To understand the energization process produced by the interaction of protons with the rotational discontinuities, we computed the trajectory of test-particles using high-cadence electric and magnetic field data generated by the hybrid simulation. The initial test-particle ensemble is randomly distributed through the simulation box with initial Maxwellian distribution function for perpendicular velocities while particles are initialized with zero velocity along the mean field. Periodic boundary conditions are imposed in x–y directions and we considered protons with the same physical parameters as in the hybrid simulation.

Figure 5 shows the tracking of two representative particles that undergo the acceleration process resulting from the interaction with discontinuities. The top panel shows the field interpolated along the particle trajectory. The middle panel shows the magnetic moment $\mu ={m}_{i}{v}_{\perp }^{2}/2{B}_{0}$, normalized to the magnetic moment evaluated with the initial proton thermal velocity (black line) and the parallel velocity (red line). Proton interaction with rotational discontinuities results in the violation of the first adiabatic invariant. The nonadiabatic particle behavior is observed once the particle crosses a discontinuity, followed by the acceleration of the particle along the magnetic field due to the field-aligned electric field inside the structure. The bottom panel of Figure 5 presents the ratio between the particle gyroradius ρ_i = v_th /Ω_ci and the radius of curvature of the magnetic field at particle location R = 1/κ, with the curvature defined as $\kappa =\parallel \hat{{\boldsymbol{b}}}\cdot {\rm{\nabla }}\hat{{\boldsymbol{b}}}\parallel$. The breaking of particle magnetic moment is observed whenever there is rapid change on (ρ_i /R) rather than the amplitude of the ratio itself. The nonadiabatic and stochastic behavior of protons occurs at discontinuities when particles experience a sharply curved magnetic field rotation at the discontinuities. This is consistent with the signatures on the perpendicular VDF (middle panel of Figure 4), showing an enhancement of particles with larger velocities each time a magnetic structure is crossed. To conclude, the proton heating and acceleration process resulting from the phase steepening of large-amplitude Alfvénic fluctuations is complex. The natural development of parallel-propagating structures that travel at around the Alfvén speed involve a bipolar electric field that accelerates particles into a mean field-aligned beam, while scattering by the magnetic field structure contribute at the same time to an enhancement of perpendicular heating. Particle that are being accelerated into the beam may resonate with the propagating structure leading to the damping of compressive fluctuations and allowing the final nonlinear steady state.

Zoom In Zoom Out Reset image size

In this Letter we made use of hybrid simulations, complemented by test-particle simulations, to investigate proton energization at the phase-steepened edges of an initial large-amplitude, constant-B Alfvénic fluctuation in a low-beta plasma. We find that dispersion leads to the initial phase steepening of the wave resulting in its rapid collapse, as predicted for dispersive Alfvénic wave packets. This process is accompanied by the formation of rotational discontinuities embedded in compressional structures characterized by an enhancement of magnetic pressure at the steepened fronts that propagate at the group speed slightly less than the Alfvén velocity. Proton perpendicular heating via pitch angle scattering and parallel acceleration take place in those localized regions, due to the interaction of protons with the parallel electric field (mainly induced by the gradients of ∣ B ∣) therein. Within the fully self-consistent hybrid simulations, it is those demagnetized protons, accelerated up to about the Alfvén speed, to ultimately mediate the damping of the parallel electric field and reduce compressible fluctuations via nonlinear wave-particle resonance. Our results provide a possible explanation for the ubiquitous presence of a stable, field-aligned proton beam, commonly observed in the Alfvénic wind, and enhanced proton heating. Even though relaxation toward a constant-B state has been predicted also by asymptotic analysis in MHD (Cohen & Kulsrud 1974), we argue that the existence of the proton beam itself is intimately related to the quenching of compressible fluctuations and to the nearly constant-B field that characterizes the Alfvénic wind. MHD simulations so far show strong coupling with compressible modes at low plasma betas that lead to parametric type instabilities, thus preventing such an asymptotic relaxation of large-amplitude fluctuations (Malara & Elaoufir 1991; Roberts & Wiltberger 1995). In any event, it will be of interest to carry out a systematic comparison between long-term evolution of a broadband Alfvénic wave packet in MHD and in hybrid models. On the other hand, our simulation displays a final degree of compressibility of fluctuations comparable to solar wind observations (Villante & Vellante 1982), although some steepened fronts survive or tend to reform. In this respect further investigations on the role of electron physics, not included within the hybrid model, is needed, since electrons may contribute significantly to the damping of compressible fluctuations via resonant interactions (Hollweg 1971).

We would like to thank the anonymous referee for constructive comments and A. Artemyev for useful discussions. This research was supported by NASA grant #80NSSC18K1211 and the HERMES NASA DRIVE Science Center grant #80NSSC20K0604. We acknowledge the Texas Advanced Computing Center (TACC) at The University of Texas at Austin for providing HPC resources that have contributed to the research results reported within this paper. URL: http://www.tacc.utexas.edu.