Characterizing the negative triangularity reactor core operating space with integrated modeling
H. S. Wilson A. O. Nelson J. McClenaghan P. Rodriguez-Fernandez J. Parisi C. Paz-Soldan


[
  \begin@twocolumnfalse
    



^1 Department of Applied Physics and Applied Mathematics, Columbia University, New York, NY 10027, USA
^2 General Atomics, San Diego, CA 92121, USA
^3 Plasma Science and Fusion Center, Massachusetts Institute of Technology, Cambridge, MA 02139, USA
^4 Princeton Plasma Physics Laboratory, Princeton, NJ 08540, USA

Characterizing the negative triangularity reactor core operating space with integrated modeling
    
H. S. Wilson^1,
A. O. Nelson^1,
J. McClenaghan^2,
P. Rodriguez-Fernandez^3
J. Parisi^4
and
C. Paz-Soldan^1

    
================================================================================================================



§ ABSTRACT 
Negative triangularity (NT) has received renewed interest as a fusion reactor regime due to its beneficial power-handling properties, including low scrape-off layer power and a larger divertor wetted area that facilitates simple divertor integration. NT experiments have also demonstrated core performance on par with positive triangularity (PT) H-mode without edge-localized modes (ELMs), encouraging further study of an NT reactor core. In this work, we use integrated modeling to scope the operating space around two NT reactor strategies. The first is the high-field, compact fusion pilot plant concept MANTA (The MANTA Collaboration et al 2024 Plasma Phys. Control. Fusion 66 105006) and the second is a low field, high aspect ratio concept based on work by Medvedev et al (Medvedev et al 2015 Nucl. Fusion 55 063013). By integrating equilibrium, core transport, and edge ballooning instability models, we establish a range of operating points with less than 50 MW scrape-off layer power and fusion power comparable to positive triangularity (PT) H-mode reactor concepts. Heating and seeded impurities are leveraged to accomplish the same fusion performance and scrape-off layer exhaust power for various pressure edge boundary conditions. Scans over these pressure edge conditions accommodate any current uncertainty of the properties of the NT edge and show that the performance of an NT reactor will be extremely dependent on the edge pressure. The high-field case is found to enable lower scrape-off layer power because it is capable of reaching high fusion powers at a relatively compact size, which allows increased separatrix density without exceeding the Greenwald density limit. Adjustments in NT shaping exhibit small changes in fusion power, with an increase in fusion power density seen at weaker NT. Infinite-n ballooning instability models indicate that an NT reactor core can reach fusion powers comparable to leading PT H-mode reactor concepts while remaining ballooning-stable. Seeded krypton is leveraged to further lower scrape-off layer power since NT does not have a requirement to remain in H-mode while still maintaining high confinement. We contextualize the NT reactor operating space by comparing to popular PT H-mode reactor concepts, and find that NT exhibits competitive ELM-free performance with these concepts for a variety of edge conditions while maintaining relatively low scrape-off layer power.

  \end@twocolumnfalse
]








§ INTRODUCTION


While great progress has been made in the fusion energy field toward power-plant relevant plasma performance and divertor technology, a major outstanding challenge remains: the coupling of a high-performance core to a realistic exhaust solution. This challenge has shifted the focus of some plasma core modeling efforts away from maximizing power output and toward optimizing power handling potential. Fusion pilot plant (FPP) tokamak concepts have been dominated by positive triangularity (PT) plasmas operating in high confinement mode (H-mode). H-mode is accessed when a PT plasma is given sufficient heating and fueling and is characterized by the formation of an edge transport barrier with high pressure gradients called a pedestal <cit.>. Given its higher confinement, H-mode is generally seen as a desirable regime for an FPP. However, the power through the scrape-off-layer PSOL must be above the L-H mode transition power PLH estimated by scaling laws <cit.> to sustain H-mode. In a PT H-mode reactor-class device, this leads to heat fluxes that will be difficult for plasma facing components to sustain without an advanced divertor <cit.>. Even if we assume we can advance divertor and material technology to sustain reactor-level H-mode heat loads, H-mode bears yet another challenge: it is accompanied by edge localized modes (ELMs) <cit.>. ELMs are instabilities that can result in large energy fluences to plasma facing materials if not mitigated in some way <cit.>. As power, current, and magnetic field are increased to reactor-relevant levels, the machine damage from ELMs is likely to be intolerable <cit.>.

There are multiple regimes that exhibit higher confinement than the traditional low-confinement mode (L-mode) with smaller ELMs than H-mode or that avoid them all together. These include I-mode, QH-mode, EDA H-mode, and quasi-continuous exhaust (QCE), among others <cit.>. While PT no-ELM or small-ELM regimes are better than PT H-mode for device longevity, they often have sensitive access conditions depending on the heating and fueling scheme or do not yet provide sufficient fusion performance improvement over L-mode to support a realistic FPP design <cit.>. 

Recently, the negative triangularity (NT) regime has resurfaced as a potential solution to the core-edge integration challenge in a reactor-class tokamak <cit.>. The “upper" and “lower" triangularities of a toroidal plasma are defined as δu,l = (Rgeo-Ru,l)/a, where Rgeo is the geometric major radius, Ru is the major radius at the highest point of the separatrix, Rl is the major radius at the lowest point of the separatrix, and a is the plasma minor radius. The average triangularity δ is the mean of the upper and lower triangularities. Because NT plasmas have Ru,l>Rgeo, NT x-points are located at a larger major radius than PT x-points. This allows for more space for a divertor and a larger divertor-wetted area, both of which are beneficial from a power-handling perspective <cit.>. Further, the NT regime has been shown to be ELM-free as long as the triangularity is sufficiently negative, even in cases where PSOL exceeds PLH by a significant margin <cit.>.  

