Core performance predictions in projected SPARC first-campaign plasmas with nonlinear CGYRO

This work characterizes the core transport physics of SPARC early-campaign plasmas using the PORTALS-CGYRO framework. Empirical modeling of SPARC plasmas with L-mode confinement indicates an ample window of breakeven (Q>1) without the need of H-mode operation. Extensive modeling of multi-channel (electron energy, ion energy and electron particle) flux-matched conditions with the nonlinear CGYRO code for turbulent transport coupled to the macroscopic plasma evolution using PORTALS reveal that the maximum fusion performance to be attained will be highly dependent on the near-edge pressure. Stiff core transport conditions are found, particularly when fusion gain approaches unity, and predicted density peaking is found to be in line with empirical databases of particle source-free H-modes. Impurity optimization is identified as a potential avenue to increase fusion performance while enabling core-edge integration. Extensive validation of the quasilinear TGLF model builds confidence in reduced-model predictions. The implications of projecting L-mode performance to high-performance and burning-plasma devices is discussed, together with the importance of predicting edge conditions.


Introduction
Upcoming burning plasma experiments in both ITER [2] and SPARC [3,4] are designed to achieve maximum fusion power and fusion gain (Q) in H-mode operation, with their baseline projections of Q = 10 (ITER Baseline Scenario, IBS) and Q = 11 (SPARC Primary Reference Discharge, PRD) necessitating a high energy confinement regime to sustain burning plasma temperatures at modest input power.
While operating plasmas in H-mode improves the energy confinement thanks to the formation of an edge transport barrier-thus reducing the power requirements to sustain the confined plasma energy-there are a number of unfortunate aspects that are associated with the formation of an edge pedestal.Operation in H-mode typically comes with associated edge localized modes (ELMs) that require either mitigation or suppression to avoid their detrimental effects on plasma facing components, which will be a pressing challenge for burning plasma experiments and fusion pilot plants [5].Even though there are variations to H-mode operation that are ELM-free [6], the predictive capabilities to project access to those regimes and to predict confinement times in reactor-relevant conditions are very limited.Similarly, I-mode regimes [7] (also characterized by an energy edge transport barrier, but without a particle transport barrier) that exhibit improved confinement relative to L-mode plasmas have limited predictive capabilities, as the underlying turbulence and transport mechanisms are subject of on-going research (e.g.[8,9]).Most predictive work -either with the use of empirical confinement laws such as τ 98y2 E [2] or edge pedestal models such as EPED [10]-is focused on type-I ELMy H-modes, thus only providing an upper bound to performance predictions for future devices, assuming that operation without ELMs requires, typically, a degradation of pedestal pressure from that coming from peeling-ballooning stability considerations [11,12].
While both H-mode and I-mode operation is associated with higher confinement levels and thus higher efficiency, every tokamak can operate in L-mode conditions.L-mode can be described as the foundational and more naturally occurring state of plasma confinement in tokamak devices.It is characterized by the lack of a significant edge transport barrier and has associated lower confinement levels.It is often assumed that the lower confinement in L-mode regimes is due to stronger turbulent transport, which is particularly important at the edge of the plasma, where shearing of turbulence is not strong enough (compared to H-mode regimes) to sustain a high cross-field gradient and transport barrier.
The upcoming SPARC tokamak will first operate with L-mode regimes, which are projected by empirical scalings to achieve breakeven (Q > 1) conditions.These high-performance L-modes are currently expected to be, arguably, the first time that Q > 1 is attained in magnetic confinement fusion.Because of this, the SPARC team is performing extensive investigations of several aspects of these L-mode plasmas [13][14][15][16], which focus on achieving optimal performance from the scenario.In this paper, we describe the progress so far in understanding the core turbulence characteristics of these high-performance, DT-fueled L-mode plasmas.Leveraging the PORTALS-CGYRO framework [1], here we present a comprehensive study of achievable core gradients and core performance in projected L-mode-like conditions in SPARC, and discuss the importance of edge conditions to attain fusion performance goals.

Background
SPARC [3,4] is a high-field (B = 12.2T ), medium-size (R = 1.85m), deuterium-tritium (DT) capable machine under construction by Commonwealth Fusion Systems in Devens, Massachusetts (USA), with construction to be completed by 2025 and first plasma in 2026.The accelerated timeline of SPARC is enabled thanks to novel high-temperature superconductor technology, which allows the operation at high magnetic field and a moderate size.The high magnetic field results in the capability to operate at high plasma current (I p ≤ 8.7M A) and high absolute density (⟨n e ⟩ ≈ 3 • 10 20 m −3 in the baseline H-mode regime), while remaining reasonably far from stability boundaries (q 95 ≈ 3.4, f G = n e /n G ≈ 0.37, β N ∼ 1.0) [3,4,12,17].These conditions are favourable for high energy confinement and fusion energy production, despite being smaller than a number of tokamaks that have been built and operated, such as JET [18], TFTR [19] and JT60-U [20].Projections of performance in SPARC indicate a maximum fusion gain in the range of Q = 8 − 11 [3,12,21] operating with a mix of 50-50 DT fuel, provided that the pedestal remains at EPED-predicted [11] levels -i.e. at ELMy H-mode levels-with enough ELM control to ensure machine survivability.Performance prediction in other regimes, such as L-mode, I-mode and other H-mode variations with somewhat degraded pedestal performance, is currently ongoing, although with less confident physics-based projections due to the lack of reliable predictive tools of edge and near-separatrix conditions.
Before targeting these high-performance H-mode discharges, however, SPARC will demonstrate breakeven (Q > 1) in a DT pulse as the conclusion of its first operational campaign.As well as demonstrating breakeven as soon as possible, the first campaign will need to minimize total neutron production so as to not preclude vessel entry for maintenance and upgrades.With these goals in mind, the SPARC team investigated several high-performance L-modes as potential breakeven candidate scenarios.
L-modes are desirable for several reasons.Unlike H-modes, L-modes do not have an input power threshold.Therefore -if SPARC L-modes are found to have sufficient confinement-Q > 1 can be demonstrated with less fusion power and neutron production, while also reducing the number of ion cyclotron resonance heating (ICRH) antennas which must be commissioned to fullpower in the first campaign.L-modes are also naturally ELM-free, which both avoids the need to commission the ELM mitigation system and reduces the risk of damage to plasma facing components (PFCs).L-modes are also more reproducible, have improved density control and less impurity accumulation Performance in SPARC first campaign than H-modes.However, L-modes also have significantly lower energy confinement than H-modes -and so studying whether we can achieve sufficient performance to reach Q > 1 in L-modes is crucially important.
Initial low-fidelity scoping of the SPARC L-mode operational space using a POPCON analysis1 [22] indicated that Q > 1 can be achieved for a range of densities and heating powers, assuming an energy confinement given by the ITER-89P scaling [23].To reduce the time required for scenario development, a family of candidate scenarios was developed for parameters which are relaxed slightly relative to the PRD.These scenarios are for single-null geometries2 at slightly reduced elongation (areal elongation κ A = 1.7 vs κ A = 1.75 in the PRD) and plasma current (I p = 8.5MA vs I p = 8.7MA in the PRD), for densities ranging from ⟨n e ⟩ ≈ 1.0 • 10 20 m −3 (near the predicted LOC-SOC transition [24]) to ⟨n e ⟩ ≈ 2.1 • 10 20 m −3 [13].These scenarios were projected by empirical scalings to achieve Q > 1 while staying within operational limits such as those for disruption avoidance and PFC protection [14].This initial scoping was used to guide a series of higher-fidelity studies, investigating physics such as RF coupling and power exhaust.In this paper, we present high-fidelity core transport modelling for these high-performance L-modes -investigating in unprecedented detail which of the candidate scenarios has the highest probability of achieving Q > 1, as well as suggesting further optimization of the scenario to maximize performance.

