Laboratory measurements of the physics of auroral electron acceleration by Alfvén waves

While the aurora has attracted attention for millennia, important questions remain unanswered. Foremost is how auroral electrons are accelerated before colliding with the ionosphere and producing auroral light. Powerful Alfvén waves are often found traveling Earthward above auroras with sufficient energy to generate auroras, but there has been no direct measurement of the processes by which Alfvén waves transfer their energy to auroral electrons. Here, we show laboratory measurements of the resonant transfer of energy from Alfvén waves to electrons under conditions relevant to the auroral zone. Experiments are performed by launching Alfvén waves and simultaneously recording the electron velocity distribution. Numerical simulations and analytical theory support that the measured energy transfer process produces accelerated electrons capable of reaching auroral energies. The experiments, theory, and simulations demonstrate a clear causal relationship between Alfvén waves and accelerated electrons that directly cause auroras.


Supplementary Methods 1: Reproducing Auroral Magnetospheric Conditions in the Laboratory.
While absolute parameters of the LAPD vary dramatically from those of the auroral magnetosphere [1], the experimental values of B 0 , n e , and T e were chosen so that magnetic pressure dominates over plasma thermal pressure as it does in the auroral magnetosphere. Specifically, the plasma parameters in the acceleration region of the auroral magnetosphere lead to an electron plasma beta that is less than the electron-to-ion mass ratio, β e < m e /m i . Under these extremely low-β e conditions, the electron thermal velocity is slower than the Alfvén velocity, which can be expressed by v te /v A = (β e m i /m e ) 1/2 < 1, and so any resonant acceleration of electrons by the Alfvén waves is expected to occur in the suprathermal tail of the electron velocity distribution. The experiments presented here were performed in hydrogen plasma with m i /m e = 1836 and electron plasma beta β e = 6.69 × 10 −5 , leading to the ratio of the electron thermal to the Alfvén velocity of v te /v A = 0.35, as shown in Supplementary Table 1.
Thermal electrons in the LAPD are much more collisional than those of the auroral magnetosphere, with a Braginskii thermal electron collision rate of ν te = 4.4 × 10 6 s −1 . Through the linear dispersion relation for inertial Alfvén waves that includes the effect of collisions [2,3], given by Supplementary Equation (1), a non-zero collision rate causes an increase in Alfvén wave damping and a slight decrease in the wave frequency. In our experiment, the ratio of the thermal electron collision frequency to the angular frequency of the Alfvén waves is ν te /ω = 0.6, indicating that the electrons in the thermal core of the distribution are moderately collisional. Despite the effects of thermal collisions on the inertial Alfvén wave propagation and damping, the electrons that are resonantly accelerated by E z of the inertial Alfvén wave have resonant parallel velocities well above the electron thermal velocity, v r /v te = ω/k z v te > 1, and therefore these suprathermal electrons are much less strongly affected by Coulomb collisions. The collision frequency is inversely proportional to the third power of the particle velocity, so the effective collision frequency for resonant electrons relative to the thermal electron collision frequency scales as ν re /ν te = (v te /v r ) 3 . For the inertial Alfvén waves in this experiment, we calculate that ω/k z v te = v r /v te = 2.4, giving a normalized collision frequency for the resonant electrons of ν re /ω = 0.04, so that the effect of collisions in the LAPD on the electron velocity distribution g e1 (v z ) at the resonant velocity is minimal.  Analysis of the linearized cold plasma two-fluid system of equations, in the inertial regime v te /v A < 1, shows the existence of a plane wave solution for inertial Alfvén waves with the dispersion relation [2,3] where ω ci is the ion gyrofrequency, k 2 ⊥ = k 2 x + k 2 y is the squared magnitude of the wave vector perpendicular to B 0ẑ , δ e is the electron skin depth, and ν te is the Braginskii collision rate for thermal electrons. Thermal collisions introduce damping. Both finite frequency ω/ω ci and finite perpendicular wavenumber k ⊥ δ e act to decrease the phase speed of inertial Alfvén waves below v A . In the relevant limit k z ≪ k ⊥ for our experiment, the group velocity is largely aligned with B 0 , so that whatever wave structure exists in the perpendicular x-y plane is primarily translated along B 0 . However, the perpendicular structure of inertial Alfvén waves is not inconsequential: the dependence of the parallel phase velocity ω/k z and parallel group velocity ∂ω/∂k z on k ⊥ leads to dispersive behavior. We model the inertial Alfvén waves in this experiment as a superposition of plane waves with different k ⊥ values. Following Fourier decomposition, measurements of the Alfvén wave are found to be consistent with the dispersion relation above.
The Sigma antenna we constructed for this experiment excites Alfvén waves by oscillating currents in insulated copper antenna elements that are immersed in the plasma. Antenna elements are connected by low-impedance feeds to external amplifiers. For these experiments, amplifiers drive the antenna at 1.177 MHz. Using three elements spaced by 2.5 cm alongx allows the antenna to impose structure on the Alfvén wave fields perpendicular to B 0 . At the (x, y) location of g e (v z ) measurements, the Alfvén wave is mostly polarized so that B ⊥ ≈ B yŷ . Consequently, the proceeding overview of the wave analysis will focus primarily on B y measurements. However, the full analysis uses both B y and B x whenever appropriate, e.g. in the calculation of E z .
Electromagnetic fields of the launched wave B ⊥ and E ⊥ perpendicular to B 0 are measured using Elsässer probes [4]. Each probe has B-dot coils to record changes in B x and B y as well as double probes to record E x and E y . In the analysis presented here, the Elsässer variables [5] are not calculated from the measured wave fields. Instead, E ⊥ and B ⊥ measurements from the two Elsässer probes are analyzed to determine experimentally the phase speed v ph = ω/k z of the launched wave modes and compared with predictions from the inertial Alfvén wave dispersion relation in Supplementary Equation (1).
A time slice of the composite B y (x, y, z = 1.9 m, t) measured by Elsässer probe 1, called B y1 for brevity, is shown in Supplementary Fig. 1(a). Similarly, we call B y (x, y, z = 4.5 m, t) measured by Elsässer probe 2 B y2 . The amplitude of B y1 is 0.015 mT. A Fourier transform over x and y givesB y1 (k x , k y , t). The absolute value |B y1 | of the transformed data from Supplementary Fig. 1(a) is plotted in Supplementary Fig. 1(b). This shows the concentration of wave power is highest in the plane wave modes (k x , k y ) = (±50, ±19) m −1 . The values for phase velocity, parallel wave vector k z , and parallel wavelength λ z for the plane wave mode (k x , k y ) = (50, 19) m −1 given below are typical for the waves produced at this frequency under these plasma conditions, although it should be noted that we do not approximate that wave as consisting of a single (k x , k y ) plane wave mode here in the analysis of the experimental measurements. All plane wave modes with measured signal above the noise floor are included whenever appropriate in the analysis of the experimental measurements.
The phase delay for a given mode (k x , k y ), attributable to the time of flight of the wave between the two Elsässer probes, is found by comparing the time seriesB y1 andB y2 . The delay for (k x , k y ) = (50, 19) m −1 is shown in Supplementary Fig. 1(c). Elsässer probe 2, farther from the antenna, records the wave arriving later in time and with diminished amplitude. Phase fronts arriving at probe 2 are delayed by 0.77 ± 0.03 µs. Given the probe spacing of z 2 − z 1 = 2.56 ± 0.05 m, the experimentally determined parallel phase speed for this mode is v ph = (3.32 ± 0.14) × 10 6 m s −1 . This experimental value, along with the definition of the parallel phase speed v ph = ω/k z and the known ω driven by the antenna amplifiers, gives an experimental value of k z = 2.23 ± 0.09 m −1 and a parallel wavelength of λ z = 2.82 ± 0.11 m.
To verify that the launched wave mode with (k x , k y ) = (50, 19) m −1 is an inertial Alfvén wave, the experimental values of v ph , k z and λ z are compared with predictions from Supplementary Equation (1) using v A = 3.4×10 6 m s −1 , ω ci = 1.6×10 7 rad s −1 , δ e = 5×10 −3 m, and ν te = 4.4×10 6 s −1 . The predicted values are v ph = (2.9±0.3)×10 6 m s −1 , k z = 2.5±0.3 m −1 , and λ z = 2.5 ± 0.3 m. Most of the uncertainty in these values comes from n e0 . The experimental and predicted values are consistent, supporting the conclusion that the launched wave is an inertial Alfvén wave. A similar analysis can be repeated for each mode (k x , k y ) with wave power above the noise floor.
Related analysis propagates each plane wave mode (k x , k y ) inB y1 to the z-location of Elsässer probe 2 using the multiplicative factor e ikz(z 2 −z 1 ) and k z given by Supplementary Equation (1). The result is inverse transformed, and if the wave as a whole is successfully described by a superposition of forward propagating inertial Alfvén plane waves, it should recreate B y2 . Supplementary Fig. 2 shows a snapshot of the results from this procedure. The measured B y1 , when propagated to the location of probe E2, successfully reproduces the amplitude, x-y structure, and temporal evolution of B y2 . Notably, sharper features measured in B y1 , corresponding to higher values of k ⊥ , are diminished in both the measured and modeled B y2 . Consistent with this observation, Supplementary Equation (1) predicts higher k ⊥ modes are more damped. Success modeling B y2 using B y1 and Supplementary Equation (1) supports the conclusion that the entire wave pattern measured by Elsässer probes 1 and 2 is a superposition of forward propagating inertial Alfvén waves.
A critical part of the experiment is knowing E z at the location of g e (v z ) measurements. E z is predicted to be small, E z /E x ≈ 0.002. This means an attempt to use a double probe to measure E z would require an impossible near-perfect alignment of the probe withẑ. Even a 1 • tilt of the probe towardx would project enough E x onto the measurement to ruin the effort. Because a direct measurement of E z is not possible with current measurement capabilities, Ampére's law is used instead to calculate E z from measurements of B x and B y . After Fourier transforming in space and time, theẑ component of Ampére's law giveŝ where c is the speed of light, K zz ≈ −ω 2 pe /ω 2 is the (z, z) entry of the plasma dielectric tensor, and ω pe is the electron plasma frequency. An accurate calculation of E z at the location of g e (v z ) measurements uses nearby values of B x and B y . Fortunately, Elsässer probe 2 is less than a meter away inẑ. Success described earlier in this section in recreating B y2 using B y1 measurements from 2.6ẑ m away demonstrates that Elsässer probe 2 data is sufficiently nearby to determine the wave fields at the location of g e (v z ) measurements.

