Plasma oscillations below 100 kHz in a high-power, magnetically shielded Hall-effect thruster are characterized experimentally. A time-resolved laser-induced fluorescence diagnostic was used to measure the evolution of the axial ion velocity distribution along the discharge channel centerline. A method was developed to correct for artificial broadening of the distribution due to finite spatial resolution, enabling accurate ion temperature measurements in the acceleration region. Time-dependent ion heating behavior is revealed, which varies with axial location. Electrostatic, collisionless particle-pushing calculations were employed to simulate the effects of a shifting acceleration front on ion motion. It is found that ions exchange energy with these oscillations, which accounts for a portion of the ion velocity fluctuations observed in the thruster plume. Phasing relationships between ion dynamics and discharge current oscillations are discussed in the context of prior Hall thruster experiments.

## I. INTRODUCTION

Hall-effect thrusters generate thrust by accelerating ions through an electrostatic field in an annular, crossed-field plasma discharge.^{1} These devices have been in use for several decades, but until recently, deep-space applications were limited by rapid erosion of the discharge channel walls. This obstacle has been overcome by magnetic shielding (MS), which was discovered through the observation of a zero-erosion state in Aerojet’s BPT-4000 life test^{2} and subsequently understood and developed through simulations and experiments at the Jet Propulsion Laboratory (JPL).^{3,4} Several unanswered questions remain regarding the fundamental operation of Hall thrusters, especially at these new operating regimes. These include the causes of magnetic pole erosion in MS thrusters,^{5–7} in addition to anomalously high electron mobility^{8} and the interactions between these effects and various oscillations and instabilities.^{9} Such phenomena can be investigated with time-resolved laser-induced fluorescence (LIF), an experimental technique that provides non-invasive measurements of the local evolution of the ion velocity distribution function (IVDF).

Hall thrusters commonly exhibit a variety of oscillations and instabilities. Among the most prominent are large-amplitude, low-frequency ionization oscillations known as the “breathing mode.”^{10} Additionally, azimuthal modes and ion acoustic turbulence exist in the hollow cathode plume,^{11–13} along with high-frequency electron–cyclotron drift waves in the discharge channel.^{14,15} This array of oscillations influences the distribution of particle energies in non-classical ways, which leads to enhanced cross-field electron transport^{16} and anomalous ion heating.^{17–19} High-power MS Hall thrusters demonstrate a unique collection of oscillatory behavior, the nature of which can be highly dependent on operating conditions, such as discharge voltage and magnetic field strength.^{20} Since these effects can lead to the production of high-energy ions and increased erosion of the pole covers (e.g., see Refs. 6, 17, 21, and 22), it is essential to characterize the plasma oscillations for the faithful simulation of thruster operation and erosion.

This work presents results from time-resolved LIF measurements taken on the Hall-effect rocket with magnetic shielding (HERMeS),^{23} a 12.5 kW magnetically shielded Hall thruster being developed in concert by NASA and Aerojet Rocketdyne for the Advanced Electric Propulsion System (AEPS).^{24} In this study, HERMeS was operated at discharge voltages of 300 and 600 V, both with discharge current regulated at 20.8 A. In the 300 V configuration, the discharge current signal is characterized by aperiodic fluctuations with a wide-band frequency spectrum between 10 and 40 kHz, whereas at the high-power, 600 V configuration, the thruster exhibits large, nearly periodic current oscillations at 55 kHz. Figure 1 shows Fourier spectra of discharge current traces from both conditions.

Time-resolved LIF measurements have demonstrated distinctly different characteristics of IVDF behavior in the presence of 300 and 600 V oscillation modes.^{25} The 300 V mode witnesses aperiodic oscillations primarily in the breadth of the distribution, associated with a higher effective ion temperature at certain phases. In contrast, the 600 V periodic mode is associated primarily with motion of the mean velocity of the ion fluid, with the acceleration region appearing to shift axially throughout the oscillation cycle. The 600 V configuration also displays large variations in ion beam energy downstream of the acceleration region, in the near-field plume. At the peaks of these oscillations, the ions acquire kinetic energies exceeding the discharge voltage, suggesting that effects other than electrostatic acceleration may be relevant to ion motion.

Plasma simulations of HERMeS operation on a magnetic field-aligned mesh have been carried out with Hall2De. This code employs a multi-fluid approach in the discharge channel and transitions to a particle-in-cell treatment of ions when they become collisionless within the acceleration region.^{17,26} These computations are validated with experimental wear tests^{27} and used in the life qualification process for long-duration science missions. While the elusive character of the oscillations in question often prohibits their simulation from first principles, Hall2De studies have demonstrated that in the high-voltage configuration, simulating the breathing mode by forcing oscillations in the anomalous transport profile with amplitudes inferred from LIF data^{28} results in a realistic discharge current oscillation amplitude and realistic erosion rates for the front pole covers.^{17}

