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 reveals 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.

## I. INTRODUCTION

Upcoming burning-plasma experiments in both ITER^{1} and SPARC^{2,3} 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.^{4} Even though there are variations to H-mode operation that are ELM-free,^{5} 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^{6} (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 the subject of ongoing research (e.g., Refs. 7 and 8). Most predictive work—either with the use of empirical confinement laws such as $ \tau E 98 y 2$ (Ref. 1) or edge pedestal models such as EPED^{9}—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.^{10,11}

While both H-mode and I-mode operation are 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,^{12–15} 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,^{16} 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.

## II. BACKGROUND

SPARC^{2,3} is a high-field ( $ B = 12.2 T$), medium-size ( $ R = 1.85 \u2009 m$), 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 \u2264 8.7 \u2009 MA$) and high absolute density ( $ \u27e8 n e \u27e9 \u2248 3 \xd7 10 20 \u2009 m \u2212 3$ in the baseline H-mode regime), while remaining reasonably far from stability boundaries ( $ q 95 \u2248 3.4 , \u2009 f G = n e / n G \u2248 0.37 , \u2009 \beta N \u223c 1.0$).^{2,3,11,17} These conditions are favorable 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 \u2212 11$ (Refs. 2, 11, and 21) operating with a mix of 50–50 DT fuel, provided that the pedestal remains at EPED-predicted^{10} 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 full-power 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 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 analysis^{22,23} 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.^{24} To reduce the time required for scenario development, a family of candidate scenarios was developed for parameters that are relaxed slightly relative to the PRD. These scenarios are for single-null geometries at slightly reduced elongation (areal elongation $ \kappa A = 1.7$ vs $ \kappa A = 1.75$ in the PRD) and plasma current ( $ I p = 8.5 \u2009 MA$ vs $ I p = 8.7 \u2009 MA$ in the PRD), for densities ranging from $ \u27e8 n e \u27e9 \u2248 1.0 \xd7 10 20 m \u2212 3$ (near the predicted LOC-SOC transition^{25}) to $ \u27e8 n e \u27e9 \u2248 2.1 \xd7 10 20 m \u2212 3$.^{12} 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.^{13} 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 modeling 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.

### A. 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^{16} is a novel, open-source^{26} toolset developed to accelerate the numerical convergence of time-independent transport solvers. The speedup achieved by the utilization of surrogate-based optimization techniques^{27} in PORTALS has enabled some of the highest fidelity predictions of core transport to date. Direct nonlinear gyrokinetic simulations with the CGYRO^{28} code have been used to predict steady-state (multi-channel flux-matched) plasmas in SPARC PRD,^{21} ITER IBS,^{29} DIII-D ITER Similar Shape,^{16,30} ARC,^{31} and JET H, D and T Ohmic plasmas.^{32,33}

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^{28} code. CGYRO is an initial-value, 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^{34} and is optimized for scalability and performance on modern high-performance computing architectures,^{35} 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^{36} and uses a saturation rule^{37} 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,^{38} but found to be negligible compared to turbulence as modeled by both CGYRO and TGLF ( $ < 1.2 %$ in all flux-matched radii).

## III. 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.7 \u2009 MA$), magnetic field ( $ B T = 12.2 T$), and shaping ( $ \kappa sep = 1.97 , \u2009 \delta 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 \xd7 10 \u2212 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. 2, 3, and 11). Auxiliary input power absorbed from ICRH is considered to be below $ P ICRH < 15 \u2009 MW$, 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, $ \Gamma e = 0$.

### A. On parameter selection and assumptions

