This work focuses on a single illustrative operating condition for the unshielded H6 Hall thruster. The H6 is a laboratory model 6-kW thruster (see Figure 1) with a design operating point of 300 V and 20 mg/s (where 1 mg/s ≈ 1 A) on xenon. In this work, we focus on a 600 V and 10 mg/s condition that exhibits several simultaneous modes, most prominently an m = 1 cathode spoke, a shared m = 0 mode between the channel and cathode regions, and mode-hopping between m = 3 and m = 4 in the channel.

Recordings took place in the University of Michigan Plasmadynamics and Electric Propulsion Laboratory’s Large Vacuum Test Facility (LVTF) circa 2011. The LVTF is a cylindrical chamber 9 m long and 6 m in diameter and at the time was maintained at high vacuum by seven TM-1200 cryopumps with a combined pumping speed of 210,000 L/s on xenon.

A Photron SA5 FASTCAM high-speed camera placed 6.5 m axially downstream of the thruster exhaust imaged the plume through a quartz window with 152 × 192 resolution in monochrome (B&W) at a framerate of 175 kHz. The fastest mode observed was a 80 kHz cathode spoke, just below the 87.5 kHz camera Nyquist frequency. The camera used a Nikon ED AF Nikkor 80–200 mm lens at its maximum aperture of f/2.8. The bright, central cathode saturates several pixels, and the algorithms largely ignore these pixels.

In preparation for video processing, any high-speed video dataset should be thought of as a 3D matrix of pixel measurements, p(t, x, y), with dimensions of time, t, and Cartesian space, x and y, and with lengths N _t , N _x , and N _y , respectively (McDonald and Gallimore, 2013).

Before applying mode analysis algorithms, several preprocessing steps should first be considered. This section outlines two prominent steps: 1) spatial identification and normalization and 2) amplitude normalization. Other preprocessing steps, not covered here, include masking and filtering which are useful in isolating specific spatial regions or dynamics, respectively. Please note that Fourier analysis requires spatial identification as a precursor to converting to polar coordinates. The other methods discussed here are useful, but not required, for processing and analysis.

Mode dynamics in Hall thruster plasmas are typically characterized by oscillating waves with slowly evolving amplitudes and frequencies. The easiest way to detect these modes is to apply FFT and Welch-averaged FFT (Welch, 1967) algorithms to diagnostic measurements of signal pixels of the high-speed video. Figure 4 shows an example of this applied to three signals from our dataset: the discharge current flowing from the anode to cathode, a high-speed pixel centered in the channel, and a high-speed pixel adjacent to the cathode.

Figure 4A shows a 1 ms time window of the three signals and that they have similar oscillatory behavior. Each signal has been AC coupled and normalized by its standard deviation for ease of comparison.

Figure 4B shows the Welch-averaged FFT of each signal over 100 ms, and the resulting power spectrum reveals prominent frequency peaks with each. Most notable is a broad spectrum peak at 20 kHz that is common to all three signals. Additional peaks at 7.5 kHz and 80 kHz are unique to the channel and cathode measurements, respectively; this uniqueness suggests that their dynamics are isolated to their respective regions. The broader width of the 20 kHz peak suggests that its frequency is more erratic than the two narrower peaks. The remainder of this paper will identify the modes associated with these three peaks and few less pronounced peaks.

Waves around the channel and cathode of a Hall thruster are primarily azimuthal (Choueiri, 2001) with the approximate Fourier form p t , θ ∼ e i m θ − ω t (3) where m is the integer azimuthal mode number, θ is the azimuthal angle, and ω = 2πf is the angular frequency. After scaling the video’s coordinates and amplitudes (Section 2.2), the next step is convert the 3D video, p(t, x, y), in Cartesian coordinates to a 2D video, p(t, θ), in polar coordinates so that it matches Eq. 3. To do this, we first isolate the channel region (i.e. 0.9 < r/r ₀ < 1.1) with radial masking and discard the rest. The radial dependence within this narrow region is assumed constant, and the radial coordinate is therefore dropped.

Next, we convert the unmasked x and y coordinates within the channel to the azimuthal coordinate with a four-quadrant inverse tangent function, θ(x, y) = arctan2(y, x). The resulting data is azimuthally binned and averaged to provide a result with uniform azimuthal spacing.

