# A comparison of Fourier and POD mode decomposition methods for high-speed
Hall thruster video

J. W. Brooks & M. S. McDonald  Naval Research Laboratory, Washington DC
<EMAIL_ADDRESS>&<EMAIL_ADDRESS>A. A. Kaptanoglu
University of Washington, Seattle WA<EMAIL_ADDRESS>
###### Abstract

Hall thrusters are susceptible to large-amplitude plasma oscillations that
impact thruster performance and lifetime and are also difficult to model.
High-speed cameras are a popular tool to study these dynamics due to their
spatial resolution and are a popular, nonintrusive complement to in-situ
probes. High-speed video of thruster oscillations can be isolated (decomposed)
into coherent structures (modes) with algorithms that help us better
understand the evolution and interactions of each. This work provides an
introduction, comparison, and step-by-step tutorial on established Fourier and
newer Proper Orthogonal Decomposition (POD) algorithms as applied to high-
speed video of the unshielded H6 6-kW laboratory model Hall thruster. From
this dataset, both sets of algorithms identify and characterize $m=0$ and
$m>0$ modes in the discharge channel and cathode regions of the thruster
plume, as well as mode hopping between the $m=3$ and $m=4$ rotating spokes in
the channel. The Fourier methods are ideal for characterizing linear modal
structures and also provide intuitive dispersion relationships. By contrast,
the POD method tailors a basis set using energy minimization techniques that
better captures the nonlinear nature of these structures and with a simpler
implementation. Together, the Fourier and POD methods provide a more complete
toolkit for studying Hall thruster plasma instabilities and mode dynamics.
Specifically, we recommend first applying POD first to quickly identify the
nature and location of global dynamics and then using Fourier methods to
isolate dispersion plots and other wave-based physics.

Keywords: Hall thruster, high-speed video, mode decomposition, singular value
decomposition

##  1 Introduction

The Hall thruster (HT) is a standard spacecraft electric propulsion system
that uses crossed electric and magnetic fields
$\left(\mathbf{E}\times\mathbf{B}\right)$ to ionize and accelerate propellant
[1, 2]. Their high specific impulse and technological maturity make them ideal
for long duration satellite station keeping and time-insensitive missions,
such as orbit raising, with thousands currently in orbit and more planned.
However, HTs experience anomalous electron transport across their magnetic
field that is not sufficiently understood to permit fully predictive thruster
models, motivating continued improvement of diagnostics for model validation.
This lack of validation is especially important for increased qualification by
simulation for new thruster designs at ever-higher power, such as the N30 or
X3 nested Hall thrusters, where full life qualification in a ground test
facility would be prohibitively expensive for 10s-100s kW thrusters [3, 4].

Anomalous electron transport has been strongly linked to plasma oscillations
in HTs [5] in both the thruster discharge channel [6, 7] and cathode regions
[8]. Long-wavelength azimuthal oscillations, the focus of this work, are
typically characterized in a Fourier representation with integer mode numbers,
$m$, including the azimuthally uniform $m=0$ “breathing” modes in HT channels
[9, 10] and cathodes [11, 12], $m>0$ azimuthally rotating “spoke” modes in the
channel, and a $m=1$ counter rotating “anti-drift” modes around the
cathode[13]. The rotating spoke channel modes, first observed in an early Hall
accelerator by Janes and Lowder [14], are now known to be ubiquitous in
unshielded thrusters[15] though they may be less common in magnetically
shielded thrusters [13, 16]. Cathode $m=1$ anti-drift modes have been
previously observed in both shielded and unshielded versions of the H6 Hall
thruster [13] (the unshielded H6 is studied in this work). All of these modes
are historically characterized by one or more in-situ diagnostics, but the
localized nature of these diagnostics both perturb the plasma and make it
difficult to relate individual measurements to global mode structures.

High-speed imaging (HSI) has become a popular complementary diagnostic to in
situ probes and has proven itself well suited to characterizing mode dynamics
[6, 17, 13]. Its popularity is due to its ease of use, nonintrusive nature,
high speed (100s of kHz), and fine spatial resolution of order millimeters per
pixel. HSI was first used on HTs to relate image brightness to plasma density
oscillations measured by electrical probes [18]. This and later work has led
to the use of pixel light intensity as a proxy for discharge current [19].

To isolate and characterize individual mode dynamics from HSI, various post-
processing algorithms have been developed. Early work subtracted the time-
average from each pixel to reveal multiple $m>1$ modes in the H6 Hall thruster
channel [15]. Later, Fourier techniques were used with azimuthal binning to
isolate the frequency spectrum associated with individual channel modes [20].
Fourier methods have also been used to provide azimuthal dispersion plots
within the CHT (cylindrical Hall thruster) [6] and H6 [17] thruster channels.
A phase-based Fourier analysis was developed to isolate the spatial structure
of the $m=1$ cathode mode on the H6 thruster [13] and $m>0$ channel and
cathode modes on the HERMeS thruster [16]. An alternative Fourier-based
visualization method, Cross-Spectral-Density (CSD), was developed on the CHT
thruster [21].