Supplementary Methods 3: Measuring the Reduced Electron Velocity Distribution Function.
The Whistler Wave Absorption Diagnostic (WWAD) provides high-precision measurements of the suprathermal tails of the parallel electron velocity distribution g e (v z ). It operates by sending a small-amplitude probe wave through a short length of plasma, and the measured damping is used to calculate g e (v z ). This technique, called wave absorption, can be implemented using any wave absorbed by resonant electrons [6,7,8,9,10]. In the overdense plasma of the LAPD (ω pe > |ω ce |), the whistler wave is convenient because it is the only propagating wave mode at frequencies just below |ω ce | [10]. As a result, wave data at these frequencies can be unambiguously interpreted as whistler waves. The whistler wave is absorbed by Doppler-shifted cyclotron resonant electrons. The resonance condition is where v z is the velocity of resonant electrons, k wzr is the real component of the whistler wave vector, and ω w is the whistler wave frequency. During each 10 µs measurement of g e (v z ), ω w is chirped to scan v z . Chirping the whistler frequency from |0.75ω ce | to |0.95ω ce | scans v z from 1.8 v te to 14 v te , or in energy units, between 5 eV and 1000 eV.
The central physical principle of this diagnostic is that damping of the whistler wave is proportional to the number of resonant electrons encountered by the wave. Measured damping of the whistler wave between the transmitting and receiving whistler probes can be used to experimentally determine the reduced electron velocity distribution g e (v z ) using where ω pe is the electron plasma frequency and k wzi is the measured damping. The rightmost expression in parentheses shows explicitly the assumption that whistler damping is exponential by relating it to the transmitted and received wave amplitudes, A T and A R , and the length of plasma L plasma between the probes. Measured time of flight of the whistler wave between the two probes gives values for k wzr that agree with the whistler wave dispersion relation [10,11]. Because every quantity on the right hand side of Supplementary Equation (4) is measured or known, g e (v z ) can be calculated from experimentally known quantities.
The relative abundance of thermal electrons causes whistler waves with frequencies resonant with thermal velocities to be completely absorbed. As a result, the thermal core of g e (v z ) is opaque to this diagnostic technique. However, because the resonant acceleration of electrons by inertial Alfvén waves affects suprathermal electrons, the WWAD measures the region of g e (v z ) where resonant electron acceleration by inertial Alfvén waves takes place. Although ω w reaches frequencies corresponding to 5 eV, whistler signal does not actually propagate to the receiving antenna until frequencies corresponding to 15 eV or higher. Measurements of g e (v z > 0) are available starting at 15 eV, while g e (v z < 0) is measured at energies above 20 eV. The asymmetry in the cutoff of measurements occurs because the suprathermal tails of g e (v z ) are themselves asymmetric. Primary electrons from the cathode cause the g e (v z < 0) tail to be greater than g e (v z > 0), and so g e (v z < 0) becomes opaque at higher energies than g e (v z > 0). Discharge current is emitted between the LAPD's cathode and anode seen in Fig. 1(a), so although g e (v z ) is asymmetric, plasma outside the cathode and anode is currentfree.
WWAD measurements in the LAPD of the suprathermal tails of g e (v z ) are consistent with separate Langmuir probe measurements and kinetic theory. Prior measurements of g e (v z ) were used to find values of n e and T e that agree with simultaneous Langmuir probe measurements [10]. Additionally, WWAD measurements of g e (v z ) resolved in Alfvén wave phase isolated perturbations to g e (v z ) caused by the Alfvén wave that are well-described by linear kinetic theory [12].
Supplementary Methods 4: Applying the Field-Particle Correlation Technique to Experimental Measurements.
The field-particle correlation technique is an innovative analysis method that uses co-located (in configuration space) electromagnetic field and particle velocity distribution function measurements to determine the rate of energy transfer between the electromagnetic fields and particles [13,14,15]. The Maxwell-Boltzmann equations govern the dynamics and evolution of a kinetic plasma. Under the weakly collisional conditions of the auroral magnetosphere, we can drop the collision term in the Boltzmann equation, which is unnecessary to describe the collisionless transfer of energy between fields and the accelerated suprathermal electrons, to obtain the Vlasov equation. Because the linear electron cyclotron frequency f ce = 4760 MHz in the experiment is more than three orders of magnitude greater than the Alfvén wave driving frequency f = 1.177 MHz, we can average the Vlasov equation over the electron gyrophase to obtain the gyroaveraged Vlasov equation. Subsequently, we integrate over the perpendicular component of velocity to obtain the evolution equation for the reduced parallel velocity distribution g e (x, y, z, v z , t) [11], Note the Supplementary Equation (5) is formally dependent on the guiding center coordinates (X, Y ) in the plane perpendicular to the strong axial magnetic field, but because the electron gyroradius, ρ e = v te /ω ce = 0.004 cm, is much less than the accuracy of our spatial position measurements, we can approximate the guiding center position as the spatial position in the perpendicular plane, (X, Y ) ≃ (x, y).
To diagnose the energy transfer between fields and electrons in this experiment, we define the reduced parallel phase-space energy density, w e (x, y, z, v z ) ≡ m e v 2 z g e (x, y, z, v z )/2, for the electrons. Note that the standard form of the field-particle correlation technique is formulated using the full electron phase-space energy density, m e v 2 f e /2, defined using the full distribution function f e (v), rather than the reduced distribution function g e (v z ). The velocity-space signature of the net resonant energy transfer caused by E z is the same for both cases, however, because the v 2 x and v 2 y contributions to the field-particle correlation sum to zero when integrated over v z [14]. Therefore, we can apply our intuition from the previous literature on field-particle correlations to the reduced parallel case here, where our experimental measurements do no separately resolve v x or v y . Multiplying the reduced parallel Vlasov equation by m e v 2 z /2 yields an equation for the evolution of the reduced parallel phase-space energy density w e (v z ), dependent on the parallel electric field E z and the reduced parallel distribution function g e (v z ).
To determine the dominant contribution to the rate of change of w e (v z ) = w e (x 0 , y 0 , z 0 , v z ) averaged over the Alfvén wave period T at the single point of measurement (x 0 , y 0 , z 0 ), we can Fourier transform the distribution function in time over the Alfvén wave period. The discrete Fourier transform of the distribution function over T for g e (v z , t j ) sampled at t j = j∆t for j = 0, 1, . . . , N − 1 and ∆t = T /N is given bŷ where ω n = nω, n = −N/2, . . . , N/2, and the Alfvén wave frequency is ω = 2π/T . The time evolution of each Fourier harmonic is given by The first few Fourier harmonics from the measured g e (v z , t) are well ordered, such that g e0 (v z , t j ) ≫ g e1 (v z , t j ) and g e1 (v z , t j ) ≫ g en (v z , t j ) for n ≥ 2.
Expanding the electron distribution function into its Fourier harmonics g e (v z , t) = g e0 (v z , t j )+ g e1 (v z , t j ) + · · · , substituting into Supplementary Equation (6), and keeping only the lowest order nonlinear term, we obtain an equation for the evolution of the reduced parallel phase-space energy density w e (v z ), At the point of measurement (x 0 , y 0 , z 0 ), this equation tells us that the rate of change of parallel phase-space energy density dominantly depends on the parallel electric field E z along with the constant g e0 (v z ) and fundamental mode g e1 (v z , t) contributions to the Fourier expansion.
Since E z (t) and g e1 (v z , t) both oscillate at the Alfvén wave frequency ω, when this equation is averaged over the Alfvén wave period T , the first two terms on the right-hand side of Supplementary Equation (9) contribute nothing to the net electron energization rate. The third term, a product of two quantities oscillating at the same frequency ω, results in oscillations at frequency 2ω as well as a constant term. This constant term, the amplitude of which depends on the relative phases of E z (t) and g e1 (v z , t), is the dominant contribution to the time-averaged rate of energy transfer between electrons and the parallel electric field E z of the Alfvén wave. Therefore, to isolate the dominant contribution to the net energy transfer rate between the Alfvén wave and the electrons, the correlation is computed, consisting of an unnormalized time average over the full Alfvén wave period. Note that the integration of the reduced parallel correlation over parallel velocity yields the rate of work done by the parallel electric field on the electrons, j ez E z = dv z C Ez (v z ).
The resulting correlation C Ez (v z ) yields the velocity-space signature of energy transfer to the electrons as a function of the parallel velocity v z . The correlation is positive in regions of velocity space where the parallel phase-space energy density increases due to work on the electrons by the parallel electric field, and is negative where the parallel phase-space energy density decreases. Previous numerical [15] and observational [16] applications of the field-particle correlation technique to weakly collisional heliospheric plasma turbulence has demonstrated that energization of particles by the Landau resonance leads to a characteristic bi-polar velocityspace signature, with a decrease in w e (v z ) at velocities just below the resonant parallel phase velocity, and an increase just above the resonant velocity. This is the characteristic signature of electron energization that we are seeking in our laboratory experiments of auroral electron acceleration by Alfvén waves.

