Transport Barriers in magnetized plasmas- general theory with dynamical constraints

A fundamental dynamical constraint—that fluctuation induced charge-weighted particle flux must vanish- can prevent instabilities from accessing the free energy in the strong gradients characteristic of Transport Barriers (TBs). Density gradients, when large enough, lead to a violation of the constraint and hence preclude unstable modes and turbulent transport. This mechanism, then, broadens the class of configurations (in magnetized plasmas) where these high confinement states can be formed and sustained. The need for velocity shear, the conventional agent for TB formation, is obviated. The most important ramifications of the constraint is to permit a charting out of the domains conducive to TB formation and hence to optimally confined fusion worthy states; the detailed investigation is conducted through new analytic methods and extensive gyrokinetic simulations.


Introduction
This paper is an investigation into the fundamental physics of Transport Barriers (TBs) in magnetically confined plasmas [1][2][3][4][5].Over a distance of centimeters, experimental TBs routinely sustain temperature differences exceeding those from the core to edge of the Sun.Crucially, the rates of heat loss in these enormous gradients are surprisingly small.This is the spectacular property that is relied upon in many proposals 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.
to produce fusion gain in terrestrial plasmas.By conducting an enquiry into the physical mechanisms responsible for the existence of such a desirable but 'unwarranted' thermodynamic state, we plan to chart out parameter regimes that are most hospitable to TBs with the highest confinement.
In the conventional gyrokinetic stability analysis, the TB formation is attributed to the suppression of turbulent transport associated with the dominant instability-the coupled Ion Temperature Gradient and Trapped Electron Mode ITG/TEM [6][7][8][9][10]; the suppressing agent is usually the E × B (velocity) shear.As experiments progress to fusion conditions, it was/is feared that the paucity of velocity shear could eliminate/weaken transport suppression.Fortunately, however, there are many experiments that do observe TBs precisely in the low shear regime [11][12][13].What mechanism, then, causes and maintains this fusion relevant low transport state?

The charge-weighted flux constraint
The answer, amazingly, lies in a very general property of strongly magnetized plasmas; there is a fundamental physical constraint that gyrokinetic fluctuations must obey to exist: the fluctuation driven net charged-weighted flux must be zero (s is the species index) ∑ s q s Γ rs = 0, where q s is the charge, and Γ rs is the usual radial transport flux of particles (i.e. the flux averaged over space scales larger than the fluctuations).For linear eigenmodes, this is the conventional quasi-linear flux.This 'flux constraint' (FC) is satisfied by low frequency plasma fluctuations that obey Maxwell's equations to the lowest order in the gyrokinetic expansion.
A small digression on the origin of the FC, epitomized in equation (1), is in order.It could be derived, for instance, from: a) collective momentum conservation after summing over fluctuation interactions in space, as in Coulomb collisional transport, b) local frame invariance of local fluctuations, c) primarily electrodynamics considerations together with their invariances.(See [14][15][16][17][18]).
Although this dynamical constraint has been known for some time, it is for the first time that its enormous power and effectiveness will be explored in the context of TB formation.We will focus on the suppression of two important instabilities in the gyrokinetic system-the coupled Ion Temperature Gradient and Trapped Electron Mode ITG/TEM-since these instabilities are typically the limiting factor for confinement.If these instabilities can be tamed, drastic improvements to confinement can be achieved.
A simple translation of the content of equation ( 1) reads: if it is not possible for a class of fluctuations to satisfy the FC no instability is possible, no matter how much free energy may be available in the equilibrium gradients.In essence, the violation of the FC will not let plasma relax the gradients via instability, and configurations with large gradients (TBs) will be stable!Tracing the amazing TB phenomenon to these fundamental roots is not only beautiful but it also greatly clarifies the conditions necessary for these crucially important, thermodynamically counter-intuitive states to be sustainable.The fluctuation dynamics are therefore much more clearly understood as a consequence of two very different dynamics operating simultaneously: the aforementioned constraint on the total chargeweighted flux, and the well-known relationship of fluctuation growth with the plasma free energy.A major accomplishment of this work lies in distinguishing these two different dynamics, their separate roles, and their very different characteristics.
Oftentimes, we are able to analytically derive approximate expressions for conditions when the constraint is soluble/insoluble.These are simple and useful: specifically, boundaries, beyond which the constraint cannot be satisfied for any growing fluctuation, and hence regions with no instabilities are identified.Extensive gyrokinetic simulations (conducted via the gyrokinetic code GENE [19]) show that these analytical boundaries explain a great deal of the observed behavior, for both linear instabilities, and also, nonlinear fluctuations.The theory-simulation combination makes the central role of the constraint 'visible' for the first time.
The invocation of the FC, thus, brings the gross stability of the very complex gyrokinetic system to the realm of (partial) analytical tractability.Combined with detailed simulations, we, then, have a powerful tool to deeply investigate fluctuations.

Free energy and instability
Whenever there is free energy (in equilibrium gradients in magnetically confined plasmas, for instance), growing instabilities can arise if they reduce the equilibrium free energy.For the gyrokinetic system, one can derive an evolution equation for the fluctuation-induced free energy [14,[20][21][22][23][24][25], where δf is the perturbed distribution function, and ) are, respectively, the particle and heat fluxes, and C is the collision operator.For simplicity of presentation, we have assumed k ⊥ ρ is small in equation ( 2), but this is not necessary.
In equation (2), fluctuation free energy grows because the free energy in the macro (equilibrium) scale is reduced by the same amount (due to quasi-linear fluxes), minus collisional entropy production (as in Boltzmann's H theorem).The Left Hand Side (LHS) is the rate of change of the free energy in the fluctuations.(As an aside, note that this expression can be shown to be closely related to the familiar definition U-TS for small fluctuations, where U is the kinetic plus field energy, and S is the entropy ∼f ln f for small fluctuations about a Maxwellian, as in the gyrokinetic ordering).The RHS is the rate of change of the equilibrium free energy, and has a classic structure: a product of thermodynamic forces dT s /dx and dn s /dx times the corresponding thermodynamic fluxes Q s and Γ s , respectively.Free energy in the equilibrium is transferred to the fluctuations, but the collisional term changes entropy (increasing it).
The growth rate of an instability is determined by how efficiently the instability causes this transfer.As expected, for the gyrokinetic system, the free energy transfer can be greatly increased by familiar factors such as 'bad' magnetic curvature drifts of ions and electrons, resonances in phase space, etc.Such specific physics is contained in the very general equation (2), although some manipulations are needed to render it more apparent.For example, in simplified models, it can be shown explicitly that 'bad' curvature affects equation (2) to give higher growth rates; its manifest effects will be seen in our calculations.
If the free energy dynamics were operating unchecked, steeper gradients 'would' engender faster instabilities for linear modes.And nonlinearly, larger fluctuation amplitudes should result from larger gradients, giving larger fluxes.This behavior is indeed seen in a truly vast variety of physical systems, including in realms far outside of magnetized plasma dynamics.
Fluctuations in TBs, however, behave otherwise!This behavior is both theoretically striking and of great practical importance: it is crucial to many approaches to attain thermonuclear energy gain.That this 'deviant' behavior is an ineluctable consequence of the FC equation (1), constitutes the theme story of this paper.
This story will be built by investigating, analytically as well as by detailed simulations, the historically most pertinent instability for TB formation, the ITG/TEM, whose suppression by velocity shear forms the conventional foundational physics of the phenomenon.But if TB formation is to happen without velocity shear, we must show how the FC could suppress ITG/TEM despite all conventional free energy drives.In the process of demonstration, we will develop powerful new results and insights.The distinct manifestations of the two dynamics (free energy driven plus the dictates of the FC), clearly revealed in simulations, provides us with a deeper understanding of the physics of gyrokinetic fluctuations.

Glimpses of what constrained dynamics can do
It is readily seen that the FC analysis, when applied to the simplest linear ITG/TEM model, is equivalent to the traditional dispersion relation (see [14] for details).The fundamental power of the FC is revealed, however, when we examine fluctuations in their full generality; the injection of the FC into the picture will yield results that have never been obtained through dispersion relations.Some examples are: a) Analytic derivations of bounds on the stability (in realistic geometry) where the dispersion relation is, analytically, quite intractable.The bounds are verified by simulations with surprising accuracy.b) Analytic descriptions of complex mode structure changes that have direct implications for nonlinear heat fluxes.Except for invoking the solubility of the FC, no comprehensible description could be constructed especially in realistic geometries where a kinetic dispersion relation is intractable.c) Qualitative insights into the stability behavior when gradients are strong.What, for instance, are the relevant variables that control the system behavior.Attempts to do this without using the dual dynamics approach fall subject to multiple contradictions-this will be explicitly demonstrated in Sec.XIII where the ITG/TEM model includes trapped electrons.d) The striking 'adaptability' of fluctuations.Like many similar complex systems, the gyrokinetic system adapts to produce the most unstable mode to decrease the free energy.
Simulations show that the ITG/TEM usually 'finds a way' to tap free energy until very close to the boundary for which the FC becomes insoluble.In other words, the gyrokinetic system does in fact behave just like the vast variety of physical systems with many degrees of freedom.The apparent violation of ubiquitous thermodynamic behavior in TBs is actually just a failure to appreciate the presence of the FC.The extra understanding accruing from this fundamental statistical property of mechanical systems with many degrees of freedom is one of the distinguishing features of this work and will be an essential element in building what we call the statistical mechanical ansatz that will guide our attempts to understand the essential behavior of the system.Adaptivity becomes apparent once the FC is understood, so that the fluctuation behavior as solubility boundaries are approached can be seen as an obvious attempt to remain unstable.
We end this section with several general explanatory remarks before we begin to present detailed calculations.Our calculations will be limited to a deep exploration of the ITG/TEM mode using both analytical methods and extensive simulations via the well known gyrokinetic code GENE.
The crucial insight that animates this work is that there exist equilibrium conditions for which it is impossible for the fluctuation-induced (charge-weighted) flux to satisfy the FC.Since the constraint has no particular knowledge of eigenmodes, the FC analysis is valid for all fluctuations, eigenmodes or not.An important conceptual ramification of our approach is that an inability to find an unstable solution could originate in circumstances that are much broader than the eigenmode problem could ever predict.
Using the new methods of constrained dynamics, deriving general conditions when no gyrokinetic fluctuations can grow turns out to be enormously easier than a corresponding eigenmode calculation.When analytic bounds can be calculated, they are remarkably close to what extensive gyrokinetic simulations find.These stability predictions often defy expectations from 'commonsense' heuristic viewpoints.They also predict completely new types of mode behavior, confirmed by extensive numerical simulations, that have not been previously appreciated.And finally, simulations show that these considerations apply to nonlinear fluctuations and transport just as for linear instabilities.
Exploring conditions for which the fluctuation-driven transport cannot satisfy FC will be our route to figuring out what makes TBs possible when, for example, velocity shear is low; the latter is, of course, a matter of great practical importance.
The reader will likely find that this rather different mode of probing stability is unfamiliar, primarily because we transcend here the methods of conventional plasma instability theory.There are good reasons, however, to stick to and rely on the FC-based stability analysis: 1) simulations in complex geometries follow criterion derived by simple methods 2) results can be qualitatively traced to fundamental statistical mechanical concepts that are far more general than gyrokinetic stability theory.Rather, those general considerations are applied to gyrokinetics as a specific context of interest.
It is of essence to emphasize that the FC pertains only to the particle and not the thermal flux, thus valorizing density gradients (as the diagonal thermodynamic drive for the particle flux).Consequently, one may then expect that there will exist parameter regimes where increasing equilibrium density gradients will eventually compel the driven fluxes to be non-zero violating the FC.Since no growing fluctuations could exist, we could get a stable equilibrium state with high gradients, the TB.
Perhaps the most telling result of FC analysis is that equilibrium density gradients above a threshold emerge as a critical ingredient to TB formation (at least when ITG/TEM is the controlling fluctuation).
Since equilibrium gradients are also a measure of the free energy drive (for instabilities), it might appear somewhat paradoxical that a system with gradients greater than a threshold can become fluctuation-free.This reflects, in fact, the very power of the deeper physics which tells us that no matter how large the free energy, the insolubility of the FC guarantees an essentially instability-free state.The TB, thus, is a perfectly warranted high gradient equilibrium!

