Reducing transport via extreme flux-surface triangularity

Based on a gyrokinetic analysis of and extrapolation from TCV discharges with large negative and positive triangularity δ, the potential of extreme |δ| in reducing turbulent transport is assessed. Linearly, both positive and negative δ can exert a stabilizing influence, with substantial sensitivity to the radial wavenumber kx . Nonlinear fluxes are reduced at extreme δ in a trapped-electron-mode regime, whereas low-amplitude ion-temperature-gradient turbulence is boosted by large negative δ. Focusing on the former case, nonlinear fluxes exceed quasilinear ones at negative δ, a trend that reverses as δ > 0. A change in saturation efficiency is the cause of these features: the zonal-flow residual is boosted at δ > 0, reducing fluxes compared with the linear drive as δ is increased, and a shift towards larger zonal-flow scales occurs with increasing δ due to finite-kx modes weakening with δ.


Introduction
Non-circular flux-surface cross-sections, achieved by plasma shaping, can have a substantial impact on plasma dynamics.Typical applications in tokamaks are magneto-hydrodynamic (MHD) stability [1,2] and control of microinstabilities and turbulence [3][4][5], while in stellarators, such shaping can also be one aspect of optimization for e.g.reduced neoclassical transport [6,7].
Original Content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence.Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.
One particular topic pertaining to shaped tokamaks is negative flux-surface triangularity, which entails promising confinement and power-handling properties [8].Investigations into the microinstabilities and turbulence characteristics of negative-triangularity (NT) configurations (see [9] for a review) have led to a number of findings.The trappedelectron mode (TEM) -much of the experimental work on NT has focused on this regime, often with mixed densitygradient and temperature-gradient drive-can be stabilized when the triangularity δ < 0, as the toroidal precession drift is stronger for deeply trapped electrons, reducing TEM growth [10].Relatedly, the critical electron temperature gradient is increased for δ < 0, particularly near the plasma edge, where |δ| tends to be larger; no substantial difference in profile stiffness is observed between NT and positive triangularity (PT) [11].In present-day medium-sized devices, the argument has been made that radially local TEM turbulence studies are more accurate at δ > 0 and can lead to transport overprediction in NT configurations [12].This, however, is unlikely to affect simulation of larger-scale devices, where it has been theorized that reducing TEM turbulence even in reactor-relevant regimes dominated by ion-temperature-gradient (ITG) turbulence may be beneficial [9]; note, however, the possible benefit of proximity to regime boundaries, where turbulent fluxes can be nonlinearly suppressed [13].
Using artificial geometry and profiles, [14] reports stabilization at δ < 0 not only for TEM but also ITG turbulence.This finding is consistent with the modification of local magnetic shear [15], which for low elongation κ ≈ 1 is reduced in magnitude by increasing δ, causing triangularity to act as a destabilizing force; as κ increases, the impact of local shear on δ is lessened, and δ increasing the perpendicular wavevector k ⊥ stabilizes ITG modes.This applies to linear growth rates and nonlinear ITG fluxes alike [4].While these ITG studies have employed kinetic electrons, the use of an adiabatic response for computational efficiency has allowed [16] to identify a locally optimal parameter point in κ-δ space at κ = 1.6 and δ = −0.625,where ITG activity is substantially reduced.Similarly, [17] has suggested substantial benefits of large values of |δ| → 1, beyond the values presently straightforwardly accessible in experiment [18].
Validation of gyrofluid-based quasilinear fluxes against experimental transport [19] confirms the applicability of mixing-length models for moderate magnitudes of δ, whereas for extreme δ, it remains to be ascertained whether such quasilinear models can successfully predict fluxes.
In the present work, the ramifications of extreme δ on instability and turbulence are assessed based on Ohmically heated, limited NT and PT discharges in the TCV tokamak.It is to be noted in this context that practical limitations may exist preventing present-day devices from achieving such δ, related e.g. to MHD stability or operational considerations.For the purpose of this study, such concerns are set aside for future work, focusing solely on the impact of extreme δ on microturbulence.
The remainder of this paper is organized as follows.In section 2, background is provided on the experimental TCV discharges studied here, including a transport analysis for later comparison with simulations.Linear and nonlinear gyrokinetic analysis in nominal geometry follows in section 3 to allow for comparison with experiment and as an anchor for the subsequent studies of section 4, which utilize parameterized geometry in order to be able to vary δ independently.There, instability and transport are assessed for different δ, as is mode localization, and the role of zonal flows in setting transport levels and bridging the gap between quasilinear and nonlinear fluxes is determined.Section 5 contains a summary of key results.

