Theory and Simulations of Compressional and Global Alfvén Eigenmode Stability in Spherical Tokamaks

Controlled thermonuclear fusion is currently being researched as a future source of sustainable commercial energy. Fusion is the most energy dense method of power generation aside from matter-antimatter annihilation. For example, roughly seven million kilograms of oil must be combusted to produce the same amount of energy as one kilogram of fuel for fusion. By comparison, roughly six kilograms of fuel would be required for the same energy output in a modern nuclear fission power plant. The benefits of fusion energy over existing methods of power production, as well as the outstanding scientific, technological, and engineering challenges have been written about extensively and will not be repeated here (see Ref. 1 and references therein). Instead, a very brief, high level background will be given in order to place the work of this thesis into context for a less specialized audience and motivate its contribution towards this goal.

Nuclear fusion occurs when two atomic nuclei are merged into a heavier nucleus. The difference in total mass of the initial vs final particles is released as energy in accordance with Einstein’s relation $\Delta E=\Delta mc^{2}$. The nuclear binding energy curve in Fig. 1.1 describes how tightly bound each nucleus is. A nuclear reaction which starts at a state of lower binding energy and finishes with higher binding energy will release energy $(\Delta E>0)$, while the opposite requires energy to occur $(\Delta E<0)$. Since the nuclear binding energy curve peaks at an isotope of iron, light nuclei, such as isotopes of hydrogen, helium, etc. release energy when fusion occurs, whereas all isotopes heavier than iron release energy during the opposite process: nuclear fission.

Nuclear fusion occurs due to the strong nuclear force, which attracts nuclei to one another. However, this force only acts at extremely short ranges at the scale length of an atomic nucleus: $10^{-15}$ m. In order to induce fusion, two nuclei must be brought extremely close together. At larger distances, the force between charge neutral atomic nuclei is repulsive due to valence electron-electron forces (either electrostatic repulsion or chemical forces which lead to the formation of molecules but prevent overlap of the atomic nuclei). If the nuclei have instead been ionized, electrostatic repulsion exists due to the positively charged protons instead of the valence electrons. The general strategy to induce fusion is to heat ions to very high temperatures (greater than 10 million degrees Celsius) such that the the positively charged nuclei have a non-negligible chance of overcoming the Coulomb barrier and fusing.

The likelihood that a fusion reaction will occur is characterized by the reactivity, which depends both on the ion temperature and specific pair of nuclear isotopes as shown in Fig. 1.2 for a few common fusion reactions. Those most relevant to this thesis are

Here, D stands for deuterium, an isotope of hydrogen with one proton (p) and one neutron (n), and T represents tritium, an isotope of hydrogen with two neutrons. He³ and He⁴ are isotopes of helium with three and four total nucleons, respectively. He⁴ is referred to as an $\alpha$ particle for historical reasons. Since DT fusion has the largest reactivity, it is the fuel of choice for future fusion reactors. However, tritium is radioactive and poses health risks if inhaled or ingested (for instance, through environmental contamination), so its use in modern fusion experiments is uncommon due to additional costs associated with safety procedures. Instead, pure deuterium fuel is prevalent in experiments, which is not radioactive and therefore poses no health risks.^***Although DD reactions generate tritium as a fusion product, the rate of fusion in modern experiments is so low that the amount of tritium produced is negligible compared to the quantities needed for a DT reactor and consequently does not constitute a hazard.

In addition to the temperature requirements, the plasma (ions and their stripped electrons) must also be confined to a certain volume so that the ions have a chance of colliding and initiating a fusion reaction. Stars such as the sun which are powered by fusion confine their plasma gravitationally – this is not an option for terrestrial fusion reactors since it requires astronomical amounts of mass and space. While there are many methods of confinement, the one that has historically received the most research attention is toroidal magnetic confinement. In this approach, a strong helical magnetic field is established in toroidal geometry such that the charged particles will (ideally) spiral around the field lines indefinitely due to the Lorentz force. In a tokamak,^†††Tokamak is an English transliteration of a Russian acronym for “toroidal chamber with axial magnetic field.” the magnetic field is set up by a combination of external coils and the poloidal field generated by an inductively-driven toroidal current in the plasma. This configuration is shown in Fig. 1.3. Unfortunately, the plasma is not confined as well as this idealized scheme would imply. Several spontaneous mechanisms exist in reality which degrade the plasma’s confinement, ranging from macroscopic instabilities to microscopic turbulent processes. The plasma also emits electromagnetic radiation due to the acceleration of charged particles. Consequently, the plasma will “leak” energy at a certain rate, cooling the plasma and reducing the rate of fusion reactions unless constant heating power is provided to maintain its temperature.

Consider again the DT fusion reaction in Eq. 1.1c. Due to momentum conservation, 14.1 MeV of the energy released is carried by the neutron as kinetic energy and 3.5 MeV is carried by the $\alpha$ particle. Charge neutral, the neutron is unaffected by the magnetic field and rapidly leaves the confinement region. In a power plant, the neutrons will be caught in an exterior blanket, converting their kinetic energy into heat which can be used to generate electricity through conventional methods (e.g. boil water to spin a turbine that drives an electric generator via induction). In contrast, the $\alpha$ particle is charged and can be confined by the magnetic field with the rest of the plasma. Ideally, the $\alpha$ particles will transfer their energy to the background fuel before leaving the system as cold ash, replenishing the plasma’s stored energy as it is lost due to the mechanisms mentioned previously. Ignition occurs when the plasma heating by the fusion products balances the heat loss processes, without the need for auxiliary heating from external sources (to be discussed in Sec. 1.2). The goal for a future commercial fusion reactor is to reach ignition.