Figure 5 illustrates the azimuthal and radial binning for a single time snapshot (t = 22.9 ms). Figure 5A shows the radial mask isolating a narrow region within the channel and also the edges of N _θ = 100 azimuthal bins; Figure 5B shows the unbinned values (black) overlayed by the bin-averaged result (red). This shows a clear m = 3 structure and also a few non-ideal features: a non-uniform spacing between peaks, a steepened waveform, and a weak fourth peak at 3π/4. This process is repeated for each time step within p(t, x, y), and the result is the 2D dataset, p(t, θ), within the channel.

In this example, we applied the binning process to the channel region. However, this same procedure could be applied to other radial regions, including the region around the cathode (Jorns and Hofer, 2014) or at radial slices between the cathode and channel. Optionally, multiple radial slices could be made to provide a 3D dataset, p(t, θ, r).

With p(t, θ) solved and matching the assumed Fourier form (Eq. 3), we can next identify the azimuthal modes within the channel. To do this, we apply a 2D FFT, F t , θ , in both time and θ to identify the mode numbers and their characteristic frequencies. This provides the 2D complex matrix, P ω , m = F t , θ p t , θ , (4) with dimensions in angular frequency, ω, and azimuthal mode number, m, and with dimensional lengths, N _t and N _θ , respectively. Note that the coordinates ω = 2πf and m range between their negative and positive Nyquist frequencies (−87.5 kHz ≤ f < 87.5 kHz and −N _θ /2 ≤ m < N _θ /2). Due to symmetry, the negative frequencies can be truncated.

The absolute value of this result, P ( ω , m ) , is the azimuthal dispersion relationship within the channel and is shown in Figure 6A. It reveals a series of evenly-spaced discrete modes (2 ≤ m ≤ 5) and what appears to be a continuous wave. Both the continuous wave and the discrete modes are propagating clockwise (negative θ) and therefore have negative wavenumbers. By convention, we present the azimuthal modes and wave numbers as positive.

Figure 6B shows several slices of the dispersion relationship at m = 0, 2, 3, and 4 and more clearly identifies the peaks originally observed in Figure 4B. From this plot, the m = 0 mode is the broad-spectrum peak at 20 kHz, the m = 3 mode is the narrow peak at 7.5 kHz, and the m = 4 mode is the peak at 12.5 kHz. In addition, a weak m = 2 mode is observed at 2.5 kHz. Figure 6B shows that all three of the indicated m > 0 modes have a roughly uniform spacing of 5 kHz. Figure 6A also identifies an m = 6 mode at roughly double the m = 3 frequency which makes it a harmonic of the m = 3 mode (Yamada et al., 2010). To observe the 80 kHz peak in Figure 4B, the above analysis could be applied to the region around the cathode instead of inside the channel.

To capture the time evolution of each mode, we apply 1D FFT, F θ , to each time step in the 2D dataset, p(t, θ), along the θ dimension. The result is the 2D complex matrix, P t , m = F θ p t , θ (5) with dimensions in time, t, and azimuthal mode number, m. Depending on the FFT algorithm, p(t, m) then needs to be multiplied by a constant to return the correct amplitude, typically 2/N _θ . Figure 7 plots the real (cosine) component, the imaginary (sine) component, and the amplitude of p(t, m) at m = 0, 3, and 4 and shows the time evolution of each. Figure 7A reveals the m = 0 breathing mode to have a mostly consistent amplitude around 10%–20% of the average channel brightness. Figures 7B, C shows the m = 3 and m = 4 modes, respectively, with amplitudes between 3% and 5%. These figures also clearly identify repeated mode hopping between these two modes. In the m = 3 and 4 figures, the real component leads the real component by roughly 90°, indicating that the modes are rotating clockwise. The m = 0 mode is not rotating, and therefore its imaginary component is zero.

