Available energy of trapped electrons in Miller tokamak equilibria

Available energy (Æ), which quantifies the maximum amount of thermal energy that may be liberated and converted into instabilities and turbulence, has shown to be a useful metric for predicting saturated energy fluxes in trapped-electron-mode-driven turbulence. Here, we calculate and investigate the Æ in the analytical tokamak equilibria introduced by Miller et al. (Phys. Plasmas, vol. 5, issue, 4, 1998, pp. 973–978). The Æ of trapped electrons reproduces various trends also observed in experiments; negative shear, increasing Shafranov shift, vertical elongation and negative triangularity can all be stabilising, as indicated by a reduction in Æ, although it is strongly dependent on the chosen equilibrium. Comparing Æ with saturated energy flux estimates from the TGLF (trapped gyro-Landau fluid) model, we find fairly good correspondence, showcasing that Æ can be useful to predict trends. We go on to investigate Æ and find that negative triangularity is especially beneficial in vertically elongated configurations with positive shear or low gradients. Furthermore, we extract a gradient-threshold-like quantity from Æ and find that it behaves similarly to gyrokinetic gradient thresholds: it tends to increase linearly with magnetic shear, and negative triangularity leads to an especially high threshold. We next optimise the device geometry for minimal Æ and find that the optimum is strongly dependent on equilibrium parameters, for example, magnetic shear or pressure gradient. Investigating the competing effects of increasing the density gradient, the pressure gradient, and decreasing the shear, we find regimes that have steep gradients yet low Æ, and that such a regime is inaccessible in negative-triangularity tokamaks.


Introduction
Energy transport in tokamaks and stellarators is largely dominated by turbulent energy losses, which severely degrade the energy confinement in these devices.A detailed understanding of how various parameters characterising the plasma and the magnetic field geometry, such as magnetic shear and the pressure gradient, affect the turbulent transport properties would be helpful in comprehending and mitigating this.The standard method to assess the turbulence properties of any given tokamak is to perform nonlinear gyrokinetic simulations.However, such simulations are computationally expensive because of the very disparate time-and length scales characterising the turbulence and the transport.Thus, it would be beneficial to find a reduced model capable of predicting the level of turbulent transport by simpler means.
In a recent publication, it was shown that the available energy (AE) of trapped electrons can serve as such a reduced model (Mackenbach et al. 2022), at least for turbulence driven by the plasma density gradient.Any plasma possesses a maximum amount of thermal energy that can be converted into instabilities and turbulence (Gardner 1963).
This "available" energy can be calculated by performing a Gardner restacking of the plasma distribution function f , in which phase-space volume elements are rearranged in a manner that respects Liouville's theorem (Kolmes et al. 2020;Kolmes & Fisch 2020).The restacking of f that minimises the thermal energy results in a "ground state" distribution function f g , and the AE is defined as the difference in thermal energy between f and f g .If one imposes the additional constraint that adiabatic invariants be conserved in the restacking process, the AE becomes relevant to magnetically confined plasmas (Helander 2017(Helander , 2020)).In fusion plasmas, the magnetic moment µ is generally conserved for all species, and the parallel adiabatic invariant J = mv ∥ dℓ is conserved for magnetically trapped electrons.
A significant portion of the electrons are trapped and can contribute to turbulence through trapped electron modes (TEMs).The AE of trapped electrons correlates with the turbulent energy flux for such TEM-driven turbulence over several orders of magnitude in saturated energy fluxes (Mackenbach et al. 2022).This correlation is expressible as a simple power law, where the saturated energy flux, Q sat , was found to be related to the available energy, which we denote by A in formulas, via approximately (1.1) This relation was found to hold for both a tokamak and stellarators, and for various values of the density gradient.Aside from this relationship, other links have been found by Kolmes & Fisch (2022) where quasi-linear plateauing is shown to be related to a concept closely connected to AE, highlighting other links to transport physics.In any case, to gain a deeper understanding, it is of interest to derive an explicit expression of AE in tokamak geometry, in order to investigate the dependence of it on various geometrical and plasma parameters.This is our aim in the present paper, where we compute AE for the family of tokamak equilibria constructed by Miller et al. (1998).The starting point is the following explicit expression for AE in a flux tube of any omnigenous equilbrium (Helander 2020;Mackenbach et al. 2023a), including that of a tokamak, Here, L is the total length of a field-line completing one poloidal turn, B 0 is some reference magnetic field strength, z = H/T 0 is the particle energy normalised by the temperature, λ = µB 0 /H is the pitch angle, and ∆ψ t and ∆α C denote the size of the flux-tube in the radial and binormal directions respectively (we have parameterised the radial coordinate by means of the toroidal flux ψ t and the binormal by means of the Clebsch angle α C ).We furthermore sum over all magnetic wells with a certain value λ.The hatted quantities in the integrand denote normalised frequencies, with ωα being the normalised bounce-averaged drift precession frequency, ωT * the normalised electron diamagnetic drift frequency, and ĝ1/2 the normalised bounce time.They are explicitly defined as where we have denoted the ratio between the gradients by η = (d ln T /dψ t )/(d ln n/dψ t ).
Finally, R[x] = (x + |x|)/2 is the ramp function.Using the above expressions, we shall find the AE of trapped electrons in any Miller tokamak.

Theory
2.1.The available energy in any omnigenous system We first note that the integral over z can be rewritten into a convenient form.We define two functions that are independent of z, namely, (2.1) With these functions, the integral over the normalised energy z reduces to the following form; (2.2) This integral can be solved analytically, and its functional form depends on the signs of c 0 and c 1 , resulting in four different conditions.The easiest case to evaluate is the case where c 0 < 0 and c 1 > 0. In this case, the argument of the ramp function is always negative, and hence the integral reduces to zero.The second case is when the argument of the ramp function is always positive, which occurs whenever c 0 ⩾ 0 and c 1 ⩽ 0. The integral then reduces to the following form, (

2.3)
There are two cases left to consider.First, we inspect the case where the argument of the ramp function is positive for low z but becomes negative for high z, that is, c 0 ⩾ 0 and c 1 > 0. The unique point where the argument of the ramp function vanishes is the following, (2.4) Thus, the integral becomes (2.5) This integral can be expressed in terms of the error function, erf(x) = 2/ √ π x 0 exp(−t 2 )dt, The final case is that where the argument of the ramp function is negative for low z but becomes positive for high z, that is, c 0 < 0 and c 1 ⩽ 0. The integral then becomes Note that I z ⩾ 0, ∀(c 0 , c 1 ) ∈ R 2 , which can also be seen in Fig. 1.The AE can now be found by executing the integral over the remaining coordinate Note that this expression is completely general; no approximations have been made in executing these integrals, aside from the preceding assumption of omnigeneity.It is also interesting to note that from this expression one can see that there are no tokamak configurations with vanishing AE, at least in leading order near the axis.This conclusion can most readily be drawn by investigating the expression for ω α from Connor et al. (1983).Here, one can find that there is always a zero crossing for ω α (with no pressure gradient), which implies that c 0 and c 1 must change sign.As such, the available energy must be non-zero (as either I(c 0 , c 1 ) or I(−c 0 , −c 1 ) must be non-zero).Formally, this corresponds to the fact that such a zero crossing implies that the device does not have the so-called maximum-J property, which is required for the linear stability of trapped electron modes (Proll et al. 2012).† To make further progress in solving Eq. (2.8), one requires the function ωα (λ), which in turn requires a specification of the equilibrium.In this paper, we will use the local construction of the equilibrium, employing a formalism developed by C. Mercier & N. Luc (1974).