Physics-based modeling techniques
To study core performance, this work leverages high-fidelity simulations using PORTALS-CGYRO and medium-fidelity simulations using PORTALS-TGLF.PORTALS [1] is a novel, open-source [25] toolset developed to accelerate the numerical convergence of time-independent transport solvers.The speedup achieved by the utilization of surrogate-based optimization techniques [26] in PORTALS has enabled some of the highest fidelity predictions of core transport to date.Direct nonlinear gyrokinetic simulations with the CGYRO [27] code have been used to predict steady-state (multi-channel flux-matched) plasmas in SPARC PRD [21], ITER IBS [28], DIII-D ITER Similar Shape [1,29], ARC [30] and JET H, D and T Ohmic plasmas [31,32].
The quality of PORTALS predictions depends on the model employed to simulate gradient-driven turbulent transport.Drift-wave-type turbulence is the dominant core transport mechanism in the plasmas of interest, and is modeled here using the GPU-accelerated CGYRO [27] code.CGYRO is an initialvalue, Eulerian, multi-species, spectral gyrokinetic solver designed specifically for electromagnetic, collisional, and multi-scale tokamak turbulence simulations.It solves the coupled gyrokinetic-Maxwell system of equations in the local δf limit [33] and is optimized for scalability and performance on modern high-performance computing architectures [34], including GPU-accelerated systems.CGYRO has been used to study a wide range of tokamak plasmas and is used in this work to predict the turbulent transport in flux-matched SPARC plasmas, to provide a high-fidelity prediction of the core temperatures and density.On the other hand, the Trapped Gyro-Landau Fluid (TGLF) quasilinear model of turbulent transport is used for the evaluation of performance at much reduced cost.TGLF solves a linear system of gyro-fluid equations [35] and uses a saturation rule [36] to estimate the transport levels associated with the typical drift-wave-type micro-instabilities in the core of tokamaks.TGLF is benchmarked against CGYRO predictions in this work to provide confidence in the reduced-model predictions.In both the PORTALS-CGYRO and PORTALS-TGLF simulations, neoclassical transport is calculated with the NEO code [37], but found to be negligible compared to turbulence as modeled by both CGYRO and TGLF (< 1.2% in all flux-matched radii).

Core turbulence modeling
To produce physics-based predictions (in contrast with the empirically-based or data-driven predictions used for the initial scoping) of core performance, we start with the definition of a potential parameter space of operation of early-campaign SPARC plasmas.With the goal of understanding the maximum fusion performance in SPARC L-modes, the scenarios presented in this section will be characterized by being operated with the full plasma current (I p = 8.7M A), magnetic field (B T = 12.2T ) and shaping (κ sep = 1.97, δ sep = 0.54).These parameters are assumed to be the same as for the PRD scenario, and assumed to provide an upper bound on expected performance, as plasma current, magnetic field and shaping parameters are expected to increase confinement time.Similarly, impurity mix is assumed to be the same as for the PRD: fuel dilution f DT = 0.85, tungsten concentration f W = 1.5 • 10 −5 , He-3 minority concentration f min = 0.05 and effective ion charge Z eff = 1.5 (for further information about these choices, the reader is referred to Refs [3,4,12]).Auxiliary input power absorbed from ICRH is considered to be below P ICRH < 15M W , consistent with commissioning plans for the machine and its actuator systems.For conservatism and because of the lack of reliable predictive models, no core fueling is assumed and therefore particle flux-matching is equivalent to finding the null-flux condition, Γ e = 0.

On parameter selection and assumptions
In the parameter scans presented in following sections, three parameters are chosen as free variables and varied to shed light on the expected performance of SPARC early-campaign plasmas: total input power (P in = P ICRH + P OH ), edge electron density (n edge ) and edge temperature (T edge ), the latter equal for ions and electrons.In the context of core transport predictions, "edge" refers to the ρ tor = 0.9 (square root of normalized toroidal flux) flux surface, which corresponds to r/a ≈ 0.943 (normalized half-width of midplane intercept) and ψ n ≈ 0.922 (normalized poloidal flux) for the plasmas of interest.We must note that the only true controllable parameter from an operational point of view is the ICRH launched power.Operational density is somewhat controllable, but the ability to reach desirable edge density levels will depend on the balance of fueling and transport near the separatrix and the scrape-off layer, which remains an open question and out of the scope of this work focused on core transport predictions.Similarly, edge temperature is a highly uncertain parameter that will be determined by the balance of transport, radiation losses and energy flows in the near-separatrix region.Due to the difficulties predicting the edge temperature, the approach taken here is to assume a potential range of edge temperatures (based on past experiments, as will be described in Section 3.2) and explore the sensitivity of fusion gain and fusion power on the achievable edge pressure.
In summary, the parameter space is defined by the following ranges: T edge = 0.9 − 2.3keV , n edge = 0.77 − 2.05 • 10 20 m −3 and P in = 8.3 − 18.3M W .There will be an exploratory, very high density case (n edge = 4.0 • 10 20 m −3 ) that was run at low T edge and high P in .To make this case survive, the tungsten concentration was halved, otherwise the plasma would radiatively collapse with the nominal tungsten content.These ranges can be visualized in Figure 4a) and b) (to be discussed in following sections), where the 12 different combinations of (n edge , T edge , P in ) are shown.
During the parameter scans, the plasma geometry will be assumed to be fixed, including the internal flux surfaces.This assumption is motivated by: 1) the prohibitive computational cost of performing predictions of scenarios with different geometry (as the PORTALS surrogate re-utilization technique would not be applicable), and 2) the expected small effect of internal equilibrium on transport, given the low normalized pressure, β N ≲ 0.5, and expected small differences in Shafranov shift for all cases studied.The power deposition is assumed to remain unchanged as well, both the spatial location of power delivered to ions and electrons, and the ICRH power partition between species.
To study the sensitivity of the ICRH power deposition so that we can make the parametric scans easier while retaining any relevant physics, scans using interpretive TRANSP [38] with various potential kinetic profiles were performed.TRANSP is run with the TORIC full-wave code [39] coupled with the bounceaveraged Fokker-Planck FPPMOD model [40] to calculate the self-consistent response of the minority species to the wave field.Figure 1 shows the results of TRANSP scans with various volume-averaged densities, volume-averaged temperatures, and temperature ratios.All the rest of parameters are kept fixed and are similar to the TRANSP simulations presented for the PRD [12], although interpretively.The results show that the relative amount of ICRH power absorbed by the ions remains within the range P ICRH,i /P ICRH ∼ 0.65 − 0.75 for the cases studied here, and therefore we considered it is a safe assumption to keep the power deposition fixed.Similarly, the change in spatial deposition is not significant, with the strongest effect at the inner most location simulated, r/a < 0.35 (motivated by the guidelines presented in Ref [1]).While the change in power deposition at that location is not completely negligible (∼ 25% increase in power when doubling the density from ⟨n e ⟩ = 1.0•10 20 m −3 to 2.0 • 10 20 m −3 ), the very stiff transport expected in the core will result in a very small change in the core temperature with variations in transported power.
We must note that the changes in both the ICRH partition and the deposition location are dominated by the changes in energy exchange power from the minority species, not the actual ICRH direct absorption by the minority or the Performance in SPARC first campaign plasma species.Because the temperature of the background minority species (He-3) is assumed to be equal to that of the background ions when no ICRH power is applied, and because both the electron and ion temperatures are fixed throughout the interpretive TRANSP simulations, there is artificial channeling of ICRH power from/to electrons to/from ions when the difference between T e and T i is too high for a given density choice.This suggests that an important future update to the PORTALS framework will be to allow for the auxiliary power delivery to evolve during the convergence process, as the temperature profiles are updated.
While the previous analysis has focused on the ICRH power, the plasmas of interest will also have a significant fraction of Ohmic heating.This work assumes fixed Ohmic power (with a representative database value of 3.3M W ), despite its potential to change with density, temperature and impurity content.It is because of this that the power scans are, in actuality, performed with respect to the total absorbed power, P in , which will be varied between 8.3M W (5M W ICRH + 3.3M W of Ohmic) and 18.3M W (15M W ICRH + 3.3M W of Ohmic).The limitation of running with fixed Ohmic power is a consequence of the capabilities of the current PORTALS framework, and will be improved in future work.