In contrast to the success of erosion modeling at 600 V, classical models have consistently underpredicted erosion rates at the 300 V condition. Recent work, however, has demonstrated evidence for the growth of lower hybrid instabilities in the inner front pole region with the potential to heat ions anomalously.^{18,19} During preliminary simulations with Hall2De in which ion heating was enhanced by simple scaling with the lower hybrid frequency, significantly improved agreement with the measured erosion rates was achieved.^{17} These findings further implicate wave-based heating mechanisms in the production of high-energy ions. Time-resolved, non-invasive measurements are necessary to experimentally explore the existence of these instabilities.

We present further in-depth analysis of time-resolved LIF data from HERMeS in order to more fully characterize the interactions between discharge oscillations and ion kinetics throughout the thruster. We develop a method to correct for spatial averaging of LIF measurements in the acceleration region to obtain accurate ion temperatures. This procedure enables investigation of the characteristics of ion heating in the low-power configuration. Additionally, ion energy oscillations are jointly investigated by pairing the LIF results with quasi-electrostatic particle-pushing simulations. These calculations demonstrate oscillation-particle energy exchange in the acceleration region. These approaches are used to study phasing relationships between discharge current and fluid properties, which are relevant to ion heating and anomalous transport effects.

This article is organized as follows: Section II describes the experimental methods and analysis approach in detail, including techniques for the deconvolution of spatially averaged velocity distributions and the simulation of ion acceleration in the presence of low-frequency oscillations. Analysis of IVDF oscillations for the 300 and 600 V operating conditions is presented in Secs. III A and III B, respectively. Section III C then discusses the ion simulation results, followed by a summary of notable results in Sec. IV.

## II. METHODS

### A. Time-resolved laser-induced fluorescence

Laser-induced fluorescence (LIF) is a non-invasive plasma diagnostic useful for providing kinetic descriptions of ion motion in plasma sources, such as Hall thrusters. LIF uses a laser to excite metastable electronic states in plasma ions, which decay and fluoresce. Due to the Doppler shift in the reference frames of moving ions, tuning the laser through a range of wavelengths around the targeted transition and monitoring the fluorescent light intensity allows for a mapping from wavelength to ion velocity, where the intensity as a function of velocity approximates the IVDF up to a normalizing factor.^{29}

A ubiquitous drawback of optical plasma diagnostics is poor signal-to-noise ratio (SNR) due to background light and unwanted laser scattering. This is mitigated in typical LIF setups by modulating the laser amplitude at a known frequency with mechanical chopping or an acousto-optic modulator (AOM). The fluorescence signal is then separated from the background collected light with phase-sensitive detection (PSD) via a lock-in amplifier. The characteristic integration time associated with the PSD low-pass filter limits the time resolution of the LIF technique but averages over background noise.

Higher time resolution can be achieved by choosing a faster PSD filter and averaging measurements over a repeatable process to increase relative signal intensity. Techniques, such as triggered ensemble averaging, may be used to study time-variation that is periodic or quasiperiodic.^{30} However, Hall thrusters can exhibit aperiodic oscillations for which these methods are unable to lock onto the phase, for example, those found in HERMeS at 300 V (see Ref. 25). Improved SNR for time-resolved LIF diagnostics in such an aperiodic context has been achieved through the transfer-function averaging technique developed by Lobbia and Durot.^{31}

This approach improves SNR by treating the thruster as a linear time-invariant (LTI) system and measuring average transfer functions to estimate the LIF signal from a reference signal, provided that the Fourier spectrum is relatively constant. It is assumed that for a given laser wavelength and position, the Fourier transform $F(\omega )$ of the true LIF signal $f(t)$ is related to the reference signal $I(\omega )$ (typically the discharge current in this context) by a complex transfer function $H(\omega )$ such that

where $\omega $ is the angular frequency. $H(\omega )$ is assumed to be constant in time but is a nontrivial function of frequency. By collecting the LIF signal for a long time period ($\u223c$30 s) and breaking the measurement interval into a large number of subintervals, the approximate transfer function is calculated in the Fourier domain for each laser wavelength by averaging over these subintervals. The signal at each wavelength is then deduced from the same reference signal via Eq. (1). This method improves the time resolution dramatically and has been shown to agree well with other, more limited methods for time-resolved LIF, such as phase-sensitive boxcar averaging.^{31}