Construction of local equilibria
Equilibria are constructed by finding a radially local solution to the Grad-Shafranov equation, and this solution allows us to find ωα .We highlight the essential components of this derivation, which essentially follows the steps taken by Miller et al. (1998), and a thorough overview is given by Candy (2009) The Mercier-Luc formalism requires the shape of the flux surface, the poloidal field B p on that flux surface, the gradients of the pressure p(ψ), and the toroidal field function f (ψ) = RB ϕ on the flux surface, where R is the major radial coordinate, B ϕ is the toroidal component of the magnetic field, and ψ is the poloidal flux.We parameterize the flux surface as R s = R s (l) and Z s = Z s (l), where l measures the poloidal arclength along the flux surface.It is also useful to define a tangential angle u, which measures the angle between the unit vector in the major radial direction e R and the vector tangential to the flux surface e l clockwise, thus (2.9b) † This correspondence between the maximum-J property and AE is shown in (Helander 2017), and can also be understood from Eq. (2.8).A device is said to be maximum-J if ∂ ψ J < 0 for all particles, which implies ωα > 0 for all λ.For η < 2/3, Eq. (2.1) implies that c0 < 0 for a radially decreasing density profile, and c1 > 0, thus the integrand of the AE reduces to zero since Iz = 0.
With this definition, the angle u can be calculated by du/dl = −1/R c , where R c (l) is the radius of curvature of the poloidal cross section, and the negative sign arises because the poloidal arclength is measured clockwise.We go on to introduce a radial-like expansion variable ρ which is zero on the given flux surface, in terms of which the cylindrical coordinates become (2.10b) The metric tensor in these coordinates has non-zero components only on the diagonal (which is to be expected as we ensured orthogonality in the construction), where we use the convention and substitute into the Grad-Shafranov equation, which in leading order reduces to This allows one to find the radial variation of the poloidal magnetic field by using (Helander & Sigmar 2005) From this equation we can immediately see that ψ 1 /R s = B p,s , with B p,s being the poloidal field on the flux-surface as indicated by the subscript.As such, the poloidal field strength can be written as The toroidal field is found from its definition where B ϕ,s = f (ψ 0 )/R s .The total magnetic field strength is also readily derived (2.18) Note that the derivatives ∂ ρ b p , ∂ ρ b ϕ , and ∂ ρ b are given in square brackets.The radial variation of the poloidal line element is readily found from the metric tensor, In these equations f ′ (ψ 0 ) is treated as a free parameter, but it is difficult to ascertain if the chosen value of this parameter is realistic.It is more convenient, however, to specify the magnetic shear, which is related to f ′ (ψ 0 ).This can be made explicit by investigating the safety factor (2.20) Taking the derivative of the safety factor with respect to ψ, one finds an equation describing this relationship, (2.21) We also wish to relate the arclength along a magnetic field line to the poloidal arclength.These quantities are related as Finally, the poloidal coordinate can be expressed in terms of the poloidal angle θ rather than the poloidal arclength by and the total arclength thus becomes (2.24)

Non-dimensionalisation and available energy
We proceed to make the various functions dimensionless as in Roach et al. (1995), and in doing so we will introduce various dimensionless constants which will be useful for the remainder of the analysis.We assume that we have been given the dependencies of the various functions in terms of the minor radial coordinate r, which in turn relates to the major radial coordinate R 0 through the inverse aspect ratio of the flux surface in question ϵ = r/R 0 .Furthermore, we define our reference field B 0 through the relation f (ψ 0 ) = B 0 R 0 .Let us now define various dimensionless functions of interest, (2.25f ) One also needs to relate ψ to r, which can be done by investigating the poloidal field as in Eq. (2.14)We proceed to define a dimensionless pressure gradient, analogous to the α parameter used in s-α geometry, (2.30) Note that this dimensionless pressure gradient is not related to the Clebsch angle.The pressure gradient can in turn be used to define a dimensionless toroidal current density (2.31) We go on to define the shear s in the following manner which can be substituted into Eq.(2.21) to relate the shear to f ′ R 0 as where we have defined the geometric constants C 1 to C 4 as dθ, (2.34b) dθ.
(2.34d ) The radial derivatives of the magnetic field become These expressions are the same as Roach et al. (1995, Eq. (14)), where the differences in sign arise because the sign convention for Rc is different here and ρ has the opposite sign.Finally, we express the total magnetic field length as Bp,s dθ. (2.36) We now turn our attention to the precession frequency, which we calculate from (1.3a).To simplify the calculation slightly, we note that the operator ∆ψ t ∂ ψt ≈ ∆ψ∂ ψ to leading order around smallness of the radial coordinate ρ, as we can approximate ∆ψ t ≈ ∆ψ∂ ψ ψ t .Using this identity, we find the same expression as in Roach et al. (1995), (2.37) where we define the bounce averaging operator in angular brackets as (2.38) We rewrite the precession frequency as Next, we investigate the Jacobian ĝ1/2 , which is the normalised bounce time, and find that it is equal to We rescale it with a factor ϵ 1/2 to acount for the fact that in smallness of ϵ the integrand of the bounce time goes as 1/ √ ϵ.Therefore, we define (2.41) The AE now becomes where the prefactor 1/ϵ to the integral deliberately not cancelled against the √ ϵ, so that the integral in brackets is to lowest order independent of ϵ, as the integration range scales as ϵ.With the above expression, we go on to define a dimensionless AE.We take steps in accordance with (Mackenbach et al. 2022), and calculate the fraction of the total thermal energy that is available.The thermal energy of a plasma in a flux tube can be calculated by expanding around ψ t = ψ t,0 and α C = α C,0 and retaining only the constant terms, resulting in We then define the available energy as a fraction of the thermal energy as (2.44) Simplifying the expression using ∆ψ = ∆r∂ r ψ, one finds that (2.45) ∆r measures the length-scale over which energy is available, i.e. a typical length-scale over which gradients can be flattened.We take this to be proportional to the correlation length, typically found to be the gyroradius.Therefore, let us set where ρ g is the gyroradius, and C r some function of order O(ρ 0 g ).This function is not known a priori, and may vary.For example, if there are large radial streamers present in the system C r may be significantly increased.The dimensionless AE now becomes (2.47) This expression has various scalings which are of interest.First, we see that reducing the aspect ratio for fixed ρ g /R 0 is beneficial since it leads to fewer trapped particles.Note that, in the limit of a large aspect ratio, the trapping fraction scales as √ ϵ, which is the same dependency found here.A reduction in the expansion parameter ρ g /R 0 (at fixed ϵ) is also found to help decrease AE.
As a final step, we introduce the dimensionless density gradient with which c 0 and c 1 reduce to an especially simple form (2.49) Importantly let us make note of the fact that the radial coordinate r may have different conventions.In previous investigations (Mackenbach et al. 2022(Mackenbach et al. , 2023a) ) the radial coordinated was defined via the square root of the toroidal flux . (2.51) In the aforementioned investigations ∆r eff was chosen as ∆r eff = ρ, resulting in a good correlation with turbulent energy fluxes.Therefore, we choose C r such that ∆r eff = ρ, which means that and we shall use this choice of C r from here on.The parameters (κ, δ) vary from the plot on the left to the plot on the right, as (2/3, 0.9), (2/3, −0.9), (3/2, 0.9), and (3/2, −0.9).All plots have R 0 = 1 and ϵ = 1/3.