TCV experiments and transport modeling
Discharges attempting to maximize the magnitude of the edge δ were performed on TCV, resulting in the flux surface shapes and pressure (gradient) profiles shown in figures 1 and 2, respectively, with the normalized electron pressure β and the Cross-sections at constant toroidal angle of flux surfaces ρtor = 0.1 − 0.9 for the NT discharge (blue) and the PT discharge (red).The radial position 0.7 chosen for subsequent analyses is marked as a slightly bolder line.Electron temperature (solid) and density (dashed) profiles for the NT discharge (blue) and the PT discharge (red).The ion temperature is assumed to be T i0 = 0.133T e0 across the domain, implying identical normalized gradients ω T (dotted) and ωn (dash-dotted) for both species.The location ρtor = 0.7 is marked as a vertical black line.collision frequency ν ei (in units of the ion sound speed c s divided by the major radius R 0 ), as derived from the profile data, reported in figure 3, where key geometric quantities are also included.Electron temperature and density profiles were obtained from the Thomson scattering diagnostic [20], whereas ion temperature and charge state are based on the diamagnetic loop [21] due to a lack of reliable charge exchange recombination data.Throughout the radius, experimental errors are close to 10% of the corresponding density or temperature value; thus, the core temperatures between Figure 3. Left: profiles of the electron β (in units of 0.1%, solid) and the collision frequency ν ei (dashed) for the NT discharge (blue) and the PT discharge (red).Right: Profiles of safety factor q (solid), elongation κ (dashed), and triangularity δ (dotted) for the NT discharge (blue) and the PT discharge (red).The location where turbulence simulations are performed is marked as vertical black lines.
the two discharges are identical within error bars at less NT heating power, while the NT core density exceeds that of the PT discharge somewhat outside of error bars.The normalized gradients are defined as ω n,T = R 0 /L n,T , where L n,T refers to the density or temperature scale length, computed based on the radial coordinate ρ tor .Note that the color convention employed in these figures, where blue corresponds to negative and red to positive δ, is applied throughout this work unless stated otherwise.
The NT discharge #73564 and the PT discharge #73589 were matched closely, achieving a separatrix triangularity of approximately ±0.6, whereas the magnitudes at a normalized toroidal flux position ρ tor = 0.7 decrease to δ ≈ ±0.3.Subsequent instability and turbulence analyses pertain to this radius, centered around time 1.25 s of either discharge.
Even for the very flexible poloidal-field coil set of TCV, developing these high triangularities presented considerable challenges.In particular, the radial plasma position had to be carefully tuned to maintain the plasma limited on the inside wall alone and to prevent the appearance of unwanted X-points.In the NT case, in addition, a narrow path had to be found to avoid both low-rational-order MHD instabilities and axisymmetric (vertical) instabilities.Although the goal was to match the stored energy between the two discharges, the energy was approximately 20% higher in the NT case in spite of the input power (Ohmic only) being 23% lower, suggesting somewhat better core confinement in the NT case.
Across most of the device, the NT and PT profiles match reasonably well, including at the 0.7 position.While in the following, a brief comparison between experimental heat fluxes and those obtained from turbulence simulations is performed, the intent is not to validate the turbulence model, and any in-depth comparison is left for future study.
To facilitate comparison with experimental heat fluxes as a baseline for the applicability of the subsequent studies performed in this work, interpretative transport simulations were performed with ASTRA [22], where the temperature and density transport equations were solved to match the experimental profiles at any given time, accounting for sources.For this purpose, a simplified equilibrium was reconstructed Shown are the electron heat flux as a solid and the ion heat flux as a dashed line.Error bars indicate variability within a 10 ms time interval.Given the slightly higher core density for NT in figure 2, the NT confinement is slightly better but similar to the PT confinement.from the plasma boundary using LIUQE [23], and the TCV CU automatic module was used in solving the transport equation for the poloidal flux.Based on these ASTRA simulations, the heat flux in figure 4 was calculated, with error bars obtained by averaging over simulations within a 10 ms time interval during the discharge.Here, the normalized radial coordinate r/a is used; note that r/a = 0.68 corresponds to ρ tor = 0.7, which will be the focus of subsequent analysis, where after assessing the linear instability picture, the ASTRA fluxes will be compared to results from nonlinear gyrokinetic simulations.
As will be shown, variation of input gradients within error bars results in large changes to simulated fluxes, due to proximity to critical gradients and turbulence regime transitions.Experimentally, this property may be the result of the Ohmic heating and possibly of the push to strong shaping.As a consequence of the sensitive simulated fluxes, quantitative comparison with experiment merely provides an anchor, and investigations presented here focus on the underlying physics of the turbulence.
Table 1.Physical input parameters for which gyrokinetic simulations are performed for NT and PT numerical and parametric equilibria.In order, these parameters are the normalized pressure β, the electron-ion collision frequency ν ei , the effective charge Z eff , normalized density gradient ωn and temperature gradient ω T , the safety factor q 0 , the magnetic shear ŝ, and the device minor radius a 0 .

