PSFC/JA-23-52 NSTX-U research advancing the physics of spherical tokamak J W Berkery¹, P O Adebayo-Ige², H Al Khawaldeh³, G Avdeeva⁴, S-G Baek⁵, S Banerjee¹, K Barada⁶, D J Battaglia⁷, R E Bell¹, E Belli⁴, E V Belova¹, N Bertelli¹, N Bisai⁸, P T Bonoli⁵, M D Boyer⁷, J Butt⁹, J Candy⁴, C S Chang¹, C F Clauser⁵, L D Corona Rivera¹, M Curie⁹, P C de Vries¹⁰, R Diab⁵, A Diallo¹, J Dominski¹, V N Duarte¹, E D Emdee¹, N M Ferraro¹, R Fitzpatrick¹¹, E L Foley¹², E Fredrickson¹, M E Galante¹², K F Gan², S Gerhardt¹, R Goldston¹, W Guttenfelder¹,¹³, R Hager¹, M O Hanson⁷, S C Jardin¹, T G Jenkins¹⁴, S M Kaye¹, A Khodak¹, J Kinsey¹⁵, A Kleiner¹, E Kolemen⁹, S Ku¹, M Lampert¹, B Leard³, B P LeBlanc¹, J B Lestz⁴, F M Levinton¹², C Liu¹, T Looby⁷, R Lunsford¹, T Macwan⁶, R Maingi¹, J McClenaghan⁴, J E Menard¹, S Munaretto¹, M Ono¹, A Pajares⁴, J Parisi¹, J-K Park¹,¹⁶, M S Parsons¹, B S Patel¹⁷, Y V Petrov¹⁵, M Podestà¹, F Poli¹, M Porcelli⁹, T Rafiq³, S A Sabbagh¹⁸, Á Sánchez Villar¹, E Schuster³, J Schwartz¹,⁸, A Sharma¹, S Shiraiwa¹, P Sinha¹, D Smith¹⁹, S Smith⁴, V A Soukhanovskii²⁰, G Staebler²¹, E Startsev¹, B Stratton¹, K E Thome⁴, W Tierens²¹, M Tobin¹⁸, I U Uzun-Kaymak¹², B Van Compernolle⁴, J Wai⁷, W Wang¹, W Wehner⁴, A Welander⁴, J Yang¹, V Zamkovska¹⁸, X Zhang²², X L Zhu²³, S Zweben¹ ¹ Princeton Plasma Physics Laboratory, Princeton, NJ, United States of America ² University of Tennessee, Knoxville, TN, United States of America ³ Lehigh University, Bethlehem, PA, United States of America ⁴ General Atomics, San Diego, CA, United States of America ⁵ Massachusetts Institute of Technology, Cambridge, MA, United States of America ⁶ University of California at Los Angeles, Los Angeles, CA, United States of America ⁷ Commonwealth Fusion Systems, Devens, MA, United States of America ⁸ Institute for Plasma Research, Bhat, Gandhinagar, India ⁹ Princeton University, Princeton, NJ, United States of America ¹⁰ ITER Organization, St. Paul Lez Durance, France ¹¹ University of Texas at Austin, Austin, TX, United States of America ¹² Nova Photonics, Princeton, NJ, United States of America ¹³ Type One Energy, Madison, WI, United States of America ¹⁴ Tech-X Corporation, Boulder, CO, United States of America ¹⁵ CompX, Del Mar, CA, United States of America ¹⁶ Seoul National University, Seoul, Republic of Korea ¹⁷ United Kingdom Atomic Energy Authority, Abingdon, United Kingdom ¹⁸ Columbia University, New York, NY, United States of America ¹⁹ University of Wisconsin, Madison, WI, United States of America ²⁰ Lawrence Livermore National Laboratory, Livermore, CA, United States of America ²¹ Oak Ridge National Laboratory, Oak Ridge, TN, United States of America ²² Tokamak Energy Ltd., Milton Park, United Kingdom ²³ Dalian University of Technology, Dalian, China November 2023 Plasma Science and Fusion Center Massachusetts Institute of Technology Cambridge MA 02139 USA This work was supported by the U.S. Department of Energy under contract numbers DE-AC02-09CH11466 (PPPL), DE-SC0008309 (UT, Knoxville), DESC0013977 (Lehigh), DE-SC0021385 (Lehigh), DE-SC0021113 (GA), DE-FG02-95ER54309 (GA), DE-FG02-91ER54109 (MIT), DE-SC0021156 (UT, Austin), DESC0021625 (Nova Photonics), DE-SC0021311 (Columbia), and DE-AC52-07NA27344 (LLNL). Additionally: EP/R034737/1 (UKAEA). A part of this research used resources of the National Energy Research Scientific Computing Center (NERSC), which is supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231. A part of this work was supported by the Department of Energy early career research program. Reproduction, translation, publication, use and disposal, in whole or in part, by or for the United States government is permitted. Submitted to Nuclear Fusion E-mail: [email protected] submitted November 22, 2023 Abstract. The objectives of NSTX-U research are to reinforce the advantages of STs while addressing the challenges. To extend confinement physics of low-A, high beta plasmas to lower collisionality levels, understanding of the transport mechanisms that set confinement performance and pedestal profiles is being advanced through gyrokinetic simulations, reduced model development, and comparison to NSTX experiment, as well as improved simulation of RF heating. To develop stable non-inductive scenarios needed for steady-state operation, various performance-limiting modes of instability were studied, including MHD, tearing modes, and energetic particle instabilities. Predictive tools were developed, covering disruptions, runaway electrons, equilibrium reconstruction, and control tools. To develop power and particle handling techniques to optimize plasma exhaust in high performance scenarios, innovative lithium-based solutions are being developed to handle the very high heat flux levels that the increased heating power and compact geometry of NSTX-U will produce, and will be seen in future STs. Predictive capabilities accounting for plasma phenomena, like edge harmonic oscillations, ELMs, and blobs, are being tested and improved. In these ways, NSTX-U researchers are advancing the physics understanding of ST plasmas to maximize the benefit that will be gained from further NSTX-U experiments and to increase confidence in projections to future devices. Keywords: Spherical tokamak, magnetic confinement fusion, NSTX, NSTX-U Submitted to: Nucl. Fusion Spherical tokamak (ST) concepts are currently being designed for fusion pilot plants [1], both in the public and private sectors. The National Spherical Torus Experiment (NSTX [2]) has historically provided much of the physics basis for the ST concept for fusion energy production, and it will continue to do so when the upgraded device (NSTX-U [3]) returns to operation. In the interim, NSTX-U researchers are advancing the physics understanding of ST plasmas to maximize the benefit that will be gained from further NSTX-U experiments and to increase confidence in projections to future devices [4, 5, 3, 6, 7]. STs have certain advantages: assuming component costs scale with size, their more compact size means they can provide net electricity more economically [8]. Their low aspect ratio, A, improves stability with favorable average curvature, and the high elongation plasmas that are possible also improve confinement and stability, enabling high β (the ratio of plasma pressure to magnetic pressure). There are challenges as well: managing the high heat flux, and start-up and sustainment of the plasma without space in the center column for an induction coil. The objectives of NSTX-U research are to reinforce the advantages while addressing the challenges: (i) to extend confinement physics of low-A, high β plasmas to the lower collisionality levels relevant to burning plasma regimes, (ii) to develop stable, low-disruptivity, large bootstrap fraction, non-inductive scenarios needed for steady-state operation, and (iii) to develop power and particle handling techniques to optimize plasma exhaust in high performance scenarios. The mission need for NSTX-U remains strong in an age of increased private interest in STs. When NSTX-U returns to operation, it will have the capabilities 2 MA of plasma current, 1 T toroidal field, 15 MW of neutral beam heating and 6 MW of high harmonic fast wave heating, to produce plasmas with high normalized beta (> 5), high bootstrap current fraction (> 0.7), and low collisionality (ν*e < 0.01). NSTX-U will be the best device to assess the collisionality scaling of confinement in an auxiliary heated ST, to pursue 100% non-inductive scenarios, and with an energetic particle phase space overlapping with fusion pilot plants. The outline of the paper is as follows. Section 2 discusses progress in the first objective listed above, including the topics of pedestal structure, core transport, impurity transport, and RF heating. The second objective is considered in Sec. 3, covering stability, disruptions, and machine protection. Finally, in Sec. 4, objective three is reviewed, with lithium research, divertor heat flux, and plasma blobs being the main elements. Understanding of the transport mechanisms that set confinement performance and pedestal profiles is being advanced through gyrokinetic simulations, reduced model development, comparison to NSTX experiment (and forthcoming, lower collisionality NSTX-U experiments), as well as assessment of transport of energetic particles from neutral beam injection (NBI) or range of frequency (RF) heating. 2.1. Pedestal structure 2.1.1. Progress towards pedestal pressure structure model Pedestal performance is a significant source of uncertainty in ST reactor design, and the pedestal properties of high confinement (H-mode) plasmas in NSTX have been seen to deviate from the standard model whereby ideal magnetohydrodynamic (MHD) peeling-ballooning macrostability limits the pedestal pressure and kinetic ballooning mode (KBM) microstability limits the pressure gradient. A long-standing problem for ST pedestal stability prediction has been the reliable modeling of peeling-ballooning (P-B) stability boundaries. Unlike in large aspect ratio devices, ideal P-B modes are often predicted stable for ST discharges with unstable edge localized modes (ELMs). In simulations with the state of the art extended-MHD nonlinear M3D-C1 code [9, 10], ELMing discharges in NSTX were seen to be limited by resistive current-driven peeling modes, with considerable sensitivity to plasma resistivity, whereas non-ELMing wide-pedestal discharges were located near the ideal pressure-driven ballooning threshold (see Fig. 1) [11]. Observations similar to the wide-pedestal H-mode have been made in enhanced-pedestal H-modes, i.e. these discharges are also limited by ideal pressure-driven modes. However, the linear stability picture is not as clear as for the other NSTX discharges, as nonlinear stability is found to play a role in stabilizing in the non-ELMing nature of these plasmas. It is now being investigated whether the impact of resistivity on P-B stability is a result of spherical tokamak geometry (aspect ratio, shaping, etc.) or profile alterations due to lithium coating in NSTX. The model thus enables higher fidelity predictions for ELM thresholds and presents a valuable basis in the quest for a predictive model for ELMs in low-aspect ratio tokamaks. Gyrokinetic analysis predicts that a variety of NSTX H-modes, from those with narrow pedestals and ELMs to wide ELM-free cases, are within 10% of KBM stability thresholds across the entire pedestal. This indicates KBM remains a viable candidate for constraining the maximum pressure gradient at low aspect ratio. Using gyrokinetic simulations to predict the onset of KBMs and other gyrokinetic instabilities, a new gyrokinetic linear threshold model [12] reproduces the NSTX experimental [13] width-height scaling (Fig. 2), which deviates significantly from standard aspect-ratio devices [14]. To predict NSTX(/-U) and future ST pedestal width and height, the complimentary gyrokinetic KBM and resistive P-B stability constraints will be combined. 2.1.2. Progress towards pedestal transport validation and predictive model The experimentally inferred ratio of electron particle to heat diffusivity (from the SOLPS code [15]) is smaller than predicted by just KBM. This is consistent with the existence of both microtearing modes (MTM) and electron temperature gradient (ETG) modes, which are expected to transport primarily electron heat flux, and are predicted to be unstable. Nonlinear electron-scale ETG simulations predict a range of electron heat flux (0-2 MW) depending on local parameters (such as Te gradient, for surfaces with r/a between ~0.90-0.98). Combined with neoclassical (NC) ion thermal transport (~ 1 MW, from the NEO code [16]), transport from ETG + NC accounts for 50-75% of the power flow across all surfaces for the wide pedestal and progressively less for the narrow pedestal. Nonlinear ion-scale MTM simulations predict electron heat flux comparable to that from ETG. Together, ETG + MTM + NC accounts for the total power flow in a wide-pedestal discharge in which the largest deviations from KBM transport ratios are observed. To support evolution towards a predictive capability, additional simulations were used to develop a reduced ETG pedestal transport model that reproduces many of the dependencies with driving gradients and equilibrium parameters, unifying the NSTX results with those from recently published analysis from the higher aspect ratio DIII-D tokamak [17]. The results of this reduced model are shown in Fig. 3. Linear gyrokinetics with the CGYRO code [18] predicts that MTMs become dominant over ion temperature gradient (ITG) modes as the density is increased near the pedestal region and are predicted to be the dominant instability in the core in H-mode plasmas at similar densities (normalized to the Greenwald density, nGW = Ip/πa², where Ip is the plasma current and a the minor radius) between NSTX(/-U) and DIII-D [19]. This underscores the importance of considering the density in applications of MTMs to STs. Global electromagnetic simulations of NSTX discharge 132588 using the GENE code [20], including carbon impurities which can't be neglected, confirm local simulation results that experimental profiles are in the KBM stable region, and that the most unstable mode is the MTM [21]. For KBM to become dominant, a 30% increase in beta would be required. MTM instabilities were found to grow near both core and edge rational safety factor locations. 2.2. Core transport Global confinement in NSTX H-modes has been observed to have scalings distinct from those in conventional aspect ratio devices. Thermal transport of core ions is typically found to be neoclassical in NSTX, while electron thermal losses dominate. Various gyrokinetic studies have been performed to explain this phenomenology. 2.2.1. Gyrokinetic validation / stability studies In addition to KBM modes that can set the profile shapes, for high-beta NSTX discharges and NSTX-U projections that span over an order of magnitude variation in collisionality, linear gyrokinetic simulations of ion-gyroradius-scale micro-instabilities also show a complex mix of MTMs, trapped electron modes (TEMs), and hybrid TEM/KBM modes. ITG instabilities, however, are typically stable in the NSTX discharges. In most cases, simulations showed that modes that can produce ion transport were found to have growth rates smaller than the flow shear rate or thresholds much higher than the experimental NSTX gradients, consistent with the observed neoclassical ion thermal transport [22]. The analysis suggests ITG instabilities are unlikely to contribute significant anomalous thermal losses in high-beta, lower collisionality NSTX-U scenarios, and that KBMs are likely to play a role limiting the confinement by sitting near the threshold in most cases. This is summarized in Fig. 4, where both ITG (stars) and KBM (circles) thresholds from various discharges are presented. Another CGYRO scan of aspect ratio and elongation found that while ETG or ITG were dominant for DIII-D, MTMs dominate in the NSTX range of parameters. Additionally, the high elongation and triangularity shaping capabilities of STs have also been shown, through calculations with the XGC code [23, 24], to reduce the linear growth rates of both MTMs and KBMs [25]. Equilibrium E × B sheared flows (electric field crossed with magnetic field) have sometimes been found to significantly suppress electrostatic ion scale transport in both gyrokinetic simulations and in experiment. For an NSTX case it was found that for sufficiently large magnetic shear, ŝ, MTMs are linearly stabilized with increasing ballooning angle θ₀ = kx0/(kyŝ) where kx0 is the radial wavenumber at the outboard midplane. The E × B sheared flow acts to move modes through θ₀, resulting in mode suppression. MTMs typically drive electron heat transport and nonlinear simulations using CGYRO found the MTM heat fluxes were suppressed when including E × B sheared flows matching experimental heat diffusivities. Meanwhile, nonlinear global gyrokinetic simulations of MTMs with the GTS code [26] are being carried out for experimental NSTX plasmas. As in previous local simulations, electron heat transport was found to be dominant, but additionally it was found that the safety factor q = 1 surface can play the role of transport barrier as turbulence propagates radially both inward, and equilibrium E × B sheared flow, while stabilizing MTM, was found to possibly destabilize an Alfvénic Kelvin-Helmholtz instability. 2.2.2. Reduced model development Gyrokinetic theory, while an excellent approach to properly describe the microturbulent transport, is computationally expensive, so reduced models are critical to achieve real time profile predictions. However, these models must be validated against gyrokinetic simulations; several efforts are proceeding. Turbulent and neoclassical heat transport have been calculated for NSTX with the flux-matching TGYRO [27] (TGLF [28] + NEO [16]) solver, predicting experimental ion and electron temperature profiles, Ti and Te [29]. For a low confinement (L-mode) plasma, linear stability analysis and scans of temperature gradients identified the low-k and high-k unstable modes driving the turbulent electron and ion heat transport at the outer core region ρ > 0.4, where k is the wavenumber and ρ is the square root of normalized toroidal flux. In an H-mode plasma ETG modes drive the turbulent electron heat transport while the low-k modes are suppressed and the ion heat transport is predominantly neoclassical (see Fig. 5). Comparison with linear gyrokinetic stability analysis shows close agreement of the real frequencies of unstable modes between TGLF and CGYRO gyrokinetic simulations, but higher growth rates are predicted by TGLF, especially for the H-mode case. An improvement in understanding of electron thermal transport is made possible by using a Multi-Mode-Model (MMM) [30] of anomalous transport, which includes a new physics-based electron temperature gradient (ETG) model [31] as well as an upgraded MTM model. The ETG model has been verified through its comparison with CGYRO simulations. The magnitude and associated real frequency of the most unstable ETG mode present in the kyρs spectrum agree with the results obtained from CGYRO [32]. The MMM model reproduces the experimentally measured thermal power of 2.0 MW [33] and produces electron temperature profiles that are consistent with the NSTX experimental data in a low collisionality discharge where ETG diffusivity is larger than MTM, as shown in Fig. 6 [32]. CGYRO was used to study linear and nonlinear ETG modes in NSTX(/-U) and compared with reduced ETG models to better determine their applicability [33]. Linear simulations determine the ETG critical gradient, contrasting it with the standard tokamak scaling formula [34]. This formula appears insufficient to describe the threshold when applied to STs. Within the deep core (r/a < 0.3), where the normalized temperature gradient is typically less than the experimental value, no ETGs manifest. In the transport region (r/a ~ 0.6 − 0.8), ETG stability varies across cases. Nonlinear simulations that calculate electron thermal transport align well with NSTX experimental data within uncertainties. The reduced ETG model [31], displays strong agreement with experimental measurement and CGYRO simulations. This, and the successes of TGYRO and MMM, instills confidence in their utility for updated profile prediction of future NSTX-U discharges [35]. 2.3. Impurity transport Besides transport of main ion species and electrons, transport of impurities and energetic ions is also important to understand. Resonant magnetic perturbations (RMPs) are sometimes applied in tokamaks to suppress ELMs. However, RMPs often result in a decrease in the plasma density, also termed density pump-out. The role of neoclassical transport in density-pumpout can be analyzed by using a coupling of the M3D-C1 and 3D-NEO codes [36]. The neoclassical fluxes of impurities evaluated for the NSTX discharge 130670 shows that, if an impurity resides in the low collisionality regime it will experience significant changes in 3D perturbed flux with RMPs (compared to the axisymmetric flux), but is weakly affected by RMPs in the Pfirsch-Schlüter regime [37]. Since the atomic number (Z) of impurity species affects the collisionality, the transport of low-Z impurities is greatly altered with RMPs but has a smaller effect for higher Z impurities. 2.4. RF Heating Simulation capabilities for RF heating are expanding greatly in general. The use of the Petra-M finite element method platform [38] allowed for simulations with very high-fidelity for the first time of the full 3D NSTX-U torus including realistic antenna geometry in the high harmonic fast wave (HHFW) heating regime [39]. More recently, RF sheath effects have also been implemented in Petra-M [40, 41]. RF sheaths can form during ion cyclotron RF operation and contaminate the plasma with impurities [40, 42] so it is very important to have a quantitative analysis for the future NSTX-U HHFW experimental campaign. A recently developed coupling between the Vorpal [43] and UEDGE [44] codes was employed to simulate NSTX-relevant scenarios with experimentally realistic plasma profiles and HHFW antenna parameters [45]. It was found that the ponderomotive contribution to the plasma parallel momentum can be significant for representative RF input power fluxes. Resonant wave-filament interactions have been shown to be a loss mechanism for HHFW power [46] and can explain the large fraction of RF power flows along the field lines in NSTX that hit the divertor, rather than heating the core plasma [47, 48]. A numerical study of the impact of the edge density fluctuations and filaments on the HHFW wave propagation have been recently started as well [49]. A density perturbation was applied to 3D HHFW simulations with Petra-M, and the fluctuations were seen to have a strong impact on the power absorption, spreading it vertically and scattering it into multiple regions. A large percent of HHFW injected power can be lost in the plasma boundary on NSTX [50]. New 1D particle-in-cell simulations are conducted to investigate a possible role of parametric decay instabilities (PDI) [51]. The present study identifies the PDI excitation near the lower hybrid resonance layer where the thermal effects becomes non-negligible, and can form a power loss channel. The modeling also suggests that parasitic electrostatic wave excitation, such as Ion Bernstein waves and hot ion plasma waves, can lead to anisotropic perpendicular ion heating, in line with the previous ion temperature measurement on NSTX during HHFW experiments [52]. When this resonance layer is removed in the simulation domain by raising the density, PDIs are suppressed despite a large increase in the input power, indicating that a presence of the lower hybrid resonance density layer can be a key parameter in determining the convective growth. This can be tested experimentally by moving the plasma toward the antenna or by gas puffing in the vicinity of the antenna. Machine learning based surrogate models are being employed to accelerate RF modeling codes [53]. In particular, fast and accurate predictions of HHFW heating in NSTX/NSTX-U have been obtained via training surrogates on full-wave solutions. High-fidelity predictions (i.e. R² ∈ [0.94 − 0.97]) are obtained for both electron and deuterium power deposition radial profiles. Figure 7 shows examples of (a) high, (b) average, and (c) low accuracy predictions of deuterium power deposition profiles together with the corresponding key simulation input parameters. The predictions show excellent agreement with the ground truth, denoting the capability of the surrogates to capture the main physics of HHFW absorption. The low accuracy case shown, while still capable of accurately capturing the physics, represents one of the worst performing cases in the surrogate models. The reason is due to an increase in the complexity of the power deposition profiles (multi-pass absorption) in this region of parameter space. The achieved average inference times on the order of microseconds are necessary to incorporate a HHFW model into integrated frameworks for both time-dependent simulations and real-time control applications. NSTX-U HHFW heating and current drive modeling with the GENRAY [54] and CQL3D [55] heating and current drive codes will also be used to inform the development of high performance scenarios using integrated modeling frameworks described later in section 3.3. Finally, transport of energetic particles (EPs) can be complicated when multiple sources of these particles are simultaneously employed, for example by neutral beam injection (NBI) and RF heating. Modifications of the NBI distribution on NSTX(/-U) as RF was added to the simulations was assessed with the ORBIT code [56], where increased losses of NB ions that are displaced to loss orbits by the RF field were found. The modified NBI distribution can subsequently have an effect on Alfvénic instabilities, either stabilizing or destabilizing, depending on how the RF tail distribution aligns and overlaps with Alfvén eigenmode resonances in the EP phase space. HHFW+NBI analysis was carried out with the Monte-Carlo MCGO code [57] as well [58]. Preliminary calculations show that RF produces larger width banana orbits, thus shifting lost ion deposition somewhat from the lower divertor region to the outer midplane. Though the total number of lost ions doesn't increase much, because of the high energy of the RF-accelerated lost ions, the heat load does increase. In summary of the first objective of NSTX-U research, to extend confinement and stability physics basis at low aspect ratio and high beta to lower collisionality, progress has been made on developing predictive pedestal structure and transport models, validating reduced core transport models, investigating impurity transport, and RF modeling to understand HHFW absorption and the impact on energetic particle stability. In order to achieve high performance discharges in NSTX-U and project low disruptivity operation to future steady state STs, various performance-limiting modes of instability were studied and predictive tools were developed. 3.1. Stability 3.1.1. Low-frequency modes, tearing modes, and islands One performance-limiting mechanism in NSTX was that the central electron temperature tended to remain largely unchanged as the external heating power increased above a certain level. Various explanations have been proposed for this behavior, including Alfven eigenmodes either directly influencing electron thermal transport [59], or modifying the power deposition from neutral beams [60]. Resistive 3D MHD simulations with M3D-C1 have now shown another possible explanation, that low toroidal mode number n, pressure-driven ideal MHD instabilities, though non-disruptive, acted to break magnetic surfaces in the plasma core, thereby flattening the electron temperature profiles above a critical beta value [61, 62]. Accounting for this effect for future devices may project somewhat lower performance, but it may be a more accurate prediction, and also motivates studies to potentially avoid these modes. Low frequency (< 50 kHz) magnetohydrodynamic (MHD) n = 1 and 2 mode activity was often present in the first 400 ms of NSTX discharges and may have inhibited even higher performance. Observations of these modes in relation to the equilibrium evolution, and a 3D study with M3D-C1 showed that they formed in the core when the safety factor at the magnetic axis, q0, reached a mode rational value, and then slowly moved outwards with the mode rational surface while remaining locked to the core, flattening the core Te and rotation profiles. When the rational surface moved sufficiently far from the core, the core would unlock from the mode. When the q0 evolution is fast enough to cross multiple rational surfaces before the core unlocks from the early mode, a characteristic “bifurcation” in the magnetic spectrogram can be observed (see Fig. 8). This understanding of these modes is the first step towards either avoiding them to increase performance, or at least maintaining a current evolution that allows the plasma to push through these early modes [63]. Electron thermal transport in NSTX was routinely characterized as anomalous [64]. Operating with reversed magnetic shear (RMS) was shown to suppress the anomalous transport through the development of electron internal transport barriers when the shear was sufficiently negative [65, 66]. By utilizing TRANSP, RMS operating scenarios were developed in NSTX-U [67]. It was found that procedures similar to those used in NSTX can be followed in NSTX-U to establish RMS. This includes creating as large a plasma as possible with a fast current ramp and early neutral beam injection to heat electrons and slow the current diffusion for the edge to the core. Additionally, all such scenarios have shear levels beyond the threshold found in NSTX for transport suppression. Early reversed shear profiles above q0 = 1 are potentially prone to resistive MHD instabilities known as double tearing modes, however. An example of such a mode from NSTX [65] has now been investigated with M3D-C1, where it was found that sometimes these modes can lead to flux-surface breakup that may persist throughout the discharge. Steady state operation will require fully non-inductive current, but simulations have shown that for large aspect ratio devices, a significant reduction of bootstrap current over large radial extent is attributable to magnetic islands. A gyrokinetic study of the effect of MHD islands on self-driven current in NSTX has now shown, however, that though tearing modes (and islands) are more prominent in high poloidal β regimes, the bootstrap current reduction is less sensitive to island width in this regime, and is therefore much less significant in STs [68]. One such model uses TRANSP code [69] to calculate the parameters of the modified Rutherford equation to predict the growth of magnetic islands [70]. An application of the model to the triggerless neoclassical tearing modes (NTMs) observed in NSTX showed that the model prediction agrees with the measurement when the EP term is added to the modified Rutherford equation [71]. The TRANSP “kick” model [72, 73] is used to calculate the input EP pressure and density profile self consistently, with the EP transport due to magnetic islands considered. Tearing modes may also be triggered by a sawtooth instability, which induces loss of EPs. An improved model of EP transport by sawteeth has now been implemented for TRANSP. Contrary to previous models already implemented in TRANSP, the new model resolves EP transport in phase space based on the workflow implemented for the kick model. The model has been tested on NSTX-U and JET discharges [74, 75], and it will enable a better quantitative understanding and prediction of plasma discharges. Finally, NTMs could be triggered by resonant magnetic perturbations (RMPs) in NSTX, and a theoretical simulation of that process showed that in the presence of a frequency mismatch between the rotation frequency of seed island and that of the RMP, the critical RMP amplitude needed to trigger an NTM oscillates as the duration of the RMP pulse is varied [76]. The critical amplitude is minimized when the RMP pulse duration is such that seed island chain executes a half-integer number of rotations with respect to the pulse. 3.1.2. Energetic particle stability Non-resonant n = 1 fishbone modes driven by EPs in plasmas with qmin slightly larger than one can lead to degradation of plasma confinement. These modes were studied with M3D-C1 by extending the code to include kinetic effects of both thermal ions and energetic ions, and the perturbation on magnetic flux surfaces as well as the transport of energetic particles can be calculated [77]. It was found that the thermal ion kinetic effects can cause an increase of the frequencies of the modes, and that Landau damping can provide an additional stabilization effect [78]. A method of analytic predictive fast ion transport modeling has been extended to NSTX-U regimes. Unlike in conventional tokamaks, the fast ion effective dynamical friction (drag) rate in STs is typically of the same order of the fast ion effective scattering rate [79]. A new theoretical formulation [80] for marginally unstable modes has reported that in those scenarios, resonances of Alfvénic eigenmodes tend to shift and split their most active regions of energy exchange with the fast ion sub-population, thereby effectively extending the instability drive range. It has been shown within a reduced analytic framework that the energetic particle distribution evolves according to a quasilinear diffusion equation when the instabilities are near marginal stability. Moreover, this quasilinear system preserves the exact mode saturation levels that are present in the fully nonlinear kinetic theory [81]. This analytic extension allows for a self-consistent, fully predictive and numerically efficient fast ion transport module within whole device modeling which is numerically realized through the Resonance-Broadened Quasilinear (RBQ) code [82]. Toroidal Alfvén eigenmode (TAE) avalanches in NSTX were always accompanied by simultaneous rapid frequency-chirping and large amplitude bursting of multiple modes, and the avalanches led to significant neutron rate drops resulting from fast-ion energy loss and loss of beam ions [83]. Using M3D-C1, the down-chirping frequency decrease has been effectively replicated in nonlinear simulations. The chirping rate was found to be closely related to the continuum damping of the TAEs and the shear of the toroidal rotation. In the nonlinear stage, a substantial portion of EPs can be transported from the core to the outer region. This transport leads to increased interaction with the perturbed fields of the mode, resulting in a further amplification of the mode. The larger experimental phenomenon [84] has now been reproduced by non-linear, multiple wave number simulation results, which identified wave–wave nonlinear coupling among different modes as an important ingredient for the onset of TAE avalanches, during which there is a resonant interaction between different modes and energetic particles [85]. Though the frequencies are somewhat larger, and the time scale is shorter, the simulations in Fig. 9 qualitatively reproduce the experimental TAE avalanche behavior. Finally, Global Alfvén eigenmodes (GAEs) in NSTX driven by energetic particles have been linked to enhanced electron transport and may be expected to affect the performance of NSTX-U and future STs where super-Alfvénic fast ions might be present. Studies of the amplitude and phase of bursts of GAEs in NSTX indicate that the first harmonic can non-linearly drive weakly damped modes in the second harmonic frequency band [86]. It is further seen that these short wavelength Alfvénic modes can be strongly toroidally localized. Fully nonlinear simulations of the evolution of unstable GAEs in NSTX-U demonstrate the importance of including all toroidal harmonics for accurate prediction of saturation amplitudes, and show a significant redistribution of resonant fast ions and modification of the beam distribution even by relatively small amplitude modes [87]. 3.2. Disruptions and machine protection 3.2.1. RWM stability predictions Resistive wall modes (RWMs) are global MHD modes that appeared at high β in NSTX, and their stability was understood through years of work on kinetic modifications to ideal MHD stability [88, 89]. Complex physics calculations of RWM stability [90, 91] have been distilled into reduced models for implementation in disruption avoidance schemes [92], and machine learning techniques guided by physics theory, have been tested as a supporting tool [93]. This work has been developed further by using the novel technique of counterfactuals to assess hypothetical magnetohydrodynamic activity or normalized β levels that would have prevented the RWM from going unstable (in this case, staying under the no-wall ideal MHD limit of βN ~ 4) in experimentally unstable discharges, motivating the usage of the counterfactual technique to simulate real-time control [94]. 3.2.2. Disruption event characterization and forecasting The reduced kinetic RWM stability model is just one of many physics modules implemented in the Disruption Event Characterization and Forecasting (DECAF™) code [95], developed to understand the chains of events leading to disruption in many machines, including NSTX(/-U) [96]; further examples of modules include the following. First, standardization of the determination of disruption timing between machines will allow comparable analysis of multi-machine databases, and DECAF has now been updated to utilize abnormal plasma current and vertical position signals to accomplish this [97]. An example DECAF analysis of the NSTX-U 2016 database is shown in Fig. 10. This plot shows the disruption triggering event in an operational space plot, giving a different detail than disruptivity diagrams [98]. Second, vertical stability detection approaches have themselves been compared between real-time reconstructions and magnetic probe measurements, resulting in critical metrics and thresholds for predicting vertical displacement events. Third, automatic detection of H to L-mode back transitions based on Te, deuterium-alpha radiation, Dα, stored energy, and energy confinement time signals has been implemented to study the correlation and causality of confinement back transitions with disruptions. Finally, density limits can lead to disruption, and the trivially calculated Greenwald limit has been further tested against the NSTX and Mega-Amp Spherical Tokamak (MAST), and MAST-U databases [99, 100]. The physics of the density limit has been recently explored theoretically through local phenomenon, such as power balance at magnetic islands [101] or turbulent transport at the plasma boundary [102, 103]. The former has been tested as a disruption indicator in NSTX-U discharges [95, 100], but so far was not found to be as reliable as the Greenwald limit. The latter was tested vs. some MAST-U discharges and found to be potentially useful for spherical tokamak disruption prediction [100]. 3.2.3. Runaway electrons Runaway electron (RE) generation presents a possibly large issue for future devices due to the potential for damage to the structures. In tokamaks RE can be generated after disruptions but also during the early stages of the discharge start-up, when higher electric fields are needed to increase the current through cold plasmas. NSTX contributed unique low aspect ratio data to a multi-machine database for RE generation and confinement during startup [104]. The study found that there is not an operation range that avoids the generation of RE during tokamak start-up, but that under certain circumstances (e.g. low-density operation), more REs are generated such that they are detected and/or can do damage. The multi-machine comparison demonstrated that the thermalization of the startup REs and the suppression of secondary generation of REs is sensitive to the dynamics of the density, the electric field and thus also the temperature during the start-up phase. It was found that larger devices have typically higher temperatures during the start-up hence slowing down the RE generation process. Data from NSTX at a low aspect ratio was critical for demonstrating the drift-orbit losses of REs are not the primary loss mechanism for startup REs [104]. 3.3. Equilibrium reconstruction, integrated modeling, and control Tools are being developed for NSTX-U for fast prediction and optimization of plasma scenarios, including neural networks trained on the EFIT equilibrium reconstruction code [105] as well as prediction of non-rigid plasma response for shape control [106]. These can be used in forward-mode for simulation of plasma scenarios before they are run, or reconstruction mode for real-time equilibrium reconstruction, and have shown some performance improvements over the existing real-time EFIT. The OMFIT integrated modeling workflow [107] has been set up for NSTX and NSTX-U, including obtaining the full kinetic equilibrium reconstruction, analysis of power balance fluxes, and prediction of heat plasma profiles based on a reduced turbulence model. The integrated workflow addresses tasks of kinetic profiles fitting integrated with the transport code TRANSP and EFIT to provide a self-consistent equilibrium reconstruction with kinetic constraints of total pressure, including beam pressure, and total current, including bootstrap current. The latest version of EFIT, EFIT-AI, can be used in this workflow to construct higher resolution equilibrium reconstructions [108]. GSevolve [109], a free boundary equilibrium code that evolves the Grad-Shafranov equilibrium including current and pressure profiles, is being prepared to be connected to the same control software as the experiment, the NSTX-U plasma control system [110]. The inclusion of high resolution current profile simulation provides a close match to interpretative experimental analysis from TRANSP. NSTX-U discharges with challenging dynamics, including vertical oscillations [111] and full vertical displacement events have been simulated successfully. The ability to simulate discharges from breakdown to termination is included. GSevolve simulations can improve NSTX-U control and performance by, for example, controlling and increasing elongation in the plasma current rampup. Tight simultaneous regulation of several plasma parameters, both scalars and profiles, during the tokamak discharge is important for maintaining steady-state, non-disruptive scenarios. Controllers with feedforward and feedback capability, based on both off-line and on-line optimization are being designed for plasma-scenario regulation in NSTX-U. A linear quadratic integral regulator has been synthesized for simultaneous regulation of the q profile and normalized beta βN [112], while a Model Predictive Controller (MPC) has been developed for simultaneous control of multiple scalar quantities, including q0, βN, and internal inductance [113]. Moreover, a hybrid MPC for dual q-profile and stored energy regulation has been proposed to explicitly incorporate the pulse-width-modulation constraints imposed by the NSTX-U NBI system [114]. Finally, a method of significantly reducing the computational time of a model-based feedforward-control optimization scheme for q-profile shaping in NSTX-U has also been devised, whereby analytical cost function gradients are used in place of numerical gradients in order to ensure a faster convergence to the optimized actuator trajectories [115]. In summary of the second objective of NSTX-U research, to develop non-inductive operation at high-performance and low-disruptivity, progress has been made in understanding performance-limiting MHD, NTMs and disruptions, developing predictions for mechanisms hypothesized to flatten central electron temperature, developing predictive models for burning plasma relevant energetic particle dynamics, and developing tools for integrated modeling tools, scenario design and experimental planning. The increased heating power and compact geometry of NSTX-U will produce very high heat flux levels that will be seen in future high performance STs as well. Nominal heat fluxes are up to 8 MW/m², however in specific conditions NSTX-U could generate more than an order of magnitude higher [116]. Innovative solutions to handle these heat fluxes and predictive capabilities accounting for plasma phenomena must be developed, tested, and improved. 4.1. Lithium The NSTX(/-U) program of research has pioneered the study of lithium as a renewable surface at the plasma-material interface that protects the underlying solid substrate, improves confinement through particle pumping, and enhances power exhaust capability, including heat flux redistribution due to vapor shielding from evaporated lithium. In addition to the glow discharge boronizations employed by the majority of worldwide fusion research devices, NSTX-U also plans to utilize a set of LIThium EvaporatoR (LITER) [117] units to provide depositional lithium coatings on the plasma facing surfaces. These coatings have been shown to increase confinement, reduce recycling, and suppress ELM activity [118]. However, the cessation of ELMs generates a concomitant core impurity buildup leading to radiative discharge collapse. To mitigate this effect an impurity granule injector (IGI) [119] will also be fielded. The IGI provides horizontal injection of sub-mm impurity granules which provide a localized density perturbation stimulating an ELM and subsequently mitigating the core impurity burden. By tuning the pacing frequency of the granule injection, a balance can thus be struck between the increased performance provided by the evaporative lithium coatings and the impurity buildup resulting from the elimination of natural ELMs. 4.1.1. Lithium vapor box The fact that NSTX-U can produce very high heat fluxes makes it an excellent candidate for a future test facility for lithium divertor solutions. One such solution is to dissipate the heat flux in a simple configuration called a lithium vapor box (LVB) that evaporates lithium in the private flux region [120]. Divertor detachment is desired to protect plasma facing components (PFCs), but can create a highly radiating X-point and upstream ionization could lead to non-negligible lithium content, so a LVB [121] baffled design has been studied with the SOLPS-ITER code using an NSTX-U magnetic equilibrium. Private flux region gas puffing was shown to reduce the upstream lithium contamination significantly more than common flux region puffing due to better access to the separatrix [116] (Fig. 11). Lithium flows primarily along the separatrix due to the hotter field line having a stronger thermal force which pushes the lithium upstream. The effect of the fuel puff for a given intensity depends strongly on the recycling coefficients assumed for the various plasma facing components in the simulation. Test stand experiments of lithium vapor boxes are being implemented [122], and such a concept could be considered for future use on NSTX-U. Additionally, many lithium vapor-shielding high-density divertor concepts rely on non-coronal lithium radiation that can be affected by radiation transport. The non-coronal lithium radiated power loss enhancement is due to fast impurity transport and atomic processes like charge exchange with neutrals, and comes about when compared with the radiated power loss of a plasma in coronal ionization-recombination equilibrium without plasma transport. However, when plasma opacity to its own radiation is considered, additional corrections to ionization and recombination rates of all lithium charges states, and additional corrections to line and continuum radiation may result in significantly modified lithium radiated power loss. Work is now underway to clarify the role of lithium radiation transport on radiated power and charge state distribution in these divertors [123]. 4.1.2. Flowing liquid Li plasma facing components Another approach is to place liquid lithium directly at the strike point. In one such design [124], a porous wall is used to stabilize the liquid metal surface, while an MHD drive is used to push the liquid metal flow underneath the porous surface. This same concept has now also been optimized to work as an evaporating surface for the LVB [121]. Analytical and numerical models for liquid lithium PFCs [125] allow parametric studies of the design variants and optimization, as well as detailed mapping of the flow and heat transfer distributions. The plasma heat flux distribution is obtained using SOLPS analysis, while numerical simulation of the plasma facing components uses a version of the 3D computational fluid dynamics code Ansys CFX. It can predict the temperature distribution on the liquid metal plasma facing component, in different design conditions. The resulting distribution can then be used to establish a two-way coupling between an analytical model and SOLPS. The code was modified to simulate MHD flows at high Hartmann numbers. Finally, the optimized design of the flowing liquid lithium device with porous walls has now been used to create conceptual replacement tiles for the NSTX-U divertor. Finally, many challenges remain for lithium solutions to be viable for fusion reactors. For example, a complete system that re-collects the evaporated lithium [126] will need to be tested on a long-pulse device with hot walls, and both hydrogenics and impurities will need to be extracted [127]. The liquid lithium program in NSTX-U is geared to close some of the gaps to a reactor, and the designs are coupled to reactor designs of liquid lithium PFCs [128]. 4.2. Divertor heat flux 4.2.1. 3D calculations of heat deposition Heat fluxes on NSTX-U solid PFCs (divertor tiles) have previously been analyzed in 3D and with time variation with the HEAT code [129]. The code has now been updated with an ion-gyro orbit module, and the gyro-orbit effects were shown to alternately enhance PFC performance by smearing magnetic shadows and degrade performance when narrow regions on edges and corners are loaded with high heat fluxes [130]. 4.2.2. Reduced model for scrape-off-layer A reduced model for scrape-off layer (SOL) plasma transport was developed and tested against NSTX(/-U) experimental observations and high fidelity interpretive simulations [131]. The model addresses core-edge physics problems such as edge ionization sources and neutral densities by coupling to the kinetic neutral transport code DEGAS-2 [132], and fast ion confinement with a novel numerical algorithm for integrating the stochastic differential equation for pitch angle scattering. The model can additionally use Langmuir probe data to reconstruct the calculation of the upstream density and electron temperature and can be used in TRANSP to constrain recycling coefficients, thereby providing an efficient and fast coupling between core and edge to support experimental planning and discharge scenario development. 4.2.3. Edge harmonic oscillations, turbulence, and edge localized modes Other plasma phenomena can influence heat flux as well, and it is important to understand these to accurately predict the performance of PFCs for NSTX-U, where the predictions can be tested and then confidently projected to future, higher power STs. One of these phenomena is the edge harmonic oscillation (EHO), which was analyzed and found to be beneficial by decreasing the divertor peak heat flux when background edge fluctuations were low, and increasing the heat flux width (more with larger EHO frequency) [133]. However, when background edge fluctuation levels were high, the EHO was found to increase the divertor peak heat flux. The divertor peak heat flux decreases with the frequency of EHO n = 1 mode, while the heat flux width increases with the frequency of the EHO n = 1 mode. High plasma edge turbulence in NSTX H-modes was found to significantly increase the divertor heat flux width to levels comparable with the divertor heat flux in L mode, and possibly up to three times larger than the multi-machine scaling results [134]. ELMs may cause problems for future STs by releasing bursts of energy from the plasma, but they were found in NSTX to cause heat flux striations that were beneficial (to the divertor, at least) by significantly increasing the heat flux width and decreasing the divertor peak heat flux during the ELMing period vs. the inter-ELM period, also to levels similar to L-mode. A characteristic signature for the effect of small ELMs on the lower outer divertor heat flux is in NSTX is shown in Fig. 12. After the L-H transition at ~ 0.153 s (red dashed line), the divertor heat flux width shrunk, which caused the peak heat flux, qpeak, to increase from 1.2 MW/m² to ~ 3.5 MW/m² (Fig. 12(d)). After 0.22 s the ELMs appear, and the deposited power on the lower divertor increased by < 50% with small ELMs (these are similar to type V ELMs published in Ref. [135]). However, qpeak significantly decreased from 3 MW/m² (ELM free), to < 1.5 MW/m² (small ELMs), a similar level to L-mode discharges. 4.3. Blobs ELM filaments transport energy and particles out of the confined plasma region to the SOL and eventually to the PFCs. These structures are similar to the intermittent “blobs” [136] appearing in the background SOL turbulence. Gas-puff imaging data from NSTX has been analyzed and while ELM filaments and blobs were found to be of similar size and shape, the ELM filaments internally rotated about three times faster and the angular velocity increased with the distance from the separatrix [137, 138]. Novel data analysis techniques were needed to be developed to assess the angular velocity on a frame-by-frame time resolution [139]. The understanding of filamentary dynamics is important to develop the ability to predict and mitigate the impact of ELMs on plasma facing components. The blobs were also studied and in NSTX were shown to decrease in quantity with NBI power, be less prevalent in H-mode than L-mode, and depend most strongly on the poloidal turbulence velocity [140]. Temporally, blobs came in pulses that lasted on the order of 25 µs, while the wait time between pulses was typically about 1 ms [141]. A theory for the mechanism of blob formation based on velocity shear breaking radially elongated streamers was satisfactorily tested against NSTX data [142, 143]. A model of the dynamics of blobs which accounts for interaction of electromagnetic perturbations with the conducting wall has also been tested for NSTX cases [144]. The effect of applied n = 1 and 3 non-axisymmetric magnetic perturbations in NSTX on edge turbulence was studied as well, but it was found that they did not affect the average poloidal size, autocorrelation time, or relative fluctuation levels of the turbulence. Finally, in summary of the third objective of NSTX-U research, to develop and evaluate conventional and innovative power and particle handling techniques, progress has been made in developing liquid/vapor lithium PFC concepts and plans for future NSTX-U missions, developing models for SOL transport and PFC heat loads, and characterizing those heat loads and edge transients. 5. Outlook The NSTX-U recovery project has been re-baselined [145]. A major task is the replacement of the central magnetic assembly, and at the present time vendors are under contract and the bundle fabrication is ongoing. All other major components of the recovery project have been delivered and machine reassembly tasks have started. The first plasma in the new phase of NSTX-U operation is now expected in late 2025. In the meantime research on NSTX(/-U) continues, using data from previous run campaigns or developing tools and models for analysis of future experiments. When it returns to operation NSTX-U will be well poised to fill the gaps in physics understanding between current spherical tokamaks and future ST fusion pilot plants. Table 1. Acronyms used in the paper - CQS: Current quench start - DCS: Disruptive current spike - DIII-D: Tokamak at General Atomics - EHO: Edge harmonic oscillation - ELM: Edge localized mode - EP: Energetic particle - ETG: Electron temperature gradient mode - GAE: Global Alfvén eigenmode - HHFW: High harmonic fast wave - IPR: Plasma current not meeting the request - ITG: Ion temperature gradient mode - KBM: Kinetic ballooning mode - LVB: Lithium vapor box - MAST: Mega Ampere Spherical Tokamak - MHD: Magnetohydrodynamics - MPC: Model predictive control - MTM: Micro-tearing mode - NBI: Neutral beam injection - NC: Neoclassical - NSTX-U: National Spherical Torus Experiment - Upgrade - NTM: Neoclassical tearing mode - P-B: Peeling-ballooning mode - PDI: Parametric decay instabilities - PFC: Plasma facing components - RE: Runaway electrons - RF: Range-of-frequency waves - RMP: Resonant magnetic perturbation - RMS: Reversed magnetic shear - RWM: Resistive wall mode - SOL: Scrape-off layer - ST: Spherical tokamak - TAE: Toroidal Alfvén eigenmode - TEM: Trapped electron mode - VDE: Vertical displacement event Table 2. Codes used in the paper - Ansys CFX: A computational fluid dynamics code - CGYRO: A collisional, electromagnetic, multiscale gyrokinetic code [18] - CQL3D: Collisional/quasiLinear 3D code [55] - DECAF: Disruption event characterization and forecasting code [95] - DEGAS-2: A kinetic neutral transport code [132] - EFIT: An equilibrium reconstruction code [105] - GENE: Gyrokinetic electromagnetic numerical experiment code [20] - GENRAY: General ray-tracing code [54] - GSevolve: A free boundary equilibrium code [109] - GTS: Gyrokinetic tokamak simulation code [26] - HEAT: Heat flux engineering analysis toolkit [129] - M3D-C1: A non-linear, 3D, resistive MHD code [9, 10] - MCGO: Monte Carlo guiding center orbit code [57] - MMM: Multi-Mode-Model of anomalous transport [30] - NEO: A multi-species drift-kinetic solver [16] - OMFIT: One modeling framework for integrated tasks [107] - ORBIT: A guiding center code [56] - PETRA-M: A 3D full-wave code [38] - RBQ: Resonance-broadened quasilinear code [82] - SOLPS: A multi-fluid (ions, electrons, neutrals) edge modelling code [15] - TGLF: A gyro-Landau-fluid model code [28] - TGYRO: A parallel transport manager that combines TGLF and NEO [27] - TRANSP: A plasma transport code [69] - UEDGE: A time-dependent 2D plasma fluid equation code [44] - Vorpal: An electromagnetic / plasma wave solver [43] - XGC: X-point included gyrokinetic code [23, 24] Acknowledgments This work was supported by the U.S. Department of Energy under contract numbers DE-AC02-09CH11466 (PPPL), DE-SC0008309 (UT, Knoxville), DE-SC0013977 (Lehigh), DE-SC0021385 (Lehigh), DE-SC0021113 (GA), DE-FG02-95ER54309 (GA), DE-FG02-91ER54109 (MIT), DE-SC0021156 (UT, Austin), DE-SC0021625 (Nova Photonics), DE-SC0021311 (Columbia), and DE-AC52-07NA27344 (LLNL). Additionally: EP/R034737/1 (UKAEA). A part of this research used resources of the National Energy Research Scientific Computing Center (NERSC), which is supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231. A part of this work was supported by the Department of Energy early career research program. Part of the data analysis was performed using the OMFIT integrated modeling framework [107]. The United States Government retains a non-exclusive, paid-up, irrevocable, worldwide license to publish or reproduce the published form of this manuscript, or allow others to do so, for United States Government purposes. Data Availability Statement The digital data for this paper can be found in: https://dataspace.princeton.edu/handle/88435/dsp01j6731612k