Miller geometry
Finally, we choose our equilibrium to be of the type discussed in Miller et al. (1998).The key step is to parameterise the flux surface as a standard D-shaped tokamak in terms of the poloidal angle θ, (2.53b) Here, R 0 (r) is the centre of the flux surface, κ(r) is the elongation, and δ(r) is the triangularity.An important feature of this parameterisation is that it is up-down symmetric, which can be seen by invariance under (Z s , θ) → −(Z s , θ).The poloidal magnetic field can then be calculated by (2.26), and the equilibrium is fully specified by the following set of 9 parameters; [ϵ, κ, δ, s κ , s δ , ∂ r R 0 , q, s, α], where s κ = r∂ r ln κ and s δ = r∂ r arcsin(δ).Henceforth we shall refer to this set of numbers which determines the local geometry as a "Miller vector", (2.54) Cross-sections are plotted in Fig. 2, to serve as a reference for the various shapes mentioned in the following sections.We furthermore recognise that it is possible to calculate the toroidal flux enclosed by a poloidal cross-section, and one may retrieve analytical expressions by expanding it around the smallness of ϵ, (2.55) with J n (x) being the n th Bessel function of the first kind.In terms of an effective r, we then have (2.56) and we find that the factor C r becomes (2.57) For shaped equilibria (i.e.κ ̸ = 1, δ ̸ = 0, or ϵ → 1), C r will differ from unity and one should keep this important caveat in mind.
2.5.An analytical limit: large aspect ratio s-α tokamak We proceed to investigate a limiting case of Miller geometries; namely that of a large aspect ratio tokamak with circular flux surfaces and a steep local pressure gradient, which we shall henceforth refer to as the s-α limit, and this calculation is equivalent to analyses given in Connor et al. (1983); Roach et al. (1995).This will serve as a computationally efficient model in such geometries, and will furthermore be used as a benchmark for the more general calculation of the AE.The algebraic details of this derivation are given in Appendix A, and here we highlight the central steps.It is convenient to express λ as a trapping parameter k 2 , where the deeply trapped particles have k = 0 and the barely trapped particles have k = 1.This mapping is given by λ = 1+ϵ(1−2k 2 ), so the magnetic field may be written as (2.58) One can now perform the bounce-averaging integrals required for Eq.(2.37) in the s-α limit, resulting in where we define where K(k) and E(k) are complete elliptic integrals of the first and second kind, respectively.The normalised bounce time, as given in (2.41) is equal to (2.61) Finally, from Eq. (2.56) we see that C r = 1 in this limit.The AE now becomes a straightforward integral of known functions over k which can efficiently be computed numerically.

Numerical results
Two codes have been constructed: one that computes the integral of (2.62) using standard integration routines, and a numerical routine that computes both the precession frequencies and the AE as given in (2.44), both of which are computationally cheap (fractions of a CPU second per evaluation).First, we shall verify the relationship between AE and turbulent transport.Next, we shall investigate the results obtained for the s-α circular tokamak, after which we shall investigate how AE varies in Miller geometries as a function of various parameters.The code used to generate these results is freely available on GitHub †.The bounce-integrals required in Eq. (2.44) are evaluated using numerical methods detailed in Mackenbach et al. (2023b).Finally, we take the prefactor ρ g /R 0 to be unity in all plots presented below, so when converting to a real device, one should multiply the AE by a factor (ρ g /R 0 ) 2 .