Instability and turbulence in nominal geometry
Based on these profiles, as part of the EUROfusion RT07 campaign, flux-tube simulations at ρ tor = 0.7 are performed with the gyrokinetic [24] turbulence code Gene [25,26] based on numerical equilibria obtained from CHEASE [27], first to ascertain the linear stability properties and then to examine turbulent heat fluxes.Note that at this position, |δ| ≈ 0.3 is substantially less than at the separatrix.For the NT/PT discharges #73564/#73589, the physical input parameters are defined and listed in table 1.Note that Z eff is used as part of a Landau collision operator with collisions among and between any species included; the primary constituent of this Z eff is a carbon impurity, which is not included as an explicit species.Simulations retain shear magnetic fluctuations but ignore parallel magnetic fluctuations.While experiments use 75% deuterium and 25% hydrogen, for computational efficiency, no hydrogen is included in simulations.Numerical input parameters are listed in table 2, which are the result of both linear and nonlinear convergence testing (where doubling resolutions, hyperdiffusivity, or box sizes does not significantly alter fluxes).As evident from figure 5, the PT and NT discharges are dominated by TEMs and ITG modes, respectively.Particularly in the ITG-dominated case, modes at finite radial wavenumbers k x (which, like the toroidal wavenumber k y , is normalized to the inverse of the ion sound gyroradius ρ s ) can attain significant growth rates γ.A so-called flipped NT scenario, where the NT geometry is combined with the density/temperature profile parameters (including ν ei and Z eff ) from the PT discharge, highlights the ability of finite-k x modes to affect the system when δ < 0, as finite-k x modes now grow at comparable rates to those at k x = 0; this change in balance of different k x coincides with a more mixed TEM-ITG state.The process underlying these changes is described, e.g. in [17], where ITG localization is a consequence of field-line curvature in combination with polarization (which is related to the squared perpendicular wavenumber k 2 ⊥ ): at negative triangularity, these quantities exhibit wide wells along the field line comparatively insensitive to k x , leading to mode localization and growth experiencing little impact from k x .Another scenario, where the Linear growth-rate spectra in realistic geometry.The linestyle encodes at which kx modes are centered: solid for 0, dotted for 0.1, dashed for 0.2, and dash-dotted for 0.4.Lines with/without crosses indicate negative/positive frequencies, corresponding to TEM/ITG.The PT discharge (red) has kx = 0 dominate across ky, whereas the NT discharge (blue) has near-identical contributions from all kx.A flipped scenario (pink) with NT geometry and otherwise PT parameters approaches the TEM-ITG transition and has reduced kx dependence; further changing this flipped scenario by using the collision settings from the NT discharge (black) leads to γ close to that of the NT discharge.
flipped case is further modified by switching in the NT collision frequency and effective charge, is also shown in figure 5, showing that the difference between NT and PT profile gradients alone cannot account for most of the changes between discharges.
The weak dependence of γ on k x for NT is compatible with highly efficient energy transfer from unstable to stable modes via zonal flows, leading to low saturated fluctuation amplitudes [29][30][31][32].Indeed, when performing nonlinear simulations, the heat flux-as shown in figure 6-for the NT discharge attains a very low amplitude in the quasi-stationary state, whereas the PT flux lies approximately a factor of two above the corresponding ASTRA value.While the present effort does not aim to undertake validation, given the similar non-δ input parameters, sensitivity of linear growth rates to moderate scenario changes, as well as the apparent Dimits regime [33]-i.e. the regime between the linear instability and nonlinear turbulence thresholds-of the NT simulation, it can be speculated that changes of physical input parameters such as profile gradients within experimental error bars will allow for matching observed fluxes in simulations.
Given the local experimental δ values of approximately ±0.3, these discharges in themselves are unable to provide the basis for assessing the potential of extreme δ with respect to turbulent transport.Thus, extrapolation to larger |δ| is required, for which the numerical equilibrium of the NT and PT discharges is converted into a parameterized geometry.

Extreme-triangularity simulations
Parameterized magnetic equilibria, most prominently the Miller geometry [34,35], allow for independent variation of individual shaping parameters such as δ.While different strategies exist regarding the conversion of a numerical equilibrium into parameterized Miller-like geometry, the approach used here is that of [36].Note that the fitting mechanism entailed therein slightly modifies some of the geometric quantities from their values in the numerical equilibrium, such that the full new set of resultant geometric parameters obtained for the NT/PT discharges #73564/#73589 reads δ = −0.303/0.280,q 0 = 1.57/1.50,ŝ = 1.16/1.38,a 0 = 0.267/0.272,along with the Shafranov shift α MHD = 0.146/0.164,the inverse aspect ratio of the ρ tor = 0.7 flux surface ϵ t = 0.182/0.200,the triangularity shear s δ = −0.460/0.493, the elongation κ = 1.17/1.12and its shear s κ = 0.026/0.143,the squareness ζ = 0.006/ − 0.017 and its shear s ζ = −0.049/− 0.033, and the flux-surface shifts dR/dr = −0.274/0.025and dZ/dr = −0.004/− 0.006.While such geometries miss certain features such as differences between upper and lower δ, they capture most of the relevant shaping properties of tokamak plasmas.Note that, consistent with typical profiles, the sign of s δ follows the sign of δ, whereas in the following analysis, δ will be varied across a large range while keeping other parameters, including s δ , fixed.However, when switching in s δ = −0.460for a nonlinear simulation otherwise matching the PT scenario, at δ = −0.6,only a small quantitative difference in fluxes is observed relative to s δ = 0.493, and all findings as reported throughout the remainder of this paper still hold for this data point.Furthermore, given that [37] reports strong sensitivity of ITG-driven fluxes to s δ in some regions of parameter space, the triangularity shear was set to a more extreme-but for large |δ| itself potentially relevant-s δ = 0.8 (for δ = 0.6 in the PT/TEM-dominated case) and s δ = −0.8(for δ = −0.6 in the NT/ITG-dominated case).Either test resulted in less than 10% change of ion or electron fluxes relative to the respective nominal s δ cases.Note that this finding does not contradict what was reported in [37]: there, substantial regions of parameter space did not show much impact of s δ , and the situation studied was an ITG scenario with adiabatic electrons.
The linear as well as nonlinear results obtained with these new equilibria and reported next, when matching the experimental δ, differ by typically multiple 10%.This supports the notion that the earlier simulations using numerical equilibria (relative to which the Millerized geometry can be seen as a perturbation) are near threshold and sensitive to moderate changes in input parameters.Aside from this consideration and the fact that the Millerized geometry and associated simulations are producing results relevant to the TCV experiments, the following analyses focus on δ-dependent instability and saturation physics rather than constituting quantitative predictions of immediately implementable scenarios.