The Fourier methods remain the most established method for isolating and
analyzing Hall thruster plasma oscillations, but they have several
disadvantages. First, their linear sine/cosine bases are not ideal for
nonlinear features. Second, Fourier methods require several preprocessing
steps for high-speed video; this added complexity requires additional
computational overhead and makes it more difficult to study multi-dimensional
dynamics.

An alternative to Fourier methods are Singular Value Decomposition (SVD) based
algorithms [22]. The chief advantage of SVD is that is tailors custom bases
for each dataset based on the energy contribution of the coherent dynamics
instead of using a presumed basis (e.g. Fourier’s sines and cosines). This
facilitates improved characterization and reconstruction of nonlinear dynamics
without making any physical assumptions of the system. SVD also requires
minimal preprocessing compared with Fourier methods.

The most notable and likely simplest SVD algorithm is Proper Orthogonal
Decomposition (POD) which is extensively used in fluid mechanics [23, 24, 22].
A classic POD application is to decompose fluid vortex shedding around a body
(e.g. a cylinder) into discrete modes [25, 26]. POD has also been used in
plasma physics to characterize plasma oscillations [27, 28, 29, 30, 31] but is
often referred to as Biorthogonal Decomposition (BD). Very recently, POD has
been used to characterize axial and azimuthal modes in Hall thruster high-
speed video [32] and azimuthal modes in a hollow cathode plume [33]. More
sophisticated SVD algorithms exist, most notably DMD (Dynamic Mode
Decomposition) [34, 35] and have applications in active control, linear
dynamics, and plasma physics [36, 37, 38].

The goal of this work is to compare Fourier and POD techniques as applied to
Hall thruster high-speed imaging in a tutorial format. To this end, this work
analyzes a high-speed video recording of the unshielded H6 Hall thruster plume
and provides a step-by-step explanation of algorithm’s implementation and
results. This work starts by introducing the H6 thruster, high-speed video
dataset, and video preprocessing in Section 2. Section 3 discusses several
established Fourier mode analysis methods and, when applied to the H6 dataset,
identifies simultaneous cathode and channel modes in addition to mode hopping
within the channel. Section 4 introduces the SVD algorithm and its most common
implementation: Proper Orthogonal Decomposition (POD). When applied to the H6
dataset, POD identifies the same mode behavior and provide several
improvements over the Fourier methods. When used together, POD and Fourier
methods provide a more complete toolkit for identifying and the isolating
global mode dynamics. The code and dataset used for this work are available
online [39].

##  2 Experimental setup

This section covers the experimental hardware (thruster, facility, high-speed
camera), the HSI dataset, and common HSI prepossessing techniques.

###  2.1 Hardware and dataset

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
$\approx$ 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.

Figure 1:  a) The unshielded H6 Hall thruster. b) End-on-view of the thruster
and the plasma in its channel (annular ring) and the cathode (center dot).
Figure reproduced here with permission from the author [40].

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 x 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 [17].

###  2.2 Video preprocessing

Before applying mode analysis algorithms, several preprocessing steps should
first be considered. This section outlines two prominent steps: 1) spatial
centering 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 all of these steps are optional for POD.

####  2.2.1  Spatial identification and scaling

This step identifies and the scales the spatial geometries in preparation for
converting to polar coordinates, a requirement for Fourier analysis. This step
is not required for POD. To identify the thruster channel origin and radius,
we fit an annular Gaussian function

$\displaystyle\eqalign{G&(x,y;x_{0},y_{0},r_{0},a,w,G_{0})\\\
&=a\,exp\left(-\frac{1}{2}\left(\frac{r(x,y;x_{0},y_{0})-r_{0}}{w}\right)^{2}\right)+G_{0}}$
(3)

to the time-average, $\bar{p}(x,y)$, of the high-speed video. In this
equation, $r_{0}$ is the radius of the channel, the radius at each pixel is

$r(x,y;x_{0},y_{0})=\sqrt{(x-x_{0})^{2}+(y-y_{0})^{2}},$ (4)

$x_{0}$ and $y_{0}$ are the center (origin) of the channel, $a$ is the
amplitude, $w$ is related to the channel width, and $G_{0}$ is an offset.