Comparison with tglf
Our first course of action is comparing AE with turbulent energy-flux calculations in tokamak geometries, to verify its relation to turbulent transport in such geometries.At the moment, nonlinear gyrokinetic simulations are computationally too expensive for detailed parameters scans, and therefore we instead employ the quasi-linear tglf (trapped gyro-Landau fluid) code (Staebler et al. 2007;Staebler & Kinsey 2010).Some key differences between the two models are highlighted before any comparison is made.tglf computes the linear eigenmodes of a variety of instabilities, ion and electron temperature gradient (ITG, ETG) modes, electromagnetic kinetic ballooning (KB) modes, as well as trapped-ion and trapped-electron modes (TIM, TEM), and then applies a quasilinear saturation rule to accurately fit the fluxes from nonlinear gyrokinetic simulations.For quasi-neutrality purposes, tglf requires the inclusion of at least one ion species.These are fundamental differences to the formulation of the AE described in this work, which only accounts for the AE of trapped electrons.Therefore, when setting up tglf, care was taken to ensure the modelled turbulent energy-fluxes were as much as possible due to instabilities dominated by trapped electrons, using settings analogous to those used in recent gyrokinetic simulations in a similar regime (Proll et al. 2022).Given the lack of collisions in this regime, the expected dominant instabilities should be of the collisionless trapped-electron mode (CTEM) variety.However, some other instabilities can also arise from interactions with the ion population.Thus, to ensure that the dominant instabilities in the tglf simulations were as relevant as possible for our comparison, only contributions from modes propagating in the electron-diamagnetic direction were included, which excludes e.g. the ubiquitious mode (Coppi & Pegoraro 1977), which propagates in the ion direction.Furthermore, we find that, for the scenarios considered in this work, adding an equally large electron temperature gradient to the density gradient, i.e. taking η = 1, significantly decreased the amount of non-TEM modes dominant in tglf simulations, and as such we set η to unity for the comparison.The recent SAT2 (Staebler et al. 2021) quasilinear saturation rule for tglf was used, as it includes the impact of plasma shaping on the quasilinear saturation (Staebler et al. 2020).Although tglf also uses a Miller parameterisation of the local equilibrium, we note that it does not use the same normalisation as Roach et al. (1995) followed in this work, and care has been taken to convert between the two.We finally stress that the current model for C r is a fairly simple model, and that prediction can be refined using a more sophisticated model.This can, for example, be done by using some fitting function for C r , where one finds the best-fit parameters which minimise the error between the energy flux and the prediction of AE. † Install the code via https://github.com/RalfMackenbach/AE-MillerMiller vector [ϵ, κ, δ, s κ , s δ , ∂ r R 0 , q, s, α] of [1/3, κ, δ, 0, 0, 0, 2, 0, 0], whereas the middle column has [1/3, κ, δ, 0, 0, 0, 2, 1, 0], and the right-most column has as Miller vector of [1/3, 1, 1/2, 0, 0, 0, 2, s, α].The first column has ωn = 3, the second and third column have ωn = 6, and η = 1 in all plots.The white masked-out regions have a dominant instability which is not in the electron direction, and as such we filter them out.Finally, note that the colour bars are the same scale in each column.
For the comparison we use the gyro-Bohm normalised energy fluxes computed by tglf, where Q e is the electron energy flux from tglf, and Q GB is the gyro-Bohm energy flux.This is compared to the estimate of the gyro-Bohm normalised energy flux from AE Mackenbach et al. (2022Mackenbach et al. ( , 2023a)), which is where the constant of proportionality is taken from the fit presented there we was found C Q ≈ 1.0 • 10 3 .With such a power law, a linear correlation between Q A and Q e from nonlinear gyrokinetic simulations was found for pure density gradient-driven TEMs, which is different from the current comparison in which both the electron temperature and the density gradient drive the TEM (η = 1).The data points in the comparison are chosen in order to verify that tglf reproduces some trends that will be discussed in following sections.
A comparison in the (κ, δ) and (s, α) planes is displayed in Fig. 3.One can see that there is good correspondence in trends: decreasing the magnetic shear and/or increasing the pressure gradient helps in reducing the energy flux, as does increasing the elongation.However, there are also differences visible between the two models for the energy flux, which are evident in the (s, α)-plot.A clear discrepancy can be seen at high shear values (s ≈ 3), where the tglf energy flux drops and the AE estimate does not, and the AE furthermore overestimates the energy flux at high shear.In the (κ, δ)-plots the trends are well captured by AE, with some differences.To further investigate the relationship between the two estimates of the energy flux, all the simulation data shown in Fig. 3 have been combined in a scatter plot shown in Fig. 4. In order to check consistency with previous findings, we have furthermore included the data points of Mackenbach et al. (2022), which are nonlinear simulations in general geometries.Here, we see that there is a linear relationship for most of the data (we have added a red line with the expected linear relationship), although there exist data points that deviate more significantly from the linear relationship.There are various reasons why such a discrepancy may occur: • There may be other instabilities present (though not necessarily dominant) that are not captured by the AE of trapped electrons, such as the ubiquitous mode, or the universal instability (Landreman et al. 2015;Helander & Plunk 2015;Romanelli 1989;Costello et al. 2023).More generally, if there are instabilities present that do not derive their energy from trapped electrons, the current AE-model is no longer expected to be an accurate measure.
• The AE length-scale C r may vary more significantly for certain choices of equilibrium parameters and the current choice given in Eq. (2.57) may not be accurate.
• Recall that AE can be interpreted as an upper bound on the amount of energy that can be released.If the portion of the AE that resides in stabilising modes deviates markedly • The tglf's quasilinear saturated fluxes in both the (κ, δ) and (s, α) planes show occasional extreme outliers for small changes in input.tglf has been extensively verified against a wide variety of nonlinear gyrokinetic simulations (although further validation for negative triangularity is currently being pursued), but the regime explored in this work is not the typical input space and could require separate verification.
• Although not present in the current set of simulation data, the AE of trapped electrons will certainly cease to be an accurate model in situations where the trapped electrons play no role, such as in the case of a pure ion temperature gradient, and no gradients in of electron temperature or density.
We stress that the scatter does mean that predictions may be faulty if one lies within the scatter of the fit.However, seeing that general trends are well captured by AE, it may serve as a useful estimate for transport and trends at low computational cost (AE calculations are roughly a factor 50 faster than the presented tglf calculations).

s-α geometry
We now shift our attention to the behaviour of AE on the various free parameters found in tokamak equilibria.Recalling that we have derived two AE expressions, one for any Miller geometry and one for the large-aspect ratio limit, let us start by investigating the latter.A plot of the AE calculated from Eq. (2.62) is given in Fig. 5 as a function of magnetic shear and pressure gradient.We note that the ranges for s and α are not meant to represent realistically attainable values here, instead, we are more interested in the general structure of the AE over the domain.There are several interesting features visible in the figure .Even in this simplest model, the available energy exhibits rich structure over the s-α plane.More precisely, AE is large when s and α are comparable, s ∼ α, and is otherwise much smaller, particularly when the absolute value of one of these quantities is large.These findings are consistent with previous investigations (Rosenbluth & Sloan 1971;Dagazian & Paris 1982;Connor et al. 1983;Kessel et al. 1994;Strait et al. 1997;Rettig et al. 1997;Kinsey et al. 2006).It is also interesting to note that the precise reduction in AE depends on the drive: for a pure electron temperature gradient, significant Figure 6: Dependence of AE on magnetic shear s and pressure gradient α.In subplot (a), the Miller vector [ϵ, κ, δ, s κ , s δ , ∂ r R 0 , q, s, α] is set to [1/3, 3/2, 1/2, 0, 0, 0, 2, s, α], and other plots have the same vector with one parameter changed.In the subplot (b) the safety factor q is reduced from 2 to 1, for (c) the elongation κ decreases from 3/2 to 1/2, and in (d) the sign of triangularity δ is changed from 1/2 to −1/2.All plots have ωn = 1 and η = 0.
positive shear is more helpful in reducing AE, while AE driven by a pure density gradient benefits more from negative shear.
Since Eq. (2.62) can be integrated numerically to high precision, it serves as a useful benchmark for the more general AE of (2.44).Accordingly, we have compared the AE in the large-aspect-ratio limit with circular flux surfaces using a code that solves Eq. (2.44).This comparison is shown in Appendix B, and we find that the codes agree.