On the choice of boundary condition
As described earlier, both the edge density and edge temperature are assumed, imposed parameters in the core modeling.The edge density assumption is such that the operating volume-averaged Greenwald fraction is scanned in the range f G = 0.11 − 0.26.This work assumes that fueling systems will be capable of sustaining such edge density 3 .
The determination of edge temperature and its validity is more subtle.Transport predictions from the separatrix in L-modes, while possible [41,42], are more complex, uncertain and have a somewhat weaker physics basis, as L-mode edge transport remains an active area of research and experimentation and the validity of current transport models is still under investigation.Radiation and impurity effects, the dynamics of neutrals, the effect of rotation and sheared flows, charge exchange processes, and the influence of equilibrium geometry (strongly shaped plasmas in double-null) near the separatrix introduce complexities that are difficult to predict for breakeven plasma conditions.Furthermore, in addition to these, performance predictions that span the entire confined plasma must be coupled to scrape-off layer and divertor simulations and the difficulty to perform nonlinear gyrokinetic simulations at the edge (large normalized logarithmic gradients lead to small simulation timesteps and these locations require higher radial resolution) drove the decision of using r/a = 0.943 as the boundary condition for the kinetic profiles, as a tunable parameter to explore the operational space.We chose T edge = 1.6 ± 0.7 keV (as a reminder, edge in this work refers to r/a = 0.943) as the nominal assumption, motivated by a database study of non-H-mode plasmas in Alcator C-Mod.A set of ∼2000 discharges in Alcator C-Mod were chosen by their stationarity.Time windows of at least 250 ms were extracted from stored Alcator C-Mod data where the variation in the magnetic field, plasma current, and line-averaged density was less than 5% during the time window.If ICRF power was present, it was required to vary no more than 7.5% during the selected time window.These requirements ensured relatively steady discharge conditions while not being overly restrictive.From the set of discharges satisfying these criteria, only plasmas with q 95 < 4.5 and without significant edge density gradient were considered so as to only explore non-H-mode-like kinetic profiles (relatively few discharges with high density gradient and q 95 < 4.5 met the above stationarity conditions in Alcator C-Mod).Outliers (plasmas with unphysically high temperature gradients due to bad fitting or bad experimental data) were removed.As a result, a total of ∼600 discharges were used to derive distributions of achievable T e,edge (as measured by Thomson Scattering [43]) and associated a/L Te profiles.From this database, plasmas with T e,edge > 600 eV were chosen to be "high-performance" non-H-mode plasmas in Alcator C-Mod (this down-selection resulted in 61 discharges).The distribution (mean and standard deviation) of a/L Te in the edge region, r/a = 0.94 − 1.0, is shown in the top panel of Figure 2.
The bottom panel of Figure 2 shows the distribution of edge temperatures when gradients are radially integrated from an assumed separatrix temperature of 200eV at r/a = 1.0 (derived from a two-point model), which was chosen as a potential operational boundary condition for the SPARC plasmas of interest.This workflow therefore suggested -under the assumption of a constant normalized gradient scale length, which implies similarity in the transport physics in that region between Alcator C-Mod and SPARC-that T edge = 1.6 ± 0.7 keV (from the integration of "highest"-performing plasmas in the Alcator C-Mod database, T edge > 600 eV) is a reasonable choice for the edge temperature in this exploratory study.The validity of this edge assumption will be the subject of ongoing studies but the approach in this work was to attempt to identify a plausible set of L/I-mode edge profiles that were derived from existing experiment.We note however, that recent work by Hatch et al. suggests that the edge region of I-mode is consistent with an electron-temperaturegradient mode limited edge.Further investigation into these findings is the subject of ongoing work, and could provide a more physically robust manner for predicting near-edge profiles in SPARC.
We acknowledge the various caveats that this workflow introduces.First, the choice of normalized gradients scale lengths found in present-day tokamaks as a reference for SPARC could be subject to criticism, as it may as well be that SPARC (e.g. with higher plasma current and magnetic field among other parameters) could achieve higher gradients and therefore edge performance than those found in, e.g., Alcator C-Mod.Secondly, the database used to derive the edge temperature does not differentiate between L-or I-mode regimes.Power threshold and regime access considerations are not being taken into account.Finally, the choice of a separatrix temperature of 200eV (and equal to electrons and ions) is subject to variations depending on the operational regime (particularly the density and input power choice) and could be prone to optimization if it is indeed found to be a critical parameter to determine overall plasma performance.While all these are valid limitations for these predictions, the work presented here will discuss the implications of these choices and the potential avenues for improvement in future work (particularly as discussed in Section 6).

