Simulations of driven breathing modes of a magnetically shielded Hall thruster

The operation of a 5 kW-class magnetically shielded Hall effect thruster with sinusoidal modulation of the discharge voltage is investigated through simulations with a 2D axisymmetric hybrid (particle-in-cell/fluid) code. The dynamic response of the thruster for different modulation amplitudes and frequencies is presented and discussed. The analysis of partial efficiencies contributing to thrust efficiency allows identifying counteracting effects limiting net gains in performance figures. Voltage modulation enhances the amplitude of plasma oscillations and can effectively control their frequency when the modulation frequency is close to that of the natural breathing mode (BM) of the thruster. The 2D plasma solution reveals that the dynamics of the ionization cycle are governed by the electron temperature response, enabling a driven BM at the modulation frequency. For modulation frequencies far from the natural BM one, voltage modulation fails to control the plasma production via the electron temperature, and the natural BM of the thruster is recovered. High order dynamic mode decomposition applied to the 2D plasma solution permits analyzing the complex spatio-temporal behavior of the plasma discharge oscillations, revealing the main characteristics of natural and externally driven modes.


Introduction
The Hall effect thruster (HET) is a mature electric propulsion technology, with relatively high efficiency and thrust. However, in spite of its successful flight history, there currently exist significant efforts to increase HET's efficiency and lifetime in response to new, on-going, ambitious space missions, including the deployment of near-Earth satellite megaconstellations for fast and far-reaching telecommunications, * Author to whom any correspondence should be addressed.
Original Content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. or the colonization of the Moon and Mars. Magnetic shielding (MS) of the channel walls [1] is an effective technique that is enabling the design of the next HET generation, featuring enhanced performances and operational lifetimes compared to classical HET designs such as the Stationary Plasma Thruster (SPT) [2,3], which relies on a magnetic lens type topology inside the channel.
The plasma discharge in both conventional and MS-HETs is characterized by the presence of the so-called breathing mode (BM), a low frequency (of the order of tens of kHz) discharge oscillation, which was first identified in the 1970s [4], and that has been typically linked to a predator-prey type ionization instability, involving an axial displacement of the ionization front inside the thruster chamber [5][6][7]. While the underlying physical mechanisms behind the onset and growth of the BM are still topics of present debate [8][9][10][11][12], several numerical and experimental studies have shown that, depending on the operating parameters, the BM can yield significant oscillations in the discharge current, I d , with amplitudes of the order of its DC value, which can greatly affect the thruster performance [13][14][15][16][17][18][19][20]. Moreover, the external anode-to-cathode circuit is known to have a significant impact on the discharge oscillations. Barral and Miedzik [21], analyzed resistor-inductorcapacitor (RLC) networks connecting anode and cathode and proportional-integral-derivative closed-loop control of the discharge voltage V d or of the applied magnetic field intensity as active means for achieving nearly oscillation-free operation. Wei et al [22] have reported experimental measurements showing a noticeable effect of the amplitude of I d oscillation on the thruster performance by changing only the capacitor value of an RLC filter unit connected between anode and cathode. Therefore, the assessment of the impact of the BM on the discharge performance continues to be of central importance in HET research.
Recent works have explored the so-called externally driven BM operation of the thruster by applying some external modulation to V d . These studies are motivated by different reasons. First, adding a sinusoidal voltage component to the anode potential with a frequency close to that of the natural (i.e. unmodulated) BM of the thruster has been proven effective to obtain time-coherent breathing oscillations, which have permitted the successful validation of time-resolved laser induced fluorescence diagnostic techniques with application to HET discharges [23][24][25][26][27][28].
Second, a novel power processing unit (PPU) concept for HETs relying on a pulsating drive of the anode voltage through a pulsating boost chopper circuit demonstrated a good performance and reduced size and weight compared to conventional PPUs for a 1 kW-class HET, thus making it attractive for high power applications [29]. A recent study has proven the pulsating boost chopper power supply (chopper PS) effective for controlling the frequency of the driven discharge oscillations for the case of a 1 kW-class anode-layer-type HET, which can greatly simplify electromagnetic interference (EMI) mitigation strategies on modern spacecrafts [30].
Moreover, several numerical and experimental studies of a modulated HET discharge have reported enhanced performance. Tamida et al [31] applied a chopper PS on a conventional HET prototype with a frequency close to that of the natural BM and with a modulation amplitude ranging between 10% and 50% of the average value of V d . Operating at 200 W, measurements indicate up to a 30% relative increase in thrust efficiency with respect to the unmodulated operation, the latter featuring a thrust efficiency of about 20%. Yamamoto et al [29] reported 23% and 16% relative increases on the thrust efficiency and the thrust-to-power ratio, respectively, with respect to the unmodulated discharge at V d = 150 V (featuring values of 22% and 50 mN kW −1 ), for a 1 kW-class HET operated with a chopper PS, when V d is modulated at a frequency close to that of the natural BM. The enhanced thrust performance is attributed to the synchronization between the predatorprey oscillation mechanism and the PPU, which optimizes the plasma production and the ion acceleration processes. Measurements in the plasma plume of the same thruster reported a negligible effect of the modulation on the plume divergence and a wider ion energy distribution function in pulsating operation [32].
Romadanov et al [33] studied the discharge and plasma plume characteristics of a 200 W cylindrical Hall effect thruster (CHT), operating with a sinusoidal V d modulation, through experimental measurements and numerical simulations using a time-dependent 1D axial fluid model of the discharge including neutrals, electrons and singly-charged ions. The time-averaged voltage was set to 220 V, and the peak-topeak modulation amplitude was varied from 8 to 50 V, while the modulation frequency was varied from 6 to 18 kHz, close to the 13 kHz natural BM frequency of the thruster. This type of modulation was proven effective to externally drive the BM resulting in enhanced I d oscillations. The amplitude of the I d oscillations and its root mean square (RMS) value were shown to depend on the modulation amplitude. Two distinct, linear and non-linear, response regimes were identified. Moreover, a resonant behavior of the discharge was observed, with maximum I d oscillation amplitude and RMS value at modulation frequencies close to that of the natural BM. The resonant frequency was also shown to depend on the modulation amplitude, ranging from 10 to 13 kHz. RMS values of the propellant utilization and current efficiencies are found to increase in modulated operation. Close to the resonant frequency, and at large modulation amplitude (half-peak values around 9%-11% of the time-averaged value) the authors report a total increase of up to a 20% in the product of the RMS values of the propellant utilization and current efficiencies compared to unmodulated operation. Enhanced oscillation and phase alignment of the ion velocity and ion density were postulated as one of the main reasons behind the observed behavior. The beam divergence is found almost unaffected by the external modulations. An experimental study [34] carried out with the same thruster showed that the rotating spoke instability, which could contribute to electron turbulent transport across the magnetic field, can be suppressed in the non-linear regime in externally driven BM operation, thus suggesting that voltage modulation can increase current efficiency.
A theoretical analysis demonstrated that for modulated operation, an increase of thrust can be obtained if the ion beam current and energy oscillate in-phase [35]. In [36], timedependent measurements of the ion current and ion energy in the plume of the CHTpm2, a permanent magnet version with magnetic shield of the low power (100-200 W) CHT mentioned above, operating with sinusoidal V d modulation, revealed a decrease in the phase angle between these two signals when the modulation frequency was close to that of the natural BM, which is optimal for high thrust production. However, recent experiments and 1D axial hybrid numerical simulations of the modulated CHTpm2 have reported lower than expected thrust gains due to increased plume divergence and lower ion energy amplitude close to the resonant frequency of the modulation [37]. Moreover, thrust gains are counteracted by gains in discharge power, thus the anticipated net increase in thrust efficiency is not achieved.
Therefore, further analysis of the complex spatio-temporal behavior of main plasma magnitudes inside the thruster chamber is required to advance in the understanding of the modulated HET response. In particular, the dynamics of the ionization and acceleration regions, which are highly coupled in a classical HET discharge, play an essential role in determining the spectral characteristics and phase-frequency relationships of relevant plasma magnitudes. Accurate numerical simulations are of particular importance in this context, since experimental measurements of plasma profiles inside thruster chamber are scarce and challenging. Reduced 1D axial simulation models used in previous studies [33,37] reproduced the main features of the experimental results, although they present several important limitations. First, they have not been able to reproduce accurately the natural BM of the thruster; second, they miss important 2D effects and are not directly applicable to complex MS topologies of new low erosion MS-HET designs; third, they rely on approximate models for plasmawall interaction, which determine plasma losses to the walls and thus the discharge performance; and fourth, they do not include double ionization collisions, whose contribution is generally non-negligible.
This work presents 2D numerical simulations of the dynamic response of a MS-HET operating with an external sinusoidal modulation of V d . The simulations are carried out with the code HYPHEN [38,39], a 2D axisymmetric multi-thruster simulator with a hybrid formulation featuring a particle-in-cell model for heavy species (including neutrals, singly and doubly charged ions), and a magnetized driftdiffusive fluid model for electrons. In a recent work [40], HYPHEN capabilities were extended to the simulation of MS-HET discharges, and partial code validation was made against the HT5k prototype [41,42], in terms of discharge current and thrust for five different operation points, without any external modulation and using xenon as propellant. Taking as reference the simulation results for one of these unmodulated operation points, a parametric investigation of the modulated response of the thruster for different modulation amplitudes and frequencies is presented. The central goal is to assess the effect of modulation in this new, high-efficiency, 5 kW-class MS-HET prototype, which greatly differs from HET designs analyzed in previous modulation studies, and to advance in the understanding of the externally driven BM operation, which could be advantageous for high-power applications. The simulation results presented here are shown to be in good qualitative agreement with those reported in previous numerical and experimental studies, thus confirming that modulated MS-HET discharges also follow main trends found for conventional HETs. The 2D solution for the electron temperature, plasma density, neutral density and ionization production term inside the thruster chamber is analyzed to explain the impact of modulation in the discharge performance, and the transition between the externally driven and the natural BM is investigated. To the authors' knowledge, there are no previous full 2D numerical studies on a modulated MS-HET discharge.
The document is organized as follows. Section 2 summarizes the simulation model and settings and presents the natural BM of the thruster characterizing the unmodulated discharge. Section 3 discusses the effect of voltage modulation on main performance figures when compared to the unmodulated case. Parametric studies in terms of modulation amplitude and frequency are presented in sections 3.1 and 3.2, respectively. Section 4 analyzes the 2D plasma solution to get a deeper insight on the dynamics of the modulated response. Section 5 studies the transition between the natural and the externally driven BM. Finally, conclusions are drawn in section 6.

