Gyrokinetic simulation of turbulence and transport in the SPARC tokamak

The world fusion program is poised to enter a new era. With ITER on the horizon and a host of new fusion experiments proposed to reach burning plasma conditions, there is an increased need to apply our most advanced computational models to predict device performance and even influence their design and operation. Years of plasma theory and computation, along with a world-wide effort in the validation of the gyrokinetic model,^1–8 have provided confidence that gyrokinetics represents an accurate model for transport and turbulence in the core of fusion devices. Commonwealth Fusion Systems, a Cambridge, MA based startup company, has proposed an ambitious plan to design, construct, and operate SPARC, a high field (12.2 T), mid-sized (major radius, R = 1.85 m, and minor radius, a = 0.57 m), tokamak by the mid 2020s with the explicit mission of demonstrating Q = 2 with 50 MW of fusion power generated.⁹ Recent work utilizing the ITER design rules and integrated modeling using the TRANSP code¹⁰ confirms that SPARC should achieve a plasma gain (Q) in the range of $9–11$ when operating the Primary Reference Discharge (PRD),¹¹ a low q95 (safety factor at 95% of minor radius $∼3.3$⁠), inductive ELM-y H-mode scenario, similar to the ITER baseline scenario. The SPARC PRD is projected to operate with dimensionless parameters already obtained in current devices.⁹ However, these parameters have only been obtained simultaneously in a handful of tokamak discharges performed to date. Nevertheless, due to the dimensionless similarity, it is anticipated that core transport will exhibit characteristics that are in many ways similar to current devices. SPARC modeling efforts to date were based on 0D modeling and 1.5D modeling that uses physics-based transport models, such as TGLF SAT-1¹² for its calculation of turbulent transport. In order to probe the nature of the turbulence and transport in SPARC and provide some insight into the accuracy of the 1.5D modeling, we present the analysis of linear and nonlinear gyrokinetic simulations for the SPARC PRD using the CGYRO code.^13,14 This work builds off of recent results by utilizing the kinetic profiles obtained from TRANSP modeling¹¹ as the starting point for the gyrokinetic analysis. This analysis does not predict the performance of SPARC directly, as self-consistent profile prediction using gyrokinetics will require extreme computational resources. Instead, it represents an investigation of SPARC conditions with gyrokinetics and helps to answer two questions: (1) Does the turbulence and transport in SPARC look similar or dissimilar to current day fusion devices and (2) What physics needs to be captured by gyrokinetic simulation to enable accurate prediction of plasma profiles and performance. This paper represents one of only a handful of attempts at nonlinear gyrokinetic simulation of a burning plasma performed to date. Compared with previous publications,^15,16 which were performed on projected ITER profiles, this work simulates positions across a large fraction of the plasma minor radius and investigates the potential role of a wide range of turbulence that has been observed in current day fusion devices. The remainder of this paper is organized as follows: Sec. II describes the simulation setups used throughout the paper. Section III covers linear and nonlinear ion-scale simulations of SPARC that attempt to match the power balance ion heat fluxes. Section IV describes sensitivity scans performed around the ion heat flux matched conditions, and Sec. V probes the importance of electromagnetic turbulence, rotation, and electron-scale turbulence for the prediction of SPARC performance. Finally, Sec. VI will discuss conclusions and future work.