Instability and transport
Per figure 7, similarly to the numerical-equilibrium situation, the dominant modes for the PT parameters are TEMs, with much weaker dependence on k x when δ < 0, whereas the NT scenario is dominated by ITGs across all δ, with weak k x dependence throughout and dominant k x ̸ = 0 branches at δ > 0.Here and hereafter, the PT and NT labels, still respectively encoded by red and blue in figures, refer to the profiles and Miller parameters of the two discharges other than δ, such that e.g. the PT curve can extend from negative to positive δ.
Eigenvalue calculations (not shown here) at k x = 0 for the PT profiles produce, across δ, multiple subdominant ITG branches at substantial fractions of the growth rate of the dominant TEM, with the relative importance of these ITG branches increased at δ < 0. At low k y or high k x , ITG can become dominant.Based on variations in driving gradients, the TEM is largely independent of both ω Te and ω Ti (implying a primary density-gradient drive) except for substantially negative δ, where ω Te becomes a driver.This property is partly due to the similar magnitudes of the density and electron temperature gradients at this radius.The mixed nature of the TEM drive is confirmed by crossphases with respect to the electrostatic potential Φ (not shown; similar analyses have been presented e.g. in [38]); these show the majority of the drive involves the The chosen ky = 0.8 is representative of the linearly and nonlinearly dominant physical scales.For PT profiles (red), TEM (frequency ω < 0, crosses) dominates, whereas ITG is the fastest-growing mode for NT profiles (blue).The mapping of linestyles to kx is solid for 0, dotted for 0.1, dashed for 0.2, and dash-dotted for 0.4.Dotted vertical lines indicate the nominal δ of the PT/NT equilibria.perpendicular temperature fluctuation T ⊥e,trapped of the trapped electrons, pointing at ∇B rather than curvature drive and highlighting mode properties characteristic of ∇T-TEM.In particular, these phases rule out the universal instability [39].
For NT parameters, ITG is the linearly dominant instability at all δ and k x ; positive δ exhibits lower growth rates, with dominant branches at k x ̸ = 0 as opposed to the dominant k x = 0 ITG at δ < 0. The cause for these properties is the same as discussed in section 3, relating to curvature and polarization (and, by extension, the local shear).As expected, the ITGs are sensitive (not shown) to ω Ti but not in immediate proximity to the critical gradient, in keeping with the interpretation that the low flux for the NT profiles in figure 6 is the result of proximity to the nonlinear critical gradient and Dimits regime.
The overall lower growth rates for the NT profiles, again consistent with nominal geometry, arise primarily from the slightly reduced TEM gradient drive.Note, however, that ITG modes and TEMs do not react equally at different δ and different k x .Furthermore, the present work does not take into account the potential impact of fast ions, and instabilities driven by such species can be substantially affected by triangularity, as well [40].
A distinct view on TEM activity is offered through the lens of the TEM available energy E A , which describes the energy available to TEM instability and turbulence in a given magnetic equilibrium and profile [41,42].E A does not contain any dynamical information, but has demonstrated some predictive capability with respect to growth-rate and flux scalings across different (tokamak and stellarator) configurations.In figure 8, the available energy, in arbitrary units, is compared with linear growth rates γ at k x = 0 as δ is varied for the PT discharge.Clearly, γ and E A exhibit mostly opposing trends.However, this analysis does not take into account possible changes in turbulent correlations lengths that can strongly affect the efficacy of E A .It is noted here that correlation lengths, as calculated from nonlinear simulations, increase substantially from strongly negative to strongly PT, and that a nonlinear correction to E A scales as the squared correlation length [43], implying that the trend in E A can be reversed and matched qualitatively with the growth rates.However, such an analysis is beyond the scope of the present effort, and the results presented here primarily highlight the fact that an available-energy calculation not accounting for nonlinear changes can fail to predict linear growth rates.
Considering E A is formulated for collisionless TEMs (although a collisional extension may be possible [44]), simulations are repeated without collisions, then additionally without electromagnetic effects, squareness, and ω Ti , all of which are neglected in the E A model, and lastly the ion temperature is set to match the electron temperature in parallel with these other modifications (note that ions are treated adiabatically in the available-energy derivation).While omitting collisions results in a substantial boost of growth rates, the other modifications cause only minimal to moderate changes, and most importantly, none of them reverse the qualitative trends.
Turning to properties of turbulence, nonlinear simulations in the Millerized geometries are performed for different δ.As figure 9 shows, for the TEM-dominated PT profiles, δ has little impact for moderate triangularity, with a slight overall trend opposing that of the linear growth rate.The inability of δ to affect fluxes at moderate δ may be a consequence [15] of the low κ = 1.12.However, a substantial reduction in electron electrostatic flux Q es e is seen for both signs of δ as its magnitude exceeds 0.5.The stronger reduction of Q es e at δ < 0 is partially offset by an increase in ion flux Q es i , consistent with the more mixed TEM-ITG state expected here from the linear analysis.Electromagnetic flutter transport (not shown) is negligible.The situation differs for the ITG-dominated NT profiles, which have lower fluxes overall but are not as marginal as in the case of the numerical equilibrium.Regardless, when increasing the ion temperature gradient by one third, Q es i and Q es e respectively increase by 260% and 180%, consistent with close proximity to the ITG turbulence threshold and confirming strong sensitivity to the driving gradient.Returning to the nominal gradients, Q es e and Q es i are comparable in magnitude and not very sensitive to δ, except at strongly negative triangularity, where an increase in ion flux is observed, which is not expected based on γ alone-as will be discussed later in section 4.2, linear physics can fully account for this feature.Note that the ITG trends in both the PT and NT cases are in apparent opposition to those reported in [17]; there, however, a different, less realistic treatment of the magnetic geometry was used, changing which their trends reverse [45] and match more closely with those reported here for δ < 0.
Flux spectra as presented in figure 10 reveal more about how the turbulence changes with δ; note the different color code in this figure.In both cases, albeit at different relative magnitudes, two overlapping flux structures are seen: a broader higher-k y peak dominated by electron flux Q es e that is attributable to TEM turbulence and a narrower low-k y peak dominated by ion flux Q es i arising from ITG turbulence.Changing δ, while slightly adjusting the TEM peak location, primarily affects the balance between ITG and TEM in accordance with the linear drive, to the point where ion and electron fluxes are comparable at δ = 0.6 in the NT case.The relatively higher k y associated with the TEM peak is expected based on the ratio of the temperature and density gradients [46].Not shown here is the parallel envelope of the flux, which at negative δ attains more structure, suggesting a more slablike nature of the ITG.The fact that the ITG flux for NT profiles begins increasing as δ is lowered below zero may be the result of the slab branch overtaking the higher-k y toroidal ITG.As a consequence, it is concluded that optimization based on triangularity may need to account for both slab and toroidal ITG.
One of the open questions regarding how triangularity affects profile evolution is whether δ impacts profile stiffness.For the present purpose, the term stiffness refers to the slope of the flux over the total gradient ω T + ω n summed over ions and electrons, rather than over each individual gradient separately.This definition is used here as reduction of the individual gradients invariably results in turbulence regime transitions, obfuscating any individual threshold gradients.The focus here lies on this stiffness near the turbulence threshold.
Stiffness is large when small changes in a normalized gradient result in large changes in flux.Focusing on the PT profiles, figure 11 shows how fluxes at three representative δ vary as ω T and ω n are lowered in synchronicity from their nominal values.Comparing δ = 0 and δ > 0, PT appears to reduce stiffness, whereas the picture for negative triangularity is more complex: fluxes are stiffer than at other δ nearer the nominal gradients, but as the electron flux approaches the ion flux (and thus the regime transition between TEM and ITG is neared), the stiffness is markedly reduced.Ion profile stiffness follows different and similarly nonlinear trends.Note that these findings expand on the picture presented in [11], where no difference in stiffness was observed between positive and negative δ. Stiffness is sensitive to the local situation, and may change substantially as e.g.criticality or a regime transition is approached.