Supplementary Methods 5: Gyrokinetic Simulations of Electron Acceleration by Inertial Alfvén Waves.
The gyrokinetic simulation code AstroGK [17] is used to evolve the electron velocity distribution function self-consistently as an inertial Alfvén wave propagates through the periodic simulation domain. Plasma parameters and wave properties are chosen to match the LAPD experiments, with β i = 1.67 × 10 −5 and T i /T e = 0.25 and perpendicular wavenumber of the wave k ⊥ ρ i = 0.05 with k z /k ⊥ ≪ 1. Linear gyrokinetic theory predicts a normalized complex wave frequency of ω GK ≡ ω/k z v A = (0.9708, −7.195 × 10 −4 ), leading to a very weak normalized damping rate of γ/ω = 7.41 × 10 −4 . We specify collisional coefficients for the ions and electrons in the collision operator [18,19] in terms of the wave frequency, as ν i /ω = ν e /ω = 4.2 × 10 −3 , yielding a regime of weak collisionality. The numerical resolution of this 3D-2V gyrokinetic simulation is (n x , n y , n z , n λ , n ε , n s ) = (10, 10, 256, 64, 32, 2), where velocity space coordinates are λ = v 2 ⊥ /v 2 and ε = v 2 /2, uniform Maxwellian equilibria are chosen for ions and electrons, and a realistic mass ratio is specified, m i /m e = 1836. The domain is a periodic box of size L 2 ⊥ × L z , elongated along the straight, uniform mean magnetic field B 0 = B 0ẑ , where all quantities may be rescaled to any parallel dimension satisfying L z /L ⊥ ≫ 1, and with L ⊥ = 40πρ i .
A plane, inertial Alfvén wave with wavevector (k x /k ⊥0 , k y /k ⊥0 , k z /k z0 ) = (1, 0, 8) is initialized with the exact kinetic eigenfunction from linear gyrokinetic theory [20]. To eliminate any transient behavior arising from the initialization that does not agree with the properties of the inertial Alfvén wave mode, we evolved the simulation linearly for five periods with enhanced collision frequencies ν i = ν e = 4.2×10 −2 k z v A , as successfully used in previous studies using eigenfunction initialization of wave modes [21,22].
After the elimination of the transient, the simulation is restarted nonlinearly with the lower collisionality ν s /ω = 4.2 × 10 −3 . At a single point in the simulation domain, the perturbed complementary gyrokinetic electron distribution function g e [17] and the parallel electric field E z are output at uniform time intervals. The parallel field-particle correlation on the gyrotropic velocity space C Ez (v z , v ⊥ ) is computed over the three wave periods 10 ≤ t/T ≤ 13 in the simulation, and the result is plotted in Figure 3(b). Integration over the perpendicular velocity coordinate yields the reduced parallel field-particle correlation C Ez (v z ), a quantity that can be compared directly to the results of the electron acceleration experiments on the LAPD.
Note that, because of the periodic domain, the gyrokinetic simulation results are relevant to measurements made at an asymptotically large distance z/λ ≫ 1 from source of the waves. This condition does not hold true for the experimental measurements, which were measured at a distance z/λ ≃ 2. Furthermore, gyrokinetic theory is valid in the limit ω/ω ci ≪ 1, so the decrease in the real frequency due to the effects of the cyclotron resonance as ω/ω ci → 1, as clear from Supplementary Equation (1), are not included in the gyrokinetic simulation. The resonant parallel velocity in the gyrokinetic simulation is given by ω GK /k z v te = 2.71, as evident in Fig. 3(b) and (c). Using the linear Vlasov-Maxwell dispersion relation to account for the finite frequency value ω/ω ci = 0.45 from the experiment, a complex frequency value of ω VM = (0.863, −1.50 × 10 −3 ) is obtained. The finite frequency effects lead to a reduction in the parallel phase velocity, and thereby a reduction in the resonant parallel electron velocity to ω VM /k z v te = 2.4. The modeling of the electron velocity distribution using Liouville mapping, as described in Supplementary Methods 7, uses this more physically accurate value of ω VM in prescribing the inertial Alfvén wave.