To isolate the spatial structure associated with each mode, we first apply 1D Fourier analysis, P f , x , y = F t p t , x , y , (6) to each pixel in the 3D video, p(t, x, y), with respect to time. Next, we index the resulting matrix, p(f, x, y), at the frequencies associated with each mode as identified in Figure 6. Figure 8 shows the resulting real (cosine) and imaginary (sine) component of each mode, which reveals the spatial extent and phase of each. The number of peaks and troughs of each structure accurately corresponds to its mode number. This approach is a computationally efficient alternative to bandpass-filtering each pixel as done previously (McDonald and Gallimore, 2013). These plots are also nearly identical to plots in previous works (Jorns and Hofer, 2014; Baird, 2020) with the difference being that the past works plot the phase, atan (imag./real), of each mode instead of the real and imaginary components.

For better visualization, each pixel in p(t, x, y) can optionally be normalized by its standard deviation before applying Eq. 6. Figure 9 shows an example of this for the m = 3 mode. This figure reveals improved spatial detail within the channel and also suggests mode structure outside the channel as well.

In this section, we applied the Fourier mode decomposition and identification methods to our H6 high-speed video dataset. First, the FFT of individual pixels was able to identify the presence of coherent mode dynamics in the channel and cathode regions. After radially and azimuthally binning the data within the channel, a 2D FFT provided an azimuthal dispersion plot and helped relate mode numbers to their characteristic frequency. A 1D FFT was then used to isolate the temporal evolution of each mode and most notably identified mode hopping between the m = 3 and m = 4 modes in the channel. Finally, a 1D FFT was applied to the original video, and plotting the frequencies associated with each provided the mode’s spatial structure. While these methods were only applied to the channel region in this work, this process can be applied to any radial region including around the cathode.

The SVD algorithm is the foundation for POD and for more advanced algorithms. This section discusses the SVD algorithm and its implementation.

The most notable feature of SVD is that it does not presume spatial or temporal bases (e.g., Fourier sines and cosines) when decomposing coherent structures. Instead, SVD provides a tailored, orthonormal basis for a particular video dataset by minimizing the L ₂ norm (error) between the original video and the video reconstructed from the new SVD bases (Brunton and Kutz, 2019). SVD also orders each basis from the highest to lowest contribution (amplitude) so that the strongest coherent structure appears first, then the second strongest, etc. In contrast, Fourier orders its bases by increasing frequency or wavenumber. A downside to SVD’s ordering is that post-processing (either visually or algorithmically) is often required to associate each identified SVD structure with its Fourier counterpart. Examples of this are provided later in this paper. An advantage of SVD is that it does not require any preprocessing, most notably conversion to a particular coordinate system (e.g., converting the Cartesian video to polar). More details on the SVD algorithm can be found in the extensive literature on the subject (Golub and Reinsch, 1970; Golub and Loan, 1996; Woolfe et al., 2008; Brunton and Kutz, 2019).

In order to apply the SVD algorithm to a video dataset, p(t, x, y), we first convert it from 3D to 2D. We do this by stacking the data associated with the two spatial dimensions into a single spatial dimension, z = stack(x, y) to get the 2D p(t, z) dataset with dimensions in time and space and with sizes N _t and N _z = N _x N _y , respectively. This process is the 3D equivalent of stacking columns in a 2D matrix to get a 1D array. The Python data structure xarray has built-in functions (stack and unstack) that make stacking and unstacking very convenient. By convention, we also transpose p(t, z) to p(z, t) so that the spatial dimension is ordered before the temporal dimension. For our high-speed video dataset, N _z > N _t , so p (z, t) is a “tall-skinny” matrix.

When applied to p(z, t), the SVD algorithm outputs three real matrices: U(z, n), Σ(n, n), and V(t,n) ^T . Here, we are using n to represent the SVD basis numbers to distinguish them from the Fourier mode numbers, m. First, U(z, n) is a non-square, 2D matrix associated with the spatial bases (called the topos, i.e., “shapes”) with dimensions of space, z, and mode number, n, and lengths N _z and N _n = N _t , respectively. Σ(n, n) is a square 2D diagonal matrix associated with the basis energies (these diagonal values are referred to as the singular values) with dimensions n by n each with length N _n . Finally, the transposed matrix, V(t,n) ^T , is a square 2D matrix that contains the temporal evolution of the bases (called the chronos, i.e., “times”) with dimensions of n and t, each with length, N _n = N _t . For example, many SVD algorithms return Σ(n, n) as 1D array of its diagonal elements.