With these fit parameters solved, we next center and normalize the spatial
coordinates to the channel radius (i.e. $x_{norm}=(x-x_{0})/r_{0}$ and
$y_{norm}=(y-y_{0})/r_{0}$ ). These normalized coordinates can then be
optionally multiplied by the dimensioned channel radius to provide actual
units to $x$ and $y$. For this work, we remain with normalized coordinates.

Figure 2 shows the results of this step as applied to the time-averaged video,
$\bar{p}(x,y)$, of our H6 dataset. In this figure, the two regions with the
most plasma dynamics, the annular channel and center cathode, are clearly
visible. The pixel intensity has raw integer units based on the camera’s 12
bit depth. The pixels associated with the cathode’s center were allowed to
saturate to better capture dynamics in the channel, and the pixels adjacent to
the saturated pixels still capture the cathode’s dynamics. The channel radius,
$r_{0}$, as identified by the Gaussian fit, is indicated with a black dashed
line. The channel edges, indicated with white dashed lines, are identified as
approximately at $(r_{0}\pm w)/r_{0}$. The $x$ and $y$ coordinates have been
centered and normalized to the channel radius as described above.

Figure 2:  An annular Gaussian function is fit to the time-averaged video and
identifies the channel origin and radius. The video’s coordinates are then
centered at the origin and normalized by the channel radius.

Alternatives to the annular Gaussian fit have been used in previous Hall
thruster work. McDonald [17] discussed both the Kasa [41] and Taubin [42]
methods and recommended Taubin for cases where the entire channel annulus is
not visible. Another option previously used [21] is a circle detection
algorithm called a Hough transform. While no one method is obviously superior
over the others, we recommend using 1) the Gaussian fit presented here because
the solved parameters are directly relatable to physical dimensions or 2) the
Hough transform as it is available prewritten in many programming languages.

####  2.2.2  Video amplitude scaling

This step scales the video’s arbitrary intensity measurements to a more
meaningful range. Common practice is to assume that oscillations in a video’s
amplitude is roughly linear with the plasma density oscillations [18] and that
camera measurements are linear with light emission [43]. As we cannot scale
the data to a definitive physical value, this section instead discusses
several data normalization methods to provide better physical intuition of the
oscillations. This step is optional for both Fourier and POD methods.

Before normalizing, the first step is to subtract the time-averaged image
(also known as AC coupling) from each frame of the raw video dataset to better
isolate the oscillations.

Next, we normalize the video amplitude in one of two ways: pixel-wise
normalization or channel-average normalization. The first method divides each
pixel by its standard deviation in time, and this has the advantage of making
the mode dynamics easier to visualize after mode decomposition. Unfortunately,
this method artificially amplifies the oscillations at different spatial
locations and makes quantifying global mode amplitudes untenable. As an
alternative, the second option divides each pixel by the average brightness
within the channel, which allows the modes to be scaled as a percentage of the
average channel brightness. While both methods are used in this work, channel-
average normalization is default.

Figure 3 shows the results of AC coupling and channel normalization at two
separate instances in time. At $t=22.9$ ms (Figure 3a), the normalized video
snapshot reveals a dominant 3-lobed azimuthal wave ($m=3$ mode) in the
channel. At $t=24.8$ ms (Figure 3b), a 4-lobed azimuthal wave ($m=4$ mode) is
dominant. The amplitudes of both modes are around 10% of the average channel
brightness. Both snapshots also show an $m=1$ azimuthal wave around the
cathode.

Figure 3:  Two time snapshots, that have been AC coupled and normalized by the
average channel brightness, reveal prominent mode structures. a) At $t=22.9$
ms, an $m=3$ mode is dominant in the channel. b) At $t=24.8$ ms, an $m=4$ mode
is dominant in the channel. Both snapshots show an $m=1$ mode around the
cathode.

##  3 Fourier-based methods

Fourier based algorithms are the most established methods for mode
decomposition and identification in HT high-speed video [44, 17, 13, 21]. This
is because Fourier series’ bases are periodic sines and cosines and are
therefore ideal for characterizing wave-like oscillations, such as plasma
waves. In this section, we use Fourier analysis to characterize simultaneous
$m=0$ and $m>0$ modes associated with the thruster discharge channel and
cathode in addition to mode hopping between the $m=3$ and $m=4$ modes in the
channel.

###  3.1 Detecting modes

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 [45]
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 4:  Three signals are analyzed (a pixel in the thruster channel
($r/r_{0}=1$), a pixel adjacent to the thruster cathode ($r/r_{0}=0.1$), and
the discharge current), and their power spectrum reveals the existence of
several prominent modes (peaks). a) The three time-series signals over 1 ms.
Each has been subtracted by their mean and divided by their standard
deviation. b) Their power spectrums, calculated over 100 ms.

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.

###  3.2 Azimuthal and radial binning

Waves around the channel and cathode of a Hall thruster are primarily
azimuthal [5] with the approximate Fourier form