Miller geometry
We now leave the realm of the s-α limit and venture into shaped, finite-aspect-ratio equilibria.Our first step is to investigate the dependence on magnetic shear and pressure gradient for a range of different Miller vectors, and the results are shown in Fig. 6.Here we see similar trends as in section 3.2: negative shear and large α tend to be especially stabilising for a pure density gradient.However, it is also clear that the magnitude and precise contours depend strongly on the chosen Miller vector, as defined in Eq. (2.54).For example, it can be seen that lowering the safety factor is stabilising, since AE is reduced over a large region of the s-α plane as one compares subfigure (a) to (b).In subfigure (c) the elongation has been reduced produce a "comet"-type configuration (κ < 1, i.e. a horizontally elongated tokamak, see Fig. 2), which can increases the magnitude of the AE, Figure 7: The effect of the geometry on AE.In plot (a) the Miller vector [ϵ, κ, δ, s κ , s δ , ∂ r R 0 , q, s, α] is set to [1/3, κ, δ, 0, 0, 0, 2, 1, 0] with ωn = 1 and η = 0, and all other subfigures have the same Miller vector with one parameter changed.The contours shown in (b) have a higher inverse aspect ratio as ϵ is increased from 1/3 to 2/3, in (c) the pressure gradient α is increased from 0 to 1/2, and for (d) we have decreased the shear from 1 to 0. and the stabilising effects of s and α become less pronounced.Finally, in subfigure (d) the sign of the triangularity has been reversed to become negative.Although the shape of the contours remains largely unchanged, the peak in AE is changed to higher α and lower s, indicating that negative triangularity can be particularly beneficial in high-shear discharges with a modest value for α.In a more general sense, when changing any of the parameters significantly, one should expect that the precise shape and magnitude of the contours will change.
With this important caveat in mind, let us investigate the influence of geometry on the AE.To do so, we display the dependence on κ and δ for various Miller vectors in Fig. 7. Several interesting general trends can be observed.First, note that increasing the elongation beyond κ = 1 generally decreases the AE in all Miller vectors considered here, although the precise effect depends on the triangularity.Second, we see that it is not true in general that positive or negative triangularity is always stabilising; it depends on the other Miller parameters.Third, we see that tokamaks with κ < 1 and δ < 0, often referred to as (negative) comet cross sections tokamaks (Kesner et al. 1995), show a reduction in AE in plots (b) and (d), at least for the pure density gradient considered here.This is perhaps unsurprising, since such tokamaks are close to having the maximum-J property as shown by Miller et al. (1989).Since AE measures deviations from the maximum-J property, it is thus expected that these configurations perform well in terms of AE.
Investigating the plots in detail, in plot (a) one sees that negative triangularity is beneficial for κ > 1 as can be seen by the reduction in AE.In the following sections, we shall see that this is a consequence of the positive magnetic shear chosen.Next, note that doubling the inverse aspect ratio, as is done when going from (a) to (b), has a stabilising effect.Naively, one would expect that doubling the inverse aspect ratio would increase the AE by roughly a factor √ 2 ≈ 1.4, due to the factor √ ϵ in Eq. (2.45).However, going from plot (a) to (b) we see a decrease of the maximum AE by some 15%.This is likely due to the fact that, in a small-aspect-ratio device, magnetic field lines spend most of their time (or more precisely, arc-length) on the inboard side of the tokamak (Helander & Sigmar 2005).There, ω λ tends to be opposite to the drift wave and therefore these orbits do not contribute to the AE for a pure density gradient.It is also interesting to note that negative triangularity no longer exhibits a reduction in AE as the aspect ratio is significantly decreased, in accordance with the findings of Balestri et al. (2023).
Going from plot (a) to (c) the pressure gradient is increased from α = 0 to 1/2.With this introduction of pressure gradient, it can be seen that positive triangularity shows a decrease in AE, where negative triangularity does not.Finally, plot (d) has the magnetic shear reduced from s = 1 to s = 0 as compared to (a), which drastically changes the picture.Most importantly, we see that the lack of this positive magnetic shear results in negative triangularity no longer being stablising.We find that the results change somewhat if one instead imposes a pure electron temperature gradient (not shown here), though the basic trends remain intact.
All in all, we conclude from these results that the AE is very sensitive to equilibrium parameters, including quantities not investigated here such as q, s κ , and ∂ r R 0 .This sensitivity is perhaps reassuring: gyrokinetic turbulence has long been known to be strongly dependent on equilibrium parameters and even slight nudges can drastically change the picture (a sentiment perhaps best captured by the old Dutch expression wie het kleine niet eert, is het grote niet weerd ).We seem to reproduce a similar sensitivity in this AE-model for trapped electrons.This sensitivity becomes especially clear when investigating the dependence of AE on triangularity, which we shall discuss in the next section.

When is negative triangularity beneficial?
As hinted at in the previous section, it is not possible to make a general statement about the effect of negative triangularity on AE; its possible benefit depends strongly on other parameters describing the equilibrium.We can however find trends, and in order to do so we define the following fraction where δ = ±0.1 is chosen to represent a typical experimentally realizable range of parameters.This fraction can be interpreted as the factor by which the AE changes upon switching from positive to negative triangularity, where ∆ < 1 implies a reduction in AE.We present an investigation of ∆ and its dependencies in Fig. 8.We see two clear trends that seem to be robust for tokamaks with κ > 1.Firstly, as noted in the previous sections, in plots (a) and (b) we see that negative triangularity tends to be especially stabilising for configurations with significant positive shear.Similar conclusions were made by Merlo & Jenko (2023), who found that the turbulent energy flux in gyrokinetic simulations follows the same trend for TEM-driven turbulence: only for sufficiently high positive shear is a decrease in energy flux found at negative triangularity.Increasing α tends to push the ∆ = 1 line (in the plot this is the log ∆ = 0 line) to even higher values of shear, implying that a significant pressure gradient may make negative triangularity less desirable.Secondly, in plots (c) and (d) we note that negative triangularity can be beneficial in situations where the gradient is small, such as in the core.The dependence on η is non-trivial; at small density gradients a nonzero value of η can make negative triangularity beneficial.As in the previous sections, the results here depend on the Miller vector and are not meant to serve as a quantitative measure for core and edge transport.However, we have found that the presented trends tend to be robust as long as κ > 1 and thus do have qualitative value.We finally note that a more comprehensive model of the effect of negative triangularity should likely take collisions, impurities, and global effects into account (Merlo et al. 2019(Merlo et al. , 2021)).
From these results we infer that negative triangularity is expected to be especially beneficial in the core of the plasma, where gradients are necessarily small.It is not clear if the benefit extends to the edge: only with significant positive shear does negative triangularity become beneficial here as well.One should also keep in mind that ∆ measures the effect of going to negative triangularity while keeping all other parameters fixed.A more complete investigation would, for example, compare experimental equilibria with positive and negative triangularity, or use a global MHD-equilibrium code to find consistent profiles.We do not attempt such an investigation here, but we note that our mathematical framework would readily allow for such a comparison.We finally remark that the above results may seem counter-intuitive as negative triangularity is often thought to automatically imply TEM stabilisation, since the bounce points of most trapped particles reside on the inboard side of the torus, where the magnetic curvature should be favourable.Consequently, it is often argued that the bounce-averaged drift is such that TEMs are stabilised.Upon calculation of (2.37), we find no such stabilisation however, as explained further in Appendix C.