This work expands on the analysis of Ref. 25, in which the experimental setup below is described in more detail. We carried out time-resolved LIF measurements using the transfer function averaging approach on the HERMeS Technology Demonstration Unit 2 (TDU-2) in JPL’s 3-m diameter by 10-m Owens vacuum facility, with background pressure between 10 and 14 $\mu $Torr. The LIF diode laser was tuned near the Xe II transition at 834.953 nm in vacuum and was modulated at frequency $fAOM\u22481.8$ MHz. Lock-in amplification was carried out in post-processing, with time constants $\tau PSD=1.0\mu $s for the 300 V oscillations and $\tau PSD=0.7\mu $s for the 600 V oscillations. This configuration allows minimum attenuation of signals below 100 kHz. The coordinate system used in this paper denotes the axial dimension as the $z$-axis, with the origin in the plane of the anode. All distances are normalized to the discharge channel length $Lch$. Five points along the channel centerline were probed in both voltage configurations, at $z\u2212$positions $[0.86,1.11,1.16,1.21,1.50]$ for 300 V and $[0.82,0.96,1.01,1.06,1.50]$ for 600 V.

### B. Data analysis

Post-processing with the average transfer-function approach yields collected fluorescence signals as a function of laser wavelength (corresponding to axial ion velocity), time, and position. Relevant plasma parameters, such as mean velocity and ion temperature, can be calculated by taking successive velocity moments of the distribution. We fit these discrete distributions with a bi-Maxwellian distribution function of the form

The use of a second Gaussian curve in this model allocates additional degrees of freedom for fitting small deviations from equilibrium, for example, the presence of a low-energy tail to the distribution in Fig. 2, and does not necessarily represent a separate ion population. A nonlinear least-squares algorithm is used for this fitting, with bounds on parameters specified manually by inspecting the distributions. Curve-fitting serves the dual purpose of smoothing noise as well as providing a functional model, which can be treated analytically.

The mean velocity and effective temperature can be straightforwardly extracted from the curve fit parameters. For a single-Maxwellian distribution in the z-direction [Eq. (2) with $A2=0$], the mean velocity and effective temperature are $v\xafz=v1$ and $Ti=mivT12/2$. These results can be extended to the bi-Maxwellian case by analytically taking the first and second moments of Eq. (2). Analytical expressions for these moments are found in the Appendix, and examples of fits to time-resolved IVDF snapshots from the data are shown in Fig. 2. In this context, the *effective* temperature $Ti$ quantifies the variance of the ion distribution in units of energy, representing an extension of the equilibrium thermodynamic property.

It is important to note that $Ti$ is a 1D quantity specifying the ion energy spread along the interrogation axis. This does not necessarily correspond to an isotropic, equilibrium temperature, and the word “heating” in this paper refers to an increase in this energy spread, without necessarily implying thermalization. Time-averaged^{32} and time-resolved^{25} LIF measurements in HERMeS along two orthogonal axes rotated 45$\xb0$ from axial and radial have produced evidence of anisotropic $Ti$ both on the channel centerline and in the beam edges. Radial IVDF measurements would be of significant interest to investigate ion heating in the direction perpendicular to the electrostatic acceleration, but these data have not been acquired due to lack of line of sight access in the channel.

Figure 2 also demonstrates an example of an artifact from the data-taking process filtered by the curve-fitting algorithm: for some wavelengths and times, the signal value appears negative far into the wings of the distribution. The negative values are likely caused by high-frequency ringing from the transfer-function conversion in the Fourier domain. Numerical methods exist for dealing with these spurious negative values before implementing quadrature,^{31} but analytical curve fits constitute a more computationally simple approach for obtaining bulk plasma properties from noisy distributions. It must be noted, however, that removing these negative values may slightly perturb the effective temperatures measured by the curve fits from their true values.

### C. Deconvolution of artificial broadening

Both atomic physics effects and finite spatiotemporal resolution can artificially broaden the apparent velocity distribution. Such phenomena are typically symmetric; therefore, the mean velocity is left unchanged; however, these effects can lead to overestimates of the effective ion temperature. Perturbations may be partially mitigated through careful experimental setup, but further analysis is often required to obtain accurate ion properties.

Atomic broadening effects include Zeeman splitting and hyperfine structure. While these splitting effects can be significant in general, both contribute only minutely to this analysis. We minimized the effect of Zeeman splitting by maintaining $\pi $-polarized laser injection relative to the magnetic field.^{7,33} Likewise, we preserved a sufficiently low laser injection power so as to reduce hyperfine structure and laser saturation effects. According to a simple model,^{32} these perturbations contributed less than $10%$ error for the observed ion temperatures. For these reasons, we neglect such atomic broadening and focus on spatial lineshape averaging.

If the true distribution varies significantly with position or time, such as in the presence of oscillations or large spatial gradients, the collected distribution is effectively a convolution of the true IVDF with the spatiotemporal collection window. This convolution has the effect of smoothing and broadening the collected distribution. Operating with increased time resolution, as was done in this study, removes the effects of temporal averaging over the low-frequency oscillations of interest. However, spatial averaging over the interrogation volume still non-negligibly contributes to the apparent width of the distribution. This distortion is especially large in the acceleration region, where the mean velocity varies strongly with position. Based on time-averaged IVDFs collected in an earlier study on HERMeS in Ref. 34, the mean ion velocity within the interrogation region can vary by up to 5 and 3 km/s for $Vd=600$ and 300 V, respectively.