All of the simulations described in this paper were performed using the continuum gyrokinetic code, CGYRO.^13,14 CGYRO is a spectral, modern gyrokinetic code that has a wide user base and was designed specifically for more efficient execution of near-edge/pedestal and multi-scale gyrokinetic simulation. Linear and nonlinear gyrokinetic simulations were all performed at three radial locations through this paper, $ρ=$ 0.4, 0.6, and 0.8 (square root of normalized toroidal flux) to probe core SPARC turbulence that exists from approximately the q = 1 surface to the plasma near-edge. All simulations used inputs obtained from the integrated modeling efforts described in Ref. 11, and a table of simulation inputs can be found in the  Appendix. Simulations utilized realistic geometries as represented by Miller parameterization of the flux surface shape,¹⁷ the Sugama collision operator,¹⁸ and a total of 6 gyrokinetic species [D, T, a lumped, equivalent impurity (Z = 9, A = 18), W (partially ionized, Z = 50, A = 184), He³, and electrons]. It is worth noting that these simulations did not include alpha particles. This choice was made due to the extremely low densities at the radial locations simulated (outside of $ρ=0.4$⁠) and the small effect that this population has on the plasma effective charge (Z_eff). Unless specifically noted otherwise in the text, all simulations were electromagnetic, evolving both the perturbed electrostatic and vector potentials $δϕ$ and $δA||$⁠. The nominal simulation domain used for nonlinear, ion-scale simulations was ∼77 × 64 ρ_s with $∼300$ radial modes and 16 toroidal modes spanning from $kθρs=$ 0.098 to $∼1.5$⁠, where $kθ$ is the poloidal wavenumber of the turbulence and ρ_s is the Larmor radius evaluated with the ion sound speed (c_s). Therefore, these simulations are capable of capturing long wavelength turbulence, such as ion temperature gradient driven modes (ITG), trapped electron modes (TEM), and micro-tearing modes (MTM), while ignoring contributions from high-k TEM and electron temperature gradient driven (ETG) turbulence. Electron-scale simulations are also covered in this paper in Sec. IV. These simulations used simulation domains of approximately 12 × 8 ρ_s represented by $∼360$ radial modes and 54 toroidal modes spanning from $kθρs∼0.75$ to 41.0 to enable the simulation of high-k TEM and ETG turbulence, including radially extended structures, such as ETG streamers. All simulation resolutions were represented by $Nξ=16$ (pitch angle grid points), $Nθ=24$ (poloidal grid points), and $Nenergy=8$ (energy grid points). Ion-scale simulations were typically run for $∼1200, ∼900, ∼550 a/cs$ at $ρ=0.4$⁠, 0.6, and 0.8, respectively, and typical time averages with typical time averages of 400, 300, and 300 $a/cs$⁠. The estimated uncertainty in the mean heat fluxes in this paper is $±10%$ as determined from the simulation time history. The ion-scale simulation domain was chosen to minimize the computation cost of each simulation while still maintaining sufficient accuracy of the results. Nonlinear convergences tests were performed in box size and radial modes to ensure convergence in these parameters. Results from these tests suggest that increasing box size to ∼100 × 100 ρ_s or increasing the radial resolution by 50% does not lead to a statistically meaningful change in the simulated heat fluxes. Linear spot checks were performed to ensure that other resolutions were numerically converged. The power balance heat fluxes quoted in this paper were the result of TRANSP calculations described in the integrated modeling in this Ref. 11.

As a first step of this analysis, linear CGYRO simulations were performed at the three radial locations of interest, $ρ=0.4,0.6,0.8$⁠. Scans were performed that span both the ion and electron-scales from $kθρs=$ 0.1 to $∼kθρs=$ 30.0 (⁠$kθρe∼0.5$⁠). Linear simulations were performed with both electrostatic and electromagnetic fluctuations evolved to demonstrate the validity of using an electrostatic approximation. The results of this linear analysis are plotted in Fig. 1. The results plotted in Fig. 1, along with a comparison of SPARC PRD linear growth rates from TGLF-SAT1 and CGYRO, were previously described in Ref. 11. The reader is referred to this reference for additional information on the code comparisons, as they are out of the scope of this manuscript.

Microtearing turbulence is shown to be unstable at long wavelengths at both $ρ=0.4,0.6$⁠, while ITG is the dominant instability at slightly higher wavenumbers (⁠$kθρs∼0.5$⁠) independent of the radial location studied. The role of electromagnetic turbulence on the linear growth rate appears to be minimal at both $ρ=0.6$ and 0.8 but results in a clear reduction in the ITG linear growth rates in the range $kθρs∼0.5−1.0$ at $ρ=0.4$⁠. High-k (⁠$kθρs>1.0$⁠) TEM is unstable at intermediate wavenumbers with ETG unstable at some level at all radial locations. As is commonly observed in inductive scenarios on present tokamaks, the low- and high-k growth rates increase as one moves closer to the edge. However, it is notable that the high-k TEM and ETG growth rates found at all of the radial locations are not large when compared with the low-k growth rates and, therefore, the impact of the electron-scale turbulence and multi-scale interactions is not expected to be large in these conditions for the base profiles.^19,20 Multi-scale simulations are outside the scope of this paper. However, electron-scale investigations presented in Sec. IV likely support the conclusion that ETG plays a minor role in the SPARC PRD.

With the observation from linear simulations that ITG is the dominant low-k instability at all radial locations, nonlinear, ion-scale simulations of the SPARC PRD were performed at $ρ=0.4,0.6,$ and 0.8. The simulation of the input profiles obtained from TRANSP predictive modeling demonstrates that the simulated heat fluxes are in disagreement with the power flows predicted via TRANSP. Given the stiff nature of transport and the differences between TGLF-SAT1 and CGYRO, this discrepancy is perhaps unsurprising. To resolve the discrepancy between power balance ion and electron heat fluxes and CGYRO, we performed multi-point scans of the ITG drive term, $a/LTi$⁠. The results from this $a/LTi$ scan are presented in Fig. 2, and this section provides a more in-depth description of the simulation results first described in Ref. 11.

