Searching for Chaotic Behavior in the Ion Current Waveforms of a Hall Effect Thruster

Abstract

The dynamics of the ion flux expelled by a 0.5 kW-class Hall thruster supplied with krypton was examined in a wide range of discharge voltages. A homemade Faraday probe installed onto a rotary arm was used for reconstructing angular profiles of the plasma plume 0.5 m downstream of the thruster exit plane. The time dependence of the ion current was measured along the thruster axis. For investigating the signal dynamics, a Fourier approach as well as methods of nonlinear time series analysis like bifurcation diagrams and recurrence plot techniques were applied, which are of interest for chaotic behavior identification. Along with the well-known breathing mode (10—30 kHz), other characteristic groups of oscillations were also detected. The bifurcation diagram revealed a drastic transition between large and small amplitude oscillating regimes while varying the discharge voltage from 550 to 700 V. In parallel to this transition, recurrent plots display a qualitative change from a periodic (or quasi periodic) oscillating regime to much less predictable dynamics.

Introduction

Laboratory plasmas and associated diagnostic systems have a broad range of applications, from small-scale plasma propulsion devices such as Hall effect thrusters to large thermonuclear reactors like ITER or DEMO [1, 2]. The initial goal of this work was to investigate the possible emergence of chaotic behavior in a Hall effect thruster, using ion current time series. The plasma analysis methods for this purpose can also be applied to high-temperature plasma devices. For example, the magnetic topology at the tokamak plasma edge can exhibit some chaotic behavior [3]. More recently, recurrence plots has been employed at the Joint European Torus (JET) to study the non-linear dynamics of Edge Localized Modes (ELMs) [4], which are tokamak instabilities characterized by a periodic crash of the pedestal temperature and which show some similarities with the breathing mode observed in Hall effect thrusters.

In this paper, the dynamics of the ion flux expelled by a 0.5 kW-class Hall thruster supplied with krypton was examined in a wide range of discharge voltages. In Sect. 2, the experimental setup is introduced. A homemade Faraday probe installed onto a rotary arm was used for reconstructing angular profiles of the plasma plume 0.5 m downstream of the thruster exit plane. In section III, the ion current time dynamics is analysed in details, using a Fourier approach as well as methods of nonlinear time series analysis like bifurcation diagrams and recurrence plot techniques, which are of interest for chaotic behavior identification. Finally, in Section IV some conclusions and perspectives are given, based on the observations and analysis of the ion current dynamics.

Experimental Setup

Hall Effect Thruster

A High Impulse Krypton Hall Effect Thruster (HIKHET) targeted for operation with krypton was designed from scratch in the Institute of Plasma Physics and Laser Microfusion (IPPLM) in Warsaw. With an average diameter of the discharge channel of 42 mm, this thruster belongs to the 0.5 kW-class of engines. In the present study, it was used for ion current dynamics examination. All the details about the HIKHET’s construction were already presented in our previous papers [5,6,7]. Therefore, to avoid redundancy, only a picture of the thruster installed in the IPPLM’s vacuum chamber (Fig. 1) and several photos of the plasma plume as taken during the thruster operation are reproduced here.

Fig. 1

HIKHET prototype dedicated to Kr operation

Ion Diagnostic Probe

During the project funded by ESA (Contract No. 4000122415/17/NL/GE), three electric probes have been designed in IPPLM for investigating the angular profile and temporal behavior of the ion beam generated by the thruster: two Faraday cups and a planar probe with a guarding ring, also known as Faraday probe [8,9,10]. For the needs of the present study, only ion signals recorded with the Faraday probe (FP) were used. Nevertheless, the consistency of the ion waveforms and profiles as captured by the FP and the other probes had been checked before (see Fig. 2). The wings in the ion current profiles around ± 60° are due to charge-exchange collisions at the edge of the plasma plume, increasing slightly the registered ion current at these positions.

Fig. 2

Angular profile of the ion current density collected with the thruster operating at a discharge voltage of 600 V and a coil current ratio I_inn / I_out = 7.0 A / 4.1 A (the anode mass flow rate was left at 0.8 mg/s with 0.17 mg/s from the cathode)