High-fidelity modeling setup
The PORTALS technique for steady-state modeling presents an important advantage (apart from higher nominal computational efficiency) to standard methods: the possibility to re-utilize surrogates when predicting scenarios with variations in certain input parameters.As explained in detail in Ref. [1], parameters that only affect the macroscopic evolution but not the local gradient-driven turbulence simulations (such as input power) and parameters whose effect is already included in nominal surrogates (such as density and temperature boundary conditions affecting core collisionality, temperature ratios and normalized pressure) are suitable for surrogate re-utilization.To this end, we exploited this capability to study 12 scenarios that required only a total of 71 profile evaluations to achieve flux-matched conditions in electron energy, ion energy and electron particle fluxes, thus resulting in the prediction of electron temperature, ion temperature and electron density.PORTALS is currently implemented as an electron-solver, thus evolving the electron density and adjusting the densities of the different ion species by keeping a constant, radially uniform concentration.This limits the capability of the framework to account for non-trace, multi-species particle transport (such as the effect of ash accumulation or radial variations of the D-T fraction).Future developments of PORTALS will include the capability to evolve the ion densities self-consistently.
Six radial locations are chosen to describe the logarithmic gradient profiles, r/a = [0.35,0.55, 0.75, 0.875, 0.9, 0.943].This choice is motivated by the database analysis presented in Ref. [1] that suggested that 5 locations are enough to study inductive, on-axis heated discharges, and that the inclusion of additional radial locations between r/a = 0.35 and the magnetic axis results in minimal changes to the total fusion power in SPARC.The additional location at r/a = 0.943 was added to ensure proper description of L-mode-like gradient scale length profiles, which are expected to have more structure than H-modelike profiles.During the process to achieve steady-state conditions, target fluxes were updated to account for the change in radiation losses (Bremsstrahlung, line and synchrotron radiation), collisional energy exchange between electrons and ions and alpha power to electrons and ions.The models employed for these target fluxes are identical to those in Ref. [21].
With the total of 6 radial locations per profile and with the goal of achieving convergence in the three described transport channels, this comprehensive study required 426 (71 profiles × 6 local radii) flux-tube nonlinear gyrokinetic simulations, which were performed with the CGYRO [27] code.Advances in GPU acceleration of CGYRO allowed us to perform this study with a total computational cost of ∼80,000 GPU hours (∼20,000 node hours), performed on the GPU partition of the NERSC Perlmutter supercomputer, equipped with NVIDIA A100 (40GB) GPUs.
All gyrokinetic simulations described in this paper were performed with the CGYRO code [27].CGYRO is a local, Eulerian gyrokinetic code that is optimized for multi-scale and near-edge simulation of tokamak discharges.The resolutions used in this work were relatively standard for the simulation of low-k (ion-scale) turbulence typically found in the core of tokamaks -such as microtearing (MTM), ion temperature gradient (ITG) and trapped electron mode (TEM) turbulence.High-k TEM and electron temperature gradient (ETG) mode contributions were ignored in this work motivated in part by the large (greater than 1) target values of Q i /Q e and the so-called fingerprint argument presented by Kotschereuther et al. [44].All simulations utilized relatively high fidelity physics including geometry represented by the Miller Extended Harmonic (MXH) formulation [45], fully electromagnetic turbulence (δϕ, δA ∥ , and δB ∥ ), Sugama collisions [46], and 4 gyrokinetic species: D, T, lumped impurity (Z = 4, A=8), and electrons.The lumped impurity species is meant to represent an average impurity that includes contributions from He-3, W, and other low-Z species.Rotation was not included in these simulations as no reliable estimate of the profile was available.Simulation domains varied slightly over the simulated radial regions (r/a = 0.35, 0.55, 0.75, 0.875, 0.9, 0.943) but nominal box sizes were [L x , L y ] ∼ 125 × 120ρ s utilizing 24 complex toroidal modes (N α ) and 512 radial modes (N r ) across the domain to represent turbulence from k θ ρ s = 0.0 to 1.22.Simulations also used 8 energy grid points (N u ), 24 pitch angles (N ξ ), and 24 points in theta (N θ ).Slightly higher radial and energy resolutions of approximately N r = 630 and N u = 12 were used at r/a = 0.9 and 0.942.
We must emphasize that the results presented here are not results of surrogate modeling of gyrokinetic turbulence, but of direct nonlinear CGYRO modeling.Surrogates are just used to accelerate multi-channel (Q e , Q i , Γ e ) convergence, but the steady-state predictions (T e , T i , n e ) are a full-physics result (of course, under the constraints of ion-scale nonlinear gyrokinetics as embedded in CGYRO) and must be interpreted as such.
Figure 3 shows an overview of the PORTALS process for the prediction of the 12 core scenarios.On the left, the multi-channel, multi-radius residual (i.e. the L 1 , scalarized residual that accounts for electron energy, ion energy and electron particle flux residuals at all radii simulated) is plotted vs. the iteration number.The 12 converged cases are numbered in green, with vertical dashed lines indicating the transition from one scenario to another, reutilizing flux surrogates.On the right, the transport and target flux profiles (including the three channels) for the case with the lowest residual per scenario simulated are shown.Convergence is achieved for each case when the residual is sufficiently low, and confirmed by visual inspection and expert judgement.Details on the residual definitions and on the convergence metrics are described in Ref. [1].Future work with PORTALS will focus on automation of convergence metrics and more robust convergence criteria.However, considering the stiff and critical gradient behavior of the plasmas studied here, exact flux-matching (e.g. in case #2) is not required for performance predictions, as extensively discussed in Refs.[1,21].

Nonlinear gyrokinetic results
Figure 4 summarizes the predictions of the three kinetic profiles.Electron temperature is found to always be higher than the ion temperature 4 .Gradient profiles show evidence of stiff inner core (r/a < 0.7) and non-stiff edge conditions (r/a = 0.7 − 0.943), which is reminiscent of experimental observations of L-mode plasmas [47].The case with very high density, that resulted in a volume-averaged Greenwald fraction of f G ≈ 0.47 required a negative electron density gradient to achieve the null-flux condition in the plasma core, hence the slightly hollow electron density profile at mid-radius.
The observation of core transport stiffness from Figure 4 is further explored in Figures 5 and 6  plasma becomes stiffer, mostly as a consequence of the increased gyro-Bohm unit, Q GB ∝ T 5/2 e .For similar power (Q in M W ) to be exhausted via turbulent transport, the lower gyro-Bohm normalized energy flux (lower Q/Q GB when Q GB is higher) requires a lower logarithmic gradient, which flattens the gradient profile.This brings up an interesting corollary: under the assumption of local, flux-tube gyrokinetic turbulence, the parameter that may separate L-mode from H-mode is the absolute temperatures and densities, which enter through the gyro-Bohm normalization in transport solvers.This means that, high-temperature L-modes -as those expected in high-confinement, burning-plasma machines-may behave more similarly to present-day Hmodes.However, as it will expanded throughout this section, for the energy confinement time to be at levels predicted by L-mode scaling laws in highperformance plasmas, we require a high edge pressure, which may mean the formation of an edge transport barrier might be required.
Figures 5c) and d), and Figures 5e) and f) show the exploration of the effect of the input power (8.3MW and 18.3M W ) at high (T edge = 2.3keV ) and low (T edge = 0.9keV ) edge temperatures respectively, both at low edge density (n edge = 0.77 • 10 20 m −3 ).At high edge temperatures, the plasma is stiffer, and the edge logarithmic gradient is not very sensitive to the input power.In fact, in these stiff conditions, the higher input power (18.3MW v.s.8.3M W ) resulted in a ∼ 15% lower fusion gain because the fusion power did not increase by the equivalent amount (only ∼ 90%).On the contrary, at low edge temperatures, the increase in input power leads to a ∼ 5% increase in fusion gain.We must emphasize, however, that these results assume that the same edge temperature is achieved at the two different power levels.
Similarly, Figures 6a) and b) demonstrate the almost invariability of the temperature profile when density is varied at high edge temperature (T edge = 2.3keV ).At low edge temperature (T edge = 0.9keV ), however, the temperature profile is more sensitive to density, with the plasma becoming colder as density is increased, as shown in Figures 6c) and d).In both cases the fusion gain is strongly affected by the density, increasing by a factor of ∼ 3 from n edge = 0.77 • 10 20 m −3 to n edge = 2.05 • 10 20 m −3 .We must emphasize here as well that these results assume that the same edge temperature is achieved at the same input power regardless of the edge density, thus assuming a higher edge pressure at higher density for the same input power.Figures 6e) and f) show the density profiles and their gradients predicted in both scans.During these density scans, there is significant variation in the target fluxes per species, with the energy exchange power between electrons and ions and the radiation losses being the most sensitive to the density variation, as shown in Figures 6g)  and c,d) low (0.9keV ) edge temperatures, both at high input power, P in = 18.3MW . Temperature profiles are shown at the top and logarithmic gradients at the bottom.e,f) compilation of density profiles and their gradients for the two scans.g,h) electron and ion sources (input power, alpha power), sinks (radiation losses) and exchange power (electron to ion energy exchange).and h).With the assumption of constant impurity concentration, the radiation losses increase with density, consistent with the model implemented in PORTALS and described in detail in Ref. [21].The trend of the electron to ion energy exchange power is more complex, with a non-monotonic radial variation and including a change of direction (when T i > T e in the outer part of the simulated plasma region) at the highest densities.
Figure 7 summarizes the fusion gain predictions for the 12 scenarios considered, as a function of the three controllable parameters.When all three are plotted independently (Figure 7a), the fusion gain is seen to increase with edge density and edge temperature, as expected, and to either decrease or increase with input power, depending on the edge temperature (as discussed in Figure 5).However, when the edge pressure is used (Figure 7b), the fusion gain is seen to increase linearly with the pressure, organizing all results in roughly the same linear trend.The effect of input power is seen most clearly at the highest edge pressure, a plasma in which operation with 8.3M W of heating power would result in Q ≈ 1.2 (breaking even), while operation with 18.3M W would result in Q ≈ 0.8.The reader is reminded that this work, focused on the core physics, does not make any claim about whether the edge pressure of p edge ≈ 150kP a is achievable with only 8.3M W of input power.Further work is needed to relate these core predictions to the edge conditions.
Figures 7b), c) and d) have colorbars that indicate two metrics of importance to understand the validity of these predictions and the integration of core and edge physics.H-factors (with respect to the τ 89p E [23] and the τ 97L E [48] L-mode confinement scaling laws) in the vicinity of H ≈ 1.2−1.3might be needed to achieve breakeven, which correspond to edge pressures ∼ 35% that of H-mode as predicted by EPED [10].Generally, for the same power and boundary conditions, it is observed that τ 97L E predicts higher (∼ 25%) H-factors than τ 89p E , which means that the later scaling law is more conservative and requires lower edge pressures for nominal confinement (H = 1.0).
We must note that, even if L-mode H factors in excess of unity seem to be required for breakeven in SPARC according to the work presented here, the regression studies that lead to the construction of scaling laws often have large errors (usually σ = ±15%).Furthermore, machines whose data was used to construct the scalings often achieve energy confinement levels higher (and lower) than the predicted scalings, owing to limitations of reproducing complex plasma dynamics with just a few global parameters.Therefore, the results presented here should not be taken as overly pessimistic, as H ≈ 1.2−1.3might be achievable in SPARC, particularly as we learn more about the edge physics and the control of core confinement in this new class of high-performance, burning-plasma machines.
An important factor that affects the predicted performance and that can explain, to some level, the differences between the predictions presented here and those presented in, e.g., Refs.[3,4], is the temperature ratio.Figure 8c shows the temperature ratio as a function of effective collisionality, ν eff = 0.02 keV , indicating that the ratio between volume-averaged temperatures is found to be T i /T e = 0.8 − 1.0 and a strong function of collisionality.Only those cases at high collisionality (those assumed to have low edge temperature and/or high edge density) display a temperature ratio close to unity, while the rest (and most) of cases have T i /T e ≈ 0.85.This has an important effect when doing empirical modeling, as an operational point in POPCONs with the same confinement time and volume-averaged electron temperature will produce less fusion power if the actual ratio is below unity.
Empirical scalings and POPCON modeling of burning plasmas often assume T i /T e = 1.0, an assumption that needs to be revisited.
Interestingly, the steady-state ion to electron heat flux ratio remains high, above Q i /Q e ≳ 1.6 at all locations for all cases studied here, as shown in Figure 8a, regardless of whether there is more electron heating (P OH +P ICRH,e ) than ion heating (P ICRH,i ).In fact, the heat flux ratio seems to be independent of the input power ratio, as shown in Figure 8b, with electron-heating dominated plasmas having as high or even higher ion to electron heat flux ratio.This is consistent with pilot plant case studies [49] that pointed to the dominance of energy transport through the ions even in electron-heating dominated plasmas as a consequence of radiation, energy exchange and the ITG nature of the turbulence.Even though the capability to include the turbulent energy exchange power as part of the energy balance has recently been implemented in PORTALS [1], the simulations presented here were performed prior to this addition and therefore do not account for the turbulent energy exchange power.Motivated by recent results in low-collisionality, ITG-dominated plasmas [50], the turbulent exchange power was evaluated with CGYRO for a characteristic radial location (r/a = 0.75) at low and high collisionalities for the SPARC L-mode database presented here.The results indicate that, while at high collisionality the effect is small (< 6% of the total exchange power), at low collisionality the turbulent exchange power can be as high as 23% of the total exchange power.This result motivates further studies to evaluate the effect of the turbulent exchange power in the energy balance, but it is left for future work.
Density peaking is an important parameter to determine fusion performance, and in SPARC it will be determined by transport (balance of turbulent and neoclassical diffusion and convection), given the lack of significant core fueling.Figure 9 shows the density peaking predicted in this work with PORTALS-CGYRO (which also accounts for neoclassical transport using the NEO code) as a function of the effective collisionality.In this figure, the predictions are compared to the database of source-free plasmas in Alcator C-Mod, ASDEX Upgrade and JET that was used to derive the Angioni scaling of density peaking [51,52].Even though the database belonged to H-modes in the three machines, the density peaking predicted by CGYRO and NEO for this set of SPARC plasmas is found to be in remarkable agreement over a wide range of collisionalities (note that the x-axis is in logarithmic scale).