Quasilinear fluxes
One method commonly employed in the context of transport modeling and turbulence optimization is the quasilinear, or mixing-length, model (see, e.g.[47], where nonlinear TEM fluxes were reproduced through quasilinear modeling).Here, the approach described in [48] is used, while accounting only for the most unstable eigenmode and retaining a single k x connection on either side of k x = 0 in the evaluation of the perpendicular wavenumber k ⊥ , resulting, for species j, in the quasilinear heat flux Here, k refers to wavenumbers k y , while w jk = Q lin j /Φ 2 is the quasilinear weight, which accounts for the Φ-T j cross phase, and θ is the ballooning angle.
The resultant fluxes, fixed to their nonlinear counterparts at δ = 0 by choosing an appropriate model constant C k , are contrasted with those from nonlinear simulations in figure 12, for a scan over δ.While for the NT profiles, the ITG-based dominant fluxes are captured well, the trend in the nonlinear electron fluxes (due to the dominant TEM) is not reproduced by the quasilinear model.Given the changes in relative growth rate of the k x ̸ = 0 modes compared with the dominant k x , another Figure 10.Spectra of electron (solid lines) and ion (dashed lines) electrostatic heat fluxes as functions of toroidal wavenumber ky for the PT (left) and NT (right) profiles.Here, the colors correspond not to the profiles but to three individual δ values: blue to −0.6, black to 0, and red to 0.6.Aside from the formation of a low-ky peak primarily in the ion flux for δ < 0, a slight shift of the bulk towards higher wavenumber is seen as δ is varied from −0.6 to 0.6.model is shown in the same figure that includes the finite-k x modes in the wavenumber sum in equation ( 1); evidently, this does not alter the excellent performance of the model for the NT profiles, whereas it does somewhat improve the predictions for the PT profiles, at least at positive δ.Thus, quasilinear modeling at substantial triangularity requires accounting for finite-k x contributions.The model still neglects subdominant eigenmodes; however, based on their properties as discussed in this section, they are unlikely to affect the outcome, as the only subdominant modes seen in this scenario are of ITG type.
Adding to the utility of quasilinear models beyond fast prediction is the fact that in comparing their output with that of nonlinear simulations, the physical origin of observed trends becomes apparent.Here, it is evident that, for the PT profiles, while the more abrupt stabilization at the most negative δ can be traced to linear stabilization, the nonlinear flux reduction at positive δ must result from a different mechanism; as will be Figure 12.Nonlinear (solid lines) and quasilinear electron and ion heat fluxes for the PT (red) and NT (blue) profiles, respectively.The model marked by dashed lines corresponds to that in equation ( 1), whereas the dotted lines are obtained by expanding the meaning of the k sum in the model to include radial wavenumbers kx.This latter approach achieves better but still not predictive performance for the PT profiles.discussed in section 4.4, zonal flows are the cause.Similarly, for the NT profiles, the saturation physics appears unaltered throughout the δ range, as linear physics can account for the nonlinear behavior in its entirety.