The FP collector and the ring (Fig. 3, right) were made of graphite (to reduce the secondary electron emission with respect to typical metals) and separated by a distance of 0.5 mm. Both parts were polarized at the same potential. The shield of the probe, also made of graphite, was left at the floating potential [9].

Fig. 3

Three electric probes mounted on the diagnostic rotary arm

Experimental Facility

Throughout the experimental session, the thruster was installed inside a cylindrical ~ 2 m³ IPPLM’s vacuum chamber equipped with a three-stage pumping system, resulting in an average pumping speed of the facility close to 18 m³/s for krypton (and 14 m³/s for xenon). Because measurements were performed while keeping a mass flow rate of 0.8 mg/s, the ultimate pressure inside the vacuum chamber varied between 1.8 and 1.9 × 10^–5 mbar (when the gas flow was switched off, the background pressure dropped below 5.0 × 10^–8 mbar). The pressure was measured by a Ionivac Oerlikon-Leybold ITR 90 vacuum gauge corrected for krypton.

For the purpose of performance test, the thruster was anchored to a thrust balance (manufactured by the MECARTEX Swiss company) that allows for thrust measurement with 2% accuracy. A dedicated rotary system was installed to drive the plasma diagnostic probes inside the chamber [7]. During the experimental session, the FP was installed on an aluminium arm at a distance of 0.494 m from the thruster exit surface. In order to repel electrons and collect only ions from the plasma beam, the FP collector and the ring were polarized at the voltage of -15 V with respect to the laboratory ground. The collector was biased by a homemade power unit regulated by adapting the number of alkaline batteries, while the ring was biased by a homemade electronically regulated module. The ion current was recorded using a 12 bit GaGe data acquisition card with the sampling frequency set to 50 MHz and the acquisition time to 20 ms.

Discharge Current Oscillations

Although HET is also called stationary plasma thruster and operates steadily (almost) as long as power and propellant are supplied, the plasma instabilities induce a modulation of the discharge current in a very broad frequency range. The most common type of oscillations, so-called breathing mode (BM), appears within the range of 10 – 60 kHz. The mechanism of these large amplitude low frequency volumetric oscillations has been investigated with simple analytical models of prey-predator type for the neutral gas ionization (0D model) [11,12,13,14], as well as with more sophisticated theories [15,16,17]. The BM appears in numerical simulations even within 1D description [16, 18, 19] and are unambiguously reproduced by 2D models, e. g. [20]. Over time, kinetic [21], hydrodynamic [22] and hybrid [20] models have been developed.

It is supposed that electron transport can be the key contributor to the stabilization of discharge oscillations [19]. So far, many articles have been written on HET discharge current oscillations. Without going deeper into details, it is worth emphasizing that transitions between global (including BM) and local modes are also commonly observed [19, 23]. It was reported that these transitions can cause an increase of the discharge current mean value or/and of the oscillation amplitude, that can be deleterious for the thruster performance [24]. The transitions may be triggered by the intentional variation in such parameters like B-field magnitude, mass flow rate and discharge voltage and thus, to some extent, controlled, as well as by the uncontrolled change of these parameters e.g. due to magnetic circuit temperature growth.

Thruster Performance

The main characteristics of the thruster performance as a function of the applied discharge voltage ${U}_{d}$ are depicted in Fig. 4. The anode efficiency $\eta$ has a minimum at 550 V and reaches its maximal values for the highest voltages. The current utilization, i.e. the ion current to discharge current ratio ${I}_{i}/{I}_{d}$ also displays a minimum at 550 V, but interestingly gives the same maximal value at both ends of the examined voltage range. The thrust $T$ (and so the specific impulse ${I}_{sp}$) as well as the discharge power ${P}_{d}$ are systematically increasing with the discharge voltage. The ion beam divergence $\alpha$ remains almost constant with respect to ${U}_{d}$. The frequency $f$ of the main BM peak, as it will be shown in the power spectral density (PSD) plots, has a local minimum at 500 V. The unambiguous drop in thruster performance around 550 V precedes a drastic transition in the ion current dynamics, that is illustrated in the next section with the bifurcation diagrams, power spectra and recurrence plots.

Fig. 4