On the effect of impurity content
Previous results correspond to the PRD nominal choice of f DT = 0.85 and Z eff = 1.5, which is achieved in the simulations by assuming a "lumped" low-Z impurity that satisfies the quasineutrality and Z eff constraints (accounting for the 5% ICRH minority concentration).Here we explore the effect of an increased effective ion charge on the predictions, by either controllable (divertor impurity seeding) or uncontrollable sources.For this study, we select a case at high input power (P in = 18.3MW ) and in middle of the T − n database (case circled in Figure 4b, with T edge = 1.6keV and n edge = 1.2 • 10 20 m −3 ).We Fig. 9 Electron density peaking as a function of the effective collisionality for the 12 cases studied here (red stars) and for a database of source-free Alcator C-Mod, ASDEX Upgrade and JET plasmas [51,52] (purple).Cases from the database with core fueling from neutral beams have been shaded.Error bar on predicted peaking accounts for the edge density profile.The lower bound is obtained with a density profile flat in the r/a = 0.943 − 1.0 region, and the upper bound is obtained with a density profile equal to zero (thus with a discontinuity at r/a = 0.943) in that region.
study two additional impurity levels: Z eff = 2.0 (f DT = 0.79) and Z eff = 2.5 (f DT = 0.73), which complement the nominal assumption of Z eff = 1.5 (f DT = 0.85).These higher impurity content plasmas are achieved by increasing the concentration of the "lumped" impurity until the desired Z eff is reached.Impurity charge (Z = 9) and tungsten concentration (f W = 1.5•10 −5 ) are assumed fixed.
Nonlinear gyrokinetic profile predictions, at constant edge conditions, display an increase in performance as the impurity content increases, as shown in Figure 10.Not only higher ion temperatures (Figure 10b) are achieved via turbulence stabilization -i.e. higher steady-state normalized logarithmic ion temperature gradients (Figure 10e), for the same input power-despite higher radiation losses (Figure 10f), but also higher fusion power (Figure 10c) despite the lower fuel density (Figure 10a).The nominal Z eff = 1.5 case had Q ∼ 0.29, and increasing the impurity content led to Q ∼ 0.33 and Q ∼ 0.4 for Z eff = 2.0 and Z eff = 2.5 respectively, i.e. a 38% increase in fusion gain from the nominal to the Z eff = 2.5 case.This increase in performance via impurity-induced suppression of core turbulence is consistent with past experimental and modeling studies [53][54][55][56][57].To understand the turbulence stabilization effect, Figure 10d) shows a comparison between the steady-state ion energy flux for the nominal Z eff = 1.5 case and the fluxes (coming from turbulent CGYRO and neoclassical NEO transport) that result in increasing the impurity content to Z eff = 2.0 and Z eff = 2.5.As the target flux is very similar, the gradients need to peak more to achieve similar levels as in the Z eff = 1.5 case (blue), resulting in the increased performance in flux-matched conditions.For clarification, fluxes of the Z eff = 2.0 and Z eff = 2.5 cases in Figure 10d) are not in steady-state, but represent the iteration #0 of the PORTALS-CGYRO process, which corresponds to the gradients of the nominal (converged) Z eff = 1.5 case but with the differing impurity content.
These results are very promising and suggest that a path for higher performance may exist that also satisfies at the same time core-edge integration constraints (divertor heat flux limits and associated impurity seeding needs).The further exploration of such operational scenarios and the investigation of whether the lower turbulent transport levels are due solely to the effect of the effective ion charge on the ion temperature critical gradient or on dilution effects (which can lead to interesting optimization studies of the puffed impurity of choice) will be subject of future work.