There are a few clear observations from Fig. 2 that should be noted. The ratio of the ion to electron heat fluxes (⁠$Qi/Qe$⁠) calculated by TRANSP is in the range 2.0 to 3.0 over the radial region studied. These conditions are, therefore, unlike those obtained in current day, dominantly electron-heated plasmas, which tend to exhibit a $Qi/Qe<1.0$⁠. In fact, $Qi/Qe$ ratios much greater than 1.0 are generally found in predominately Neutral Beam Injection (NBI) heated plasmas, which tend to be dominated by low-k turbulence.²¹ ITG is clearly playing an important role in determining the transport in the SPARC core as relatively small changes in input gradients lead to substantial changes in the simulated heat fluxes at most radial locations. At $ρ=0.4$⁠, a significant increase (⁠$∼40%$⁠) in $a/LTi$ is needed to match the TRANSP simulated heat fluxes. In contrast, extremely stiff transport is observed at $ρ=0.6$⁠. At this radial location, once the value of $a/LTi$ is increased to the point that ITG turbulence is unstable, small changes in the inputs result in dramatic changes in the heat fluxes. This is demonstrated by the fact that a 5% increase from the input value of $a/LTi$ can take the simulated ion heat flux from zero to approximately 5× the experimental value. It is clear from the simulation that only a $∼1.5%$ increase in $a/LTi$ would likely be required to match the TRANSP power flows. The strong sensitivity of the heat fluxes to input gradients is a result of a strong Dimit's shift present at this radial location.²² At the nominal input value of $a/LTi$ (0%), ITG is linearly unstable and robustly non-zero heat fluxes are predicted early in the nonlinear simulation (see Fig. 3). However, the zonal flows grow throughout the simulation, leading to the suppression of the finite-n turbulence and a slow decrease to near zero heat fluxes as demonstrated in Fig. 3.

It is not clear whether extending these simulations would eventually lead to another phase of finite turbulence, but such investigation is out of the scope of this work. At the outermost simulated radii, $ρ=0.8$⁠, the response to $a/LTi$ is notably less stiff. Although the most unstable linear modes at this location were in fact ITG modes, it will be demonstrated through sensitivity scans that turbulence at $ρ=0.8$ is also near threshold for long wavelength TEM and likely exhibits some mixed mode characteristics. However, a 35% decrease in the ITG drive still results in nearly a 5× decrease in the predicted ion heat flux. Additional heating power applied to the SPARC PRD will result in relatively little change in the core ion temperature profile since small changes in gradients can easily exhausted the additional heat flux associated with increased heating. The presence of stiff transport in SPARC is an observation that would likely extend to other high performance plasmas, such as ITER's baseline operating scenario. In this work, an attempt was not made to match simulated heat fluxes to the exact power balance values. The intention of this paper is to demonstrate the nature of the turbulence and transport in the SPARC core without focusing on an exact quantitative reproduction of the TRANSP heat fluxes. Therefore, the sensitivity scans performed in Sec. IV will be performed around the conditions closest to being ion heat flux matched (+45%, +5%, and −35% $a/LTi$ at $ρ=0.4$⁠, 0.6, and 0.8, respectively).

The electron particle fluxes obtained from the scan of $a/LTi$ are plotted in the right panel of Fig. 2. Since SPARC is operated without neutral beams and with a very high line-integrated density leading to a negligible core source from edge gas fueling, the steady state electron particle flux should be approximately zero. However, we can see that CGYRO simulation predicts notably non-zero values of the electron particle flux for the conditions best matching the ion heat fluxes as shown in Fig. 2. Outward particle fluxes are predicted at $ρ=$ 0.4 with stronger, inward particle fluxes predicted at the outer radii simulated. In all simulations described in this paper, the deuterium and tritium particle fluxes are found to be approximately 1/2 of the electron particle flux value (⁠$ΓD≈ΓT≈0.5Γe$⁠). Naively the inward particle fluxes might be thought to imply that more peaked density profiles would be expected in the SPARC PRD. However, the interconnected nature of all the transport channels and drive terms makes this conclusion much less clear. This underscores the need for self-consistent profile prediction, where all kinetic profiles are modified in a manner such that all transport channels are accurately reproduced. Such work is out of the scope of this paper. However, in Sec. III, we will perform scans of multiple inputs around the approximately ion heat flux-matched conditions to understand the sensitivity of the nonlinear results of scans in input parameters.

The completion of ion heat flux-matched nonlinear gyrokinetic simulations provides an opportunity to probe low- and high-Z trace impurity transport in the SPARC PRD. These studies are of importance as the accumulation of impurities in a burning plasma can lead to excessive radiative losses, fuel dilution, and overall decrease in the fusion performance. Building off of ion heat flux-matched conditions described above, we performed six additional ion-scale simulations. Each of the ion heat fluxed matched conditions performed at $ρ=$ 0.4, 0.6, 0.8 was re-run only changing the gradient scale length of one of the target impurity species [the lumped impurity (Z = 9, A = 18) and partially ionized W species]. These impurity species were chosen because they were already present in the original simulation, and they represent reasonable proxies for two first wall materials being considered for SPARC [C (Z = 6) and W]. The particle fluxes obtained from each pair of simulations (each pair consists of two identical runs with different $a/Lnz$ values used) were then normalized to the impurity density and plotted vs $∇n/n$⁠. This results in the following expression: $Γ/n=−D∇n/n+V$⁠. The impurity diffusion and convection can, therefore, be determined by a linear fit to the data. This approach has been used in many previous publications.^21,23 The results from these simulations are plotted in Fig. 4.