Normalized values of the thrust T (specific impulse I_sp), anode efficiency η, current utilization ${\mathrm{I}}_{\mathrm{i}}/{\mathrm{I}}_{\mathrm{d}}$, ion beam divergence α, frequency of the main BM peak f, and discharge power P_d as a function of the applied discharge voltage U_d

Ion Current Data Analysis

The initial goal of this work was to investigate the possible emergence of chaotic behavior using ion current time series. The signals were collected only in a thermally stable state of the thruster. Keeping the main operating parameters unchanged, i.e. a krypton mass flow rate of 0.8 mg/s and a coils current ratio of I_inn / I_out = 7.0 A / 4.0 A (leading to a magnetic field of up to ~ 500 G in the thruster channel), the discharge voltage was varied in the range 300–700 V. Other coils current ratios are not considered in this paper and have been reported in other works [25, 26]. A comparison of thruster performance with krypton and xenon can be found in [27, 28].

The first attempts to classify the discharge current–voltage (I–V) characteristics by regimes relevant to specific features of the oscillations date back to the 70 s, where Tilinin selected six main patterns [29]. Gascon et al., in the search for dynamic characteristics, refers to that in [30] and [31].

To reconstruct basic characteristics (like invariants) of dynamic systems, it is sufficient to capture time series representing a separate observable and apply the Takens embedding theorem [32, 33]. Even though the discharge current could be considered to be such an observable as it is strongly correlated with ion local oscillations in the plasma [34], it also includes an electron component. Therefore, we decided to use the ion current density as measured in a far-field of the thruster (0.5 m downstream of the thruster exit plane). Despite the fact that the ion current density measured in the near-field of the thruster by a small Faraday probe has already been already with the global discharge current [34], in the present study a systematic examination of the ion current density fluctuations was performed in a wide range of discharge voltages.

Bifurcation Diagram

The changes in the ion current behavior can be represented by the frequency of occurrence of signal instantaneous amplitudes. A histogram of local extrema, so-called bifurcation diagram, that displays the evolution of the ion current amplitudes as captured by the FP for different discharge voltages is presented in Fig. 5(a). Blue dots represent local maxima, while red ones denote local minima of the ion current. A smoothing procedure based on a moving average method was applied to the current waveforms, using a time window of 0.2 μs. The window size of the data set used to search for minima and maxima was set to 100 μs.

Fig. 5

Bifurcation diagram with the discharge voltage as a control parameter (a), density plot (b) and pictures corresponding to the changing voltage (c)

The standard bifurcation diagram can be enriched with information about population of extrema of a given value identified within relevant time series, see Fig. 5(b). It may provide important information about how often the system visits its specific states. In this analysis, the voltage resolution was limited to 50 V and the ion current resolution was set to 0.035 mA. Additionally, photos showing the behavior of the plasma beam for different applied voltages are presented in Fig. 5(c).

It can be concluded from Fig. 5 that the amplitude of the ion current oscillations gradually increases with the growing voltage from 300 to 550 V. Usually, bifurcation diagrams indicate values of a control parameter (in our case the discharge voltage) for which dramatic changes of a displayed function take place. Here, such a catastrophic transition in the ion current character corresponds to 600 V. After some threshold voltage between 550 and 600 V, the ion current amplitude collapses and does not show an upward trend when increasing the voltage further. A more complex behavior is observed at 500 V, for which the instantaneous amplitudes of the oscillating current gather themselves into several separate groups of different densities. This ion current behavior is reflected in Fig. 5(c), where the shape of the beam continuously expands up to 550 V, while for higher voltages the plasma beam collapses, the plasma seems to change color and the broadening of the beam with increasing voltage is no longer visible.

Power Spectrum

While the BM have usually been studied by looking only at the discharge current, in this work pure oscillations of the ion current, undisturbed by the electron current, were examined. The Welch method was used to determine the power spectral density (PSD) of the ion current time series (normalized to zero mean value and unit variance), splitting up 2 ms long signals into ten segments with 50% overlap. The PSD for all voltages from 300 to 700 V are shown in Figs. 6, 7, and 8 and summarized in the 2D diagram shown in Fig. 9. The main peak visible at a frequency of around 30 kHz (and its harmonics) in each spectrum corresponds to BM.

Fig. 6