Supplementary Methods 6: Laplace Transform Solution of Field-Particle Correlation at Finite Distance.
Assume an inertial Alfvén wave exists with an electric field parallel to B 0 with the form E z = E z0 e ik·x−iωt , where k is complex, ω is real, and both satisfy the inertial Alfvén wave dispersion relation. E z is the only wave field that explicitly enters the reduced gyroaveraged equation for the evolution of g e in Supplementary Equation (5). The distribution g e (x, y, z, v z , t) is separated into a uniform static background g e0 (v z ) and a small perturbation caused by the inertial Alfvén wave the varies in space and time g e1 (x, y, z, v z , t). The linearized form of Supplementary Equation (5) is To account for the finite interaction length of electrons with the Alfvén waves in this experiment, Supplementary Equation (11) is solved using a Laplace transform over z. Because the duration of the inertial Alfvén wave burst is much longer than the lifetime of electrons traversing the LAPD, a Fourier transform is applied over time. The result of applying these transforms to Supplementary Equation (11) is where a tilde is used to denote the Laplace and Fourier transformed versions of the quantities g e1 and E z , and p is the Laplace transform conjugate variable of z. The boundary condition at the antenna g e1 (x, y, z = 0, v z , t) is introduced in Supplementary Equation (12) by the Laplace transform and acted on by the Fourier transform.
Supplementary Equation (12) can be solved algebraically forg e1 , and the real space solution g e1 is produced by carrying out inverse Laplace and Fourier transforms. The result is The first line on the right hand side of Supplementary Equation (13) is the boundary condition propagated ballistically to (x, y, z). The second line is the linear wave-particle interaction. As expected, the numerator of the wave-particle interaction is zero when the interaction length z = 0. The resonant denominator allows the wave particle interaction to grow for electrons at or near the inertial Alfvén wave phase velocity v z = ω/k z .
The solution for g e1 discussed so far only applies to electrons originating at the Alfvén wave antenna and traveling from there to the location of g e measurements. In other words, the solution described above is valid for g e1 (v z > 0). Following a similar procedure, a solution is produced for g e1 (v z < 0). (13), a boundary condition g e1 (x, y, z = 0, v z , t) is needed. While measurements of the perturbation to the distribution function are not performed at the antenna, a reasonable approximation can be generated by requiring that the boundary condition be consistent with the fluid properties of the inertial Alfvén wave. Specifically, g e1 (x, y, z = 0, v z , t) must include the parallel current of the Alfvén wave required by two fluid theory J z = in e e 2 E z /(m e ω) [23]. To this end, the boundary condition is given the form