Quasilinear gyro-fluid results
The TGLF quasilinear model for core turbulence has been tested against the database of 528 nonlinear gyrokinetic results that were presented in this work.With a much reduced computational cost, TGLF estimates transport fluxes by solving a linear system of gyro-fluid equations [35] as an eigenvalue problem and uses a saturation rule [36] to calculate fluxes from the growth rates and real frequencies of unstable micro-instabilities, together with the eigenmode solution and associated quasilinear weights.The saturation rule in TGLF has evolved over the years, as the community has gathered more understanding of the dynamics of turbulence in the plasma core.Here we test the three latest saturation rules, SAT1 [58,59], SAT2 [60] and SAT3 [61], with electromagnetic effects (δϕ, δA ∥ , and δB ∥ ) and only capturing long-wavelength turbulence (k θ ρ s ≤ 1.22) for consistency with the CGYRO simulations.As a reminder, the SAT1, SAT2 and SAT3 saturation rules were originally constructed to reproduce the cross-scale turbulence coupling, the 3D dependence (poloidal angle, radial and poloidal wavenumber) of the saturated intensity and the isotope effect, respectively, with their regression coefficients fitted with databases of a few tens of nonlinear CGYRO simulations.
Figure 11 shows a comparison of CGYRO and TGLF predictions for each of the three transport fluxes under consideration (electron energy, ion energy and electron particle fluxes) at the 6 flux surfaces considered in this work, totalling 528 local simulations.Remarkably, TGLF is capable of reproducing the transport fluxes well, spanning over four orders of magnitude.This is in spite of the fact that TGLF is fitted to only a few tens of CGYRO simulations per saturation rule, in different parameter spaces than this database of SPARC plasmas, and with just ∼10 CPU-seconds per evaluation (vs ∼200 GPU-hours for nonlinear CGYRO).Generally, it is observed that SAT3 performs the best among the saturation rules tested in this work.As seen in Figure 11d, SAT3 performs on average better than SAT1 and SAT2.The strongest differences between CGYRO and TGLF SAT3 are observed at the edge, where TGLF is generally seen to overpredict electron and ion energy fluxes.Unfortunately, SAT3 performs the worst of the three at reproducing particle fluxes.Interestingly, SAT2 is found to perform worse (highest root mean square error) than SAT1 and SAT3 at reproducing energy fluxes at the plasma edge (positions r/a = 0.9 and 0.943).This is somewhat in contrast to the results presented in Refs.[41,42], where SAT2 was found to be in good agreement with ASDEX Upgrade experiments and capable of reproducing the parametric dependencies of the L-mode scaling law well.
A comparison of gradient-driven transport fluxes such as the one just described in Figure 11 is not entirely fair, as the critical gradient and stiff nature of turbulence in the core of tokamak plasmas may result in very different predicted fluxes for very similar gradients.Therefore, a flux-matched comparison of predicted kinetic profiles can provide more important insights.PORTALS-TGLF (SAT3) was run for all cases in the database, and the results of comparing the predicted profiles and performance with PORTALS-CGYRO are shown in Figure 12.
Generally, TGLF SAT3 is seen to reproduce the predicted temperature profiles well (Figures 12a and b), with the exception of the very edge, r/a = 0.943, where CGYRO predicted higher gradients for both electron and ion temperatures.This is consistent with the results in Figure 11 that showed TGLF SAT3 overpredicting the transport fluxes at the edge, which leads to lower gradients in steady-state conditions when the plasma has a non-stiff edge (low density and temperature boundary conditions).In global terms, ion temperature peaking is seen to be underpredicted by TGLF SAT3, while electron temperature peaking is overpredicted (Figure 12g).The volume-averaged temperatures are seen to be slightly underpredicted by TGLF SAT3 (Figure 12h), which is consistent with the overprediction of the transport fluxes at the edge.As per the electron density gradients (Figure 12c), the agreement is not as good as for the temperature gradients, with more scattered behavior and the appearance of hollowed profiles with TGLF SAT3, particularly at r/a = 0.9 and 0.87.This leads generally to an underprediction of the density peaking with TGLF (Figure 12f), consistent with other studies with quasilinear profile predictions [21].Example density profiles showing disagreement (lower peaking with the quasilinear model) between CGYRO and TGLF are depicted in Figure 13.
In terms of performance, energy confinement time is underpredicted by TGLF, a consequence of the underprediction of both the volume averaged temperatures and density peaking.The lower temperatures and densities lead to lower fusion power, for which CGYRO predicts generally higher performance.These results are, however, encouraging, as TGLF SAT3 is capable of predicting performance reasonably well at a tiny fraction of the computational expense, which will allow the further exploration of the parameter space and optimization studies at much reduced cost.
TGLF SAT3 is capable of predicting the increase of performance with the increase in impurity content, as shown in Figure 13.In CGYRO, it was observed that Q increased by ∼ 38% from the nominal Z eff = 1.5 case to the Z eff = 2.5 case, while in TGLF Q increased by ∼ 30%.However, the physics reasons behind the improvement in performance are different.In TGLF, the change in edge temperature gradients is not as strong as in CGYRO, and the improvement in performance is mostly due to the increase in near-axis (r/a = 0.35) temperature gradients.An additional case with Z eff = 3.0 was also studied with TGLF, and it is seen to predict a roll-over in performance, which is consistent with the dilution of fusion fuel becoming too strong to be compensated by the turbulence stabilization, as f DT = 0.67 for Z eff = 3.0.
It is important to note that predictions with TGLF spanning to short wavelengths (TGLF default k θ ρ s ≤ 24.0 instead of k θ ρ s ≤ 1.22 used here for consistency with CGYRO) did show differences in standalone fluxes and in kinetic profile predictions, particularly with electron energy and particle fluxes having non-negligible contributions at k θ ρ s ∼ 2 − 3 at the larger radii, even in situations with max (γ/k θ ρ s ) high−k < max (γ/k θ ρ s ) low−k , a zonal flow mixing criterion usually employed to assess the relative influence of multi-scale effects [62].While this could suggest the need to run nonlinear gyrokinetic simulations that also include shorter wavelengths than used in this study (to capture mid and high wavenumber TEM and ETG modes, and cross-scale interactions), further work is needed to verify the capabilities of TGLF to capture cross-scale interactions and mid-wavenumber flux contributions correctly, for which very limited studies exist.Due to the very large computational resources required for such an exercise, this is out of the scope of this work but will be the subject of future investigations.