TEM fluctuation localization
To determine how the changes in geometry resulting from sign variations in δ affect instability and turbulence, figure 13 shows key geometric quantities, linear eigenfunctions (the dominant mode is ITG at k y = 0.1 and TEM at k y = 0.8), and turbulent fluctuations as functions of the parallel coordinate z; here and hereafter, non-zonal fluctuations are taken to be fluctuations with k y > 0. The squared perpendicular wavenumber k 2 ⊥ , a measure of the polarization, increases away from the outboard midplane z = 0.As described in [17], peaks in polarization can confine the ITG eigenfunction, as seen here for k y = 0.1, where Φ(δ > 0) is substantially narrower than at δ < 0. The curvature components K x and K y -describing the toroidal and radial variation of the background magnetic field, respectively, in combination with the parallel variation [49]are divided by the background magnetic field, corresponding to the secular and periodic contributions to the curvature drift [50], respectively.The latter shows somewhat larger regions of unfavorable curvature particularly for δ < 0, but this is not reflected in the k y = 0.8 TEM eigenfunctions.Rather, the latter appear to follow the k 2 ⊥ physics closer towards the outboard, where δ > 0 has lower values (up to |z| ≲ π /2).The nonlinear picture may be complicated by the fact that at δ = −0.6, a nonnegligible contribution to Φ may arise from ITG-like fluctuations; however, based on the fluxes shown earlier in figure 9, it may be concluded that nonlinearly, the local shear allows fluctuations to broaden for δ < 0. Note that the zonal Φ(k y = 0) was removed from the nonlinear data for the purpose of this comparison.