To evaluate Supplementary Equation
The first moment of this term reproduces the required J z of the inertial Alfvén wave.
To this point, collisions have been omitted from the discussion of g e1 . Collisions are introduced using a velocity-dependent Krook collision operator [2,11]. This operator adds a collision term to the right hand side of Supplementary Equation (11) with the form −ν(v z )g e1 where ν(v z ) is the velocity-dependent Coulomb collision rate. The result of the collision term is ω → ω + iν(v z ) in Supplementary Equation (13) and Supplementary Equation (14), except in the Fourier kernals e iωt and e −iωt and the implicit ω in E z . By adding collisions, this generalizes our solution to the collisionless Vlasov equation to be a solution to the weakly collisional Boltzmann equation.
Alfvén wave antennas are often inefficient because Alfvén wavelengths are usually much longer than the antennas that can be installed through typical vacuum ports. This experiment is no exception. As a result, some of the power injected into the plasma by the antenna is likely not deposited into the inertial Alfvén wave mode, and some portion of the measured g e1 will be due to non-Alfvénic oscillations. In earlier work [11], non-Alfvénic perturbations were modeled using the velocity perturbation of electrons traversing a voltage drop across the antenna. This effect is included as a homogeneous solution to Supplementary Equation (11) in the particular solution for g e1 . Because the amplitude and phase of the voltage drop across the antenna are not measured, they are allowed to be adjusted to minimize the χ 2 difference between the measured and modeled g e1 . The amplitude and phase of the voltage drop across the Alfvén wave antenna are the only free parameters in the entire model for g e1 . The voltage drop across the antenna is typically 1 volt or less. Of all the contributions to the modeled g e1 , the largest is the waveparticle interaction in Supplementary Equation (13). The wave-particle interaction requires a finite interaction length to develop, which ultimately explains why the resonant signature in C Ez evolves near the antenna and the bipolar signature of resonant electron acceleration asymptotically approaches the resonant velocity as distance from the antenna increases.
The solution for g e1 in Supplementary Equation (13) is for a single inertial Alfvén plane wave mode. As discussed in Supplementary Methods 2, the inertial Alfvén waves generated in this experiment are a superposition of multiple plane waves. Because of this, a realistic model requires all the g e1 's, each attributable to a separate inertial Alfvén plane wave mode, to be summed. The result is used for comparing with the measured g e1 as well as for predictions of C Ez in this experiment.