A model calculation: the ITG instability with adiabatic electrons
The ITG mode is a major instability which is often the dominant source of energy loss in tokamaks.With the approximation of adiabatic electrons, its physics is simple enough that solubility of the FC can be tested analytically.We call this mode the ITG ae .Though limited in its applicability, it qualitatively captures all features of the interplay of the FC, free energy balance and adaptivity that pertain to the full electron case: how the FC leads to stabilization, how it compels serious modifications of eigenfunction structure, how the gyrokinetic system responds adaptively to the FC and to the spatial structure of curvature, and how such adaptivity implies that only the FC results in stability for all but the most extreme geometries.And in many TBs, cases with full electron dynamics sometimes quantitatively approach the ITG ae solubility bounds section 13.The relative simplicity of the ITG ae makes it an ideal model to explore.
It is well known that the stability of the slab-like ITG ae is controlled by the parameter η = (1/L T )/(1/L n ), where 1/L n and 1/L T are the usual inverse gradient scale lengths of density and temperature.The mode is unstable only when η exceeds some critical value η crit .So for a given temperature gradient, the slab ITG ae can be stabilized by a sufficient density gradient.
We will show that the stabilization of the slab ITG ae by low η is a consequence of the insolubility of the FC.Once this connection is recognized, a large variety of fluctuation behaviors can be simply understood, even aspects that were previously not.
In toroidal geometries, the curvature-driven ITG ae is both more unstable and more complex.In the next section, we derive a very simple and general criterion for the amount of density gradient beyond which the FC is insoluble.Since the density gradient is the primary cause of stability, it is convenient to define a more suitable normalized density gradient parameter (i.e.other than η): the Fraction of density gradient in the Pressure gradient.F P = 1/(1 + η) will be the parameter of choice in this study.We will soon derive that, for ITG ae , the FC is insoluble when the density gradient is large enough to make Note that the density gradient magnitudes corresponding to equation ( 4) sometimes do actually arise in ITBs and pedestals, and, therefore, are experimentally relevant.Such density gradients are crucial to allowing the existence of the TB if velocity shear is low, as we will see.
One must emphasize here that equation ( 4) is not just another way of writing the conventional η criterion for linear stability; it represents altogether different physics.
In fact, when examined in the context of our past experience with linear theory, this result is extraordinarily surprising.Let us contrast the simplicity of equation ( 4) with well known critical temperature gradient criterion in terms of 1/L Tcrit ; the latter is found by decreasing the temperature gradient, and scanning over all k θ ρ i , until the maximum γ vanishes.In toroidal geometries 1/L Tcrit varies by a factor 2-4 or more depending upon the safety factor, magnetic shear, the temperature ratio T i /T e , etc.The fact that the constrained dynamics lead to a stability criterion that involves only the fractional density gradient and has no cognizance of the multidimensional parameter space that determines the conventional stability attests to its simplicity, and, as we will verify shortly, universality.This is the first example where simulation behavior cleanly distinguishes between considerations of the FC and of energetics.Energetics crucially determine γ and 1/L Tcrit , but not the critical F P .
Let us put this in an even larger physical perspective.Consider the following thought experiment: we maintain a strong temperature gradient at a constant value, and increase the density gradient.Eventually, stability results.From the general perspective of statistical mechanics, this behavior is qualitatively different than in a vast variety of systems.Even when temperature gradient magnitudes are extremely high, and when factors like the curvature drive are highly destabilizing, and despite there being an extremely large number of degrees of freedom (and hence, a very large number of potential mode structures), none of these succeeds in tapping the strong gradient free energy.
It will be a misconception to think that the adiabatic electron assumption, by itself, offers an explanation for this unexpected result.Since adiabatic electrons carry no density fluxes, the ITG ae cannot relax density gradients; density gradients, thus, are not a drive for the instability.This merely implies that increasing the density gradient should not increase the growth rate.In no way, can it explain why density gradients drive the growth rate to zero.
What explains this amazing phenomenon is the insolubility of the FC (equation ( 4)-No fluctuations can grow (in any geometry for any plasma parameters) when FC is violated-no relaxation routes are accessible.
We will now reinforce the power of the FC route to stability by presenting gyrokinetic simulations for different and realistic geometries and parameters.In particular, we show • the factors that strongly affect the free energy equation may have little impact upon the FC solubility.• the response of fluctuations to each dynamic is cleanly distinguished.
To the best of our knowledge the results below are new.Simulations were carried out using a standard magnetic geometry (so-called 'Miller geometry' with representative shaping) for different temperature gradients, geometric and other equilibrium parameters.In all cases, we fix the temperature gradient at some value and scan the density gradient until stability is found.The growth rate of the most unstable mode (scanning k θ ρ i ) to resolve the peak γ, is plotted.
We start with a standard case in figure 1(a) (Cyclone Base Case with elongation and triangularity).We choose a wide range of temperature gradients (R/L T ) up to very high values (R/L Ti = 10, 20, 40 and 80, where R is the torus major radius).The magnitude of the equilibrium gradient obviously strongly affects the free energy equation, but not the solubility bound in terms of F P .The simulation results reflect this fact: growth rates vary by an order of magnitude, but the critical F P varies by ∼5% (a hundred times less) and is always quite close to the analytic upper bound for the solubility of the FC.
Without an understanding based on the solubility of the FC, this result would be quite perplexing.
We now display stability curves for a variety of magnetic geometry parameters: magnetic shear, safety factor q, and Shafranov shift.As we traverse this space, the critical temperature gradient R/L Tcrit varies by ∼270% (figure 1(b)) while the critical F P varies only by 5% (figure 1(c)), about two orders of magnitude less.And to boot, the stabilizing F P is always close to the analytic upper bound imposed by the FC.Finally, in figure 1(d), we show that the stability curves for representative stellarator geometries from Wendelstein 7X [26] and the proposed NCSX [27], are similar to that of tokamaks.(These stellarator results were for VMEC geometries for NCSX and W7X courtesy of M. Zarnstorff, P. Xanthopoulos, I. McKinney and M.J. Pueschel.) Vast experience with ITG has told us that the temperature ratio T i /T e strongly affects the growth rate and critical temperature gradient.For the results shown in figure 1(e), R/L Tcrit varies by 400% as T i /T e is varied between 1/4 and 4. The growth rate (normalized to the ion transit frequency v thi /R varies by ∼500%.This is because T i /T e affects the relative magnitude of the adiabatic electron contribution to the free energy in equation (2).But since adiabatic electrons contribute no particle flux, we expect no T i /T e effect upon the upper bound of F P .Simulations find exactly that-the critical F p varies by ∼4% (again about two orders of magnitude less than the variation in γ or R/L Tcrit ) and is close to the computed solubility limit.
Finally, we consider the ETG, which is well-known to be nearly isomorphic to the ITG ae .We compare a case with finite Debye length, λ D ρ e = 1/2, to one with λ D ρ e = 0 in figure 1(f ) Somewhat similarly to T i /T e , the finite Debye length gives an additional contribution to the free energy equation ∼δE 2 that lowers the growth rate by a factor ∼2, but this does not contribute to the particle flux, and the FC bound is unaffected.The critical F P is affected by only ∼10%, far less than the growth rate.
Let us summarize the simulation results reported above.They illuminate, for the first time, a crucial point about the deep underlying physics of the gyrokinetic system: effects that strongly affect the free energetic dynamics may not, often do not impact the one parameter solubility limit due to the FC, equation (4).The latter is what controls the stabilization of the ITG ae by density gradients; the predicted threshold is totally insensitive to numerous effects that strongly affect free energy dynamics.
Since the solubility of the FC emerges as the key dynamic that commands the simultaneous existence of steep gradients and weak instabilities, and hence of TB formation, it is this dynamic that we will continue exploring and testing in this study.
The simulations also highlight the 'adaptability' of the gyrokinetic system.When temperature gradients are high, the modes become stable only quite close to the solubility boundary of the FC, despite huge variations in parameters.One does intuitively expect that a system with many degrees of freedom would adapt -'find a way' to tap strong equilibrium free energy to relax it if that were at all possible.The gyrokinetic system, just like other similarly complex physical systems, does exactly that but for a 'hard' constraint in the dynamics that precludes accessibility in some regions of the equilibrium parameter space.Even then, the system remains unstable until very close to the FC boundary.
This latter picture is further reinforced by simulation results that follow, which are in many ways much more striking than the ones in figure 1.But first, we indicate how the solubility criterion for the FC can be obtained analytically.The growth rate of the most unstable mode (maximum for all k θ ρ i ) is plotted (a) vs F P for a standard geometry like the cyclone base case) (b) for four diverse geometries vs R/L T showing that the critical temperature gradient varies by ∼270% (c) for those same geometries, the critical F P varies by ∼5% (d) for two stellarator geometries (e) for an ITB geometry, for many values of T i /Te (f ) for a low aspect ratio A geometry, an ETG for two different λ Debye .

Mean field theory for computing FC bounds: Simplified Kinetic Model
The complexity of the gyrokinetic system suggests developing a mean field theory as is, very effectively, done in other branches of physics.Interestingly, the Simplified Kinetic Model (SKiM), constructed here, leads to exactly the same solubility bounds as the somewhat more accurate version in appendix B.
We begin by quantifying bounds on the solubility of the FC in this simple system, to arrive at analytic expressions.We can then use gyrokinetic simulations for realistic parameters to examine how well these bounds describe the behavior of the system.We find that bounds derived using SKiM have a remarkable agreement with realistic simulations.The ultimate justification for this model is the high degree of agreement that its predictions have with simulations.However, it adds value beyond simulations, particularly in elucidating basic conceptual issues.
The mean field theory on the gyrokinetic equation (in ballooning coordinates) will be constructed by a statistical averaging on the variable l, the distance along a field line; l characterizes the conventional 'large number of degrees of freedom' of this system.It is straightforward to systematically expand the gyrokinetic equation in large ω and small k ⊥ ρ i to derive a dispersion relation that will contain certain eigenfunction averages: These have a transparent interpretation as an eigenfunction average ion drift frequency < ω di >, parallel wave number < k ∥ >, and perpendicular wavenumber < k ⊥ >.We can consider these eigenfunction averages to encapsulate many degrees of freedom in a single average quantity, conceptually like in a mean field theory.
The gyrokinetic distribution function in a geometry with constant k ∥ , k ⊥ and < ω di > is easy to solve (it is equivalent to a straight field line geometry with constant curvature, which can be Fourier transformed in space).In terms of the eigenfunction averages defined above, the non-adiabatic part of the distribution function (s is the ion species index) is where ω ⋆ s is diamagnetic frequency, and is the sum of a contribution from density gradients ω ⋆ ns and temperature gradients ω ⋆ Ts (mv 2 /2T − 3/2), and C is some measure of a collision frequency.
The ITG ae dispersion relation, obtained by using the usual quasineutrality condition must be solved numerically, even for a single species.We can compare the analytic SKiM dispersion relation to simulation results only after using the eigenfunction provided by simulations.There is good qualitative, often, semi-quantitative agreement between SKiM and the simulation growth rates.We can derive more elaborate eigenfunction averages than equations ( 5)- (7), to lead to better agreement with simulations for given simulation eigenfunctions.But the intuitively clear expressions in equations ( 5)-( 7) will suffice for our present purposes; they encapsulate in a tractable way the average effect of the dynamics in many degrees of freedom along a field line.
The particle flux from the non-adiabatic response equation ( 8) may be computed (for a single ion species) as for an instability with growth rate γ > 0 and real frequency ω r .
(Note if γ ⩽ 0 additional terms arise due to the Landau contour.)We have made one small simplifying assumption above: that gradients are steep, as in a TB, so that ω ⋆ s >> ω.This is quite well satisfied for ITG ae and the ITG/TEM discussed in section 13.
It is, now, trivial to note that if 1/L n > 3/(2L T ), the RHS of equation ( 10) is positive definite; consequently the flux condition cannot be satisfied.This condition, translating as F P > 0.6 guarantees stability of all fluctuations (for any < k ⊥ > , < k ∥ >, < ω d >), and any real frequency ω r and γ > 0).It is equally evident that stability could be possible for F P < 0.6 because of the term proportional to v 2 ; the latter contribution, however, will depend on the mode details.
The stability criterion F P > 0.6 is to be contrasted with its typical counterpart from a dispersion relation.The latter is complicated enough that numerical solution is required, even for this simplified model, and there are many well-known dependencies of the eigenvalue upon a host of parameters.On the other hand, the criterion F P > 0.6 is easily obtained from very brief analytical arguments, is far simpler in character, takes no cognizance of geometrical complexity or fluctuation attributes, and is true for all fluctuations, eigenmodes or not; it is universal, in particular, it holds for fluctuations with any temporal or spatial structures, any spectrum in ω and k y .
And the criterion also predicts the behavior of far more complex simulations in realistic geometry as shown above and below.
A similar criterion for instability has been derived from a dispersion relation like equation (9) using a Nyquist analysis [28].The present analysis is simpler, and more crucially, reveals that the stability criterion follows only from the solubility of the FC, and not from the dispersion relation as a whole.One consequence is that the criterion does not depend upon energetic considerations.Also, fluctuations must satisfy the FC in realistic geometry; it is not just a feature of the model geometry here; FC considerations have great generality.
And there is a deeper point.The criterion for solubility of the FC found here applies to any fluctuation, whether or not it satisfies the dispersion relation.Therefore, the insolubility of the FC cannot be due to considerations that pertain specifically to eigenmodes whose study constitutes the bulk of conventional stability theory!The reasons why the particle flux can be always positive (in the direction down the density gradient) for sufficiently large density gradients, is to be found in fundamental physical considerations: Since the diagonal fluxes driven by thermodynamic forces are positive definite (see the discussion in section 12), a sufficiently large density gradient will, eventually, cause the LHS of equation (10) to be nonzero; the FC will, then, be insoluble and the mode will be stable.
The FC insolubility criterion is a very robust result, takes no cognizance of plasma, magnetic or fluctuation properties.This is, perhaps, why the insolubility bound of SKiM agrees so well with simulations in much more complicated geometries: the insolubility criterion is independent of many of the details where SKiM is different from the more complex system.The results displayed in figure 1 illustrate this 'universal' nature of the stability limit.SKiM model is equivalent to a straight field line system which is, after all, a valid gyrokinetic system, it is just in a simple geometry.
Before closing this section, let us consider several details that are pertinent to particular results in figure 1.Note that the parameter T i /T e does affect the dispersion relation, equation (9), via the adiabatic electron response.It also affects the free energy equation through such a term.Since the dispersion relation pertains to eigenmodes, the eigenmodes will depend upon T i /T e .It is well known that large T i /T e is stabilizing for the ITG ae and larger T i /T e reduces the growth rate and has a strong effect upon the critical temperature gradient.But T i /T e is totally absent in the FC, equation (10), as adiabatic electrons do not carry any flux.The fact that this particular attribute of the eigenmode is absent in the insolubility criterion is a particular manifestation of the general property we emphasized earlier-the FC bound samples all fluctuations, eigenmodes or not.As the results in figure 1(e) show, this analytical FC-based approach predicts the behavior of the simulations quite well.
In a similar vein, for the SKiM dispersion relation for ETG, Debye length λ D enters as an additional term ∼λ 2 D k 2 ϕ.This term, just like large T i /T e , decreases free energy (the growth rate).However, it causes no particle flux, so the F P bound is hardly affected (figure 1(f )).
In the following sections, we will consider a more nuanced bound on the FC (derived from SKiM along similar lines), and its comparison to simulations.This comparison is in many ways more striking evidence of the validity of the approach above.And it has more practical relevance to nonlinear simulation results.But first, we consider in more detail one of the basic conceptual issues raised by the results described above.

Interpretations of FC bounds in conjunction with the simulation results: a statistical mechanical ansatz
The criterion F P > 0.6 defines the absolute limit for stability.There is nothing in the logic of the argument that implies that simulations must find stability only for F P close to 0.6.
In fact, it is easily shown by numerically solving either the 0D dispersion relation or the FC solubility criterion that many values of < k ⊥ >, < k ∥ >, and < ω d > lead to stability (FC insolubility) at F P values much less than 0.6.And this can be true by a margin of 50% or more.Similarly, simulation results for low values of k y ρ i modes attain stability for F P values that are significantly less than 0.6.
Nonetheless, the remarkable fact is that simulations, when scanned over a wide range of k y ρ i , consistently find stability when F P ∼ = 0.6.
We believe that the likely explanation for this surprising simulation behavior must be found, at least partially, in the dynamic complexity of the gyrokinetic system.A system with a large number of degrees of freedom, and a large number of eigenmodes, contrives to find some mode that has values of < k ⊥ >, < k ∥ > and < ω d > close to the ones that approach the absolute maximum F P = 0.6.And possibly, the number of possible unstable eigenmodes is particularly large when the temperature gradients are high, in other words, the system 'finds a way' to be maximally unstable.
The preceding statement ought to be viewed as a metaphor describing a class of behavior displayed in many complex systems.Such systems often manifest 'self-organization' that is widely acknowledged as behaving in an 'adaptive' manner, that is, organized structures in the system autonomously respond to changing or challenging conditions.Adaptation in complex systems is a very highly studied phenomenon, we refer the reader to e.g.[29,30].
Adaptivity equips the gyrokinetic system to find a subset of modes from the spectrum of modes (with a wide range of < k ⊥ > , < k ∥ > and < ω d >) that remain unstable until the absolute bound on F P is reached.
As we proceed with our enquiry, we will attempt to organize our findings in a somewhat systematic manner.This organization will be done under the banner of what we call a 'statistical mechanical ansatz' -its first element is: • Gyrokinetic systems have apparent adaptivity for any specified plasma conditions, there will be some class of fluctuations that will remain unstable until close to the upper bound on FC solubility.
Upper bounds on stability, therefore, have not only the virtue of simplicity but are very meaningful indicators of the overall stability of very complex systems.Like all descriptions in statistical mechanics, the current approach forgoes an attempt to describe details of the system (the specific eigenmodes, for instance), in favor of getting at the most important pertinent, general properties-the stability boundary.
In the following, we will see that simulations show that the adaptability of gyrokinetic instabilities applies to more than just the FC.It also applies to energetic effects such as stabilizing curvature.So the ansatz above can be generalized to be more inclusive.
The preceding ansatz might also explain why an approximate theory such as the SKiM model is able to predict the behavior of realistic simulations so accurately.The bounds derived from SKiM and gyrokinetic simulations agree rather well (within 10%).This is considerably more accurate than the typical level of agreement of SKiM predicted growth rates with simulation growth rates using equations ( 5)-( 7) and the dispersion relation equation (9).Thus the utility of SKiM to predict the characteristics of a specific eigenmode is limited; it, however, excels in predicting the general (independent of the nature of fluctuations, magnetic geometry, and detailed plasma parameters) phenomenon like the upper bound on FC solubility.The SKiM upper bounds are in excellent agreement with the simulation results in figure 1, calculated, for example, for complex magnetic geometries.Some additional analyses on this subject are given in appendix B.
The FC will shortly be derived and examined under more specific conditions; it will no longer be universal.Agreement with simulations, however, becomes even more striking reinforcing our faith in the validity and power of these concepts.