Zonal-flow saturation
Zonal flows [51] are fluctuations of the electrostatic potential Φ that vary only across flux surfaces, thus producing radial electric fields and, in consequence, a shear flow.They may saturate linear instabilities, most commonly ITG and ∇n-TEM [46]; when plasmas are driven by instabilities and dynamics are not dominated by collisions, this process typically occurs by catalysis of energy flow from unstable to stable eigenmodes [52][53][54], rather than by an inertial forward cascade.The analyses of zonal-flow behavior presented in this section pertain exclusively to the PT profiles.
A common but imperfect-and not necessarily wellconditioned-measure of zonal-flow strength in saturation is the shearing rate ω s = ⟨d 2 Φ /dx 2 ⟩, where the average is over the flux surface.Alternatively, one may measure the zonal Φ(k y = 0) and relate its magnitude to the non-zonal Φ(k y > 0), averaged over all finite k y .Both resultant quantities are shown in figure 14 as functions of δ; their opposite scaling with δ is explained readily by a reduction of the characteristic k x of the Figure 14.Shearing rate ωs (dashed line) and zonal-to-non-zonal ratio Φ Z /Φ NZ (solid line) versus triangularity, averaged over the quasi-stationary state.Here, Φ Z and Φ NZ are obtained by respectively averaging over all ky = 0 and ky > 0 fluctuations in a given turbulence simulation.The former decrases while the latter increases with δ, pointing at a substantial shift in characteristic zonal-flow kx.zonal flow with triangularity, i.e. a spectral down-shift as δ increases.This is directly related to the corresponding aforementioned increase in correlation lengths: as for δ < 0 the linear drive of modes with k x ̸ = 0 is very strong (and potentially even dominant), it is expected that turbulent fluctuations (e.g. in Φ(k y > 0)) with k x ̸ = 0 grow to larger amplitudes here than in the δ > 0 simulations.In turn, this causes zonal flows to be driven more strongly at smaller radial scales, as energy transfer is typically local in k-space [55,56].
A priori, it is unknown whether nonlinear fluxes should follow the (inverse) trend of ω s or Φ Z /Φ NZ .The earlier comparison of quasilinear and nonlinear fluxes in figure 12 appears to require (presuming they constitute the primary saturation mechanism in this mixed turbulence regime) that zonal flows provide stronger flux suppression at positive δ, consistent with the Φ Z /Φ NZ curve.Obeying this ratio (or, alternatively the scaling of k x Φ) rather than the shearing rate is in keeping with the stable-mode-saturation paradigm.The question why the zonal flow exhibits different strength and scale for positive versus negative δ, however, requires further analysis.Zonal flows are typically seen as in a state of balance between nonlinear drive and linear damping, both of which are quantified next.
The strength of the nonlinear drive of zonal flows, measure by growth rates of secondary instability [57], differs only mildly as δ is varied, see figure 15.Given that this standard secondary-instability analysis does not retain any k x connections of the streamer or sideband (compare [58]), the secondary γ at the higher zonal-flow k x is larger by a factor of four as expected; note that the two k x shown here are representative of the characteristic zonal-flow k x at low and high δ.As the secondary-instability analysis depends on the structure of the frozen streamer, any changes here are the result of the different eigenfunction widths for positive and negative δ.The secondary growth rates relatively minimally favor larger-scale zonal flows at higher triangularity, but this difference in secondary-instability drive is too weak to substantially affect the outcome when compared with the difference in primary-instability drive between k x .More importantly, the δ > 0 zonal-flow drive is slightly weaker, opposing the trend in zonal-flow strength observed in the quasi-stationary state.Thus, changes with δ in the ability of the zonal flow to maintain itself against linear damping are the likely candidate to explain the nonlinear simulations.
Figure 16 reports the residuals Φ res of zonal flows after initial collisionless decay [59].To this end, linear simulations at k y = 0 are initialized with a finite-amplitude zonal flow and evolved until the zonal flow has decayed to its residual, see e.g.[60].For these simulations, the collision operator was removed to ensure stable residuals; note that any collisional decay of the residual affects the system equally at each δ, as ν ei is held constant throughout the nonlinear scan.Clearly, little difference can be found between different k x , and thus the scale disparity implied by figure 14 cannot be the result of differential ability of the zonal flows to maintain a residual.Relatedly, A theoretical prediction is shown as a dotted line.Little variation among kx is seen, while residuals are boosted as δ is increased from negative to positive.the frequency of the geodesic acoustic mode (not shown) as extracted from the same simulations during the zonal-flow decay is largely independent of δ.More importantly, however, Φ res increases by a factor of nearly 2 between δ = −0.8 and 0.8.This sufficiently affects the zonal-flow balance to account for the differences between quasilinear and nonlinear fluxes in figure 9.
A theory of zonal-flow residuals that captures various geometric effects has been proposed in [61].Comparison with the dotted fit curve shows good agreement for the present scenario, where due to specific limitations of the theory, non-δ Miller parameters have been absorbed into a scalar free parameter, and the fit results in f nd = 3.0 for where S is the shaping function.Another theoretical approach to zonal-flow residuals at finite triangularity [62] finds that the residual for δ < 0 is lower due to an expanded particle trapping region in negative triangularity, although the theory requires fixed eddy size, which is not the case here.Their approximate expression for the residual, which does not account for all Miller parameters used here, is qualitatively consistent with the present data, although the δ effect is stronger in the present TCV case than in their Gene simulations.While the dependence of the characteristic k x of the linear drive on δ can account for spectral changes in the zonal flow, the increased zonal-flow amplitude at positive δ is the result of a larger residual Φ res .It may be possible to include this effect in a predictive manner into a quasilinear transport model to improve its performance in recovering nonlinear fluxes.One path is provided by the inclusion of the triplet correlation time [53], where linear complex frequencies are combined to provide a measure of saturation efficiency.By incorporating zonal-flow damping in this quantity, the resonance between unstable and stable mode will be broken more strongly at δ < 0 (where Φ res is lower), and resultant quasilinear fluxes will increase at this triangularity.Note, however, that in a collisional regime such as the present one, collisional damping may be too strong for a quasilinear model based on the triplet correlation time to be quantitatively applicable.Future study will need to determine at what ν ei such a model will become inapplicable.

Conclusions
TCV discharges with moderately large triangularity δ of either negative or positive sign in the outer core were examined by means of linear and nonlinear gyrokinetic simulations.Due to close proximity to critical gradients and resultant substantial profile stiffness, simulated fluxes are highly sensitive to moderate gradient variations, making a quantitative comparison of transport difficult; note, however, that such moderate variations in gradients produce fluxes not inconsistent with experiment.Converting the numerical magnetic equilibrium into a parameterized form resulted in moderate quantitative changes and allowed for extrapolation to larger magnitudes of negative and positive δ.While the ITG-dominated NT discharge exhibits only moderate dependence of growth rates on δ and nonlinear flux scalings captured well by quasilinear modeling, the TEM-dominated PT discharge offers a more complicated picture.There, irrespective of collisionality, growth rates at δ < 0 are low despite a large reservoir of available energy, whereas at δ > 0, nonlinear fluxes are reduced below quasilinear estimates.Although some of this nonlinear-quasilinear discrepancy can be eliminated by accounting for instability at finite k x when composing quasilinear fluxes, the considerable remainder is the result of much stronger zonal flows at δ > 0. The latter are driven at nearly equal rates across δ; at the same time, the zonal-flow residuals, measuring how well the system can maintain zonal flows once produced, rise substantially as δ is increased from negative to positive, explaining the lower nonlinear (compared to quasilinear) fluxes there.Finally, the profile stiffness was evaluated for the TEM-dominated PT scenario at different δ, and it was found that, possibly due to proximity to boundaries between ITG, ∇T-TEM, and ∇n-TEM regimes, stiffness changes between and even within triangularities, depending on the specific gradient settings.
Overall, the analyses and results presented here show that substantial benefit in terms of flux reduction can arise from extreme flux-surface triangularity of either sign, notwithstanding potentially prohibitive other considerations in achieving such geometries.Furthermore, whether negative or positive δ performs better in terms of instability and transport depends on not only the instability regime but also multiple other aspects of a given scenario.Finally, the ability to quantitatively predict transport based on quasilinear models may hinge on taking into account zonal-flow residuals when comparing configurations across large ranges of δ.
In future investigations, where the benefit of extreme triangularity is to be assessed under reactor-relevant conditions (where heating will be from fusion alphas rather than Ohmic and discharges will be diverted, while the gyroradius scale will be more strongly separated from the device scale), ITG turbulence will be the focus, with possible modifications due to the presence of fast ions.The orbits of such fast ions are critically impacted by triangularity [63], which may impact turbulence and potentially H-mode access.Such a scenario will likely require larger elongation [16], in addition to careful consideration of triangularity profile effects.