PSD of the ion current waveform for discharge voltages of 300 V (red), 350 V (blue) and 400 V (green)

Fig. 7

PSD of the ion current waveform for discharge voltages of 450 V (red), 500 V (blue) and 550 V (green)

Fig. 8

PSD of the ion current waveform for discharge voltages of 600 V (red), 650 V (blue) and 700 V (green) (Color figure online)

Fig. 9

PSD as a function of the discharge voltage

For 300 V the BM peak is very fuzzy, perhaps due to a significant noise distortion, and the widest among the recorded voltages below 600 V. The PSD is also poor in harmonics. Between 300 and 500 V, the number and intensity of harmonics increases with the voltage. For 500 V, additional peaks around the main BM peak and its harmonics are clearly visible, that may suggest a quasi-periodic process. A dramatic change occurs in the PSD between 550 and 600 V, with the collapse of the BM harmonics, as clearly visible in Fig. 9.

For discharge voltages of 600 V and higher, the BM peak becomes extremely broad while its harmonics magnitude appears very low or inseparable. On the other hand, the intensity of the peak at around 200 – 300 kHz grows to a value that is comparable to the one of the BM peak. Due to its temporal scale which is commensurable with the time-of-flight of ions through the acceleration zone, the relevant oscillations were categorized as transit time oscillations [24]. Interestingly, the intensity of this peak increases for the highest as well as for the lowest values of the discharge voltage, while it disappears in the range 450–550 V.

Recurrence Plot

To gain more insight into the ion current dynamics, one of the main graphical tools for time series analysis, so-called recurrence plot (RP), was used. The RPs are a simply achievable way for the diagnosis of dynamical systems. They display important information about time scales which are not easily obtainable by other methods. Often, characteristic patterns which appear in the RPs can be used to classify dynamic processes [35]. For example, the simple pattern that distinguishes periodic systems is a set of continuous lines parallel to the diagonal. However, in the case of noise, the RP displays randomly placed points without any specific structure, whereas specific patterns appear for systems with more complex dynamics (including chaotic systems) [36].

In the analysis of chaotic systems, the Takens embedding theorem implies that the information about an attractor can be unfolded from only one time series $x\left(t\right)$. For that purpose, a time delay $\tau$ is used to create a multivariate vector in a d-dimensional phase space $y\left(t\right)=\left(x\left(t\right), x\left(t+\tau \right), \dots , x\left(t+(d-1)\tau \right)\right)$. The integer d-value is called the embedding dimension (ED). Following the algorithm described in [38], a method based on the autocorrelation function was used for a reasonable time-delay selection. Consequently, the time delay corresponding to the first zero of the autocorrelation function was chosen, showing how fast the information between $x\left(t\right)$ and $x\left(t+\tau \right)$ is lost.

The minimum ED that allows reconstructing characteristics of the attractor is determined by using a method based on the identification of false nearest neighbors (FNN). This method is described in details in [37] and is briefly summarized in Appendix for clarity.

The obtained recurrence plots are presented in Figs. 10, 11 and 12. The condition for identifying two points $y\left({t}_{i}\right)$, $y\left({t}_{j}\right)$ of the phase space as neighbors (marked by a black dot in the RP) was chosen as:

$${R}_{d}\left(i,j\right)<0.2 L$$

(1)

where $L$ is the average distance between all points visited in the phase space. The value of 0.2 was chosen empirically, in order to obtain a sufficient statistics in the RP.

Fig. 10

Recurrence plots for three sets of voltages: 300 V, 350 V and 400 V. $\uptau$ indicates the time delay, d the embedding dimension and R the radius for which the nearest points in the multidimensional space are located

Figure 10 shows the obtained recurrence plots for the three lowest voltages. When analyzing the RP for 300 V, the patterns seem to indicate a periodic process with quite regular disturbances, as evidenced by the appearance of lighter and darker areas. In the case of 350 V, it can be seen that the periodicity is frequently interrupted. It is interesting to note that this loss of stability creates a regular map. The white areas (areas where there is no predictability) are fairly well-defined on the horizontal and vertical bands. On the other hand, in the case of 400 V, these white areas are moved towards the graph borders. Thus, the loss of regular behavior occurs when the signal is shifted by very large but specific delays. This reveals again some regularity.