Multiplying these three matrices together, U z , n Σ n , n V t , n T = | | | u 1 u 2 … u N | | | σ 1 ⋱ σ N | | | v 1 v 2 … v N | | | T = σ 1 u 1 v 1 T + σ 2 u 2 v 2 T + ⋯ + σ N u N v N T ≈ p z , t , (7) reconstructs p(z, t) with remarkable accuracy, owing to the orthonormality of the u _n and v _n and the ordering of bases from strongest to weakest. Related to this, another major advantage of SVD over Fourier is its ability to reconstruct data (like this high-speed video) with fewer bases. In this equation, u _n , σ _n , and v _n , are the topo, singular value (energy), and chrono associated with the nth basis, which are ordered from highest contribution (energy) to lowest. In Eq. 7, the multiplication of u n v n T is an outer product. As a final note, we refer to σ _n as the energy because our measurement is light intensity. However, it is common in the literature to measure the velocity or magnetic field and therefore refer to σ n 2 as the energy.

The POD method, also referred to as Biorthogonal Decomposition (BD), has been used extensively for plasma physics datasets (Dudok de Wit et al., 1994) but has only been very recently introduced to Hall thruster and hollow cathode high-speed video (Désangles et al., 2020; Becatti et al., 2021). The POD algorithm is nearly identical to the SVD algorithm with the main distinction being that the U, Σ, and V ^T matrices are often trimmed to retain a smaller set of bases (e.g., n < 30).

After applying the POD algorithm (Section 4.1) to our dataset, we inspect each of the U, Σ, and V ^T matrices to identify and interpret the dynamics associated with each mode. First, the energy for each basis, σ _n , is plotted in Figure 10. This highlights that the lower numbered bases capture the majority of the dynamics of the video and are therefore the focus of this section.

To visualize the spatial bases (topos), we unstack the z dimension in U(z, n) to get U(x, y, n). Figure 11 shows the first 10 topos, u _n , and each reveals an m = 0 or m > 0 azimuthal structure in the channel or around the cathode. The bases associated with m > 0 rotating modes each have a near-duplicate (but rotated) topo that represents an effective sine-cosine pairing. For example, the n = 2 and 3 bases are effectively the Fourier sine and cosine (real and imaginary) components of the m = 3 rotating mode in the channel and therefore well matches the modes in Figure 8. Note that the POD bases associated with the m = 2 Fourier mode is not shown in Figure 11 despite being present (as shown by Figure 8). This is because its energy contribution is lower than the first 10 POD bases.

Each of these topos has an associated time evolution, i.e., chronos, associated with it. Several select chronos, σ n v n T , are plotted in Figure 12. Referring back to Figure 11, we see that their corresponding topos are roughly equivalent to the m = 0, 3, and 4 Fourier modes in the channel. Therefore, it is not a surprise that their chronos in Figure 12 are nearly identical to the Fourier-solved mode evolution shown in Figure 7.

Finally, the Welch-averaged FFT of select POD chronos are shown in Figure 13. Note that only a single basis from a sine-cosine pair is presented as their power spectrums of each are nearly identical. Many of the peaks in Figure 13 are the same peaks as identified in the raw data (Figure 4B) and the Fourier mode analysis (Figure 6B). A new structure is revealed by the n = 4 basis, related to the Fourier m = 1 cathode mode, which has two peaks around 38 and 80 kHz. From this analysis alone, it is not clear if the m = 1 cathode mode truly has two characteristic frequencies, if this basis is the combination of two Fourier modes, or if there is another explanation.

Comparing the results of the POD bases (Figures 11–13) to the Fourier modes (Figures 6–8) reveals a near one-to-one mapping as is indicated in Table 1. This is not too surprising as there are already established scenarios in which POD modes reduce to Fourier modes (Holmes et al., 2012).

The primary advantage of the SVD-based methods is that they do not presume a universal basis, Fourier or otherwise, and instead create tailored bases which optimally capture particular features of the dynamics. In addition, these methods require less preconditioning of the video, e.g., radial binning and converting the pixel coordinates to polar, and therefore they can provide detailed images of the spatial structures of each basis or mode (Figure 11). The most basic SVD method, POD, is easier to implement than the Fourier methods, does not require a linear assumption, and produces very similar results. More advanced SVD methods, such as DMD, have potential for Hall thruster modal analysis but require more development.