With a few assumptions, it is possible to analytically deconvolve the effect of spatial averaging from the collected IVDF to correct for artificial broadening. We model spatial averaging as blurring of the IVDF by a Gaussian kernel, a common approach in image processing. The circular cross section of the optical collection beam is thus assumed to have a Gaussian intensity profile, which reduces the blurring to a straightforward convolution integral. The key parameter is the kernel width $\sigma $, which was estimated using an adjustable aperture and an optical power meter to measure the profile diameter of a laser focused by the optical collection apparatus in the reverse direction. As shown schematically in Fig. 3, the value of $\sigma $ used must take into account the projection angle of the collection beam onto the laser injection axis. In the present work, fluorescence was collected at an angle of $\theta =20\xb0$ from normal to the laser axis; therefore, the measured radius of the beam must be multiplied by a factor $sec\theta $ to represent the true averaging scale along the beam axis.

The measurement function at an axial position $z$ is of the form

The blurring process can be simulated by allowing $z\u2032$ to range over all values. We now make the crucial assumption that within the collection volume, the true IVDF $f(z,vz)$ shifts in mean velocity at different positions but maintains a fairly constant lineshape. This allows for the linear approximation of the IVDF across the collection volume with a first-order Taylor series in velocity space, according to the mean velocity gradient $\u2202v\xafz/\u2202z$ at each position. Defining $\tau =(z\u2032\u2212z)\u2202v\xafz/\u2202z$, the measured ion velocity distribution can then be written as the convolution integral

where

Working backward to obtain the true distribution from $fmeas$ can be achieved with numerical deconvolution algorithms, but this difficult process tends to amplify noise and depends on heuristically tuned hyperparameters. For the present case, we exploit the convenient properties of Gaussian functions to yield an analytical result consistent with the bi-Maxwellian fitting model.

Consider a single-Maxwellian distribution with amplitude $A$, mean velocity $v\xafz$, and thermal velocity $vT$. For this distribution, Eq. (4) yields another Maxwellian with unchanged mean velocity but with thermal speed

and amplitude

Inverting these relations yields the true parameters $A$ and $vT$. Due to the linearity of the convolution operator, the corrections can be applied independently to both distributions in a Bi-Maxwellian fit. The deconvolved temperature can then be calculated from the corrected parameter values as described in the Appendix.

The analytical deconvolution approach was used to correct the effective ion temperature measurements in the 300 V operating condition. Velocity gradients were estimated numerically from previous time-averaged LIF measurements in Ref. 34. A corrected temperature trace is shown in Fig. 4 at a position of 1.16 $Lch$, in the center of the acceleration region. The algorithm generally produced corrected temperatures 1–4 eV below the spatially averaged values, as expected. However, the process is dependent on the simplifying assumptions made. Equation (6) broke down in one position (1.11 $Lch$), with the correction factor under the radical exceeding $vT\u2217$. This is likely due to the lineshape changing rapidly with position in this location (Fig. 5), which invalidates the Taylor approximation made here. In order to approximate the temperatures at 1.11 $Lch$, a single-Maxwellian fit was instead used; however, this still led to unreasonably low temperatures sometimes, which are excluded from the results. The deconvolution procedure was not applied to the 600 V data since the highly time-dependent spatial velocity gradient is not well-approximated by time-averaged data in this configuration.

### D. Quasi-electrostatic ion ensemble model

In the 600 V operating condition, phenomena of interest include both the presence of high-energy ions in the plume and interactions between discharge oscillations and the ion velocity distribution. As an aid for interpreting the LIF data, we developed a simple computational model for the simulation of collisionless ion trajectories through the oscillating accelerating potential. This framework ignores particle interactions and integrates ion motion through a time-dependent electric field $Ez(z,t)$, which is imposed (not calculated self-consistently) in these simulations. These simplified dynamics should serve as a suitable approximation downstream of the ionization region, where particles are largely collisionless. This quasi-electrostatic approach provides a computationally inexpensive means for isolating the portion of the ion dynamics, which corresponds to single-particle acceleration in a slowly varying electric field.

We carried out 1D particle-pushing calculations for individual ions, as well as for an ensemble of ions. In the latter case, a continuous ion flux derived from LIF data was supplied at the upstream computational domain boundary. We randomly generated initial particle velocities via inverse transform sampling applied to an experimental IVDF snapshot. Applying the electrostatic force from the oscillating potential drop to an ensemble of noninteracting ions demonstrated the evolution of the IVDF throughout the acceleration region. We ignored bulk volume ionization and implemented a standard second-order leapfrog integrator for numerical stability. After reaching a steady state, particle counts at each velocity were computed using binned averaging over several simulation time steps. Binning particle velocities results in approximate simulated IVDFs with resolution in both position and time. We compare the results of these simulations with time-resolved IVDF measurements in Sec. III C.