Figure 11 shows the RP for a medium range of voltages. The patterns in the RP for 450 V are similar to the ones for 400 V, except that instead of a white frame, there are white and black stripes with a distinct white band in the center of the graph. For 500 V, laminar diagonal patterns are broken with vertical and horizontal stripes, displaying a very regular loss of predictability. This may indicate a certain irregularity of these regular processes resulting from chaotic behavior. It is worth noting that in the power spectrum of this signal, besides the leading frequency of 27.6 kHZ (and its harmonics) which represents BM oscillations, there is a peak of about 13.2 kHz (with harmonics) that seems to be incommensurable with the BM one (27.6 / 13.2 = 2.09…). In this sense the signal could be considered as quasi-periodic. The weaken diagonal patterns between stronger ones are likely representative for this lower frequency which, by its value, could be depictive for spoke oscillations.

Fig. 11

Recurrence plots for three sets of voltages: 450 V, 500 V and 550 V. $\uptau$ indicates time delay, d embedding dimension and R radius for which the nearest points in the multidimensional space are located

For 550 V, the pattern indicates almost perfect harmonic oscillations. Parallel lines are clearly visible in relation to the diagonal, proving the periodicity of the process.

The last three RPs are depicted in Fig. 12. At 600 V, the pattern is characterized by noise where some regular behavior appears sporadically. At 650 V, frequent irregular patterns reappear. At 700 V, the structure of the lattice is very visible.

Fig. 12

Recurrence plots for three sets of voltages: 600 V, 650 V and 700 V. $\uptau$ indicates the time delay, d the embedding dimension and R the radius for which the nearest points in the multidimensional space are located

A recurrence quantification analysis (RQA) is performed based on the RP, using three typical quantities: the determinism (DET), the laminarity (LAM) and the Shannon entropy (ENTR). The DET corresponds to the percentage of RP points that form diagonal lines:

$$DET=\frac{\sum_{l={l}_{min}}^{N}l\cdot P\left(l\right)}{\sum_{l=1}^{N}l\cdot P\left(l\right)}$$

(2)

where $P\left(l\right)$ is the histogram of the lengths $l$ of the diagonal lines and ${l}_{min}$ is a threshold value to consider a line (${l}_{min}=2$ here).

The LAM is the percentage of recurrence points which forms vertical lines:

$$LAM=\frac{\sum_{v={v}_{min}}^{N}v\cdot P\left(v\right)}{\sum_{v=1}^{N}v\cdot P\left(v\right)}$$

(3)

where $P\left(v\right)$ is the histogram of the lengths $v$ of the vertical lines and ${v}_{min}$ is a threshold value to consider a line (${v}_{min}=2$ here).

The ENTR is the Shannon entropy of the probability distribution of the diagonal lines lengths $p\left(l\right)$:

$$ENTR=-\sum_{l={l}_{min}}^{N}p\left(l\right)\cdot \mathrm{ln}p\left(l\right)$$

(4)

These quantities are presented on Fig. 13 as a function of the discharge voltage. Three different regions can be observed: a first zone $300 V\le {U}_{d}\le 400 V$ of increase of DET, LAM and ENTR with ${\mathrm{U}}_{\mathrm{d}}$, followed by a saturation zone for $450 V\le {U}_{d}\le 550 V$ and finally a collapse for $600 V\le {U}_{d}\le 700 V$. It can be noted that the dynamics is strongly correlated with the intensity of the BM harmonics, as observed in the PSD in Fig. 9.

Fig. 13

Recurrence Quantification Analysis (RQA) of the RP showing the determinism (DET), the laminarity (LAM) and the Shannon entropy (ENTR) as a function of the discharge voltage ${\mathrm{U}}_{\mathrm{d}}$

Conclusion and Future Work