It is generally observed that both species exhibit relatively small values of impurity diffusion (⁠$∼1.0$ m²/s or less) that peaks near $ρ=0.6$⁠. The pinch direction changes from inward to outward for the tungsten transport at large major radius, while the low-Z convection remains inward at all radii. The implication for a steady state profiles of these impurities is shown in the plot of the impurity peaking factor, V/D. These results imply a slightly hollow W profile near the edge, which is peaking near the core. In contrast, the predicted low-Z profile is slightly peaked throughout the profile. Despite the predicted peaking, the overall peaking factors found here correspond to $a/Lnz$ values in the range $∼−1.0$ and 3.5 and, therefore, do not imply that any dangerous impurity accumulation would be expected in this radial region. Since it is well known that neoclassical impurity transport can play an important role in the peaking of impurities, the NEO code²⁴ was also employed to evaluate the neoclassical transport coefficients at $ρ=0.4,0.6,0.8$⁠. The neoclassical diffusion and convection at these locations are found to be negligibly small (⁠$<5%$⁠) compared with the turbulent contributions, consistent with observations for low- and mid-Z impurities many present-day tokamaks outside of the deep core.^23,25 However, the relatively small neoclassical contributions are in contrast with recent W impurity transport studies on JET.⁴⁶ The reason for this discrepancy is likely due to significant differences in the discharges studied, as the JET results focus on a strongly beam driven (relatively high rotation), hybrid scenario that exhibits relatively larger electron density gradients. Both the high levels of rotation and larger density gradients found in the JET discharge can significantly enhance W neoclassical transport contributions. The JET conditions are in stark contrast to the intrinsically rotating, standard ELM-y H-mode conditions of the SPARC PRD (note that no rotation was assumed for the base simulations), where the electron density gradients are small over most of the profile (particularly relative to the temperature gradients). More exhaustive studies into the differences between JET and SPARC impurity transport are out of the scope of this paper but may be the subject of future work.

To better understand the nature of the transport in these conditions, we performed six nonlinear simulations at each radial location (18 total simulations). These additional simulations were in the form of ±20% changes in the turbulence drive terms, $a/LTe, a/Ln$⁠, and ν from their base values obtained from the TRANSP simulation output, which used TGLF-SAT1 to predict the kinetic profiles. These quantities are the normalized electron temperature gradient scale length, the normalized electron density scale length, and the collision frequency. The collision frequency was scaled in a manner such that all collision rates (e–e, e–i, etc.) were modified by the same percentage. As noted in Sec. III, these scans were performed around the $a/LTi$ values that were found to best match the power balance ion heat flux in Fig. 2. In Figs. 5 and 6, we plot the results of the these scans. The objective of these scans is to demonstrate the sensitivity of the results to changes in drive terms, to provide some insight into the types of turbulence responsible for particle and heat transport, and to provide an indication of the steady state profiles that may be predicted if nonlinear CGYRO simulation were used to predict the self-consistent steady state profiles in the SPARC PRD.

Figure 5 plots the sensitivity of the simulated heat fluxes to changes in inputs. The rows of Fig. 5 from top to bottom plot the results at $ρ=0.4,0.6,0.8$⁠. At $ρ=0.4$⁠, it is found that only moderate changes in the ion and electron heat fluxes result from $±20%$ changes in the drive terms. The changes include a small but finite increase/decrease in the predicted ion and electron heat fluxes with increased/decreased $a/Ln$⁠, a slight decrease in fluxes at increased $a/LTe$⁠, and flat or slightly decreased fluxes with increased ν. The ratio of ion to electron heat flux does not show significant variation over the scan, moving only from a value of approximately 5.4 to 4.7, suggesting that the dominant turbulence type (ITG) is largely unaffected by changes in the inputs at this location. These $Qi/Qe$ ratios also are found to significantly exceed the power balance values, which were $∼2.0$ at this radial location. At this time, it is unclear how this discrepancy would be resolved, perhaps requiring multi-parameter sensitivity studies to gain insight.

At $ρ=0.6$⁠, variations in ν do not show a clear trend within these parameter variations. Likely most of the observed variation is explained via uncertainties in the simulated heat fluxes due to the intermittent nature of the near marginal turbulence found at this radial location. Variations of $a/LTe$ and $a/Ln$ show much different behavior. The reduction in these input values led to essentially a complete stabilization of the turbulence, resulting in approximately 0 ion or electron heat flux. This behavior arises due to the stiff, near marginal state at which turbulence exists at this radial location. Small variations in the simulation inputs can tilt the balance of finite-n fluctuations and zonal flows such that the turbulent state collapses into a zonal flow dominated condition. The behavior of the simulations results at $ρ=0.8$ is straightforward as this condition exists farther from a marginal turbulent state. A clear trend of increased/decreased ion and electron heat fluxes with increased/decreased values of $a/Ln$ is observed, while the behavior $a/LTe$ exhibits a increase in fluxes at increased values, potentially indicating the onset of $a/LTe$ driven TEM at this location. The $Qi/Qe$ ratios at this location show a larger variation during the parameter scans. There is a clear difference observed with $a/Ln$ variation, taking the $Qi/Qe$ ratio from approximately 3.0 to 1.5 over the $±20%$ scan, suggesting a significant change in the turbulence during the scan.