A remaining milestone before ignition is “scientific breakeven,” where the fusion power generated is greater than the total heating power required to sustain the reaction. It is characterized by the gain factor $Q=P_{\text{fusion}}/P_{\text{heat,ext}}$. Hence, $Q>1$ is the condition for the power generated by fusion to exceed the power required to heat the plasma to fusion temperatures (breakeven). Ignition corresponds to $Q=\infty$ since it requires no external heating – the reaction is self-sustaining. To date, the highest achieved fusion gain was $Q=0.62$ by the Joint European Torus (JET) in 1997 with DT fuel.³ Shortly thereafter, the Japanese tokamak JT-60 reached plasma conditions with DD fuel that would correspond to $Q=1.25$ if achieved with DT fuel.⁴ A massive international collaboration is currently underway to construct the ITER tokamak in southern France, which is being designed to achieve a transient peak value of $Q=10$ and sustain $Q=5$ for several minutes. Its plasma volume will be ten times larger than the largest existing tokamak, JET. Construction is scheduled to finish in 2025, with initial experiments beginning soon thereafter, and high performance DT experiments commencing around 2035.

In current experiments, the plasma is heated primarily by external sources due to the low amount of fusion power generated. Some heating is provided by the resistive dissipation of the inductively-driven plasma current. Additional heating is required to reach the high temperatures desired for fusion experiments. There are many methods to inject electromagnetic waves with external antennas in order to heat a resonant sub-population of the plasma. This resonant population subsequently heats the rest of the plasma via collisions.

Another very common method which is intimately connected to the subject of this thesis is neutral beam injection (NBI). In this scheme, large linear devices external to the tokamak chamber electrostatically accelerate ions to high energies much larger than the desired temperature of the plasma. Just before being launched into the tokamak, the ions are sent through a volume of neutral gas in order to neutralize the energetic ions through charge exchange collisions, allowing the energetic particles to penetrate into the plasma.

Once within the plasma, charge exchange will ionize most of the energetic particles near the plasma core (if they entered as ions instead of neutrals, they would be immediately confined to the edge by the magnetic field and therefore ineffectively heat the plasma). The energetic beam ions subsequently transfer their energy to the background plasma through collisions.⁵ Fusion born $\alpha$ particles heat the background plasma in a similar fashion. Fusion plasmas often have temperatures of the order $T\sim 1-10$ keV,^‡‡‡Here and elsewhere, “temperatures” quoted as energies include implicit multiplication by the Boltzmann constant $k_{B}=1.38\cdot 10^{-23}$ m² kg s^-2 K^-1, as is standard in the plasma physics community. whereas NBI involves energies of $\mathcal{E}_{beam}\sim 40-100$ keV in modern experiments. In ITER, 1 MeV beam ions and 3.5 MeV $\alpha$ particles will also be present. These ions satisfy the ordering $v_{th,i}\ll v_{\text{EP}}\ll v_{th,e}$ where $v_{th,i}$ and $v_{th,e}$ are the background (thermal) ion and electron temperatures, respectively. Hence, they are referred to as “energetic particles” (EP) or “fast ions” in order to distinguish them from the Maxwellian background.^§§§Runaway electrons can also become energetic particles in tokamaks, satisfying $v_{th,e}\ll v_{\text{EP}}$, though their dynamics are somewhat different from those of fast ions and are not studied in this thesis.

The steady state velocity distribution for a constant source of ions injected with velocity $v_{0}$ has previously been calculated analytically and is known as the slowing down distribution⁶

Here, $\theta(x)$ is the heaviside step function, defined piecewise as 1 for $x>0$ and 0 otherwise. Hence it cuts off the distribution at velocities above the injection velocity. In the denominator, $v_{c}$ is the critical velocity, defined as $m_{b}v_{c}^{2}/2=14.8A_{b}T_{e}/A_{i}^{2/3}$. Here, $A_{i}$ and $A_{b}$ are the atomic numbers of the thermal ions and neutral beam ions, respectively, and $m_{b}$ is the mass of the beam ions. When $v_{0}>v_{c}$, the injected ions will transfer most of their heat to the thermal electrons and the rest to the thermal ions, which then equilibrate through their own collisions. This is the usual case for neutral beam injection in modern experiments, and will also be true for energetic particles in ITER. The slowing down distribution is a valid description of fusion products (which are always born at a specific energy) and also NBI ions, with one caveat. Consider the case of deuterium NBI. Due to the neutralization process described in the beam generation, neutral molecular deuterium D₂ and D₃ also form in addition to neutral D atoms. Since these species all have the same energy (proportional to the beam voltage) but different masses, they enter the plasma at different energies. Hence, NBI distributions are more accurately described as a weighted sum over slowing down distributions with injection velocities $v_{0}$, $v_{0}/\sqrt{2}$, and $v_{0}/\sqrt{3}$.

Because the beams are injected into the plasma at a specific location and with a specific orientation relative to the background magnetic field, the distribution function of their fast ions will be anisotropic in velocity space. The peak value of the pitch of the NBI distribution in the core can be estimated as $v_{\parallel}/v\approx R_{\text{tan}}/R_{0}$, where $R_{\text{tan}}$ is the value of the major radius where the beam line is tangent to the magnetic field (known as the tangency radius), and $R_{0}$ is the radius of the magnetic axis (unique distance where the poloidal magnetic field vanishes). Further discussion of suitable models of the distribution function will be discussed in Sec. 2.2 and Sec. 4.2.

Although they comprise a minority of the plasma by density $(n_{\text{EP}}\ll n_{e})$, the energetic particle pressure can be comparable to the thermal plasma pressure due to their large energy even in current devices. Therefore, they require special consideration in fusion plasmas because they introduce qualitatively different physical behavior from what is expected due to the background ions and electrons. Fundamentally, energetic particles must be treated as a kinetic species. Fast ion orbit widths and Larmor radii are much larger than those of the thermal species. Importantly, the typically non-Maxwellian fast ion distributions can be a source of free energy to destabilize plasma waves that would otherwise be stable. As will be discussed in detail in Sec. 2.3, waves can be driven unstable or further stabilized by fast ions depending on the sign of gradients in the fast ion distribution. This fact was well known at the time of Mikhailovskii’s 1975 review of instabilities in a non-Maxwellian plasma.⁷ Instabilities driven by fusion products were first identified by Kolesnichenko in 1967⁸ and treated specifically for shear Alfvén waves driven by neutral beam ions in a tokamak by Rosenbluth in 1975.⁹