Importantly, experiments have shown that a plasma with NT shaping exhibits improved confinement over PT L-mode plasmas with otherwise similar parameters (current, magnetic field, density, and auxiliary heating) in both DIII-D <cit.> and TCV <cit.>. This improved confinement is attained without entering H-mode, so PSOL is not required to be greater than PLH as it would be in PT H-mode, which allows for the use of techniques like seeding noble gas impurities in NT to further lower PSOL while maintaining plasma performance <cit.>. Recent experiments on DIII-D have also extended the observed operating space in NT to reactor-relevant levels in non-dimensional parameters <cit.>. In a diverted configuration, DIII-D NT plasmas simultaneously exhibited βN > 3, fGr > 1, and q95 < 3 with H98y2 > 1 where βN is the normalized plasma beta, fGr is the Greenwald fraction <cit.>, q95 is the safety factor at ψN = 0.95, and H98y2 is normalized confinement time from the ITBP H98y2 energy confinement scaling law <cit.>. Additionally, recent gyrokinetic simulation work has found the reduction in turbulent transport in NT to be independent of machine size <cit.>, further encouraging the use of NT in a reactor-class device.

While a few NT FPP designs have been proposed at varying levels of fidelity <cit.>, these studies primarily focused on one design point and the performance trade-offs between various input parameters have not yet been fully established. To provide greater context for these trade-offs, we evaluate the performance of a reactor-class NT tokamak around two published operating points: the MANTA design <cit.> and the larger design from Medvedev et al <cit.>. We accomplish this by using the STEP code <cit.>, which facilitates predictive integrated modeling through easy and self-consistent data transfer between various equilibrium, stability, and transport codes. The two NT design points studied in this work differ most notably in size, magnetic field, and current, as described in section <ref>. Specifically, we investigate changes to Pfus and PSOL that result from changes in temperature and density profiles, auxiliary power Paux, seeded impurity fraction fimp, triangularity δ, and major radius Rmaj. The core effects from changes in toroidal magnetic field Bt, volume, and plasma current IP will be investigated through comparison of the smaller volume, higher magnetic field, lower current MANTA design and the larger volume, lower magnetic field, higher current Medvedev design. Of particular interest to us is the effect of the edge pressure boundary condition on fusion performance in NT, as a full characterization of the NT edge is currently absent from the literature. 

In section <ref>, we introduce the high-field (MANTA-like <cit.>) and high-volume (Medvedev-like <cit.>) base cases around which we analyze NT reactor performance. We describe the integrated modeling workflow used with the STEP code <cit.> for the majority of simulations in this work. The density profile changes from a reactor-relevant particle source, i.e., a particle source localized toward the plasma edge, are investigated. We determine that the Angioni scaling <cit.> predicts a density peaking similar to that which would be evolved to in TGYRO from a near-edge particle source. Thus we implement a psuedo-evolving scheme for density for simplicity in subsequent scans, as described in section <ref>. In section <ref>, we discuss Pfus density changes from δ and Rmaj scans in a high-field core, and find that geometry is not as leveraging as the electron pressure at ρ = 0.8 pe,0.8. In section <ref>, we discuss the NT edge and its present uncertainty. Due to the lack of a physics-based predictive model of an NT edge, we scan various temperature and density edge boundary conditions (Te,0.8 and ne,0.8, respectively) and find that Pfus and PSOL are both highly dependent on both Te,0.8 and ne,0.8. We use infinite-n ballooning stability codes on the region beyond TGYRO evolution (ρ = 0.8 to ρ = 1.0) to determine that a H98y2≈ 1 (confinement is what would be expected from a PT H-mode plasma with similar parameters) and fGr≈ 1 high-field NT operating point is likely feasible from a ballooning stability standpoint for a variety of potential pedestal widths. In section <ref>, we evaluate the impact of Paux on Te,0.8, and find that a relatively high temperature boundary condition is required for sufficient Pfus and cannot be overcome by additional Paux. We compare the performance of NT reactor-like core scenarios at various T,e,0.8 and Paux to other published FPP concepts. In section <ref>, we determine that at levels of intrinsic impurities (helium, tungsten) assumed by other FPP designs, there is minimal effect on fusion power when compared to the effect of seeded impurities utilized in this work (krypton). Including krypton at various concentrations results in a nearly linear downward trend in PSOL from increased impurity fraction and a potential for using impurity fraction to optimize Pfus.






§ SIMULATION SETUP AND METHODS





 §.§ Establishing a high-volume and a high-field NT reactor base case



There are currently two main strategies to reach reactor-relevant performance in tokamaks: the high-volume approach and the high-field approach.

For PT H-mode, two notable FPP designs utilizing the high-volume approach are ARIES-ACT2 <cit.> and EU-DEMO <cit.>. Meanwhile, the high-field approach is well-represented by ARC-class devices <cit.> and generally by the SPARC project <cit.>. For reference, basic parameters describing these concepts are displayed in table <ref>. In this work, we initialize a base case for both strategies applied to NT, using work by Medvedev et al <cit.> and the MANTA collaboration et al <cit.> as starting points for the high-volume and high-field strategies, respectively. This enables us to evaluate advantages and disadvantages of both strategies in the NT operating space.

MANTA (Modular, Adjustable, Negative Triangularity ARC) is a high-field (Bt≈ 11 T) NT FPP design. It is compliant with requirements laid out in the National Academy of Science, Engineering and Medicine’s (NASEM) report “Bringing Fusion to the U.S. Grid” <cit.>, made possible in part by utilizing a FLiBe liquid immersion blanket, demountable HTS magnet joints, seeded krypton, and an NT core that exhibits sufficiently high confinement without ELMs. The simple divertor design of MANTA requires that PSOL be below 40 MW for a separatrix density of 0.9 × 10^20/m^3 <cit.>. Fusion power is additionally constrained to be within 400-500 MW when using 40 MW of auxiliary heating to meet the NASEM net electric goal of ≥ 50 MWe <cit.>. 