In the parameter scans presented in the remainder of the paper, 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 $ \rho tor = 0.9$ (square root of normalized toroidal flux) flux surface, which corresponds to $ r / a \u2248 0.943$ (normalized half-width of midplane intercept) and $ \psi n \u2248 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 Sec. III B) 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 \u2212 2.3 \u2009 keV , \u2009 n edge = 0.77 \u2212 2.05 \xd7 10 20 \u2009 m \u2212 3$, and $ P in = 8.3 \u2212 18.3 \u2009 MW$. There will be an exploratory, very high-density case ( $ n edge = 4.0 \xd7 10 20 \u2009 m \u2212 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 Figs. 4(a) and 4(b) (to be discussed in Sec. IV), 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, $ \beta N \u2272 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^{39} with various potential kinetic profiles were performed. TRANSP is run with the TORIC full-wave code^{40} coupled with the bounce-averaged Fokker–Planck FPPMOD model^{41} 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,^{11} 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 \u223c 0.65 \u2212 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. 16). While the change in power deposition at that location is not completely negligible ( $ \u223c 25 %$ increase in power when doubling the density from $ \u27e8 n e \u27e9 = 1.0 \xd7 10 20 \u2009 m \u2212 3$ to $ 2.0 \xd7 10 20 \u2009 m \u2212 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 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*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.

_{i}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.3 \u2009 MW$), 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.3 \u2009 MW$ (5 MW ICRH + $ 3.3 \u2009 MW$ of Ohmic) and $ 18.3 \u2009 MW$ (15 MW ICRH + $ 3.3 \u2009 MW$ 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.

### B. 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 \u2013 0.26$. This work assumes that fueling systems will be capable of sustaining such edge density. We must note that the edge-separatrix region, $ r / a = 0.943 \u2013 1.0$, has no effect on core performance from the perspective of flux-tube local gyrokinetic turbulence and the assumptions employed in this work.

The determination of edge temperature and its validity is more subtle. Transport predictions from the separatrix in L-modes, while possible,^{42,43} 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 time-steps 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 \xb1 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^{44}) and associated $ a / L T e$ 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 T e$ in the edge region, $ r / a = 0.94 \u2212 1.0$, is shown in the top panel of Fig. 2.

The bottom panel of Fig. 2 shows the distribution of edge temperatures when gradients are radially integrated from an assumed separatrix temperature of 200 eV 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 \xb1 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 temperature gradient mode limited edge.^{75} Furthermore, 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 gradient 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. Second, 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 200 eV (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 Sec. VI).

### C. 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. 16, 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. 16 that suggested that five 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-mode-like 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 six 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^{28} 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 (40 GB) GPUs.

All gyrokinetic simulations described in this paper were performed with the CGYRO code.^{28}^{,} 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 micro-tearing (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.^{45} All simulations utilized relatively high-fidelity physics including geometry represented by the Miller Extended Harmonic (MXH) formulation,^{46} fully electromagnetic turbulence ( $ \delta \varphi , \u2009 \delta A \u2225$, and $ \delta B \u2225$), Sugama collisions,^{47} and four 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 ] \u223c 125 \xd7 120 \rho s$ utilizing 24 complex toroidal modes ( $ N \alpha $) and 512 radial modes (*N _{r}*) across the domain to represent turbulence from $ k \theta \rho s = 0.0$ to 1.22. Simulations also used 8 energy grid points (

*N*), 24 pitch angles ( $ N \xi $), and 24 points in theta ( $ N \theta $). 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.

_{u}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*, $ \Gamma e$) convergence, but the steady-state predictions (

_{i}*T*,

_{e}*T*,

_{i}*n*) 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.

_{e}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 judgment. Details on the residual definitions and on the convergence metrics are described in Ref. 16. 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. 16 and 21.

## IV. 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. In this work, the only assumption that could somewhat change the temperature ratio is the assumed $ T i / T e$ at $ r / a = 0.943$, which has only been given the value of 1.0. Future work could explore sensitivity to this assumption, but we do not expect a strong deviation from 1.0 given the high edge collisionality. Gradient profiles show evidence of stiff inner core ( $ r / a < 0.7$) and non-stiff edge conditions ( $ r / a = 0.7 \u2212 0.943$), which is reminiscent of experimental observations of L-mode plasmas.^{48} The case with very high density, that resulted in a volume-averaged Greenwald fraction of $ f G \u2248 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 Fig. 4 is further explored in Figs. 5 and 6. Figures 5(a) and 5(b) show the exploration of effect of the edge temperature at low input power ( $ P in = 8.3 \u2009 MW$) and low edge density ( $ n edge = 0.77 \xd7 10 20 \u2009 m \u2212 3$), indicating a clear flattening of the logarithmic gradient as the boundary condition is increased. At high edge temperatures, the plasma becomes stiffer, mostly as a consequence of the increased gyro-Bohm unit, $ Q G B \u221d T e 5 / 2$. For similar power (*Q* in *MW*) to be exhausted via turbulent transport, the lower gyro-Bohm normalized energy flux (lower $ Q / Q G B$ 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 H-modes. However, as it will expanded throughout this section, for the energy confinement time to be at levels predicted by L-mode scaling laws in high-performance plasmas, we require a high edge pressure, which may mean the formation of an edge transport barrier might be required.