In some cases, fast-ion-driven instabilities can have benign effects on plasma performance. In those circumstances, they need not be avoided and can even serve as effective diagnostics for background plasma properties,¹⁰ pellet injection,^{11, 12} or the fast ion distribution.¹³ The interaction between energetic particles and plasma waves is most frequently studied within the context of fast ion transport. Fast ions can drive waves which consequently lead to fast ion redistribution or loss from the system, impacting their ability to heat the background plasma and potentially damaging the device.¹⁴ Theoretically, judicious excitation of waves can also be used to more effectively transfer the energy of fusion products to the background plasma using the $\alpha$ channeling scheme.^{15, 16, 17} The motivation for the work of this thesis, to be described in more detail in Sec. 1.6, is a situation where instabilities excited by fast ions can modify the electron temperature profile without any fast ion transport.¹⁸ In general, the study of fast-ion-driven instabilities is concerned with predictability and control. In order to avoid their deleterious effects or harness advantageous ones, we must be able to reliably predict their excitation and self-consistent interaction with the fast ions and background plasma.

Effective energy exchange (driving or damping the wave) between a wave and fast ions requires a resonant wave-particle interaction. Conceptually, resonance occurs when periodic motion of a particle is synchronized with the periodic fluctuation of a wave such that the particle experiences a net force due to the wave during each orbit. More precisely, this condition can be stated as $\oint\bm{v}\cdot\delta\bm{E}\,dt\neq 0$, where the time integration is over a periodic orbit. Integration of an analogous expression over the full fast ion distribution yields the total energy exchange. In a tokamak, resonance occurs when the following general condition is satisfied:

Here, $\left\langle\dots\right\rangle$ denotes orbit averaging over time, $\omega$ is the wave frequency, $\left\langle\omega_{\phi}\right\rangle$ is the average toroidal orbital frequency, $\left\langle\omega_{\theta}\right\rangle$ is the average poloidal orbital frequency, and $\left\langle\omega_{cs}\right\rangle$ is the average cyclotron frequency of the species of interest: $\omega_{cs}=\left|q_{s}\right|B/m_{s}$. The frequencies must be related by integers $n$, $p$, and $\ell$. The toroidal mode number of the wave defines $n=k_{\phi}R_{0}$, while $p$ can be different from the wave’s poloidal mode number $m$ due to particle drift motion. The integer $\ell$ is the cyclotron coefficient, typically equal to zero for low frequency waves and non-zero for waves in the ion cyclotron frequency range. Eq. 1.3 implies a “slow instability,” e.g. with $\gamma\ll\omega_{\text{bounce}}$, since it describes a net wave-particle interaction over a full orbit (global resonance) as opposed to an instantaneous/transient resonant interaction which would describe a “fast instability” (local resonance). Further discussion of the local vs global resonance conditions can be found in Ref. 19, Ref. 20, and Ref. 21. The resonance condition in Eq. 1.3 can be equivalently written in terms of wave vectors and velocities instead of frequencies:

Here, $k_{\parallel}$ and $k_{\perp}$ are projections of the wave vector $\bm{k}$ parallel and perpendicular to the background magnetic field, while $v_{\parallel,\text{res}}$ and $v_{\text{Dr}}$ are the same projections for the fast ion guiding center velocity. These equivalent conditions are known as the general Doppler shifted cyclotron resonance condition. Note that here and for the rest of the paper, the “Doppler shift” refers to the shift in the resonance due to a particle’s parallel and drift motion, not the bulk rotation of the plasma. The Landau resonance corresponds to $\ell=0$, the “ordinary” cyclotron resonance has $\ell=1$, and the “anomalous” cyclotron resonance has $\ell=-1$. For sub-cyclotron frequencies, and in the usual case where $\left|\left\langle k_{\parallel}v_{\parallel}\right\rangle\right|\gtrsim\left|\left\langle k_{\perp}v_{\text{Dr}}\right\rangle\right|$, counter-propagating modes $(k_{\parallel}<0)$ can only satisfy the ordinary cyclotron resonance, while co-propagating modes can interact through the Landau or anomalous cyclotron resonances, depending on their frequency. In this work, co- and cntr-propagation are defined relative to the direction of plasma current and beam injection.

With the understanding that fast ions can resonantly destabilize plasma waves, we will now introduce the specific types of waves that are the subject of this thesis.

Only a small subset of the rich ecosystem of plasmas waves is relevant to this thesis. The rest are discussed in textbooks such as Ref. 22. Of interest to this work are the waves which exist in magnetohydrodynamics (MHD). MHD is a model which describes the thermal plasma as a charged fluid. It is applicable when (1) the system is sufficiently collisional $(\lambda_{\text{mfp}}\ll L)$ such that the electrons and ions have Maxwellian distributions, (2) the plasma is macroscopically quasineutral $(\lambda_{D}\ll L)$, and (3) the plasma is non-relativistic such that wave phase velocities are small compared to the speed of light $(v_{p}=\omega/k\ll c)$.

In toroidal geometry, the MHD waves discussed above become more complicated. There are still two main branches of the dispersion: the compressional and shear waves, but the eigenfrequencies and eigenfunctions (mode structures) will depend on details of the toroidal geometry and equilibrium plasma profiles. Tokamaks host a vibrant zoo of Alfvén eigenmodes, depicted in Fig. 1.5, which may be driven unstable by phase space gradients in fast ions generated by heating in the ion cyclotron range of frequencies (ICRF), neutral beam injection, or fusion products. These waves may be excited at frequencies spanning many orders of magnitude: from EPMs and beta-induced Alfvén acoustic modes (BAAE)^{32, 33} at very low frequencies – tens of kHz, and gap modes including toroidicity-induced Alfvén eigenmodes (TAEs)^{34, 35} with $f\sim v_{A}/4\pi qR\approx 100$ kHz, to moderate frequency sub-cyclotron modes (CAEs, GAEs) $f\sim f_{ci}/2\approx 1$ MHz, and even above the cyclotron frequency (ion cyclotron emission, or ICE).^{36, 37} A number of helpful reviews exist on the topics of energetic particles and fast-ion-driven instabilities in tokamaks and related devices which can provide much more detail than there is room for in this brief thesis introduction.^{38, 39, 35, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49} There are also two textbooks^{50, 51} devoted exclusively to the derivation of Alfvén waves in various theoretical frameworks. The two instabilities addressed in depth in this thesis will now be introduced in more detail.