Supplementary Methods 7: Liouville Mapping of the Electron Velocity Distribution.
We can model the effect of the inertial Alfvén wave on the electron velocity distribution by a mapping of phase space through the wave fields according to Liouville's theorem, a technique we denote here as Liouville mapping [24,25]. In this technique, we build up each point in the electron velocity distribution at the spatial point z obs > 0 and time t obs > 0 of observation by mapping backwards in time according to the Vlasov equation (where the trajectory of a small volume in 3D-3V phase space at velocity v i is determined by its velocity and the Lorentz force law). This backwards time integration is continued until the small volume of phase space reaches the time t = 0 at which the wave was initially launched at z = 0 from the antenna, or until the particle reaches the antenna at z = 0. At the end of this time integration, we obtain the velocity of that point in phase space, v f . In the absence of collisions, the relevant case for the Vlasov equation, we can invoke Liouville's theorem to set the phase space density at v i at the observation time and position equal to the phase space density at v f of the "initial distribution" at the end of the backwards time integration.
Therefore, one simply needs to determine the form of the "initial distribution" from physical considerations. We take our initial conditions for the electrons at t = 0 to be an isotropic Maxwellian distribution with T e = 4 eV everywhere. We consider that the antenna turns on at t = 0, launching an inertial Alfvén wave that propagates from z = 0 in the positive z direction with the appropriate parallel phase velocity given by the Vlasov-Maxwell dispersion relation, ω VM = ω/k z v A = 0.863. See discussion at the end of Supplementary Methods 5 about the difference between the complex linear inertial Alfvén wave frequencies using the linear gyrokinetic dispersion relation [20] or the linear Vlasov-Maxwell dispersion relation [26]. For the Liouville mapping results here, we choose on the Vlasov-Maxwell approach. For each volume of phase space (at velocity v i at the start of the backwards integration) that remains at z > 0 for all t ≥ 0 as it is mapped back in time, we simply need to integrate back in time to t = 0. In this case, the phase-space density at the position and time of observation at velocity v i is simply set equal to the phase-space density of the initial Maxwellian equilibrium at velocity On the other hand, for volumes of phase space (at velocity v i at the start of the backwards integration) that reach the antenna at z ≤ 0 before reaching t = 0 when mapped back in time, we require a more sophisticated treatment of the "initial distribution". We take the antenna to be a physical barrier at z < 0, so that we can only sample the electron distribution at z ≥ 0. Therefore, we integrate the electron phase space trajectories backwards in time only until that phase-space volume reaches z = 0. We adopt an idealized approximation that the wave appears at z = 0 and can only effect the electron distribution after integration over a finite distance along z > 0. Following the full solution of the linearized kinetic equation for an inertial Alfvén wave generated by an antenna outlined in Schroeder et. al (2017) [11], we account for two specific "near-antenna" effects that can influence the electron velocity distribution measured at the probe: (i) a "fluid" component to the boundary condition that arises from the varying electrostatic potential φ(t) of the Alfvén wave driven by the antenna; and (ii) a "homogeneous" term in the solution that yields a boundary condition describing the non-Alfvénic response of the plasma in the sheath region at the antenna surface.
The fluid contribution to the z = 0 boundary condition yields a standard Boltzmann correction to the Maxwellian distribution due to the electrostatic potential φ(t) of the Alfvén wave at the time t z=0 when the phase-space volume reaches the antenna. The electrostatic potential at the antenna is determined by integrating from beyond the leading edge of the wave at z * back to z = 0, φ(z = 0, t z=0 ) = − 0 z * E z (z, t z=0 )dz. Because the wave launched from the antenna at z = 0 has no length of interaction over which it can modify the electron velocity distribution, the only effect on the distribution at z = 0 is the varying electrostatic potential needed to launch the wave. Therefore, in this case, the phase-space density at the position and time of observation at velocity v i is simply set equal to the phase-space density of the Maxwellian equilibrium at velocity v f modified by the electrostatic potential, The homogeneous boundary conditions adds a perturbation δf h to f obs (v i , t obs ) due to the varying potential change through the plasma sheath given by ∂v z e iω(z obs /vz−t obs )+iαs (16) where V s0 is the amplitude of the plasma sheath potential and α s is the phase of the sheath potential relative to that of the Alfvén wave [11]. The exponent gives the phase of the sheath potential ballistically propagated to the position z obs and time t obs at which the electron velocity distribution is observed. The model shown in Fig. 3(f)-(h) uses V s0 = 0.28 V and α s = −π/2.
Because the resonant acceleration of the electrons is governed by the parallel component of the electric field E z , we do not model the effect of the perpendicular electric field on the electron distribution (which would lead to a conservative oscillation of the distribution in the perpendicular direction that averages out over one period). Rather, we model only the parallel electric field at z ≥ 0 by The Liouville mapping code was validated by comparison to published results in Kletzing (1994) [24]. The dimensionless parameters of the Liouville mapping calculation are β i = 1.67 × 10 −5 , m i /m e = 1836, v te /v A = 0.35, perpendicular wavenumber k x δ e = 0.29, parallel wavenumber k z δ e = 9.63 × 10 −3 , wave frequency ω/k z v A = 0.863, and Alfvén wave amplitude E ′ z = cE z /v A B 0 = 2.4 × 10 −7 . In Fig. 3(f)-(h), the electron velocity distribution function f e and E z needed to compute field-particle correlation C Ez (v ⊥ , v z ) are determined at z obs /δ e = 1386, or z obs /λ z = 2.12, with the correlation taken over 100 uniformly spaced samples in time over the time range 2.13 ≤ t/T ≤ 3.13, where T is the period of the Alfvén wave. The electron velocity distributions plotted in Fig. 4 are taken (a) at z obs /δ e = 1386 (equivalently z obs /λ z = 2.12 or z obs = 5.27 m), and t obs /T = 3.13 and (b) at z obs /δ e = 6930 (equivalently z obs /λ z = 10.61 or z obs = 26.35 m) and t obs /T = 11.63. The nonlinear trapping width [27] in parallel velocity, shown by the dashed vertical gray lines in Fig. 4, is given by ∆v z /v te = ± eφ/T e = ±0. 19, showing good agreement with the extent in v z of the effect of the wave on the electron velocity distribution.