Compared to MANTA, the design outlined in Medvedev et al. is a lower field (B_t≈ 6 T) larger major radius (R_maj≈ 7 m) NT tokamak <cit.>. The primary focus of <cit.> was the MHD stability of a theoretical NT reactor-class design. As such, only the pressure profile is reported in <cit.>; the density and temperature profiles are not explicitly constrained. To facilitate a direct comparison of a high-volume NT case with a high-field NT case, we use the PRO-create module in OMFIT <cit.> to initialize similar profiles in the high-volume case to those in MANTA. Due to the high volume, the separatrix density had to be lowered significantly to avoid exceeding the Greenwald limit <cit.>, and consequently the ne,0.8 scans performed in section <ref> are over lower values than in the high-field case.



All scans in this work are done around either the high-volume or a high-field operating points with H98y2≈ 1 and fGr≈ 1 outlined in table <ref>. Parameters taken directly from references <cit.> and <cit.> are in bold. Figure <ref> shows the equilibrium cross sections of both reference cases in R and Z coordinates. As will be discussed in section <ref>, there is significant uncertainty in predicting the edge condition for an NT FPP. To partially accommodate this uncertainty, we restrict the edge of our high-volume and high-field base cases via upper bounds on the normalized parameters H98y2 and fGr because they are both largely affected by the edge pressure condition in any plasma. Reactor-class fusion devices may be able to operate at fGr>1 because they will exhibit higher power densities than current devices <cit.>. Supporting this idea, non-ELMing NT plasmas have accessed H98y2 > 1 simultaneously with fGr > 1 on DIII-D <cit.>. For these reasons, we chose to establish the high-volume and high-field base cases with normalized confinement time and density H_98y2≈ 1 and fGr≈ 1, respectively. Obtaining H98y2≈ 1 and fGr≈ 1 is accomplished by varying the temperature and density at ρ = 0.8 until a converged solution is found, where the method to find a converged solution is described below in subsection <ref>. 





 §.§ Integrated modeling with the STEP code



The integrated modeling in this work was performed with the STEP (Stability, Transport, Edge, Pedestal) code <cit.>. STEP enables self-consistent iteration between various OMFIT <cit.> modules. We utilize CHEASE <cit.> for equilibrium calculations, TGYRO <cit.> with TGLF <cit.> for transport, and BALOO <cit.> for infinite-n ballooning instability. It is of note that bootstrap current was not included in the simulations in this work except in the discussion of ballooning stability in the edge in subsection <ref>, but the inclusion of bootstrap current is not expected to alter the results of the core significantly. TORIC <cit.> was used on the MANTA core scenario to provide heat deposition profiles which were passed to the STEP workflow through the CHEF <cit.> module. A diagram of this workflow is shown in figure <ref>. TORIC was not used to solve for heating in the high-volume case. Heating for the high-volume case was instead copied from the heat deposition profiles solved for in MANTA and scaled as needed. TGYRO is a physics-based transport solver that uses NEO for neoclassical transport calculations and TGLF for turbulent transport calculations with moderate computational cost.

The modules mentioned thus far are primarily made for and/or trained on PT plasmas, and as such may not capture all NT effects. While the geometry change is taken into account in terms of surface area and volume, we note that TGLF is not a full gyrokinetic model, so there may be gyrokinetic effects of NT shaping that are not accounted for. However, gyrokinetic analysis done in other work has implied that the improved confinement of NT is likely due to gradients beyond ρ = 0.8 <cit.>, where TGYRO evolution does not extend to in this work. Gradients in this region are subject to stronger triangularity shaping, as toroidal plasmas become more circular farther from the edge. For example, in the high-field base case, triangularity increases from about -0.5 to -0.3 from ρ = 1.0 to ρ = 0.8. Even so, reference <cit.> found that even though the effect is stronger at high radii, NT exhibits reduced transport over PT at low radii as well. Additionally, TGYRO/TGLF has been shown to reliably recreate DIII-D NT shots with various saturation rules <cit.>, increasing the confidence with which we can apply these models to a reactor concept. However, we primarily use the SAT-2 saturation rule because it has been shown to better match experiment than SAT-0 at high powers and includes geometry effects <cit.>. Unless otherwise stated, all simulations in this work evolve the temperature profiles from ρ = 0.8 to ρ = 0 with TGYRO/TGLF, using the ICRH heating profiles that were simulated for MANTA in TORIC with CQL3D<cit.> (with scaling as needed) and holding all density profiles constant with electron density profiles at the Angioni peaking pkAngioni <cit.> and impurity density profiles scaled with electron density. Krypton is also included at a fraction of 0.001 as a seeded impurity with an otherwise 50/50 D/T fuel mix in all simulations unless otherwise noted. The use of impurities is elaborated on in section <ref>. 

Only simulations fully converged in TGYRO are shown in this work. The definition of “full" convergence for purposes of this work is described in <ref>.  





 §.§ Density psuedo-evolution using density peaking scaling law


Most currently operating tokamaks are fueled using neutral beam injection (NBI). NBI allows on-axis particle sourcing, which can increase density peaking and overall aides the tokamak's performance <cit.>. However, due to the relatively large size and high-field of a reactor-class tokamak, fueling by NBI may not be practical in an FPP <cit.>. Instead, the particle sources in a reactor are likely to be confined to the edge region, outside of ρ = 0.8. 