With the preceding physics background in mind, we will now review the specific experimental observations that motivate this work, and identify the problems that the thesis aims to address with a theoretical approach.

The National Spherical Torus Experiment (NSTX) and its recent upgrade (NSTX-U),¹⁰¹ shown in Fig. 1.8, are low aspect ratio $(R/a=1.3)$ tokamaks, also known as a “spherical tokamaks” (STs). Spherical tokamaks are currently being researched as a potentially more cost-effective route to a fusion energy pilot plant^{102, 103, 104, 105} than conventional large tokamaks $(R/a=2.5-4)$, since construction cost scales with device size. NSTX and the similar spherical tokamak MAST in the UK both found strong inverse confinement scaling with normalized electron collisionality $\nu_{e}^{*}$ ($\propto n_{e}/T_{e}^{2}$) – indicating that plasma confinement may be more favorable at high temperatures in spherical tokamaks than in conventional ones.^{106, 107, 108} This may be explained by magnetic field lines spending relatively more time in the “good curvature” region^{109, 110} and strong $\bm{E}\times\bm{B}$ rotation shear.^{111, 112, 113} Spherical tokamaks are also capable of achieving high $\beta\approx 10-40\%$, compared to $\beta\approx 3-10\%$ for conventional aspect ratio tokamaks,¹⁰² allowing more efficient plasma confinement.

NSTX has a major radius of 0.85 m, and a minor radius of 0.68 m. It achieves densities on the order of $10^{19}$ m^-3 and temperatures of approximately 1 keV. The upgrade from NSTX to NSTX-U included the installation of a new central magnetic and an additional neutral beam source. Consequently, the maximum toroidal field and plasma current were doubled from 0.5 T to 1 T and 1 MA to 2 MA, respectively, and discharge duration was increased from 1 s to 5 s.

The NSTX-U deuterium neutral beams can operate at a maximum voltage of 90 kV, with total power doubled from 6 MW to 12 MW during the upgrade. Importantly, the new neutral beam source was installed in a different geometry from the original beam sources, as shown in Fig. 1.9. Each beam has three beam lines with slightly different injection angles. The original beam lines inject more radially, with tangency radii of $R_{\text{tan}}=0.5\text{ m, }0.6\text{ m, and }0.7\text{ m}$, whereas the new beam sources lead to a more field-aligned (tangential) energetic particle population, due to tangency radii of $R_{\text{tan}}=1.1\text{ m, }1.2\text{ m, and }1.3\text{ m}$ (the magnetic axis is located near $R_{0}\approx 1.05$ m due to the Shafranov shift). Combination of these beam sources provides substantial flexibility in tailoring the fast ion distribution in phase space in order to study fast-ion-driven instabilities.

The primary diagnostics used to measure CAEs and GAEs on NSTX(-U) are magnetic Mirnov coils⁵² and a reflectometer array.^{114, 68, 95} The 10 Mirnov coils are toroidally spaced non-uniformly, with sufficient separation to detect $|n|\leq 18$, and detect frequencies up to 5 MHz ($\approx 1-2f_{ci}$ in NSTX). They provide very precise measurements of the unstable mode frequencies and their time evolution throughout the discharge, as well as a measurement of the edge polarization of the fluctuations. The reflectometer array measures plasma displacements at 16 locations, bunched mostly in the edge/pedestal region, with a few points extending to near the magnetic axis. These internal measurements can be used to calculate radial mode structure of the modes across most of the plasma minor radius.⁹⁵ Both diagnostics can be used to measure the mode amplitudes.

Previous comparison of the frequencies of the observed modes in the NSTX H-mode discharge # 141398 and the most unstable modes in HYM modeling show a close match in frequency for each toroidal harmonic.⁶³ Comparison of the mode structures inferred from reflectometry measurements and present in HYM simulations show qualitative similarities. Preliminary analysis of high-$k$ scattering measurements^{115, 116} has shown possible signatures of CAE to KAW mode conversion in NSTX,¹¹⁷ though further analysis is needed for conclusive identification.

NSTX-U provides an excellent laboratory for studying energetic-particle-driven instabilities. Due to its low magnetic field and high beam power, it is capable of operating across a wide range of $v_{\text{EP}}/v_{A}$ and $\beta_{\text{EP}}/\beta_{\text{tot}}$, which are key parameters controlling the activity of these instabilities.¹¹⁸ Moreover, dimensionless parameters for the fast ions in NSTX-U generated from neutral beam injection can be comparable to those of fusion alphas in large burning plasmas such as ITER, connecting these studies to future burning reactors. The focus of this thesis will be advancing the theoretical understanding of the global and compressional Alfvén eigenmode stability, motivated by two significant experimental observations anomalous electron temperature profile flattening in NSTX and robust GAE suppression with the off-axis beam sources in NTSX-U.

The primary motivation of this thesis is the anomalously flat electron temperature profile observed in NSTX during H-mode neutral-beam-heated discharges, which has been linked to the presence of high frequency Alfvén activity, identified as a mixture of CAEs and GAEs.^{66, 90, 91, 67, 52, 68, 95, 96, 69} While the ion temperature matches predictions from the global transport codes TRANSP,¹¹⁹ indicating neoclassical ion confinement, the electron temperature deviates from these descriptions at high beam power. As the beam power is increased, the electron temperature radially broadens while the central electron temperature stagnates or even decreases. The flattening occurs in both L-mode and H-mode discharges, most dramatically at higher beam power, presenting opportunity for its further investigation on NSTX-U with its doubled capacity of neutral beam heating to a maximum of 12 MW. The observed broadening of the electron temperature profiles contrasts with the expectation that increasing beam power will result in an increased on-axis temperature with negligible impact on its radial profile. Anomalously low electron temperature limits fusion performance (when not operating in a hot ion mode), and could imperil future spherical tokamak development.