Gradient-threshold like behaviour
Our next step is to investigate the dependence of AE on the gradient strength ωn .From Eq. (2.8), one can show that there are two distinct scalings (Mackenbach et al. 2023a).In a strongly driven regime, one finds that the AE scales linearly with the gradient strength ωn .For a weakly driven regime one can expand around small ωn , and one finds that the AE scales with the gradient strength as (3.4) These scalings are reminiscent of gradient-threshold (or critical gradient) type behaviour (Dimits et al. 2000).Gradient thresholds are signified by a sudden decrease in energy flux when decreasing the gradient below some threshold value.The aforementioned scaling behaviour of the AE is displayed in Fig. 9 which similarly shows a rapid decrease below some threshold value.In plot (b) we estimate a critical threshold-like quantity from AE, by fitting a straight line to the strongly driven regime, i.e. we find the best-fit parameters a 0 and a 1 in the formula with ωn ≫ 1.The gradient threshold, denoted by ωc , is then defined as the interception with the abscissa, hence Figure 10: The gradient threshold as a function of equilibrium parameters.In plot (a) the Miller vector [ϵ, κ, δ, s κ , s δ , ∂ r R 0 , q, s, α] is [1/3, 3/2, 1/2, 0, 0, 0, 2, s, α], and plot (b) has [1/3, κ, δ, 0, 0, 0, 2, 1, 0].In all plots η = 0.
One could, of course, use different definitions for ωc , e.g. one could define the intersection point between the two straight lines on the log-log plot of Fig. 9 as ωc .However, we have found that the definition of Eq. (3.6) has several benefits: it is computationally cheaper, less prone to numerical noise, and seems to behave more smoothly.Other attempted definitions show the same trends.
We illustrate how ωc varies as a function of various equilibrium parameters in Fig. 10.Note that subplot (a) in Fig. 10 has the same Miller vector as Fig. 6 (a), and subplot (b) in Fig. 10 has the same Miller vector as Fig. 7 (a).Focussing on plot (a), an interesting trend is that increasing shear tends to increase ωc linearly, and ωc tends to plateau for low shear to some value.This is similar to the findings of Jenko et al. (2001), though their investigation focusses on electron-temperature gradient turbulence.It is also interesting to note that, in addition to the reduction in AE in the negative-triangularity configuration, it also benefits from a high critical gradient, which is in line with the findings of Merlo et al. (2015).This effect becomes even more pronounced as one increases the shear, which furthermore reduces the AE in the negative triangularity configuration.This implies that negative triangularity may be beneficial in a different sense: since the critical gradient estimated from AE is higher in negative triangularity geometries, the profiles may be able to sustain much higher gradients and thus higher core density/temperature.

Tokamak optimisation
In this section we aim to find AE-optimised tokamaks for a certain set of equilibrium parameters, at fixed gradients (ω n = 1 and η=0).To this end, we choose to optimise over κ and δ while keeping all other parameters fixed.In order to find somewhat realistic solutions, we restrict ourselves to a bounded optimisation space, namely The SHGO algorithm from Endres et al. (2018) is ideally suited for finding the global minimum in this low-dimensional bounded parameter space and is also available in scipy.
Finally, we shall vary magnetic shear and α, and investigate its effect on the global minimum found.
The results are displayed in Fig. 11, where the optimal values of AE, κ, and δ values are displayed as a function of s and α.For a visual aid of the shape of the cross sections, we refer to Fig. 2. It can be seen that both the optimal triangularity and elongation tend In all plots ωn = 1 and η = 0.
to be in the corners of the optimisation domain, and hence one should expect that these results are strongly dependent on this domain.Firstly, we see that vertically elongated tokamaks tend to be beneficial for all parameters considered here.It is furthermore interesting to note that the negative triangularity solution tends to be optimal whenever there is significant shear and the pressure gradient is not too large, which is consistent with the findings of Section 3.4.
From this plot, an important conclusion can be drawn: there is no such thing as a single "optimal" solution.The global minimum depends sensitively on other equilibrium parameters, such as shear and pressure gradient, which are, in turn, determined by the profiles of the safety factor, density, and temperature.Therefore, if one is interested in finding an AE-optimised tokamak, one should take care when choosing the profiles.One could also choose to let the profiles be part of the optimisation by describing them with some number of free parameters and constraints (e.g. one could use a fixed number of Fourier modes on top of a profile and optimise for the mode amplitudes).In reality, the profiles are themselves set by equilibrium conditions, making a self-consistent optimisation highly non-trivial.A more consistent investigation could perhaps solve this by coupling the current AE-model to a transport solver, which would calculate selfconsistent profiles.