$p(t,\theta)\sim e^{i(m\theta-\omega t)}$ (5)

where $m$ is the integer azimuthal mode number, $\theta$ is the azimuthal
angle, and $\omega=2\pi 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,\theta)$, in
polar coordinates so that it matches Eq. 5. To do this, we first isolate the
channel region (i.e. $0.9<r/r_{0}<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,
$\theta(x,y)=\text{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_{\theta}=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$\pi$/4. This process is repeated for each time step
within $p(t,x,y)$, and the result is the 2D dataset, $p(t,\theta)$, within the
channel.

Figure 5:  The azimuthal binning process for a single instant in time is
shown. a) A radial mask isolates a narrow radial region within the channel,
and the $N_{\theta}$=100 azimuthal bin boundaries are overlayed. b) The
average of each azimuthal bin (red) overlays the raw (black) data.

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 [13] or at radial slices between the
cathode and channel. Optionally, multiple radial slices could be made to
provide a 3D dataset, $p(t,\theta,r)$.

###  3.3 Azimuthal mode identification

With $p(t,\theta)$ solved and matching the assumed Fourier form (Eq. 5), we
can next identify the azimuthal modes within the channel. To do this, we apply
a 2D FFT, $\mathcal{F}_{t,\theta}$, in both time and $\theta$ to identify the
mode numbers and their characteristic frequencies. This provides the 2D
complex matrix,

$\displaystyle
P(\omega,m)=\mathcal{F}_{t,\theta}\bigl{\\{}p(t,\theta)\bigr{\\}},$ (6)

with dimensions in angular frequency, $\omega$, and azimuthal mode number,
$m$, and with dimensional lengths, $N_{t}$ and $N_{\theta}$, respectively.
Note that the coordinates $\omega=2\pi f$ and $m$ range between their negative
and positive Nyquist frequencies (-87.5 kHz $\leq f<$ 87.5 kHz and
$-N_{\theta}/2\leq m<N_{\theta}/2$). Due to symmetry, the negative frequencies
can be truncated.

The absolute value of this result, $\bigl{|}P(\omega,m)\bigr{|}$, is the
azimuthal dispersion relationship within the channel and is shown in Figure
6a. It reveals a series of evenly-spaced discrete modes ($2\leq m\leq 5$) and
what appears to be a continuous wave. Both the continuous wave and the
discrete modes are propagating clockwise (negative $\theta$) 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 [46]. 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.

Figure 6:  a) The azimuthal dispersion plot within the channel. A continuous
wave and series of azimuthal mode numbers are identified and travel counter-
clockwise (i.e. are negative). b) The power spectrum of select mode numbers
helps to identify their characteristic frequency.

###  3.4 Mode evolution

To capture the time evolution of each mode, we apply 1D FFT,
$\mathcal{F}_{\theta}$, to each time step in the 2D dataset, $p(t,\theta)$,
along the $\theta$ dimension. The result is the 2D complex matrix,

$\displaystyle P(t,m)=\mathcal{F}_{\theta}\bigl{\\{}p(t,\theta)\bigr{\\}}$ (7)

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_{\theta}$. 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% to 20% of the average channel brightness. Figures 7b and 7c shows the
$m=3$ and $m=4$ modes, respectively, with amplitudes between 3% to 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 degrees, indicating that the modes are rotating clockwise. The
$m=0$ mode is not rotating, and therefore its imaginary component is zero.

Figure 7:  Time evolution of the $m=0$, 3 and 4 modes as captured by the 1D
FFT in $\theta$. Repeated mode hopping is observed between the $m=3$ and $m=4$
modes. The signals are normalized by the average channel brightness.

###  3.5 The modes’ spatial structures

To isolate the spatial structure associated with each mode, we first apply 1D
Fourier analysis,

$\displaystyle P(f,x,y)=\mathcal{F}_{t}\bigl{\\{}p(t,x,y)\bigr{\\}},$ (8)

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 [17]. These plots are also nearly identical to plots in previous
works [13, 16] with the difference being that the past works plot the phase,
atan(imag./real), of each mode instead of the real and imaginary components.

Figure 8:  Fourier analysis (Eq. 6) reveals the spatial structure of several
dominant modes in the channel. The amplitude is spectral density (au) with
each subfigure having a different scaling.

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.

Figure 9:  Pre-normalizing each pixel in $p(t,x,y)$ by its standard deviation
instead of by the average channel brightness provides better visualization of
the mode structure. Amplitude is spectral density (au). These plots are
directly comparable to the $m=3$ subplots in Figure 8 where the channel
normalization is used.

###  3.6 Fourier methods conclusion

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 first 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.

##  4 SVD-based methods