The inferred electron diffusion profile required to explain this observation is unusual and quite large in magnitude. Moreover, since the temperature profile is so flat, it is unlikely that microturbulence could be the source of this anomalous diffusion. Local gyrokinetic simulations of turbulence in the “core flat region” ($r/a\leq 0.15$ – where gradients are absent) do not reproduce this level of anomalous transport.^{112, 120} Conversely, high frequency Alfvén activity is often observed in discharges with the anomalous electron temperature flattening. This was demonstrated by Stutman et al.in an experiment where an H-mode plasma was reproduced with 2 MW, 4 MW, and 6 MW of injected beam power in successive discharges.¹⁸ As the beam power increased, the amplitude and number of distinct high frequency modes increased in conjunction with flattening of the temperature profile and an increase in the peak electron diffusivity by an order of magnitude, as depicted in Fig. 1.10. In addition to dedicated experiments designed to investigate this phenomenon, a detailed database of shots with substantial Alfvénic mode activity has been compiled by Fredrickson¹¹⁸ and extended by Tang¹²¹ to catalog details about CAEs and GAEs. Statistical analysis of this database further strengthens the link between between strong CAE/GAE activity and $T_{e}$ flattening. For these reasons, the CAEs and GAEs are suspected to be the primary cause of the unexplained electron energy transport.

In response to observations of anomalously large ratios of $T_{i}/T_{e}$ in early NSTX operations, it was also suggested that the simultaneous excitation of many CAEs could generate broadband turbulence, leading to efficient heating of the thermal ions through stochastic diffusion.^{122, 123} In that case, further experimental analysis ruled this out as a viable mechanism for NSTX conditions based on the observed mode amplitudes.¹²⁴ Moreover, recent analytic studies by Kolesnichenko suggest that the CAE/GAE-induced energy transport could explain anomalously high confinement observed in certain JET discharges¹²⁵ –- introducing the possibility of harnessing these mechanisms for improved performance. Predictive capabilities and control of the $T_{e}$ profile with neutral beam heating is vital to achieving and exploring high performance plasmas in STs such as NSTX-U and future ST-FNSF designs. In order to probe the aforementioned favorable confinement scaling in regimes that can not be accessed by other devices, STs must be able to achieve their target $T_{e}$ profiles, which is jeopardized by their spontaneously flattening in the presence of CAEs and GAEs. The importance of this topic is emphasized in the NSTX-U 5 year plan:¹²⁶ “The successful development and implementation of an energetic particle model for electron thermal transport is essential to achieve the high priority goal of $T_{e}$ and $T_{i}$ profile predictions.”

Two theoretical mechanisms have been previously proposed to explain how the CAEs and GAEs could modify the electron temperature profile. The first involves the stochastization of electron orbits induced by the presence of many modes of sufficient amplitude. Test particle simulations using the guiding center particle code ORBIT ¹²⁷ have shown that there is a sharp increase of two orders of magnitude in the electron diffusion coefficient if a threshold in the number of unstable GAEs is surpassed, provided that the mode amplitudes are sufficiently large.¹²⁸ It is likely that a qualitatively similar effect occurs for sufficiently many CAEs, though this has not been confirmed. A strong dependence on the mode amplitude was also reported, as well as sensitivity to $\delta E_{\parallel}$ fluctuations. The second process is an energy-channeling mechanism where a core-localized CAE or GAE mode converts to a KAW at the Alfvén resonance location. Since the KAW efficiently damps on thermal electrons, this process modifies the effective beam energy deposition profile, redirecting neutral beam power from the core to the edge through the Poynting flux. This possibility has been studied analytically by Kolesnichenko et al.in the case of GAEs^{129, 130, 131, 132} and has been confirmed numerically for CAEs by Belova et al.^{133, 63}

Both of these mechanisms – energy channeling and orbit stochastization – have been shown to have an effect in numerical calculations. However, there is a quantitative gap between the levels of effective transport they predict when mode amplitudes are scaled to experimental magnitudes and the amount necessary to explain the experimental $T_{e}$ profile flattening. While this phenomena is currently unique to NSTX, it could potentially be relevant to ITER where both the neutral beam ions and fusion alphas will be super-Alfvénic and hence have the potential to excite CAEs/GAEs.⁹⁷ Solving this problem can be decomposed into two parts: 1) for a given plasma and beam scenario, which modes will be unstable? 2) for a given spectrum of modes, how will they affect the electron energy transport? My thesis focuses primarily on this first part, aiming to improve understanding of the CAE/GAE stability properties.

The second experimental observation motivating the further theoretical development of CAE/GAE stability properties is the efficient suppression of GAEs with off-axis beam injection, recently discovered in NSTX-U. Early NSTX-U operations found that the original NSTX neutral beam sources excited a broad spectrum of high frequency Alfvén eigenmode activity (mostly GAEs), just as they did in NSTX. When additional beam power was supplied by the new, off-axis beam sources, all instabilities in this frequency range rapidly vanished.⁸⁹ Three examples are shown in Fig. 1.11, though this was an ubiquitous phenomenon present in at least 100 unique discharge time windows.⁹⁶ At first glance, this is a surprising finding. Typically, higher beam power will further destabilize Alfvén eigenmodes since it supplies the system with more energetic particles. Conversely, suppressing instabilities while increasing the plasma heating is the best case scenario.

Subsequent HYM modeling of an NSTX-U discharge reproduced both the excitation of GAEs with the original beam sources and also their complete stabilization with the addition of the new beam source.¹³⁴ Further simulations found that the GAE suppression can be achieved with only 7% of the total beam ions being supplied by the new beams, much lower than the 25% of total beam ions supplied in the modeled experimental discharge. A more complete theoretical understanding of this very efficient stabilization mechanism could contribute to the development of additional phase space engineering techniques for control of fast ion instabilities.^{16, 135, 136} In particular, theoretical advancements enabling the control of CAEs/GAEs will facilitate the investigation of their role in the anomalous electron energy transport.

The main goals of this thesis are to advance the theoretical understanding of CAE and GAE stability properties in application to the anomalous electron temperature profile flattening that they are associated with. Both analytic theory and numerical simulations are employed towards this goal. Each chapter in this thesis is written to be mostly self-contained, so there is an intentional degree of redundancy in some areas in order to remind readers of previous results and relevant background for each section. Appendices appear immediately after each chapter. The thesis is outlined as follows.