The flux constraint for low k modes, and implications for nonlinear simulations
We begin this section with several nonlinear simulation results for ITG ae that will motivate the need to create a framework specifically for low k modes.This will uncover arguably the most extraordinary consequence of the FC.
In figure 2(a), we plot the nonlinear heat flux (normalized to its value at F P = 0.2) versus F P .Nonlinear simulations were carried out for very diverse geometries and parameters (see also table 1), and the main results can be summarized as follows • Using Miller geometry, we vary safety factor q (1.4 ⩽ q ⩽ 4), ŝ (−0.5 ⩽ ŝ ⩽ 2), Shafranov parameter α (0.18 ⩽ α ⩽ 6), and local inverse aspect ratio ϵ (0.1 ⩽ ϵ ⩽ 0.3).• Gradients vary substantially as well (15 ⩽ R/L T ⩽ 50 and 1 ⩽ T i /T e ⩽ 3).) • There is an almost 'universal' behavior for these cases: The nonlinear heat flux for all cases falls by about two orders of magnitude for F P ∼0.5.This is significantly before F P = 0.6.• The drop in heat flux with F P is much faster than the growth rates; in the examples in figure 1, the maximum growth rates typically fell by only a factor ∼2, at F P ∼0.5, not the two orders of magnitude found in nonlinear simulations.• The nonlinear drop is very close to an analytic bound for solubility of the FC that we will derive for low k y modes.
We now turn our attention to understanding this behavior.
It is well known that nonlinear heat fluxes for ITG ae are dominated by lower k y ρ i fluctuations.The new analytic bound for such modes gives a good approximation for the F P value at which flux drops 2 orders of magnitude or more.
A reduction in the heat flux by about two orders of magnitude makes a TB possible, even if the flux is not zero.This new lower FC bound implies that a lower density gradient may be enough for a TB.With impurities, the F P bound becomes ever lower i.e, even less density gradient is enough.With experimentally relevant impurity levels, in fact, a substantial proportion of observed TBs satisfy this new bound.Hence TBs are possible, at least insofar as ITG ae is concerned, even without velocity shear; insolubility of the FC assures it.
Because of its deep relevance, the linear bound for solubility of the FC for low k modes will be derived via SKiM in section 9.In this section, for extensive comparison to simulations, we will simply read it from figure 2(b) where regions of solubility are demarcated in the < k ⊥ ρ i >-F P plot.Notice that though the maximal bound F P = 0.6 was 'universal', the new bound depends upon the eigenfunction averaged < k ⊥ ρ i >.For low values of < k ⊥ ρ i >, stability results when F P > 1/3.The maximum soluble F P ∼ = 0.53 occurs at < k ⊥ ρ i > ∼ = 1.3.(In all these cases, Z eff = 1; the impurity modifications will be described later in sections 10, 11, and 14).
The behavior of low k y ρ i eigenmodes as F P exceeds ∼1/3 is a dramatic example of adaptivity in the gyrokinetic system.Consider an example for ITB geometry in figure 3(a) The parameters are similar to a JET campaign with slightly reversed shear, and pellets to increase density gradients.It is also similar to some high β poloidal cases on DIII-D with beam fueling.The eigenmode shape dramatically changes as F P is increased, so that the eigenfunction average < k ⊥ ρ i > is increased to stay in the soluble (=unstable) range, computed analytically.The simulation eigenfunction has < k ⊥ ρ i > that clearly tracks the analytic bound.Starting from a conventional form (a 'bump' localized at the outboard midplane) at F P = 0.28 < 1/3, the eigenfunctions change, even drastically, as F P becomes larger (figure 3(b)).It has a considerable amplitude at higher θ for F P = 0.35, and ends up becoming broader and broader as we approach F P = 0.46; to remain unstable, the mode must shift to higher < k ⊥ ρ i >.Increasing F P beyond 0.46 results in stability.However, before becoming stable, the eigenfunction has undergone enormous evolution to remain in the soluble region of the FC, by increasing the eigenfunction average < k ⊥ ρ i >.

This qualitative behavior is essentially universal for ITG ae in diverse geometries and parameters. The controlling influence of FC solubility on the mode dynamics as F P is increased, emerges as the fundamental reason for the universal behavior.
One of the more spectacular consequences of the eigenfunction having a richer content in larger k ⊥ is its effect on the nonlinear heat flux.Larger k ⊥ implies a reduction in the size of turbulent eddies; a smaller eddy 'step size' would, in turn, reduce the turbulent transport.Also, increases in k increase the nonlinear coupling to stable modes in the system, and this should decrease the turbulence amplitudes, and hence the turbulent heat flux.A common dimensional estimate of the heat diffusivity from linear eigenmodes is the 'mixing length' rule ⊥ which must go down for larger < k ⊥ > The Since simulations indicate that the FC bound for low k modes has considerable nonlinear significance, let us turn to other remarkable consequences.
In the four diverse cases in figure 4, the average < k ⊥ ρ i > from all the unstable simulation eigenfunctions stays in the soluble region.There are two important aspects of these results: 1) Most obviously, the modes become stable as the analytic boundary is approached.This is qualitatively like the results in figure 4. 2) As the boundary is approached, the eigenfunction structure changes very strongly to increase < k ⊥ ρ i >.The eigenfunctions are given in figure 5, for F P = 0 and for F P very close to the stability limit found in the simulations.
The cases with F P = 0 have a very familiar structure: a bump on the outboard midplane, in the region of 'bad' curvature.This also is in a region of low k ⊥ (θ).For F P near the solubility boundary, the structure changes radically, and is concentrated in a region of much larger k ⊥ (θ).
Note that each case has eigenfunctions that differ substantially in detail near the stability limit.Nonetheless, the eigenfunction average < k ⊥ ρ i > tracks the analytic curve quite well.
Let us take stock of what we have learnt: • The eigenfunction undergoes very severe structural modifications to remain unstable.• There is absolutely no apparent reason for such extreme changes, other than to stay in the soluble region of the FC • But for the availability of the analytic solubility curve, and the eigenfunction average < k ⊥ ρ i >, these modifications would seem inexplicable and even random.• At high F P , each geometry's eigenfunctions are both quite complex and uniquely different.The only common factor among them is that their < k ⊥ ρ i > is rather near the analytic solubility boundary for the F P value.
And perhaps most important of all, the eigenfunctions are incomprehensibly complex defying understanding.Nevertheless, their behavior becomes coherent in a reduced description in terms of eigenfunction averages.This is the hallmark of a system when a statistical mechanical approach is justified.Since drastic modification of the eigenfunctions near the stability boundary has profound implications for the nonlinear heat flux, we have extended the probe by including finite k x = k y ŝη 0 ; different values of η 0 (0, π/2 and π)providing an additional degree of freedom.One could now expect that the solubility limit can be approached more closely when there are multiple η 0 to choose from.Actually, there are two contrasting regions.For F P ∼0, finite k x modes are almost always The simulation eigenfunctions [Reψ (Imψ)normalized to be unity (zero)and at the origin] for the four cases in figure 4 for F P = 0 and F P near the stability limit.The value of < k ⊥ ρ i > (θ) for the geometry is also shown (the dashed line).Near the solubility limit, the eigenfunction structure shifts to reside in regions of larger < k ⊥ ρ i >.The eigenfunction structure is often quite complex, and differs strongly from case to case.Nonetheless, the eigenfunction averaged < k ⊥ ρ i > (θ) tracks the analytic expression.Simulation eigenfunction (Reψ , Imψ ) for representative cases, for F P near the stability limit.We include cases where kx is not zero on the outer midplane, i.e. η 0 is not zero.The value of < k ⊥ ρ i > (θ) for the geometry is also shown (the dashed line).Near the solubility limit, the eigenfunction structure shifts to reside in regions of larger < k ⊥ ρ i >.The eigenfunction structure is often quite complex, and differs strongly from case to case.Nonetheless, the eigenfunction averaged < k ⊥ ρ i > (θ) tracks the analytic expression.
more stable than k x = 0.However, as the analytic FC bound in figure 3(b) is approached, we find that some modes with k x ̸ = 0 can be more unstable, and solubility limit can be approached more closely.
In the following, for each configuration and F P value, we perform simulations for η 0 = 0, π/2, and π.Among these, we plot the value of < k ⊥ ρ i > that is closest to the bound, or, below it.If it is in fact below the curve, that is, it shows a disagreement with the analytic bound.So plotting the lowest < k ⊥ ρ i > from multiple η 0 is a more stringent test of the accuracy of that bound.
Results are plotted in figure 8 for multiple cases.The minimum value of < k ⊥ ρ i > is usually distinctly closer to the boundary between soluble and insoluble regions for the FC.In a few cases, the point penetrates very slightly into the insoluble region.One must keep in mind that the analytically computed expression is an expansion to the lowest order in a parameter, so there are some corrections to it.The results strongly support the analytic bound.
In figure 8, we have labeled points as hollow circles when their D mix < 2% of the maximum D mix for that scan.Such low D mix occurs only when k ⊥ has increased by several folds.Most of those points that extend slightly into the analytic insoluble region have very low D mix , so they are quite feeble modes.
In figure 8, the eigenfunctions, just before stabilization, display enormous complexity.Representative eigenfunctions are shown in figure 6.All cases have < k ⊥ ρ i > quite close to the analytic boundary (including for k x ̸ = 0), but their mode structure is too complex to understand otherwise.
Let us revisit the nonlinear relevance of the analytic bound.Compare the trends in D mix to the trends in the nonlinear heat flux (figure 7).Between F P = 0.2 and F P = 0.53, the reductions in heat flux are a reasonably good match to the trends in D mix .The reduction in γ, however, is usually considerably weaker than in the nonlinear χ.This is evidence that the increase in < k ⊥ > is playing an important role in the reduction of the nonlinear heat flux.And as we have seen above, the increase in < k ⊥ > is due to the FC solubility boundary plus the adaptivity of the gyrokinetic system.
Finally, in figures 7(h) and (i), the variation of nonlinear heat flux with F P for T i /T e = (1, 2 and 3) is displayed.We normalize to gyroBohm units in the ion temperature.We can consider the results as a thought experiment where we keep the ion temperature and gradients constant, and only change the T e for the adiabatic electrons.The FC solubility bound is unaffected by this, but, the energetic equation is strongly affected by the stabilizing effect of higher T i /T e .This is manifest as a strong reduction in heat flux by almost an order of magnitude.Despite such a major change, the F P value at which flux is reduced by roughly two orders of magnitude is almost unchanged.This correlated perfectly with the simulation results for D mix .This is yet more evidence that the gyrokinetic system is ruled by two separate dynamics, the free energy equation and the FC, and both linear and nonlinear simulations distinguish between these quite cleanly.
The simulation results in this section, even more than previous sections, impel a statistical mechanical approach.Due to the complexity and major variations in eigenfunctions (particularly near the solubility limit of the FC), the conventional stability theory becomes both intractable and opaque.However, when analyzed as a mean field theory (dealing with statistically averaged quantities), the gyrokinetic system becomes much more transparent.The large number of degrees of freedom that allows a statistical description also imbues the system with strong adaptability; stability behavior becomes predictable-the modes become stable only near the solubility limit of the FC.
In other words, if one chooses a description in terms of eigenfunction averages, instead of eigenmodes, and posits the system with adaptability, the behavior of the system is vastly simplified, and becomes quite tractable; simulations strongly support this approach.
We, thus, add another element to the statistical mechanical ansatz:  1, along with the analytically computed solubility boundary.The magenta crosses are the < k ⊥ ρ i > for η 0 = 0, and the solid diamonds are the minimum < k ⊥ ρ i > among the three η 0 simulated (for unstable modes).Hollow circles are also the minimum among the three η 0 , but these are the cases that are only feebly unstable (D mix < 0.02 of the maximum of the scan).As expected, the addition of more degrees of freedom (multiple η 0 ) allows some modes to approach the solubility boundary more closely.Sometimes, feeble modes can exist slightly beyond the boundary, so the analytic theory is not perfect.But the overall agreement is very good.Despite feeble modes crossing the boundary, the overall agreement is very good, especially considering that approximations must be used to obtain the analytic results.
• In situations with large free energy (large gradients), the insolubility of the FC is often the operative dynamic that leads to stability.• As is observed in very many systems with a large number of degrees of freedom, the gyrokinetic system manifests apparent adaptive behavior, and 'finds a way' to dissipate free energy until a 'hard' dynamical constraint prevents it-here, the FC.
After analyzing the simulation and theoretical results in the following sections, we will add new elements to the ansatz.

Graphical visualization of the two interacting dynamics-the constraint and free energy
Exposing physics through a simultaneous graphical representation of the FC and of the free energy (like the curvature drive) is very elucidating.Consider the SKiM dispersion relation equation ( 9), which we can write as D(ω) = 0.It can be shown that Im(D(ω)) = 0 and Re(iωD(ω)) = 0 represent, respectively, the FC and the free energy equation.Each of these equations corresponds to a curve in the upper half plane (ω r , γ).We gain a great deal of insight by plotting each of these two curves; the simultaneous solution to the two is equivalent to the dispersion relation.
In figure 9, drawn for ITG ae in an ITB, we shade the region with positive particle flux (in the direction driven by the density gradient) orange.The region of negative particle flux is white.The red boundary between these is where the FC is satisfied.The free energy balance is the blue curve.
As F P is increased, the FC curve undergoes a profound change.As expected, the region of positive particle flux grows, since a greater particle flux is driven by the density gradient (thermodynamic forces drive their respective thermodynamic fluxes).Also, the region of negative particle flux shrinks.Hence, the curve where the FC is satisfied becomes smaller and smaller.Eventually, for large enough F P , the entire upper half plane has positive flux: the FC has no solution, and there is no instability.
Notice that the free energy curve does not change qualitatively during this process.Stability for higher F P comes because the red curve shrinks(and disappears); the solution will have progressively lower γ, until no positive value γ solves the FC.
For a very different sequence, consider progressively decreasing the curvature drift, going from positive to negative values.This will also stabilize the mode.See figure 10.We keep F P = 0.2, well away from the solubility limit.As we see, the free energy curve shrinks until it has a solution only near γ∼0.The FC, on the other hand, always has a considerable extent at high γ.The reason for stability in this case is clearly that the free energy balance equation cannot be satisfied for strongly negative (stabilizing) curvature.
Hence, these plots make it very clear that stability can result from two very different physical origins.This is the essence of the dual dynamics approach.These same types of plots will be very clarifying for cases with non-adiabatic electrons.
We now turn to the stabilizing effects of curvature, i.e. energetic stabilization.