The electron particle fluxes obtained from the sensitivity scans are plotted in Fig. 6. These results demonstrate the significant sensitivity of the simulated particle fluxes to rather modest changes in the inputs. At $ρ=0.4$⁠, the sensitivity is less dramatic with changes in $a/LTe$ and $a/Ln$ representing the largest knobs for changing the overall particle flux. Increased $a/Ln$ or decreased $a/LTe$ brings the predicted fluxes much closer to a steady state value of 0. The particle fluxes at $ρ=0.6$ show significant variation over the scan. The $±20%$ change in collision frequency at this location can change the predicted particle flux level by a factor of 3. A larger variation in the flux is seen with $a/Ln$ variations. However, it should be noted that the reduced $a/Ln$ values led to the stabilization of the turbulence as noted in the discussion of the heat flux sensitivities. At $ρ=0.8$⁠, a strong $a/Ln$ sensitivity is found, with a 8× increase in the inward flux when going from −20% to +20% increase in the value of $a/Ln$⁠. The sensitivity of the results to changes in $a/LTe$ and ν at this location is weak. Combined with the results plotted in Fig. 5, it would appear that a 20% reduction in $a/Ln$ might allow for both particle and heat fluxes to be approximately matched. The simple interpretation of Fig. 6 suggests that the flattening of the density profile is anticipated at $ρ=0.8$ with increased density peaking anticipated at $ρ=0.4$⁠. However, it is important to emphasize that this analysis underscores the coupled nature of the transport problem and motivates the self-consistent modification of all the kinetic profiles to obtain conditions that match target heat and particle fluxes.

One objective of the simulations performed in this paper was to understand the physics fidelity needed to perform gyrokinetic simulations that represent the turbulent state in the PRD of the SPARC tokamak. Due primarily to computational limitations, gyrokinetic simulations of plasmas are often performed with reduced physics fidelity to enable more tractable simulations. This generally involves reducing the temporal and spatial scales captured in the simulations and/or disabling specific physical processes thought to play a less important role. In this section, we examine the impact of electromagnetic fluctuations, rotation, and electron-scale turbulence on the SPARC gyrokinetic simulations.

In order to investigate the effect of including electromagnetic (EM) turbulence on the simulated heat fluxes, we performed a set of electrostatic runs to complement the existing ion heat flux-matched, electromagnetic simulations described in Sec. III. These simulations were identical with the exception that only $δϕ$ was evolved in the gyrokinetic equation, as opposed to keeping both the $δϕ$ and $δA||$ components as was done in all of the simulations previously described. Recall that some attempt to understand the electromagnetic nature of the turbulence was discussed in Sec. III and is plotted in Fig. 1. In Sec. III, it was found that the ITG linear growth rates at $ρ=0.6$ and 0.8 were essentially unchanged by enabling electromagnetic effects. This is in contrast to the non-negligible reduction in the ITG growth rates due to electromagnetic effects observed at $ρ=0.4$⁠. Additionally, it was shown that microtearing modes are linear unstable at the lowest simulated k values at both $ρ=0.4$ and 0.6 (Fig. 1). The electron heat fluxes calculated from nonlinear, ion-scale simulations of the approximately ion heat-fluxed matched conditions are shown in Fig. 7.

We plot Q_e since the inclusion of electromagnetic turbulence might be expected to contribute primarily to the changes in the electron heat flux, and although ion heat fluxes are not plotted here, they show the same qualitative conclusions. Importantly, it is found that the direct contribution to the electron heat flux from the $δA||$ component of the fluctuations is small (⁠$<10%$⁠) and negative at all radial locations. Small negative contributions to the heat flux have been demonstrated to be a characteristic of electromagnetic ITG previously.²⁶ The lack of heat flux driven directly by $A||$ components of the fluctuations suggests that, despite their linear instability at $ρ=0.4$ and 0.6, microtearing modes are playing a negligible role in these plasma conditions and that the effects observed are primarily due to the stabilization of the ITG. This finding is largely in agreement with previous work that demonstrated a negligible role of MTM transport in the presence of ITG.²⁷ The linear results in Fig. 1 indicate a reduction in the ITG growth rates at $ρ=0.4$ when comparing electrostatic and electromagnetic simulations. At this radial location, the stabilization of the ITG resulting from going from electrostatic to electromagnetic nonlinear simulation reduces the heat flux level by approximately an order of magnitude. However, we note that the ion to electron heat flux ratio is essentially unchanged between the two simulation types, suggesting that there was no notable change in the turbulence type, simply the overall level was modified. The large decrease in ITG transport resulting from the inclusion of electromagnetic effects, despite a relatively modest reductions in growth rate, has been observed previously²⁸ and was recently explained via a decrease in energy transfer efficiency between stable modes, unstable modes, and zonal flows.²⁹ At $ρ=0.6$ and 0.8, the simulated heat fluxes are nearly identical independent of the simulation type. This is largely consistent with the unchanged growth rates shown in Fig. 1. For the simulation at radial locations where the electromagnetic and electrostatic simulations yield similar results (⁠$ρ=0.6$ and 0.8), there does not appear to be a significant difference in the electron particle fluxes. The overall conclusion of this analysis and its implications for the simulation of the SPARC performance is that (1) microtearing plays a negligible role for these particular conditions and (2) electromagnetic contributions must be included in SPARC simulations due to the potential stabilization of the ITG and the resulting influence on the simulated heat fluxes.