In Chapter 2, analytic conditions for net fast ion drive are derived for beam-driven, sub-cyclotron CAEs and GAEs. Both co- and counter-propagating modes are investigated, driven by the ordinary and anomalous Doppler-shifted cyclotron resonance with fast ions. Whereas prior results were restricted to vanishingly narrow distributions in velocity space, broad parameter regimes are identified in this work which enable an analytic treatment for realistic fast ion distributions generated by neutral beam injection. The simple, approximate conditions derived in these regimes for beam distributions of realistic width compare well to the numerical evaluation of the full analytic expressions for fast ion drive. Moreover, previous results in the very narrow beam case are corrected and generalized to retain all terms in $\omega/\omega_{ci}$ and $\left|k_{\parallel}/k_{\perp}\right|$, which are often assumed to be small parameters but can significantly modify the conditions of drive and damping when they are non-negligible. Favorable agreement is demonstrated between the approximate stability criterion, simulation results, and a large database of NSTX observations of cntr-GAEs.

In Chapter 3, a similar analytic approach is taken for co-propagating CAEs and GAEs driven by the Landau resonance. Approximations applicable to realistic neutral beam distributions and mode characteristics observed in spherical tokamaks enable the derivation of marginal stability conditions for these modes. Such conditions successfully reproduce the stability boundaries found from numerical integration of the exact expression for local fast ion drive/damping. Coupling between the CAE and GAE branches of the dispersion due to finite $\omega/\omega_{ci}$ and $\left|k_{\parallel}/k_{\perp}\right|$ is retained and found to be responsible for the existence of the GAE instability via this resonance. Encouraging agreement is demonstrated between the approximate stability criterion, simulation results, and a database of NSTX observations of co-CAEs.

In Chapter 4, a comprehensive numerical study is presented in order to investigate CAE/GAE stability properties for a wide range of beam parameters in realistic NSTX conditions. Linear simulations are performed with the hybrid MHD-kinetic initial value code HYM in order to capture the general Doppler-shifted cyclotron resonance that drives the modes. The simulations reveal that unstable GAEs are more ubiquitous than unstable CAEs, consistent with experimental observations, as they are excited at lower beam energies and generally have larger growth rates. The local analytic theory derived in Chapter 2 and Chapter 3 is used to explain key features of the simulation results, including the preferential excitation of different modes based on beam injection geometry and the growth rate dependence on the beam injection velocity, critical velocity, and degree of velocity space anisotropy. Drive due to velocity space anisotropy is capable of explaining most trends theoretically, though it is found that gradients with respect to $p_{\phi}$ can be responsible for a substantial fraction of the fast ion drive for co-propagating modes. The background damping rate is inferred from simulations and estimated analytically for relevant sources not present in the simulation model, indicating that co-CAEs are closer to marginal stability than modes driven by the cyclotron resonances.

In Chapter 5, the numerical discovery of strong energetic particle modifications to GAEs in HYM simulations of NSTX-like plasmas is presented and investigated. Key parameters defining the fast ion distribution function – the injection velocity and injection geometry – are varied in order to study their influence on the characteristics of the excited modes. It is found that the frequency of the most unstable mode changes significantly and continuously with beam parameters, in accordance with the Doppler-shifted cyclotron resonances which drive the modes, and depending most substantially on the injection velocity. This unexpected result is present for both counter-propagating GAEs, which are routinely excited in NSTX, and high frequency co-GAEs, which have not been previously studied. Large changes in frequency without clear corresponding changes in mode structure are signatures of an energetic particle mode, referred to here as an energetic-particle-modified GAE (EP-GAE). Additional simulations conducted for a fixed MHD equilibrium demonstrate that the GAE frequency shift cannot be explained by the equilibrium changes due to energetic particle effects.

Lastly, a summary of the main results and discussion of future research directions is given in Chapter 6.

The analysis of this chapter focuses on fast ions interacting with CAEs/GAEs through the ordinary or anomalous cyclotron resonances. Drive/damping due to the Landau resonance is treated in Chapter 3. General expressions for the growth rate of these instabilities were originally derived for mono-energetic beam^{137, 138} and bi-Maxwellian¹³⁹ distributions, as well as for an arbitrary distribution⁷ in a uniform plasma. These derivations were later extended and applied to NBI-driven CAEs/GAEs in various experimental conditions dating back to the TFTR era^{19, 140} and continuing in more recent years with applications to JET¹⁴¹ and NSTX.^{90, 92} The recent studies on NBI-driven modes had two key limitations. First, they did not correctly treat the cutoff at the injection energy, an approach suitable for shifted Maxwellians generated by heating in the ion cyclotron range of frequencies (ICRF), but not for slowing down distributions from NBI. Second, they assumed a delta function in pitch for tractability, which is unrealistic considering the more broad distributions present in experiments, as inferred from Monte Carlo codes such as the NUBEAM ¹⁴² module in TRANSP.¹¹⁹ Prior studies also assume $k_{\parallel}\ll k_{\perp}$ and $\omega\ll\omega_{ci}$ as simplifying approximations, whereas the modes excited in spherical tokamaks such as NSTX may have frequencies approaching $\omega\lesssim\omega_{ci}$ and $k_{\parallel}\sim k_{\perp}$.

Choice of parameter regimes to study has been informed by prior and ongoing numerical modeling of CAEs/GAEs with the 3D hybrid MHD-kinetic initial value code HYM.^{63, 93, 134} The simulation model couples a single fluid thermal plasma to a minority species of full orbit kinetic beam ions and also includes the contributions of the large beam current to the equilibrium self-consistently.¹⁴³