Figure 1 .
Figure 1.Cross-sections at constant toroidal angle of flux surfaces ρtor = 0.1 − 0.9 for the NT discharge (blue) and the PT discharge (red).The radial position 0.7 chosen for subsequent analyses is marked as a slightly bolder line.

Figure 2 .
Figure 2.Electron temperature (solid) and density (dashed) profiles for the NT discharge (blue) and the PT discharge (red).The ion temperature is assumed to be T i0 = 0.133T e0 across the domain, implying identical normalized gradients ω T (dotted) and ωn (dash-dotted) for both species.The location ρtor = 0.7 is marked as a vertical black line.

Figure 4 .
Figure 4. Heat-flux profiles based on ASTRA calculations.Red lines correspond to the PT discharge, blue lines to the NT discharge.Shown are the electron heat flux as a solid and the ion heat flux as a dashed line.Error bars indicate variability within a 10 ms time interval.Given the slightly higher core density for NT in figure2, the NT confinement is slightly better but similar to the PT confinement.

Figure 5 .
Figure5.Linear growth-rate spectra in realistic geometry.The linestyle encodes at which kx modes are centered: solid for 0, dotted for 0.1, dashed for 0.2, and dash-dotted for 0.4.Lines with/without crosses indicate negative/positive frequencies, corresponding to TEM/ITG.The PT discharge (red) has kx = 0 dominate across ky, whereas the NT discharge (blue) has near-identical contributions from all kx.A flipped scenario (pink) with NT geometry and otherwise PT parameters approaches the TEM-ITG transition and has reduced kx dependence; further changing this flipped scenario by using the collision settings from the NT discharge (black) leads to γ close to that of the NT discharge.

Figure 6 .
Figure 6.Electrostatic heat-flux powers as functions of time for realistic-geometry NT (blue) and PT (red).Solid lines correspond to electron, dotted lines to ion fluxes; note that the ions and electrons produce very similar saturated fluxes in the NT case.Electromagnetic fluxes (not shown) are at much lower levels.

Figure 7 .
Figure 7. Linear growth rates in Miller geometry as functions of δ.The chosen ky = 0.8 is representative of the linearly and nonlinearly dominant physical scales.For PT profiles (red), TEM (frequency ω < 0, crosses) dominates, whereas ITG is the fastest-growing mode for NT profiles (blue).The mapping of linestyles to kx is solid for 0, dotted for 0.1, dashed for 0.2, and dash-dotted for 0.4.Dotted vertical lines indicate the nominal δ of the PT/NT equilibria.

Figure 8 .
Figure 8. Collisional (black solid line) and collisionless (red solid line) growth rates at ky = 0.8, for Miller geometry, neither of which match the qualitative behavior of the available energy E A (black dashed line, arbitrary units).The trend in E A holds over a wide range of elongations κ (not shown).Included for comparison are simulations without various other effects relating to geometry, finite β, and ion dynamics (blue and pink solid curves), in keeping with assumptions made in deriving E A .

Figure 9 .
Figure 9. Electron (solid lines) and ion (dashed lines) electrostatic heat fluxes for the PT (red diamonds) and NT (blue squares) profiles as functions of δ.Dotted vertical lines indicate the nominal triangularity for each scenario.The TEM turbulence in the PT case is reduced at large |δ|, whereas the ITG fluxes in the NT case increase for strongly negative δ.

Figure 11 .
Figure 11.Impact of (combined density and temperature) gradient reduction on electron (solid lines) and ion (dotted lines) fluxes for three different δ = −0.6 (blue), 0 (black) and 0.6 (red) for the PT profiles; at a multiplier of 1, the nominal fluxes are reproduced.Stiffness as obtained from the slope of these curves can change within and between triangularities.

Figure 15 .
Figure 15.Growth rates of the secondary instability producing the zonal flow in the PT scenario.Two zonal-flow wavenumbers are considered, kx = 0.058 (black diamonds) and kx = 0.23 (red squares).Only a small change is seen with δ.

Figure 16 .
Figure 16.Collisionless zonal-flow residuals for different radial wavenumbers kx (colors) as functions of δ for the PT geometry.A theoretical prediction is shown as a dotted line.Little variation among kx is seen, while residuals are boosted as δ is increased from negative to positive.

Table 2 .
Numerical input parameters for which gyrokinetic simulations are converged for NT and PT numerical and parametric equilibria.