To assess the effect on density peaking pk of a particle source in this region, we scan a Gaussian electron source from ρ = 0.7 to ρ = 0.9 in a MANTA-like core scenario. In these scans, TGYRO/TGLF with SAT-0 converges the density profile to one that has no more than a 2% deviation from the peaking predicted by the Angioni scaling <cit.> pkAngioni, as shown in figure <ref>. SAT-0 was used in this case because it is known to be easier to converge particle flux and heat flux simultaneously. Given the results of our particle source scans we suspect that TGYRO/TGLF would predict a similar peaking as is expected from the Angioni scaling. Though it is a challenge to converge particle flux and heat flux simultaneously with SAT-2, we chose to use SAT-2 in the remainder of this work because it includes additional geometry effects that SAT-0 and SAT-1 do not. We therefore use the pkAngioni prediction for density and TGYRO/TGLF with SAT-2 for temperature profile prediction in all subsequent parameter scans for ease. This density peaking prediction is dependent on collisionality, NBI source, and normalized plasma pressure β. As both the collisionality and β are dependent on temperature, the peaking prediction changes as the temperature profile evolves in TGYRO. As a result, we pseudo-evolve the density between iterations by adjusting the core density after running TGYRO/TGLF and repeating until we find a converged solution with density peaking within the Angioni prediction. The reference “Angioni" profile in figure <ref> is that which was predicted using this procedure. 



It is of note that the Angioni scaling only uses a dataset of AUG and JET H-mode observations. However, it was found that collisionality was the most statistically significant parameter in the analysis <cit.>. The normalized collisionality of MANTA is 0.4 at ρ = 0.9 and drops monotonically to 0.02 in the core. Note that the collisionality is high in the edge, which is not allowed in PT H-mode reactor core scenarios. However, NT does not have the same current-limiting pedestal demands as PT H-mode <cit.>.













§ IMPACT OF THE NT PRESSURE BOUNDARY CONDITION VARIATION ON FUSION PERFORMANCE


A significant remaining uncertainty in NT performance prediction is establishing the proper edge pressure condition. This is because the NT edge region is unique in that it is not a true “L-mode" or “H-mode" edge. Experiments have shown that NT can have steeper gradients than PT L-mode in the region outside of ψ_N = 0.9, forming a small “pedestal" while remaining ELM-free <cit.>. However, they are still often able to recover the same pressure as PT H-mode in the core <cit.>. Because the NT edge has yet to be characterized to the extent of the PT H-mode edge, which still itself retains significant uncertainty during extrapolations to an FPP, there is more uncertainty around what pressures and pressure gradients could potentially be obtained in an NT reactor edge. 



To capture effects related to the edge boundary condition, we employ a “brute force" characterization of the edge in both the high-field and high-volume configurations. This characterization is performed by scanning over four electron temperature and four electron density boundary conditions at ρ = 0.8, for a total of 16 simulations. We use the psuedo-density evolving method described in section <ref> and evolve temperature from ρ = 0 to ρ = 0.8 with TGYRO/TGLF and SAT-2. The resultant electron density and converged temperature profiles are shown in figure <ref> with boundary condition values shown by the dotted gray lines. These values can also be seen from the x and y axis in figure <ref>. In figure <ref> and <ref>, H98y2 is given by the colorbar. The green region in figure <ref> is that which was evolved in TGYRO/TGLF. Note that as the edge value increases in both electron temperature and density, H98y2 increases as well. The highest pressure boundary condition for the high-field configuration tested in this work exhibited a confinement time 20% over that expected from the τ98y2 scaling <cit.>, while the highest for the high-volume configuration was 5% below. Attaining H98y2 > 1 for both cases is difficult at low pressure boundary conditions.



In both the high-field and high-volume cases, fusion performance is seen to depend heavily on the pressure boundary condition at ρ = 0.8. This can be seen in figures <ref> and <ref> for the high-field and high-volume cases, respectively. The same conclusion for high-field was drawn in reference <cit.> from high fidelity modeling of PT L-mode operation in SPARC. Investigation of the edge physics in tokamak plasmas requires further work in both cases.



In both figure <ref> and figure <ref>, the top plot shows Pfus increasing with increased electron density at ρ = 0.8 (ne,0.8) for three distinct temperature values at ρ = 0.8 (Te,0.8). Note that Pfus also increases with increased Te,0.8. The bottom plot in both figures shows the same trend for PSOL for both cases except for at the lowest temperature in the high-volume case (blue line in figure  <ref>). All points shown are converged using the psuedo-evolving scheme for density and TYGRO/TGLF for temperature as described in section <ref>. The green circled points in figure <ref> and figure <ref> indicate simulations in which the Greenwald fraction (calculated with volume-averaged density) exceeded unity. The red circled points in figure <ref> are those in which the infinite-n ballooning stability limit was exceeded in the region from ρ = 0.8 to ρ = 1.0 for a pedestal width of 0.1 in ψN. The scans over ne,0.8 are at much smaller values for the high-volume case than for the high-field case due to exceeding the Greenwald limit at higher densities. This results in higher PSOL for the high-volume case as well, due to the lower separatrix density required. Therefore, a significant benefit of the high-field case from a power handling standpoint is the ability to employ higher separatrix density over higher volume cases without exceeding the Greenwald limit. This is even with the high-volume case employing higher IP, because the minor radius is more than twice that of the high-field case, which reduces the Greenwald limit significantly. Note also that the high-field case displays higher Pfus, but also converged with higher Te,0.8 than the high-volume case. This is once again likely attributable to the higher ne,0.8 employed in the high-field case. The mutual dependence of edge temperature and density and their significant effect on Pfus and PSOL motivates additional work in characterizing the NT edge boundary condition in future experiments and modeling.





 §.§ Ballooning stability in the region outside of ρ = 0.8



In addition to the impact of pressure boundary condition on core fusion performance, another uncertainty related the edge in an NT FPP is the specific conditions under which an NT plasma remains ELM-free. The leading experimental explanation of ELM-free performance at significant negative triangularity is the closure of the second infinite-n ballooning stability region, which restricts pressure gradient growth <cit.>. In these experiments, where pressure gradients are confined to the first stability region, no ELMs are observed. However, it is of note that DIII-D NT experiments have also suggested that there is a gradient limiting mechanism that precedes ballooning instability <cit.>. Thus ballooning instability is likely only sufficient as an upper bound on the pressure gradient; it is possible that the pressure gradient will be even lower, which would lower global performance. This further justifies the scans of the edge pressure boundary condition performed earlier in this section. In this work, we use ballooning stability as a proxy for access to the ELM-free state, assuming that if the normalized pressure gradient remains in the first stability region the plasma will not generate ELMs, as suggested by reference <cit.>. 
Gradients at the core-edge boundary (ρ = 0.8) in this work were maintained well below the infinite-n ballooning stability limit in all scans.   