Existence of solutions with high gradients yet low AE
In this section, we investigate how this AE model may relate to the suppression of TEMs when the density gradient is increased.To do so, we note several interesting properties that arise as one increases this gradient.First, the normalised pressure gradient α scales linearly with the density gradient (assuming a constant ratio of the poloidal magnetic field pressure to the thermal pressure, which e.g.occurs if one is operating at a fixed β-limit).The shear depends on the pressure gradient, as such a gradient drives the bootstrap current, which in turn changes the rotational transform profile.The bootstrap current density has an off-axis maximum in realistic scenarios, and such an off-axis maximum can locally lower the shear.This is most readily seen by inspecting the expression for shear in a large-aspect-ratio, circular tokamak, which depends on the current density profile j(r) as (3.8)where ȷ measures the average current density inside the radius r.From this expression, it is clear that for current density profiles that peak at r = 0, the shear is always positive.An off-axis maximum, supplied by the bootstrap current, can cause a locally lower shear.Hence, as one raises ωn one simultaneously increases α and decreases s.To estimate the magnitude of the effect of the bootstrap current on the shear, we note that the bootstrap current is proportional to the density and temperature gradients, and thus to the pressure gradient j b ≈ j b,0 α(r). (3.9) This is an approximation since the different transport coefficients relating the bootstrap current to the various gradients are not identical (Helander & Sigmar 2005), but we ignore this minor complication.We furthermore write the total current density as j = j b + j e , where j e is the equilibrium current, and assume j b ≪ j e = j e,0 ȷ(r).To first order in the smallness of the bootstrap current, (3.8) then gives Finally, following Miyamoto (2005) we estimate the ratio j b,0 /j e,0 as where β p is the local ratio of the thermal pressure over the poloidal magnetic field pressure, and the angular brackets denote a volume average.We shall take j b,0 /j e,0 to be on the order of 10%, implying that the shear may change as ds/dα ∼ s/10.Finally, one can relate the pressure gradient to ωn as where η i = ∂ r ln T i /∂ r ln n, with T i (r) being the ion temperature.We assume that the factor ϵβ p (1 + η + η i ) ∼ 0.1, so that dω n /dα ∼ 10.We illustrate the competing effects of the density gradient, pressure gradient, and shear in Fig. 12.In subfigures (a) and (b) we see various iso-contours of the AE in (ω n , s, α)-space, where (a) has positive triangularity and (b) has negative triangularity.It is especially interesting to note that in subfigure (a) there are paths in parameter space in which ωn increases but the AE decreases.These paths generally require that, as the density gradient increases, the pressure gradient should also increase and the shear should decrease.As we have argued, these trends are indeed found in tokamak discharges.One such path is indicated in subplot (a) as a blue line.Importantly, the blue line has which is the right order of magnitude for both ds/dα and dω n /dα.Subfigure (b) exhibits drastically different features.Planes of constant AE tend to lie parallel to planes of constant ωn , indicating that not much stabilisation is possible by changing the shear or the pressure gradient: the AE rises when ωn is increased.In subfigure (b), we again plot a line along the direction of increasing α and decreasing magnetic shear in red.Finally, note that for s δ we have used the estimate from Miller et al. (1998), s δ ≈ δ/ √ 1 − δ 2 .In subfigure (c) we display the AE along the blue and red lines given in subfigures (a) and (b) as a function of the density gradient.Note that the positive-triangularity case exhibits a distinct maximum, with low AE both to the left and right of the peak.One could interpret the existence of the latter as two distinct low-transport regimes; one with low gradients, and one with high gradients (which also has decreased magnetic shear and increased α).It is furthermore interesting to note that the negative-triangularity tokamak rises to far higher values in terms of AE and does not seem to drop back down to low levels along the chosen domain.Hence one could perhaps conclude that reaching a low-transport state with high gradients is not feasible in a negative-triangularity discharge.This is in line with findings of Saarelma et al. (2021) and Nelson et al. (2022), where the H-mode was found to be inaccessible in negative-triangularity tokamaks on basis of the ballooning instability, though the physical reason is of course different.This rise in AE in negative triangularity is perhaps unsurprising given that we have found that negative triangularity is stabilising in cases with significant positive shear, a weak pressure gradient, and a slight density gradient, exemplified in Figs. 8 and 11.Since, along the chosen path shear decreases and α increases with increasing density gradient, which is opposite to what is stabilising for negative-triangularity tokamaks, we see a sharp increase in AE.It may be feasible, however, to have a significant reduction in transport by tailoring the qprofile in such a way that negative triangularity becomes favorable, which likely implies significant positive shear.With such a reduction in AE, one could perhaps enjoy much improved transport whilst staying in an L-mode like regime.The parameters described in Marinoni et al. (2019) do seem to meet such requirements, especially near the edge where the reduction in transport seems greatest as compared to the positive triangularity case.A more comprehensive investigation, which shall be undertaken in a future publication, would self-consistently calculate the bootstrap current which would give precise paths in (ω n , α, s)-space.However, given the nature of the iso-contours in this three-dimensional space, we expect the observed trends to be robust, as long as the path has the correct general dependencies (i.e.decreasing shear and increasing α with increasing density gradient).