Discussion
Sections 4 and 5 have presented and discussed the results of the ion-scale nonlinear gyrokinetic and quasilinear gyrofluid modeling of SPARC core plasmas when four parameters are varied: the edge (r/a = 0.943) density, the edge temperature, the input power and the impurity content.While there are still unknowns and the assumed edge conditions are not fully understood, the results presented here provide a few take-away lessons that will be useful for the design of experiments and operation of SPARC and other future burning-plasma machines.
It is found that L-modes characterized by energy confinement times in line with empirical scaling laws (i.e.H ∼ 1.0; where H is the confinement factor relative to empirical L-mode confinement time) in high-performance plasmas require moderate (∼35% that of H-mode) pressures.In these hot, dense conditions, turbulence is very stiff, as a consequence of the high gyro-Bohm unit and dominance of ITG turbulence, generally expected in reactorlike plasmas [49].As such, these plasmas behave more closely to present-day H-modes than to L-modes, which could also explain the very good agreement with the database of H-mode density peakings.In conditions closer to H ∼ 1.0, the normalized logarithmic temperature gradient profiles become radially flat and stiffer, requiring a high edge pressure to maintain the energy confinement time as predicted by the scalings.In general, this opens the question of the validity of L-mode confinement scaling laws in high-temperature, high-density plasmas (i.e.reactor-like and burning plasmas), where the core behaves more similarly to H-modes.This apparent contradiction may be clarified by looking at core gyrokinetic turbulent transport as a dimensionless, local process, where normalized quantities such as logarithmic gradients normalized to the plasma minor radius, collisionality or temperature ratios are the ones that determine the gyro-Bohm normalized turbulent transport levels.In this sense, as long as the assumptions of flux-tube gyrokinetic theory (well attained in burning plasmas) are valid, the core transport dynamics are independent of the edge conditions.If L-modes in high-performing devices like SPARC are empirically found to achieve high core temperatures and densities, the behavior of the core may be more similar to H-modes, where the high temperature is achieved by means of an edge transport barrier.The lack of neutral beam injectors and limited flexibility to control the current density profile in SPARC may mean that the formation of internal transport barriers should not be, generally, expected.It is also important to note that the increased stiffness found in these high-performance "L-mode" plasmas, which potentially separates them from present-day machines, also implies that the edge transport similarity assumed when using the Alcator C-Mod database of edge gradient profiles needs to be revisited.Future work will include the characterization of the edge transport in SPARC L-modes via multi-machine databases that span a wide range of edge collisionalities and input powers.
The requirement of high edge pressure may be a challenge for the integration of core performance and edge and divertor constraints, particularly if a transport-regulated edge barrier, r/a = 0.943 − 1.0, is needed.As a result, the sustainment of a high enough core energy confinement time in "L-mode" conditions may still require high input power as a result of the high magnetic field that contributes to a low E × B shearing rate [41].This high power to be transported through r/a = 0.943 − 1.0 to sustain the high gradient would be challenging for managing a reasonable divertor heat flux, particularly during early operations, when the divertor conditions and flux mitigation will still be somewhat unknown and solutions to the heat flux and loading challenges will still be under development.Given the present uncertainties and risks that we will only be able to retire after first experiments have been performed, the SPARC team is devising contingency plans that include potentially pivoting to H-mode operation (including at reduced magnetic field, 8T [63]) if needed in order to get the Q > 1 milestone.
The attainment of a high edge pressure may be achieved by means of access to I-mode regimes [64], which have potential advantages to H-modes, such as ELM avoidance, low impurity accumulation and better density control.I-mode edge profiles were almost certainly the main contributors to the edge profiles for the high edge temperatures (T edge > 0.6keV ) shown in the database of Alcator C-Mod plasmas in Figure 2. The rather large ion to electron heat flux ratio seen in Figure 8 (even at low Greenwald fractions, below what could be estimated to be the transition to the low-density branch for mode transitions [65]) may also be a sign of the potential for I-mode-like regimes, as there is evidence for a critical ion heat flux required for L to I transitions [66].The potential for I-mode access is being examined and will be the subject of a follow-on publication.The high Q i /Q e ratio also suggests that T i > T e may be observed at the separatrix, which could further lead to higher fusion power at the same edge pressure.However, its impact at the boundary condition for the core simulations (r/a = 0.943) has not been explored yet and it will be part of future work.The modeling of transport from separatrix to the core, coupling fueling, edge E × B shear and integration with scrape-off layer and divertor physics models (a truly self-consistent integrated model) will also be part of future work.
Interestingly, this work has found that increasing the impurity content in these plasmas, even if the fusion fuel becomes more diluted, may lead to higher fusion power and performance, under the assumption of a fixed edge boundary condition.This is consistent with previous gyrokinetic modeling and experimental validation work [55].This work so far has only explored the effect of increasing the effective ion charge with a fixed low-Z, "lumped" impurity, but future work will explore the effect of other impurity choices and their concentration.This is very encouraging, as lower edge pressures may be allowable if higher fusion power is attained by means of core turbulence stabilization via higher impurity content.Other potential phenomena, such as the effect of fast ions on turbulence [67], have not been addressed in this work, which could also favor high-performance cores with lower edge pressures.Toroidal rotation has also been assumed to be zero throughout the core, as well as the particle fueling.Both effects are generally expected to increase performance (the effect of rotation on standalone turbulence fluxes was explored in Ref. [68] for the SPARC PRD scenario), but how much of an increase still remains an open question, until reliable first-principles or physics-based models of edge particle Performance in SPARC first campaign fueling and residual momentum stresses are available in a fully predict-first, physics-based fashion.

Conclusions
This work has presented a high-fidelity core (r/a < 0.943) transport study of, arguably, unprecedented fidelity in magnetically confined plasmas.Thanks to surrogate-acceleration of transport solvers enabled by PORTALS, nonlinear gyrokinetics can now be used to predict operational scenarios in future burning plasmas and in this work a large database of scenarios was predicted at a moderate computational expense.
Results of the parametric exploration of SPARC early-campaign scenarios project near-breakeven conditions in L-mode-like plasmas.Breakeven is predicted when the edge pressure becomes ∼ 35% that of H-mode, but the actual fusion gain will be determined by how much input power is required to sustain such edge pressures.Higher impurity content is, somewhat unexpectedly, predicted to increase fusion power, opening an interesting path for optimization and transport physics investigations.Enabled by an extensive database of 528 local turbulence simulations, quasilinear results with TGLF show promise for the use of reduced models to widely explore the parameter space because, although some differences appear (particularly in the density profile predictions and edge temperature gradients), TGLF is able to capture the parametric dependencies well, together with the overall plasma performance, at a fraction of the cost.
This work has pointed out the need to understand better the edge conditions in L-mode plasmas, as well as the validity of energy confinement scaling laws to predict the performance of high-temperature and high-density L-modes.The first plasmas in upcoming machines such as SPARC and ITER will be L-mode-like.Furthermore, the challenge of core-edge integration, ELM mitigation and divertor heat flux control has recently led fusion researchers to explore the potential of non-H-mode operation directly in fusion power plant studies (e.g.[69][70][71]).It is only by the careful validation [72] of predictive simulations against experimental data that we will be able to understand the physics of these plasmas and to optimize the performance of such devices [73].More importantly, L-mode operation in high-performance devices could be more similar to H-mode operation than to present-day L-modes, consequence of the high gyro-Bohm transport in high temperature, high density core plasmas.This could have important implications for the design and operation of future fusion power plants.The validation of this modeling insight in the first campaigns of SPARC will certainly be a key milestone in the path to economic fusion energy.
The upcoming years, as we witness SPARC's first experimental operations, promise not only a test-bed for our predictions but also a fertile ground for learning and improving our predictive capabilities.This cycle of prediction, validation, and refinement is crucial as we navigate from achieving burning plasma conditions to envisioning fusion power plants.The insights from validating high-fidelity simulations such as the ones presented in this work against real-world data will be invaluable, guiding the optimization of future fusion devices.As we refine these models, their improved predictive power will be instrumental in the design, operation, and optimization of fusion power plants, marking a significant leap towards a future powered by clean, sustainable fusion energy.