Curvature and adaptability for ITG ae
Configurations favorable for TB formation, often, have strongly stabilizing curvature at most positions of the poloidal angle.This is well known [31] for tokamaks, especially, where TBs are often found under conditions of negative magnetic shear or strong Shafranov shift.In light of the eigenfunction adaptability, we must examine if such energetic effects are sufficient for stability or if some help from the constraint will still be necessary.
We observed in section 6 the remarkably adaptive behavior of ITG ae : as the FC constraint tends towards insolubility, the mode structure changes, often radically, to remain unstable by staying in the soluble region of the FC.Response to curvature stabilization turns out to be similarly adaptivesimulations show that ITG ae mode 'fights' the influence of stabilizing curvature by residing, preferentially, in the 'bad' curvature region.Such an adaptation circumvents curvature stabilization, and the role of FC dynamics for stabilization will remain crucial; it will persist even when non-adiabatic electron effects are included in the ITG/TEM because mode adaptation to prevent stabilization by curvature will still pertain.We first consider an equilibrium sequence similar to the ITB case A (with parameters given in table 2), but make the magnetic shear much more negative.Recall that case A is roughly representative of a JET campaign with controlled density profiles in the ITB region using pellet injection.Also, quite similar parameters are often attained in high β poloidal cases on DIII-D.We consider these parameters and also make the shear more negative yet.Cases #3 and #4 correspond, roughly, to the actual experimental ITBs on JET and many cases on DIII-D.Cases #5 and #6 have much higher negative shear.However, they are roughly representative of some other experimental discharges with strongly reversed shear (e.g.NCS shots on DIII-D and strongly reversed shear shots on JT60U).For completeness, cases #1 and #2 are included and are roughly representative of situations before an ITB, where α is not high and ŝ is not very low.
The measure of the resulting curvature drive in the gyrokinetic equation -ω d -is plotted as a function of poloidal angle θ in figure 11(a) (for the standard ballooning coordinate description).Positive values are destabilizing and negative values are stabilizing.For experimental ITBs on JET, the curvature is very preponderantly negative.Since shear is even more negative, it becomes extraordinarily stabilizing, except for a small region of 'bad' curvature around the outboard mid-plane.
The parametric scan of the simulations shows a distinct separation between energetic and FC dynamics.The D mix (figure 11(b) is reduced an order of magnitude by the energetic effects (from #2 to #6), but the stability point in F P is hardly affected.
The plots in figure 11 are also an example of the strong adaptability of these modes to avoid outright stabilization by overwhelming 'good' curvature.We turn to this aspect now and use the eigenfunction average curvature, equation (7).
Despite the enormously stabilizing curvature over most of θ, the eigenfunction, by strongly concentrating in the bad curvature region, still remains unstable; the eigenfunctionaveraged curvature equation (7) remains positive, although it does decrease somewhat (figure 11(c).The contrast in figures 11(a) and (c) is quite amazing!Without such adaptability, the hugely negative curvature would have directly stabilized the modes, and density gradients would have been unnecessary.We also show in figure 11(c) what would happen if the mode had not adapted.We do this via a thought experiment: we compute the average curvature for different equilibria but using the eigenfunction from case #1.In this hypothetical situation, the average curvature would drop to sharply negative values.Using the SKiM dispersion relation and typical values for the other parameters, we plot the growth rate for various average curvatures.Sufficiently negative eigenfunction averaged curvatures stabilize the mode, as seen in figure 11(d), for all F P .Such values would have arisen if the mode structure had not adapted to stay predominantly in the bad curvature region.
Notice that SKiM finds that somewhat before the point of complete stabilization by curvature, the critical F P is unrelated to the constraint prediction: one stability point occurs at low F P ∼0, and it also occurs for other F P that are much less than the solubility boundary of the constraint.Depending upon the particular value of negative average curvature, the critical F P can be anything from 0 up to the constraint value.In other words, curvature stabilization is through the free energy equation, and so the stability point is unrelated to the FC solubility bound.
Simulations find that stabilization independent of F P rarely happens for TB level gradients (R/L T > 15 or so).The simulation results in figure 11 are typical of many cases we have examined.Because of adaptability to stabilizing curvature, only the FC can give stability for most ITB parameters.Hence, (sufficient) density gradients are usually a crucial ingredient for suppressing ITG modes in ITBs, especially for those with sharp gradients.This is, at least qualitatively, in full agreement with the JET experimental results for cases with low velocity shear; density gradients were, indeed, crucial to the ITB in the experiment.In fact, simulations find that the crucial role of density gradients persists even for much more negative shear than is typically found in experiments.
To emphasize: the mode adapts to avoid curvature stabilization; it is only stabilized by the FC.
As F P increases, eigenmode structure in TBs must adapt to both the curvature and to the FC solubility to remain unstable.
Interestingly, as the FC solubility boundary is approached, the mode must undergo two metamorphisms simultaneously: It must increase k ⊥ to remain unstable and at the same time strongly peak in whatever bad curvature region that is available.The mode tries to do exactly that, as shown in simulations figure 12.The eigenfunction adapts to have positive eigenfunction average curvature except for cases with truly extreme curvature and large F P near the solubility limit.This simultaneous combination exceeds the adaptive capacity of the system.At larger F P , the region of destabilizing curvature region is not near θ = 0, but is located at θ > π (this region is not shown in figure 11(a)).The eigenfunction average curvature stays positive even as the mode structure changes to increase k ⊥ , except very near the solubility limit for the cases with the most extreme curvature (case #5 and case#6).
In the supplement, even a more extreme equilibrium sequence is simulated, with α reaching extremely high values rarely seen experimentally.At α∼21, the mode becomes predominantly curvature stabilized, with a growth rate dependence on F P qualitatively like to figure 11(a).Other aspects of this sequence confirm the results of this section.

Derivation of the flux constraint for low k modes
Examining equation (10), the mode number appears in two places: through < ω di > and through the arguments of Bessel function J 0 .In the limit of vanishing k, the FC becomes insoluble for F P > 1/3 (it differs from the limiting value F P > 0.6 due to the contribution of the second term in the square brackets).However, simulations show: 1) low k instabilities regularly exceed F P = 1/3, and 2) instabilities occur at < k ⊥ > that cannot be considered small.But almost always, < ω di > remains small compared to < k ∥ v th >.We therefore analytically derive a new FC bound by taking < ω di > small compared to < k ∥ v th >, but allow < k ⊥ > ρ i to have any value.By neglecting < ω di > in the denominator, one can integrate the resonant denominator in v ∥ for low γ.We are particularly interested in the region near marginal stability, where a delta function approximation to the resonant denominator suffices.The integral in v ∥ becomes trivial, and equation (10) becomes Examining equation (11), the value of ω r that allows the largest solubility domain in 1/L n is for ω r = 0. (This is where the temperature gradient term is minimum.Note that the exponential term factors out.)Hence ω r = 0 sets the solubility limit, and the boundary is given by, and k ⊥norm = k ⊥ v thi /Ω i .This was used to obtain the analytic curve in the previous section.Note that equation (12), an approximation of equation (10), is not as accurate as the original expression.And indeed, simulations do find points that lie slightly over the boundary in the nominally insoluble regime.One should note various ways the approximation could break down: • Even if the expansion is fairly accurate for the FC bound, there will be next order corrections, so that the boundary can be somewhat shifted from the computed value • The expression for < k ⊥ > is intuitive and can be derived in the limit of small < k ⊥ >, but it is only an approximation; there could be corrections.• Trapped ions are not included in the kinetic treatment, and as ω and γ become small, these could affect the results.• Because of any combination of these effects, even if the overall theory is reasonably accurate, residual, weak instabilities are potentially possible beyond the computed solubility bound.

The effect of impurities on the solubility of the flux constraint for low k modes
Since our FC demands the charge-weighted turbulent flux to be zero, impurities, because of their higher charge, could have a disproportionately large effect on the solubility condition.We will, therefore, generalize equation ( 12) to include impurity species treating the impurity terms exactly in the same way as the main ions.Equation ( 12) modifies to ) where v ths is the thermal speed of species s and ρ s = v ths /Ω s .Graphs of the computed solubility bound for sample parameters are shown in figure 13, where ρ i is the main ion gyroradius.We see that impurities shift the solubility boundary to lower F P .Many experimental cases for pedestals and ITBs find that the impurity gradients 1/L n are steeper than the electron values [32,33].This shifts the curve to even lower F P , i.e. it is a stabilizing effect (in the opposite case where the impurity gradient is more shallow than the electron value, this shifts the curve to higher F P values).
For common impurities and relevant impurity levels, e.g.Z eff = 2 ∼ 3, a significantly lower density gradient is needed to stabilize ITG modes.This brings significantly more experimental TB into the range where the ITG is stable without velocity shear.
Before we show, through simulations, the heightened importance of impurity stabilization, let us contrast the true operative physics here with the usual invocation regarding impurities: the so-called dilution.
In the context of TBs, conventional approaches attribute stabilization by impurities to dilution.In fact, for the free energy equation and curvature-driven modes, stabilization via dilution is a reasonable conclusion because it can be shown that higher Z impurities contribute much less free energy than the bulk ions they displace.The effect is then proportional to the impurity to ion density fraction.
Impurities affect the FC strongly, way more than what could be explained by the effects of dilution.Quite unlike dilution, the shift in F P (figure 13) arises because impurities contribute a disproportionately large charge flux.If the impurities had been inert, then the FC boundary would be unchanged, since in equation ( 12), the numerator and denominator would be reduced by the same impurity dilution factor.Impurity dilution, then, would result in the same solubility boundary as without impurities.Instead, quite the opposite happens: the impurities contribute more strongly because their denominator is smaller, since their velocities are smaller (by ∼ √ m bulk /m impurity ).Hence, they give a larger contribution to the FC relative to their charge fraction in the equilibriumexactly opposite to how they affect the free energy equation.
Hence, for ITG stabilization by density gradients, impurities are not a dilution effect at all; they stabilize by substantial contributions to charge flux.
Note that it may be straightforward to model energetic particles as very high temperature ionic species in equation (13).And with appropriate gradients of density or temperature, they, like impurities, could be stabilizing through their contributions to ionic charge flux making it harder for FC to be satisfied.Recent investigations have demonstrated that energetic particles can have a stabilizing effect [34][35][36].However, this interesting deviation is beyond the scope of this already crowded enquiry.We will limit ourselves to only the thermal species.
We now show how strongly the simulations agree with the analytically computed bounds for the FC with impurities.

Impurity effects on FC-theory and simulations
In this section, we find that simulations corroborate the analytical predictions: when impurities are included, the (main plasma) density gradient needed to stabilize the ITG ae is reduced.For experimental parameters, this reduction can be up to a factor of ∼2.This brings a far larger number of experimental TBs into the regime where the ITG ae is stable without velocity shear.And in TBs, it is frequently the case that impurity gradient scale lengths are shorter than the electron density scale length further reducing the needed density gradient.
These results, on impurities being facilitators of TB formation, could be of great practical importance.They, indeed, constitute important theoretical/scientific evidence of the fundamental processes that control such formation.The analytically predicted solubility bounds on the FC are found to agree with simulations fortifying the validity of the FC-based approach.
In figure 14, simulations, with diverse geometries and temperature gradients, track the analytic bounds, except that quite weak instabilities (with D mix < 2% of the maximum value for the scan, denoted by an open circle) persist into the nominally insoluble region.As for the plots in previous sections, three values of the k x were simulated (k x = k y ŝη 0 ), for η 0 = 0, π/2 and π, and the lowest value of < k ⊥ > ρ i (which is the one that most challenges the analytic boundary).Simulations also track the analytic bounds predicted for different impurity levels and impurity gradients, with the same parameters as in figure 13.
It is striking that, despite tremendous changes in geometry and temperature gradient, the behavior of all but these very weak modes respects the FC bound of low k modes with the impurity modifications, validating the analytical theory.Consistent with the theory, the impurities cause stability for a lower F P , and this tracks the analytic curve: in the insoluble domain, only feeble modes persist (D mix < 2% of the maximum value).
These residual instabilities have little impact on the behavior of D mix for the low k modes, as seen in figures 15(a) and (b).Despite the tremendous differences in geometry and gradients R/L T , there is a dramatic and rapid reduction in D mix near the analytic solubility bound for the FC.Specifically, the cases with impurities drop at the analytically predicted F P value for their respective impurity parameters, validating the analytic bounds.The differences in parameters are manifest in somewhat different trajectories, but the endpoint of stability (or near stability) is consistent with the predictions.
Nonlinear simulations of the heat flux follow the same behavior as D mix (see figures 15(c) and (d)), and the flux drops by about two orders of magnitude at the analytically predicted F P bound.This occurs for both ITB-like and pedestal configurations with much steeper gradients.
We examine the effects of T i /T e in figures 15(e) and (f ).We compare the nonlinear fluxes for cases with T i /T e = 2 and 3 and Z eff = 1 and the ones with impurities (fluxes are in gyroBohm normalized units in the ion temperature).Recall that T i /T e > 1 does not affect the FC bound, but it does have a stabilizing effect in the free energy equation.The cases with T i /T e = 2 and 3 have far lower flux than the case with Z eff = 1 and T i /T e = 1 at low F P .This is the expected energetic effect.But the FC solubility limit for all these cases is the same, so they all drop to low values at the same F P .Now, let us turn on the impurities.For low F P = 0.2 < 1/3, the FC is not very restrictive.The impurity cases have a larger flux than the cases T i /T e = 2 and 3.But as F P increases, this reverses.The presence of impurities implies a lower FC bound, and so the flux drops below the cases with T i /T e > 1.Once again, the FC dominates the behavior at large F P near the solubility limit and is more important than energetic effects.
Our methods (a combination of analytics and simulations) again confirm the power of the FC dynamics and show that the impurities can play an outsize, even, crucial role in TB formation.

Insolubility of the FC and statistical mechanics-generalization to non-adiabatic electrons
The main results from our exploration of the theoretical structure of the gyrokinetic system, are strongly supported by extensive gyrokinetic simulations (using the GENE code [19]).Let us step back and examine the conceptual foundations.
Since FC is derived for essentially all perturbations (and not just eigenmodes), its origin and its operational manifestation must be anchored in physics concepts more general than that guide the conventional eigenmode stability analysis.Since the FC is a (charged) particle FC, it naturally valorizes the density gradients.Let us consider this in more detail.
In the so-called Onsager region, the thermodynamic forces are linearly related to the thermodynamic fluxes.It is well known that this classic region does not, generally, pertain to transport from gyrokinetic instabilites.But for the simpler problem of calculating the charge flux for a given fluctuation, methods can be used that are completely analogous to the Onsager region.Let us try to understand the FC solubility through such concepts.
In the classic Onsager theory, the gradient dQ i /dx of a physical quantity Q i represents a thermodynamic force.The flux Γ i of Q i is given by where the matrix M ij often has off-diagonal components.(i.e.temperature gradients could drive particle flux, for instance).The second law of thermodynamics demands that this matrix must have positive definite eigenvalues.For our purposes, the important corollary is that the diagonal components must be positive-thermodynamic forces drive their diagonal fluxes down the gradients.
In the linear gyrokinetic equation, the non-adiabatic part of the distribution function is driven by the gradients dQ i /dx; there is also an additional term proportional to frequency.And crucially, in the context of showing insolubility of the FC to general fluctuations, that frequency can be taken to be arbitrary.In other words, it is not considered to be an eigenfrequency, and the fluctuation structure also need not be an eigenmode.
In the presence of a general fluctuating potential ϕ, the induced quasilinear flux, then, has an additional term (due to the finite frequency) modifying equation ( 14) to, For example, in the simple ITG ae with a single species contributing, the ion particle flux will be In general, A and B have similar magnitudes but typically A < 0, and necessarily B > 0. So, for some density gradient that is large enough, the ion particle flux will always be positive, and the FC is insoluble.This insolubility is unrelated to the amount of free energy in the equilibrium, which is to say, how strong the gradient is.The requisite (stabilizing) density gradient is determined by the temperature gradient; the solubility condition, thus, can be expressed in terms of a ratio like η or The fact that diagonal terms (like B) in the response matrix are positive definite is a consequence of entropy considerations, qualitatively like the classic case, but different in some important details.For example, the classic considerations do not apply to highly non-equilibrium situations like perturbations that are increasing exponentially, or, where collisions are zero so there is no entropy production.Nonetheless, it can be shown to be true that the diagonal components of the gyrokinetic response matrix G are positive, even for exponentiating fluctuations and without collisions.For this, we use the free energy equation equation ( 2), which describes entropy considerations for the gyrokinetic system.See appendix A.
The generality of the preceding arguments transcends the simplified model ITG ae .For example, even with non-adiabatic electrons, and electromagnetic fluctuations, the matrix Q ij will have qualitatively the same properties.
The FC might be readily satisfied by having the ion and electron charge fluxes cancel each other.That is, even though the equilibrium density gradient will drive a particle flux of each species, the net charge flux will vanish, and the FC is soluble.So in the general case, there need not be a regime for TBs.
However, if one of the species (say electrons) is sufficiently close to adiabatic (causing very little contribution to the charge flux) and the ionic species is not, then a large enough density gradient will lead to a net charge flux, i.e. the FC will be insoluble.
To the best of our knowledge, no one has recognized/investigated that electrons tend to be near adiabatic for the instabilities that must be quenched for TB formation.We will do this in the following sections.
The next section will begin by modifying the FC when electrons are non-adiabatic.
It is time to add an additional item to the evolving statistical mechanical ansatz: • In situations with large free energy (large gradients), the insolubility of the FC is often the operative dynamic that leads to stability.• As is observed in many systems with a large number of degrees of freedom, the gyrokinetic system manifests apparent adaptive behavior, and 'finds a way' to dissipate free energy until a 'hard' dynamical constraint prevents it-here, the FC.This apparent adaptivity also applies to the curvature drive-the eigenmode structure avoids regions of stabilizing curvature, even when these are strongly predominant • The basic trends of stability due to the FC can be understood from the very general statistical mechanical propensity for thermodynamic forces to drive the corresponding (diagonal) thermodynamic fluxes down the gradient.When only one of the species is sufficiently close to adiabatic, this propensity eventually makes the FC insoluble as density gradients are increased.

The flux constraint including non-adiabatic electrons
In the preceding section we discussed how the basic concepts behind the working of FC would generalize to non-adiabatic electrons.Here we consider the full range of pertinent dynamics in detail.Non-adiabatic electron response often plays a crucial role in gyrokinetic instabilities.For the ITG/TEM, this arises primarily from trapped electrons, e.g. it can lead to the TEM.But passing electrons can also play a role.Since non-adiabatic electrons contribute charge flux in the opposite direction, the ion (plus impurities) charge flux will effectively diminish; violation of FC, therefore, will become more difficult.We will see, however, that even with non-adiabatic electrons, the FC may still be potent enough to suppress instabilities and allow TB formation.We find that FC violation can still take place in diverse parameter regimes that are similar to experimental TBs.Electron contribution, however, modifies the solubility conditions.
It is possible that when the non-adiabatic electron charge flux is large enough, density gradients do not, necessarily, lead to insolubility.This situation occurs for many typical tokamak geometries, in fact.In such geometries, the stabilization seen in the previous sections is either greatly reduced or eliminated because the FC remains soluble for much larger density gradients.
Within the conventional framework of dispersion relations, this is interpreted as the onset of electron instabilities, e.g. the TEM that can access new free energy from density gradients.
However, the dominant factor is not so much an increase in the free energy as that the constraint becomes soluble.We show this for the two non-TB equilibria in table 1 (case non-ITB A and B).We can include trapped electrons in the SKiM model.Details are given below; here let us discuss the conceptually important results.We use simulation eigenfunctions to evaluate the parameters in SKiM.We plot the E-FC plots for each case at F P = 0.5 in figure 16.We compare to the plots with the trapped electron terms neglected-the latter being ITG ae .Although the trapped electrons do indeed quantitatively increase the free energy significantly, the qualitative difference is in the FC.Without trapped electrons, the FC is insoluble except for γ∼0.With trapped electrons, the FC is soluble in a large region, up to high γ.A strong instability, then, becomes possible.
In fact, even if the free energy curve were the same as the ITG ae case, the fact that the constraint was rendered soluble by non-adiabatic electrons would have produced an instability with a high growth rate.
When the FC does not constrain the system to a low growth, there is nothing to prevent the relaxation of free energy.Then Figure 16.E-FC plot for cases trapped electrons compared to adiabatic electrons.The red curve is again the solution of the FC, and the blue curve is the free energy balance.The eigenvalue lies at the intersection.The electron flux makes the FC soluble for high growth rates, when it is insoluble for adiabatic electrons.This allows a strong instability for the full electron cases.It also makes the free energy curve go to higher γ, but the qualitative change is to the FC.The free energy curves indicate that the ITG drive contributes a substantial fraction of the ITG/TEM free energy.In fact, even without the boost in free energy from trapped particles, there would still be a substantial instability just due to the FC curve modified by trapped electrons.
the gyrokinetic system behaves like the vast majority of physical systems: instabilities arise to relax the free energy.Strong gradients will drive strong hybrid instabilities of ITG and TEM, whatever name one might wish to call it.
Irrespective of the details of the instability, we expect that a highly adaptive system like gyrokinetics will produce a strong instability once the FC becomes soluble.For example, the ITG ae results showed an extreme adaptability of the eigenmode structure in order to produce an instability that satisfies the FC and avoids stabilizing curvature.With such adaptability, the inclusion of non-adiabatic electrons that make the FC soluble at high γ should be expected to lead to strong instability, whatever the specific details.Just as for the ITG ae , we expect that the gyrokinetic system will display extreme adaptability to avoid curvature (energetic) stabilization, by changing its mode structure to avoid stabilizing curvature.We find that this is, in fact, the case, when trapped electrons are added.
In order to restore the conditions for TB formation, then, we must probe into mechanisms/conditions that would still violate FC even when the non-adiabatic electron contribution of negative charge flux is pushing towards solubility.How do we do this?
We naturally turn to the density gradients and the FC constraint to ensure stability so that TB formation (the driving practical aim of this paper) may take place in regimes of weak velocity shear.In fact, such stable regimes have been discovered well before this work, in experiments and in simulations of tokamak geometries with low or negative magnetic shear and/or high Shafranov shift.We will have to provide a framework to understand how when the entire electron dynamics is included.It is not an exaggeration to state that it is only by resorting to the constrained dynamics that we may make sense of the surprising behavior found in the simulations.Conventional interpretations are, at best, inadequate.
We will now go past ITG ae and investigate the gyrokinetic system with non-adiabatic electrons.We expect the generality of the concepts that allowed us to explore the ITG ae , will help us understand TB formation with full electron dynamics; the latter is not expected to vitiate the basic statistical mechanical principles.
We reemphasize that the basic intellectual drive for this study comes from a desire to create a deep and detailed framework to understand and interpret this amazing phenomenon of TBs.Equally important is the corollary that this understanding could guide us to access even more promising TB formations.
Anticipating results in this section, the scope of this effort could be summarized as: • Understand the fundamental origin of TB formation from processes unrelated to velocity shear • Show that the simulation behavior in such cases is better explained by a collection of insights that we call the statistical mechanical ansatz • The simulation behavior has the characteristics of stabilization due to the FC • The energetic stabilization-free energy considerations do not play as dominant a role as in conventional stability theory • The realization that stabilization is due to the FC implies different mode behavior and different requirements for stabilization as compared to energetic stabilization • In particular, density gradients are a crucial aspect of the stabilization.Impurities and impurity gradients are just as relevant; they can lead to TBs in unexpected regimes.
The FC dynamics valorizes the density gradients as a stabilizing force.In the framework of conventional stability theory, this is a bit odd since density gradients are just as much a source of free energy as, for instance, temperature gradients.It could happen that some particular class of instabilities may not access their free energy, but it is quite inconceivable they emerge as stabilizers of all fluctuations when they are sufficiently large-that is, they prevent access to all other sources of free energy.This is probably the most profound scientific phenomenon that this study uncovers and dwells on.
Stabilization by energetics, on the other hand, requires extreme geometries (at least for the tokamak).These tend to be very challenging from the viewpoint of global ideal MHD stability and other practical considerations.This is, therefore, the first practical advantage that the new understanding brings-with sufficiently strong density gradients TBs could be formed in less extreme and challenging magnetic geometries; the range of accessibility could be vastly increased.
In the stellarator context, a practical consequence of this understanding is to develop better informed optimization strategies for the magnetic geometry, and, to give a different interpretation than has sometimes been adopted.
An additional practical implication is that stabilization due to the FC means that density gradients and/or impurity gradients can make it easier to initiate TBs from a geometry that did not initially have one.
We turn now to simulation results for tokamaks.This displays the characteristics of the statistical mechanical ansatz with non-adiabatic electrons quite well (the stellarator has similar general considerations, but will not be considered here due to space).

Simulation results with non-adiabatic electrons
For quite some time, it has been realized that tokamak geometries with low or negative magnetic shear and high Shafranov shift lead to weak ITG/TEM, and hence, TBs can result.By relating this class of TB formation to the general dynamics delineated above, we can recognize and explain various qualitative trends for the first time.These trends are often directly relevant to experiments.
To review, four pivotal, general behaviors uncovered in previous sections, were: • In situations with large free energy (large gradients), the insolubility of the FC is often the operative dynamic that leads to stability.• The basic trends of stability due to the FC can be understood from the very general statistical mechanical propensity for thermodynamic forces to drive the corresponding thermodynamic fluxes.In cases where one species is close to adiabatic, this propensity eventually makes the FC insoluble as density gradients are increased.• As is observed in very many systems with a large number of degrees of freedom, the gyrokinetic system manifests apparent adaptive behavior, and 'finds a way' to dissipate free energy until a 'hard' dynamical constraint prevents it-here, the FC.• The apparent adaptivity was explored in the context of the curvature drive that could be stabilizing or destabilizing.We found that in order to stay unstable, the eigenmode structure avoids regions of stabilizing curvature, even when these are strongly predominant.
Let us now view the results of electromagnetic gyrokinetic simulations with full electron dynamics in diverse geometries and equilibria.(Note that beta values consistent with α are used.This has been found to suppress passing electron effects for the ITG/TEM.We have analyzed this extensively, but due to length, will discuss the details elsewhere.Summarizing: at zero beta, non-adiabatic passing electrons can sometimes lead to a considerably enlarged instability domain, which can be understood in terms of the FC.But even low values of beta, and especially realistic ones, suffice to reduce these so that trapped electron effects dominate, as discussed below.) In figure 17, we plot D mix as a function of F P for two different equilibria (for ITB and pedestal gradients) doing a parameter scan in ŝ and α (all these cases have electron collisionality ν * ∼0.05-0.1;higher collisionality produces a similar result but stability is reached for a less extreme ŝ or α).
The figures show that an initial increment in ŝ or α brings the mode close to F P stability point of the ITG ae (the vertical dashed line).Further, much larger increases do not significantly improve the stability point in F P .As can be seen in figure 18, the curvature structure becomes extraordinarily heavily weighted toward stabilizing curvature (except for a small region of destabilizing curvature) and continues to become more stabilizing.
But eventually, the huge preponderance of good curvature hardly improves the stability in F P In all cases, the requisite F P to cause stability approaches the solubility bound for the ITG ae .This behavior is roughly generic.How does this result square with the prevailing picture before now?18.In all cases, increases in ŝ or α eventually reduce the threshold F P needed for very low D mix down to the limit of the ITGae, but even large increases beyond that do not materially reduce that F P .
The prevailing point of view would, roughly, go like this: the large region of stabilizing curvature of TB geometries (negative ŝ or large Shafranov shift for tokamaks) leads to an orbit-averaged stabilizing curvature for the trapped electrons.This stabilizes the TEM mode.Then, some combination of density gradients and stabilizing curvature stabilizes the ITG.We'll call this the 'conventional' picture of ITG/TEM stability in a TB.
This picture is inconsistent with the displayed simulation results.To see why, we must first recognize the fact that the ITG and TEM are inevitably tightly coupled together.The ITG mode produces an electrostatic potential fluctuation, and the trapped electrons react to this with a bounce average response in exactly the same way as for a 'pure' TEM, including the bounce averaged curvature effects.And similarly, a TEM mode produces E × B fluctuations that convect the ion temperature, and this couples that free energy into the fluctuation, just as with an ITG ae .So the ITG and TEM are not like two 'elementary particles', each with an independent existence, rather, each is more like a quark in the sense that it is only found in nature bound to another quark.The ITG/TEM is nearly always a composite mode (though for various parameters one aspect might be predominant).
Consequently, an extraordinarily stabilizing curvature for a TEM should affect the energetics of the ITG/TEM similarly.If curvature stabilizes TEM, ITG/TEM should be stabilized too-consequently it should be more stable than the ITG alone; and this relative stabilization would increase as the curvatures become more extreme.Simulation results are completely the contrary: the incremental stabilizing effect on the critical F P gets weaker as the curvature becomes more extreme.In fact the F P stability point approaches the adiabatic ITG ae most closely when the stabilizing trapped electron curvature is overwhelmingly large, even when it is far larger than the destabilizing curvature in the tiny region near θ = 0.
Fortunately, the new dynamic explored in this paper does explain the simulations: just like what we found in ITG ae , the Figure 18.ω d (θ) for the cases in figure 17.The curvatures become very strongly preponderantly negative, but this does not affect the F P bound found in simulations in figure 17. mode structure adapts, to avoid interacting with the stabilizing regions of phase space.Specifically, its spatial structure 'decouples' from the trapped electron response for particles with highly stabilizing orbit-averaged curvature, so it does not interact much with such particles.The mode might only couple strongly with the minority fraction of trapped electrons that have destabilizing orbit average curvature.Because there is little coupling to most trapped electrons, the 'effective' trapped electron fraction is considerably smaller than the actual fraction of trapped electrons.So the non-adiabatic electron response is small.Hence, the electrons produce little charge flux, pushing the F P limit for ITG/TEM close to that of ITG ae .
The FC based stabilization is qualitatively different from the conventional stability analysis controlled by energetic factors; consequently, different parametric behaviors will be predicted.FC predicts that increasing density gradients lead to stabilization; whereas we will see that the conventional picture of curvature stabilization requires extreme geometries, and in this case, density gradients play an insignificant role.
In short, the essential physics of stabilization (as revealed in the FC context) is the same as that of ITG ae ; note that the same behavior vis-a-vis curvature was found in section 8.

Analysis using SKiM including trapped electrons-connectivity of ITG and TEM
When trapped electrons are included, the dispersion relation in the mean field theory SKiM, has an additional term on the right hand side of equation (17).
Qualitatively, this term is familiar from early investigations of the ITG and TEM modes.It contains an integral over energy, curvature resonances, and an effective collision term.
In SKiM, two new quantities are introduced for trapped electrons, which are defined and discussed in appendix B: 1) an eigenmode average of the orbit averaged curvature < ω d orbit e >, appearing in the resonant denominator.It is similar to the ions, but the weighting involves orbit averages of ϕ rather than ϕ itself.And 2) an overall strength of the coupling to the trapped particles, the 'effective' trapped electron fraction < f trap >.The expression, independent of curvature, depends only upon the magnitude of orbit averages of ϕ (relative to ϕ itself).It describes the eigenfunction averaged strength of coupling to trapped particles, irrespective of whether the trapped response is stabilizing or destabilizing as regards curvature.For some eigenmode structures, it can be considerably less than the actual fraction of trapped particles in the geometry.
We now use equation (17) to show the qualitative difference in the predictions of the conventional picture and the one guided by the statistical mechanical ansatz.
Let us begin by examining the notion that the stabilization of the TEM can explain the behavior of the ITG/TEM; we will find that it is definitely not so!Reasonable parameters can allow the ITG/TEM to be quite unstable even though the TEM is stabilized by curvature.This is because the non-adiabatic electron response allows the FC to be satisfied for large F P , by having an electron charge flux to balance the ions.This electron flux requires an expenditure of free energy, however, due to the stabilizing orbit averaged curvature.We'll see that the ITG drive can supply more than enough energy for this and have enough left over to attain a large growth rate.
Consider a representative example.In figure 19, we consider typical parameters corresponding to the convention mode of thinking, where the average trapped electron curvature is quite stabilizing.
We compare the stability of the TEM, the ITG ae , and the ITG/TEM.We display two kinds of graphs (all with R/L Te = 20, as in a TB): figure 19(a) is a γ − F p plot while the other three are E-FC plots in the ω − γ plane.Recall these plots from section 7: they show the curve where the FC is satisfied (red), and the curve where the free energy balance is satisfied (blue).The orange region is where the FC is not satisfied because of too much ion flux, the white region is where the FC is not satisfied because of the opposite situation, and the boundary between them where the charge flux vanishes is the red curve.
There are three separate regimes to study: 1.The TEM-specific regime: set the ion temperature gradient to zero and keep only density and electron temperature gradients.Now as in the conventional scenario: we adjust the electron curvature to be sufficiently negative (TEM is always stable for all density gradients, we use R/L Te = 20).We also assume a strong coupling to trapped electrons: < f trap > is a typical physical fraction of trapped particles.2. ITG ae specific regime: set < f trap >= 0 and use the same R/L Ti = 20.3. ITG/TEM hybrid regime: combine the ITG and TEM, including the ion temperature gradient in the TEM case with trapped electrons (R/L Te = R/L Ti = 20).
All these cases are run at the same k ⊥ and k ∥ , at representative values typical of the eigenfunctions above (k ⊥ ρ i ∼0.5 and k ∥ = 0.71 and < ω di = 0.24).
Growth rates from equation ( 17) are shown in figure 19.As expected, the ITG ae is stabilized by high values of F P .The TEM is always stable.But the ITG/TEM is far more unstable than the ITG ae , up to F P values much greater than the ITG ae .
Stabilizing the (pure) TEM does not cause the ITG/TEM to be stable; it remains much more unstable than ITG ae for all F P .
This surprising behavior becomes comprehensible when we contrast the effects of trapped particles on the FC, and on free energy dynamics.We plot the E-FC curves for these three cases at F P = 0.5.The ITG ae (figure 19(b)) is stable (all orange) because no point in the upper half plane satisfies the FC with adiabatic electrons.But there is considerable free energy, i.e, growth rates would be quite large due to free energy balance alone.The TEM (figure 19(c)) is stable primarily because the free energy is so low: the stabilizing curvature has so much reduced it that balance is only achieved for low growth rates in rather tiny regions.These regions are so shrunken that they do not intersect the FC for γ > 0, so there is no instability.And distinct from ITG ae , the FC could, in principle, be satisfied up to much higher growth; the electron flux due to the non-adiabatic electron response can easily balance the ion charge flux; stability is due to the lack of free energy Now consider the ITG/TEM in the last panel.The trapped electrons (the TEM part) provide the particle flux to balance the ion flux (to satisfy the FC) and ITG part supplies the freeenergy drive.Thus large γ is allowed even at F P = 0.5.The intersection of the red and blue curves yields the unstable eigenfrequency of the dispersion relation.
Note that the trapped electron response does, indeed, reduce the free energy available to the perturbation: the free energy (blue) curve in the ITG/TEM is limited to lower γ compared to that of the ITG alone and by a large margin for ω r ≲ 1.
Figure 19.The Energy-FC plot in the complex ω (ωr-γ) plane.The fallacy of understanding the ITG/TEM as decoupled ITG and TEM (a) γ in SkiM for representative parameters with large trapped particle fraction < ftrap > and sufficient negative average electron curvature < ω d orbit e > to stabilize the pure TEM mode (with ion temperature gradient equal zero.)However, the coupled ITG/TEM (with R/L Ti = R/L Te ) is not stabilized by F P even though the ITG and TEM is stable (the ITGae is stabilized by F P as usual).The reason for this is revealed by the E-FC plots for F P = 0.5 (b) the ITGae is stabilized by the constraint as always (c) quite the contrary, the TEM is stabilized mainly energetically: the free energy curve is so close to γ = 0 that it does not intersect the FC curve, so, there is no unstable curve.However the FC alone would allow a substantial growth rate.For the ITG/TEM in (d) the free energy curve, now supplemented by the free energy of ion temperature gradients, allows high growth rates, unlike the TEM.The FC curve also allows high growth rates, since electron particle fluxes can balance the ions.Hence strong instability results for the ITG/TEM even though the corresponding ITG and TEM are stable.