Within the particle-pushing simulations, we tried two approaches for implementing the electric field oscillations at the 600 V operating condition. In the first case, conservation of energy was used to approximate the time-averaged electric field from time-averaged LIF measurements,^{34} and then, oscillations were simulated by shifting this electric field profile back and forth axially. We estimated the shape of the time-averaged potential profile via single-particle conservation of energy based on the LIF mean velocity $u(z)$, yielding $\varphi (z)=\varphi 0\u2212miu(z)2/(2e)$, where $\varphi 0$ is the total potential drop. In this case, the (normalized) waveform of the axial shifting was chosen to be a smoothed square wave of form $w(t)=tan\u22121\u2061(cos(\omega t\u2212\delta )/\gamma )$, which served as a simple analytical intermediate between purely sinusoidal and square-wave oscillations to match the upstream ion behavior. This expression produces a pure sinusoid as $\gamma $ grows large, whereas for small $\gamma $, the behavior approaches a square wave. For this work, we chose $\gamma =0.1$ and tuned the phase $\delta $ to match the upstream LIF data.

The second, preferred approach for implementing the electric field oscillations was to use $Ez(z,t)$ extracted from previous Hall2De simulations of HERMeS;^{5} this provided a higher-fidelity estimate that accounted for realistic variations in the potential drop magnitude as the location of the acceleration region oscillated. Readers may refer to the literature (e.g., Refs. 5, 6, 17, and 26) for detailed descriptions of the physics incorporated in the 2D (r–z) hybrid multifluid-PIC Hall2De simulations. The Hall2De model requires an experimentally calibrated anomalous transport coefficient, parameterized as the 2D anomalous collision frequency profile $\nu AN(z,r,t)$, in order to account for non-classical mobility of electrons across the radial magnetic field. Previous simulations found that the global discharge oscillations observed in HERMeS could be accurately replicated by sinusoidally shifting this anomalous transport profile in time.^{6} The amplitude and shape of the oscillations were validated through detailed comparisons between the simulated and measured time-dependent discharge current and ion velocity oscillations (Fig. 18 of Ref. 17). Other aspects of the HERMeS Hall2De simulations have also been extensively validated;^{3,5,6,17,23,25,27,28} therefore, the simulated $Ez(z,t)$ from Hall2De may be used in our simpler 1D particle-pushing simulations as a close approximation of the actual electric field in the thruster.

The weaker oscillations at $Vd=300$ V did not warrant the same forcing in Hall2De. Alternatively, our 1D quasi-electrostatic model of these oscillations simply used the method described above to generate an $Ez(z,t)$ profile derived from time-averaged LIF data.^{34} For the 300 V oscillations, we axially shifted the potential profile from time-averaged LIF with the waveform of a measured discharge current fluctuation.

## III. RESULTS

### A. 300 V IVDF oscillations

Figure 6 displays time averages of the ion temperature plotted over axial distance from the anode along the HERMeS channel centerline. The temperature profile for time-averaged LIF is simulated by first averaging the time-resolved IVDFs over the collection period before extracting temperature information, whereas for the mean time-resolved LIF profile, the temperature is first computed and then averaged over time. The spatial averaging correction technique is used to produce the deconvolved profile, where it must be noted that the second position ($z=1.11Lch$) encountered numerical difficulties as described in Sec. II C. In contrast to the other data points, a Maxwellian fluid approximation was applied to this point, shown as hollow in Fig. 6.

The mean deconvolved $Ti$ profile shows the importance of correcting for spatiotemporal averaging during the LIF collection process for precise measurements of ion temperatures. As expected, the error introduced by finite resolution is the largest in the acceleration region ($z=1.16Lch$), where both oscillation amplitudes and velocity gradients are the largest. The correction at this location is >3 eV, whereas far upstream and downstream of the acceleration zone, there is negligible error due to the relatively flat velocity profile at those positions.

Figure 5 shows instantaneous LIF traces for four different phases of the aperiodic discharge current oscillation in two locations of interest. The distributions exhibit a slight axial shift, but the width of the lineshape varies in a different manner in the two positions despite their adjacency. Upstream, the effective temperature increases due to the presence of a high-energy tail as the discharge current peaks. However, the downstream position observes increased $Ti$ at opposite phases, consisting of a low-velocity tail relative to the mean, which appears during $Id$ troughs.