In this work, we investigated the dynamics of ion current waveforms from HIKHET as collected by a Faraday probe, for discharge voltages going from 300 to 700 V. Apart from the analysis of the power spectra, bifurcation diagrams and recurrence plots with autocorrelation functions and false nearest neighbors statistics were also examined. The PSD and bifurcation diagram showed an increase of the BM oscillations amplitude (~ 20-30 kHz) and of the harmonics intensity with the discharge voltage up to 550 V, after which a sudden collapse of the ion current amplitude is observed, while a less intense but higher frequency oscillations (~ 200–300 kHz) compatible with the ion transit time is present for the highest voltages. In this experiment, the thruster performance are affected by the BM intensity, as the anode efficiency and current utilization decrease with increasing voltage up to 550 V, while they increase sharply (as well as the specific impulse) at higher voltages. The RPs and RQA allowed to observe three different voltage regions: a first zone $300 V\le {U}_{d}\le 400 V$ associated with an increase of the determinism, laminarity and Shannon entropy with increasing voltage, followed by a saturation zone for $450 V\le {U}_{d}\le 550 V$ and finally a collapse of these three quantities for$600 V\le {U}_{d}\le 700 V$. The observed dynamics is strongly correlated with the intensity of the BM harmonics, and a detailed observation of some RPs unveils patterns indicating regular loss of predictability of the system (e.g. at${U}_{d}=500 V$), which could be an indicator of chaotic behavior. The study of chaos indicators in the ion current waveforms is an approach that could lead to a deeper understanding of relevant physical processes. Controlling discharge current instabilities (e.g. with chaos control methods) will be crucial to improve the performances of Hall thrusters in the future.

In the present study, the observations made in the bifurcation diagram and recurrence plots seem compatible with the presence of a chaotic behavior, although a strong evidence validating the hypothesis of deterministic chaos in the ion current has yet to be established. That is why further work will include other methods of nonlinear time series analysis e.g. to estimate the Lyapunov exponents of the system. However, it will be necessary to clean up the signals from noise with dedicated methods relying on notions of deterministic chaos theory. On this aspect, the main issue is to avoid distortion of the nonlinear dynamics during the noise reduction process.

Appendix

The main idea of searching for FNN is that, when the chosen embedding dimension is smaller than the actual dimension of the attractor, i.e. if the attractor is tried to be reconstructed in too small phase space, let say d-dimensional, some of the points that appear very close to each other (due to projection from higher dimension) could be separated by a much greater distance in d + 1 dimensions. Such points are called false neighbors. The distance between two points $y\left({t}_{i}\right)$ and $y\left({t}_{j}\right)$ in d dimensions ${R}_{d}\left(i ,j\right)$ is defined using the Euclidian norm:

$${R}_{d+1}^{2}\left(i,j\right)={R}_{d}^{2}\left(i,j\right)+{\left(x\left({t}_{i}+d\tau \right)-x\left({t}_{j}+d\tau \right)\right)}^{2}$$

(5)

The goal is therefore to find the smallest embedding dimension that minimizes the percentage of FNNs. In practice, 100 nearest neighbors for each point at the reconstructed trajectory in d-dimensional phase space are checked. Two conditions are used to identify FNNs when moving from the phase space in d dimensions to d + 1 dimensions. The first one assumes that the distance should not increase by a ratio greater than the given threshold value ${R}_{tol}=15$, such that:

$$\frac{{R}_{d+1}^{2}\left(i,j\right)-{R}_{d}^{2}\left(i,j\right)}{{R}_{d}^{2}\left(i,j\right)}>{R}_{tol}^{2}$$

(6)

However, this condition only allows to identify false neighbors that are very close to each other in d dimensions, but can fail with a limited dataset subject to perturbative noise. In this case, the distance between the nearest neighbors does not always satisfies ${R}_{d}\ll {R}_{A}$, with ${R}_{A}$ the size of the attractor. For points with ${R}_{d} \simeq {R}_{A}$,the second condition is:

$${R}_{d+1}\left(i,j\right)>{A}_{tol}{R}_{A}$$

(7)

where ${A}_{tol}=2$ is used hereafter. The plots of the autocorrelation functions and the FNN percentages versus embedding dimension, for all the examined discharge voltages, can be found in Appendix in Figs. 14 and 15 respectively. From Figs. 15, the embedding dimension d = 4 seems to be appropriate for studying the dynamic of the system [36].

Fig. 14

Autocorrelation functions applied to the ion current signals collected by the FP

Fig. 15

Percentage of FNN vs. ED, applied to the ion current signals collected by the FP