However the greatly increased region where the FC is soluble allows a substantial instability.
In other words, the energetic cost of exciting the stabilizing trapped electrons to balance the FC is one that the ITG/TEM instability can and does pay.
This example shows the fallacy of considering the TEM as somehow independent from the ITG.One must understand the behavior of the coupled system; it is really a hybrid mode!

Using SKiM to understand the stabilization of ITG/TEM seen in simulations: the statistical mechanical ansatz vs the conventional picture
In the conventional picture, we expect that sufficiently large negative curvature (< ω d orbit e >) will stabilize the ITG/TEM.This is, indeed, what is found in SKiM (see figure 20(a)).For a certain small range of negative orbit average curvature, growth rate dependence upon F P is similar to the simulations: stabilization occurs near F P ∼1/2.But the simulations, the stability point in F P varies strongly as the curvature is varied.Even a 20% change in the negative curvature causes either complete stability for all F P , or, an instability up to large F P ∼0.6 − 0.7.
This strong variation of stabilizing F P with negative curvature is the typical behavior found by SKiM when there is a large trapped particle fraction.For example, as a further test, we reduce the ion curvature to zero (i.e. a slab-like ion mode) in figure 20(b).Again, stabilization at F P ∼1/2 is only possible for a small range of negative < ω dorbit >, and the mode becomes unstable for small changes beyond this.Contrasting mode behavior in the conventional picture with behavior from the statistical mechanical ansatz, using SKiM.(a) In the conventional picture < ftrap >= 0 remains large and the electron curvature is made more negative.Stabilization results for any value of F P depending upon how negative the curvature is.This is quite unlike the simulation results for progressively more negative curvature (b) This qualitative behavior is generic energetic stabilization as in the conventional picture: a similar scan with the ion average curvature set to zero behaves similarly (c) The growth rate in SKiM for a parameter scan like that in the statistical mechanical ansatz, where < ftrap >= 0 is reduced by eigenfunction adaptation to the progressively more negative curvature.For either positive or negative electron curvature, stability is reached a particular F P near the FC solubility limit for ITGae.(d) Growth rates from SKiM compared to simulation, where all eigenfunction averages are taken from the simulation eigenfunction and the SKiM γ is from its dispersion relation equation (18).The agreement is fairly close.Hence SKiM is a valid model to understand the simulations.This is quite the typical behavior in simulations (see figure 17) which find, that over an enormous range of curvature (figure 18), F P ∼1/2 always produces stability.
To understand simulations, let us appeal to adaptivity, a major element of what we have called the statistical mechanical ansatz: the eigenfunction adapts to avoid curvature stabilization by changing its structure.Simulations find that adaptation largely decouples from the trapped electron orbits with strongly negative curvature.By decoupling from most trapped electrons, the < f trap > becomes small.
Going back to SKiM, let us track γ as < f trap > is lowered (figure 20(c)) for two values of < ω d orbit e >= 0.05, −0.05,As we would expect, for low < f trap > the stability becomes similar to the ITG ae .The average < ω d orbit e > is taken from simulation eigenfunctions; it is often slightly positive because the eigenfunction decouples from the trapped electrons with the most stabilizing curvature.Notice that the stability point stays near F P ∼1/2 for low < f trap >.This qualitative behavior is the same if the sign of the electron curvature is changed, i.e. it is slightly negative.This happens because the non-adiabatic response is vanishing for either sign of curvature; the FC stabilization, energetic stabilization, is not strongly affected by curvature drive.Note the similarity of these results to the ITG ae : curvature affects the magnitude of the growth rate (here by a factor of ∼2) but the stability point in F P is not strongly affected.
Finally, we can evaluate the coefficients in SKiM for the actual eigenfunctions calculated in the simulations.When we do, SKiM gives good qualitative agreement with the simulation growth rates.These simulation eigenfunctions all give small < f trap > as we will see in the next section, and frequently < ω d orbit e > is small and positive (destabilizing).This, qualitatively, is just like the ITG ae in section 8, where the mode structure adapted to keep the average curvature drive positive; stabilization was due to the FC.
Let us now turn to E-FC plots.When parameters appearing in SKiM are averaged over the simulation eigenfunction, the stabilization with F P follows the ITG ae path.The ITB case C is shown in figure 21.Like the ITG ae , stabilization results because the FC curve shrinks to have low γ as F P is increased.And also like the ITG ae , the free energy curve alone, however, goes to large γ.On the other hand, plots of the conventional picture with curvature stabilization and stabilizing electron curvature look exactly the opposite (figure 20(c)).As we artificially increase the negative electron curvature, the free energy curve shrinks to be very close to γ∼0, but the FC can go to high γ.
These E-FC plots are qualitatively like those for the ITG ae , in the respective cases of FC stabilization and energetic (curvature) stabilization.
This difference in underlying physics would manifest experimentally: if the FC embodied physics controls the underlying physics, density gradients play a dominant role, whereas if energetics controls stability, one will need extremely strong curvature, and density gradients are not required.
In summary, the results with electron dynamics are qualitatively the same as for the ITG ae in section 8. From the statistical mechanical point of view, the addition of non-adiabatic electrons does not introduce any fundamental new element.
It is just that the solubility of the FC now depends upon the electron dynamics, so stabilization requires the non-adiabatic electron response to be small.This, in turn, requires that the eigenfunction must adapt to give such a small response, despite the fact that the equilibrium has a large fraction of trapped particles.This adaptation results from the eigenfunction avoiding the regions of stabilizing curvature, as found for the ITG ae .
The same fundamental dynamic operates for both the ITG ae and the ITG/TEM: the eigenfunction adapts to avoid energetic stabilization in the case of stabilizing curvature.But having done so, it is only the insolubility of the FC that can impose stability.