The natural BM response to constant discharge voltage
The present study uses the same simulation model and main settings detailed in [40], which are briefly outlined here. Figure 1 presents a sketch of the simulation domain, which corresponds to the cylindrical axisymmetric half meridian plane of the thruster, including the thruster chamber, with length L c and width H c , and the near plume region, with axial and radial extensions equal to 6L c and 6H c , respectively. The electric circuit connecting the anode and the centrallymounted cathode is also depicted in figure 1. Both the anode wall and cathode exit planes correspond to boundaries of the simulation domain. A discharge voltage V d is imposed between the anode and the cathode. The reference ϕ = 0 for the electric potential is set at the cathode boundary, so the anode potential is ϕ = V d . The discharge current I d flowing from cathode to anode through the external circuit is indicated in figure 1 according to the flow direction of positive charges.
A prescribed xenon neutral mass flowṁ A is injected through the anode wall. The xenon neutral mass flow injected through the central cathode isṁ C = 0.075ṁ A . Simulations include singly-charged and doubly-charged ions generated through ionization collisions in the bulk plasma. Plasma recombination at the thruster walls contributes to neutral density. Charge-exchange collisions generating radiallyexpanding, slow ions and fast neutrals in the near plume are not considered in this study, since their effect on the discharge performance was found negligible for the present case.
The code version for HET simulations consists of three main modules coupled within a time-marching sequential loop: the Ion module (I-module), which follows a Lagrangian approach for simulating the neutral and ion species; the Electron module (E-module), which applies quasineutrality and solves a fluid model for the magnetized electron population; and the Sheath module (S-module), which provides the proper coupling between the quasineutral plasma bulk, and the thruster walls. While the I-module operates on a structured mesh of the simulation domain, the E-module uses an unstructured magnetic field aligned mesh [43], defined by the externally applied magnetic field B to limit the numerical diffusion arising from the strong anisotropic transport of magnetized electrons. The details of both meshes are provided in [40] and omitted here.
Every simulation timestep, the I-module takes as inputs the externally applied magnetic field B, the electric field E = −∇ϕ, and the electron temperature T e and performs the following tasks: (i) the propagation of macroparticles one timestep ∆t forward, according to the electromagnetic fields acting on them; (ii) the injection of new macroparticles into the domain and the removal of existing ones; (iii) the interaction of macroparticles with the thruster walls, such as neutral reflection and ion recombination; (iv) the generation of new ion macroparticles due to the ionization of neutrals; and (v) the computation, through a particle-to-mesh weighting process, of the macroscopic properties characterizing each heavy species, including the neutral density n n , the plasma density n e = s̸ =e Z s n s and the net ion current density j i = s̸ =e eZ s n s u s , with Z s , n s and u s the charge number, particle density and macroscopic velocity of the sth ion population, including singly and doubly charged ions. Further details can be found in [38,44,45]. The E-module, taking these heavyspecies magnitudes as inputs, solves a quasineutral, driftdiffusion fluid model for the magnetized electron population, including the plasma current conservation, a generalized Ohm's law solving for ϕ and the electron current density vector j e = −en e u e , and the electron energy equation solving for T e with a Fourier's law for the electron heat flux vector q e . The electron model is further detailed in [38,46,47], and it has been recently upgraded for the simulation of MS-HET discharges in [40].
HYPHEN E-module models the azimuth-averaged, wavebased electron anomalous transport through a turbulent electron collisionality ν t = α t ω ce , where α t (z, r) is a phenomenological function representing the electron local turbulence level and ω ce is the electron cyclotron frequency [6,[48][49][50]. In [40], five different operation points were considered, defined by a pair (V d ,ṁ A ). The former ranged from 300 to 400 V, and the latter from 10 to 14 mg s −1 . For each operation point, only two experimental measurements were available: the time-averaged values of the thrust and the discharge currentF andĪ d , respectively. Therefore, to complete the electron model, an axial stepout type function for α t (z, r), which has provided good fittings in previous studies [13,[51][52][53], was chosen, featuring two fitting parameters only: α t1 and α t2 (> α t1 ) acting approximately inside and outside the thruster chamber, respectively. For each operation point, the pair (α t1 , α t2 ) was tuned to reproduce the experimental values (Ī d ,F) with relative errors below 5%. The obtained turbulent fittings yielded numerical results in good agreement with experimental measurements of plasma magnitudes for a similar prototype [54,55], and turned out to be also good reproducing the main oscillations of I d , which correspond to the natural BM of the thruster, as described later. Previous studies have successfully developed more advanced empirically-calibrated models for α t (z, r), informed by local plasma properties [56][57][58]. Moreover, noninvasive time-resolved measurements have recently shown a direct correlation between the anomalous collision frequency and breathing cycle [59]. Despite current remarkable efforts, the precise characterization and modeling of the anomalous electron transport still remains an open problem in HET research, subject to several challenges that block the development of predictive simulation tools, including the nonuniqueness of empirically-calibrated more elaborated models [56,60]. Here, the experimental data availability for the MS-HET prototype limits the tuning of more complex turbulent models.
The present work takes as reference the simulation case 1 of [40]. For simplicity, the RLC filter unit connecting anode and cathode in [40] is not included in this work (since its effect on I d was proven negligible). Therefore, this unmodulated simulation case features V d = 300 V (constant with time) andṁ A = 14 mg s −1 . For this case, the fitting process yields (α t1 , α t2 ) = (0.8, 8)%. The oscillations of I d are shown in figure 2(a). The time-averaged value of the discharge current isĪ d = 15 A, and it features a half-peak amplitude ∆I d of about 6.7% ofĪ d , consistent with experiments. Figure 2(b) shows the normalized frequency spectrum of I d .
The dominant frequency f BM = 20 kHz corresponds to the natural BM, which is close to the 22 kHz reported in experiments. The frequency spectrum in figure 2(b) reveals the spectral complexity of I d , with several peaks at frequencies close to f BM . This result suggests that the natural discharge oscillations can be represented by a non-harmonic cluster of several modes, similar to those identified in a recent study through high order dynamic mode decomposition (HODMD) for a conventional SPT-100-like HET discharge [61]. The application of HODMD to this unmodulated reference case is left for section 4. In the rest of the paper, simulation cases including V d modulation are compared with this unmodulated reference case.

Effect of voltage modulation on main performance figures
The sinusoidal modulation applied on the anode-to-cathode discharge voltage (refer to figure 1) is defined as whereV d is the constant time-averaged discharge voltage value, ε is the relative half-peak modulation amplitude, and f s is the modulation frequency.
In order to compare with the unmodulated reference case presented in section 2, all modulated cases analyzed in this work consider the same operation point, characterized by the pair (V d ,ṁ A ) = (300 V, 14 mg s −1 ). Each modulation case is defined by the pair (f s , ε); the unmodulated reference case will be referred to as f s = 0 or ε = 0. Sections 3.1 and 3.2 present parametric studies on ε and f s , respectively. Time-averaged values of main performance figures, including a valuable set of partial efficiencies (defined below), are reported for all simulated cases.
Computing efficiencies in a highly oscillating HET discharge requires special care. Here, meaningful estimations of the thrust efficiency and a set of partial efficiencies will be obtained considering time-averaged quantities over a sufficiently large number of complete discharge oscillation cycles. The thrust efficiency is defined and factorized, in terms of the time-averaged variables, as beingP d = I d V d the time-averaged discharge power deposited into the plasma (the net power delivered through cathode electron emission is around 1%-2% ofP d and has been neglected here). The propellant utilization and the current efficiency are defined as respectively, whereṁ =ṁ A +ṁ C ,m i ∞ andĪ i∞ are the timeaveraged ion mass flow and ion beam current across the downstream plume boundary, respectively, and m i and e are the xenon atom mass and the electron charge, respectively. The voltage, divergence, and dispersion efficiencies are defined, respectively, as withP ∞ andP z∞ the time-averaged flows of total and axial plasma energy across the downstream plume boundaries, respectively, and eP ∞ /Ī i∞ the time-averaged effective energy of ions downstream. The values ofF,P ∞ andP z∞ are computed from the corresponding surface integrals at the downstream plume boundary considering all plasma species (the contributions of neutrals and electrons are found negligible in all cases). The product η cur η vol quantifies the fraction ofP d transferred to the downstream plume; η div assesses the plume divergence based on axial energy and total energy flows; and η disp quantifies the level of velocity dispersion of all plasma species (which would be one if only axial mono-energetic ions contribute to thrust and power downstream). Additionally, the charge efficiency, is equal to one if no doubly-charged ions are generated, and decreases with the production of doubly-charged ions.
To show the effect of discharge oscillations on the performance figures estimations, and to establish a fairer comparison with the results reported in [33], RMS values of I d and F, and estimations of the efficiencies above using RMS values of the involved quantities are also reported in this work. For any time series y 1 , . . . , y K ∈ R, its RMS value is defined as

Modulation amplitude parametric study
A parametric investigation on the effects of the modulation amplitude is presented in this section. For all modulated simulation cases analyzed here, the modulation frequency is set equal to that of the natural BM characterizing the unmodulated discharge presented in section 2, so that f s = f BM = 20 kHz. Results for four modulated discharges with relative half-peak amplitude ε =2.5%, 6.25%, 12.5% and 25% are presented and   Table 1 contains these data only for the unmodulated case. A monotonic growth of I d and F RMS values and their oscillation amplitudes with the modulation amplitude is found, in line with the results reported in [33]. As expected, RMS values are noticeably higher than their corresponding (arithmetic) timeaveraged values when signal oscillations are large. In [27,33], the linear and non-linear responses of a 200 W CHT prototype against the amplitude of the anode voltage modulation were identified. In the linear regime, at low modulation amplitude (up to ε ∼ 5%), the amplitude of the I d oscillations increases linearly with the modulation amplitude while the RMS values remain nearly constant. In the non-linear regime, at higher modulation amplitudes, I d oscillation amplitude and RMS values exhibit a superlinear growth with the modulation amplitude. Interestingly, these two regimes are also observed in figure 3(a).
The time-averaged and RMS values for the enhanced F shown in figure 3(b) in the non-linear regime are consistent with the results presented in [35][36][37], where it is shown that, for modulated operation, an increase in F can be obtained when oscillations of the ion beam current and energy are high and their phase shift is low. These aspects are analyzed in [36] and [37] considering modulated discharges with different modulation frequencies (close to the natural BM) at a constant modulation amplitude, an analysis which is undertaken in the next section. Figure 4 compares the total and partial efficiencies of modulated cases with the unmodulated case of table 1, identified there by an asterisk. The trends of η u and η cur with the modulation amplitude are in good qualitative agreement with the experimental measurements and numerical results reported in [33] (propellant utilization values reported there correspond to our ratio η u /η ch ).
Section 4 will show that ionization is enhanced by external modulation, thus corresponding to a larger η u . A similar behavior, but with slighter changes, is found for η cur . The increase in η cur (especially in the non-linear regime) is due to an enhancement of the ion beam current I i∞ with the voltage modulation. 1D numerical simulations and experiments reported in [33,37] indicate that the enhanced I i∞ is related to the increase in ε and the phase alignment of the ion density and energy at the near plume (i.e. end of the acceleration region). The same behavior (not shown) is found in our simulations. Finally, since η/η * < (η u η cur )/(η * u η * cur ), the product η vol η div η disp decreases with ε increasing, thus concluding that voltage modulation has a small effect on thrust efficiency due to the different trends of the partial efficiencies.
The comparison of solid lines (time-averaged values) and dashed lines (RMS values) in figure 4 indicates that using RMS values to estimate efficiencies does not seem appropriate in scenarios with high discharge oscillations, and it may lead to significant misestimation of the performance figures. This is particularly the case for η u defined in equation (3), which is reported to increase up to a 6% if estimated with RMS values, and only up to a 2% if obtained with time-averaged quantities. Discrepancy between these particular values is induced by oscillations inṁ i∞ exclusively.

Frequency parametric study
Setting the amplitude modulation to ε = 25%, this section presents simulation results for f s ranging from half to double the dominant frequency of the natural BM, f BM = 20 kHz.  I d frequency spectrum than that of the unmodulated one in figure 2: a single prominent peak at f s is found, corresponding to highly coherent I d oscillations. This result is in line with experimental studies reporting that sinusoidal voltage modulation close to the natural BM yields quasi-periodic discharge oscillations [23,27,28].
In figure 5(b) the modulated discharge exhibits a clear capacitive/inductive character for f s = 10/40 kHz, with I d ahead/behind V d , respectively. A similar behavior has been reported in experimental studies for a 1 kW-class HET operated with a chopper PS [29], and in experimental and numerical studies for the CHTpm2 prototype with sinusoidal voltage modulation [37].
As already seen in section 3.1, modulated cases feature higher I d oscillation amplitude than the unmodulated one. To further investigate this aspect, figure 6 shows the I d and F response for several frequencies between 10 and 40 kHz. Interestingly, the oscillations of I d and F exhibit a nonmonotonic behavior with f s . In fact, the modulated response exhibits a resonant-kind behavior similar to that reported in experiments and 1D numerical simulations for the 200 W electromagnetic CHT [33] and its permanent magnet version CHTpm2 [36,37]. This resonant response of the thruster is characterized by enhanced I d and F oscillations for an interval of f s close to f BM .
In [33], the resonant modulation frequency f r , is defined as the modulation frequency for which the I d oscillation amplitude and its RMS value are maximum. The value of f r is shown to depend on the modulation amplitude, being slightly lower than f BM for high modulation amplitude (ε = 9%), and approximately equal to f BM for low modulation amplitude (ε < 5%). At resonance, for both the electromagnetic and the permanent magnet versions of the 200 W CHT, it is found ∆I d /I d ∼ ±100% for ε = 9% [33] and ε = 18.2% [37], respectively. For the latter,Ī d andF are found to increase a 6% and 4% with respect to the unmodulated case, respectively [37].
Here, for ε = 25%, the maximum oscillation amplitude of I d is ∆I d /I d ∼ ±50% at f s = 18 kHz, while RMS values are maximum at f s = 17 kHz. Therefore, there exists a resonance region centered at a resonant frequency f r ∼17-18 kHz, slightly lower than f BM = 20 kHz, with a width of 2-3 kHz, shaded in blue in figure 6. Within this resonance regionĪ d andF are shown to increase a 3% and 5% with respect to the unmodulated reference case, respectively. Interestingly, doubling and halving f s with respect to f BM leads to roughly the same decrease in I d oscillation amplitude and RMS value. At lower modulation amplitudes, f r is found to approach f BM . In particular, for ε = 6.25%, we find f r = f BM . These results reveal the prominence of the natural BM of the thruster over the external modulation on the low-frequency discharge oscillations for low modulation amplitude. The transition between the natural and externally driven BM at lower modulation amplitude is detailed in section 5. Figure 7 shows the changes in efficiencies for modulated cases with respect to f s = 0 versus the modulation frequency. Slight changes (of about 3%-4% maximum) are found in all efficiency figures, in line with the results in [37]. The nonmonotonic behavior of the thrust efficiency η close to resonance is also in good agreement with numerical and experimental results reported there. The maximum η gains with respect to the unmodulated reference case occur within the resonance region and are limited to ∼2% only, as a consequence of counteracting contributions of all the partial efficiencies in equation (4), as described next.
The maximum of η u is close to the resonance frequency, as in [33,37]. The current and voltage efficiencies, η cur and η vol , do not exhibit here clear trends with f s . Experimental studies in [33] do not report a clear trend of η cur either, while η cur is found to remain nearly constant in [37]. All modulated cases feature a smaller η disp than the unmodulated case, with a minimum in the resonance region, thus indicating a larger velocity dispersion in the ion velocity distribution function (VDF) for those cases.
For a given species, the time-averaged axial flux-VDF of particles leaving the domain through the downstream boundary is defined as with F(v r , v θ , v z ) the species VDF. Figure 8(a) shows F z for singly-charged ions at the downstream boundary at the channel midline for f s = 0, 10, 18 and 40 kHz. For the unmodulated case, F z presents a single dominant peak at 290 eV and a much lower ion velocity dispersion (consistent with a higher η disp ) than for the modulated cases. For all modulated cases, the ion-flux VDF exhibits two main peaks near 225 eV and 375 eV, the minimum and maximum values of the modulated V d (t). Interestingly, the dominant peak is the one with the highest energy. Compared to cases at f s = 10 and 40 kHz, the results for f s = 18 kHz show a lower velocity dispersion of the singly-charged ion population, in line with [35]. A similar double peak behavior for modulated cases is found for doublycharged ions.
The behavior of η disp is also related to that of η ch (figures 7(e) and (g)), since the presence of doubly-charged ions contributes to the velocity dispersion of the ion population. The lower η ch values near resonance indicate that a larger fraction of doubly-charged ions is produced at these conditions. This is consistent with the higher T e around f r reported by experiments [37]; this aspect is further analyzed in section 4. The decrease in η div within the resonance region indicates higher plume divergence, in line with results in [37], and, along with η disp decrement, partially compensates the gains in η u η cur .
The increase in the product η cur η vol η div ≈P zi∞ /P d is also limited, whereP zi∞ corresponds to the dominant ion contribution toP z∞ (refer to equation (4)). This aspect is further analyzed as follows. For sinusoidally oscillating quantities,P d andP zi∞ can be expressed as where E zi∞ = eP zi∞ /I i∞ is the average axial energy per ion crossing the downstream plume boundary; and φ d , φ ∞ ∈ [−180, 180] deg are phase shifts between I d and V d , and between I i∞ and E zi∞ , respectively, referred to as discharge phase and ion phase. A negative phase shift means that the first signal of the pair (the currents, in this case) is ahead of the other signal. Assuming that time-averaged values remain nearly constant with f s , and that φ d , φ ∞ ∈ (−90, 90) deg (i.e. cos φ d , cos φ ∞ > 0), bothP d andP zi∞ increase for higher oscillation amplitude and phase alignment of the involved signals.   Figures 9(b) and (c) present the time response of I i∞ and E zi∞ , respectively, while the one corresponding to I d is in figure 5(b). Two main facts are revealed by these results. First, the behavior of I i∞ is found similar to that of I d : it exhibits a non-monotonic behavior in terms of oscillation amplitude and time-averaged values, with both quantities peaking near resonance; and second, the absolute value of both φ d and φ ∞ are shown to decrease near resonance, indicating a phase alignment of signals I d and V d , and I i∞ and E zi∞ , respectively. Indeed, both curves φ d , φ ∞ (f s ) are nearly coincident, and full phasing (i.e. φ = 0) is found for f s ∼ 25 kHz, thus above the resonance region (and f BM ). The results for φ d (f s ) confirm the capacitive/inductive character of the discharge commented above for f s lower/higher than 25 kHz. A similar behavior is reported in [37]. These findings indicate that near resonance bothP zi∞ andP d increase, thus limiting the product η cur η vol η div ≈P zi∞ /P d and, consequently, the gains in η induced by modulation.
Previous studies for a voltage modulated thruster [35][36][37] have shown that high oscillation amplitude and low phase shift between ion current and energy is optimal for thrust production. Here, we find maximum F near resonance (refer to figure 6(b)), where I i∞ oscillations are maximum and φ ∞ is low. While average values of E zi∞ (about 251-253 eV) are found to barely change with f s (maximum change is about 0.8% with respect to the unmodulated case) its oscillation amplitude is shown to decrease at resonance with respect to its level at low f s (see black and red curves in figure 9(c) corresponding to f s = 10 and 18 kHz, respectively). This fact, along with the higher plume divergence discussed above, limits thrust production, in line with [37]. Furthermore, the combination of the sub-optimal F with the aforementioned enhanced P d yields limited η gains for modulated resonant discharges.
Finally, figure 8(b) shows F z for singly-charged ions at the downstream boundary at the thruster channel midline computed at the time instants with maximum and minimum E zi∞ for a given discharge oscillation cycle. Results are presented for the unmodulated reference case described in section 2 (blue curves) and for the modulated case with f s = 18 kHz (red curves). The axial flux-VDFs present two peaks at energies consistent with the time response of E zi∞ in figure 9. Interestingly, while both the high and low-energy peaks in the unmodulated case are similarly populated, the high-energy peak clearly dominates in the modulated case. This result indicates that voltage modulation near resonance enhances and optimizes the acceleration of the ion population, thus favoring F production, as commented above. These aspects are further analyzed next.

Further insights on the 2D plasma response
In this section we analyze the time-dependent response of the plasma variables in order to discuss the predator-prey character of the driven modes. First we will consider a 0D analysis with global variables, then we will look at 1D and 2D behaviors.

Global response
In order to define globally-averaged values of the plasma variables we consider the domain in figure 1 bounded by the red line, which covers the portion of the simulation domain where most ionization takes place, extending axially up to z/L c =1.4, thus including the thruster chamber and the position of T e and E z peaks downstream [40]. Spatially-averaged values within this region of any magnitude ζ are denoted asζ.
The ion continuity equation states where is the net ionization production term, including single and double ionization between electrons and heavy species (ions and neutrals) with, for each ionization process κ (i.e. Xe → Xe + , Xe → Xe 2+ and Xe + → Xe 2+ ): n κ the heavy species involved, ∆Z κ the heavy species charge number jump in the process, and R κ (T e ) the corresponding ionization rate, which is a non-linear function of the electron temperature. Ion current oscillations are therefore governed by oscillations of S e , which in turn depend on the dynamics of T e , n n and n e . Current continuity in the quasineutral plasma relates the produced ion and electron current densities, The surface integral of the electric current (i.e. j = j i + j e ) collected at the anode wall provides I d , while the surface integral of the ion current collected at the downstream plume boundary provides I i∞ .
On the other hand, T e is determined from the inertialess electron energy equation for an isotropic pressure tensor, (13) where P ′ ′ e = −5T e j e /(2e) + q e is the electron energy flux gathering the enthalpy and heat fluxes; P ′ ′ ′ elec = j e · E is the work of the electric field over fluid electrons per unit volume and time, corresponding to the main energy source term for electrons; and P ′ ′ ′ inel accounts for the power losses per unit volume from inelastic collisions (e.g. excitation and ionization). Figure 10 shows the phase shifts and the time response of several plasma variables for the same simulations than figures 5 and 9. Interestingly, the results for the modulated discharge found here for the MS-HET are in good qualitative agreement with those reported in [37] for the CHTpm2. This fact suggests that the same fundamental phenomena drive the modulated response in both prototypes, and indicates that anode voltage modulation seems a robust technique to control the frequency and the phasing between current and voltage across different HET designs.
For modulated cases, a quasi-periodic time response with dominant frequency equal to f s is found for all magnitudes in figure 10, with significantly larger oscillation amplitudes than for the unmodulated case, in line with the behavior of I d . The ionization of neutrals is triggered by a combination of the rise ofT e and the level of neutral replenishment, yielding an increase in plasma density until the neutral population is depleted andT e decreases. The ionization then drops until neutrals are replenished andT e rises again, leading to another cycle. For all cases,ñ n is ahead ofñ e , the phase lag between signals ranging from 100 • to 150 • . This behavior is typical of a predator-prey type ionization instability [5,6], and suggests that the main discharge oscillations in modulated cases correspond to an externally driven BM (this aspect is further assessed later). For low and high f s , far from the resonant region,ñ n -V d phase shift exhibits an asymptotic behavior, typical of resonant oscillators. A similar trend is found forñ e -V d andS e -V d phase shifts.
For f s = 18 kHz,S e exhibits much larger oscillations yielding higher neutral depletion and plasma density peaks. The enhanced averaged plasma production (it increases up to a 6% from the unmodulated case) is in line with the increase in η u with respect to the unmodulated discharge observed in figure 7(a). Near resonance, theS e -V d phasing ( figure 10) has a central impact on the discharge performance, as described next. First, it is responsible for both the discharge and ion phasing (i.e. the decrease of |φ d | and |φ ∞ |) discussed in section 3.2, since (1) S e is the source of both ion and electron currents contributing to I d and I i∞ (refer to equations (10) and (12)), and (2) V d ultimately determines E zi∞ in figure 9(c). The similar phase-frequency characteristic found for φ d and φ ∞ in figure 9(a) is a consequence of the coupling between the plasma generation and acceleration in the HETs, which limits the performance improvement in the modulated discharge (especially in terms of η). On the one hand, both F and P zi∞ are maximized ifS e and V d are maximum at the same time every discharge cycle. On the other hand, the higher F (or P zi∞ ) is accompanied by higher P d .
Considering only modulated cases, theS e -V d phasing near resonance also contributes to the lower ion velocity dispersion observed in figure 8(a) for f s = 18 kHz (red line): if the peaks of bothS e and V d occur closer within the cycle, then the majority of the ions will go across a similar potential fall when accelerating downstream (∂V d /∂t ≈ 0 near its maximum). Furthermore, this phasing is behind the fact that the high energy peak of F z in figure 8(b) is much more populated than the lower one.
The electron temperature response is determined by the energy balance for electrons in equation (13), including P ′ ′ ′ elec as the main energy source, and is found to oscillate nearly in-phase with V d for all f s values. This is consistent with the fact that (1) MS of the channel walls limits power losses to the walls and avoids T e saturation [62]; and (2) the period of the modulation (∼10 −4 s) is much longer than the electron residence time (∼10 −6 -10 −7 s), so that j e quickly adapts to V d changes (with f s ). The largerT e peak near resonance (it increases from about 20 to 24 eV from the unmodulated case to the case with f s = 18 kHz) favors double ionization collisions, and is thus in line with the slight decrease in η ch and η disp , as indicated in section 3.2.
For f s near f BM , T e is found to govern S e (and n n ) dynamics and, therefore, is the main factor responsible for the control of the dominant frequency of the modulated discharge oscillations. As commented above,T e oscillates nearly in-phase with V d for all f s values ( figure 10). On the other hand,ñ n oscillates nearly out-of-phase with V d for f s = 10 kHz, and their phase shift decreases (in absolute value) for higher f s . The rise inS e depends on that ofT e and the level of neutral replenishment in the chamber. For f s = 10 kHz, neutrals have longer time available to replenish the chamber beforeT e begins to rise, and they trigger the ionization before the rise ofT e . For f s = 18 and 40 kHz, neutrals have shorter time to replenish the chamber, and a lower value ofñ n is found at the initial rise ofS e , which is now triggered mainly by the rise ofT e . Near resonance, at f s = 18 kHz, the ionization is enhanced by an optimal combination ofñ n andT e , yielding the maximum peak ofS e andñ e , and the highest level of neutral depletion. Figure 11 shows the axial-temporal contour maps over several modulation cycles of radially-averaged plasma magnitudes (denoted with circumflex accent). For f s = 10 kHz, the ionization front, identified with the peak ofŜ e , is formed on every cycle at z/L c ≈ 0.8, and moves upstream towards the anode along the cycle. The larger neutral replenishment in the chamber in this case triggers the ionization close to the chamber exit, whereT e is higher, and the ionization front is formed well beforeT e reaches its maximum value. Then, the ionization front travels upstream, as neutrals are consumed (and the electron temperature increases), and stays near the anode a significant part of the cycle, untilT e , which oscillates nearly in-phase with V d , decreases with time and the ionization front vanishes. For the case with f s = 18 kHz the ionization front appears at a distance from the anode similar to the case with f s = 10 kHz, while for f s = 40 kHz the front forms closer to the anode (at x/L c ≈ 0.7) due to the lower level of downstream neutral replenishment. For both f s = 18 and 40 kHz, the ionization front is triggered by the rise in the electron temperature. Near resonance (e.g. at f s = 18 kHz), an optimal phasing between neutral replenishment and electron temperature is found every cycle, yielding the highest ionization production, as commented above. The upstream motion of the ionization front towards the anode is much more evident for f s = 18 than for f s = 40 kHz, and in both cases the ionization extinguishes when the electron temperature drops within the cycle. These results confirm that, for modulated cases with f s close to f BM , the dynamics of T e are the main factor responsible for slowing down or speeding up the inherent ionization cycle of the unmodulated discharge, yielding an externally driven BM at f s .

Local response
To have a full 2D picture of the dominant spatio-temporal modes we apply the HODMD [63,64] data-driven technique, which expands a spatio-temporal dataset into dynamicsrelevant modes identifying involved frequencies, amplitudes and growth rates, and that allows to separate transient behaviors from asymptotic oscillations. HODMD extends standard DMD [65] capabilities for the application to non-linear system with high spectral complexity, and overcomes the limitations of standard DMD when dealing with noisy data. A recent study [61] has already demonstrated the capabilities of HODMD to identify and characterize the natural BM and the ion transit time [7,13,66] from HYPHEN simulations of a SPT-100-like HET.
Given a set of spatio-temporal data composed of K snapshots of dimension N, i.e. x 1 , . . . , x K ∈ R N , each of them representing the values of a physical variable of interest x at the N spatial mesh points involved in the numerical simulations for a given time instant k, i.e. x k = x(t k ), HODMD aims at decomposing each snapshot x k as where ψ m are the complex spatial modes (normalized with RMS-norm), Ω m their corresponding complex frequencies and a m > 0 are their real amplitudes. Here, only purely oscillating modes are considered for the analysis and a m assesses their dynamical relevance. The accuracy of the HODMD reconstruction of the original data, x k , depends on the number M of modes considered and a tunable parameter d that defines the number of 'delayed' snapshots to be used in the analysis [63]. In order to identify an adequate value for M and d, a similar procedure to the one in [61] is followed: M and d are varied until the relative RMS error in the reconstruction is found small enough (between 5% and 15%); once this is achieved, it is checked that small variations of the resulting d parameter do not affect the output HODMD modes. Results shown next are for M = 40 (including complex conjugates) and d =600. Figure 12 shows the normalized amplitude spectrum resulting from the application of HODMD to the spatio-temporal data of S e (z, r, t) within the red-bounded region of figure 1 for unmodulated and modulated cases. For the unmodulated case, the largest-amplitude mode corresponds to the natural BM frequency f BM = 20 kHz, but there are non-negligible contributions of several modes with frequencies close to f BM , and which are not harmonics of f BM . In particular, the amplitudes of the second and the fifth-largest modes, at 29.5 and 70.2 kHz, respectively, are just 18% and 42% lower than that of the largest one. Therefore, the natural BM cannot be reconstructed within a selected accuracy by a simple spatio-temporal sinusoidal term (nor the superposition of several sinusoidal harmonics). Instead, a wide non-harmonic BM cluster is found, which confirms the rich spectral complexity of the natural BM for this thruster, in line with the I d time response and frequency spectrum shown in figure 2.
Modulated cases show a simpler spectrum, with a more prominent role of the mode at f s . This result is consistent with the quasi-periodic I d response shown in figure 5 for these cases. Yet, cases with f s = 10 and 18 kHz show non-negligible contributions of two to three harmonic modes, and several non-harmonic secondary modes with frequencies close to f s .
Interestingly, for f s = 40 kHz (> f BM ), the mode-amplitude hierarchy contains modes with frequencies close to f BM . Figure 13 shows the spatial structure of the dominant HODMD mode of S e , T e , n n , and n e , in terms of its magnitude and phase angle within the region of analysis depicted in figure 1. The first three rows of figure 13 show results for the unmodulated reference case. 2D contour maps of mode magnitude and phase angle are depicted in the first and second row. The third row shows 1D axial profiles of magnitude (black lines and left y-axis) and phase angle (blue lines and right y-axis). While dashed lines correspond to values along the thruster channel midline, solid lines, corresponding to radially averaged values, are included to emphasize the spatial non-uniform 2D character of the modes. The reference for the phase angle (i.e. 0 • ) is set to the phase of the S e mode at the midpoint of the anode wall. Results for the modulated discharge with f s = 18 kHz are found qualitatively similar to those of the unmodulated discharge. Therefore, 2D contour maps are omitted, and only 1D axial profiles of magnitude and phase angle are reported in the fourth row of figure 13, in the same fashion as for the unmodulated case. The information in figure 13 must be interpreted as follows: the magnitude plots indicate where the oscillation takes place; the phase plots describe the standing or progressive structure of the wave and allow identifying the phase relation between the different plasma variables.
Magnitude 2D maps and 1D axial profiles in the first column of figure 13 show that, for both the unmodulated reference case and the modulated case with f s = 18 kHz, the region where S e oscillates extends from the anode to z/L ≈ 0.8. This result is consistent with the spatio-temporal contour maps ofŜ e depicted in figure 11(b) for f s = 18 kHz, showing that the ionization front, identified with the peak ofŜ e , forms at z/L ≈ 0.8 and travels upstream towards the anode during the modulation cycle. In fact, in this region S e exhibits an upwards progressive-wave structure (i.e. traveling to the left in the diagrams), indicated by the positive slope of the 1D axial phase profiles (blue lines in the third and fourth rows of the first column of figure 13). Its characteristic axial phase velocity is 2.68 km s −1 for the unmodulated reference case and 2.63 km s −1 for the modulated case with f s = 18 kHz. The main oscillations of n n and n e (third and fourth columns of figure 13) are also found between the anode and z/L ≈ 0.8, indicating that the ionization region is contained there. In contrast, T e is found to oscillate downstream the ionization region, mainly outside the thruster chamber, where most of the ion acceleration takes place. A similar behavior is found for E z (not shown). T e has essentially a constant phase through the analyzed region, indicating a dominant standing-wave character. The 1D axial profiles along the channel midline for the phase of n n and n e indicate that the former leads by 75 • -135 • within the ionization region, and that both present a dominant standing-wave character.
The spatial non-uniformity of the modes in the thruster channel is especially evident for the case of n n . The radial nonuniformity of S e is mainly due to the ion recombination at the thruster walls, which significantly contributes to n n , and it is partially responsible for that of n e . For the n n mode, radially averaged 1D profiles of magnitude and phase along the channel (solid lines) greatly differ from those along the thruster channel midline (dashed lines). Interestingly, the radially averaged n n mode corresponds to an apparent 1D axial wave travelling downstream along the thruster channel from the anode wall (i.e. solid blue lines with negative slopes). This behavior is consistent withn n maps in figure 11, in which the downstream travelling of neutrals is more evident for cases with f s = 10 and 18 kHz. The phase characteristics of this axial wave allows to identify two different regions. On the one hand, near the anode, the characteristic axial phase velocity is mainly a result of the competition between neutral injection and ionization: for the unmodulated case, which features lower ionization near the anode, the HODMD mode for n n recovers the expected convective axial wave dominated by neutral injection from the anode, and the phase velocity is 0.31 km s −1 , close to macroscopic axial injection velocity for neutrals, set to 0.3 km s −1 ; for the modulated case with f s = 18 kHz, the phase velocity is of the same order of magnitude but smaller, 0.15 km s −1 , due to the enhanced ionization and the phasing between S e and n n , which decrease the effective neutral replenishment rate. On the other hand, further downstream the anode wall, for z/L c > 0.2-0.3, neutral fluxes from ion recombination along the lateral walls yield an apparent axial phase velocity of 1.76 km s −1 (unmodulated case) and 1.51 km s −1 (modulated case), about five to six times larger than the macroscopic axial velocity of neutrals.

Transition from driven to natural BM
Voltage modulation at relatively high modulation amplitudes (e.g. ε = 25%) was found to yield limited gains in thrust efficiency, but, it has been proven to be an effective strategy to control the frequency of the discharge oscillations for f s ≈ f BM . While this feature can help reduce EMI on modern spacecrafts, the enhanced I d oscillations obtained in the non-linear regime ( figure 3(a)) may pose issues to other subsystems. Therefore, anode modulation at relatively low amplitudes within the linear regime (e.g. ε = 5% − 10%) seems an interesting compromise allowing to control the frequency of the oscillations while keeping low current oscillation amplitudes. However, operating at low ε may result in a loss of control over the discharge oscillations and the reappearance of a dominant natural BM if f s is sufficiently far from f BM [30].
For a given modulation amplitude, the natural BM is expected to dominate over driven oscillations if f s far enough from f BM .
The transition between the driven and the natural BM of the thruster is here investigated for a relatively low modulation amplitude within the linear regime, set to ε = 6.25%. Two limit simulation cases are run with f s = 3 kHz and 110 kHz, one order of magnitude lower and larger than f BM , respectively. Figure 14 shows the normalized frequency spectrum of I d and main plasma magnitudes for f s = 3 kHz and 110 kHz. The results for I d reveal that the natural BM of the thruster at f BM = 20 kHz is present in the response as the dominant mode (for f s = 3 kHz) or a co-dominant one (for f s = 110 kHz). Therefore, voltage modulation fails to control the main oscillations of the plasma discharge for both f s = 3 kHz and 110 kHz. Moreover, I d response deviates from a quasi-periodic one at f s , characterizing modulated discharges with f s ≈ f BM . Instead, a higher spectral complexity is found, similar to that of the unmodulated case.
As discussed in section 4, when f s is close to f BM , T e is found to govern the dynamics of S e , thus enabling the control of the discharge oscillation frequency. Here, for f s far from f BM , a quasi-periodic time response at f s is still found forT e . However, the main oscillations ofS e occur now at f BM , thus uncoupled from T e dynamics, and its response presents a higher spectral complexity, similar to that of I d , with a secondary mode at f s . Similar results are found forñ n andñ e . Therefore, the natural ionization cycle of the discharge at f BM , governed by the dynamics of heavy species (i.e. neutrals and ions), which, contrary to electrons, cannot adapt to V d changes when f s is far from f BM , is recovered, and it is only partially modulated by V d through T e dynamics.
For both modulation cases with f s = 3 and 110 kHz, the HODMD modes corresponding to f BM = 20 kHz for S e , T e , n n , n e present a similar structure and spatio-temporal characteristics as those of the natural BM of the thruster shown in figure 13. For the same plasma magnitudes, figures 15(a) and (b) show the 3 kHz and 110 kHz modes for modulated cases with f s = 3 kHz and f s = 110 kHz, respectively. In addition to their frequency, several aspects reveal that these modes are induced by the voltage modulation and do not present the features of a predator-prey type ionization instability typically associated to the BM characterizing HET operation.
First, the roughly flat axial profile of the phase of S e corresponds to global oscillations inside the thruster channel, and the progressive-wave structure travelling upwards towards the anode, which characterizes both the natural and driven BM in figure 13 is not present. Second, the spatial distribution of the magnitude of the 110 kHz n n mode reveals that the main neutral density oscillations at that frequency take place only within a very localized region close to the chamber walls. These n n oscillations are due to slow neutrals from ion recombination at the walls at f s = 110 kHz, which become ionized before reaching the center of the channel. Finally, even though n n oscillations take place in the bulk plasma for the 3 kHz mode, n e oscillations are found to lead those of n n in most of the domain, contrary to the case of a BM. This result for the f s = 3 kHz case seems to indicate that, throughout the slow modulation cycle, the growth in plasma density within the channel is limited by the T e decay and ion recombination at the walls, which later promotes a neutral density peak inside the thruster vessel.

Conclusions
The dynamic response of an MS-HET with sinusoidally modulated discharge voltage has been investigated through 2D (axial-radial) simulations of the discharge, then assessed through Fourier analysis and HODMD techniques. Parametric analyses for several modulation amplitudes and frequencies show a good qualitative agreement with experimental measurements and 1D numerical results reported in previous studies for CHT prototypes in the low power range (<1 kW), suggesting that the same fundamental phenomena drive the modulated response across different designs. The ubiquitous presence in these ExB devices, including also conventional and anode-layer-type HETs [30,31], of the BM oscillations of the discharge current, attributed to the same fundamental predator-prey type ionization instability in which the electron temperature solution inside the vessel plays an essential role, could explain this interesting result. The role of plasma-wall losses could be another reason behind the similar modulated response across these prototypes: on the one hand the surfaceto-volume ratio of the MS-HET is comparable to those of the low-power CHT and HETs of the previous modulation studies; on the other hand, the MS topology considered here (which to some extent is much more similar to that of the CHTpm2 than to that of a conventional HET) significantly limits plasma-wall interaction.
For constant modulation frequency, time-averaged values and oscillation amplitudes of I d and F increase with the amplitude of modulation and a linear and a non-linear regime are observed in the response for low and high modulation amplitude, respectively. For highly oscillating discharges, the computation of efficiency figures based on time-average quantities is shown to avoid misestimations induced by RMS values. Within the non-linear regime, a resonant-like response is found for frequencies close to the natural BM. The analysis of partial efficiencies has revealed that anode modulation yields only limited gains in thruster performance (about 2% in thrust efficiency) due to counteracting effects inherent to the coupling between plasma production and acceleration processes characterizing HET operation.
The 2D solution for relevant plasma variables reveals that, for modulation frequencies close to the natural BM, the electron temperature governs the ionization production dynamics, and is therefore the main factor responsible for the control of the dominant frequency of the modulated discharge oscillations, enabling the externally driven BM induced by anode voltage modulation. The complex spatio-temporal structure of the 2D plasma solution is analyzed through the HODMD technique, which has revealed the travelling character of the ionization production wave inside the thruster channel, common to the natural and driven BM. Ion recombination at the lateral walls of the thruster channel significantly contributes to neutral replenishment in the thruster vessel and gives rise to an apparent 1D axial neutral density wave travelling downstream with phase velocity about five to six times larger than the macroscopic axial velocity of neutrals.
For low modulation amplitudes within the linear regime, the natural BM of the thruster reappears as a dominant or co-dominant mode of the discharge current oscillations when the modulation frequency is far from the natural BM one. The transition from a driven BM to the natural one occurs when the electron temperature, which oscillates at the modulation frequency, loses control over the dynamics of the ionization production, and the natural ionization cycle of the discharge, governed by heavy species dynamics, is recovered. The main oscillations of the plasma production inside the thruster vessel then correspond to the natural BM, and secondary modes at the modulation frequency, induced by the discharge voltage through the electron temperature dynamics, do not exhibit the features of a predator-prey type ionization instability.
Discharge voltage modulation has been proven an effective technique to control the frequency of the discharge oscillations in a range close to the natural BM of the thruster without loss in performance, thus presenting interesting applications for EMI mitigation and for plasma diagnostics. The analysis of the effects of anode modulation for different operating conditions covering the nominal operation range of the thruster would permit evaluating its impact on the performance map and I-V characteristics.

Data availability statement
The data that support the findings of this study are available from the corresponding author upon reasonable request. National Research Agency), under Grant No. PID2019-108034RB-I00/AEI/10.13039/501100011033. This research has been carried out on the same configuration of the HT5k prototype of SITAEL used in [40], where numerical results were partially validated with experimental data. The present analysis is a theoretical exploratory study and no experimental validation has been undertaken so far. Therefore, the conclusions of this paper pertain exclusively to its authors.
This work is part of the PhD thesis of J Perales-Díaz, under development, and will be included as an independent chapter within it.