This phenomenon is more clearly observed by plotting the temperature oscillations directly, as shown in Fig. 7, along with mean velocity and current traces. Such phase plots demonstrate that upstream of the acceleration zone, the ion temperature oscillates in-phase with the discharge current, while at axial locations beyond 1.11 $Lch$, the phasing flips and temperature is out of phase with $Id$. The observed broadening is consistent with the axial motion of the acceleration front such that higher temperatures in both positions occur at phases nearer to the center of the potential hill. However, more complicated ion heating dynamics may be needed to explain the full variation, as these oscillations continue far upstream and downstream of the acceleration zone, where the velocity profile is largely constant.

The oscillations observed in the 300 V operating condition are relatively unique in the context of other time-resolved plasma measurements in Hall thrusters. This may be due to the difficulty of obtaining accurate time-resolved data for aperiodic oscillation modes. One time-resolved LIF study by Fabris *et al*.^{35} demonstrated the appearance of a consistent low-velocity tail at discharge current peaks in a much smaller, 400 W Hall thruster oscillating at 23 kHz. However, their experiment did not witness the opposite oscillations upstream of the acceleration region, as is apparent for HERMeS. Their result is likely attributable to ions born at lower potentials due to a time-dependent ionization region overlap, consistent with the observation of much larger mean velocity oscillations than those found in HERMeS at this configuration. HERMeS ions born at downstream potentials would have velocities less than 1 km/s, appearing as a distinct slow peak in the distribution at 1.16 $Lch$. This is in contrast to the bulk broadening observed in Fig. 7(c). The lack of evidence for separately ionized populations in this larger thruster suggests that the present oscillations may have an altogether different causal relationship with time-dependent ion heating than that found in other Hall thruster literature studies.

### B. 600 V IVDF oscillations

In the 600 V operating condition, the HERMeS acceleration region is located further upstream; therefore, measurements were made at different axial locations to span the potential drop. As noted in Sec. II C, the large mean velocity oscillations in this configuration prevent the spatial averaging correction approach from being used to obtain precise $Ti$ measurements due to the spatial sparsity of the present dataset. However, the lack of temperature offset correction does not alter phasing relationships used to determine the qualitative character of the IVDF oscillations throughout the thruster.

Instantaneous LIF profiles from 600 V operation are shown in Fig. 8. Figure 8(a) displays IVDFs in the center of the time-averaged potential hill ($z=1.01Lch$), which quickly jump from a fast 20 km/s population at discharge current troughs to a slow 3–5 km/s population at times of maximum discharge current. Occasionally, both populations are present in the same frame, but this is most likely due to the diagnostic’s inability to resolve the rapid transitions between the two. Such frames are omitted from calculations of $Ti$. In reality, at the 600 V operating condition, there was never a slow ion population distinct from the main population.

Further downstream in the near-field plume [Fig. 8(b), $z=1.50Lch$], the distributions display a more gradual axial shifting motion. At discharge current troughs, the mean velocity of the ion population contains axially directed kinetic energy in excess of 600 eV, despite having only fallen through a maximum instantaneous potential drop of approximately 560 eV, as predicted by Hall2De. The acceleration of these short-lived, high-energy ions is further investigated computationally in Subsection III C.

Figure 9 compares the phasing of ion fluid properties throughout the thruster. Bimodal distributions due to limited temporal resolution are omitted or treated as separate populations (dashed line) depending on the context. Unlike for the 300 V oscillations, 600 V ion temperatures vary in-phase with the discharge current at all locations. The waveform of the mean velocity oscillation changes dramatically over short distances from nearly sinusoidal upstream to a square wave within the acceleration region. Further downstream, the waveform becomes distorted with sharper peaks and troughs and lags the discharge current as the ions enter the plume.

The upstream oscillations between fast and slow ion populations are likely driven by the approximately binary condition of whether the ions in that location have yet fallen through a strongly shifting potential profile. If the instantaneous potential drop is sufficiently narrow, this would filter the mean velocity curve to yield a square wave in the center of the acceleration region. Such a square wave is observed in the mean velocities at 1.01 $Lch$, as shown in Fig. 9(c). In the downstream positions, however, all of the ions have already been accelerated; therefore, the mean velocity oscillations observed in Fig. 9(e) must be due to another mechanism. A 1D model shows that ion energy exchange with discharge oscillations partially accounts for this behavior, as explored in Sec. III C.

The oscillations in the 600 V configuration appear similar in character to Hall thruster breathing oscillations observed in other contexts, although at an unusually high frequency for a large device, such as HERMeS. Young *et al*.^{30} observed similar behavior in a BHT-600, which oscillated at nearly the same frequency and displayed a comparable downstream phase lag between discharge current and mean velocities. Likewise, similar axial motion of the IVDF is obtained from other time-resolved LIF^{36} and retarding potential analyzer (RPA)^{37} measurements of Hall thrusters at a wide range of designs and power levels.

### C. Ion acceleration simulations