Pure TEM in the framework of the statistical mechanical ansatz
Studying the behavior of pure TEM (no ion temperature gradients) is theoretically interesting because it is an even more striking (than ITG ae and the ITG/TEM) manifestation of what we have called the statistical mechanical ansatz; it predicts that there are regimes where the TEM is stabilized by density gradients, and simulations find the predicted regime.
Ideally, any new theoretical understanding should predict some new and unexpected regime.Within the conventional picture, TEMs are driven unstable by the free energy in density gradients.So to find a regime where the opposite happens is indeed striking.
Since the TEM has considerably less free energy drive than the ITG/TEM, it can be stabilized by curvature more easily.But just before that happens, there is a regime where: • The eigenfunction adapts to decouple from the trapped electrons that sample highly stabilizing curvature • Hence the non-adiabatic response is small • Increasing density gradients, then, will tend to result in a much higher ion charge flux than the balancing electron charge flux • Consequently, the FC will become insoluble: the TEM will be stabilized by increasing density gradients, irrespective of energetic considerations The results of simulations, performed with geometries in figure 18 (where the TEM is not yet stable but is becoming weak), are shown in figure 23.
When density gradients are increased, the TEM at first becomes unstable.But then, with increasing density gradients, it becomes weaker or even stable.This behavior has been found for all the ITB geometries considered above.This situation is, surely, incomprehensible in previous pictures of the TEM, where only the free energy drive is included in the physical interpretation.This behavior (frames (a) and (b) in figure 23) occurs in all the TB geometries we have tested (when close enough to stabilization) and persists when plasma parameters are changed: impurities are included or/and when T i /T e > 1 figure 23(c In the second case, the growth rates increase with increasing density gradient.But for sufficiently negative curvature, the mode is stable for all density gradients.This is the characteristic behavior for stabilization by energetic considerations that we also see in the ITG ae , which is supposedly a very different type of mode.But when examined under the lens of the dual equations of the FC and free energy balance, the behavior of the ITG ae and TEM is quite similar.
To nail down the origin of stabilization, we examine the E-FC contour plots.First consider the case with low < f trap >= 0.2.Near the peak growth rate at F P = 0.28 the energy curves and FC curves are similar and intersect near their peak growth rate.However, at F P = 0.6 the situation is much different.The energy curve has gotten considerably larger both in maximum height and its horizontal extent.As expected for a TEM, increasing the density gradient has increased the free energy.The FC curve has gotten much lower, however, and is limited to low growth rates.The energy and FC curves have pulled apart so much, in fact, that there is no intersection at a positive growth rate.In short, the FC is stabilizing this mode.
This detailed study of pure TEM is very edifying-it shows how the mode behavior that is in striking violation of the conventional stability notions becomes perfectly understandable in terms of the FC based larger statistical mechanical perspective.We must remember that stability considerations must respect the simultaneous dictates of free energy and the FC; the FC could give stability despite very strong free energy drive (from whatever source) just as for the ITG ae .
When one species is nearly adiabatic (here the electrons), the fundamental tendency for thermodynamic forces to drive from SKiM showing that such qualitative behavior does not result from stabilizing electron curvature and large < ftrap >.Hence, only the statistical mechanical ansatz produces this behavior.Next we elucidate this regime using E-FC plots (f ) for a case in (d) at modest density gradient F P = 0.28, showing instability (g) At large density gradient F P = 0.6, the free energy is indeed larger than at F P = 0.28: free energy balance extends to higher γ, and over a larger extent in ωr However, the FC shrinks, as predicted by the statistical mechanical ansatz.The shrinkage is sufficient that no intersection between the FC and free energy balance occurs, so there is no instability.thermodynamic fluxes causes the FC to become insoluble as the density gradient is increased.This occurs in the simulations because the adaptive eigenfunction decoupled from the highly stabilizing regions of phases space, just like the behavior of the ITG ae The stability framework developed in this paper is general, and is by no means limited to a given set of modes.The statistical mechanical ansatz should apply to all, and with the same consequences for TB equilibrium as the ones in figure 18.

Summing up the lessons of pure TEM analysis
The results in the TEM section, very strongly, drive home the point: the FC is true for all fluctuations independent of what they are, what drives them, and how strongly.And if the plasma dynamics ensures that one species is nearly adiabatic (contributes little charge flux), there will be regimes where the FC insolubility will dictate stability.The resulting stability persists no matter what equilibrium gradients contribute to the free energy.
And, from a practical point of view, the FC-stabilized regimes become accessible, for instance, to curvature structures that are not nearly as extreme as those that are required for energetic stabilization.Creating TBs becomes easier in that sense.

Impurities and nonlinear heat fluxes for ITG/TEM
In earlier sections, we showed that the addition of impurities significantly affected the FC solubility for the ITG ae , and altered the F P threshold for stability.In the equilibrium scans above, we found that the eigenfunction adapts to reduce the non-adiabatic electron response.Hence, we expect that the addition of impurities to cases with electron dynamics would alter the F P value where D mix is strongly reduced, like the ITG ae .This is indeed the case.
In the scan of ŝ for ITB-like parameters (specifically: figure 17(a), let us include carbon impurities with Z eff = 2.We find that as ŝ becomes more negative, D mix is, now, strongly reduced at the commensurate F P of the ITG ae FC solubility limit (figure 24(a)).This F P is distinctly different from the expected bound when Z eff = 1.In figure 24(b), we consider a pedestal equilibrium with large enough α that the F P behavior was close to the ITG ae (specifically: pedestal case D, with α = 6 in figure 17(d).With impurities, the D mix is strongly reduced at F P similar to the solubility limit of the ITG ae .Now let us consider an equilibrium where the mode behavior is not quite as stable as in the ITG ae limit, i.e. when the non-adiabatic electron response is not small and plays a significant role.Specifically, consider the ITB case A with ŝ = −0.3 in figure 17(a), where the drop in D mix is not as sharp as for adiabatic electrons.What do impurities do in this equilibrium?impurities shift the curves to the right (figure 24(c)) qualitatively like the adiabatic electron case.And it makes the D mix substantially smaller than when there are no impurities, also qualitatively like the ITG ae .But the non-adiabatic electrons still allow instabilities up to significantly larger F P than for the corresponding ITG ae .
Nonetheless, the observed stabilization is very likely still due to the FC.As a test, we increase R/L T considerably, from 15 in figure 24(c) to R/L T = 20, 30 and 50 in figure 24(d).This increases the D mix by factors of 2-4, but the F P value for a precipitous drop in D mix changes only by ∼10 − 15%.This behavior is characteristic of the FC stabilization; the drop in D mix is still due to the FC.It is just that the non-adiabatic electron contribution to FC, canceling a part of the ion plus impurity current, drives the system towards FC solubility (towards destabilization).
Next, we show that the nonlinear simulations for the general ITG/TEM yield results similar to linear D mix just as for the ITG ae .We display two equilibria where the D mix showed a strong drop at F P close to the bound for the ITG ae : figures 25(a) and (b), represent, respectively, an ITB and a pedestal configuration.As we see, the total electrostatic nonlinear heat flux (the sum of electrons plus ions) behaves similarly to the D mix .
In these figures, the simulations were electromagnetic, but we have plotted only the electrostatic component of the heat flux.In some of these cases, electromagnetic modes (probably MTM) were present which produced a small heat flux.Such modes have not been part of our theoretical consideration, and obviously, any F P bound for such modes is not expected to be related to the ITG ae bound.Such modes, and their role in TB physics, will be considered in the future.Now, let us consider nonlinear results for configurations where the F P stabilization was not quite as strong as for ITG ae .For the pedestal case A (figure 17(d) with α = 6), the nonlinear behavior is again qualitatively similar to D mix , but the drop in the nonlinear flux is even faster (figure 25(c)).Finally, the ITB case A, also did not show stabilization as strongly as the ITG ae .With impurities included, the drop in the total electrostatic nonlinear heat flux follows D mix , but the drop is more rapid (especially in the case with strong impurity gradients).
A common feature, revealed by simulations, is that the nonlinear flux drops faster than the estimated D mix .Hence, the strong reductions in the D mix imply even stronger reductions in the nonlinear heat flux (most relevant to experimental confinement).

Summary-conclusions
We will now conclude our somewhat long enquiry into the physics of highly non-thermodynamic states like the TBs that are, by definition, endowed with strong gradients.States of high confinement are possible only when the turbulent transport originating by plasma instabilities can be minimized.Energetically, a plasma with high gradients should be violently unstable due to the enormous amount of free energy that can feed instabilities.This is the typical behavior of a vast variety of physical systems.And yet, TBs (both internal and edge) are often created in almost all tokamaks, stellarators, and other magnetic geometries.Simulation results including impurities in cases with full electron dynamics and electromagnetic effects (a) For Z eff = 2, as ŝ is made more negative, the stabilization occurs at the predicted FC bound for the ITGae including impurities (the thick dashed lines) rather than without impurities (the light dashed line).(b) A pedestal case, showing stabilization at the analytically predicted FC bound for the respective impurity parameters (c) ITB case A, where non-adiabatic electrons are small but large enough to extend the FC solubility limit.Nonetheless, qualitative trends with impurities have similarities to the ITGae (d) This same case varying the temperature gradient.The trends with F P are similar, verifying that the stabilization in the simulations is from FC insolubility and not from energetic effects.This paper is a detailed and comprehensive study, using analytical as well as computational tools, to figure out how states with strong gradients (therefore large free energy) manage to avoid feeding fluctuations that, in turn, could destroy them (by wiping out the gradients).
Before now, this situation has been understood via the properties of specific instabilities, the ITG and TEM, each with their unique stability properties.And when these unique and complex instabilities are coupled together, fortuitously, complex simulations find that they happen to have stability domains arising from their unique properties; such domains are, in fact, accessible to experiments.The accomplishment of this work is that this situation can be viewed very differently, and in the process, understood more clearly.Both these instabilities (and likely many others) are subject to basic dynamical relations in the gyrokinetic system: the restraining influence of the FC, and the destabilizing influence of the free energy balance equation.Since both the ITG and TEM are fluctuations of a dynamical system with very many degrees of freedom; the latter is expected to behave adaptively to tap free energy wherever possible.When viewed through the lens of the dual dynamical equations, the ITG ae , the ITG/TEM, and the TEM behave essentially the same: their stabilization can be understood via what we have called the statistical mechanical ansatz.When the FC becomes insoluble for growing modes, stabilization results.And because stabilization usually derives ultimately from the FC, basic trends follow from the essential physical character of the latter.Since thermodynamic forces drive their corresponding thermodynamic fluxes (down the gradients), when one species is close to adiabatic, the other will necessarily drive enough flux at high enough density gradients so that the FC has no solution.Let us review our deliberations in more detail.The most important new concept that we have developed and exploited in this paper originates in a characteristic property of fluctuations in a gyrokinetic system: the system will sustain only those classes of fluctuations for which the chargeweighted (radial) particle flux vanishes.The FC, though known in the literature, was never fully harnessed; we believe that it is the first time that its enormous stabilizing potential has been recognized and marshaled to develop what could qualify to be a 'theory' of fluctuations in a TB.
The essential theme in the constrained (rather than conventional) dynamics is that the Fluctuation behavior (stability) must be determined as a consequence of these two very different dynamics operating simultaneously: the FC and the conventional free-energy.We have concentrated on investigating the stability of the canonical ITG/TEM system.We have studied it extensively: in its various limiting cases, using different geometries, shaping, and plasma parameters (all with strong gradients).We were able to distinguish the workings of these two different dynamics, their separate roles, and their very different characteristics.In particular, we were able to identify what aspects of our extensive simulation results could be attributed to what dynamics.
In the process, we encountered an extremely interesting aspect of linear gyrokinetic instabilities-their adaptability; the corresponding eigenmodes change their structure so as to remain maximally unstable when subjected to stabilizing conditions (like stabilizing curvature) that will reduce free energy.For example, the fluctuations minimize interacting with stabilizing curvature, to a remarkable degree, even when such regions are both highly predominant and very deep.In principle, that would prevent the realization of TB-like states unless some other mechanism were to drive the system back to stability.In short, the gyrokinetic system manifests apparent adaptive behavior, and 'finds a way' to dissipate free energy until a 'hard' dynamical constrain prevents it-here, the FC.
Understanding adaptability is, theoretically, rather interesting.Understanding adaptability in a gyrokinetic system, however, is an absolute must for building a comprehensive theory of TBs.It is the adaptability that renders the standard free energy based stability analysis even more incapable of explaining the why of experimentally observed TBs (especially those without velocity shear).And it is the adaptability that allows the system to remain unstable until very close to the bounds set by the FC.Dramatic structural changes of the fluctuation can be understood as a response to the FC, after the latter is understood.
The structural changes in the eigenfunctions become comprehensible only in a statistically averaged (reduced) description.Contrary to previous understandings of the gyrokinetic system that are mode specific, the current approach transcends mode type and explains stabilization trends found in simulations for all the investigated modes ( ITG ae , ITG/TEM, or TEM).
Let us get back to the FC, the principal actor in this paper.Since the constraint pertains only to particle flux, the density gradients get, automatically, valorized in the theory.The simplest manifestation of the FC violation is seen in the suppression/elimination of fluctuations when density gradients exceed a certain threshold.Of course, the FC violation takes place when the electron flux cannot balance the combined flux of ions and impurities.By investigating the ITG/TEM dynamics in a large variety of situations (different geometries, different plasma parameters) we were able to demonstrate the universal working of the constraint; there is a large enough F P (the normalized density gradient fraction) beyond which the mode is either fully stable or stable enough that the turbulent transport is drastically lowered.
Although tokamaks have been the main focus of this paper, the same physics applies to stellarators.Notably, recent results on W7X have demonstrated the key role of density gradients in enabling high ion thermal confinement [37,38].In stellarators there are additional tactics (beyond those employed in tokamaks) available to realize the TB geometry.We expect that the framework outlined in this paper may find fruitful applications to optimizing stellarators for turbulent transport.
We will not repeat here all the 'results' discussed in detail in the text.We will, however, list the essence of our findings: • It is the constrained free energy dynamics-the conventional free energy considerations subjected to the zero FC-that determines the gross stability of gyrokinetic fluctuations.• In most cases relevant to TBs (high gradients), the FC is the crucial and efficient stabilizer-it succeeds in reducing turbulent transport, for instance, in much less demanding geometries than would be required by the free-energy route.• The insolubility of the FC at some level of F P seems to be possible in substantial regions of plasma-parameter space suggesting an equally large parameter space in which TBs are possible • The importance of the FC in suppressing fluctuations is even more impressive in the context of their adaptability, which lets them counter other stabilizing influences and remain maximally unstable.Thus reducing free energy, say by extreme geometries, may not stabilize the system as efficiently as desired.
• The FC is violated whenever the negative Flux ( = chargeweighted flux) of electrons cannot balance the positive Flux of ions plus impurities.So if the electron response to fluctuations is more adiabatic than that of ions, it is relatively easier for the FC to be insoluble when there are density gradients.(For ETG, as usual, the roles of electrons and ions should be reversed.)• Since the Flux is charge-weighted, impurities by virtue of their higher charge can play an outsize role in creating conditions hospitable for TBs.They can be the deciding factor even when the electron and ion Fluxes were to balance; this provides us a powerful knob to induce TBs.The action of impurities, thus, can be much more potent than the mere "dilution" of the main species.• The strongest knob for inducing stabilization by the FC is, however, the density gradient.
The reader can find in the text a very large survey of the parameter space for which detailed simulations were performed.
The aim was to chart out the existence of and the access route to territories where TBs (without velocity shear) could reside.
In order to be able to scan parameters easily to elucidate the essential physics, we conducted these simulations in simplified Miller-like geometries and controlled scans were highly educational as they were necessary for building confidence in the FC physics.The large number of variations that confirm the theoretical framework can be taken as showing that it has a substantial domain of validity.The next obvious step is to carry out this combination of analysis-simulations for reconstructed equilibria of experimentally observed TBs.It is likely, then, that we would be well equipped to find/perfect recipes to engineer these high confinement configurations in current and future machines.
it does not change the result that the response matrix Q ij must be positive definite, so the diagonal elements are positive.So the gyrokinetic response matrix, defined for an arbitrary fluctuation, has the same properties as the classic Onsager matrix, and for essentially similar reasons.

Appendix B. A more advanced version of SKiM
Here we explore why so simple a model as SKiM is so remarkably accurate, and why it can predict what simulations find in a far more complicated geometry.We will also derive a somewhat more complete FC relation that has the same solubility bounds.
The FC bound in SKiM follows from the fact that the perturbed distribution function can be written as a resonant denominator multiplying the driving term in the gyrokinetic system that arises from equilibrium gradients.The latter is the part due to [ω ⋆ ns + ω ⋆ Ts (mv 2 /2T − 3/2)]∼[1/L n − 3/(2L T ) + (v/v th ) 2 /L T )].The resonant denominator gives a weighting, in velocity space, of the plasma response to the gradient driving term.
For the flux to be zero, the resonance must weight both the positive and negative regions of the driving term equally.This becomes very challenging as F P → 0.6, since then, the constant in the driving term 1/L n − 3/(2L T ) become close to zero and slightly negative.Hence, over most of velocity space, the driving term [1/L n − 3/(2L T ) + (v/v th ) 2 /L T )] is positive.To solve the FC, the resonant plasma response must have a high concentration it the small region v/v th ∼0.
And to attain such high levels of concentration as F P → 0.6, the resonance must be progressively sharper, so the growth rate γ → 0.
The requisite concentration can happen by some combination of parallel resonances ∼v ∥ k ∥ and perpendicular resonances ∼ω d .
In the case of a more realistic geometry than SKiM, the plasma response would be different in detail, but as long as the plasma response could still concentrate near v/v th ∼0 the maximum value of F P would be the same.Hence, SKiM would correctly predict the solubility bound, despite having a different resonance structure.
Let us further consider the case of the FC bound for low k y perturbations where ω d is small compared to k ∥ v ∥ .In the resonant denominator, the plasma response is concentrated in a region where v ∥ ∼0.At low k ⊥ ρ i , this is no tendency to concentrate in the perpendicular direction.Hence the FC bond is considerably lower than in the case above: F P = 1/3.However, if k ⊥ ρ i becomes larger, the Bessel function J 0 reduces the response at high k ⊥ ρ i .This is another way of saying that the response becomes relatively more concentrated at low v ⊥ .This allows F P to increase beyond F P = 1/3.But the Bessel function does not concentrate as strongly as a resonance might.So the maximum F P is less than 0.6, and the SKiM calculations above find it is about 0.53.
Therefore, for a more accurate model than SKiM, as long as the plasma response is due to parallel ion motion could concentrate near v ∥ ∼0, and as long as Bessel function effects reduced the response in the perpendicular direction, the maximum value of F P would be similar to SKiM.Naturally, the plasma response could be quite different in detail.
It can be greatly simplified if we take k ⊥ , ω d and v ∥ to be constant.Introducing the Fourier transform in l, ϕ(l) = ´ϕk ∥ e ik ∥ l , we obtain for the FC: This is just the SKiM FC, equation (10), for a given k ∥ integrated over the values in the spectrum, weighted by |ϕ k ∥ | 2 .In SKiM, there is only a single average k ∥ , which is a spectrum average (the SKiM expression is equivalent to . Equation (B3) is more realistic, in that it has a full spectrum of k ∥ contributions, weighted by the Fourier amplitudes |ϕ k ∥ | 2 .Crucially, the FC solubility bounds for equation (B3) are identically the same as SKiM.When the SKiM FC is insoluble, it means the flux is positive definite for all k ∥ .Then the integral of such positive definite contributions weighted by the positive quantity |ϕ k ∥ | 2 is also positive definite, and the converse.This also applies to the FC bounds for the low k y case as well.
For any choice of ϕ(l), a reasonable estimate for the average k ⊥ and ω d are the SKiM values, that is an average in l weighted by |ϕ(l)| 2 .
Even without the approximations leading to equation (B3), the basic properties of the plasma response will give a similar result.For example, if the response is concentrated to low v ∥ , and k ⊥ ρ i is small, the FC bound will be ∼1/3, and it will increase due to the J 0 as k ⊥ ρ i increases.In other words, even the exact response will have similar properties to SKiM insofar as the solubility of the FC.The actual response will include numerous complex effects such as trapped ion orbits, etc.This could have a major effect on a dispersion relation based upon equation (B1).But it is not so surprising that the maximum bound on solubility will be similar to that from SKiM.

Appendix C. SKiM with trapped electrons
Now we define the parameters < ω d orbit e > and < f trap > used throughout this paper.These can be derived in various simplifying limits, but they are sufficiently intuitive that we simply present them here.In SKiM, the parameters for trapped particles are: where <> orbit is the orbit average or 'bounce average', i.e. the time average over a bounce orbit (for mathematical simplicity we take <> orbit to vanish outside the trapped region, to make equations (C2) and (C1) more transparent to write).The dl ′ dΩ v is a weighting over phase space, where dl ′ is the configuration space weight, and the relevant velocity coordinate for trapped particles is the velocity pitch angle, with a weight proportional to the solid angle, dΩ v .The average trapped electron curvature < ω d orbit e > is evidently the orbit averaged curvature < ω d,e > orbit weighted by the magnitude of the orbit averaged response of that orbit < ϕ > orbit .It is the obvious generalization of the ion expression < ω di >= ´dl ′ |ϕ| 2 ω d / ´dl ′ |ϕ| 2 to a trapped species, where all quantities are orbit averaged.
The effective trapped particle fraction f trap averages the magnitude of the response ∼ < ϕ > 2 orbit over phase space and compares it to ϕ 2 .If the bounce average response were as large as possible < ϕ > 2 orbit = ϕ 2 , this expression would be the actual physical fraction of trapped particles.But because the orbit averaged response can be considerably less than this, the 'effective' trapped fraction f trap for a particular eigenfunction can be considerably less than the actual fraction.
Heuristically, the magnitude of the non-adiabatic electron response from trapped particles is ∼f trap , whereas < ω d orbit e > gives whether the curvature acting on that response is energetically stabilizing or destabilizing (destabilizing is positive).
Values obtained from equation (C1) and (C2) for several cases we considered are shown in figures 26 and 27.

Figure 1 .
Figure 1.The growth rate of the most unstable mode (maximum for all k θ ρ i ) is plotted (a) vs F P for a standard geometry like the cyclone base case) (b) for four diverse geometries vs R/L T showing that the critical temperature gradient varies by ∼270% (c) for those same geometries, the critical F P varies by ∼5% (d) for two stellarator geometries (e) for an ITB geometry, for many values of T i /Te (f ) for a low aspect ratio A geometry, an ETG for two different λ Debye .

Figure 2 .
Figure 2. (a) Heat flux vs F P for ten diverse nonlinear simulations (b) Plot of the soluble and insoluble regions of the FC for low k modes: the F P solubility limit depends upon k ⊥ ρ i .

Figure 3 .
Figure 3. (a) The solubility regions together with the simulation results (the open circles)for the < k ⊥ ρ i > from the eigenfunction.The unstable eigenfunctions track the analytically computed boundary to stay in the soluble region, by increasing < k ⊥ ρ i > as the boundary is approached, (b) The eigenfunction for various F P values (denoted by arrows in (a)).To increase the < k ⊥ ρ i >, the eigenfunctions progressively broaden to larger θ.

Figure 4 .
Figure 4. < k ⊥ ρ i > from simulation eigenfunction vs F P for representative cases (the open circles).The unstable eigenfunctions track the analytically computed boundary to stay in the soluble region, by increasing < k ⊥ ρ i > as the boundary is approached.

Figure 5 .
Figure 5.The simulation eigenfunctions [Reψ (Imψ)normalized to be unity (zero)and at the origin] for the four cases in figure4for F P = 0 and F P near the stability limit.The value of < k ⊥ ρ i > (θ) for the geometry is also shown (the dashed line).Near the solubility limit, the eigenfunction structure shifts to reside in regions of larger < k ⊥ ρ i >.The eigenfunction structure is often quite complex, and differs strongly from case to case.Nonetheless, the eigenfunction averaged < k ⊥ ρ i > (θ) tracks the analytic expression.

Figure 6 .
Figure 6.Simulation eigenfunction (Reψ , Imψ ) for representative cases, for F P near the stability limit.We include cases where kx is not zero on the outer midplane, i.e. η 0 is not zero.The value of < k ⊥ ρ i > (θ) for the geometry is also shown (the dashed line).Near the solubility limit, the eigenfunction structure shifts to reside in regions of larger < k ⊥ ρ i >.The eigenfunction structure is often quite complex, and differs strongly from case to case.Nonetheless, the eigenfunction averaged < k ⊥ ρ i > (θ) tracks the analytic expression.

Figure 7 .
Figure 7. D mix (the maximum for all all three η 0 ) and the growth rate (also the maximum for all three η 0 ) is compared to the nonlinear heat flux for many diverse cases.For clarity, results are normalized to unity at F P = 0.The D mix follows the nonlinear result considerably more closely than γ, indicating that the increase in < k ⊥ ρ i > plays an important role in the nonlinear flux reduction.In (h) and (i), nonlinear results for three values of the temperature ratio are shown.Despite the fact that the transport level is strongly affected, D mix reaches low values for nearly the same value of F P .

Figure 8 .
Figure8.The minimum value of < k ⊥ ρ i > is plotted for equilibrium from table 1, along with the analytically computed solubility boundary.The magenta crosses are the < k ⊥ ρ i > for η 0 = 0, and the solid diamonds are the minimum < k ⊥ ρ i > among the three η 0 simulated (for unstable modes).Hollow circles are also the minimum among the three η 0 , but these are the cases that are only feebly unstable (D mix < 0.02 of the maximum of the scan).As expected, the addition of more degrees of freedom (multiple η 0 ) allows some modes to approach the solubility boundary more closely.Sometimes, feeble modes can exist slightly beyond the boundary, so the analytic theory is not perfect.But the overall agreement is very good.Despite feeble modes crossing the boundary, the overall agreement is very good, especially considering that approximations must be used to obtain the analytic results.

Figure 9 .
Figure 9.The complex ω (ωr-γ) plane with various plots scanning a range of Fp from 0 to 0.44.The blue contour denotes the line satisfying the energy equation and the red contour denotes the line satisfying zero charge flux.Eigenvalues lie at the intersection of the two.Dark orange regions denote areas with positive charge flux and light regions denote areas of negative charge flux.Note the progressively lower line satisfying the zero charge flux constraint, forcing growth rates to small values and indicating access to the transport barrier regime.

Figure 10 .
Figure 10.The complex ω (ωr-γ) plane with various plots scanning a range of eigenmode-averaged drift frequencies.The blue contour denotes the line satisfying the energy equation and the red contour denotes the line satisfying zero charge flux.Eigenvalues lie at the intersection of the two.Dark orange regions denote areas with positive charge flux and light regions denote areas of negative charge flux.Note the progressively lower line satisfying the zero charge flux constraint, forcing growth rates to small values and indicating access to the transport barrier regime.

Figure 11 .
Figure 11.(a) The curvature ω d i(θ) for the equilibrium sequences shown, for decreasing ŝ.It becomes overwhelmingly stabilizing at strongly negative ŝ (b) D mix from simulations for k θ ρ i = 0.35; the maximum D mix for three values of η 0 are shown (c) the eigenfunction averaged curvature for F P = 0 for both the true eigenfunction and the 'unadapted' eigenfunction for case 1; the simulation eigenfunction adapts to maintain destabilizing curvature, and if it had not, average curvature would become very stabilizing (d) Growth rates from SKiM for parameters representative of these eigenfunctions, but for various < ω d >.Without adaptation, the curvature would stabilize the mode independent of the F P bound.

Figure 12 .
Figure12.Plots of the eigenfunction averaged < k ⊥ ρ i > and average curvature < ω d > for cases 1-6.The < k ⊥ ρ i > follows the analytic bound for all cases.The average curvature remains positive except when the eigenfunction is called upon to simultaneously adapt to the energy and the FC while the curvature is extremely preponderantly negative.

Figure 13 .
Figure13.FC solubility boundaries including carbon impurity in a deuterium plasma.The case with Z eff = 2 has the same density gradient scale length for all species, whereas the case Z eff = 2 1/L nZ = 3/L nd has three times steeper impurity gradient than deuterium.The F P is defined in terms of the electron density gradient, F P = 1/Lne/(1/Lne + 1/L Ti ) which is the most readily measured experimental quantity.

Figure 14 .
Figure 14.The simulated < k ⊥ > ρ i vs F P for several equilibrium, including impurities, plotted with the analytic FC solubility bound.The impurity parameters and analytic predictions are identical to those in fig(13) (recall this is for carbon impurity).Green is Z eff = 1, purple is Z eff = 2 and red is Z eff = 2 but the impurity density gradient scale length is 1/3 of the deuterium gradient.

Figure 15 .
Figure 15.The simulated D mix and nonlinear heat flux for diverse cases with impurities.Figures ((a) (c) (e)) [((b) (d) (f ))] plot these quantities for ITB [pedestal-like] configurations.Both D mix and nonlinear heat fluxes become appropriately small for all cases as the analytic bound in Fp is approached.Notice that the temperature ratio, although it strongly affects the magnitude of the heat flux, has very little effect on the value of Fp at which the flux becomes very small-totally as predicted by the theory.

Figure 17 .
Figure 17.D mix from simulations with full electron dynamics and electromagnet effects, in scans of ŝ and α.The corresponding curvatures are shown in figure18.In all cases, increases in ŝ or α eventually reduce the threshold F P needed for very low D mix down to the limit of the ITGae, but even large increases beyond that do not materially reduce that F P .

Figure 20 .
Figure 20.Contrasting mode behavior in the conventional picture with behavior from the statistical mechanical ansatz, using SKiM.(a) In the conventional picture < ftrap >= 0 remains large and the electron curvature is made more negative.Stabilization results for any value of F P depending upon how negative the curvature is.This is quite unlike the simulation results for progressively more negative curvature (b) This qualitative behavior is generic energetic stabilization as in the conventional picture: a similar scan with the ion average curvature set to zero behaves similarly (c) The growth rate in SKiM for a parameter scan like that in the statistical mechanical ansatz, where < ftrap >= 0 is reduced by eigenfunction adaptation to the progressively more negative curvature.For either positive or negative electron curvature, stability is reached a particular F P near the FC solubility limit for ITGae.(d) Growth rates from SKiM compared to simulation, where all eigenfunction averages are taken from the simulation eigenfunction and the SKiM γ is from its dispersion relation equation(18).The agreement is fairly close.Hence SKiM is a valid model to understand the simulations.

Figure 21 .
Figure 21.E-FC plot for the representative TB scan (ITB case C in (figure 20(d)), with eigenfunction averaged parameters computed from the simulation eigenfunction.As F P is increased in (a)-(d), the FC curve shrinks while the free energy curve does not.The stabilization seen in the simulations is clearly from the flux constraint.

Figure 22 .
Figure 22.E-FC plot for a case like the conventional picture of curvature stabilization (figure 20), for constant F P = 0.5.The curvature becomes more negative in (a)-(d).The free energy curve shrinks till the mode becomes stable.Stability is primarily from the free energy balance forcing γ to be near zero.
). Guided by SKiM, two possible scenarios of TEM stabilization could emerge • reduce < f trap > maintaining slightly positive < ω d orbit e > • the conventional picture without adaptation: < f trap >= 0.45 and < ω d orbit e > is scanned to progressively more stabilizing values Results are shown in figure 23.For the first case, SKiM and simulations are alike: the TEM is initially destabilized by increasing density gradients and is, then, stabilized, reflecting the characteristic behavior of FC stabilization.

Figure 23 .
Figure 23.The regime where density gradients ultimately stabilize the TEM, an unexpected regime predicted by the statistical mechanical ansatz.(a) Simulations of ITB case A with varying ŝ Increasing density gradients first destabilize, and then stabilize, the TEM, as predicted (b) Similar behavior from ITB case D with varying α (c) Similar behavior arises with impurities and T i /Te not equal to one (d) Results from SKiM for the TEM and representative eigenfunction averaged parameters, showing that at small < ftrap > this behavior manifests (e) Resultsfrom SKiM showing that such qualitative behavior does not result from stabilizing electron curvature and large < ftrap >.Hence, only the statistical mechanical ansatz produces this behavior.Next we elucidate this regime using E-FC plots (f ) for a case in (d) at modest density gradient F P = 0.28, showing instability (g) At large density gradient F P = 0.6, the free energy is indeed larger than at F P = 0.28: free energy balance extends to higher γ, and over a larger extent in ωr However, the FC shrinks, as predicted by the statistical mechanical ansatz.The shrinkage is sufficient that no intersection between the FC and free energy balance occurs, so there is no instability.

Figure 24 .
Figure 24.Simulation results including impurities in cases with full electron dynamics and electromagnetic effects (a) For Z eff = 2, as ŝ is made more negative, the stabilization occurs at the predicted FC bound for the ITGae including impurities (the thick dashed lines) rather than without impurities (the light dashed line).(b) A pedestal case, showing stabilization at the analytically predicted FC bound for the respective impurity parameters (c) ITB case A, where non-adiabatic electrons are small but large enough to extend the FC solubility limit.Nonetheless, qualitative trends with impurities have similarities to the ITGae (d) This same case varying the temperature gradient.The trends with F P are similar, verifying that the stabilization in the simulations is from FC insolubility and not from energetic effects.

Figure 25 .
Figure 25.Nonlinear simulation results including full electron dynamics and electromagnetic effects compared to D mix .The D mix is scaled to the nonlinear result at low F P to see if the trends are similar as F P is increased.First, cases without impurities.When non-adiabatic electrons are weak in (a) and (b) the trends of the nonlinear heat flux match the trends in D mix rather closely.(c) A pedestal case where non-adiabatic electron effects are non-negligible, the trends are qualitatively quite similar but not as close quantitiatively.(d) An ITB case with non-negligible non-adiabatic electron effects, and where impurities are included.The trends in nonlinear heat flux match the trends in D mix fairly well.

Figure 26 .
Figure 26.Evolution of eigenfunction parameters as F P increases for magnetic shear scans in figures 17 and 18.

Figure 27 .
Figure 27.Evolution of eigenfunction parameters as F P increases for α scans in figures 17 and 18.

Table 1 .
Parameters used for the simulations.

Table 2 .
Parameters for simulations.