The derivation presented in this chapter corrects and builds on prior work by providing a local expression for the fast ion drive due to an anisotropic beam-like distribution interacting via the ordinary and anomalous cyclotron resonances. The effect of finite injection energy of NBI distributions is included consistently, yielding a previously overlooked instability regime. Terms to all order in $\omega/\omega_{ci}$ and $\left|k_{\parallel}/k_{\perp}\right|$ are kept for applicability to the entire possible spectrum of modes. As in previous works, full finite Larmor radius (FLR) terms are also retained. The analytic expression can be integrated numerically for any chosen parameters in order to determine if the full fast ion distribution is net driving or damping. More interestingly, it is found that when the beam is sufficiently wide in velocity space, such as realistic distributions resulting from NBI, the integral can be evaluated approximately in terms of elementary functions, yielding compact conditions for net fast ion drive/damping that depend only on a small set of parameters describing the fast ion and mode parameters. Such expressions grant new insights into the spectrum of CAEs and GAEs that may be excited by a given fast ion distribution, as well as providing intuition for interpreting experimental observations and simulation results. Since damping sources such as electron Landau and continuum damping are not addressed in this chapter, the net fast ion drive conditions derived here should be considered as necessary but not sufficient conditions for instability.

The chapter is structured as follows. The dispersion relations, resonance condition, and model fast ion distribution function used in this chapter are described in Sec. 2.2. In Sec. 2.3, the local analytic expression for the CAE and GAE growth rates is adapted from Ref. 7 and applied to the fast ion distribution of interest. Approximations are applied to this expression in Sec. 2.4 in order to derive useful instability criteria for the cases of a very narrow beam width in velocity space (Sec. 2.4.1) and a beam with realistic width (Sec. 2.4.2) when FLR effects are small (Sec. 2.4.2.1) and large (Sec. 2.4.2.2). The derived conditions are also compared against the numerically calculated growth rates for realistic parameter values in Sec. 2.4. In Sec. 2.5, the dependence of the fast ion drive/damping on the mode properties ($\omega/\omega_{ci}$ and $\left|k_{\parallel}/k_{\perp}\right|$) is presented and compared against conclusions drawn from the approximate stability boundaries. A comparison of the approximate stability conditions against a database of cntr-GAE activity in NSTX and simulation results is shown in Sec. 2.6. Lastly, a summary of the main results and discussion of their significance is given in Sec. 2.7. The majority of the content of this chapter has been peer-reviewed and published in Ref. 144.

One goal of this work is to extend previous derivations to include finite $\omega/\omega_{ci}$ and $\left|k_{\parallel}/k_{\perp}\right|$ effects in the stability calculation, since experimental observations and modeling of NSTX suggests that these quantities may not always be small. Experimental observations often show CAEs with frequencies from $\omega/\omega_{ci}=0.3$ to exceeding the cyclotron frequency. GAEs are observed with somewhat lower frequencies of $\omega/\omega_{ci}\approx 0.1-0.5$. While $k_{\perp}$ can not be measured accurately on NSTX due to limited poloidal coil resolution, it can be calculated for the most unstable modes excited in simulations,⁶³ which show that $\left|k_{\parallel}/k_{\perp}\right|\approx 1$ is not uncommon, and can even reach $\left|k_{\parallel}/k_{\perp}\right|>3$ in some cases. This motivates using the full, unsimplified dispersion relations in uniform geometry when numerically calculating the growth rate, instead of using the common $\omega/\omega_{ci}\ll 1$ and $\left|k_{\parallel}/k_{\perp}\right|\ll 1$ assumptions found in previous works. The more complicated eigenmode equations in nonuniform toroidal systems^{77, 78, 53, 140, 19} have been derived in the past but are too complicated for our purposes.

Define $\bar{\omega}=\omega/\omega_{ci0}$, $N=kv_{A}/\omega$, $A=(1-\bar{\omega}^{2})^{-1}$, and also $F^{2}=k_{\parallel}^{2}/k^{2}$, $G=1+F^{2}$. Here, $\omega_{ci0}$ is the on-axis ion cyclotron frequency. Then in uniform geometry, the local dispersion in the MHD limits of $E_{\parallel}\ll E_{\perp}$ and $\omega\ll\left|\omega_{ce}\right|,\omega_{pe}$ is readily given by¹⁴⁵

The “$-$” solution corresponds to the compressional Alfvén wave (CAW), while the “$+$” solution corresponds to the shear Alfvén wave (SAW). The coupled dispersion in Eq. 2.1 will be used in the full analytic expression for fast ion drive. Notably, it can modify the polarization of the two modes, which in turn changes how the finite Larmor radius (FLR) effects from the fast ions contribute to the growth rate (see Eq. 2.19). Its low frequency approximations are $\omega\approx kv_{A}$ for CAWs and $\omega\approx\left|k_{\parallel}\right|v_{A}$ for SAWs. Throughout this chapter, CAW/CAE and SAW/GAE will be used interchangeably, where CAW and SAW formally refer to the solutions in a uniform slab, while CAE and GAE refer to their analogues in nonuniform and bounded geometries. Net energy transfer between a mode and the fast ions requires a sub-population of particles obeying the Doppler-shifted cyclotron resonance.

Here, $\left\langle\dots\right\rangle$ denotes poloidal orbit averaging and $\ell$ is an integer cyclotron resonance coefficient. Two resonances are studied in detail in this chapter for the sub-cyclotron modes: the $\ell=1$ ordinary cyclotron resonance and $\ell=-1$ anomalous cyclotron resonance. Orbit averaging in Eq. 2.2 is required to satisfy the global resonance condition, as opposed to the local resonance, which describes a net synchronization condition between the wave and particle on average over its orbit, even while not being in constant resonance at all points in time. This resonance condition is applicable so long as the growth rate of the mode is sufficiently smaller than the inverse particle transit time, which is satisfied by these modes according to HYM simulations.