Fig. 1
Fig. 1 TRANSP interpretive analysis of ICRH power delivery to plasma.Left column depicts the change in the fraction of ICRH power delivered to ions as a function of density, temperature and temperature ratio.Right column shows the change in the ratio of power delivered within the r/a = [0.15,0.25, 0.35, 0.45, 0.55, 0.75] flux surfaces as a function of the same parameters.Shaded regions indicate the range of variation of the parameters in the predictive results from this work.When not varied, ⟨ne⟩ = 1.5 • 10 20 m −3 , ⟨Te⟩ = 4.0keV and T i /Te = 1.0 are assumed.

Fig. 2 (
Fig. 2 (top) Distribution (mean and standard deviation) of normalized logarithmic gradient in database of Alcator C-Mod non-H-mode plasmas in the edge region, r/a = 0.94 − 1.0.Distributions are separated by the value of edge temperature.(bottom) Distribution of edge temperatures when the gradients are radially integrated from a boundary condition of 200eV at r/a = 1.0 (dashed is the integration of the mean, solid is the average of the integration of upper and lower bounds.)

Fig. 3
Fig. 3 Summary of PORTALS-CGYRO multi-channel convergence.(left) L 1 residual normalized to the radial-average of the target flux per iteration, with the 12 converged cases numbered in green.(right) Transport and target flux profiles for the 12 converged cases, including the three channels: ion energy flux, electron energy flux and electron particle flux (transformed into energy units as a convective flux, Q Γ = 3 2 TeΓ).Shaded area indicates ±2σ as estimated for the CGYRO transport fluxes.
Figure4summarizes the predictions of the three kinetic profiles.Electron temperature is found to always be higher than the ion temperature4 .Gradient profiles show evidence of stiff inner core (r/a < 0.7) and non-stiff edge conditions (r/a = 0.7 − 0.943), which is reminiscent of experimental observations of L-mode plasmas[47].The case with very high density, that resulted in a volume-averaged Greenwald fraction of f G ≈ 0.47 required a negative electron density gradient to achieve the null-flux condition in the plasma core, hence the slightly hollow electron density profile at mid-radius.The observation of core transport stiffness from Figure4is further explored in Figures5 and 6.Figures 5a) and b) show the exploration of effect of the edge temperature at low input power (P in = 8.3M W ) and low edge density (n edge = 0.77 • 10 20 m −3 ), indicating a clear flattening of the logarithmic gradient as the boundary condition is increased.At high edge temperatures, the

Fig. 4
Fig. 4 Summary of database of PORTALS-CGYRO profile predictions.a) and b) depict the 12 scenarios chosen from combinations of the three input parameters (input power, edge temperature and edge density).c), e) and g) show the electron temperature, ion temperature and electron density prediction from the r/a = 0.943 boundary condition.c) and e) also show the deposited power density for electrons and ions respectively (only two cases studied, P in = 8.3M W and P in = 18.3MW , with fixed shape and species partition).d), f) and h) show normalized logarithmic gradients of the profiles.In all figures, red and blue are used to indicate low input power (P in = 8.3M W ) and high input power (P in = 18.3MW ) respectively, and different line styles for the temperature boundary conditions, for better visualization.In b), point with P in = 18.3MW , n edge = 1.2 • 10 20 m −3 and T edge = 1.6keV is highlighted with a black circle, as it is the scenario chosen for studies of the impurity effect in Section 4.1.

Fig. 5
Fig. 5 Exploration of core transport dynamics (predicted electron and ion temperatures together with power deposition profile at the top and normalized logarithmic gradients at the bottom) at fixed edge density, n edge = 0.77 • 10 20 m −3 .a,b) Scan of edge temperature at low input power (P in = 8.3M W ), c,d) scan of input power at high edge temperature (T edge = 2.3keV ), and e,f) scan of input power at low edge temperature (T edge = 0.9keV ).

Fig. 6
Fig. 6 Scan of edge density at a,b) high (2.3keV ) and c,d) low (0.9keV ) edge temperatures, both at high input power, P in = 18.3MW . Temperature profiles are shown at the top and logarithmic gradients at the bottom.e,f) compilation of density profiles and their gradients for the two scans.g,h) electron and ion sources (input power, alpha power), sinks (radiation losses) and exchange power (electron to ion energy exchange).

Fig. 7
Fig. 7 Study of fusion gain for the 12 scenarios considered.a) Fusion gain as a function of edge density, edge temperature and input power as separate free parameters.b) Fusion gain as a function of edge pressure and input power, with colorbar indicating the level of energy confinement over the τ 89p E confinement law (H factor).c) H factor (both for τ 89p E and τ 97L E scaling laws) as a function of edge pressure and input power, with colorbar indicating fusion gain.d) Fusion gain as a function of edge pressure and input power, with colorbar indicating the ratio of edge pressure over the EPED-predicted H-mode pressure at each density.

Fig. 8
Fig. 8 a) Radial profile of the ratio between ion and electron fluxes (transport flux in solid and input power flux in dashed).b) Transport flux ratio vs input power flux ratio at midradius (r/a = 0.5).c) Ion to electron temperature ratio on-axis and volume averages, as a function of the effective collisionality.

Fig. 10
Fig. 10 Study of the effect of impurity content on performance.a) Electron and fuel density profiles, b) electron and ion temperature profiles, c) fusion power density profile, d) ion energy flux comparison between converged nominal case and transient (iteration #0) higher Z eff cases, e) normalized logarithmic ion temperature gradient, and f) radiated power profile.

Fig. 11
Fig. 11 Validation of TGLF saturation rules against nonlinear CGYRO for the database of 14 scenarios (12 nominal scenarios + 2 variations of Z eff ), resulting in 528 local nonlinear CGYRO simulations to be compared to TGLF.a) Electron energy, b) ion energy and c) electron particle fluxes at the 6 flux surfaces for TGLF SAT3 and nonlinear CGYRO.d) Relative change in root mean square error of the three transport fluxes at each location when changing the saturation rule from SAT3 to SAT1 and SAT2.Horizontal dashed lines indicate the average (radius and channel) relative change in root mean square error.For those points above zero, SAT3 performs better.

Fig. 12
Fig. 12 Validation of TGLF SAT3 against nonlinear CGYRO in multi-channel, flux-matched conditions for both models.Plots provide comparisons between CGYRO (x-axis) and TGLF SAT3 (y-axis) of the following quantities: a), b) and c) local normalized logarithmic gradients of the electron temperature, ion temperature and electron density in steady state conditions, d) fusion gain, e) energy confinement time, f) density peaking (both standard definition, n e,0.2 /⟨ne⟩, and n e,0 /n e,edge ), g) electron and ion temperature peaking, and h) volumeaveraged electron and ion temperatures.

Fig. 13
Fig.13Study of the validity of TGLF to predict the effect of impurity content on performance: profiles of electron temperature, ion temperature and electron density profiles and their respective logarithmic gradients.Fusion gain vs Z eff is also shown.Z eff = 3.0 case is only available with TGLF as of the time of writing.