An example of the infinite-n ballooning stability calculated by BALOO is given in figure <ref> for the high-field case with fGr≈ 1 and H98y2≈ 1 and a pedestal width of Δped = 0.1. Here the red region indicates the ballooning unstable region while the normalized pressure gradient and the total pressure in the edge are shown by the blue profiles in figures <ref> and <ref>, respectively. Since edge profile prediction is outside the scope of this work, we extrapolate from the end of the TGYRO evolution region (ρ=0.8) to ρ = 1 with an H-mode-like  pedestal shape. In figure <ref> the pressure gradient remains in the first stability region even at the steepest point. Similar analysis of the high-volume case reveals a significantly larger gap between the first stability limit and the normalized pressure gradient, as the edge is instead limited by the imposed f_Gr=1 constraint. Should Greenwald fractions greater than unity be achievable in an NT FPP, it may be possible to raise the edge pressure further in machines with larger Rmaj before destabilizing ideal ballooning modes <cit.>. 

For ease of comparison with calculations presented in <cit.>, the bootstrap current is omitted from figure <ref>. However, we note that inclusion of this effect can impact the local magnetic shear in regions of strong pressure gradients, potentially leading to a reduction of the maximum ballooning-stable pressure gradient and adding uncertainty dependent on the choice of bootstrap model. As such, results presented in this work should be treated as a theoretical upper limit on the edge pressure gradient rather than a precise prediction, as mentioned above. The reduction of the maximum achievable edge pressure gradient with the inclusion of bootstrap current is less pronounced at larger Rmaj, again suggesting that an increase in aspect ratio could be used to recover edge performance that may otherwise be limited by ballooning stability <cit.>; the uncertainty in ballooning stability due to bootstrap current is more important at low aspect ratio. 



Beyond the conventional presentation of ballooning stability presented in figure <ref>, it can be informative to examine a scan of the pedestal width and height for each given case. This is especially true in NT scenarios where the pedestal is not well-modeled using H-mode pedestal predictions like EPED <cit.> and a physics-based predictive model capable of describing the pedestal width and height has yet to be developed. In figure <ref>, we inspect the stability of potential NT edge pressure gradients resulting from various pedestal widths in the region from ρ = 0.8 to ρ = 1.0 by running gk_ped <cit.> on the high-field case. The gk_ped code is a linear gyrokinetic threshold model that generates self-consistent equilibria (including the bootstrap current, calculated with the analytical formulae from <cit.>) that vary in pedestal width, pedestal density, and pedestal temperature and then solves for the ballooning critical pedestal via the BALOO stability formalism at each point. Here we define the ballooning critical pedestal as the profile form that achieves the highest pressure gradient at a particular pedestal width before becoming ideal ballooning unstable. In figure <ref>, which shows the calculations without the inclusion of the bootstrap current, the ballooning critical pedestal is described by the relationship

    Δped = 0.35βθ,ped^1.03,

where Δped is the pedestal width and βθ,ped is the normalized poloidal pedestal pressure. The color bar in figure <ref> represents the fraction of radial locations of the pedestal half-width that are ballooning unstable. We note that, for the high-volume cases presented, the modeled scenarios presented in this work feature pedestals that lie below the ballooning critical pedestal calculated by gk_ped.

Equation <ref> can be used to describe an upper ballooning-stable limit on the NT pedestal height as a function of pedestal width, as characterized by the traditional H-mode-like  pedestal shape. However, as seen in the comparison to figure <ref>, which includes an analytic model for the bootstrap current in the edge region, this ballooning critical pedestal may over-predict the edge gradient in NT FPPs. This discrepancy is potentially compounded by observations on DIII-D that the actual edge gradient in NT plasmas lies somewhere below the infinite-n stability limit <cit.>, suggesting that an additional model is required for accurate prediction of the NT edge profile. However, we note that the location of high-field equilibrium point in the instability region in figure <ref> does not invalidate the value of pressure employed at ρ = 0.8 for the core profile modeling presented in the work. Indeed, it is possible to achieve a similar pedestal height as in figure <ref> while remaining stable to ballooning modes by going to a larger pedestal width. This highlights the need to develop a constraint for the NT pedestal width similar to the EPED model in H-mode scenarios, as the ballooning stability itself cannot fully constrain a pressure boundary condition at ρ = 0.8.

Because the region between ρ = 0.8 and ρ = 1.0 is beyond the TGYRO evolution in this work, any further analysis in this region is out of scope. However, figures <ref> and <ref> suggest that NT boundary condition at ρ = 0.8 may reach high pressures and pressure gradients, beyond those expected in a PT L-mode, without becoming ballooning unstable. This is supported by NT experiments on DIII-D that have observed increased edge pressure and pressure gradients in the NT edge without ELMs <cit.>, though we note that edge transport barrier mechanics and turbulence suppression remain an open area of research in NT. Figure <ref> illustrates that even in a high-field case with H98y2≈ 1 and fGr≈ 1, there are many ballooning-stable pedestal shapes that may be possible in NT. In particular, it is possible to achieve the same temperature and density boundary conditions at ρ=0.8 with reduced pedestal gradients by increasing the width of the NT pedestal, which is not constrained in this work. Thus the profile shown in <ref> is just one of many possibilities of an NT pedestal shape given the pressure at ρ = 0.8 and we encourage further experimental endeavours to characterize this region. In particular, physics-based constraints on the pedestal width or on the functional form of the pedestal would be valuable for improved predictive capabilities for NT FPPs, as they would enable a significant reduction in the parameter space available for stability codes like gk_ped.