The role of rotation, both intrinsic and externally driven, in setting and regulating turbulence in fusion plasmas has been well established,³⁰ and as a result, a clear link between rotation and confinement has been identified. However, there is currently much debate about the ability of gyrokinetic modeling to accurately simulate momentum transport and the resulting rotation. Due to this uncertainty in predictions of momentum transport and the resulting intrinsic rotation, the rotation profile was set to 0 for all SPARC TRANSP analysis¹¹ and the gyrokinetic simulations described above. An empirical model for the predicted rotation profile in SPARC and performed nonlinear CGYRO simulations was used to assess the potential impact on gyrokinetic based predictions. The model chosen for this work is based on the work published by Kosuga, which is described here in Ref. 31. The Prandtl number was taken to equal to 1, consistent with Ref. 32 and with theoretical expectations for ITG dominated plasmas.³³ From this model, we have obtained a potential range of Mach numbers for the SPARC PRD, which span from a value of 0.07 to 0.25. The profile of the rotation was set so that the Mach number is constant across the profile, which is consistent with data from the Alcator C-Mod tokamak.³² This observation seems relevant since Alcator C-Mod was operated without external momentum input and exhibited only intrinsically driven rotation—a situation that is identical to SPARC due to its lack of neutral beam heating. Three profiles were created with Mach numbers of 0.07, 0.16, and 0.25 as shown in Fig. 8.

Using the three different rotation profiles shown in Fig. 8, we performed a set of nine nonlinear gyrokinetic simulations, spanning both the Mach number values (0.07, 0.16, 0.25) and radial location (⁠$ρ=$ 0.4, 0.6, 0.8). The resulting ion heat fluxes obtained from these simulations are plotted in bottom frame of Fig. 8. These nonlinear simulations used the identical numerical setup to those described in Secs. 3 and 4 (ion-scale, EM turbulence) with all of the CGYRO rotation effects enabled (Mach number, rotation shearing, and ExB shearing). The results at $ρ=0.4$ and 0.8 are relatively clear. At $ρ=0.4$⁠, there appears to be very little change in the simulated heat fluxes as the Mach number is scanned from 0 to a value of 0.25, suggesting that even a relatively large change in the rotation profile will likely result in only minor changes in the predicted turbulence and the resulting performance of SPARC. In contrast, at $ρ=0.8$⁠, the effect of rotation is much more dramatic. At small values of Mach number (0.07), there is a marginal reduction in the driven heat flux but at values of Mach number of 0.16 and above the turbulence is effectively stabilized by the rotation, leading to essentially no driven ion heat flux. This is a particularly important observation in the context of the stiffness displayed in Fig. 2. While the ITG transport appears quite stiff at inner radii, the stiffness is reduced in the near-edge region at $ρ=0.8$⁠. Therefore, the effect of rotation stabilizing turbulence will be much more pronounced in this region. The stabilization of the turbulence via ExB shearing will lead to a modification of the gradients at $ρ=0.8$⁠, which will, in turn, result in significantly different predicted temperature profiles.

The results at $ρ=0.6$ are less clear. It should be recalled from Fig. 2 that this radial location was the closest to marginal stability, existing only a $∼5%$ change in $a/LTi$ from the ITG critical gradient. As a result, it should be expected that only small variations in input parameters in this condition may lead to significant changes, as the interplay between a wide range of mechanisms becomes important when the turbulence exists in such a marginal state. When rotation is added to the simulation, there is both a destabilizing effect from rotation shearing and stabilizing effect from E× B shearing with both influencing the heat flux level produced by the simulation. When the lowest Mach number profile is used in simulation in these conditions, the plasma collapses back into a state with stabilized turbulence (similar to Fig. 3). At higher values of Mach number, the simulation predicts a non-zero heat flux level that is reduced (relative to the simulation neglecting rotation) as Mach number increases. The overall conclusion of these simulations is that rotation will likely play an important role in predictions of profiles and performance in the SPARC PRD. More specifically, it appears as though significant changes in the SPARC profile due to rotation will likely occur in the outer regions of the plasma, where the turbulence is less stiff, and the effect of rotation is more dramatic. However, as discussed above, the accurate simulation of momentum transport is an area of active research and debate within the plasma physics community without a clear resolution. At this point, the prediction of SPARC profiles will continue to take a conservative approach, setting rotation to zero, with the general assumption that the inclusion of rotation will result in increased profile peaking and performance.