In this chapter, we will make the approximation of $\left|k_{\perp}v_{\text{Dr}}\right|\ll\left|k_{\parallel}v_{\parallel}\right|$. Consequently, when $\omega<\omega_{ci}$ and $\left\langle v_{\parallel}\right\rangle>0$ (co-injection), Eq. 2.2 can only be satisfied for $\ell=1$ if $k_{\parallel}<0$ (mode propagates counter to the fast ions). Likewise, $\ell=-1$ requires $k_{\parallel}>0$, corresponding to co-propagation. Due to periodicity, the drift term can be approximated for passing particles²⁰ as $\left\langle k_{\perp}v_{\text{Dr}}\right\rangle\approx s\left\langle v_{\parallel}\right\rangle/qR$ for integer $s$, though this term yields relatively small corrections due to the large values of $\left|k_{\parallel}\right|$ relevant to these modes. In this approximation, the resonance condition can be rewritten as $\omega-k_{\parallel,s}v_{\parallel,\text{res}}=\ell\left\langle\omega_{ci}\right\rangle$ with $k_{\parallel,s}=k_{\parallel}+s/qR$. Conversely, for trapped particles the drift term can be approximated as¹⁴⁶ $\left\langle k_{\perp}v_{\text{Dr}}\right\rangle\approx s\omega_{b}$. HYM simulations indicate that the $s=\pm 1$ sidebands are usually more relevant than larger $\left|s\right|$.⁶³ For quantitatively accurate growth rates, all sidebands should be summed over, as done in Ref. 19 in the limit of $\omega\gtrsim\omega_{ci}\gg\omega_{b}$, and also in Ref. 92. Practically, these procedures require complicated non-local calculations which would preclude analytic progress except in extraordinarily special cases, contrary to the purpose of this chapter, which is to derive broadly applicable instability conditions. To this end, only the primary resonance $(s=0)$ will be kept when deriving approximate stability boundaries in Sec. 2.4.

Combination of the resonance condition with approximate dispersion relations can yield relations that will be useful later on. Introduce $\left\langle\bar{\omega}_{ci}\right\rangle\equiv\left\langle\omega_{ci}\right\rangle/\omega_{ci0}$ as the average cyclotron frequency of the resonant particles, normalized to the on-axis cyclotron frequency $\omega_{ci0}$. This value is approximately 0.9, as inferred from inspection of the resonant particles in relevant HYM simulations. Then defining $v_{\parallel,\text{res}}\equiv\left\langle v_{\parallel}\right\rangle>0$ (treating co-injected particles only) and rearranging Eq. 2.2 gives

The stability calculation will be applied to a slowing down, beam-like background distribution of fast ions, motivated by theory and NUBEAM modeling of NSTX discharges.⁶³ In order to satisfy the steady state Vlasov equation, the distribution is written as a function of constants of motion $v=\sqrt{2\mathcal{E}/m_{i}}$ and $\lambda=\mu B_{0}/\mathcal{E}$ in separable form: $f_{0}(v,\lambda)=C_{f}n_{b}f_{1}(v)f_{2}(\lambda)$, defined below

The constant $C_{f}$ is for normalization. The first component $f_{1}(v)$ is a slowing down function in energy with a cutoff at the injection energy $v_{0}$ and a critical velocity $v_{c}$. The cutoff at $v=v_{0}$ is contained within $f_{\text{tail}}(v;v_{0})$, which is in general a function which rapidly goes to zero for $v>v_{0}$. For ease of calculation, this is assumed to be a step function. The second component $f_{2}(\lambda)$ is a Gaussian distribution centered on some central value $\lambda_{0}$ with width $\Delta\lambda$. The variable $\lambda$ is a trapping parameter. To lowest order in $\mu\approx\mu_{0}$, it can be re-written as $\lambda=(v_{\perp}^{2}/v^{2})(\omega_{ci0}/\omega_{ci})$. Then, assuming a tokamak-like field $B\approx B_{0}/(1+\epsilon\cos\theta)$ for $\epsilon=r/R$, passing particles will have $0<\lambda<1-\epsilon$ and trapped particles will have $1-\epsilon<\lambda<1+\epsilon$. Loosely, smaller $\lambda$ means the particle’s velocity is more field aligned, such that $\lambda$ is a complementary variable to a particle’s pitch $v_{\parallel}/v$. For analytic tractability, $\lambda_{0}$ and $\Delta\lambda$ are treated as constants in this model, ignoring any velocity dependence of these parameters which may be present, especially broadening in $\lambda$ at lower energies due to pitch angle scattering. The dependence on $p_{\phi}$, is neglected in this study for simplicity, as it is expected to be less relevant for the high frequencies of interest for these modes. The model distribution does not include the two additional energy components that are present due to molecular deuterium production in the neutral beam source, as these have a quantitative but not qualitative impact on the analysis. Such effects can be recovered by summing over three beam distributions (with injection velocities $v_{0}$, $v_{0}/\sqrt{2}$, and $v_{0}/\sqrt{3}$) with appropriate weights. Comparison between the model distribution used in this study and those calculated with the Monte Carlo code NUBEAM for NSTX and NSTX-U can be found in Fig. 5 of Ref. 143 and Fig. 4 of Ref. 134, respectively.

The NSTX operating space spanned a range of normalized injection velocity $v_{0}/v_{A}=2-6$, depending on the beam voltage (typically $60-90$ keV at $2-6$ MW) and field strength ($0.25-0.50$ T) for each discharge. The beam injection geometry $\lambda_{0}$ and width in velocity space $\Delta\lambda$ are mostly determined by the neutral beam’s geometry and collimation, yielding typical $\lambda_{0}=0.5-0.7$ and $\Delta\lambda=0.3$. For this study, $v_{c}=v_{0}/2$ is used as a characteristic value. The new beam line on NSTX-U has much more tangential injection, with $\lambda_{0}\approx 0$, and also lower $v_{0}/v_{A}=1-3$ due to higher nominal field strength. A comparison between the model fast ion distribution used in this chapter (Eq. 2.6) and a NUBEAM calculation for the well-studied H-mode discharge $\#141398$, can be found in Fig. 5 of Ref. 143.

In this section, the fast ion drive/damping is derived perturbatively in the local approximation for a two component plasma comprised of a cold bulk plasma and a minority hot ion kinetic population, and applied to the anisotropic beam distribution of interest. The formula presented here extends the results obtained in Ref. 90, 92, which focused on $\omega\ll\omega_{ci}$, $k_{\parallel}\ll k_{\perp}$, and also did not study high frequency co-propagating modes ($\ell=-1$ cyclotron resonance coefficient). In contrast, the following derivation is appropriate for all values of $\omega/\omega_{ci}$ and $\left|k_{\parallel}/k_{\perp}\right|$, which is important since mode frequencies can be on the order $\omega/\omega_{ci}\sim 0.5$ or larger, and in contrast to the common large tokamak assumption, $\left|k_{\parallel}/k_{\perp}\right|$ can be of order unity, as inferred from simulations.⁹³