§ RELATIVE EFFECT OF TRIANGULARITY AND MAJOR RADIUS ON FUSION PERFORMANCE


The importance of the boundary pressure at ρ = 0.8 on plasma performance is also seen when exploring changes in fusion performance due to geometry. While scanning triangularity and major radius, we found that changes in fusion power density due to Rmaj did not dominate over changes due to boundary electron pressure pe,0.8, as will be shown.

The motivation for studying a high-volume case is that Pfus increases with increasing volume <cit.>, so larger R_maj can be employed for lower field devices to achieve performance comparable with more compact high-field devices. In a high-field device like MANTA, compactness is prioritized because size is expected to be a major cost driver. However, there is an additional benefit of large R_maj in that it allows for a larger central solenoid, enabling increased flux swings and longer pulse lengths. Thus R_maj is an import parameter to optimize for FPP performance. In this work, we do not discuss the potential engineering benefits of increased Rmaj in an NT FPP, only the effect on transport and core performance. 








Volume generally decreases with increasingly negative triangularity, which affects fusion power output. NT has shown robust avoidance of ELMs when δ < δ_crit in experiment, where δ_crit is a critical triangularity that is device-dependent, around δ_crit∼-0.15 on DIII-D <cit.>. This δ_crit is unknown experimentally in a MANTA-like device, so it is safer to assume stronger negative triangularity to ensure the plasma is robustly ELM-avoidant. However, plasmas with stronger NT are more vertically unstable, so δ is also a parameter ripe for optimization <cit.>. 



Interestingly, volume does not directly correlate with negative triangularity in these scans, but peaks at around δ = -0.3. This is likely due to squareness not being held constant across all scans. We have plotted the last closed flux surface of the equilibrium shapes scanned in this section in figure <ref> for reference. Each color is at a different major radius while the transparency of the contours increases with increasing δ.



To illustrate the impact of Rmaj and δ on fusion performance in an NT scenario, the fusion power density is plotted against δ in figure <ref> with each distinctly colored line representing a different Rmaj. We do not focus on the obvious increase in Pfus with volume, instead investigating fusion power density to discover any underlying transport effects attributable to the change in shape. The dotted lines contain simulation points in which equilibria were initialized with the same Bt, IP, Paux, target aminor, temperature profiles, and density profiles, with only δ and Rmaj changing between each simulation. Note that this does lead to changes in magnetic shear and safety factor q as well, so there is potential performance optimization to be done at distinct δ and Rmaj owing to global stability considerations. For example, q_95 = 2 for the Rmaj = 6 m cases here, which is the lower stability limit. Additional scans that could increase IP to maintain q_95 constant at each Rmaj could expand the Rmaj optimization picture and are left to future work and should be combined with dedicated MHD stability modeling. Also note that scans of the shaping parameters lead to greater variation between the equilibrium in each simulation than in the Te,0.8 and ne,0.8 scans presented in section <ref>, leading to increased variation in equilibrium parameters that are targeted as inputs in CHEASE. In any case, an increase in the fusion power density with less negative triangularity is evident. This suggests that the optimal triangularity for operation of a NT FPP may be that which is just negative enough to avoid ELMs - further decrease in δ may result in a general decrease in the fusion power density.

In figure <ref>, the electron pressure at ρ = 0.8 is plotted against δ, highlighting some variation in edge pressure within each δ scan. In all TGYRO simulations in this work, the edge pressure at ρ = 0.9 was fixed while the pressure from ρ = 0 to ρ = 0.8 was allowed to evolve until convergence was met, resulting in some variation in edge pressure at ρ = 0.8. We tested the role of edge pressure on these results by performing two additional scans over δ at constant Rmaj = 5 m. The first is shown by the solid blue triangle markers, with Te,0.9 increased by 1 keV with all other parameters the same as the corresponding dashed Rmaj = 5 m line in figure <ref>(a). The second is shown by the solid blue square markers, with Te,0.9 decreased by 1 keV. When decreasing Te,0.9, it was more difficult to converge TGYRO in the edge for these cases. The Te,0.9 + 1 keV and Te,0.9 - 1 keV scans resulted in fusion power density changing more drastically than the change from Rmaj = 4 m to Rmaj = 6 m, but still displayed the same increase in fusion power density from increasing triangularity. Thus edge pressure has a significantly more profound effect on fusion power density than Rmaj, though there is still the benefit of increased Pfus from the increased volume at higher Rmaj. 

Though there is increased vertical stability at weaker triangularity <cit.> and increased Pfus density at weaker triangularity, operating below a critical triangularity δcrit will be necessary for ELM avoidance <cit.>. For a MANTA-like device δcrit cannot yet be experimentally verified, but it is clear that it would be beneficial to operate as close to δcrit as possible while remaining ELM-free. 




§ THE EFFECT OF THE TEMPERATURE EDGE CONDITION ON HEATING REQUIREMENTS



Given the sensitivity of core performance on the temperature and density edge condition, establishing a set of reliable actuators and controls for these parameters is paramount to the successful design of an NT FPP. One potential path to control the temperature edge is to leverage auxiliary heating power Paux. For MANTA, the full wave code TORIC was used with CQL3D<cit.> to determine 1D power deposition profiles for both ions and electrons from ICRH heating such that the total Paux≈ 40 MW<cit.>. In Paux scans in this work, heating was not calculated self-consistently but instead was simply scaled from the MANTA power deposition profiles for both the high-field and high-volume cases. In figure <ref>, scans of Te,0.8 and Paux on the high-field and high-volume core are shown. The colorbar gives Pfus in mega-watts, and selected tokamak FPP operating points are indicated by blue stars. Their corresponding parameters are given in table <ref>.