In the simulations discussed in all other sections of this paper, we chose to focus exclusively on the simulation of turbulence existing with the values of $kθρs≲1.5$⁠. However, the role of multi-scale (coupled ion and electron-scale) turbulence in fusion plasmas has been well documented in recent years^34–42 and it is not clear a priori whether or not electron-scale turbulence and the resulting multi-scale interactions will play an important role in SPARC. The existence of significant multi-scale interaction would dramatically increase the computational resources needed for accurate simulation of the SPARC tokamak and would likely make profile prediction with gyrokinetic simulation intractable for the foreseeable future. In this section, we will perform electron-scale gyrokinetic simulations at $ρ=0.4,0.6,0.8$ to attempt to understand the impact that electron-scale turbulence and estimate the role of multi-scale interactions, in setting the heat fluxes in the SPARC PRD. The setup of the simulations performed here is described in Sec. II, but as a reminder, simulations were performed such that they resolve wavenumbers in the range $∼0.7<kθρs<41.0$ with box sizes sufficiently large to capture structures such as ETG streamers. The results from nonlinear, electron-scale simulations at the three radial locations studied are shown in Fig. 9.

It is shown in the top panel of Fig. 9 that the peak high-k electron heat fluxes predicted by nonlinear simulation decrease significantly moving toward the plasma core. This observation is generally consistent with the linear growth rates plotted in Fig. 1. Simulations at $ρ=0.6$ and 0.8 display essentially two phases. To best illustrate the differences in these phases, they are labeled 1 and 2 on the time history of the simulation performed at $ρ=0.8$ and the electron density fluctuations corresponding to these two times are shown in the bottom two panels of Fig. 9. Early in the simulation, ETG streamers form. These streamers are radially elongated and due to their radial extent, drive more substantial levels of electron heat flux. In fact, at $ρ=0.8$⁠, ETG streamers appear capable for driving a level of electron heat flux (⁠$∼0.2$ MW/m²) that exceeds that found in the TRANSP simulation (level shown in Fig. 2). At $ρ=0.6$⁠, the first phase exhibiting ETG streamers is pushed later in time and temporarily saturates at a much lower value (⁠$∼0.03$ MW/m²). This value is smaller than the TRANSP simulated heat flux value by approximately 5×. In the second phase of the simulations at both $ρ=0.6,0.8$⁠, the turbulence becomes dominated by electron-scale zonal flows, which destroy the coherent radial structures of the ETG streamers and thus decrease the overall heat flux level to a negligible level. This type of behavior has been reported in several publications when performing electron-scale simulations.^19,39,41,43 At $ρ=0.4$⁠, the heat flux level is relatively constant throughout the simulation and is characterized by only small ETG streamers, driving less than 1/10 of the electron heat flux predicted by TRANSP, making its overall contribution negligible.

To understand the overall role of electron-scale turbulence (and any potential multi-scale interactions) in these conditions, we can use these electron-scale simulations in conjunction with the results in Fig. 2 and can build off of key observations from recent work. The important observations are: (1) if electron-scale simulation predicts a negligible level of electron-scale turbulence driven heat flux, it will probably not play a dominant role, even in multi-scale simulation^38,39,41 (2) if ion-scale turbulence is strongly driven (i.e., far from marginal) the impact of electron-scale turbulence and multi-scale simulations will likely also be negligible.^19,38,39 In the previous work, cases that have exhibited the most significant multi-scale interactions are cases with strongly driven ETG streamers and weakly unstable ion-scale turbulence. With these criteria in mind, at $ρ=0.4$ it appears unlikely that electron-scale turbulence plays an important role, as the electron-scale simulations display only a negligible level of heat flux. At $ρ=0.6$⁠, the ion-scale turbulence is extremely close to marginal (see Fig. 2), which is an ideal condition for multi-scale interactions to occur.¹⁹ However, as shown in Fig. 9, the level of electron heat flux, even prior to the collapse of the ETG streamers into a zonal flow dominated state, is only about 20% of the experimental level. This suggests that the role of electron-scale turbulence is still likely small at this radial location. The largest potential contribution from ETG turbulence occurs at $ρ=0.8$⁠. At this radial location, ETG streamers, before collapsing into a zonal flow dominated state, drive a level of electron-heat flux exceeding the TRANSP value. However, the ion-scale turbulence at this radial location is more strongly driven and farther way from its critical gradient than other locations. It has been demonstrated in multi-scale simulation that this type of situation reduces or even nearly eliminates ETG streamers, due to increased shearing effect generated by low-k turbulence.^19,44 The conclusion of these simulations is, therefore, that electron-scale turbulence likely does not drive a large amount of heat flux and coupled multi-scale interactions likely do not play a large role in the prediction of the SPARC PRD. This conclusion is supported by the linear rule of thumb discussed in Sec. II. However, we note that this is an active area of research and only nonlinear multi-scale simulations can likely answer this question definitively.