We first apply the electrostatic acceleration model developed in Sec. II D to the trajectory of a single ion in the thruster operating at 600 V. Figure 10 displays plots contrasting such trajectories for constant and oscillating potential profiles. In this figure, the oscillatory results are generated by using the $Ez(z,t)$ profile inferred from sinusoidally shifting the electric field inferred from LIF measurements,^{34} as described in Sec. II D. An immediate observation is that the final ion velocity in the plume is highly dependent on the phase of the oscillation at which the ion crosses through the acceleration region. This effect appears to be due to a “surfing” phenomenon in which the oscillating potential exchanges energy with the ions. When the acceleration region is shifting downstream, it moves with the direction of ions such that the particles spend more time being accelerated; this ultimately leads to a higher kinetic energy. Conversely, at times when the potential is shifting upstream against the ion motion, the ions experience the acceleration for a shorter time, diminishing their final speed. Previous studies have demonstrated the ability of ions to exchange energy in this way with ion transit time oscillations, resulting in energies greater than the discharge voltage.^{38} Our results show that surfing on lower-frequency breathing-like oscillations can also lead to significant ion energy variations.

The quasi-electrostatic model demonstrates that this potential-riding effect can lead to a peak-to-peak kinetic energy amplitude in excess of 50 eV for a sinusoidal oscillation waveform. These simulations reproduce a large portion of the $\u223c$100 eV energy oscillations observed in the plume, without considering additional oscillations in the height and width of the potential profile or collective plasma effects. Low-frequency ion energy exchange, therefore, appears to be a significant contributor to the high ion energies found in the plume, previously attributed solely to oscillations in the cathode-to-plume coupling potential.^{25} Note that these energy exchange effects are included by default in the fluid simulations, but isolating their effects in this context is useful for an interpretable exploration of the direct interaction between potential oscillations and ion kinetics.

The model was next applied to a system of many collisionless ions falling through the time-dependent Hall2De potential profile. Figure 11 shows contour maps of simulated IVDFs compared with analogous plots for the 600 V TRLIF data. Simulations are shown for both the assumption of the Hall2De electric field (middle row) as well as for the axially shifting LIF electric field (bottom row). We synchronized the oscillation phase in the acceleration zone, where both the simulation and data display a rounded square waveform while oscillating between fast and slow populations. In the plume position [Figs. 11(c) and 11(f)], the simulated IVDF waveform with the Hall2DE electric field appears quite similar qualitatively to the measured contour. In both cases, the formation of spikes in ion velocity takes place, leading to a distorted triangle waveform in the downstream position. In the simulation, these spikes are caused by the phase-dependent interaction of the shifting acceleration region with ions, as described above for the single-ion simulation. The potential profile from Hall2De includes higher-fidelity features (height and shape of the instantaneous potential profile), which are more likely to be representative of the actual physics of the thruster than the rounded square-wave potential oscillation. However, even with this different assumption for the electric field oscillation, we observe spikes and troughs in the ion velocity corresponding to the quasi-electrostatic acceleration of ions through a moving potential profile, which are qualitatively similar to the behavior of both the data and the simulation based on Hall2De.

The similarity of Figs. 11(f) and 11(c) suggests that the ion velocity oscillations in the near plume are caused by a combination of the surfing effect and oscillations in the net ion acceleration voltage. The middle position (1.06 $Lch$) witnesses the least agreement between the simulation and LIF data. Phasing matches in this position, but the experimental traces transition more quickly to the distorted waveform found downstream. It appears that due to effects not captured by the electrostatic model or the inferred Hall2De potential, the details of the transition from these upstream, square-wave-like IVDF oscillations to the downstream waveform are not fully described by this simplified model.

We carried out corresponding calculations for the 300 V case in order to explore the extent to which single-particle acceleration accounts for the subtler IVDF motion found at this operating condition. The 300 V fluctuations were emulated by forcing a potential profile from time-averaged LIF measurements to axially shift in-phase with the discharge current waveform. Figure 12 shows the results of these simulations. The mean velocities agree very well between the electrostatic model and the LIF data. This agreement demonstrates that IVDF motion for the aperiodic oscillations is also partially caused by collisionless acceleration in the time-dependent potential. Despite successful modeling of the velocity fluctuations, this quasi-electrostatic model does not account for ion heating effects. More complex mechanisms are responsible for the oscillating breadth of the distribution found at 300 V.

## IV. CONCLUSION