It is clear in both cases that Pfus increases with Te,0.8. Though Pfus also increases with Paux, it is not as pronounced as the change due to Te,0.8. Note that we cannot claim the level of Paux that will be required to maintain a certain Te,0.8, as the relationship between Paux and Te,0.8 is ultimately governed by edge physics. Instead, the two are varied independently in figure <ref> to scope potential Paux and Te,0.8 combinations and their corresponding Pfus. The relatively weak dependence of Pfus on Paux suggests the possibility of higher gain solutions at a given Te,0.8 by going to lower Paux. We ultimately find that Te,0.8 is more important than Paux in determining fusion performance, which once again motivates further investigation of the physics governing the NT edge.

Similar solutions can be found in both the high-field and high-volume cases at various Paux, but are harder to converge at Paux < 10 MW and Te,0.8 < 5.5 keV in the high-field case, requiring higher Te,0.8 than the high-volume case to reach the same Pfus. Note that ARC and EU-DEMO, representing the high-field and high-volume path for PT H-mode FPPs, respectively, have approximately the same Te,0.8 and Paux values. However, ARC has a fusion power of 525 MW while EU-DEMO has a fusion power of 1800 MW, due predominantly to its larger volume. If we compare the two points in grey squares on figure <ref>, they are also at approximately the same Te,0.8 and Paux value, but exhibit Pfus of 786 MW for the high-field case and 940 MW for the high-volume case. While there are likely other parameters at play constituting the difference between ARC and EU-DEMO performance, NT appears to approximately follow a similar trend to PT H-mode between the high-field and high-volume approaches in this case. Note also that the zone in which the converged simulations lie also differs between the high-field and high-volume case in NT. This is due to the density difference between cases, with the high-volume case able to converge at a lower ne,0.8. Referring to table <ref>, ne,0.8 for the high-field case is more than three times that of the high-volume case. Thus, a higher Te,0.8 is required for power balance in TGYRO. The high-field case also displays higher sensitivity to Te,0.8 than the high-volume case. 

The expected Te,0.8 and Paux needed for comparable Pfus to leading PT H-mode FPPs is comparable to that seen in these devices. ARC is the lowest shown here at Pfus = 525 MW and ARIES-ACT2 the highest at Pfus = 2600 MW, but note that ARIES-ACT2 has significantly higher Te,0.8 which we've seen to be very influential on Pfus.




§ IMPURITY ANALYSIS



Given the primary benefit of operating in NT is enhanced confinement without ELMs, we are interested in robustly avoiding H-mode, contrary to most other FPP concept regimes. Discarding the requirement to maintain P_SOL above the L-H power threshold grants freedom to use seeded impurities to radiate heat in the edge, lowering PSOL and thereby reducing divertor heat loads. Noble gas impurity seeding in PT L-mode experiments on DIII-D has shown enhanced confinement with low P_SOL <cit.>, similarly to NT. Employing both NT shaping and seeded impurities allows more control over Pfus and PSOL and potentially easier divertor integration than with NT alone. In a reactor, high-Z noble gas impurities such as krypton and xenon are expected to be most useful given their ionization temperatures, because they primarily cause radiation in the edge, decreasing PSOL significantly with minimal effect on Pfus. 



In all simulations mentioned thus far in this work, we included only krypton impurities at a fraction of 0.001 unless otherwise noted. Note that all impurity density profiles in this work are set to scale with the electron density profile, which is a limitation of this work. These results are likely to be affected by impurity transport causing impurity density profiles to differ significantly from the electron density profile. Tungsten is the leading candidate for plasma facing component in FPPs and so tungsten ions will likely also be in a reactor-class plasma along with helium ash, so its effect should be more closely studied in future work on the subject. 

To better assess the impacts of impurities on fusion performance in an NT FPP, we plot power density profiles from line radiation in the high-field base case with three distinct impurity combinations in figure <ref>. In this figure, the power density profiles of the additional line radiation from including krypton (Kr) and tungsten (W) is plotted by the dashed lines in context of the alpha power density (solid lines) and auxiliary power density (dotted line). Note that the auxiliary power is the same for all three impurity combinations. The dilution effect of adding helium at a fraction of 0.02 in a D, T, and Kr mix results in a decrease in Pfus of only 121 MW, or about 9% of Pfus in the D, T, and Kr only mix. The power density profiles of a D, T, Kr, and He mix are not shown in figure <ref> because they overlap the D + T + Kr + He + W profiles. Including tungsten at a fraction of 1.5× 10^-5 does not significantly affect radiated power, and the Pfus difference between the D, T, and Kr profile and the D, T, Kr, He, and W profile in figure <ref> is only 127 MW, or ∼10%. Note that the presence of tungsten will be inescapable in any reactor class device using tungsten as a plasma facing material (the leading metal candidate) and that it elicits serious concern for radiative collapse due to its high atomic number. For example, the tungsten fraction used in EU-DEMO modeling is 10^-5 <cit.> and in any reactor-class tokamak the upper limit on tungsten fraction is likely to be on the order of several 10^-5 <cit.>. H-mode has the benefit that ELMs flush impurities out of the core <cit.>, but the tungsten fraction limit is also set by maintaining PSOL > PL-H, a limit that NT does not need to adhere to. Additionally, ELMs can cause sputtering, leading to increased tungsten influx <cit.>. While impurity transport in NT is an area of ongoing research, preliminary analysis on diverted DIII-D NT experiments suggests rapid impurity transport. This is indicated by hollow impurity profiles and lower Zeff in NT plasmas than PT plasmas with similar confinement factors <cit.>. 



In figure <ref>, we plot Pfus against krypton impurity fraction fKr in the high-field case. The circles plotted are at 40 MW of auxiliary power while the triangles are at 20 MW of auxiliary power. The color bar gives PSOL. Figure <ref> shows that the same downward trend in PSOL versus fKr change is seen at Paux = 40 MW and when Paux =20 MW. This is promising for prospective direct PSOL control using impurity seeding. However, figure <ref> also shows that a small change in fKr results in a large change in PSOL, and research into whether this level of impurity control in a reactor will be possible is ongoing<cit.>. Additionally, figure <ref> shows a peaked trend in Pfus, indicating potential for Pfus optimization via impurity fraction. This improvement in Pfus from certain levels of additional impurity seeding while simultaneously decreasing PSOL has been observed in experiment in NT on DIII-D <cit.> and encourages further study of this phenomenon for implementation in reactor modeling scenarios. For a given fKr, we also see that Pfus is higher when Paux = 40 MW than when Paux=20 MW, indicating potential for Pfus control via input power. At each auxiliary power, H98y2 remained approximately constant over the scan of fKr, indicating that improvement in confinement from impurities may overcome the decrease in Pfus from dilution.