Conclusions
We have shown that it is possible to simplify the analytical expression for the AE of trapped electrons in the case of an omnigenous system, which speeds up calculations.If one furthermore employs an analytical local solution to the Grad-Shafranov equation, explicit expression of various quantities needed in the calculation of the AE (e.g., bounceaveraged drifts, bounce times) can be found as in Roach et al. (1995).Making use of an equilibrium parameterisation proposed by Miller et al. (1998), we go on to investigate how AE depends on these equilibrium parameters.Using this set-up, we observe several interesting features of the AE: (i) A comparison is made between AE and tglf.We observe a fairly good correlation between energy flux and A 3/2 , indicating that AE can be a useful measure for tokamak transport.
(ii) Increasing the magnitude of the magnetic shear or increasing the Shafranov shift tends to be stabilising as indicated by a reduction in the AE, and these trends hold for many different choices of geometry.Especially negative shear reduces the AE substantially for pure density gradients.
(iii) Vertical elongation tends to be stablising, as indicated by a reduction in AE.Negative triangularity can be stabilising, particularly in configurations with significant positive shear or small gradients, but not always.
(iv) The AE has different scalings with respect to the gradient strength in weakly and strongly driven regimes.We employ this difference in scaling to estimate a gradientthreshold like quantity, and we find that it has similar behaviours as found in criticalgradient literature; an increase in shear tends to increase this gradient-threshold and negative triangularity benefits from an especially high gradient-threshold.
(v) Using AE for shape-optimisation we show that the optimal solution is strongly dependent on pressure gradients and magnetic shear, implying that the optimisation is sensitive to the density, pressure, and q-profiles.
(vi) An investigation is presented on how AE varies as the density and pressure gradient increase consistently, while shear decreases.We find that in such scenarios one can find solutions with large gradients yet low AE.Such solutions tend to exist for positive triangularity tokamaks but not for negative triangularity tokamaks.
The results suggest that various observed trends regarding turbulent transport in tokamaks may partly be understood in terms of AE, which has a simple physical interpretation and is cheap to compute.The analytical framework can readily be extended to account for an equilibrium model which allows for other shaping and plasma parameters such as plasma rotation (Hameiri 1983;Miller et al. 1995), squareness (Turnbull et al. 1999), and up-down asymmetry (Rodrigues & Coroado 2018), though no such investigation is presented here.
Figure 14: Comparison of the AE as calculated with two different codes.Plot (a) is calculated using a code that calculates AE in the s-α limit, and (b) is calculated using the Miller code.The plots are visually indistinguishable.Calculated using q = 2, ϵ = 10 −6 , ωn = 3, and η = 1.All other parameters for Miller are set to zero, as required in the limit of the s-α tokamak.Plot (c) presents the relative error, where it can be seen that the relative error is very small for large regions of s-α space.Note that the colourbar scale in plot (c) is logarithmic.
comparison presented here we use 10 3 equidistant nodes for θ.The integral over the pitch angle is done using quadrature methods.
The comparison is shown in Fig. 14.In this figure, three different contour plots are shown; (a) is the available energy as calculated from Eq. (2.62).Plot (b) shows the result as calculated from Eq. (2.45).Finally, plot (c) shows the relative error between the two codes (more precisely, it is the difference between plot (a) and (b), divided by plot (a)).It can be seen that the error is typically quite small, with a maximal value of 1% and a mean value of 0.004%.If different parameters are chosen (safety factor, density gradient, or η), the error remains similarly small.
All plots presented in the current publication are generated using the same or even more refined numerical parameters as used here, so that we have a high degree of confidence that the presented trends are indeed physical and not numerical.Further convergence checks (increasing the resolution of θ and adjusting the tolerances of the quadrature methods) do not alter the plots presented in this publication in a visually discernible manner.
As an additional check, we highlight that a recent publication has evaluated the AE of trapped electrons in quasi-symmetric systems (which includes tokamaks) in two asymptotic limits: those of a very strong and a very weak density gradient (Rodriguez & Mackenbach 2023).It was found that the AE scales with elongation as A ∝ κ 1/4 if the density gradient is sufficiently small, and A ∝ κ −3/4 if the density gradient is sufficiently large.Importantly, this analysis assumed fixed ∂n/∂r eff and r eff /R 0 , instead of fixed ∂ r n and r/R 0 .If one properly accounts for this different definition of the radial coordinate, we find that the code reproduces the correct scaling behaviours, as may be seen in Fig. 15.We also note that in the aforementioned investigation it was found that at zero shear, more negative triangularity is found to increase the AE if the gradient is sufficiently strong and κ > 1.If the density gradient is sufficiently strong and κ < 1, the AE decreases with more negative triangularity.These trends are reproduced and can be found in Fig. 7, subplot (d).In this section, we investigate the difference in trapped particle orbits in positive and negative triangularity tokamaks.To this end, we investigate the dependence of Eq. (2.37) on δ, and we set the other components of the Miller vector equal to [ϵ, κ, δ, s κ , s δ , ∂ r R 0 , q, s, α] = [1/3, 2, δ, 0, 0, 0, 2, 0, 0].The result for a positive triangularity tokamak (δ = 0.5) is plotted in Fig. 16, where we have plotted ω λ as a function of its bounce points θ, which satisfy 1 − λ B(θ) = 0.
(C 1) Furthermore, we have shown the AE per λ, called A λ , which is the integrand of Eq. (2.45).This is done by coloring a line of constant λ (which corresponds to constant B) according to its A λ .Finally, we also display ω λ as a function of the trapping parameter k 2 which maps λ → [0, 1] according to where the subscripts max and min refer to the maximal and minimal values of the functions respectively.With this convention, k 2 = 0 corresponds to the most deeply trapped particles and k 2 = 1 to the most shallowly trapped particles.We have furthermore included a red dashed line, which delineates where ω λ changes sign, which determines stability in a purely density-gradient-driven TEM.In the figure, ω λ > 0 corresponds to instability (and associated AE).It can be seen that this positive triangularity tokamak is unstable up to roughly k 2 = 1/2, and the magnetic well is relatively narrow The same information is displayed for a tokamak which has δ = −0.5 in Fig. 17.It can be seen that the precession frequencies are unstable for a broader range of values for k 2 .The AE is furthermore weighted by the bounce-time of a particle, which can become very large at the bottom of a magnetic well in a negative triangularity tokamak.As such, the negative triangularity case (with the Miller vectors as chosen here) has higher AE than the positive triangularity case.
We have tried various numerical experiments to assess the origin of this difference.From Eq. (2.37), we note that the term involving 1 − λ B ∝ v 2 ∥ , and hence we identify this term as the curvature component of the drift.The term involving λ B ∝ v 2 ⊥ on the other hand we identify as the gradient drift.Setting the term involving the parallel velocities equal to zero results in the found trends inverting, showcasing that this drive plays an important part in determining the precession.The poloidal curvature, R −1 c , furthermore plays an important part.By setting this term equal to one in Eq. (2.37), we also find that negative triangularity is preferred over positive triangularity.Therefore, we postulate that this curvature drift plays an important part in determining stability.Importantly, the particles that experience curvature drive in negative triangularity tokamaks are the deeply trapped particles, which tend to be most unstable against the TEM with a density gradient.This is in contrast to positive triangularity tokamaks, where the most shallowly trapped particles experience significant curvature drive.These shallowly trapped particles however, are stabilised by the fact that they experience an averaged drift, and as such the curvature drive here is less deleterious.

Figure 1 :
Figure 1: Contour plot of I z as a function of c 0 and c 1 .

Figure 3 :
Figure 3: Comparison between codes, showing the correspondence between the estimate of the energy flux from AE and tglf.The top row displays the energy flux from tglf, and the bottom row displays the corresponding estimate from AE.One can see agreement in trends, though some details differ.The left column has aMiller vector [ϵ, κ, δ, s κ , s δ , ∂ r R 0 , q, s, α] of [1/3, κ, δ, 0, 0, 0, 2, 0, 0], whereas the middle column has [1/3, κ, δ, 0, 0, 0, 2, 1, 0], and the right-most column has as Miller vector of [1/3, 1, 1/2, 0, 0, 0, 2, s, α].The first column has ωn = 3, the second and third column have ωn = 6, and η = 1 in all plots.The white masked-out regions have a dominant instability which is not in the electron direction, and as such we filter them out.Finally, note that the colour bars are the same scale in each column.

Figure 4 :
Figure 4: A scatter plot showing the relation between the two estimates of the non-linear energy flux.The red dashed line is shows the expected linear relationship, x = y.The plot consists of N = 1171 points.The gray points have some transparency, and as such darker regions arise due to a high density of points.We have furthermore added simulation data from the gyrokinetic code gene (Jenko et al. 2001), also presented in Mackenbach et al. (2022), as blue markers.

Figure 5 :
Figure5: The AE of a large aspect ratio circular tokamak, as a function of magnetic shear s and pressure gradient α.The plots have been generated using q = 2. Plot (a) has ωn = 3 and η = 0, whereas plot (b) has a pure electron temperature gradient, i.e. ωn = 0 and ωn • η = 3.

Figure 9 :
Figure 9: Example of dependence of AE on gradient strength.Two scalings are found in plot (a).In plot (b) we define a gradient threshold by fitting a straight line to the strongly driven regime, and finding its ωn interception with the abscissa.

Figure 15 :
Figure 15: Available energy as a function of the elongation κ for a very weak density gradient in the left plot (ω n = 1/100 at κ = 1), and a very strong density gradient in the right plot (ω n = 100 at κ = 1).

Figure 16 :
Figure 16: The precession frequency and AE distribution for a positive triangularity tokamak

Figure 17 :
Figure 17: The precession frequency and AE distribution for a negative triangularity tokamak