We performed new analysis of time-resolved LIF measurements along the channel centerline of a high-power Hall thruster in order to investigate the relationship between low-frequency plasma oscillations and the evolution of the IVDF. For the aperiodic oscillations present at a discharge voltage of 300 V, we develop an approximate analytical method to correct for artificial IVDF broadening due to spatial averaging over the interrogation volume. Ion temperatures are shown to be non-negligibly inflated by this averaging, and the correction implied by the analytical deconvolution technique suggests artificially broadened IVDFs by up to 4 eV in some positions. Large oscillations in ion temperatures are observed after this correction, reaching an amplitude of close to 8 eV at the acceleration front. Phase coupling between oscillations in discharge current, ion temperature, and ion mean velocity is analyzed in order to explore the ion dynamics in the presence of these low-frequency oscillations. For the 55 kHz oscillations found in the 600 V operating condition, a similar approach is carried out to analyze phasing relationships between discharge current and ion fluid properties. A simplified quasi-electrostatic model is constructed for the simulation of collisionless ion trajectories through this configuration’s strongly oscillating acceleration region, which enables the investigation of ion potential-riding effects that determine the downstream evolution of the IVDFs. This model demonstrates that ion–oscillation interactions account for a significant portion of large plume kinetic energies, previously attributed solely to cathode potential oscillations.^{25}

While both configurations display complex ion behavior, the 55 kHz oscillations at 600 V are comparable with other analyses of Hall thruster IVDFs with breathing mode oscillations. The temperature oscillations at 300 V are especially unique, however, and future work is required to further explore mechanisms, which may drive the time-dependent ion heating observed at this regime. The observed heating patterns are inconsistent with ionization region overlap or collisional effects; therefore, it is possible that an anomalous ion heating process occurs in phase with the low-frequency current oscillations. Investigating the details of such a phenomenon may shed light on the observed non-classical erosion rates at 300 V as well as characterize possible thermal consequences of Hall thruster anomalous transport models. In both configurations, spatially denser time-resolved measurements would allow for more accurate temperature determination and the observation of the short length scale at which the instantaneous IVDFs vary in shape. Similar observations near the magnetic poles are necessary to fully characterize the extent to which time-dependent ion heating affects erosion.

## ACKNOWLEDGMENTS

This research was carried out at the Jet Propulsion Laboratory, California Institute of Technology, under a contract with the National Aeronautics and Space Administration (No. 80NM0018D004). The support of the joint NASA GRC and JPL development of HERMeS by NASA’s Space Technology Mission Directorate through the Solar Electric Propulsion Technology Demonstration Mission (SEP TDM) is gratefully acknowledged. P.J.R. acknowledges funding from the JPL Education Office Summer Internship Program (SIP).

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

## DATA AVAILABILITY

Subject to the constraints of United States Government export controls (ITAR and EAR), the data that support the findings of this study can be made available from the authors upon reasonable request.

### APPENDIX: TEMPERATURE OF A BI-MAXWELLIAN DISTRIBUTION

This appendix contains the analytical derivation of the effective ion temperature $Ti$ of a plasma with a bi-Maxwellian ion velocity distribution function (IVDF). The kinetic temperature of a 1D plasma is defined via the second moment of the distribution function in a reference frame moving with the mean flow velocity,

For a bi-Maxwellian plasma with IVDF in the form of Eq. (2), $v\xaf$ is first computed by dividing the first moment of the distribution by the zeroth moment,

Equation (A3) differs from the expression for temperature of a single Maxwellian in that it not only depends on the thermal velocities, but also the separation between the mean velocities of the two distributions and the overall flow speed.

## REFERENCES

*46th AIAA/ASME/SAE/ASEE Joint Propulsion Conference and Exhibit, Nashville, TN*(AIAA, 2010), AIAA 2010-6698.

*AIAA Propulsion and Energy Forum, 2018 Joint Propulsion Conference, Cincinnati, OH*(AIAA, 2018), AIAA 2018-4647.

*35th International Electric Propulsion Conference, Atlanta, GA*(ERPS, 2017), IEPC-2017-154.

*Propulsion and Energy Forum, 52nd AIAA/SAE/ASEE Joint Propulsion Conference, Salt Lake City, UT*(AIAA, 2016), AIAA 2016-4839.

*AIAA Propulsion and Energy 2020 Forum—Virtual*(AIAA, 2020), AIAA 2020-3620.

*33rd International Electric Propulsion Conference*(ERPS, 2013).

*34th International Electric Propulsion Conference*(ERPS, 2015), IEPC-2015-155.

*36th International Electric Propulsion Conference, Vienna, Austria*(ERPS, 2019).

*36th International Electric Propulsion Conference*(ERPS, 2019), IEPC-2019-651.

*AIAA Propulsion and Energy Forum, 2018 Joint Propulsion Conference, Cincinnati, OH*(AIAA, 2018), AIAA 2018-4645.

*AIAA Propulsion and Energy Forum, 2018 Joint Propulsion Conference, Cincinnati, OH*(AIAA, 2018), AIAA 2018-4723.

*36th International Electric Propulsion Conference, Vienna, Austria*(ERPS, 2019), IEPC-2019-516.