While the inclusion of radiative impurities could lead to benefits in the core performance and in the power handling properties of NT FPPs, we note that we do not include in this work any decrease in the edge temperature resulting from significant radiation just inside of the separatrix. The impurities studied above were chosen to consolidate radiation in the core region <cit.>, but any drop in the temperature boundary condition will lead to a decrease in fusion power, as discussed in section <ref>. As such, proper characterization of the role of impurities with a full impurity transport code that extends out into the SOL is needed for more robust NT FPP design, and should be the subject of future work.



§ CONCLUSION



These studies demonstrate the feasibility and flexibility of the NT reactor concept. Even with edge pressure conditions low enough to be infinite-n ballooning stable, a variety of operating points exist in which FPP-relevant fusion power is possible (∼400-500 MW for a MANTA-like device with 40 MW input power <cit.>), with opportunities for increased fusion gain by decreasing input power, increasing fueling, or optimizing seeded impurity fraction. Due to the lack of a requirement to maintain H-mode, NT also grants the freedom to increase the seeded impurity fraction. This can increase radiation in the edge to bring scrape-off layer power down to acceptable levels for simple divertor integration (<40 MW for a MANTA-like device with a separatrix density of 0.9 × 10^20/m^3 <cit.>). 

The performance of an NT reactor will be heavily dependent on the edge condition. This dependency is not unique to NT but is perhaps of greater consequence in NT designs than in PT designs due to the present limitations in modeling the NT edge. Changing the edge temperature by 1 keV was found to have a more profound effect on the fusion power density than changing the major radius by 2 meters in a high-field case, though we note that this is without adjusting scans to account for changes in global stabilization parameters such as q95. Additionally, we found the temperature at ρ = 0.8 to have a significantly larger effect on fusion power than auxiliary heating in both the high-field and high-volume cases, though the high-field case exhibited stronger dependence than the high-volume case. We reiterate that the density and temperature chosen for scans over other parameters in this work were to satisfy H98y2≈ 1 and fGr≈ 1 while meeting temperature convergence in TGYRO with density at the Angioni peaking. We have not explored the extent to which these edge parameters can be accessed, and rely heavily on experimental observations to inform this work.

Though the effect of triangularity on fusion power density is minimal, fusion power density is found to consistently increase with less negative triangularity, suggesting that the optimal triangularity for an NT FPP may be that which is just negative enough to sustain ELM-free operation. Targeting this minimum negative triangularity will also have the added benefit of reducing vertical instability concerns <cit.>. In this work, the core was modeled from ρ = 0 to ρ = 0.8, so future work should prioritize high fidelity transport modeling in the edge where the triangularity is strongest to identify any benefits attributable directly to the negative triangularity geometry, like that seen in reference <cit.>. 

For a high-field case with fGr≈ 1 and H98y2≈ 1, we demonstrated a possible model for extrapolation to ρ = 1.0 that relates pedestal width to height based on the stability found from infinite-n ballooning models. Thus for a given edge pressure boundary condition, there are a variety of pedestal shapes that remain ballooning stable. No scans of the high-volume case exceeded the infinite-n ballooning stability boundary, so an increase in Rmaj could enable an increase in fusion power on two fronts: from the increased volume as well as the increased pressure boundary condition <cit.>. 

The main advantages of NT come from the improved pressure edge condition over PT L-mode supported by experiment and the absence of a requirement to remain in H-mode allowing flexibility in impurity and heating requirements. We found that the temperature boundary condition and auxiliary power needed to reach fusion powers between 500-900 MW is below or on the order of representative PT H-mode FPP concepts. Future work to develop a physics-based model for the behavior of an NT reactor edge must be prioritized to more accurately predict fusion performance, though we have shown in this work that feasible core operation can be attained at a variety of edge conditions. To further evaluate the feasibility of integrating a given operating point with simple divertor operation, whole-device modeling like that done in reference <cit.> could illuminate the engineering trade-offs attributable to the potential plasma core trade-offs described here, such as pressure boundary condition, impurity content, geometry, and heating. Ultimately, there are likely a variety of core operating points for a NT tokamak FPP with competitive fusion power to leading PT H-mode concepts even at relatively high radiation, at low input power, and with a ballooning-stable pressure boundary and significantly reduced scrape-off layer power. 



§ ACKNOWLEDGEMENTS

Most calculations performed in this study were completed through the OMFIT framework <cit.>, and the generated equilibria are available upon request. This work was supported by the U.S. Department of Energy, Office of Science, Office of Fusion Energy Sciences under Award DE-SC0022272.





§ REFERENCES


[heading=none]





§ FULL TGYRO CONVERGENCE


We define “full" convergence to be met when the residual between the total flux calculated in the TGYRO/TGLF model (ftot) and the target flux from power balance (ftar) for each point between ρ = 0.35 and ρ = 0.8 is less than or equal to 0.02. Here the residual is defined as 

    (ftot-ftar)^2/ftot^2+ftar^2.

Converging from a residual of 0.02 to a residual of 0.00 resulted in less than a 5% change in Pfus in a few representative cases, so pushing convergence past this point is not expected to significantly change the qualitative results presented in this work. Flux matching between ρ = 0.0 and ρ = 0.35 is challenging but has ultimately been shown to have a marginal effect on output fusion power due in part to the relatively small plasma volume in the core compared to the edge <cit.>.