Matrix factorization-based methods, based on the SVD (Singular Value
Decomposition) algorithm [47], are an alternative approach to decomposing and
identifying modes from high-speed video and provide several advantages over
Fourier methods. Below, we first detail the SVD algorithm, discuss the
simplest SVD-based mode analysis method (the POD algorithm), and apply it to
the present dataset.

###  4.1 The SVD algorithm

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_{2}$ norm (error) between the
original video and the video reconstructed from the new SVD bases [22]. SVD
orders these bases from the highest to lowest contribution to the
reconstruction instead of by increasing wavenumber or frequency. A downside to
this is that the bases often need to be associated with their Fourier
counterparts through post-processing (manually or algorithmically). 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 [47, 48, 49, 22].

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=\text{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)$, $\Sigma(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.
$\Sigma(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}$. The exact dimension and format
of these three matrices may vary depending on the particular SVD algorithm
being used and shape of $p(z,t)$. For example, many SVD algorithms return
$\Sigma(n,n)$ as 1D array of its diagonal elements.

Multiplying these three matrices together,

$\begin{split}&U(z,n)\Sigma(n,n)V(t,n)^{T}\\\ &=\begin{bmatrix}|&|&&|\\\
u_{1}&u_{2}&\ldots&u_{N}\\\ |&|&&|\end{bmatrix}\begin{bmatrix}\sigma_{1}&&\\\
&\ddots&\\\ &&\sigma_{N}\end{bmatrix}\begin{bmatrix}|&|&&|\\\
v_{1}&v_{2}&\ldots&v_{N}\\\ |&|&&|\end{bmatrix}^{T}\\\
&=\sigma_{1}u_{1}v_{1}^{T}+\sigma_{2}u_{2}v_{2}^{T}+\ldots+\sigma_{N}u_{N}v_{N}^{T}\\\
&\approx p(z,t),\end{split}$ (9)

reconstructs $p(z,t)$ with remarkable accuracy (owing to the orthonormality of
the $u_{n}$ and $v_{n}$). 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}$, $\sigma_{n}$, and $v_{n}$, are the
topo, singular value (energy), and chrono associated with the
$n^{\mathrm{th}}$ basis, which are ordered from highest contribution (energy)
to lowest. In Eq. 9, the multiplication of $u_{n}v_{n}^{T}$ is an outer
product. As a final note, we refer to $\sigma_{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 $\sigma_{n}^{2}$
as the energy.

###  4.2 Proper Orthogonal Decomposition (POD)

The POD method, also referred to as Biorthogonal Decomposition (BD), has been
used extensively for plasma physics datasets [27] but has only been very
recently introduced to Hall thruster and hollow cathode high-speed video [32,
33]. The POD algorithm is nearly identical to the SVD algorithm with the main
distinction being that the $U$, $\Sigma$, 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$, $\Sigma$, and $V^{T}$ matrices to identify and interpret the
dynamics associated with each mode. First, the energy for each basis,
$\sigma_{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.

Figure 10:  The SVD algorithm orders its bases from highest to lowest energy
(i.e. their total contribution to the reconstructed video), and therefore the
lower numbered bases are more likely to contain the most prominent dynamics.

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.

Figure 11:  The POD topos for the first 10 bases show $m=0$ or $m>0$ azimuthal
structures in the channel or around the cathode. Amplitude units are
arbitrary. By inspection, topos correspond to Fourier sine/cosine pairs as
shown in Table 1

Each of these topos has an associated time evolution, i.e. chronos, associated
with it. Several select chronos, $\sigma_{n}v^{T}_{n}$, 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.

Figure 12:  The POD chronos for the $n=1$, 2, 3, 5 and 6 bases. These bases
are roughly equivalent to the time evolutions of the $m=0$, 3, and 4 Fourier
modes 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.

Figure 13:  The Welch-averaged FFT of select POD chronos helps in identifying
the presence of a coherent mode structure and their characteristic
frequencies. This figure is directly comparable to the Fourier results in
Figure 6b.

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 [50].

Table 1: The POD bases and Fourier modes show a near one-to-one mapping. POD basis, n | Fourier mode, m | Location
---|---|---
1 | 0 | Channel
2, 3 | 3, sine and cosine | Channel
4 | 0 | Cathode
4, 6 | 4, sine and cosine | Channel
7, 8 | 1, sine and cosine | Cathode

###  4.3 SVD methods conclusion

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.

##  5 Conclusions

This work provided a survey of existing Fourier and new SVD-based mode
decomposition techniques for high-speed Hall thruster video data. To highlight
the implementation of these methods, they were each applied to the same H6
dataset which contained simultaneous modes associated with the channel and
cathode. Both Fourier and POD methods were able to characterize the spatial
and temporal evolution of each mode, including mode hopping between the $m=3$
and $m=4$ channel modes. Both methods also had their various strengths and
weaknesses.

Fourier methods are well suited for characterizing linear wave dynamics and
therefore excel at identifying both discrete mode and continuous wave dynamics
in the various regions of the Hall thruster plasma. A notable feature of the
Fourier methods are their ability to naturally provide dispersion
relationships. However, Fourier methods require preprocessing, e.g. converting
the video to polar coordinates, which adds computational complexity.

The main advantage of the SVD methods (notably POD) is that they produce a
tailored set of orthogonal bases for each dataset that are ordered by their
relative amplitudes (high to low), and these bases do not presume any form
(e.g. sines and cosines in Fourier). The POD method is relatively simple to
implement as preprocessing, most notably coordinate conversion, is optional.
In addition, the POD method can naturally locate mode structures at different
spatial locations (e.g. cathode and anode) unlike Fourier which requires the
dataset to be preprocessed with a particular spatial location in mind.

Together, Fourier and POD methods provide for a more complete toolkit for
studying global Hall thruster dynamics as captured by high-speed imaging.
Moving forward, we recommend first applying POD to a dataset in order to
identify its major dynamics and the spatial locations of each. With this
information, Fourier methods can then be intelligently targeted to extract
dispersion relationships or other select modal information. The code and data
used in this work are made are available online [39] to assist future
researchers.

## Acknowledgments

This analysis was performed while JWB held an NRC Research Associateship award
at the Naval Research Laboratory (NRL), with JWB and MM supported by the NRL
Base Program. The imaging data was collected at the University of Michigan by
MM under advisor Alec Gallimore in the Plasmadynamics and Electric Propulsion
Laboratory, where the use of the Photron SA5 FASTCAM was made possible via
grant FA9550-09-1-0695 from the Air Force of Scientific Research (AFOSR). AAK
acknowledges funding from the Army Research Office (ARO W911NF-19-1-0045) and
the Air Force Office of Scientific Research (AFOSR FA9550-18-1-0200).

## References

## References

  * [1] Goebel D M and Katz I 2008 Fundamentals of Electric Propulsion: Ion and Hall Thrusters (JPL Space Science And Technology Series no 1) (NASA JPL) URL https://www.wiley.com/en-us/Fundamentals+of+Electric+Propulsion%3A+Ion+and+Hall+Thrusters-p-9780470436264
  * [2] Boeuf J P 2017 Journal of Applied Physics 121 011101 ISSN 0021-8979 URL https://aip-scitation-org.ezproxy.cul.columbia.edu/doi/10.1063/1.4972269
  * [3] Florenz R, Liu T, Gallimore A, Kamhawi H, Brown D, Hofer R and Polk J Electric propulsion of a different class: The challenges of testing for MegaWatt missions 48th AIAA/ASME/SAE/ASEE Joint Propulsion Conference & Exhibit Joint Propulsion Conferences (AIAA) URL https://arc.aiaa.org/doi/10.2514/6.2012-3942
  * [4] Hall S, Jorns B, Gallimore A, Kamhawi H, Haag T, Mackey J, Gilland J, Peterson P and Baird M High-power performance of a 100-kW class nested hall thruster 35th International Electric Propulsion Conference URL https://pepl.engin.umich.edu/pdf/IEPC-2017-228.pdf
  * [5] Choueiri E Y 2001 Physics of Plasmas 8 1411–1426 ISSN 1070-664X URL https://aip.scitation.org/doi/abs/10.1063/1.1354644
  * [6] Parker J B, Raitses Y and Fisch N J 2010 Applied Physics Letters 97 091501 ISSN 0003-6951 URL https://doi.org/10.1063/1.3486164
  * [7] McDonald M S and Gallimore A D 2011 Measurement of Cross-Field Electron Current in a Hall Thruster Due to Rotating Spoke Instabilities 47th AIAA/ASME/SAE/ASEE Joint Propulsion Conference & Exhibit (AIAA 2011-5810) URL https://arc.aiaa.org/doi/abs/10.2514/6.2011-5810
  * [8] Jorns B A, Cusson S E, Brown Z and Dale E 2020 Physics of Plasmas 27 022311 ISSN 1070-664X URL https://aip.scitation.org/doi/full/10.1063/1.5130680
  * [9] Sekerak M J, Gallimore A D, Brown D L, Hofer R R and Polk J E 32 903–917 URL https://doi.org/10.2514/1.B35709
  * [10] Dale E Two-zone hall thruster breathing mode mechanism, part I: Theory 36th International Electric Propulsion Conference URL https://pepl.engin.umich.edu/pdf/IEPC-2019-354.pdf
  * [11] Goebel D M, Jameson K K, Katz I and Mikellides I G 14 103508 ISSN 1070-664X URL https://aip.scitation.org/doi/10.1063/1.2784460
  * [12] Georgin M P, Jorns B A and Gallimore A D 26 082308 ISSN 1070-664X URL https://aip.scitation.org/doi/10.1063/1.5111552
  * [13] Jorns B A and Hofer R R 2014 Physics of Plasmas 21 053512 ISSN 1070-664X URL https://doi.org/10.1063/1.4879819
  * [14] Janes G S and Lowder R S 1966 The Physics of Fluids 9 1115–1123 ISSN 0031-9171 URL https://aip.scitation.org/doi/10.1063/1.1761810
  * [15] McDonald M S and Gallimore A D 2011 IEEE Transactions on Plasma Science 39 2952–2953 ISSN 0093-3813 URL http://umich.edu/~peplweb/pdf/IEEE-TransPlasmaSci-2161343-2011.pdf
  * [16] Baird M 2020 Investigating Newly Discovered Oscillation Modes in Magnetically Shielded Hall Effect Thrusters Utilizing High Speed Diagnostics Ph.D. thesis Western Michigan University Kalamazoo, MI URL https://scholarworks.wmich.edu/dissertations/3656
  * [17] McDonald M S and Gallimore A D 2013 Comparison of Breathing and Spoke Mode Strength in the H6 Hall Thruster Using High Speed Imaging 33rd International Electric Propulsion Conference (Washington, DC) URL http://electricrocket.org/IEPC/28xa7w8t.pdf
  * [18] Darnon F, Lyszyk M, Bouchoule A, Darnon F, Lyszyk M and Bouchoule A 1997 Optical investigation on plasma investigations of SPT thrusters 33rd Joint Propulsion Conference and Exhibit (American Institute of Aeronautics and Astronautics) URL https://arc.aiaa.org/doi/abs/10.2514/6.1997-3051
  * [19] Hara K, Sekerak M J, Boyd I D and Gallimore A D 2014 Journal of Applied Physics 115 203304 ISSN 0021-8979 URL https://aip-scitation-org.ezproxy.cul.columbia.edu/doi/10.1063/1.4879896
  * [20] McDonald M and Gallimore A 2011 Parametric investigation of the rotating spoke instability in hall thrusters 32nd International Electric Propulsion Conference (Wiesbaden, Germany: IEPC 2011-242) URL http://www.umich.edu/~peplweb/pdf/IEPC-2011-242.pdf
  * [21] Romadanov I, Raitses Y and Smolyakov A 2019 Plasma Physics Reports 45 134–146 ISSN 1562-6938 URL https://doi.org/10.1134/S1063780X19020156
  * [22] Brunton S L and Kutz J N 2019 Data-Driven Science and Engineering: Machine Learning, Dynamical Systems, and Control 1st ed (Cambridge: Cambridge University Press) ISBN 978-1-108-42209-3 URL https://doi.org/10.1017/9781108380690
  * [23] Benner P, Gugercin S and Willcox K 2015 SIAM review 57 483–531 URL https://epubs.siam.org/doi/10.1137/130932715
  * [24] Taira K, Brunton S L, Dawson S, Rowley C W, Colonius T, McKeon B J, Schmidt O T, Gordeyev S, Theofilis V and Ukeiley L S 2017 AIAA Journal 55 4013–4041 URL https://arc.aiaa.org/doi/10.2514/1.J056060
  * [25] Noack B R, Afanasiev K, Morzyński M, Tadmor G and Thiele F 2003 Journal of Fluid Mechanics 497 335–363 URL https://doi.org/10.1017/S0022112003006694
  * [26] Oudheusden B, Scarano F, Van Hinsberg N and Watt D 2005 Experiments in Fluids 39 86–98 URL https://link.springer.com/article/10.1007/s00348-005-0985-5
  * [27] Dudok de Wit T, Pecquet A, Vallet J and Lima R 1994 Physics of Plasmas 1 3288–3300 ISSN 1070-664X URL https://aip.scitation.org/doi/10.1063/1.870481
  * [28] Levesque J P, Rath N, Shiraki D, Angelini S, Bialek J, Byrne P J, DeBono B A, Hughes P E, Mauel M E, Navratil G A, Peng Q, Rhodes D J and Stoafer C C 2013 Nuclear Fusion 53 073037 ISSN 0029-5515 URL https://iopscience.iop.org/article/10.1088/0029-5515/53/7/073037
  * [29] van Milligen B P, Sánchez E, Alonso A, Pedrosa M A, Hidalgo C, de Aguilera A M and Fraguas A L 2014 Plasma Physics and Controlled Fusion 57 025005 URL https://iopscience.iop.org/article/10.1088/0741-3335/57/2/025005/meta
  * [30] Hansen C, Victor B, Morgan K, Jarboe T, Hossack A, Marklin G, Nelson B and Sutherland D 2015 Physics of Plasmas 22 056105 URL https://doi.org/10.1063/1.4919277
  * [31] Kaptanoglu A A, Morgan K D, Hansen C J and Brunton S L 2021 Phys. Rev. E 104(1) 015206 URL https://link.aps.org/doi/10.1103/PhysRevE.104.015206
  * [32] Désangles V, Shcherbanev S, Charoy T, Clément N, Deltel C, Richard P, Vincent S, Chabert P and Bourdon A 2020 Atmosphere 11 518 URL https://www.mdpi.com/2073-4433/11/5/518
  * [33] Becatti G, Goebel D M and Zuin M 129 033304 ISSN 0021-8979, 1089-7550 URL http://aip.scitation.org/doi/10.1063/5.0028566
  * [34] Schmid P J 2010 Journal of Fluid Mechanics 656 5–28 URL https://doi.org/10.1017/S0022112010001217
  * [35] Tu J H, Rowley C W, Luchtenburg D M, Brunton S L and Kutz J N 2014 Journal of Computational Dynamics 1 391–421 URL https://www.aimsciences.org/article/doi/10.3934/jcd.2014.1.391
  * [36] Taylor R, Kutz J N, Morgan K and Nelson B A 2018 Review of Scientific Instruments 89 053501 URL https://doi.org/10.1063/1.5027419
  * [37] Sasaki M, Kawachi Y, Dendy R O, Arakawa H, Kasuya N, Kin F, Yamasaki K and Inagaki S 2019 Plasma Physics and Controlled Fusion 61 112001 URL https://iopscience.iop.org/article/10.1088/1361-6587/ab471b
  * [38] Kaptanoglu A A, Morgan K D, Hansen C J and Brunton S L 2020 Physics of Plasmas 27 032108 URL https://aip.scitation.org/doi/10.1063/1.5138932
  * [39] Brooks J 2021 Code and data for a comparison of fourier and pod mode decomposition methods for high-speed hall thruster video URL http://dx.doi.org/10.5281/zenodo.5150716
  * [40] Brown D L 2009 Investigation of Low Discharge Voltage Hall Thruster Characteristics and Evaluation of Loss Mechanisms. Ph.D. thesis University of Michigan URL https://deepblue.lib.umich.edu/handle/2027.42/63660
  * [41] Kåsa I 1976 IEEE Transactions on Instrumentation and Measurement IM-25 8–14 ISSN 1557-9662 URL https://ieeexplore.ieee.org/document/6312298
  * [42] Taubin G 1991 IEEE Transactions on Pattern Analysis and Machine Intelligence 13 ISSN 1939-3539 URL https://ieeexplore.ieee.org/document/103273
  * [43] Vora P L, Farrell J E, Tietz J D and Brainard D H 1997 Linear models for digital cameras Proceedings, IS&T’s 50th Annual Conference pp 377–382 URL http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.53.7490
  * [44] McDonald M, Bellant C, St Pierre B and Gallimore A 2011 Measurement of cross-field electron current in a hall thruster due to rotating spoke instabilities 47th AIAA Joint Prop. Conf. URL https://arc.aiaa.org/doi/abs/10.2514/6.2011-5810
  * [45] Welch P 1967 IEEE Transactions on Audio and Electroacoustics 15 ISSN 1558-2582 URL https://doi.org/10.1109/TAU.1967.1161901
  * [46] Yamada T, Itoh S I, Inagaki S, Nagashima Y, Shinohara S, Kasuya N, Terasaka K, Kamataki K, Arakawa H, Yagi M, Fujisawa A and Itoh K 2010 Physics of Plasmas 17 ISSN 1070-664X URL https://aip.scitation.org/doi/10.1063/1.3429674
  * [47] Golub G H and Reinsch C 1970 Numerical Mathematics 14 403–420 URL https://link.springer.com/article/10.1007/BF02163027
  * [48] Golub G H and Loan C F V 1996 Matrix Computations 3rd ed (Baltimore: Johns Hopkins University Press) ISBN 978-0-8018-5414-9 URL https://jhupbooks.press.jhu.edu/title/matrix-computations
  * [49] Woolfe F, Liberty E, Rokhlin V and Tygert M 2008 Applied and Computational Harmonic Analysis 25 335–366 URL https://doi.org/10.1016/j.acha.2007.12.002
  * [50] Holmes P, Lumley J L, Berkooz G and Rowley C W 2012 Turbulence, coherent structures, dynamical systems and symmetry (Cambridge university press) URL https://doi.org/10.1017/CBO9780511919701