Supplementary Methods 8: Comments on Connection Between the Experiment and the Three Modeling Approaches.
A few brief notes on the connection between the LAPD lab experiment and the AstroGK simulation, Laplace transform solution using analytical theory, and Liouville mapping models are worthwhile mentioning here. First, the gyrokinetic simulation self-consistently determines the electromagnetic fields from the electron (and ion) velocity distributions, but it does not include the finite frequency effects (due to finite value of k z ρ i , which results in a finite value of ω/ω ci ), which reduce the wave frequency, in this case to 89% of the gyrokinetic wave frequency (see discussion at the end of Supplementary Methods 5). Furthermore, the gyrokinetic simulation results are also effectively in the limit of an asymptotically large distance from the antenna, so they do not account for a finite distance from the antenna. Second, the Laplace transform solution directly takes into account the finite distance from the antenna, showing that the zero crossing of C Ez (v z ) shifts away from the resonant velocity towards zero. Finally, the Liouville mapping result can be used, assuming an inertial Alfvén wave frequency that includes the finite-frequency effects (and a similarly altered value of k z for the given experimental value of the driving wave frequency) to make a direct comparison to the experimental results. Although the dimensional amplitude of the reduced parallel correlation C Ez (v z ) from the Liouville mapping result is about a factor of two smaller than the experimental determination, this is likely due to the fact that the Liouville modeling used a single-temperature Maxwellian equilibrium velocity distribution and only included the dominant Alfvén wave mode. The experimental measurements found a more complicated equilibrium velocity distribution and the analysis included additional subdominant wave modes generated by the antenna, so this likely is the root of the overall amplitude difference. But, the most important aspect of the comparison between the Liouville mapping results and the experiment is that a simple Maxwellian equilibrium model quantitatively recovers the velocity of the zero crossing and peak of C Ez (v z ), the key features that confirm the characteristic velocity-space signature of Landau-resonant electron energization. An additional benefit of the Liouville mapping modeling is that we can trivially extend the calculation to determine the perturbed electron distribution function, specifically the generation of a population of accelerated electrons, at distances beyond the physical length of the experimental chamber, as shown in Fig. 4.  Table 2: Comparison of plasma parameters. Parameters measured or estimated in the LAPD experiment compared to those predicted in the auroral acceleration region, where the key dimensionless parameter characterizing this region is the ratio of the electron thermal to Alfvén velocity, v te /v A (bold).
For the comparison of the experiment to the corresponding position in the auroral acceleration region with v A /v te = 3 at s ≃ 2.85 R E , or an altitude z ≃ 2.56 R E , the resulting plasma parameters are given in Supplementary Table 2. At this position, we take a typical Alfvén wave amplitude of δB ⊥ = 10 nT [32,33,34], which corresponds to a perpendicular electric field wave amplitude of E ⊥ ∼ 200 mV m −1 . For a perpendicular wavenumber k ⊥ δ e ∼ 1, the perpendicular wavelength is approximately λ ⊥ ∼ 10 km. The parallel wavelength is typically estimated to be a factor of 10 to 100 times longer than the perpendicular wavelength [31], and such a value is also consistent with the parallel wavelength that would be estimated by the condition of critical balance in strong Alfvénic turbulence [35]. Here we take the parallel wavelength to have a value of L = 500 km. With the value of L /L ⊥ ∼ 50 and k ⊥ δ e ∼ 1, the relation between for E /E ⊥ given by Chaston (2006) [31] yields a value E z ∼ 2 mV m −1 .