Figures 5(c)–5(f) show the exploration of the effect of the input power ( $ 8.3$ and $ 18.3 \u2009 MW$) at high ( $ T edge = 2.3 \u2009 keV$) and low ( $ T edge = 0.9 \u2009 keV$) edge temperatures, respectively, both at low edge density ( $ n edge = 0.77 \xd7 10 20 \u2009 m \u2212 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.3$ vs $ 8.3 \u2009 MW$) resulted in a $ \u223c 15 %$ lower fusion gain because the fusion power did not increase by the equivalent amount (only $ \u223c 90 %$). On the contrary, at low edge temperatures, the increase in input power leads to a $ \u223c 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, Figs. 6(a) and 6(b) demonstrate the almost invariability of the temperature profile when density is varied at high edge temperature ( $ T edge = 2.3 \u2009 keV$). At low edge temperature ( $ T edge = 0.9 \u2009 keV$), however, the temperature profile is more sensitive to density, with the plasma becoming colder as density is increased, as shown in Figs. 6(c) and 6(d). In both cases, the fusion gain is strongly affected by the density, increasing by a factor of $ \u223c 3$ from $ n edge = 0.77 \xd7 10 20 \u2009 m \u2212 3$ to $ n edge = 2.05 \xd7 10 20 \u2009 m \u2212 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 6(e) and 6(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 Figs. 6(g) and 6(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*in the outer part of the simulated plasma region) at the highest densities.

_{e}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 [Fig. 7(a)], 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 Fig. 5). However, when the edge pressure is used [Fig. 7(b)], 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.3 \u2009 MW$ of heating power would result in $ Q \u2248 1.2$ (*breaking even*), while operation with $ 18.3 \u2009 MW$ would result in $ Q \u2248 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 \u2248 150 \u2009 kPa$ is achievable with only $ 8.3 \u2009 MW$ of input power. Further work is needed to relate these core predictions to the edge conditions.

Figures 7(b)–7(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 $ \tau E 89 p$ (Ref. 24) and the $ \tau E 97 L$ (Ref. 49) L-mode confinement scaling laws) in the vicinity of $ H \u2248 1.2 \u2212 1.3$ might be needed to achieve breakeven, which correspond to edge pressures $ \u223c 35 %$ that of H-mode as predicted by EPED.^{9} Generally, for the same power and boundary conditions, it is observed that $ \tau E 97 L$ predicts higher ( $ \u223c 25 %$) H-factors than $ \tau E 89 p$, 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 $ \sigma = \xb1 15 %$). Furthermore, machines whose data were 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 \u2248 1.2 \u2212 1.3$ might 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. 2 and 3 is the temperature ratio. Figure 8(c) shows the temperature ratio as a function of effective collisionality, $ \nu eff = 0.02 \xd7 \u27e8 n e \u27e9 20 \xb7 R 0 , m \xb7 \u27e8 T e \u27e9 keV \u2212 2$, indicating that the ratio between volume-averaged temperatures is found to be $ T i / T e = 0.8 \u2212 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 \u2248 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 \u2273 1.6$ at all locations for all cases studied here, as shown in Fig. 8(a), 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 Fig. 8(b), 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^{50} 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,^{16} 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,^{51} 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.^{52,53} 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).

### A. 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.3 \u2009 MW$) and in middle of the

*T*–