In this paper, we presented a set of gyrokinetic simulations of plasma turbulence and transport in the core of the SPARC primary reference discharge (PRD). These simulations indicate that the turbulence expected in SPARC is similar to that found in many current tokamaks. Namely, it appears that long wavelength, electromagnetic ITG plays the dominant role in the SPARC core. A brief look into the impurity peaking predicted via gyrokinetic simulation indicates that dangerous impurity accumulation is not expected in the region studied. Scans were performed in the SPARC PRD to determine the sensitivity of the gyrokinetic results to changes in other turbulence relevant quantities. These simulations demonstrated that some sub-dominant TEM may be present at $ρ=0.8$ and that $±20%$ variations of other inputs may significantly modify the balance of ion and electron heat flux predicted via simulation. It was shown that particle fluxes can vary significantly during relatively modest changes in inputs and that a self-consistent modification of all kinetic profiles was needed to recover the TRANSP predicted particle and heat fluxes simultaneously. The role of electromagnetic turbulence, rotation, and electron-scale turbulence was also explored via nonlinear simulation. These simulations indicate that electromagnetic effects do indeed play an important role in the deep core (⁠$ρ=0.4$⁠) of SPARC by leading to stabilization of the ITG turbulence. The inclusion of a modeled rotation profile was shown to play a negligible role at $ρ=0.4$ while exhibiting nearly full suppression of the turbulence at $ρ=0.8$ for the studied rotation profiles. The results at $ρ=0.6$ were found to be more complicated due to the near marginal nature of the turbulence and the interplay between the stabilizing and destabilizing effects of rotation. Electron-scale turbulence and multi-scale interactions were found to likely play a negligible role with the potential exception at $ρ=0.8$⁠. However, at this radial location, it is expected that the stronger ion-scale turbulence will likely serve to suppress the overall impact of the high-k turbulence found in electron-scale simulation.

One of the motivations for this work was to determine the physics fidelity needed to accurately simulate SPARC and similar high density burning plasma conditions via gyrokinetic simulation. The nonlinear ion-scale simulation results from Sec. III B, complemented by the sensitivity scans in Sec. IV clearly demonstrate that dominant role that electromagnetic ITG turbulence plays in the core (ρ = 0.4–0.8) of the condition studied with microtearing and TEM playing likely negligible roles. The sensitivity scans demonstrate that modifications to turbulence drive terms other than the primary ITG drive, $a/LTi$ (namely, $a/LTe$ and $a/Ln$⁠) can significantly modify the predicted heat and particle fluxes. Therefore, a self-consistent, multi-channel modification of all kinetic profiles (n_e, T_e, T_i) through a framework such as TGYRO⁴⁵ is likely needed to reproduce expected levels of ion and electron heat fluxes (Q_i, Q_e), the ratio $Qi/Qe$⁠, and electron particle fluxes (γ_e). The exact value of rotation in SPARC is difficult to predict, but it is clear that rotation will increase any performance predictions, implying that the assumption of zero rotation used so far is conservative. Although it was not clear a priori, linear and nonlinear simulation suggests that high-k turbulence, in the form of ETG streamers, will likely exist in the core of these burning plasma conditions, but will play a minimal role in setting heat and particle transport. As a result, high-k turbulence likely does not need to be simulated for accurate profile prediction. The conclusion from these combined findings is that electromagnetic (evolving $δϕ$ and $δA||$ fluctuations) ion-scale (⁠$kθρs≲1.0$⁠) nonlinear gyrokinetic simulation, performed with an accurate collision model (Sugama¹⁸ operator used in this work), impurities, rotation, and realistic geometry contains sufficient physics to accurately simulate the turbulence and transport expected in the high density burning plasma conditions represented by the SPARC primary reference discharge. Since such simulations can be performed with relatively low computational cost (<100k CPU hours), SPARC conditions could be viewed as ideal for the demonstration of profile prediction using gyrokinetic simulation, as the effects of high-k/multi-scale turbulence are generally expected to be small. The results in this paper provide the first insights into the turbulence expected in SPARC and similar high density burning plasma conditions and provide motivation to proceed with gyrokinetic profile prediction in future work.

We would like to acknowledge the support of the SPARC physics team for their feedback on this work. This research 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 under Contract No. DE-AC02–05CH11231. This work was funded by Commonwealth Fusion Systems under RPP005.

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

The local input parameters used for the simulations in this paper are indicated in Table I.