*n*database [case circled in Fig. 4(b), with $ T edge = 1.6 \u2009 keV$ and $ n edge = 1.2 \xd7 10 20 m \u2212 3$]. We study two additional impurity levels: $ Z eff = 2.0$ (

*f*= 0.79) and $ Z eff = 2.5$ (

_{DT}*f*= 0.73), which complement the nominal assumption of $ Z eff = 1.5$ (

_{DT}*f*= 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 (

_{DT}*Z*= 9) and tungsten concentration ( $ f W = 1.5 \xd7 10 \u2212 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 Fig. 10. Not only higher ion temperatures [Fig. 10(b)] are achieved via turbulence stabilization—i.e., higher steady-state normalized logarithmic ion temperature gradients [Fig. 10(e)], for the same input power—despite higher radiation losses [Fig. 10(f)], but also higher fusion power [Fig. 10(c)] despite the lower fuel density [Fig. 10(a)]. The nominal $ Z eff = 1.5$ case had $ Q \u223c 0.29$, and increasing the impurity content led to $ Q \u223c 0.33$ and $ Q \u223c 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.^{54–58}

To understand the turbulence stabilization effect, Fig. 10(d) 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 Fig. 10(d) 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.

## V. 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^{36} as an eigenvalue problem and uses a saturation rule^{37} 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,^{59,60} SAT2,^{61} and SAT3,^{62} with electromagnetic effects ( $ \delta \varphi , \u2009 \delta A \u2225$, and $ \delta B \u2225$) and only capturing long-wavelength turbulence ( $ k \theta \rho s \u2264 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, totaling 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 Fig. 11(d), 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. 42 and 43 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 Fig. 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 Fig. 12.

Generally, TGLF SAT3 is seen to reproduce the predicted temperature profiles well [Figs. 12(a) and 12(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 Fig. 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 [Fig. 12(g)]. The volume-averaged temperatures are seen to be slightly underpredicted by TGLF SAT3 [Fig. 12(h)], which is consistent with the overprediction of the transport fluxes at the edge. As per the electron density gradients [Fig. 12(c)], 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[Fig. 12(f)], 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 Fig. 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 in performance with the increase in impurity content, as shown in Fig. 13. In CGYRO, it was observed that *Q* increased by $ \u223c 38 %$ from the nominal $ Z eff = 1.5$ case to the $ Z eff = 2.5$ case, while in TGLF *Q* increased by $ \u223c 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 \theta \rho s \u2264 24.0$ instead of $ k \theta \rho s \u2264 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 \theta \rho s \u223c 2 \u2013 3$ at the larger radii, even in situations with $ max ( \gamma / k \theta \rho s ) high \u2212 k < max ( \gamma / k \theta \rho s ) low \u2212 k$, a zonal flow mixing criterion usually employed to assess the relative influence of multi-scale effects.^{63} 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.

## VI. DISCUSSION

Sections IV and V have presented and discussed the results of the ion-scale nonlinear gyrokinetic and quasilinear gyro-fluid 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 \u223c 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) edge 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 reactor-like plasmas.^{50} 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 \u223c 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 \u2212 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.^{42} This high power to be transported through $ r / a = 0.943 \u2212 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, 8 *T*^{64}) 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,^{65} 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.6 \u2009 keV$) shown in the database of Alcator C-Mod plasmas in Fig. 2. The rather large ion to electron heat flux ratio seen in Fig. 8 (even at low Greenwald fractions, below what could be estimated to be the transition to the low-density branch for mode transitions^{66}) 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.^{67} 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*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}*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.^{56} 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,^{68} 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. 69 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 fueling and residual momentum stresses are available in a fully predict-first, physics-based fashion.

## VII. 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 $ \u223c 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., Refs. 70–72). It is only by the careful validation^{73} 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.^{74} 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 toward a future powered by clean, sustainable fusion energy.

## ACKNOWLEDGMENTS

The authors would like to thank M. Greenwald for the experimental data used in Fig. 9. This work was funded by Commonwealth Fusion Systems under RPP020, and used resources of the National Energy Research Scientific Computing Center (NERSC), a U.S. Department of Energy Office of Science User Facility located at Lawrence Berkeley National Laboratory, operated (Contract No. DE-AC02-05CH11231), for the PORTALS-CGYRO simulations (CGYRO from main branch with hash 016ac47, PORTALS version 1.0.0). Clusters hosted at the Massachusetts Green High Performance Computing Center were used to perform the PORTALS-TGLF simulations (TGLF from main branch with hash 037b2f0). Alcator C-Mod data used in this work were obtained under DOE Award DE-FC02-99ER54512. ChatGPT-4 and Microsoft Copilot were used to enhance parts of the manuscript for clarity and coherence.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**P. Rodriguez-Fernandez:** Conceptualization (lead); Data curation (equal); Formal analysis (lead); Investigation (lead); Methodology (lead); Software (equal); Visualization (lead); Writing – original draft (lead). **N. T. Howard:** Conceptualization (equal); Data curation (equal); Formal analysis (equal); Funding acquisition (lead); Investigation (equal); Methodology (equal); Project administration (lead); Writing – review & editing (equal). **A. Saltzman:** Data curation (supporting); Formal analysis (supporting); Supervision (supporting). **L. Shoji:** Formal analysis (supporting). **T. Body:** Formal analysis (supporting); Writing – review & editing (supporting). **D. J. Battaglia:** Formal analysis (supporting); Writing – review & editing (supporting). **J. W. Hughes:** Data curation (supporting); Formal analysis (supporting); Writing – review & editing (supporting). **J. Candy:** Software (lead); Writing – review & editing (supporting). **G. M. Staebler:** Software (lead). **A. J. Creely:** Project administration (equal); Writing – review & editing (supporting).

## DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request.

## REFERENCES

*TRANSP*