Guilhem Dif-Pradalier

Peter Beyer

Peter Manz

Raul Sanchez

Carlos Hidalgo

Rémi Dumont

Pascale Hennequin

Aspects of Turbulence Self-Organisation in Magnetised Fusion Plasmas

Keywords: Bohm scaling?, 25

2.3.3 Practical implications for fusion: perspectives for model reduction . 29 2.3.3.1 What does it suggest for transport scaling? . . . . . .

Chapter 1

Introduction "Les jeunes gens qui débutent dans la [physique], les amateurs qui se tournent vers la [physique] seront-ils longtemps encore satisfaits de travailler dans la nuit, sans pouvoir répondre à aucune interrogation sur le sens et la portée de la recherche où ils s'engagent ? [. . . ] On rencontre cependant tous ces gens, tous ces jeunes gens qui croient que tous les travaux [formels] amènent un profit à l'espèce humaine, parce qu'on leur a persuadé qu'il en va ainsi de toutes les tâches spirituelles. [. . . ] Sous prétexte que je lisais tard des livres en comprenant plus facilement qu'un ajusteur n'eût fait [les lois de la thermodynamique], je ne me prenais pas pour un homme anonyme, je croyais docilement que l'ouvrier dans la rue, le paysan dans sa ferme me devaient de la reconnaissance puisque je me consacrais d'une manière noble, pure et désintéressée à la spécialité du spirituel au profit de l'homme en général, qui comprend, parmi ses espèces, des ouvriers et des fermiers. Mes maîtres faisaient tout pour m'entretenir au sein d'une illusion si agréable pour euxmêmes. [. . . ] C'est ainsi que bien des naïfs pensent passivement que [toutes les epistémées] méritent la gratitude des hommes. Cette croyance annonce le mythe de la cléricature."

Paul Nizan

In most of everyday life, the matter which surrounds us is organised in convoluted structures, assemblies of more elementary constituents. As ambient thermal energy increases, becomes comparable or surpasses the binding energies of such assemblies the coherence of these familiar structures weakens or disappears. The elementary constituents, most of them intrinsically charged become more susceptible to the influence of other elementary constituents and in turn directly influence them as well. They also become more susceptible to the influence of external electromagnetic fields. Collective motions of great complexity arise, intrinsic to such media. Thermonuclear plasmas display such characteristics -primarily characterised by the unrelenting motion of its interacting constituents (ions, electrons and neutrals) and the excitation of a broad variety of collective dynamical modes. The general definition of a plasma is quite loose and difficult since this state of matter1 encompasses a continuous spectrum of media which increasingly deviate from the concept of gas (see Table 1.1), finally ending up formalising a distinct state of matter, of which a thermonuclear plasma can be seen as its extreme limit.

The dynamics of such a medium spans a broad range of spatially and temporally correlated different scales: from the fast electron density plasma oscillations at ω p = v th /λ D [with v th = T /m and λ D = ε 0 T /e 2 n] or the fast Larmor gyromotion ω c along a magnetic field line, to the intermediate (in the adiabatic limit) bounce frequencies ω b of Table 1.1: Intrinsic complexity of the thermonuclear plasma state as compared to gas.

Property

Gas Fusion plasmas

Electrical Conductivity

Very low Usually very high (higher than copper at thermonuclear temperatures ∼ T 3/2 ) Independently acting species (w/o chemical reactions)

One; all gas particles behave in a similar way, influenced by gravity and binary collisions

Two or three: electrons, ions, and neutrals; distinguished by their charge, they behave semi-independently, coupled through electromagnetic fields and long-range collisions; this allows for new instabilities

Velocity distribution

Maxwellian (far from the boundaries) since collisions are dominant; few hot tails of fast particles Often non-Maxwellian; collisional interactions are weak, external forcing can drive the plasma far from local equilibrium, with a significant population of fast particles Interactions Binary; three-body collisions are extremely rare Binary: within a sphere of Debye; collective outside with long-range interactions due to electric and magnetic mean fields.

particles trapped in local inhomogeneities of the magnetic field, collisional frequencies ν x,y [with x, y ∈(electrons, ions)], and down to the inverse of the energy confinement time τ -1 E which characterises the rate at which a system loses energy to its environment. Regarding ω p ω c,e ω c,i ω b,e ω turb ω b,i ν ii ν ei τ -1 E ν ee 5 10 11 10 9 10 7 10 6 10 5 10 5 10 1 1 0.2 (s -1 ) spatial scales, they range from the Fermi scale λ SF where the strong force dominates or the tunnelling length λ tunn. which allows fusion reactions to occur without beating the full Coulomb potential barrier to the electron ρ e and ion ρ i Larmor scales at which most electrostatic instabilities develop up to the system scale a and the large mean free path λ mf p characteristic of weakly collisional media.

λ SF λ tunn. λ D ρ e δ b,e ρ i δ b,i a R λ mf p 10 -15 10 -12 10 -4 10 -4 5 10 -4 10 -3 5 10 -3 1 3 10 3 (m)

Interactions occur across scales. These interactions determine the overall organisation of the system, making the problem conceptually intricate and computationally intensive. The present manuscript is concerned with the exploration of some of these organisations and with their practical incidence in the quest for controlled fusion on Earth.

∴

• "Le chaos est impossible aux yeux de la raison, car il est impossible que, l'intelligence étant éternelle, il y ait jamais eu quelque chose d'opposé aux lois de l'intelligence ; car le chaos est précisément l'opposé de toutes les lois de la nature"

Voltaire

• "Chaos should be regarded as extremely good news"

Chogyam Trungpa Rinpoche

Nuclear fusion is the elementary process, accompanied by the release or absorption of energy, by which multiple atomic particles combine to form heavier nuclei. If the energy to initiate the reaction comes from accelerating one of these particles, the process is called beam-target fusion; if both nuclei are accelerated, it is beam-beam fusion. If the nuclei are part of a plasma sustained out of thermal equilibrium, one speaks of thermonuclear fusion, which represents today's most promising approach for a controlled source of energy production on Earth. We will be focusing here on the latter.

Given the criticality of zero-carbon energy production, efficiency of fusion power production and successful operation of the forthcoming Iter experiment is hugely important. The interpretation of experimental data, optimisation of performance and safe forecasting of operation in fusion devices depend on the reliability of a hierarchy of models, from fast reduced models used in the control room to 5 or 6-dimensional first principles models, important to test fundamental ideas and calibrate the reduced models. Large gradients -large reservoirs of free energy-are inherent to the magnetic fusion configuration and sources for magneto-hydro-dynamic (MHD) instabilities and for microturbulence. The former mostly constrains accessible operational space whilst the latter controls the performance of the confined plasma. This dissertation is primarily interested in questions related to turbulence, from a perspective of basic understanding: how it self-organises, how it transports, what the organising principles are. It is largely based on 5-dimensional so-called flux-driven gyrokinetic modelling.

The gyrokinetic approach is a reduction from the fully 6-dimensional kinetic model, making use of the strong guiding magnetic field inherent to the tokamak configuration. The prevailing paradigm for gyrokinetic computations assumes separations of scales between rapidly evolving micro-scale fluctuations, typically at electron ρ e or ion ρ i Larmor scales atop temporally fixed and smooth macro-scale background gradients, typically at the scale a or R of the device. Modelling strategies which enforce this scale separation are "gradient-driven": fluctuation-fluctuation feedback on the mean is discarded. Mean gradients are fixed in time and small-amplitude perturbations are computed about this state [Sydora 1995, Kotschenreuther et al. 1995, Beer & Hammett 1996, Dimits et al. 1996, Candy & Waltz 2003, McMillan et al. 2008, Peeters et al. 2009, Gorler et al. 2011]. This scale separation allows for computational efficiency. Such frameworks may either be local (flux-tube, radial periodicity is enforced) or global (no radial periodicity enforced). They have proven fruitful in the plasma core.

One difficulty comes from the fact that mean profiles are commonly "stiff", which often occurs in situations where external drive acts on profiles on timescales which are long compared to relaxation rates through instabilities. In such situations, modest variations of mean gradients may lead to large excursions in expected fluxes. Dynamics of mean gradients may also be important to adequately assess nonlinear feedback loops (e.g. through shear) and thus quantify transport. Flux-gradient relations are indeed not always found to be single-valued (the same flux may be obtained for different values of the gradients, see e.g. Ref. [Diamond et al. 1997]), rather leading to the concept of fluxgradient landscape. Furthermore, detailed knowledge of mean gradients is marred, from an experimental standpoint by the combination of limited temporal frequency acquisition and spatial uncertainties -especially in the outer regions of the plasma. For all of the above, constraining mean gradients and asking for fluxes is an ad-hoc constraint on the accessible physics that needs to be carefully assessed, especially near marginal (nonlinear) stability. Near marginality is to be understood in a nonlinear sense, i.e. as a regime of (i) minimal available free energy or driving source which allows (ii) sustainment of a finite level of turbulence (i.e. is nonlinear in nature). At marginality, mean profiles may greatly differ from profiles that would lead to linear instability onset. Marginal (nonlinear) profile gradients may indeed be either more gradual than their linearly unstable counterparts (as e.g. in Self Organised Criticality) or steeper (as e.g. in the case of a shear-stabilised "Dimits upshift"). Near marginal regimes are a central topic of this work. They are believed to be relevant for fusion because instability growth is usually fast with respect to dissipation processes, leading to so-called 'stiffness': the clamping of profiles near turbulence onset which minimises available free energy. These marginal regimes are believed to become even more important to fusion in future devices such as Iter where the ratio of external heating over confined volume of plasma decreases as compared to current machines. As in all devices the external driving/heating power that is installed is fixed; Iter being much larger than current-day machines, the ability there to 'over-drive' or push past marginal stability will likely be more difficult than in some of current-day devices, possibly further restricting accessible operational space to the vicinity of near-marginal regimes.

A more physically consistent approach -which in particular consistently describes these near marginal regimes-has been proposed in which imposed fluxes, which are known experimentally drive the system and both fluctuations and mean gradients dynamically evolve on the same timescale, driving and back-reacting on one another. This "flux-driven" framework [Carreras et al. 1996, Sarazin & Ghendrih 1998, Garbet & Waltz 1998, Beyer et al. 2000, Chang & Ku 2008, Ku et al. 2009, Idomura et al. 2009, Sarazin et al. 2010, Korpilo et al. 2016, Grandgirard et al. 2016, Wang et al. 2018, Villard et al. 2019] relaxes assumptions of scale separation, keeps the fluctuation-fluctuation feedback on the mean (mean gradients become an output of the computation) and results in a tenfold or more increase in computational demands2 . Physically, in analogy to local models which have started investigation through refined resolution of the involved problem of (multiscale) electron and ion-scale coupling [Howard et al. 2014, Maeyama et al. 2015], flux-driven frameworks intrinsically add to gradient-driven approaches the study of (multiscale) ion-scale and equilibrium coupling.

These questions are at the heart of the current manuscript.

• In Chapter 2, we investigate turbulence self-organisation near marginal (nonlinear) stability. With questions from the above paragraph in mind, we have strived to broaden our approach to related fields of knowledge, using analogies and work from other communities to inform our understanding of near-marginal plasma behaviour.

• In Chapter 3, we show multiscale interplay between core, edge and Scrape-Off Layer, spontaneous establishment of a stable transport barrier in the outer edge (at the closed/open field line interface) and highlight the roles of turbulence spreading and boundary conditions in the understanding of edge dynamics. The plasma edge is where the more important transport bifurcations develop; steps towards a better understanding of its dynamics thus have important implications for fusion.

∴

"All good things arrive unto them that wait and don't die in the meantime"

Mark Twain

As last remarks, whilst writing the present manuscript we had to make choices. In order to provide a consistent overview of our ongoing body of research, the unifying thread of the self-organisation of turbulence has been chosen. This necessarily sets aside some of our other works more briefly touched upon, in particular some regarding fluid closures [Sarazin et al. 2009, Gillot et al. 2021a], advanced equilibria for the Vlasov system [Ogawa et al. 2016a, Ogawa et al. 2016b, Laribi et al. 2019] or the concerned with the understanding of Edge Localised Mode (ELM) dynamics and their control by Resonant Magnetic Perturbations (RMPs) in the H-mode (edge) pedestal, using reduced magneto-hydro-dynamics (MHD) [Orain et al. 2013, Bécoulet et al. 2014, Orain et al. 2014, Orain et al. 2015, Fil et al. 2015, Morales et al. 2016 Gillot et al. 2022] Abstract: Tokamaks, like many dynamical systems prone to developing instabilities display threshold behaviour. The regime of (nonlinear) near marginal stability is particularly interesting both conceptually and for practical purposes. Large, hot devices such as Iter are expected to be stiff, i.e. display a propensity to remain in the vicinity of said marginality, making this regime an important one to understand. It is characterised by the competition between secondary structures -zonal mean flows and fronts/avalanches-which organise into tertiary structures -staircases. This self-organisation, a robust outcome of proximity to marginal stability is shared by a large variety of different systems, from sub-cellular to astrophysical scales. In fusion, it was first observed in numerical computations, later leading to its experimental characterisation. Self-organisation of these competing mesoscopic structures and their impact on transport raise acute questions. Amongst possible benefits, it may provide a route to understand the transition, as the size of the plasma column increases from unfavourable 'Bohm'-like to more favourable 'gyro-Bohm' scalings for heat confinement. The question however of its detailed incidence on macroscopic transport levels remains largely uncharted. We review below some key items of near-marginal organisation: what are the expectations from a computational perspective; how it compares to (difficult) experimental measurements and what frameworks of understanding have been proposed from theory. In addition, we show possible shortcomings of leading modelling approaches in the community when organisation near marginality is not adequately considered. In particular: large and systematic underpredictions of fluxes, flows and transport with implications for extrapolation to Iter and to future devices. Interestingly however, constructive recommendations as to how reduced models must be improved near marginal stability, especially through refoundation of commonly used turbulence saturation rules may be formulated.

Forewords

Anti-diffusion ('negative viscosity') phenomena are the remarkable processes by which upgradient transport leads to mesoscopic structure formation in many physical and biological systems, including those that are chaotic or turbulent. The formation of layers or staircases -i.e. regular patterns in which steep-gradient layers (transport barriers) are situated between zones of stronger mixing -is a typical outcome of anti-diffusion. This state of inhomogeneous mixing resembles closely that of spontaneous phase separation, familiar from statistical mechanics, and so is contrary to the naively expected limiting case of homogenisation. Layered systems exhibit very different transport properties to those of familiar weakly inhomogeneous systems. Thus, understanding the physics of anti-diffusion and layering is crucial to developing an accurate description of transport in multi-phase and active fluids, metallurgy, magnetically-confined plasmas, oceanographic and atmospheric flows -which determine climate -and stellar interiors. Important fundamental questions remain unanswered concerning an underlying common mathematical foundation for the different physical phenomena. Furthermore, mathematically, the anti-diffusion problem in its starkest form is ill-posed. There is thus considerable interest in how the problem can be regularised, and in the dynamics resulting from different approaches. Computationally, the efficient and accurate representation of states with two or more highly disparate scales represents a significant challenge. Much therefore can be gained from a broad, interdisciplinary endeavour in order to foster inter-disciplinary exchanges between researchers in geophysical and astrophysical fluid dynamics, plasma physicists, researchers in active and biological matter, and soft-condensed-matter physicists.

In this spirit, the author has participated, partially organised or organised several "Festivals de Théorie" which precisely aim at fostering such interactions, participated in the 2014 KITP programme "Wave-Flow Interaction in Geophysics, Climate, Astrophysics, and Plasmas", the 2015 BIRS programme "The Mathematics of Layers and Interfaces", partially organised the 2021 KITP programme "Layering in Atmospheres, Oceans and Plasmas" and is one of the organisers of the 2024 Isaac Newton Institute "Anti-Diffusive Dynamics from Sub-Cellular to Astrophysical Scales".

2.2 Lessons learnt (or to be learnt) from other fields, ranging from sub-cellular to astrophysical scales "It would seem reasonable, therefore, to class as unnatural any turbulence experiment with artificial forcing that interferes too strongly with PV evolution. Phenomena like PV mixing, jet selfsharpening, and PV staircase formation are unlikely to be well represented in experiments that strongly violate [the material invariance of PV]"

M. E. [START_REF] Dritschel | Multiple Jets as PV Staircases: The Phillips Effect and the Resilience of Eddy-Transport Barriers[END_REF] Many systems display layering behaviour. Elucidating basic principles through interdisciplinary cross-fertilisation could play an important role in understanding, characterising the impact and predicting behaviour in magnetised plasmas.

Brief forewords on formal analogies between fusion plasmas & planetary fluids

This section provides a brief overview of the more detailed discussion from Appendix B.1.

There are clear analogies between aspects of astrophysical fluid dynamics (AFD) or of geophysical fluid dynamics (GFD) and of magnetised fusion plasmas (MFE). The staircase idea itself originated in GFD and AFD. Some of the simplest non-trivial models for each system share deep formal structures. These models are the celebrated systematic scaling theory for large-scale atmospheric motions by Charney which leads to the derivation of the quasi-geostrophic equations of GFD [Vallis 2006] or the well-known Hasegawa-Mima (HM) Table 2.1: The planetary-tokamak analogies summarised.

model [START_REF] Hasegawa & Mima | Pseudo-three-dimensional turbulence in magnetized nonuniform plasma[END_REF] of MFE. Each of these models obey the following equation:

∂ ∂t -(∇f × z) • ∇ f -∇ 2 f + log (F) Potential Vorticity (PV) = 0 (2.1)
where f = ϕ stands for the electric potential in the MFE case and for the depth fluctuation f = h in the GFD context. Similarly, the free energy F = n0 ωc,i is either the ion background density gradient (MFE) or the mean layer depth F = H0 f (GFD). All notations are those of Appendix B.1 and the more detailed analogy is summarised in Table 2.2.1. Analogies are such that Eq.(2.1) is often referred to as the Charney-Hawegawa-Mima (CHM) model.

In the presence of inhomogeneities, the CHM model Eq.(2.1) predicts the existence of a linear wave that will admit in both GFD and MFE the same structure for its dispersion relation and thus similar propagating properties. Waves in fluids result from the action of restoring forces on fluid parcels displaced from their equilibrium positions. In the GFD case the dominant restoring force comes from rotation and the associated wave is the well-known westward Rossby wave; in plasmas, rotation or gravity are negligible and the dominant force comes from compressibility (it may also come from electromagnetics -this effect and associated Alfvèn waves are ignored here). The so-called drift waves are akin to sound waves and are the MFE analogue of the Rossby waves. Often, the more general Hasegawa-Wakatani (HW) framework in MFE or the generalised Charney equations in GFD are preferred due to the possibility that they offer to describe either drift wave or interchange instabilities in MFE (i.e. turbulent transport) or their atmosphere-ocean counterpart, the barotropic instability. However, the CHM model captures essential physics of drift-Rossby wave Zonal Flow formation while involving the minimal number of relevant scalar fields [Diamond et al. 2008].

Equation (2.1) in both contexts can be compactly written as d(PV)/dt which expresses the material invariance of a quantity (PV) known in both cases as the 'potential vorticity'. This conservation of (PV) on each fluid element carries important consequences: PV being comprised of two parts, a fluctuation and a mean, it straightforwardly implies a tight connection between how dynamics that unfold at the fluctuation scales f will be leaving a footprint at the larger, equilibrium mean scales F. Staircases are precisely that: emerging (tertiary in MFE) structures borne of small-scale fluctuations and developing (storing free energy) at larger meso-and macro-scales via flows and profile corrugations.

Layering in Magnetised Plasmas

Scale selection and secondary pattern formation are critical physics issues underpinning the quest to achieve good confinement in magnetised plasmas. On the one hand, avalanches -which occur when an array of convective cells respond and overturn sequentiallyconstitute bursty, intermittent transport events, which can span the system size (L). On the other, secondary zonal flows -driven by turbulent Reynolds stresses produced by the cell pattern -tend to shear apart the large avalanches. Since the basic cells are a few ion gyroradii (ρ i ) in size, and ρ i ≪ L, this pattern competition plays out over a large dynamic range. Magnetised plasmas evolve according to a reduced form of the Vlasov/Boltzmann equation, and thus have many features in common with systems governed by potential vorticity (PV) conservation -see below, section 2.2.5 and appendix B. In particular, in colisionless plasmas, the phase space distribution is the effective PV, conserved along particle trajectories.

The E × B staircase is the remarkable outcome of this pattern competition. Such staircases are quasi-periodic arrays of transport barriers -namely domains of local gradient steepening, interspersed by zones of strong mixing. Staircases are observed in simulations and in experiments. A fundamental issue is the physics underpinning the inhomogeneous mixing and gradient sharpening, which defines the staircase. In particular, the theory should determine the causal relation between the profile steepening zones (barriers) and the spatially contiguous secondary E × B flow shear layers. One promising candidate explanation is bistable (two-scale) mixing, in which one mixing scale is emergent and related to the profile gradient scale itself (cf. the Rhines scale in geophysical fluids), thus triggering transport bifurcations. Another is heat flux jamming (which resembles traffic jams), which can arise from avalanching. The time delay between the heat flux response and changes to the local temperature gradient is critical to jamming. More generally, the resilience of staircase patterns in a background of avalanches, and the various types of long-time staircase evolution -including coarsening and migration -are as yet not well understood.

While models are developing, they have yet to confront the fully multiscale character of the problem. Furthermore, they have also yet to inform us on the quantitative impact that such structures may have on the overall transport. This is an area that we contributed to develop in MFE, which is under active investigation in MFE and beyond. There is much to gain from detailed confrontation to other systems which display layering behaviour (as made abundantly clear below) since the development of basic understanding across these various fields may bear some important and practical fruit for fusion as well -eg. (i) the discrepancies in heat flux prediction between state-of-the-art models (gyrokinetic or quasilinear), should staircase patterning be accounted for or not, the open questions of (ii) the exact role of avalanching in heat and particle transport, or (iii) the pathway to macroscopic transport barriers from a merging process of weaker, separated barriers (or said differently, assess the viability of E × B staircase states as alternative regimes of enhanced confinement, in comparison with more familiar transport barrier scenarios).

Multiphase Flows & Active Fluids

Multiphase flows and, in particular, phase separation below the critical point, have been active areas of research for more than half a century [Toner & Tu 1998, Marchetti et al. 2013, Pandit et al. 2017]. Over the past decade, spurred by advances in experimental techniques and the growth of computational facilities, there has been considerable interest in examining the suppression of such phase separation by the turbulence in a well-stirred binaryor ternary-fluid mixture. For example, turbulence-suppressed phase separation has been studied by coupling the Cahn-Hilliard equations (CH), used to study phase separation in two-phase systems, with the Navier-Stokes (NS) equations. The interfaces between the coexisting phases play an important role in regulating multiphase turbulence. Indeed, the CHNS system is arguably the simplest paradigm for the study of the competition between phase separation and mixing, a theme that is fundamental to the study of anti-diffusion.

Much theoretical work on active fluids has focused on bulk and one-component systems. However, the theoretical understanding of boundaries and interfaces, whether generated spontaneously upon transitions between active phases, or forming through segregation in multi-species bacterial suspensions, is essential for exploiting active matter and active flows with the goal of engineering new materials and devices, as well as for elucidating the role of active processes in biology. For instance, the notion of liquid-liquid phase separation of immiscible active fluids has been brought to the forefront by the discovery of membrane-free organelles and subnuclear structures. Therefore, an important open challenge is understanding and quantifying the interfacial properties of active fluids.

Active fluids exhibit spontaneous spatio-temporally chaotic flows that have been dubbed 'active turbulence'. This behaviour can occur at essentially zero Reynolds number and is distinct from the well-known inertial turbulence of high Reynolds number fluids. It is perhaps more akin to the elastic turbulence observed twenty years ago in sheared polymer solutions, where nonlinear elastic effects can destabilise laminar flow. It remains an interesting open question as to whether ideas from dynamical systems (e.g. topological entropy) will be useful for quantifying such turbulence.

Finally, another important direction stems from recent experimental demonstrations that one can use light to generate spatial patterns of active/passive regions in active suspensions of cytoskeletal filaments cross-linked by optogenetically modified motor proteins. It is then interesting to explore the use of spatially and temporally modulated activity to control active flows and the dynamics of topological defects in orientationally ordered active liquid crystals. By understanding how the spatial and temporal modulation scales of activity interact with the length and time scales of active energy injection, one can hope to achieve quantitative control of the patterns of active flow. This of course echoes the question, central in magnetised plasmas of the understanding and tailoring of triggered transport transitions through shear layering and/or shear mergers.

Layering & Mixing in Stratified Fluids

There is abundant observational evidence -oceanographic and atmospheric -for the formation of layers in single-component stably-stratified flows, where stratification results from the gradient of temperature or other scalar [Balmforth et al. 1998, Caulfield 2021, Radko 2013]. The layers may result from a variety of forcings (e.g. the Phillips mechanism by which the stirring of a stratified flow leads to layering) or through the nonlinear development of instabilities (such as the Kelvin-Helmholtz instability with either vertical or horizontal shear). Understanding the propensity of density-stratified fluids to form layers, together with the associated transport of heat and other scalars, such as pollutants, is a fundamental and critical problem in geophysical and environmental fluid dynamics.

In fluids where the density depends crucially on two components with disparate diffusivities, there is the further possibility of double-diffusive convection. The most widely studied example -one of the mechanisms driving oceanic flows -is thermohaline convection, resulting from competing heat and salt gradients and in which the diffusion of salt is much slower than that of heat. Thermohaline convection can occur in two very different regimes: in the fingering regime the temperature gradient is stabilising and the solutal gradient destabilising; conversely, in the diffusive regime, the roles of both gradients are reversed. Remarkably, in both regimes, long-lived, robust density staircases have been observed: in the fingering regime in warm seas, such as the Mediterranean, and in the diffusive regime in colder regions, particularly in the Arctic. Although various mechanisms have been proposed for double-diffusive layer formation -metastable states, horizontal intrusions, and a variety of instability mechanisms -precisely how layering occurs remains a topic of great research interest.

Double diffusive layering is also of potential importance astrophysically, particular in deep stellar interiors. Here the two competing ingredients may be gradients in temperature and composition or, alternatively, temperature and magnetic field (the magnetic pressure producing the necessary buoyancy). Unlike thermohaline convection in the oceans, however, the astrophysical parameters are extreme, impossible to simulate exactly. In order to be able to incorporate the large-scale and longtime influence of layering on transport into stellar evolution models, it is therefore vital, as for layering in single-component systems, that the physics of layering and how it depends on the various parameters is understood. Double-diffusive phenomena are also likely to play a role within the confined plasma, at the confluence where turbulent "diffusion" from both core and edge and impurity diffusion from both edge and core meet. How such processes may lead to flow onset is very much an open question.

Geophysical (GFD) & Astrophysical (AFD) Fluid Dynamics

Potential vorticity (or PV) is a generalised vorticity that is conserved along fluid trajectories in rotating GFD flow [START_REF] Dritschel | Multiple Jets as PV Staircases: The Phillips Effect and the Resilience of Eddy-Transport Barriers[END_REF], Vallis 2006]. Thus GFD, often called the 'Fluid Dynamics of Potential Vorticity', is governed by the incompressible advection of a conserved phase space density, in turn related to the advecting velocity via the Biot-Savart law. Geophysical flows are often quasi-two-dimensional, being strongly constrained by rotation; such flows tend to display anti-diffusive behaviour -as found in inverse cascades -with up-gradient transport. Inhomogeneous PV mixing leads to the formation of PV staircases, in which PV layers, manifested as zonal jets, are separated by small regions of sharp PV gradients, which act as barriers to PV mixing. Thus the system effectively decomposes into domains of strong PV mixing (i.e. jets), separated by transport barriers, in which the mixing is weak.

The structure of PV staircases and the dynamics of PV mixing thus underpin the zonal jet pattern, as observed in the atmospheres of the gas giant planets. Critical questions include the micro-physics of inhomogeneous mixing and the origin of barriers, and the role of linear and nonlinear stability criteria in constraining the PV staircase structure. Astrophysically, the precise role played by the PV staircase in determining the structure of the Jovian zonal jets is a topic of great interest and some controversy. Formal analogies between GFD and magnetised plasmas have been well emphasised over the years. A partial overview of relevant ideas for flow formation in both GFD/AFD and MFE may be found in appendix B.

Climate & Atmosphere-Ocean Interactions

It is well appreciated that while climate models are able to capture the large-scale flows, referred to as the 'general circulation', they cannot resolve small-scale processes as occur in naturally occurring layers. Broadly speaking, layers in the atmosphere and ocean [START_REF] Timmermans | Understanding Arctic Ocean Circulation: A Review of Ocean Dynamics in a Changing Climate[END_REF] are manifest in one of two dynamically distinct forms: primarily vertically stacked layers separated by rapid density variations, and primarily horizontal layers separated by strong shear flows.

Density interfaces and fronts act as barriers to transport, which can be harmful or beneficial. Atmospheric inversions can trap pollution close to the ground, inhibiting its dispersion until morning-time convection erodes the interface. The stratospheric polar vortex above Antarctica inhibits dispersion of CFCs during the southern spring, resulting in the formation of the ozone hole. By contrast, the polar vortex over the Arctic is more rapidly broken up by atmospheric waves created by flow over the underlying mountain ranges; thus the ozone hole is less pronounced. In the ocean, our present climate relies upon the transport of heat and CO 2 absorbed at the surface being carried to depth by wintertime convection. A grave concern is that increased glacier melt, particularly in the Labrador Sea, will introduce so much fresh water on the surface as to create a pycnocline with such a strong density contrast that convection cannot penetrate through it, thereby trapping heat and greenhouse gases near the surface.

These matters may appear as somewhat distant from our prime interests in fusion plasmas. This is partially true; there are however important ideas that echo some of our own questions: e.g. how to treat layering in large scale climate models and what practical implications this may have. As an obvious analogue, we just have to remind the reader that in MFE, much of scenario forecasting and experimental interpretation is performed based on reduced models, especially using quasilinear approximations. Both in climate models and plasma simulators, a most important issue is to understand how turbulent transport may be influenced by layering, and then to be able to model this. This is crucial since direct numerical simulations in the realistic parameter regimes are either simply impossible or so time consuming that parameter space cannot be adequately parsed.

Cross-cutting themes to help address key questions for MFE

There are several key questions that one would seek to address due to the observation of layering and of secondary or tertiary structures in drift-wave interchange turbulence. Many are of direct practical interest for fusion modelling, some are broad enough to have implications across the various disciplines. In particular, below are some cross-cutting themes of significant portent for MFE:

• What are the various mathematical approaches to anti-diffusive behaviour? In what ways might the problem be regularised? How might the disparate length scales that characterise layered states be encompassed in an asymptotic approach? Relatively simple PDE models in time and one spatial dimension have been devised to explain how layering may occur in a stirred stratified fluid, and recently these have been adopted to explain staircase formation in drift wave plasma turbulence. Can we identify, discuss and elucidate many different proposed mechanisms of layering, such as bistable inhomogeneous mixing, 'jamming' or cross-diffusion?

• To make progress in understanding layering in plasma turbulence across a wide range of parameters, we need to address how do layers and interfaces form and evolve in various settings. What factors control the 'strength' of the interface between layers?

What factors (waves, convection, dynamical instability) result in the erosion of the interface? How strong must an interface be to resist being eroded? Do mergers between layers occur and are we able to trigger or control access to macroscopic transport barriers through such routes (and is it even desirable)?

• Layered states have a range of disparate length scales: from the system size, down to the interlayer scale and then down to the very small interfacial scales. Modelling staircase formation thus poses a significant computational challenge. To date, computational studies mostly employ uniform domains (in physical or spectral space), simplified boundary conditions and forcings; it is important therefore to both explore techniques such as adaptive mesh refinement and to examine the influence of different types of forcing and dissipation, together with the role of geometrical and boundary constraints.

• On even more practical grounds, it is key to answer how is turbulent transport influenced by layering and then, if relevant (as indications seem to show), how can anti-diffusion and flow patterning be incorporated efficiently in reduced models of plasma turbulence.

2.3

The plasma staircase: predicted computationally et al. 2017b] the factors that control this spatial segregation, the dynamics of this global organisation as well as its incidence on confinement. The regime of near-marginality [Diamond & Hahm 1995] is not merely academicallyappealing for its rich dynamics, secondary pattern formation and mesoscale organisation. It is as well a likely operating regime for current and future magnetic confinement devices. Non-burning Iter plasmas for instance, characterised by a large ratio of plasma volume over external heating are likely to reside close to marginal stability. Fusion plasmas indeed display the property, common in dynamical systems that upon surpassing a critical threshold, an instability may promptly build up, inducing large fluxes which deplete the driving gradients and inhibit the instability. Background gradients thus hover in the vicinity of nonlinear near marginal thresholds. Large, hot devices such as Iter are expected to be stiff due to the temperature dependence of the gyroBohm heat flux scaling, making near marginality a regime which models must confront. Many strategies have been devised in modelling to mimic natural processes. All are not equivalent and different choices may critically affect the nature of computed statistical equilibria. This point is further discussed in section 2.3.2.1.

Turbulence spreading [START_REF] Garbet | Radial propagation of turbulence in tokamaks[END_REF], Hahm et al. 2004, Yi et al. 2015, Heinonen & Diamond 2019, Singh & Diamond 2020, Dif-Pradalier et al. 2022] in particular is ubiquitous in flux-driven regimes, its more visible manifestations being relaxation events through emission of fronts or avalanches [Hwa & Kardar 1992, Carreras et al. 1996, Newman et al. 1996, Sarazin & Ghendrih 1998, Idomura et al. 2009, Sarazin et al. 2010, Dif-Pradalier et al. 2010, Van Compernolle et al. 2015, Wang et al. 2020].

This fact underlies departure from a local or Fickian transport paradigm [Dif-Pradalier et al. 2010, Ida et al. 2015, Hahm & Diamond 2018]. Connection between spreading and departure from local, single-valued flux-gradient landscapes is still debated. Turbulence regulation through zonal flow formation [Rosenbluth & Hinton 1998, Diamond et al. 2005] is yet another natural trend of near-marginal drift-wave turbulence. One of the merits of E × B staircase patterning is certainly to provide, whilst segregating regions where avalanching is dominant from regions where zonal flow concentrates (forming the staircase microbarriers) a natural and dynamic means for both antagonistic trends to simultaneously exist. Salient features of the E × B staircase are represented in Figs.2.1 and 2.2 and summarised below (when details are absent, further illustration of the assertions below may be found eg. in Ref. [Dif-Pradalier et al. 2017b]):

(i). Ubiquitous signatures of staircase formation are mean profile corrugations, reflecting the propensity of an initially homogeneous near-marginal profile in a bath of driftwave (or Rossby wave) turbulence to develop inhomogeneities that will tend to selfsharpen and endure. Self-sharpening relies upon positive feedback and may occur e.g. through negative viscosity-like arguments [Phillips 1972, Balmforth et al. 1998, Hinton 1991, Ashourvan & Diamond 2016], clustering instability [Flynn et al. 2009, Kosuga et al. 2013, Kosuga et al. 2014] or electric field curvature [Itoh & Itoh 2016]. Which mechanism is operative in what region of parameter space is still unclear and will be focused on elsewhere. The mean pressure profile corrugations lead to strong localised temperature gradients (the "risers" in Fig. 2.1). The name staircase comes from the step-like idealised resulting pressure profile with, statistically, quasiconstant mesoscale step spacing ∆. Corrugations are shown in Fig. 2.2 on the mean ion temperature gradient ∇T temporally-averaged over 0.53 ms between 1030 a/c s and 1343 a/c s . The step size ∆ is discussed below, item (v);

(ii). Co-located with the mean profile corrugations are zonal mean flows (section 2.3.2.1 and Fig. ) or vorticity jumps that define lasting (from a few collision times to [at least] fractions of the energy confinement time) "valleys" in configuration-time space of hindered transport. These zonal mean flows, of typical radial extent δ flow ∼ 10 ρ s act as a set of weak or permeable transport barriers [Ghendrih et al. 2014c, Dif-Pradalier et al. 2015, Hornung et al. 2017] that regulate the turbulent heat flux.

They are thus also referred to as "microbarriers" throughout the paper and are visible in Fig. In the plasma case, avalanches and corrugations interplay: the intense mixing region downstream of a corrugation is the locus of emission of the heat avalanches and an important local free energy reservoir for the nonlocal avalanche-mediated transport.

On the other hand, the staircase pattern of shear layers acts as regularly spaced (with spacing ∼ ∆) weak elastic transport barriers: some avalanches carry a sufficient momentum to carry across these shear layers -which are disturbed or even break before reforming in the avalanche wake as e.g. around time 1760 a/c s in Fig. 2.2but most do not, halt and dump their energy upstream of the corrugations. A related behaviour has been observed in 2D flux-driven modeling of the plasma edge [Ghendrih et al. 2014a].

(iii). The E × B staircase is a near-marginal pattern that progressively weakens and disappears with increasing distance to instability threshold;

(iv). The E × B staircase is a weakly collisional pattern, as expected from the fact that its microbarriers are zonal mean flows, hence sensitive to collisional damping [Hinton & Rosenbluth 1999];

(v). Avalanches, in low to moderately heat-driven plasmas are ubiquitous. They play an important role for staircase dynamics given the near-marginal character of the latter. Avalanching is the dominant transport mechanism in-between the staircase shear Figure 2.2: Three general features of the plasma staircase are visible here: (i) the mean profile corrugations here displayed on the temperature gradient, (ii) the strong, long-lived and coherent shear flows defining "valleys" of hindered transport-the mean radial E × B shear profile is shown in Fig. 2.6 (left) and (iii) the radial transport dominated by avalanchelike events in-between the staircase steps. layers (the "treads" in Fig. 2.1), see item (vi) below. Avalanches, importantly, are statistically contained by the staircase steps (section 2.3.2.2) and rarely propagate across successive staircase shear layers and over a significant fraction of the plasma volume.

The size distribution of avalanches therefore closely follows the size distribution of ∆ (the 'tread width' or 'step size'). The most probable step size ∆ stat ∼ 40 -50 ρ s is mesoscale (section 2.3.2.2), which despite a fat tail [Dif-Pradalier et al. 2017b, Hornung et al. 2017, Wang et al. 2020, Milovanov et al. 2021a] implies a favourable mesoscale (and not a macro-scale) most probable scale length for avalanches. Interestingly, despite strongly non-Gaussian statistics the staircase mesoscale ∆ stat reintroduces a favourable gyro-Bohm like scaling for heat confinement (section 2.3.2.2);

Q turb (r) = -K(r, r ′ ) ∇T (r ′ ) dr ′ (2.2)
where K is a generalised diffusivity best described in the general case as a Lévy distribution and may adequately be approximated as a Cauchy-Lorentz distribution:

K(r, r ′ ) = Λ π A /2 (A /2) 2 + |r -r ′ | 2 (2.3)
For all positions r ′ ∈ [0, a], a being the minor radius, K(r 0 , r ′ ) defines the extension over which an action of the thermodynamic force ∇T (r ′ ) at location r ′ may extend and impact the transport processes at r 0 . Let A be that distance, defined as the typical e-folding scale length of the kernel K (|r 0 -r ′ |). Finding K is thus key to discriminating the nature of the transport interactions in the plasma: ei- the autocorrelation length of the turbulence. We typically find an e-folding scale length A ∼ 10 ℓ c , which emphasises that the transport is not local, not diffusive and not described by Fick's law. In the regions of the 'risers', heat transport is dominantly neoclassical [Dif-Pradalier et al. 2015]: these regions are microbarriers for transport. ∆ stat ∼ 40 ρ s is a natural 'nonlocal' mesoscale for the system, in addition to the microscale autocorrelation of the turbulence ℓ c ∼ 5 -8 ρ s . These two micro-and mesoscales, shown in Fig. 2.4 echo the two-scale turbulence correlation length [START_REF] Hennequin | [END_REF] recently measured in Asdex Upgrade. Remarkably, nonlocal avalanching fills-in all the space left unaffected by the immediate vicinity of the shear pattern at the staircase steps such that A ∼ = ∆: the quasi-regular staircase steps ∆ define the outer scale of avalanche propagation (further details in section 2.3.2.2);

ther local if 0 ≤ A ≤ ℓ c or genuinely nonlocal if A ≫ ℓ c , with ℓ c representing
(vii). The flux-gradient relation is often multivalued and non-monotonic near marginal stability, in contrast to oft-reported results in Ion Temperature Gradient turbulence [Dimits et al. 2000]. In particular local values of the fluxes may not bijectively relate to local values of the gradients as different gradient drives may be associated to the same flux (multivalued) and larger gradient drives may be associated to lower fluxes (non-monotonic). Observation of coherent structures, as discussed here for the E × B staircase is often indication of a breakdown of the local flux-gradient paradigm;

(viii). The whole staircase pattern dynamically evolves and meanders, river-like over intermediate to long (collisional to confinement) timescales as the staircase steps show propensity to remain at constant ambient gradient drive;

(ix). Where staircase steps (risers) nucleate does not necessarily rely upon special values of low-order safety factor q rationals, yet synergistic reinforcement of shear layers, as they meander, in proximity of a low-order q rational is sometimes observed. The connection between low-order q rationals and staircase steps is tenuous with a simplified Boltmann response. A fully kinetic response leads to an interesting situation where the overall organisation is the convolution of 2 patterns: one is associated with localised marked electron temperature corrugations in the immediate vicinity of low order q rationals, mostly due to passing electron dynamics [Dominski et al. 2012].

These fine radial scale structures seem to have a modest overall impact on transport [Rath et al. 2021]. They lie on top of a second, larger mesoscale pattern which is the globally-organised staircase observed regardless of the electron response (Boltzmann or kinetic). This latter organisation on the other hand has a marked impact on transport in the confined core of tokamaks. Interestingly, similar patterns are found as well in stallarators, with a more pronounced connection between loworder q rationals and the permeable barriers. Such conclusions have been drawn when analysing experimental heat transport and poloidal rotation data from the TJ-II heliac [van Milligen et al. 2022]. The possible stronger connection between flows/staircase steps and low-order rational ι = 1/q surfaces in stellarators than in tokamaks poses interesting open questions such as the relative contribution of magnetohydrodynamics (MHD) to transport or different possible synergies between MHD and turbulence in both configurations. This is left for future work.

The results above are robust properties of the E × B staircase, in L-mode-like plasmas, based upon tens of gyrokinetic calculations representative of varying plasma parameters in drift-wave Ion Temperature Gradient (ITG) turbulence with Boltzmann electrons, using the Gysela code [Grandgirard et al. 2007, Grandgirard et al. 2016]. Different main ion species are considered: Hydrogen, Deuterium or Helium, different collisionalities , 6], different heating mechanisms: thermal baths, gradient-driven and flux-driven and different boundary conditions. These findings have been confirmed in various independent studies [Dif-Pradalier et al. 2010, Kosuga et al. 2013, Villard et al. 2013, Kosuga et al. 2014, Ghendrih et al. 2014c, Dif-Pradalier et al. 2015, Rath et al. 2016, Ashourvan & Diamond 2016, Peeters et al. 2016, Hornung et al. 2017, Wang et al. 2018, Guo et al. 2019, Choi et al. 2019, Wang et al. 2020, Rath et al. 2021, Qi et al. 2021, Milovanov et al. 2021a, Garbet et al. 2021, Rath & Peeters 2022]. Subscripts 0 stand for initial evaluation at mid-radius.

ν ⋆ = 0 and ν ⋆,0 ∈ [0.005, 10], different plasma sizes ρ -1 ⋆,0 ∈ [128, 512], different gradient drives R/L T ∈ [2, 12] or η = L n /L T ∈ [1
Most computations of E × B staircase formation have been run with the simplified assumption of a Boltzmann electron response. This impacts turbulence organisation in the vicinity of rational q surfaces (see discussion above, last item), even for ion-scale (ITG or Trapped Electron Mode) turbulence. Of course, realistic flux-driven computations including both kinetic ion and electron responses are required to assess the generic character of staircase organisation, especially as ion and electron channels may interplay [Hornung et al. 2017]. Studies interestingly have started relaxing the adiabatic or electrostatic assumptions [Rath et al. 2021, Rath & Peeters 2022] and tend to confirm that the above picture still endures as the base state on top of which finer-scale structures develop. More work is needed to really establish this fact yet the robust observation of such patterns in varieties of experimental plasma conditions [Dif-Pradalier et al. 2015, Hornung et al. 2017, Choi et al. 2019, Qi et al. 2021] allows to gain confidence in these assertions.

Salient staircase properties, relevant for fusion plasmas

Here, we focus on some of the more interesting properties of the plasma staircase. The interested reader may find additional features, e.g. in [Dif-Pradalier et al. 2017b] that will not be reported here.

Forcing & 'distance to marginality' matter

The prevailing paradigm for gyrokinetic computations assumes separations of scales between rapidly evolving micro-scale fluctuations, typically at electron ρ e or ion ρ i Larmor scales atop temporally fixed and smooth macro-scale background gradients, typically at the scale L of the device.

Modelling strategies which enforce this scale separation are "gradient-driven": fluctuation-fluctuation feedback on the mean is discarded.

Mean gradients are fixed in time and small-amplitude perturbations are computed about this state [Sydora 1995, Kotschenreuther et al. 1995, Beer & Hammett 1996, Dimits et al. 1996, Candy & Waltz 2003, McMillan et al. 2008, Peeters et al. 2009, Gorler et al. 2011]. The scale separation allows for computational efficiency. Such frameworks may either be local (flux-tube, periodicity is enforced) or global (no periodicity enforced). They have proven fecund in the plasma core.

Experimentally however, driving power is known whilst all other quantities, including mean gradients are not. A paradigm shift has been proposed where known, imposed fluxes drive the system and both mean gradients and fluctuations dynamically evolve, driving and back-reacting on one another. This "flux-driven" framework [Carreras et al. 1996, Sarazin & Ghendrih 1998, Garbet & Waltz 1998, Beyer et al. 2000, Chang & Ku 2008, Ku et al. 2009, Idomura et al. 2009, Sarazin et al. 2010, Korpilo et al. 2016, Grandgirard et al. 2016, Wang et al. 2018, Villard et al. 2019] relaxes assumptions of scale separation, keeps the fluctuation-fluctuation feedback on the mean (mean gradients become an output of the computation) and results in a tenfold or more increase in computational demands.

Validity of gradient-driven approaches is debated, especially at the plasma edge -where scale separations are questionable, or near marginal stability [Diamond & Hahm 1995, Sanchez & Newman 2015] -i.e. between linear and nonlinear thresholds, in the socalled "Dimits shift" region [Rosenbluth & Hinton 1998, Dimits et al. 2000].

Fluxdriven approaches indeed show near marginality that distinct turbulent states and flux levels can be obtained from almost identical underlying mean plasma gradients. This fact underlies a multivalued (bistable mostly), non monotonic flux-gradient landscape [Diamond et al. 1997[START_REF] Nakata | [END_REF], Dif-Pradalier et al. 2017a, Dif-Pradalier et al. 2017b, Villard et al. 2019], with competition [Miki et al. 2012] between mean (equilibrium) flow shear and zonal (fluctuating) flow shear. As underlying shear varies, gradients dynamically evolve at almost constant flux or conversely, flux levels vary may be reached at almost constant gradients.

Connection between spreading and departure from local, single-valued flux-gradient landscapes is still debated -see e.g. the review papers [Sanchez & Newman 2015, Ida et al. 2015, Hahm & Diamond 2018]. In order to test breakdown of Fickian transport at microscales or impact of scale separations near marginal stability, we have repeatedly run Gysela either in its native flux-driven configuration or 'downgraded' it to an equivalent 'gradient-driven' framework. Comparing FD and GD approaches within the same code is especially interesting for the flux-driven approach sets no constraints upon largescale avalanche-type front propagation and allows for consistent feedback between flows, shear and mean profiles and staircase organisation clearly depends on the latter, possibly on both. In practice, the flux-driven computation is first run up to flux equilibrium. Based on the end state, the computation is then restarted with gradient-driven forcing. Several conclusions can be drawn:

(i) the staircase organisation disappears with gradient-driven forcing. Space and time characteristics of the self-organised turbulent state is compared in Figs.2.5 and 2.6 for flux-and gradient-driven regimes with identical plasma parameters, in the banana collisional regime (ν ⋆ = 0.2 at ρ = 0.5). The staircase pattern visible in Fig. 2.6 that has built up and endures in the flux-driven computation is rapidly [in less than a millisecond] damped away and ultimately lost in the gradient-driven approach. The time interval and colourbars for which data is displayed in Fig. 2.6 are identical. The temperature corrugations and strong inhomogeneities in mean shear and poloidal flow, characteristic of the staircase pattern give way in the gradient-driven framework to a more homogeneous turbulence, with fluctuations more isotropic, flow and shear patterns more random in time and evenly distributed poloidally. This observation is confirmed in Refs. [Rath et al. 2016, Peeters et al. 2016], where gradient-driven models are reported to lead to spurious secular growth of zonal flow activity near marginal stability.

(ii) our understanding of near-marginal stiffness is clearly questioned and echoes the recent observation of an increased nonlinear upshift of the temperature gradient length threshold, attributed to E × B staircase onset [Rath et al. 2016]. These questions point towards the role of staircase organisation and generally of structure formation near marginal stability in the possible breakdown of the local flux-gradient paradigm.

As such they may be seen as a generalisation of earlier heuristic models where e.g. the heat diffusivity χ = -Q/∇T is itself a nonlinear function of either zonal flow shear [Hinton 1991] or of ∇T itself.

(iii) the fact that large though highly intermittent flux-driven fluxes are observed at low temperature gradient drives (location [0]) where both the gradient-driven flux and its statistical variations remain small tends to indicate clear sensitivity of the fluxgradient relation to temporal coarse-graining, despite the fact that the present results are flux-surface averaged and thus already strongly spatially coarse-grained. This result echoes results in Refs. [START_REF] Ghendrih | Thermodynamical and microscopic properties of turbulent transport in the edge plasma[END_REF], Idomura et al. 2009]. The local flux-gradient relation is violated (becomes multivalued) in a flux-driven framework at short to intermediate (collisional) time scales;

(iv) the gradient-driven framework appears to be overly stiff. Strong stiffness endangers predictive reliability in a model as large variations in predicted fluxes can be obtained from modest or even possibly indistinguishable differences in ambient mean profiles (within statistical fluctuations or experimental error bars). Reliability of gradientdriven predictions near marginal stability is thus questioned. The flux-driven framework is less stiff and appears less sensitive to details of the ambient mean profiles.

The above findings echo significant discrepancies in heat transport predictions between flux-driven approaches and gradient-driven or quasilinear approaches (section 2.3.3 and Ref. [Gillot et al. 2022]) has been reported near marginality. These observations may have important consequences yet are still in their infancy and require further investigation.

(v) Strong stiffness indeed impedes nonlocal turbulence propagation and especially turbulence spreading [Hahm et al. 2004], an important ingredient in the plasma edge [Dif-Pradalier et al. 2017a, Dif-Pradalier et al. 2022]. These observations have implications for the so-called 'transport shortfall' problem [Scott et al. 1991, Holland et al. 2011], as discussed in chapter 3.

Local flux-gradient relations may thus appear as a poor description of flux-driven nearmarginal transport since secondary (zonal flows, avalanches) and tertiary (staircases) structures are not adequately represented, making transport processes possibly overly stiff, i.e. too sensitive to details of the mean gradients. What practical implications this may have is however still an open issue. We illustrate this point in the two sections below. Comprehensive confrontation to local gradient-driven and quasilinear predictions is shown in section 2.3.3. How staircase organisation may be seen as a possible route to either confinement improvement or to gyro-Bohm breaking is shown below, section 2.3.2.2.

Taming avalanches: Bohm or gyro-Bohm scaling?

"Life always has a fat tail."

Eugene Fama

Previous sections have established staircase step layers as efficient microbarriers for transport. We now characterise transport in-between these microbarriers and investigate its possible consequences for global confinement. The question is possibly of great significance for the taming or not of avalanches-like transport events may either lead to gyro-Bohm breaking or open routes to improved confinement.

Staircase steps contain the avalanche activity-In-between staircase steps the transport of heat or momentum dominantly occurs through avalanche-type processes. Avalanches lead to strongly non-Gaussian statistics with a marked Lévy-type non-local, non-diffusive character [Dif-Pradalier et al. 2010, Milovanov et al. 2021a]. In the absence of regulating mechanisms for avalanches a departure from favourable gyro-Bohm confinement scaling is thus expected. The E × B staircase organisation is especially relevant because its microbarriers are efficient avalanche regulators. This was noted in Ref. [Dif-Pradalier et al. 2010], Fig. 2 and is further shown in Fig. 2.7.

To illustrate this point, we isolate a 100 a/c s time window during which the staircase organisation is steady. There we compute the Lagrangian space-time autocorrelation A of the turbulent heat flux in two distinct regions of the plasma chosen such that (i) both regions are bounded by steady staircase microbarriers and (ii) the radial width ∆ between these microbarriers is significantly different in both cases. Each of the two regions is a different "tread" of the staircase as defined in Fig. 2.1; their width ∆ is practically measured as the distance between the bounding temperature corrugations. The first region shown in Fig. 2.7-(a) has a width of ∆ = 58 ρ s , the second region in Fig. 2.7-(b) a large width of ∆ = 153 ρ s .

Individually, avalanches are hard to systematically track. Statistically, the radial extent of space-time autocorrelation measures the coherence of radial heat transport. The full width at half maximum of A [black contours] is thus an estimate for the avalanche size in both cases (a) and (b). Interestingly, the full width at half maximum of A closely tracks, in both significantly different cases (a) and (b) the local step size of the staircase. This leads to the following conclusions:

• the staircase step size ∆ (easy to measure) is a good proxy for the radial extent of avalanches;

• avalanche activity is a route to gyro-Bohm breaking: in two distinct regions of the plasma and at the same time, the turbulent heat flux has at least three different radial scalings: two mesoscale ones associated to staircase organisation ∼ 50ρ s and ∼ 150ρ s , in addition to the robust local ℓ c ∼ 6ρ s scale (see Fig. 2.4);

• regardless of the step size (the 'tread width'), avalanche transport fills-in the 'treads'. Without staircase microbarriers avalanches would thus propagate unhindered across larger portions of the unstable plasma domain. As a consequence, stable staircase microbarriers are efficient regulators for the radial propagation of avalanches.

Statistics of staircase step sizes-The relative strength of beneficial microbarriers and detrimental avalanches evolves dynamically, as both interplay: avalanches are dominantly • improved confinement: favoured in the case of strong, steady and closely packed microbarriers. An interesting perspective being that these microbarriers may coalesce [Hornung et al. 2017, Ashourvan & Diamond 2016] into larger scale macrobarrier(s).

Whether observations in Ref. [START_REF] Hillesheim | [END_REF]] are a manifestation of such phenomena is yet unclear;

• degraded confinement: favoured in the case of weak or vanishing microbarriers or when avalanching dominates across a large portion of the plasma.

To investigate the relative strength of microbarriers and avalanches, microbarrier shear strength (as a measure of microbarrier permeability [Dif-Pradalier et al. 2017b, Hornung et al. 2017]) is a useful quantity. From the section above 2.3.2.2 however, statistics of the step size of the staircase ∆ is possibly a more robust observable of microbarrier permeability. Avalanches indeed do not cross staircase microbarriers without perturbing them and possibly destroying them, either way leaving a footprint on ∆. ∆ is thus a measure of the least favourable radial transport scale length in the system: larger values are indicative of degraded confinement and its systematic computation conveys information on the actual transport processes. ⋆ . The statistics of step sizes ∆ is shown in (d). The most probable microbarrier spacing is ∆ stat = 40 ± 2 ρ s . ∆ however can be large, displaying statistics of heavy-tailed distributions: a Fréchet distribution approximates it well, which is symptomatic of extreme value statistics.

Gyro-Bohm scaling through successive staircase steps-∆ can be displayed as a function of ambient gradient drive as well as of ρ ⋆ . Its statistics are shown in Fig. 2.8 and Fig. 2.9. The colourscale in subplots (a) and (b) represent how often a given binned ∆ value is found (the most encountered value is normalised to 100) as a function of ambient drive. Subplot (c) displays its behaviour with respect to ρ ⋆ . The larger ∆ values correspond to temporary destruction of staircase steps, often prior to their reconstruction elsewhere. As compared to our earlier results [Dif-Pradalier et al. 2010, Dif-Pradalier et al. 2015], the present dataset is larger, encompassing steady-state weakly driven plasmas (mostly using thermal bath boundary conditions) as well as more strongly driven (hence less quiescent) plasmas. The rich dynamics of the staircase pattern is now also embedded in the dataset, hence a large distribution of step sizes visible in subplot (d) and a slightly larger most probable step size ∆ stat = 40 ± 2 ρ s (the uncertainty reflects a 95% confidence interval). The most probable staircase step spacing also appears quite independent of the ambient drive-subplots (a) and (b).

Interestingly, despite statistical variability, the ∆ stat scale is approximately constant throughout the scanned ρ ⋆ values. This fact is strongly suggestive of a gyro-Bohm scaling for the heat transport, reintroduced in the system despite avalanching, via the mesoscale and ρ ⋆ independent ∆ stat ∼ 40 ρ s scaling of the staircase step layers. In this regime of parameters (near-marginal, low-collisional), heat transport tends to a gyro-Bohm scaling owing to the successive staircase steps. We note also that the physics of why a gyro-Bohm scaling may here be found, statistically, is significantly different from diffusive/random walk arguments of turbulent transport scaling with the local microscale ℓ c ∼ 6ρ s (see Fig. 2.4).

Extreme deviations, physical interpretation-A well defined scale ∆ stat is observed with yet significant deviations from it. The distribution of step sizes in Fig. 2.8-(d) is characteristic of fat-tailed distributions and is well described by a Fréchet distribution, a special case of Weibull (or Generalised Extreme Value) distributions with lower bound:

F(∆) = τ 1+κ σ e -τ with τ = 1 + κ ∆ -µ σ -1/κ (2.4)
with respectively a scale parameter σ = 10, a location µ = 44 and a shape κ = 0.6. As also shown in Fig. 2.8-(d) a log-normal distribution, obtained through the multiplicative product of many independent random positive variables decays too fast and inaccurately describes the right tail of steps sizes. Fréchet PDF tails naturally emerge from the products of a finite number of random correlated variables and is often used as a phenomenological description of relaxation in disordered systems [Rinne 2008].

The observation of Fréchet statistics illustrates interesting aspects of the physics of staircase organisation. Fréchet (or Weibull) statistics naturally arise in the description of wear and tear of a system that is said to fail when either wear and tear accumulates beyond an acceptable level or if a fatal shock occurs. If the occurrences of fatal shocks can be modeled by a Poisson process whose rate function is state-dependent and if a random process governs the system state, then wear and tear statistics are often well described by Weibull distributions. In Ref. [Milovanov et al. 2021a] and Fig. 2.9, we further discuss implication of Fréchet statistics in the light of the nonlinear Schrödinger equation (NLSE). This point is further addressed in section 2.5.3.

In analogy, if we postulate that the PDF of staircase step spacing emerges from disruptions of individual staircase steps, it is possible to connect the global statistics of ∆ to the statistics of 'wear and tear' of individual steps. The state of individual steps results from a random process: the organisation of turbulence. Through avalanching, a natural statedependent disruptive process continuously happens within the system, the size distribution of avalanches closely obeying Poisson statistics. In such framework, Fréchet statistics for staircase step spacing naturally emerges through combined influence of both the turbulent generation of staircase steps and the disruptive termination by avalanches.

Comparable statistics would certainly be observed if avalanches were to be understood not only as responsible for traumatic disruptions of individual staircase steps but also as active participants in the onset of individual shear layers. In the latter interpretation, both the life-time and the birth-rate [Kosuga et al. 2014] of staircase shear layers is influenced by avalanches and similar Fréchet statistics for ∆ would certainly be expected. Either way, the observation of a Fréchet distribution Eq.(2.4) would then stand for a signature of avalanches in the global organisation of near-marginal plasma turbulence.

Lastly, it could be argued that a power law description could also fit the data in Fig. 2.8-(d) rather than a stretched exponential [Ghendrih et al. 2014a]. A power law would characterise an absence of characteristic size for staircase step spacing. Clearly self-similarity at all scales is not observed but self-similarity could be argued beyond ∆ stat ∼ 40 ρ s as it is known that when accounting for finite size effects in a dataset a power law often crosses over to an exponential decay, hence leading to the curvature in the log-log plot observed in Fig. 2.8-(d). However, with the Fréchet interpretation above, there is little need to invoke finite size effects to explain the fact that ∆ statistics deviates from a power law.

Practical implications for fusion: perspectives for model reduction

The Fréchet statistics above suggest a route as to how a system described by strongly non-Gaussian statistics may display favourable gyro-Bohm like transport provided, as the system increases in size, that enough staircase microbarriers are formed. This could be expected near marginal stability. Inversely, uncontrolled/untamed avalanches at the system size may be expected in small plasmas (large ρ ⋆ ) or strongly driven systems, providing a natural route to gyro-Bohm breaking. One key question for transport prediction thus becomes "what is the distance to (nonlinear) marginality?" 2.3.3.1 What does it suggest for transport scaling?

In addition to microscale regulation via shear suppression [Biglari et al. 1989, Waltz & Kerbel 1995, Hahm & Burrell 1996], transport can thus be directly regulated at mesoscales via modulation of front propagation at the scale of the staircase step width. This emergent mesoscale near marginality is consistent with a transition from Bohm-like (avalanche-dominated) transport for small plasma sizes to gyro-Bohm like (staircase-regulated avalanching) transport for large plasma sizes ρ ⋆ = ρ i /L ≲ 1/300 [Dif-Pradalier et al. 2017b]. Above marginal stability, as staircase organisation fades, avalanching activity is less regulated and Bohm-like confinement can be observed [START_REF] Nakata | [END_REF]. Interplay between layering and spreading, permeability of shear layers (staircase-related or otherwise) to extended transport events are seldom explored though interesting angles to pursue a deeper understanding of transport [Floriani et al. 2013, Dif-Pradalier et al. 2017b, Hornung et al. 2017, Ashourvan & Diamond 2017, Wang et al. 2020]. Analysis of Gysela data predicts an unfavourable scaling of staircase microbarrier strength with decreasing ρ ⋆ [Dif-Pradalier et al. 2017b]. Recent experimental findings (see next section and Fig. 2.17-(e) in Ref. [Hornung et al. 2017]) tend to display a comparable trend. No such behaviour is visible in Fig. 2.8-(c). Further understanding the statistics of ∆ with decreasing ρ ⋆ is certainly a matter of importance to understand the ρ ⋆ scaling of heat transport. Improved confinement may reside in accessing regimes with faster decaying tails for the ∆ [avalanche] distribution, provided the staircase structure survives. An unfavourable behaviour of the tail of the avalanche distribution on the other hand may provide a natural route to gyro-Bohm breaking.

Reduced models: possible shortcomings in near marginal regimes

Staircase organisation appears to be quite robust, numerically observed with adiabatic or kinetic [Rath et al. 2021, Qi et al. 2021] electron responses, in electrostatic or electromagnetic [Rath & Peeters 2022] frameworks. Near-marginal regimes are relevant to large-scale devices such as Iter. Yet, gradient-driven approaches hinder emergence of such structures [Peeters et al. 2016, Dif-Pradalier et al. 2017b] and further differences between flux-driven (FD) and gradient-driven (GD) transport have been documented [START_REF] Nakata | [END_REF]. Whether these observations are of practical incidence in fusion-relevant configurations is non trivial. The matter is important for there are increasing requirements for fast, reduced, yet reliable models to explore the vast parameter space of magnetised plasma turbulence, interpret experimental results and forecast future large experiments such as ITER. Currently, the more advanced reduced models are based on quasilinear theory (QLT) [G. Laval & Adam 2018, A.A. Vedenov & Sagdeev 1962, Drummond & Pines 1962]. With the advent of machine learning techniques, the ubiquitous closure problem of QLT is approached through data-driven techniques that use large-scale databases of first-principles GD computations. Systematic shortcomings, if any, within reference GD strategies are thus likely to be carried over to the reduced models. Comprehensive understanding of discrepancies between first-principles FD, GD and quasilinear approaches is thus important and timely, and the topic of Ref. [Gillot et al. 2022].

To this end, we confront reference results from nonlinear FD gyrokinetics using the Gysela framework [Grandgirard et al. 2016] to state-of-the-art GD nonlinear gyrokinetics and GD quasilinear calculations, using respectively Gkw [Peeters et al. 2009] and Qua-LiKiz [Bourdelle et al. 2007, Citrin et al. 2017]. We further complement the study with a twofold confrontation with the QL transport framework of QuaLiKiz-Jetto and with the nonlinear local transport framework of Gene-Tango. In these computations, N r,eval instances of QuaLiKiz or Gene [Jenko et al. 2000] locally compute at various radii r eval,i (with 0 ≤ i ≤ N r,eval ) flux-surface-averaged transport coefficients which are passed on to the Jetto [M. [START_REF] Romanelli | JINTRAC: A System of Codes for Integrated Simulation of Tokamak Scenarios[END_REF] or Tango [Parker et al. 2018] one-dimensional integrated modelling suites and used to evolve profiles through flux-driven transport equations. After a transport timescale, new local values at each r eval,i location from the evolved profiles are fed to local QuaLiKiz or Gene and the process loops. In practice N r,eval = 50 for QuaLiKiz-Jetto and N r,eval = 9 for Gene-Tango. Main results are: (i) steady-state predictions of fluxes moderately depend on the nature of forcing well above nonlinear threshold; (ii) near marginality however, nonlinear and quasilinear GD models sizeably underpredict turbulent heat transport. Under a driving flux, profiles display 'stiffness', i.e. hover in the vicinity of their near-threshold flux-matching values. Large, hot devices such as ITER are expected to be stiff due to the temperature dependence of the gyroBohm heat flux scaling, making near marginality a regime which models must confront. Proximity to nonlinear thresholds implies additional complexity as it favours secondary pattern formation and mesoscale organisation. Despite this additional complexity, (iii) the underlying assumptions of QLT hold well across nonlinear regimes. We show that transport underprediction rather stems from the choice of closure, i.e. the nonlinear saturation rule. This work stresses the relevance of QLT for model reductions of turbulence whilst providing guidelines whereby reduced models can be improved. The nonlinear and QL approaches tested here are the current workhorse for estimating transport and confinement in turbulent fusion plasmas. This work thus has implications for present-day experimental data analysis and scenario extrapolation for fusion production. Novel saturation rules should strive to incorporate near marginal flux-driven specificities, often dubbed turbulence spreading [Hahm et al. 2005, Singh & Diamond 2020, Dif-Pradalier et al. 2022], transport nonlocality or staircase organisation.

A few words on the modelling frameworks-Models are often categorised on being either 'local' (e.g. 'flux-tube') or 'global'. In the local approximation, mean profiles are piecewise constant and boundary conditions periodic; in the global approach, both assumptions are relaxed. There are documented differences between both approaches but for the present discussion, being either local or global is secondary to the fact of being either gradient-or flux-driven. Large scale mean (equilibrium) gradients, as the main source of free energy will of course contribute to driving meso-and micro-scale dynamics. Micro-and meso-scales back react non linearly on the equilibrium profiles via turbulent fluxes. A central question is whether any scale separation (in time and space) exists between turbulence dynamics at micro-and meso-scales on the one hand and equilibrium scales on the other hand. This question is likely all the more critical close to (nonlinear) marginal stability where meso-scale dynamics is more pronounced, memory of smaller scale turbulent activity being 'stored' in meso-scale alteration of equilibrium profiles.

Importantly, in flux-driven frameworks, no scale separation between equilibrium and fluctuations is postulated which implies that the sources and sinks which drive the system out of equilibrium evolve on length scales coarser than that of the turbulence as well as on slow, adiabatic time scales. A continuum of (turbulent) micro-and meso-scales can thus feedback on meso-to macro-scales. In that respect, flux-driven approaches are necessarily global whilst the inverse is not true. One could thus argue that flux-driven and global gradient-driven approaches should render close results well above threshold. Closer to (nonlinear) turbulence threshold, near marginal regimes are precisely the regimes where global gradient-driven approaches may prove significantly different from flux-driven approaches, for instance through manifestation of Self-Organised Criticality-like phenomena. The curious reader is also referred to e.g. section 4 of Ref. [Sanchez & Newman 2015], which further illustrates these matters.

Flux-driven Gysela resolves ion Larmor radius scale turbulence and collisional transport in global tokamak geometry, spanning from r/a = 0 to r/a = 1.2, a being the minor radius of the torus. A centrally peaked heat source drives a deuterium plasma out of equilibrium. For r/a ≥ 1, a heat sink is progressively applied which allows convergence to a steady temperature profile on energy confinement times. Steady-state and coarse-grained (see below) density n, temperature T , source S and safety factor q profiles from Gysela are shown in Fig. 2.10. Together with the zonal mean shear shown in Fig. 2.11, they are the reference inputs used to initialise all the other codes. On practical grounds, the Gysela profiles need some amount of coarse graining or smoothing before serving as input for either Gkw or QuaLiKiz: this is the essence of the scale separation assumption inherent to the gradient-driven approach. In a gradient-driven framework and even more so in the local limit, profiles are indeed required to be smooth below a cut-off radial scale ℓ ≤ ℓ c = max(λ lin , λ E×B ) whose physics is not included within gradient-driven models.

Here λ lin ≳ 10 -20 ρ i denotes the typical radial extent of unstable growing modes and λ E×B ≳ 10ρ i the typical width of mean E × B shear structures [Dif-Pradalier et al. 2015]. To this end, Gysela observables are both time averaged at steady-state over 30, 000 Ω ci -which is larger than a typical linear growth time and radially smoothed through sliding windows of 60 ρ i to smear out flux-driven specificities in the profiles. Here Ω ci is the ion cyclotron frequency. Sliding averages of 20 ρ i have also been investigated and found insufficient near marginality to smear out memory of meso-scale organisation from the flux-driven Figure 2.10: Radial Gysela profiles of ion and electron temperature, heat source (left), density and safety factor (right) for both near marginal and above threshold cases. The profiles are temporally averaged over 30, 000 Ω ci and radially through sliding windows of 60 ρ i .

framework. This point is further discussed below whilst evoking the role of mean E × B shear.

Gradient-driven models on the other hand (either local or global) exploit to numerical advantage the assumption of a scale separation between a fixed (mean) equilibrium and fluctuations. Sources and sinks are also required in gradient-driven approaches to maintain the system out of equilibrium; the important point is that they depend on the local dynamics of the plasma: maintaining fixed background mean gradients thus acts so as to counter-act part of the natural dynamics of the system (which would naturally relax towards equilibrium), especially hindering spreading and mesoscale dynamics. The appealing semantic contrast between local and global may sometimes cloud the important fact that, especially near marginal stability, what matters is whether these scale separations are postulated or not. Said differently, models are better classified near marginal stability on whether scale separations are postulated or not, whether mesoscale organisation can develop or not and turbulence spreading can occur, unhindered. These questions have been documented both with Boltzmann and kinetic electron responses -see e.g. Fig. 2 in Ref. [Rath et al. 2016], Figs. 5 and 9 in Ref. [Peeters et al. 2016], Figs.10 and11 in Ref.[Dif-Pradalier et al. 2017b] and Refs. [Rath et al. 2021, Qi et al. 2021]. Inclusion of electron dynamics bears additional physics that will be specifically addressed in a coming work. In the present manuscript, we compare as a first step all approaches with Boltzmann electrons.

Gradient-driven nonlinear Gkw exists in both local and global versions; it is here run in its local (flux-tube) setting, solving a limited subset of the whole plasma volume twisting around the torus due to the magnetic shear of the background magnetic equilibrium.

Gkw is compared to Gysela at 3 different locations r/a = 0.3, 0.4 and 0.6. At each of these locations, one first computes the coarse-grained values of the Gysela profiles (as detailed above) at equilibrium, shown in Figs.2.10 and 2.11. These values are used to define the reference background Maxwellian for Gkw at each location; Gkw then evolves the perturbed distribution function with reference to this fixed background. The resolution chosen for all Gkw runs is such that resolved radial and poloidal wavenumbers extend from (k θ ρ i ) min = (k r ρ i ) min = 0.051 to (k θ ρ i ) max = 2.6 and (k r ρ i ) max is ad-justed to match the radial resolution used in Gysela, i.e. 11.3 for the 'above threshold' case and 8.1 for the 'near marginal' one. This amounts to a radial box size which spans 1/(k r ) min = 2 π × 19.6 ρ i ≈ 120 ρ i at the low field side midplane. A large radial box size of ∼ 120ρ i allows one to avoid unnatural interactions with the periodic boundaries for the Gkw computations, in consistency with the local framework, are radially-periodic. Gkw thus effectively computes at each location r/a = 0.3, 0.4 and 0.6 a nonlinear realisation in the local gradient-driven approximation of the dynamics locally expected of the coarsegrained flux-driven Gysela profiles.

Further reducing model complexity, turbulent structures are expected with the quasilinear approximation to bear memory in shape and localisation of their linear generation mechanism. They remain radially thin around a reference magnetic surface and only depend on local plasma parameters and on local gradients. This is direct consequence of the assumption which relates the distribution function fluctuations to the potential fluctuations through local equilibrium parameters. This induces, as with the gradient-driven approach, a spatial scale separation between a local turbulent behaviour, and a slower and smoother evolution of the profiles. Quasilinear QuaLiKiz is run here in 2 configurations: in its standalone version, it is inherently local; in its version coupled to the transport code Jetto, the local transport coefficients from QuaLiKiz are used as inputs for transport equations, thus allowing for the equilibrium to evolve on transport timescales whilst retaining the locality and scale separation assumptions on shorter length-and time-scales. QuaLiKiz solves the gyro-kinetic dispersion relation, here with Boltzmann electrons. For efficiency, the analytic distribution response is simplified by computing the shape of the potential fluctuations in the fluid limit [Bourdelle et al. 2007], while keeping a kinetic treatment of the wave-particle resonance. Stable modes are neglected; unstable modes are accounted for through a double power law in k θ for the turbulent intensity spectrum. The amount of numerical integrations is limited by performing the resonant velocity integration analytically. The effect of zonal flow shear is modelled by perturbative modification of this response. Extensive benchmarks with local gradient-driven gyrokinetic codes have led to refining the closure for the potential fluctuations [Casati et al. 2009, Citrin et al. 2017], calibrated to a database of local nonlinear gyrokinetic simulations akin to those presented here with Gkw.

Finally, two additional series of computations have been performed in 'near marginal' and 'above threshold' settings using the integrated scheme of Gene-Tango. This additional framework complements the Gkw and QuaLiKiz-Jetto approaches: Gene is a gradient-driven local nonlinear gyrokinetic framework akin to Gkw and the transport framework of Tango is similar to Jetto. Comparison between the combined Gene-Tango and QuaLiKiz-Jetto frameworks allows to test (i) the impact of the QL reduction with respect to nonlinear evolution. Since both Gene and QuaLiKiz have scale separations built in their framework, comparison of both approaches to Gysela also allows to specifically assess (ii) the influence of the scale separation assumption. This point is further discussed below and is found to be important near marginal stability.

All 4 approaches can handle complex geometries yet here, for the sake of simplicity, each is set to run in the same simplified toroidal magnetic geometry with circular and concentric flux surfaces

B = (B 0 R 0 /R)[re θ /qR 0 + e φ ]
, where e θ and e φ are the unit vectors in the poloidal and toroidal directions, B 0 is the magnetic field on axis, r the minor radius and R = R 0 + r cos θ the major radius.

Let us close this section with 2 concluding remarks:

• First, it is worth emphasising that the impact of the type of forcing on turbulent transport and achievable steady-states is not restricted to numerical simulations. Experimentally, it is well documented in fluids that the statistics of flow states as well as transitions between them may critically depend on the type of forcing applied -e.g. in von Kármán flows either at constant impeller speed or at constant torque [START_REF] Saint-Michel | [END_REF]]. These observations are analogous to respectively the gradient-and flux-driven cases discussed above. Not surprisingly, the difference between both regimes manifests close to transitions between two equilibrium states while disappearing away from the transition. Again, this behaviour is reminiscent of magnetised plasma dynamics near marginal (nonlinear) stability or above threshold.

• Last, comes the important question of the nonlinear saturation of turbulence in all 3 approaches. E × B shear has at least two components, a mean part and a fluctuating part which may affect and regulate the system differently, over different timescales and through possibly distinct pumping and damping mechanisms. In Gysela both components are self-consistently computed and act in concert. A radial force balance E r -V T B P + V P B T = ∇p/ne is satisfied throughout the computation [START_REF] Dif-Pradalier | [END_REF]]. This is not so for QuaLiKiz, Gkw or Gene for which the mean part is unconstrained. It is an input that can be freely imposed in local approaches without breaking constitutive orderings, echoing the fact that the above radial force balance is not required to be satisfied. Gkw and Gene thus only compute the fluctuating E × B shear as part of the nonlinear response. For QuaLiKiz in its standalone version, the mean E × B shear is an input as well and the effect of the fluctuating E × B flows on the saturation of turbulence is part of the closure scheme, thus approximated from local nonlinear gyrokinetics. In the QuaLiKiz-Jetto framework this must be slightly nuanced: (i) the transport equation leads to an evolving pressure gradient; (ii) toroidal velocity V T can either be prescribed from measurements or self-consistent from the momentum transport equation including external torque and QuaLiKiz momentum transport and (iii) poloidal velocity V P is evaluated from the neoclassical transport model Nclass [Houlberg et al. 1997]. The QuaLiKiz-Jetto interface then estimates E r through radial force balance. From there comes the perpendicular velocity shear input into QuaLiKiz at each radial grid point. Note that a similar procedure exists in the Gene-Tango interface but has not been used for the present study. Note as well that even with this consistent evaluation of radial force balance not all relevant flows are taken into account in these approaches and in particular no structuring at mesoscales can occur due to the built-in scale separation assumptions. In the case of QuaLiKiz there are also additional missing intrinsic rotation or turbulence-generated zonal flows effects.

In the present manuscript, we focus on testing consequences of constitutive assumptions in the models and therefore choose to impose in Gkw and QuaLiKiz the local values of mean E × B shear from Gysela. One could also ask for a complementary approach and e.g. only consider the self-consistent fluctuating E × B shear from nonlinear Gkw or Gene or the shear effectively allowed in QuaLiKiz through the choice of closure and ask how this would compare to Gysela. This is done in the Gene-Tango framework where the mean E × B shear from Gysela is not imposed and only the fluctuating E × B shear from Gene is included. A comprehensive discussion of these issues is beyond the scope of the present work but one would expect discrepancies between models to be especially visible near marginal stability. Indeed, the possibility to self-consistently evaluate the mean E×B shear is directly connected to the role of mesoscale organisation and only fully present in flux-driven Gysela.

To illustrate this point, let us emphasise that computed fluxes from Gkw and especially from QuaLiKiz show significant sensitivity near marginality to imposed levels of mean E × B shear (within one standard deviation) when mean E × B shear pro-files from Gysela are only smoothed over 20 ρ i (approximately twice the mean flow width or the width of mean profile corrugations [Dif-Pradalier et al. 2017b]). This sensitivity is much less pronounced with the 60 ρ i radial smoothing of the Gysela profiles reported here. This is certainly illustrative of the important role that mean E × B shear plays in FD approaches near marginality. Conversely, the fact that QuaLiKiz-Jetto and Gene-Tango provide comparable results with the 60 ρ i radial smoothing despite different imposed levels of mean E × B shows that the strong coarse graining applied to Gysela effectively minimises discrepancies with the other approaches. The near marginal flux underprediction which we report below is thus likely significant and possibly a lower-bound estimate of actual differences.

This question of a consistent evaluation of mean E×B shear/of mesoscale organisation is certainly an interesting one for prospective gradient-driven computations which would aim to provide quantitative answers for new/untested plasma configurations for which there is no a priori "ground truth" (experimental or flux-driven). To try to mitigate this problem local frameworks are thus increasingly coupled to transport equations, which provide a step towards a more self-consistent coupling between mean and fluctuations. In that respect, the lingering discrepancy between Gysela on the one hand and both QuaLiKiz-Jetto and Gene-Tango on the other hand, visible in Fig. 2.13-(c), is certainly indicative of the fact that coupled models where the turbulence model is yet based on a scale separation still miss part of the dynamics near marginal stability. This important observation certainly calls for further studies on the matter, whilst especially relaxing the assumption of Boltzmann electrons.

Two distinct regimes-Two paradigmatic simulation regimes are considered in the electrostatic regime with Boltzmann electrons. Both cases are run in the so-called "local limit" [Lin et al. 2002], at ρ * = ρ i /a ≤ 1/250, where comparison to local Gkw is fair. Coarsegrained (see above) radial profiles of normalised temperature gradients R/L T = -R∂ r T /T and zonal flow shear are plotted in Fig. 2.11 for the 'Above Threshold' (top) and 'Near Marginal' (bottom) cases. Shaded areas represent temporal standard variations. The large deviation from the mean shearing rate in the near marginal case (bottom, right axis) results from the meandering of staircases, which have already been reported to play an important role in this regime of parameters [Dif-Pradalier et al. 2017b]. The linear (black hourglass symbols) and nonlinear (red squares) thresholds are discussed in next section. In order to broadly span parameter regimes, main plasma parameters vary significantly between cases, as illustrated in Table 2.2. This choice echoes the fact that parameter space is broad and parsing it is uneasy. A comprehensive discussion of the precise impact of each parameter: ρ ⋆ , aspect ratio R/a, τ ,. . . is a daunting task. We take a first step in this direction and choose to broadly span parameters from the 'near marginal' case to the 'above threshold' case. One wishes thus to minimise the possibility that found conclusions may strongly depend on the precise corner of parameter space that is investigated. Of course the present approach will require to be complemented by further dedicated studies specially focusing Case ρ ⋆,50 R/a ν ⋆,50 (q 50 ; q 95 ) τ = T i /T e near marginal 1/250 3.2 0.24 (1.4 ; 4.0) 1 < τ < 1.3 above thresh. 1/350 6 0.02 (1.7 ; 2.8) 0.9 < τ < 1

Table 2.2: Main plasma parameters in considered cases. Subscripts 50 and 95 respectively denote parameter values estimated at locations r/a = 0.5 and r/a = 0.95. Our choice of parameters in Table 2.2 is however not totally random as we have performed with Gysela extensive aspect ratio (in the range R/a = 3-10) and ρ ⋆ scans (in the range ρ -1 ⋆ = 190 -380) and found numerically a confinement time scaling law of the form: Caschera 2019]. The 'near marginal' case has both a larger ρ ⋆ and a lower aspect ratio than the 'above threshold' case. Given the above scaling, one may expect for the 'near marginal' case a degraded confinement with respect to what would have been obtained should we have run the 'near marginal' profiles with the R/a and ρ ⋆ parameters of the 'above threshold' case. This is an important point: the parameters chosen for the 'near marginal' case are not such that they strengthen flux-driven specificities. Rather the opposite: beneficial zonal flow activity through mesoscale staircase organisation is indeed observed to be enhanced at smaller ρ ⋆ values [Dif-Pradalier et al. 2010], allowing to recover favourable gyro-Bohm-like confinement scaling by taming avalanching/spreading activity through successive staircase steps (see e.g. section 2.4 in Ref. [Dif-Pradalier et al. 2017b] and Ref. [Milovanov et al. 2021a]). Anticipating on the following sections, the discrepancy found in the near marginal regime between Gysela on the one hand and either Gkw, Qua-LiKiz, QuaLiKiz-Jetto or Gene-Tango on the other hand strengthens the necessity to further understand and incorporate near marginal physics in reduced gradient-driven or quasilinear models of turbulence.

τ c Ω ci ∝ (R/a) 0.88 ρ -2.4 ⋆ [
For each set of local values of the Gysela parameters, two thresholds are to be distinguished: linear R/L lin T represents the normalised temperature gradient above which an unstable mode grows, at vanishing E × B shear. The inclusion of self-generated flow shear Kubo numbers of order unity-QLT is valid [Diamond et al. 2010] in the low Kubo number limit K = τ jump /τ int < 1, ratio of a jumping time τ jump of particles from one turbulent eddy to the next over a nonlinear eddy-particle interaction time τ int . Kubo numbers are estimated from the Gysela flux-driven data in various ways, summarised in Tab.2. Alternatively, the slower evolution of the potential field provides a relevant correlation time τ corr jump for turbulent fluctuations, trade-off between unstable growth and nonlinear saturation. It must be computed as a Lagrangian correlation time, in the co-moving frame of the eddies to correct for the Doppler shift induced by their phase velocity. We estimate it through image registration, following the toroidal shift between 3-dimensional turbulence snapshots of φ. It is compared to turbulence-driven stochastic transport times of particles which, assuming a diffusive ansatz for E × B fluctuations, drift across eddies in about τ diff,x int = L x /⟨|ṽ E,x | 2 ⟩ 1/2 , with x = {r, θ} and L x the transverse correlation lengths computed from Gysela outputs. Interestingly, employing a Eulerian correlation time would result in severe Kubo number overestimation, locally up to factors of 25.

Three Kubo numbers, combinations of the above nonlinear times are plotted in Fig. 2.12 for both 'above threshold' and 'near marginal' regimes. The various definitions for K, coherent, provide the following picture: (i) on the basis of order unity Kubo numbers, QLT should be marginally valid. Yet, as shown below, key assumptions at the heart of the QLT reduction remain valid throughout nonlinear evolution, which strengthens the case for quasi-linear integrated modelling. Interestingly also, (ii) consistently larger K values near marginality stress the more percolative nature of transport there. Avalanching emerges as a key theme to distinguish between 'above threshold' and 'near marginal' regimes as they likely underpin the larger K values computed near marginality. It is a likely indication that incorporating avalanching and its zonal mean flow regulation may significantly alter transport predictions and improve model behaviour near marginality.

Near marginal heat flux under-prediction-Heat fluxes are computed with Qua-LiKiz, Gkw or Gene from the Gysela time-averaged steady-state profiles plotted in Figs.2.10 and 2.11 with the same degree of approximations (electrostatic & Boltzmann electrons). All codes have different normalisations, with Gysela e.g. being normalised to QGys = Q/(n 0 T 0 v T i0 ) with v T i0 = (T 0 /m i ) 1/2 the ion thermal velocity. In Fig. 2.13-(a) and (b), fluxes from all codes are in the units of Q = QGys / QgB , with QgB = ρ 2 ⋆,50 a/R and values of ρ ⋆,50 and R/a given in Tab.2.2. Consistent rescaling factors have been applied to QuaLiKiz and Gkw fluxes. Non-axisymmetric (turbulent) contributions to heat fluxes are shown in Fig. 2.13-(a) and (b). Turbulence spreading and profile corrugations, inherent to flux-driven complexity are absent or hindered in all QL or GD approaches. A fair comparison thus requests that the Gysela reference data be significantly coarse-grained before being handed over to Gkw, Gene or QuaLiKiz as input profiles. This smoothing includes the E × B shear profiles. Hence the time averaging over 30, 000 Ω ci and radial smoothing over 60 ρ i performed on the Gysela observables, as detailed above -to smear out visible flux-driven specificities in the profiles. The stabilising effect of mean zonal flow shear is accounted for in Gkw and QuaLiKiz whilst locally imposing the smoothed reference Gysela mean shear values. In Gene-Tango only the fluctuating E × B shear is included -see discussion above. Sensitivity to gradient fluctuations, inherent to gradientdriven approaches is further assessed by additional scans in R/L T and E × B shear within one temporal standard deviation of the Gysela profiles.

Without inclusion of E × B shear stabilisation [i] QuaLiKiz heat fluxes are overestimated with respect to Fig. 2.13 by over an order of magnitude (not shown here). With the inclusion of shear, [ii] at locations of low or vanishing E × B shear and despite the 60 ρ i smoothing, gradient-driven (standalone) QuaLiKiz commonly displays [subplots (a) and (b)] variations by factors in heat fluxes from one radial position to the next whilst Gkw exhibits much less sensitivity to E × B shear stabilisation. Interestingly, this large sensitivity of QuaLiKiz to shear is mitigated when called within the integrated framework of Jetto [M. [START_REF] Romanelli | JINTRAC: A System of Codes for Integrated Simulation of Tokamak Scenarios[END_REF] to allow for a flux-driven QL profile evolution, driven by a central source that mimics the one of Gysela [subplot (c)]. In the regime above threshold [iii] reasonable agreement in computed fluxes is found across fidelity hierarchy. Conversely near marginality, [iv] despite significant smoothing, heat flux discrepancies in Fig. 2.13-(b) are well outside allowed gradient sensitivity and fluctuation 'error' bars. Secular growth of zonal flows in Gkw is responsible for the observed large flux underprediction. This echoes previous observations casting concern on near marginal gradient-driven predictions [Peeters et al. 2016]. The soundness of separating fluctuations from mean is thus clearly questioned near marginality.

The conclusions above are further confirmed when comparing FD Gysela to the FD quasilinear frameworks of QuaLiKiz-Jetto and Gene-Tango, in the same two regimes. In subplot (c), profiles from QuaLiKiz-Jetto and Gene-Tango are evolved until heat fluxes match the nonlinear Gysela reference fluxes. The figure of merit now becomes how close evolved quasilinear or nonlinear profiles are at flux equilibrium with those of Gysela. Remarkable [v] profile agreement is found for both approaches in the 'above threshold' case, which in the case of QuaLiKiz also echoes the agreement in fluxes displayed in panel (a). Near marginality however, [vi] a large over-prediction of the temperature is required for both QuaLiKiz-Jetto and Gene-Tango to carry the same flux as Gysela. Interestingly, the fact that both quasilinear QuaLiKiz-Jetto and nonlinear Gene-Tango frameworks provide consistent results illustrates the fact that [vii] the QL reduction is not a priori responsible for the observed flux under-prediction in the 'near marginal' regime. A constitutive ingredient is missing in this regime, at least with Boltzmann electrons, which is due to either the local or the gradient-driven approximation. This rather clearly [viii] points towards a problem with the scale separation assumption near marginality and gives a workable route for improvement. Axis for improvement: saturation rules-Flux discrepancies between Gysela and Gkw likely stems from disregarding the feedback of fluctuations on an assumed fixed "equilibrium". This enforces local and single-valued flux-gradient relations, underestimates turbulence spreading, avalanching and mesoscale organisation. All of which contribute to transport, especially near marginality. Discrepancies with QuaLiKiz may either come from violating assumptions central to QLT -linearity of fluctuations-or by inheriting shortcomings akin to those of Gkw, through the choice of saturation rules.

To disentangle these questions, we compute in Fig. 2.14, for all three approaches, the complex argument α n of the n-th Fourier component of the heat flux [panel (a)] -a proxy for the linear cross-phase between electric potential and pressure fluctuations. They depend on the toroidal wave number n labelling each eigenmode. In Gysela, Q Gys Outlook-The success of reduced models especially hinges on the reproduction of nonlinear gyrokinetic fluxes [Citrin et al. 2017]. Flux under-prediction in the dynamically important regime of near marginal stability is thus a matter of importance. At the heart of this paper lies the fact that fkux-driven (FD) and gradient-driven (GD) models provide significantly different flux predictions with Bolzmann electrons close to marginal stability R/L T c , which underpins basic discrepancies in how nonlinear saturation of turbulence is modelled. This should foster renewed interest in ways to complete quasilinear (QL) or GD models near marginality. The robustness of linear features [Besse et al. 2011] across fidelity hierarchy and across turbulent regimes is encouraging perspective and provides constructive directions whereby reduced models could be improved.

Discrepancies between the nonlinear frameworks of Gysela and of both Gkw and Gene-Tango are strong indications that turbulence spreading and mesoscale patterning are in fact central to accurate transport predictions near marginality. Larger Kubo numbers near marginality, as shown in Fig. 2.12 reinforce this point, which is also differently stressed by recent works in the plasma edge [Singh & Diamond 2020, Dif-Pradalier et al. 2022]. The fact that QuaLiKiz behaves better near marginality than Gkw is possibly because the QL closure does not include the long-lived zonal flows that quench transport near marginality in GD approaches [Peeters et al. 2016]. In QuaLiKiz, stable modes are neglected and the turbulent intensity spectrum is fitted [Casati et al. 2009, Citrin et al. 2017] onto databases of GD nonlinear computations, similar to those presented here with Gkw. With this procedure, QL models inherit the shortcomings of the primitive gradient-driven models onto which they are adjusted. The near marginal transport shortfall in QuaLiKiz indeed largely comes from issues with the QL closure, i.e. the choice of saturation rule and not the QL reduction, per se. The QuaLiKiz-Jetto framework provides a step towards coupling fluctuations with mean dynamics. This framework yet emphasises a similar near-marginal flux under-prediction, well reproduced by the similar yet fully nonlinear framework of Gene-Tango. This is further indication that the leading-order problem is likely the assumption of a scale separation between 'equilibrium' mean scales and fluctuating scales, the latter not being able to feed back onto the former. In the case of QuaLiKiz, this scale separation is inherited from the choice of the saturation rule. This provides directions for improvement. In physical terms, near marginal regimes require description of transport below or at linear stability and of possible coupling to modes presently predicted as stable in QuaLiKiz. Importantly, it also requires to model the self-advection (spreading) of turbulent domains and the possibility of non-monotonic flux-gradient relations. Several routes can be explored. In the spirit of current frameworks, QL models could be trained on near marginal flux-driven databases such as provided by the likes of Gysela. This would likely lead to QL closures with regime-dependent turbulent intensity spectra. Alternatively to present closures, QL models could also be coupled to dynamic equations for the turbulence intensity, e.g. in the form of reaction-diffusion [Heinonen & Diamond 2019] or k -ε equations [Baschetti et al. 2021], enriching accessible nonlinear dynamics.

The present work provides a framework of understanding. Ongoing studies are concerned with further characterising the 'above threshold' and 'near marginal' regimes [Garbet et al. 2007] when kinetic features of electron dynamics are present. This is important to assess relevance for Iter extrapolations. Electron dynamics is indeed known to locally modify turbulence organisation near low order rational surfaces [Dominski et al. 2012] yet, interestingly, key features of near marginal turbulence with Boltzmann electrons (flow patterning, shear effectiveness and staircase organisation) -central here to the 'near marginal' regime-robustly endure in kinetic electron regimes [Rath et al. 2021].

Experimental evidence

We here present results from Refs. [Dif-Pradalier et al. 2015, Hornung et al. 2017], which contain the first (to our knowledge) characterisation of the E × B staircase in actual tokamak experiments. This is done by analysing turbulent fluctuations measured by ultrafast sweeping reflectometry [Clairet et al. 2011, Hornung et al. 2013] over extended radial profiles. Whilst many of its expected features are recovered under various experimental plasma conditions in Tore Supra, striking new features appear, potentially of great portent when extrapolating towards Iter-grade plasmas.

Staircases are identified through observation of their successive shear layers. The key signatures, consistently with earlier findings are: a robust pattern of quasi-regularly spaced local reductions of the radial correlation of turbulent fluctuations accompanied by a sudden variation of their tilt across these layers, no correlation to low-order safety factor q rationals and a shear flow layer width ∼ 11 ± 4ρ s that scales with the local sound gyroradius ρ s , close to expected zonal flow values [Dif-Pradalier et al. 2015, Fujisawa 2009].

Two new findings are of significance for basic understanding and future plasma experiments. Staircases seem to depend strongly on the free energy source: observed in nearmarginal ion drift-wave turbulence, their seeming disappearance in electron drift-wave dominated regimes-involving different space and time scales and kinetic properties-challenges our understanding of the basic principles of drift-wave self-organisation. Also, a striking increased permeability of transport barriers with plasma size ρ -1 ⋆ = a/ρ s is found, a being the minor radius and ρ -1 ⋆ being akin to a Reynolds number for neutral fluid turbulence. Of concern for large experiments, this observation suggests an unfavourable scaling of transport with ρ -1 ⋆ and opens new routes to understanding gyro-Bohm breaking through microbarrier permeability. These observations at last provide a relatively inexpensive and reliable means to characterise an otherwise difficult measurement: flow shear with high space and time resolution. The method has broad applicability to characterising barrier onset and dynamics, transport and its many forms of bifurcation to regimes of improved confinement.

The experimental setup

The results below are obtained in low confinement (L-mode) plasmas of the Tore Supra tokamak [Dumont et al. 2014] of major radius R 0 = 2.38 m, minor radius a = 0.72 m and magnetic field on axis B 0 = 3.8 T. Turbulent fluctuations are diagnosed with an ultrafast sweeping microwave reflectometer [Clairet et al. 2011]. The probing beam is launched from the tokamak midplane. Reflections from electron density fluctuations at successive plasma radii allow us to reconstruct radial profiles of density fluctuations with very high spatial (< 1 mm) and temporal (3 µs) resolution. High fluctuation levels characteristic of the plasma edge [Gerbaud et al. 2006] and magnetohydrodynamic (MHD) activity in the plasma coreconcentrated in the vicinity of low-order rational q surfaces-can affect this reconstruction. Staircase signatures are therefore looked for in MHD-free plasmas and in experimental conditions where fluctuation levels remain approximately constant, with values below 2%, so that the reflectometer response remains well-defined [Hornung et al. 2013]. The bulk of our observations thus corresponds radially to regions r/a ∈ [0.5, 0.8].

The shear flows at the staircase steps are detected through local contractions of the coherence γ 2 of adjacent plasma layers, indicative of local reductions in the size of turbulent structures [Nazikian et al. 2005] and defined as γ 2 (r, ∆r, F ) = ⟨P σ1σ2 ⟩ 2 /⟨P σ1σ1 ⟩⟨P σ2σ2 ⟩ with P σ1σ2 , P σ1σ1 and P σ2σ2 the frequency cross and auto power spectra computed from the complex reflectometer signal time series σ 1 (r, t) and σ 2 (r + ∆r, t), r the radial position, t the time and F the fluctuation frequency of the signal. The spectra are computed using a 256 points sliding time window with 50 % overlap. Physically, γ 2 is a cross-correlation function in the frequency domain from which a radial profile of mean coherence length ⟨L coh ⟩ F is derived, ⟨•⟩ F denoting the average over the turbulent frequencies F ≥ 15 kHz [Maget et al. 2006]. ⟨L coh ⟩ F is a proxy for the correlation length of the turbulence. Importantly, its shape is robust regardless of the exact averaging domain in frequency F or of its definition using the full width at half maximum of ⟨γ 2 (r, ∆r)⟩ F , as displayed in Fig. 2.15, or an e-folding decay length.

Staircase evidence

The E × B staircase is expected to imprint on the mean coherence ⟨γ 2 ⟩ F its quasiregular structure, at the location of its shear layers [Dif-Pradalier et al. 2015]. Quasiregularly (e) show the real part of the normalised cross-correlation function C (r, ∆r, ∆t), computed from the complex reflectometer signals σ 1 (r, t) and σ 2 (r+∆r, t+∆t) at positions r/a ∼ 0.63 and 0.67. Interestingly, the tilt reverses sign whilst crossing the S 2 coherence minimum, as is expected from the presence of a shear flow at S 2 . In order to test the radial variation of tilting asymmetries against the local minima of the mean coherence length ⟨L coh ⟩ F , we define a scalar measure of the correlation anisotropy φ(r) = i,j C (r, ∆r i , ∆t j )M ij , with M ij = {-1, 0, 1} for ∆r i ∆t j {<, =, >} 0, respectively. Positive time delay correlations as in Fig. 2.15-(d) lead to φ > 0 whereas cross-correlation functions that are tilted in the opposite direction as in Fig. 2.15-(e) lead to φ < 0. A remarkable correlation between ⟨L coh ⟩ F minima and the radial change of sign of φ(r) is evidenced in Fig. 2.15-(c). The pattern of local reductions of fluctuation correlations at S 1 -S 3 mirrors extrema of mean flow shear, providing a likely visualization of E × B staircase steps.

An intriguing feature of the E × B staircase revealed by numerical simulation is the propensity of its shear layers to meander, river-like with their radial locations remaining at roughly constant drive of the turbulence ∇p/p [Dif-Pradalier et al. 2015]. This idea is tested here, within the large uncertainties that affect the measurement of the pressure profile by comparing the mean coherence lengths measured in a reference ohmic phase to a phase where 1.6 MW of additional hybrid power is injected. The power is centrally deposited (within r/a < 0.5) resulting in a steepened pressure gradient in the central region. The estimated radial inward drift of the isovalues of ∇p/p in Fig. 2.16-(a) is about r/a ≈ 0.02, which is in the same ballpark as the observed drifts of the local minima in Fig. 2.16-(b).

Figure 2.17: Staircase properties. Quantities h s , δ flow and ∆ that characterise the E × B staircase and its microbarriers, plotted against key plasma parameters: ρ ⋆ and safety factor q.

2.4.0.3 Barrier permeability, impact on confinement Based on the above results, an extensive database search has led to identifying over 200 coherence length minima displaying a staircase signature, encompassing different plasma conditions, including Ohmic and externally heated plasmas. The bulk of staircase signatures is observed in the range 0.5 < r/a < 0.8. Below r/a ∼ 0.5, both MHD activity near low-order safety factor q rationals and measurement accessibility at operational field and density values limit the database. The outer region r/a > 0.8 is discarded as the assumption of a linear response of the reflectometer there is less well established. Consistently with predictions, staircases are observed in low to moderate collisionality regimes ν ⋆ ∈ [0.01; 4] and in the vicinity of turbulence threshold: assuming T e = T i , the measured electron density L n = -n(dn/dr) -1 ∼ 0.2 m and temperature L Te = -T e (dT e /dr) -1 ∼ 0.1 m gradient scale lengths correspond to a moderate turbulent drive η = L n /L T ∼ 2 -3. Important features of the E × B staircase can be quantitatively evaluated from the mean coherence length profiles, as illustrated in Fig. 2.17. (i)-Linearly, modes localise near low-order q rationals. The plasma there is also prone to growing magnetic islands. These effects may result in localised profile corrugations and local reductions of coherence. No clear correlation is found however between staircase step lo-calisation and q rationals evaluated by EFIT [Lao et al. 1990], suggesting as in calculations [Dif-Pradalier et al. 2015] a nonlinear turbulence-borne origin for the staircase. The drop of the mean coherence length h s at the location of the staircase shear layers is indicative of the strength of flow shear there and measures the permeability (porosity) [START_REF] Ghendrih | Shearing effects on density burst propagation in SOL plasmas[END_REF] of the successive micro-barriers [details below]. Interestingly also, no clear correlation is found, Fig. 2.17-(b) between flow shear strength h s (barrier impermeability) and low-order q rationals. A synergistic reinforcement of a staircase shear layer in proximity of a low-order q rational as noted in Ref. [Dif-Pradalier et al. 2015] is possible but does not stand out here as statistically significant-the role of magnetic fluctuations or increased zonal flow inertia with kinetic electron response may be key.

(ii)-The width δ flow of the mean coherence length minima ("flow thickness") is similarly a measure of the radial extent of the shear flows at staircase steps and likely influences barrier permeability too. We find δ flow ∼ 11 ± 4 ρ s from Fig. 2.17-(c), a value in close agreement with the predicted radial extent of shear layers in simulation [Dif-Pradalier et al. 2015] and zonal flow thickness in experiments [Fujisawa 2009]. This observation reinforces our assessment of the series of ⟨L coh ⟩ F local minima as successive layers of mean zonal flows. (iii)-Within the present E×B staircase understanding, regions in-between the shear layers have a typical radial scale ∆ (the "avalanche size") and are regions of strong avalanchemediated transport. Two characteristic scale lengths are thus predicted for the system: one ℓ c ∼ 5ρ s is associated with the local (micro) scale of the turbulence, the other ∆ is mesoscale-in-between ℓ c and the macroscale a-and a measure of the typical radial extent of avalanches [Dif-Pradalier et al. 2010, Kosuga et al. 2014]. Interestingly, the existence of these two scale lengths is quantitatively consistent with recent observations from Doppler reflectometry [START_REF] Hennequin | [END_REF]]. As opposed to Bohm scaling, unfavourable when extrapolated to Iter, so-called gyro-Bohm scaling predicts a transport at the scale of ℓ c , i.e. that scales with ρ s , based upon the assumption of a local nature for turbulence. In computations, as the system size increases (as ρ ⋆ decreases) the typical scale of avalanches is found not to increase in size and ∆ ∼ 30ρ s , statistically [Dif-Pradalier et al. 2010]. Thus, even though transport is strongly nonlocal, avalanches rarely propagate throughout the whole system as they are stopped by the narrow staircase shear layers that act as a set of semi-permeable transport barriers. Despite avalanching, the staircase microbarriers are thus expected to reintroduce a favourable gyro-Bohm-like scaling for confinement. This ideal picture is somewhat questioned in Fig. 2.17-(d) for the (circled) small values of ρ ⋆ < 1/350 as no saturation at mesoscales is observed and ∆ can be very large. (iv)-An increased permeability of the microbarriers with ρ -1 ⋆ may explain this detrimental ∆ scaling, a proxy for permeability being a combination (h s ) α (δ flow ) β , with α, β ≤ 0. Both quantities h s and δ flow [Figs.2.17-(c) and (e)] are decreasing functions of ρ -1 ⋆ , indicating a decreasing barrier strength with increasing plasma size. As a result avalanches, less hindered at small ρ ⋆ may propagate on larger distances, as illustrated in Fig. 2.17-(f). This barrier permeability dependence with ρ ⋆ should be better explained based on a weakened drive of the staircase shear layers rather than on an increased damping, collisionality at small ρ ⋆ varying weakly throughout our database. These results could offer a new route to understand the observed breakdown of gyro-Bohm scaling, of concern for future large machines as ρ Iter ⋆ < 1/500.

Multi-channel staircasing?

What organisation emerges from the competition of several instability channels, each possibly leading to a distinct staircase patterning? This question is generic, turbulence generally resulting from a mixture of interplaying free energy sources rather than from a single in- Coupled analysis of high-resolution density and electron temperature fluctuations-the latter unfortunately unavailable on Tore Supra-in the spirit of the above discussion should allow for renewed insights in fundamental aspects of drift-wave turbulence self-organisation. Overall, more experiments are required to fully document staircase regimes. This is a topic of practical interest for near-marginal regimes are relevant to large-scale devices such as Iter. Clearly, more experiments are required to fully document staircase regimes and the practical importance of mesoscale organisation of microbarriers.

Theoretical understanding -the quest for overarching principles

"Of course, Michelangelo's David was wholly contained within a block of marble, which did not, however, provide great insight into what could lie beneath."

Victor A. Coonin

An important question is certainly what is the critical parameter for staircase formation? It is known that for the onset of turbulence in a pipe via directed percolation [Pomeau 1986] -a process akin to avalanching-the Reynolds number is the relevant parameter. Interestingly, we showed that collisionality ν ef f (as expected) damps the staircase shear layers -see e.g. Fig. 8 in Ref. [Dif-Pradalier et al. 2017b]-and that fewer of such structures are unambiguously observed in the edge, i.e. past ρ = r/a ≥ 0.7, or more generally that staircasing is feature of near-marginal stability -see e.g. Fig. 4 in Ref. [Dif-Pradalier et al. 2017b]. Thus the heat drive Q is a key parameter, or more precisely Q turb /Q neo . By analogy, one might argue that Q turb /(Q neo ν ef f ), the ratio of heat flux (drive) over collisionality (damping) is an important discriminating ratio. However, other mechanisms may well be at play and if collisionality is an obvious control parameter, the rationale for the heat flux or Q turb /Q neo as the other key control parameter is less clear. Furthermore, it is interesting to note however that some flow structuration in the edge E r well of Jet [START_REF] Hillesheim | [END_REF] or diiid [Ashourvan et al. 2019] have been interpreted as possible staircases. There, collisional dissipation is high and heat fluxes are low; should these observations indeed be staircases, this certainly points towards the necessity for different or additional mechanisms. Much work is still needed to better apprehend (and then possibly control) what would the critical parameters for staircase formation be.

Several models have been proposed to account for staircase formation and are (very) briefly reviewed below. This understanding is also still in its infancy yet is critical to gain understanding of these phenomena and to build reduced models with confidence -which section 2.3.3.2 has shown to likely be important.

Time delay and (traffic-like) jams

The approach through Eqs.(2.2) and (2.3) when A ≫ ℓ c is sometimes referred to as "strong" nonlocality in contrast to "weak" nonlocal approaches where a partial differential equation keeps relating fluxes and gradients as recently emphasised through a variation on the telegraph equation [START_REF] Gürcan | Transport of radial heat flux and second sound in fusion plasmas[END_REF]. For an overview of nonlocality, see e.g. the interesting reviews [Sanchez & Newman 2015, Ida et al. 2015, Hahm & Diamond 2018]. The rationale here is the following: a change in the thermodynamic force/deviation from marginal profile ∇δT will drive a flux:

∂ t δT + ∂ x Q{δT } = 0 (2.5)
which must remain invariant under the transformations x → -x and δT → -δT . This is usually referred to as the Joint Reflexion Symmetry principle [Hwa & Kardar 1992, Diamond & Hahm 1995]. The simplest mean average flux which satisfies JRS can be written as:

Q 0 {δT } = a 2 (δT ) 2 -b ∂ x (δT ) + c ∂ 3 x (δT ) (2.6)
Injecting Eq.(2.6) into (2.5) leads to a Burgers equation for the deviation from marginal profile:

∂ t δT + a δT ∂ x δT = b ∂ 2 x δT (2.7)
which essentially states that the heat transport processes are dominated by large avalanches. However, this is not the more general form for the heat response. A change in the thermodynamic force/deviation from marginal profile ∇δT will drive a flux, which may indeed have two components: an in-phase one modelled by Q 0 in Eq.(2.6) and an out-of-phase flux coming from the time delay τ which accounts for a non-Markovian temporal response:

∂ t Q = -τ -1 (Q -Q 0 {δT }) (2.8)
Microscopically, this change is mediated by the local mixing time of the distribution function [or of PV in the GFD context]. A finite time τ is associated to this mixing activity. Eq.(2.8) states that the flux will naturally respond with a finite time delay commensurable with this mixing time τ . Eliminating the heat flux between Eqs.(2.8) and (2.5) leads to a telegraph equation for the deviation from marginality:

∂ t δT + a δT ∂ x δT = b ∂ 2 x δT -c ∂ 4 x δT -τ ∂ 2 t δT (2.9)
For sufficient time delay between temperature modulations and local heat flux, heat flux jams are predicted, leading to the growth of patterns which resemble shock trains -a point well emphasised in the context of road traffic [Flynn et al. 2009]-from which develop the staircase 'jumps' [Kosuga et al. 2013, Kosuga et al. 2014, Qi et al. 2022]. Sharp curvature of the temperature will drive E × B shear, through radial force balance. Interestingly, when corrugations are envisioned as "jamitons" [Kosuga et al. 2013, Kosuga et al. 2014], they naturally depend on the strength of the dissipative mechanisms, which echoes the fact that increased collisional dissipative processes are indeed responsible for the weakening of the corrugations. This approach emphasises the plasma as a heat engine and build on the analogy between flux jams and traffic jams. The key role of time delay is analogous to that of driver reaction time. A short reaction time maintains smooth traffic flow, and longer reaction times are more likely to trigger jams. Interestingly, long delay times are to be expected near criticality, consistent with the observation of staircases near marginal stability. This idea is highlighted in Fig. 2.19 where the turbulent heat flux and the temperature gradient at ρ = 0.3 close to the source are plotted against time. As expected the turbulent heat flux responds to the source-driven temperature gradient with a comparable dynamics. Interestingly, this response is not instantaneous but involves a quasi-constant time delay between flux and gradient. The correlation analysis between the two signals gives τ ≈ 30 a/c s ≈ 6900 ω -1 c . This approach Eq.(2.9) involving a "weak" nonlocality does not seem to hold far from the source region where the strong avalanching discussed above takes place as described through Eqs.(2.2)-(2.3). Other mechanisms must thus likely be invoked to account for staircase formation.

Interplay of 2 scales, one emergent

Let's consider the following generic setup, relevant to planetary oceans: there is a uniform (density) gradient, say along z. The fluid is stirred (or is it shaken?

). What can be expected for the distribution of salt s (or heat)?. A robust outcome is a series of well-mixed layers with sharp interfaces, at steady-state. An interesting framework of understanding has been laid out in Refs. [Phillips 1972[START_REF] Posmentier | The Generation of Salinity Finestructure by Vertical Diffusion[END_REF]]. Let's consider the one-dimensional salt diffusion equation along z: ∂ t s = -∂ z Γ s where the flux of salt Γ s is a function of the local gradient of salt: Γ s = f (∂ z s). Combining these two equations leads to:

∂ t s = - ∂Γ s ∂[∂ z s] ∂ 2 z s (2.10)
Immediately, a condition for instability reads: ∂Γs ∂[∂zs] > 0, which is often referred to as antidiffusion. Implementation of a framework which displays anti-diffusive behaviour -where the flux-gradient relation follows Fig. 2.20-leads to 'inhomogenous mixing' and to layering. The model described by Eq.(2.10) provides a good starting point to describe a mechanism for layering. There are however a few issues with this formulation [Pruzina et al. 2022]:

• the problem is mathematically ill-posed for the growth rate of fluctuations diverges as the wavenumber of the perturbation increases; • one possibility to regularise this behaviour is to introduce a time delay in the adjustment of the flux to changes in the gradient. This is much in the spirit of the previous section; it may remove the 'ultraviolet catastrophy' of high wavenumber blow-up but since the model does not provide a mechanism to arrest the steepening of the interface (dotted line in Fig. 2.20) solutions still develop discontinuities in finite time;

• it relies on specifying the buoyancy flux function, and hence cannot give a full description of the physics from first principles, nor of the intricacies of layer evolution. For example, Eq.(2.10) takes no account of the impact of the buoyancy on the velocity field even though clearly the latter is key to the global evolution of the system and the theory should at least try to provide constrains on achievable velocity field configurations;

• last, it provides no physics-based understanding as to why the flux-gradient relation should indeed behave as drawn/postulated in Fig. 2.20.

A solution to the above issues was put forth in Ref. [Balmforth et al. 1998] extending on the classic mixing length model of Prandtl and Kadomtsev to propose that the formation process is due to the interplay of mixing on two length scales ℓ [∂zs,e] , one of which is emergent -i.e. depends on both the local profile scale length ∂ z s and the turbulence intensity e:

∂ t s = ∂ ∂z ℓ [∂zs,e] e 1/2 ∂ z s (2.11)
The model is also accompanied by a dynamic equation for the energy production term, involving scales ℓ. Classic emergent scales in geophysics are the Rhines scale [Rhines 1975] of drift-Rossby turbulence, which is defined by the cross-over between the wave interaction mismatch frequency and the decorrelation rate characteristic of the inverse cascade or the Ozmidov scale [Ozmidov 1965], defined as the characteristic size of the largest eddy that is not significantly affected by buoyancy in a stably stratified fluid. More specifically, in the nonlayered initial state, the relevant lengthscale ℓ 0 is that of the turbulent eddies induced by the stirring. As a result of the stirring, regions of larger buoyancy gradient ∂ z s appear and a shorter lengthscale emerges, typically at the Ozmidov scale ℓ Oz in the problem above. A natural route that was followed by Ref. [Balmforth et al. 1998] is to propose that lengthscale ℓ in Eq.(2.11) interpolates between these two scales and thus appropriately transitions from the stirring length when stratification is weak (∂ z s ≪ e) to the Ozmidov scale when ∂ z s ≫ e:

1 ℓ [∂zs,e] 2 = 1 ℓ 2 0 + 1 ℓ 2 Oz (2.12)
The interplay of processes on two length scales (one of which is dependent upon the gradient scale) triggers self-sharpening of modulations, leading to a staircase. The two scale mixing process is amenable to representation by reduced models of the k-epsilon variety obtained in the fusion context from the starting point of the Hasegawa-Wakatani model [Ashourvan & Diamond 2016, Ashourvan & Diamond 2017, Guo et al. 2019], by wave kinetics [Garbet et al. 2021], and by mean field theory [Yan & Diamond 2022]. Formally, the staircase structure emerges from a boundary value problem due to the growth, coarsening, and saturation of modulations [Malkov & Diamond 2019].

Long-time dynamics observed from this modelling approach often features a collapse of the staircase pattern into a macro-barrier on logarithmic timescales through a succession of merger events [Ashourvan & Diamond 2017, Pruzina et al. 2022]. For long times t, a general law describing the evolution of the number of layers N as 1/N ∼ ln t suggests a self-similar structure to merger dynamics and a link to Cahn-Hilliard models of layering [START_REF] Kawakatsu | [END_REF], Balmforth & Young 2005]. Producing evidence of such dynamics in actual experiments or in first-principles numerical experiments (though challenging observations, for different reasons) would be a major step forward as it would start to address the key question of control of this structure for actual confinement purposes. It could also well be, on the other hand, that such long-time evolution is hardly accessible and thus of little practical interest owing to the intrinsic dynamic nature of the staircase (which may keep being destroyed and reformed) in an actual fusion plasma. This is however a central question for both theory and observations.

Nonlinear Schrödinger equation with subquadratic nonlinearity

Theoretically, most notions of staircase formation follow from the two-scale mixing or the jam concepts. Not all do however; the proposal in this section is such an example, based on Ref. [Milovanov et al. 2021a]. The idea here is to envision the staircase-bound jet zonal flows as a wave packet of coupled nonlinear oscillators for which a complex time-and wavenumber dependent wave function can be introduced. The nodes of the network are the zonal flow layers (the staircase steps) and the springs are the avalanches which allow for a dynamic propagation of information. Through a mean-field approximation the complex wave function obeys a sort of discrete nonlinear Schrödinger equation (NLSE) with dispersive linear term and subquadratic power nonlinearity:

iℏ ∂ψ n ∂t = E n ψ n + β|ψ n | 2s ψ n + V (ψ n+1 + ψ n-1 ), (2.13)
Description through a NLSE naturally emerges when considering a system of coupled oscillators with preservation of total energy. It is important to stress however that exponent s is usually equal to unity when it describes many-body bosonic systems. In particular, the quadratic nonlinearity of the NLSE is an exact representation of Bose-Einstein condensates in the dilute limit [Lieb & Seiringer 2002]. Here however 0 < s < 1, i.e. the nonlinearity is sub-quadratic (SQ-NLSE). Mathematical foundations of this type of equation have been earlier investigated in Refs. [Milovanov & Iomin 2014, Milovanov & Iomin 2019]. It is found in the backbone approximation that the subquadratic power leads directly to white Lévy noise, and to a Lévy-fractional Fokker-Planck equation for radial transport of test particles (via wave-particle interactions). The latter equation was originally introduced in fusion plasma based on heuristic arguments, while its derivation from first principles had not been obtained at the time (see Ref. [Milovanov & Rasmussen 2018] for full discussion). This is an important aspect of the theory: fluctuations generated from vibrations of the coupled oscillators with a sub-quadratic nonlinearity 0 < s < 1 naturally produce a noise which obeys Lévy statistics whereas the same fluctuations with a quadratic nonlinearity give rise to Gaussian statistics. As avalanches are the channel through which, in this interpretation, the staircase flow layers exchange masse, energy or momentum, having their statistics correct is absolutely key. Lévy-type statistics have been shown to fit well nearmarginal transport dynamics dominated by the interplay of avalanches and staircase shear layers in flux-driven gyrokinetics [Dif-Pradalier et al. 2010]. This is a step towards justifying the use of a SQ-NLSE to describe near-marginal plasma dynamics.

Importantly, the sub-quadratic nonlinearity allows for a 'weak localisation' of the plasma avalanches -this is the translation of the staircase layers being 'permeable' or 'leaky' transport barriers [Hornung et al. 2017, Dif-Pradalier et al. 2017b]. In a self-consistent description, the avalanches, which are driven by white Lévy noise, interact with jet zonal flows to form a system of semi-permeable barriers to radial transport. In this interpretation, the staircase motif operates as a marginally stable, dynamic fluctuating system "at the edge of its localisation-delocalisation transition" [Milovanov et al. 2021a]. Said differently, the staircase pattern is the plasma realisation of an "edge of chaos" where avalanches may both enhance or weaken the subtle & dynamic equilibrium of the permeable barriers. The fact that avalanches are weakly localised likely leads to an increased probability for large amplitude and extreme events. Interestingly, in this interpretation, the role of avalanches is dual: they are as much a driving force towards stability -a means to re-establish dynamic equilibrium of the whole pattern when destabilised-than they are a driving force towards destabilisation (which is the usual interpretation, especially when large enough).

At this transition point between localisation and delocalisation, the event-size distribution of plasma avalanches is shown to be a power-law:

w τ (∆n) ≈ ∆n -τ , with drop-off exponent τ = √ 17 + 1 2 ≈ 2.56 (2.14)
This fancy number (see Sec.V of Ref. [Milovanov et al. 2021a] for a more complete discussion) is an exact result of the NLSE model and one of the main outcomes of the theory. This result should not be overlooked. Many of the theories described above (based on either the time delay or the emergent second scale) are yet to produce quantitative predictions for important quantities such as the step spacing of the staircase, the event-size distribution of avalanches or the distribution of time delays between fluxes and gradients. The above result Eq.(2.14) is actually one of the few quantitative predictions that can be tested against first principles. One can indeed translate τ into a κ NLSE = √ 17+1 8 ≈ 0.64 parameter which describes the shape of the event-size distribution of avalanches in the entire range of ∆n variation. Earlier Gysela results found a best fit using Fréchet distributions -see Eq.(2.4) and Eq.(3) in Ref. [Dif-Pradalier et al. 2017b] which rendered κ Gysela ≈ 0.67, as shown in Fig. 2.9.

Both results predict a significant power-law tail to the event-size distribution of plasma avalanches and are remarkably close. This strongly suggests that tokamak plasmas in L-mode may produce large-amplitude and extreme avalanches with much greater a probability than what one would expect based on a Gaussian statistical paradigm. From a statistical perspective, this behaviour is compatible with dynamics of a self-organised criti-cal state [Bak et al. 1987]. Notwithstanding, the critical exponents pertaining to this state are found to be inconsistent with classic models of avalanche transport based on sandpiles and their generalizations (e.g., Ref. [Milovanov et al. 2021b] for review), suggesting that the coupled avalanche-jet zonal flow system operates on different organising principles. Also the results indicate that the plasma staircase could be both an agent to monitor the evolving plasma conditions in a tokamak -implications of this finding for "contactless" plasma diagnostics using ultrafast-sweeping reflectometry have also been briefly discussed in Ref. [Milovanov et al. 2021a]-and an interesting route forward for improved confinement (see section 2.6 below) where the intrinsic instability of macro transport barriers would give way to a dynamic array of resilient and semi-permeable weak transport barriers.

From a theoretical perspective, it should be noted that it would be of great significance if the use of SQ-NLSE could be justified on firmer grounds -i.e. if a NLSE with subquadratic nonlinearity could be derived from first principles within a justifiable corner of parameter space and under an explicit series of testable hypotheses.

(Some) open questions

"Quand les talons claquent, la cervelle se vide"

Maréchal Liautey

Let us conclude with a summary of observations/thoughts to be tested or expanded in forthcoming works.

mechanism, or rather: mechanisms? the obvious open question which looms large is that of the understanding, prediction and possibly control of how staircases form, i.e., how is 'inhomogeneous mixing' realised. Staircases, ubiquitous in nonequilibrium dynamics are intrinsically fascinating for one expects turbulent mixing to tend toward homogenisation [Rhines & Young 1982] or quasilinear flattening [Sagdeev & Galeev 1969] of profiles. This question has been somewhat addressed in the previous section 2.5, is still developing but is clearly not settled. The constitutive mechanism(s) are likely not unique: several routes can indeed be at play and lead to staircasing, possibly within the same system. For instance in MFE, avalanches may play a dynamic role and either of the 3 mechanisms discussed above may be of theoretical value; in other cases, impurities coupled to interchange turbulence may also drive staircasing. The underlying physics might be quite different then as it resembles double diffusive convection (DDC) which occurs in stellar interiors (DDC with magnetic buoyancy) [Hughes & Brummell 2021], oceans (thermohaline DDC)-especially the polar oceans [Radko 2022], and planetary atmospheres;

single or multiple barriers? plasma "jets" spontaneously form at the staircase steps, uncorrelated to low-order q rationals and are best envisioned as weak transport barriers whose strength increase with proximity to criticality. Transport is likely subdiffusive across the jets and fundamentally nonlocal, non-diffusive and avalanchemediated elsewhere. The flow pattern is key to arresting the expansion of nonlocality to macroscales and containing it to mesoscales. Depending on mesoscale flow shear strength, it may also break gyro-Bohm transport scalings. Assessing the viability of E × B staircase states (or more generally of several permeable/weak transport barriers spread out thought the bulk of the confined plasma) as alternative regimes of enhanced confinement (which in comparison with more familiar transport barrier scenarios usually involve a single, larger and stronger barrier, possibly prone to instability) is a timely question. As a matter of comparison, if the plasma staircase proves in tokamaks even remotely as important in terms of transport regulation as it is known to be in our atmosphere, this question could prove fecund;

account for mesoscales, spreading & avalanche dynamics in reduced models? transport organises in a continuum of scales yet mesoscales ∆: ρ i ≪ ∆ ≪ L appear to play a central role as staircase organisation spontaneously separates distinct transport regimes into distinct physical regions -patches of strong turbulent mixing with systematic redistribution of heat and angular momentum are dynamically bounded by sharp transitions (with different porosities) that precludes, hinders or mitigates the systematic distribution between adjacent patches. An important question becomes whereby to understand and predict the (distribution of) relevant transport (meso?)scales that the turbulent organisation would select. Flux-driven systems often sit near marginal stability. Models which aim at providing a predictive viewpoint on turbulent transport should likely endeavour to factor in aspects of mesoscale organisation, stiffness, spreading, secondary or tertiary flow organisation into their modelling framework, aspects which are often precluded or hindered [Peeters et al. 2016, Dif-Pradalier et al. 2017b]. There are indications indeed that wave-breaking and positive eddy feedback on the mean fields are critical, naturally lead to a strong dynamical irreversibility which challenges linearised wave theories and models with scale separations -see e.g. an interesting discussion in Ref. [McIntyre & Palmer 1985], p.967 and recent observation of significant discrepancies in heat transport predictions between flux-driven approaches and gradient-driven or quasilinear approaches [Gillot et al. 2022] which are reported near marginality. These observations may have important consequences for reduced modelling yet are still in their infancy and require further investigation.

In that spirit, validity of gradient-driven approaches is debated, especially at the plasma edge -where scale separations are questionable, or near marginal stability -i.e. between linear and nonlinear thresholds, in the so-called "Dimits shift" region [Rosenbluth & Hinton 1998, Dimits et al. 2000]. Flux-driven approaches indeed show near marginality that distinct turbulent states and flux levels can be obtained from almost identical underlying mean plasma gradients. This fact underlies a multivalued (bistable mostly), non monotonic flux-gradient landscape [Diamond et al. 1997 transport scaling & the Bohm-gyro-Bohm "transition": transport scalings may require reinterpretation because of mesoscale organisation. Indeed, in strongly driven systems, staircase organisation vanishes. Avalanches are not tamed/regulated at mesoscales by weakened/absent staircase shear layers and Bohm-like transport is reported. Conversely, near marginal stability staircase organisation regulates (sta-tistically) avalanches at mesoscales ∆ ∝ 40 -50ρ i . ∆ is independent of the system size ρ ⋆ which naturally leads to a favourable gyro-Bohm like ρ ⋆ -scaling for transport: as the system increases in size, the confined core radially accumulates more and more staircase shear layers in its bulk. The relevant transport scale becomes the staircase step size, precisely because the inter-step volume remains avalanche-dominated. This favourable gyro-Bohm like scaling thus emerges for entirely different reasons that the oft-invoked argument of turbulent transport being adequately represented through a diffusive (random walk) process, which step size scales with ρ i . Bohm to gyro-Bohm transition in this interpretation is not due to fundamental changes in basic transport mechanisms. It naturally appears through combination of forcing (proximity to marginal stability and strength of the staircase organisation) and plasma volume (more shear layers in the radial as size increases);

Multi-scale interaction, forcing & boundary conditions: the resilience of staircase patterns in a background of avalanches, the processes through which layering occurs -including instability drive, which can occur at a variety of scales-with implications on scale selection (e.g. ρ e , ρ i , meso-scale) and the various types of long-time staircase evolution -including coarsening and migration-are as yet not well understood.

Similarly, the impact of forcing and boundary conditions is in its infancy. Since the systems of interest here are often out-of-equilibrium, these questions are central. Here again, cross-fertilisation will help clarify the relative role of stirring (e.g. the Phillips mechanism by which the stirring of a stratified flow leads to layering) and nonlinear development of instabilities (such as the Kelvin-Helmholtz instability with either vertical or horizontal shear or drift-wave interchange in plasmas) in the onset and sustainment of layered solutions. Furthermore, there is compelling evidence in plasma physics and astrophysical fluids that significant differences in transport (as a result of fundamental differences in layering) may be driven upon modification of either boundary geometry, proximity to nonlinear instability threshold or forcing -whilst either constraining a constant flux or imposing a constant temperature or density. While models are developing, they have yet to confront the fully multiscale character of the problem and, given realistic forcing and boundaries the truly global nature of the organisation at play -e.g. the role of the plasma edge in tokamaks, at the confluence of core and SOL influence yet influencing both as well is hard to overstate. In order to be able to incorporate the central question of large-scale and long-time influence of layering on transport, it is vital that the physics of layering and how it depends on the various parameters be better understood. Such questions are of great physical, mathematical and computational interest. It is of great research interest for the geophysical and environmental fluid dynamics community as well as in the plasma community and will help inform generic problems related to anti-diffusion, motivated by a considerable variety of physical and biological phenomena.

Tailor/influence staircase organisation to one's needs? As flow/staircase structures are elusive to diagnose, we might collectively be already factoring in the effect of such organisation in some of the experimental data, like M. Jourdain would in a famous play by Molière. A relevant question to ask is whether one could find actuators to weaken/reinforce locally the spontaneous organisation and thus open access to improved performance. Staircase step meandering and merging could be a way to obtain larger/less permeable barriers or even (?) transport bifurcations in addition to the more common possibilities that are based on local tailoring the q profile or magnetic shear. Modulation experiments (which could provide an interesting perspective on the old problem of the "cold pulse") and the interplay through momentum and vorticity transport with the staircase layers are likely the way forward to build understanding. This echoes the aforementioned problem with boundary conditions as well for staircase structures emerge in reduced models as a boundary value problem due to the growth, coarsening, and saturation of modulations.

last words for experimentalists: whilst indeed focusing on dynamics rather than previous staircase experiments which were concerned primarily with existence and characterisation, such experiments would be more challenging to theory and set higher benchmarks for validation. Briefly they should target the following issues: (i) do modulations amplify or destroy pre-existing staircases? (ii) what is the connection between staircasing and what is often broadly referred to as "nonlocality" phenomena? (iii) on what profiles do staircases form: density, temperature, current,. . . ? (iv) what is (are?) the critical parameter(s) for staircase formation ( Q turb Qneo ν ef f is a possibility though it is likely not the whole story)? (v) what is the meandering/destruction/reformation dynamics of staircases -and do they always reform in the wake of a destructive event? (vi) can an emergent dynamic scale be identified which would be key to staircase formation? (vii) can the merging/collapse of staircase steps into meso-or macro-scale barriers be observed and what of its implications for confinement?

Jesting interlude

"Rien n'est plus puissant qu'une idée dont l'heure est venue" Victor Hugo Personal note: this research is both an entry door to beautiful fundamental questions regarding complex systems and a possible interesting road ahead for Society. It manages a form of reconciliation between personal inclinations: developing fundamental knowledge, wishing to be a physicist whilst at the same time being a somewhat mindful citizen "at the World", as Merleau-Ponty would say. Unfortunately, on this latter point, beyond naive enthusiasm, one can always fear (amongst other things) that if it were to become practical one day, this publicly funded research --collective, trans-generational investments--will likely be sold out the next day and capital gains cornered. This, in the present author's view is quite unacceptable:

"On est gouvernés par des lascars qui fixent le prix de la betterave et qui ne sauraient pas faire pousser des radis"

Michel Audiard

Illustrated differently, crime of Hubris is not faced solely by mythical figures, scientists or engineers, as made abundantly clear by the 'interesting' (sic) viewpoint below, which has unfortunately grown almost as pervasive as it is lethal: "Macroeconomics was born as a distinct field in the 1940's, as a part of the intellectual response to the Great Depression. The term then referred to the body of knowledge and expertise that we hoped would prevent the recurrence of that economic disaster. My thesis in this lecture is that macroeconomics in this original sense has succeeded: Its central problem of depression prevention has been solved, for all practical purposes, and has in fact been solved for many decades" Robert E. Lucas Jr, Nobel Prize in Economics 1995 American Economic Review (2003) -5 years before the great 2008 crisis. . . As a physicist, one clearly wonders, provided our community does not organise and push back on coming political decisions whether it is worthwhile to be working to contributing a way out for a deleterious political system. Though I don't fully subscribe to all of its implications, Brassens (as often) has a point: "Un anarchiste est un homme qui traverse scrupuleusement entre les clous, parce qu'il a une sainte horreur de discuter avec les agents"

Georges Brassens

And so does Mark Twain, ironically remarking that: "Prosperity is the best protector of principle"

Mark Twain

Our political organisations (unfortunately) vastly focus on the first part and forget that prosperity comes about due to a predatory behaviour on (amongst other things) our natural resources. The latter part of Twain's remark is invoked when agreeable; its tartness (and pertinence) is often lost altogether. Yet, in the end, it is men and women that produce effects on the world and "A man cannot be comfortable without his own approval"

Mark Twain

Not much room is provided within our Institutions to practically enact these thoughts. The views and opinions expressed herein do not necessarily reflect those of the Institution which employs the present author.

Chapter 3

The plasma edge: a key region under conflicting influences Abstract: Magnetically-confined thermonuclear plasmas hold the promise of sustainable and safe fusion production on Earth. For over three decades, the observation of rapid core confinement improvement upon favourable modifications of edge operating conditions has been a nagging source of puzzlement. Transport properties of driftwave turbulence in the plasma edge and interaction of the confined plasma with its material boundaries have long been recognised as essential to the resolution of this conundrum. Basic understanding of these questions has however remained largely elusive. We show here, from the primitive kinetic equations, that a narrow, localised region at the plasma-material boundary interface has a global impact on transport. It is especially central to explaining transport properties of turbulence in the outer edge and responsible for spontaneous transport barrier onset at the plasma edge. These observations strongly suggest a paradigm shift where turbulence in the confined plasma (the "dog") is not only locally driven by the local gradients but 'nonlocally' controlled by fluxes of turbulence activity, primarily though not exclusively coming from the edge (the "tail"). We anticipate these results to impact our understanding of transport barrier formation and access to high confinement regimes, highly relevant to the safe operation of large burning plasma systems such as the Iter tokamak.

Why focus on the edge?

Magnetised fusion plasmas are oft separated, as sketched in Fig. 2.1, between a hot confined core region where fusion occurs, an unconfined peripheral boundary layer (the Scrape-Off Layer or SOL) -cold and comparatively rarefied, starting at the magnetic separatrix and magnetically connected to the material boundaries-and an in-between edge region, loosely defined, set between core and separatrix. Core and SOL have been extensively studied, mostly independently, the edge usually serving as fixed boundary condition for both, its dynamics difficult to apprehend. Strict decoupling however between all three regions is increasingly being questioned (e.g. [Dif-Pradalier et al. 2017a, Caschera et al. 2018, Dif-Pradalier et al. 2021]). Tokamak plasmas are indeed prone to self-organisation, a prominent example being the spontaneous transition between the low-confinement "L-mode" to high-confinement "H-mode" [Wagner et al. 1982]. The H-mode branch of operation is one of several improved confinement states that have been experimentally discovered, revitalising the fusion programme towards ITER. Significant performance enhancement has been made possible by exploiting the spontaneous L to H transition, induced by increased radial electric field shear E ′ r . Onset of differential rotation, which scales with E ′ r and steepening of the ion pressure profile in a localised region of the L-mode edge -the so-called "pedestal"-stabilises turbulence, reduces transport and initiates a self-reinforcing feedback [Biglari et al. 1990, Kim & Diamond 2003] which locks-in the bifurcated state. Comprehensive description of this dynamics from first principles is still lacking and certainly requires the intermediate step of apprehending turbulence dynamics in the L-mode edge, from where the bifurcation originates. More generally, realistic modelling of fusion-grade plasmas must address this delicate balance and tackle dynamics in the vicinity of the closed to open field line interface.

3.1.1 A thorny, multi-scale problem where forcing and boundary conditions matter

"Far better an approximate answer to the right question, which is often vague, than an exact answer to the wrong question, which can always be made precise"

John Tukey

A tradition of works has attempted to model the separatrix and edge regions. Owing to the strength of the guiding magnetic field, plasma turbulence is often stated as quasi two-dimensional. Unstable modes helically extend along field lines whilst being at leading order radially pinned to flux surfaces. On the basis of this property, separations of scales between fluctuations and mean background are commonly performed. Such approaches, referred to as "gradient-driven" [Sydora 1995, Kotschenreuther et al. 1995, Beer & Hammett 1996, Dimits et al. 1996, Candy & Waltz 2003, McMillan et al. 2008, Peeters et al. 2009, Görler et al. 2014] are computationally efficient and have extensively been used in the plasma core (this point is further discussed in section 2.3.2.1 and Appendix C.1). Their validity however faces challenges in the edge and wanes further outwards whilst nearing the separatrix and SOL. Turbulent scales at the plasma edge indeed become comparable to free energy gradient scales and oft-assumed separations between a slow, large-scale background and fast, small-scale fluctuations become inadequate, as already noted by Kadomtsev [START_REF] Kadomtsev | [END_REF].

In the peripheral plasma, gradients are steep, with intrinsic temporal variability. Precise experimental measurements are challenging. Fluxes may vary by factors in gradient-driven frameworks upon scanning imposed mean gradients within experimental error bars. This generates considerable difficulty to predict performance of magnetic confinement devices or to assess safe operation of Iter. At statistical equilibrium however, power that has been imparted to the plasma must come out. This driving power is known. This has led to proposing a paradigm shift in modelling where known, imposed fluxes drive the system and both mean gradients and fluctuations dynamically evolve in concert, act and back-react on one another. Such a framework (further discussed as well in sections 2.3.2.1 and Appendix C.1) is referred to as "flux-driven" [Carreras et al. 1996, Sarazin & Ghendrih 1998, Garbet & Waltz 1998, Beyer et al. 2000, Ku et al. 2009, Idomura et al. 2009, Sarazin et al. 2010, Korpilo et al. 2016]. It relaxes assumptions of scale separation and comes at the price of a significant (tenfold or more) increase in computational demands.

Tokamak plasmas are well-known to display transport "stiffness", the property that upon surpassing a critical gradient turbulence activity rapidly increases, leading to large fluxes which in turn rapidly deplete the driving gradients, lowering the turbulence intensity and resulting in the fact that actual gradients hover in the vicinity of their nonlinear 'nearmarginal' threshold values. This feedback acknowledges the experimental fact that sources are finite, as taken into account in flux-driven frameworks. Stiffness may be modified, eg. by safety factor tailoring or fast-ion enhanced electromagnetic stabilisation [Citrin et al. 2013]. Several problems are associated with this stiffness property as the success of models especially hinges on the reproduction of experimental nonlinear fluxes. Experimental uncertainties on mean gradients are commonly important, well above the percent. Gradient-driven Figure 3.1: Snapshots of the electrostatic potential in three different configurations, detailed on p.70, section 3.2.1.2, at statistical equilibrium. Case-1 is the reference flux-driven configuration and features, as in experiments, a transition from closed to open field lines through introduction of a poloidally-localized toroidal limiter (separatrix is at normalized radius ρ = r/a = 1) and a wall (at ρ = 1.3) within the computational domain. Case-2 differs from Case-1 through its external sink alone: uniformly distributed in the poloidal (along θ) and of gradually increasing strength past ρ ≥ 1 (see Methods). Cases 1 and 2 are initialized identically; subsequent discrepancies in their temporal evolution is thus direct incidence of the outer boundary condition alone. Case-3 is the gradient-driven twin of Case-2: it tests for the influence of turbulence spreading on the global self-organized state. Background mean gradients for Case-3 are the converged profiles of flux-driven Case-2, other parameters being equal between both computations. Case-3 hovers throughout nonlinear evolution about the statistical state of flux-driven Case-2, effectively prescribing a scale separation between background and fluctuations. For proper comparisons, the same radial-poloidal areas are magnified for each Case: top of the machine [column (a)]; low-field side midplane (b) and bottom regions (c). Reference (circled) positions are systematically displayed in all subplots as well as in Figs. 3 and4 and supplementary Figs. 3 and4 at various radial (ρ = 0.9, 0.96 and 1.02) and poloidal (θ 1 through θ 5 ) locations. The separatrix in Case-1 is clearly apparent; for comparison, it is drawn at the same locations for Cases 2 and 3 as dotted grey lines in subplots (b). The solid blue line at r/a ∼ 0.84 represents the approximate transition (see Appendix F.0.2) between (linearly) convectively unstable core and convectively stable edge. models have been reported to display from factors to orders of magnitude variations in computed fluxes upon a few percent variation of their input gradients. This sensitivity poses considerable difficulty in assessing the performance of gradient-driven predictions in a stiff context. Furthermore, the edge of L-mode plasmas has been experimentally reported as less stiff than the plasma core [Sauter et al. 2014]. The mechanism whereby this "destiffening" occurs is yet unclear and it further contributes to questioning applicability of stiff gradient-driven frameworks to the plasma edge.

In fully bifurcated H modes, the large gradients have built up over a limited region of space -the so-called edge "pedestal". With such large gradients 1/L p , electromagnetic modes become more prominent as the effective beta β ef f = β(qR/L p )2 for instability increases. Similarly, validity of gyrokinetics [Brizard & Hahm 2007] which requires ρ/L p < 1 is sometimes quoted as less straightforward. No understanding of H-mode however may be complete without comprehensive understanding of the degraded confinement in the first place. The present study is thus concerned with elucidating onset of edge dynamics in L-mode where β ef f remains moderate and where millimetric Larmor radii and centimetric gradient lengths is reasonable ground for gyrokinetics.

Furthermore, transport processes in the edge qualitatively switch from being dominantly of conductive nature in the core (the confinement time of particles is long with respect to the confinement time of energy) to being dominantly convective in the vicinity of the SOL and of the material boundaries. This change in the nature of transport occurs as magnetic field lines open up, connect to material boundaries filaments detach from the plasma at the crossing of the magnetic separatrix and particle recycling increases. Assuming low rotation and subdominant turbulence-generated flows -a debatable assumption in some experimental situations [START_REF] Crombé | [END_REF], Bell et al. 1998] and theoretically at staircase step locations in near marginal regimes [START_REF] Dif-Pradalier | [END_REF], Dif-Pradalier et al. 2017b]-the radial force balance E r -v φ B θ + v θ B φ = ∇p/(ne) reduces to E r ∼ ∇p/(ne) < 0 at leading order in the core. Conversely, in the SOL the radial electric field is dominated by the narrow boundary layer, the so-called "sheath", which develops in the parallel at the interface between plasma and material boundary, typically over a Debye length λ D = (ε 0 T /ne 2 ) 1/2 . An essential outcome of this physics, further discussed below in section 3.2.1, consists in the establishment of a sharp potential drop across the sheath which confines the plasma bulk electrons. As a consequence, ions are accelerated from the plasma towards the wall, therefore driving at leader order a sonic ion fluid velocity at the sheath entrance [Stangeby 2000].

∆ϕ b.c. = ϕ w -ϕ se = - T e 2e ln 2π m e m i 1 + T i T e ⇒ E r ∼ -∇T e > 0 (3.1)
Other factors further complicate the physics at the plasma edge, which we do not consider here. For instance, field geometry: recent advances in positive versus negative triangularity [START_REF] Camenen | [END_REF], Camenen et al. 2007, Austin et al. 2019, Fontana et al. 2019, Marinoni et al. 2019] have stressed the potential importance of this point; there is also plasma elongation, limited versus diverted geometry, etc. Interaction with neutrals is also a central topic as neutral penetration depth and spatial distribution constrain particle sources in the edge.

The plasma edge, central to understand improved confinement and safe operation is thus at the confluence of complex and often competing influences. As no understanding of H-mode or of other improved confinement regimes may be complete without comprehensive understanding of the degraded confinement in the first place, we now concern ourselves with the depiction of edge dynamics in L-mode from, flux-driven gyrokinetics.

Some outstanding problems

"Je suis ancien combattant, militant socialiste et bistrot. C'est dire si, dans ma vie, j'en ai entendu, des conneries"

Michel Audiard

As stressed above, the L-mode edge is challenging from a modelling standpoint. Relative turbulence fluctuations levels δn/n, n being the plasma density are experimentally observed to increase from core to separatrix [Liewer 1985]. At the same time, large fractions of the outer edge are either linearly stable or marginally unstable to drift-wave interchange turbulence, questioning the origin of observed fluctuations there. The origin of turbulence activity in the edge has prompted diverse speculation, connecting to a tradition of works that have long recognised that convective motions where no impenetrable boundaries are present can expand beyond the region of convective drive. Interfacial contam-ination by turbulence is referred to as nibbling [Corrsin & Kistler 1955] or engulfment [Townsend 1976] in jet interfaces, overshooting or penetration [Brummell et al. 2002] in geophysical and astrophysical fluids or spreading [START_REF] Garbet | Radial propagation of turbulence in tokamaks[END_REF], Hahm et al. 2004] in magnetised plasmas. Subcritical transitions in parallel flows have been evidenced, belonging to the class of directed percolation [Pomeau 1986]. Similar mechanisms in plasmas are appealing [Heinonen & Diamond 2019] yet precise mechanisms whereby the L-mode edge transitions to and sustains a turbulent state are unclear. Hence its sometimes-used name of "No Man's Land" (NMsL) bearing witness to the difficult understanding of its dynamics. The incentives to clarify this problem are high for turbulence activity in the edge controls stress and electric field dynamics, affecting access to transport bifurcated states and deposition patterns on plasma facing components, highly relevant to the safe operation of large burning plasma systems.

We discuss a generic situation, based on experimental Tore Supra parameters from shot #45,511 where the following conundrum is found: experimentally, the edge is measured to be turbulent, with fluctuations increasing with proximity to the separatrix. In contrast, local analysis of the profiles predicts convective stability in the edge and increased stability with proximity to the separatrix -see appendix F for further details. In this context, linearly stable or marginally unstable "No Man's Land" has long since been an area of active speculation. This problem is especially thorny within a stiff modelling framework as modest variations in the driving gradients can predict from linear stability to fully active turbulence. Several attempts have repeatedly failed at reproducing the experimentally-observed turbulence behaviour in NMsL, under-predicting by factors the level of fluctuations and fluxes there [Scott et al. 1991, Holland et al. 2011]. Increased resolution [Waltz 2017], electromagnetic effects or kinetic effects beyond gyrokinetic limitations (unstable low-frequency drift modes couple to stable high frequency ion cyclotron modes in the so-called "cyclokinetic" approach [Deng et al. 2016]) have not been found to cure the shortage of turbulent transport in the plasma edge. This underprediction is commonly referred to as the "transport shortfall" conundrum. An interesting approach to the shortfall problem is described in Ref. [Görler et al. 2014] and Fig. 3.2. There, reconciliation of experimental and computed transport levels is attempted, up to r/a ≲ 0.85 within the gradient-driven framework of Gene. A shortfall is found at nominal parameters: computed transport levels decrease as r/a increases; ion heat fluxes are overestimated by about ∼ 60% at r/a ∼ 0.5 and underestimated by a factor ∼ 2 between 0.6 ≤ r/a ≤ 0.85. The ion temperature gradient is then scanned at all radial locations where fluxes are evaluated until computed transport levels match measured transport levels. The figure of merit becomes: "is this required gradient variation (-10%, +10%, +22%, -2%) within experimental error bars"? On the basis of a positive answer, it is sometimes quoted that the shortfall conundrum is no longer a problem.

One could argue differently on the basis that: (i) corrections (-10%, +10%, +22%, -2%) are not known a priori. They are also (ii) not systematic, non smooth and sometimes of significant amplitude (+22%), all of which makes predictions of transport levels impractical without first knowing the experimental reference values. Furthermore, (iii) this flux matching procedure is not reported beyond r/a ∼ 0.85 when increased proximity with the separatrix is expected to be increasingly difficult from a modelling standpoint, e.g. due to interplay with open field lines or with the boundaries. One could ask, not for an iterative a posteriori matching procedure but rather for a more comprehensive and predictive understanding of transport levels in the edge.

We are thus concerned for the remainder of this chapter with the following questions:

• what are the mechanisms whereby turbulence activity may exist in the edge?

• how does this turbulence redistribute and possibly interplay with the core and SOL?

• what are the implications for flow generation and stress patterns in the vicinity of the magnetic separatrix?

which translate, on more practical terms, to:

1. can we provide new perspectives on the decades-long conundrum of predicted turbulence activity "shortfall" at the plasma edge?

2. what link does the shortfall problem entertain with edge transport barrier onset?

3. what basic mechanisms underpin spontaneous barrier onset at the separatrix?

SOL-edge-core interplay through improved boundary conditions

Given the complexity of the L-mode edge and of understanding early stages of edge pedestal build-up, there is considerable merit to search for a "bare-bones first-principles" approach as the model appropriate to provide resolution of the edge transport shortfall problem.

The work reported here is based on electrostatic gyrokinetics. It resolves the whole plasma volume, from core to wall. It includes transition from closed to open field lines, in fluxdriven regimes, and thus takes a alternative approach to the common flux-tube, local, gradient-driven approximations. Massive increase in computational power in recent years has been instrumental to permit the present investigation. Our study connects to a larger tradition of works that have shown pervasive though usually highly non-trivial impact of boundary conditions and of forcing on nonlinear evolution of complex systems. Recent examples in fluids have experimentally demonstrated that conjugate forcings (either energyconstrained or temperature-constrained) produce non-unique statistical steady states as well as different bifurcations between regimes -eg. Ref. [START_REF] Saint-Michel | [END_REF]]. We highlight similar conclusions here for magnetic fusion, with substantial implications for modelling.

Penalisation to model core, edge and a simplified SOL

Transition from closed to open field lines comes from the introduction of material boundaries within the computational domain. Many techniques have been proposed to treat numerically the interface between plasmas or fluids and solid boundaries. Discussion of their relative merits is beyond the scope of the present manuscript. We focus here on the so-called "penalisation" technique, which models the interface between plasma and wall with an immersed boundary technique. Such approaches were originally developed in computational fluid dynamics [Peskin 2002, Mittal & Iaccarino 2005] and have since been introduced in plasmas [Isoardi et al. 2010, Paredes et al. 2014, Caschera et al. 2018]. Immersed boundaries methods allow to simulate viscous flows within physical confines on numerical grids that do not conform to the shape of these confines. The technique of penalisation uses masks to define both (i) the precise geometry of the immersed material boundaries and (ii) to modify, in volume, the governing equations within application of these masks -hence the name "penalisation". With this technique, boundary conditions in an usual sense still exist but are demoted to the frontier of the computational domain, far from the physically-motivated boundaries that are within the bulk of the numerical grid. Through application of the masks, physical boundaries are thus modelled as effective sources or sinks. This has several advantages: (i) complex 3-dimensional boundaries are easily treated;

(ii) boundaries can vary in time without needing to regenerate the grids;

(iii) grids can be structured or unstructured. They do not require to conform to the surface of the material body;

(iv) both plasma and wall are treated on an equal footing -the material boundary is effectively computed as if it were a plasma with altered dynamics: adjustment of the transition layer between plasma and wall is dynamic and self-consistent. The transition is progressive in a controlled (and somewhat physically-motivated) fashion, thus often less prone to numerical instabilities than more regular boundary conditions with abrupt transitions at the edge of the numerical grid;

(v) the shape of the masks (i.e. abruptness of the transition), strength with which one alters the system's dynamics when transitioning from plasma to wall (i.e. the strength of the penalisation) and the penalised solution towards which the system is progressively driven (e.g. the target wall distribution function) are controllable parameters. These three parameters actually represent the control knobs that one needs to carefully adjust to mimic interplay with a complex boundary.

Penalisation is here applied to the usual [Grandgirard et al. 2016] Gysela equations to mimic the presence of a toroidal limiter and a wall. Complete details are given in Appendix.A. Both the gyrokinetic equation for the ion guiding-center distribution function Fs and quasi-neutrality equation are modified through introduction of a series of masks M mat (r, θ), M SOL (r) and M wall (r), combinations of hyperbolic tangents of value in range [0, 1]. They are adjustable in location, stiffness and shape and are illustrated in Fig. A.1. The computational domain extends from inner core (r/a = 0) to the material boundaries (r/a = 1.3).

Penalisation of the quasi-neutrality equation

Unless stated otherwise, a Boltzmann response is now assumed for the electrons in the core, edge and SOL, both for the sake of simplicity and of numerical cost. Generalisation to a kinetic response of trapped electrons in the confined core is straightforward, provided passing electrons, electrons in the SOL and electrons intercepting a material boundary remain treated adiabatically. This is detailed in Appendix E. The case of a fully kinetic response for SOL particles or particles intercepting material boundaries is non trivial and under current active investigation [Munschy et al. 2023b, Munschy et al. 2023a]. This point is briefly touched upon in section 3.4.

At the interface between plasma and wall, a boundary layer -the "sheath" develops. Sheath physics has attracted considerable attention for about a century in various areas of physics, see e.g. Ref. [Stangeby 2000]. In the following, we will consider the simplified problem where magnetic field lines intercept material boundaries at normal angles. Nonnormal angles lead to additional complexity, onset of a magnetised pre-sheath -see e.g. recent works such as Ref. [Geraldini et al. 2019] and will not be touched upon here. The difference of mobility between heavy ions and light electrons (the typical small parameter pervasive in these regions is m i /m e ∼ √ 2 × 1836 ∼ 60.6 for Deuterium) is what fuels sheath dynamics. As fast electrons are lost to the wall, a steady-state sharp potential drop establishes in the parallel, at the interface between plasma and material boundary, typically over a Debye length λ D = (ε 0 T /ne 2 ) 1/2 . This sharp potential drop across the sheath repels -i.e. confines bulk (slow, thermal) electrons and accelerates ions from the plasma towards the wall, therefore driving at leader order a sonic ion fluid velocity at the sheath entrance.

The sheath is thus a region of positive space charge which shields the quasi-neutral plasma and allows for dynamic charge equilibration between plasma and wall.

The radial electric field throughout the SOL is dominated by sheath physics. The typical scale at which it establishes is the Debye length, which is well below the ion Larmor radius and therefore at subgrid scales for ion-scale gyrokinetic computations. Part of the problem that we are concerned with here consists in incorporating essential parts of this physics whilst not resolving microscopic sheath dynamics. In particular, one wishes to mimic the behavioural change of the electric field from confined plasma to unconfined SOL described in section 3.1.1. A simple way to do this is to specifically alter the quasi-neutrality equation Eq.(??) in the SOL and there force relaxation of the electric potential ϕ towards its expected presheath condition ΛT e /e.

In order to do that, different regions of the plasma may be selected or dismissed using mask functions which select specific regions of the numerical domain. Various mask functions are displayed in Fig. A.1. They take values between 0 and 1 with control over the abruptness of transition (typically through hyperbolic tangents). This technique is versatile: many different geometries of the boundary can be implemented, possibly dynamically, by simply choosing or altering shapes of various mask, while the underlying magnetic configuration remains unmodified. The geometry of the immersed boundary in the poloidal plane is essentially defined by mask M mat (r, θ), which mimics an axisymmetric Tore Supra-like limiter and wall. The last closed flux surface (LCFS), or magnetic separatrix (dashed line in Fig. A.1) is tangent to the limiter. Combination M SOL (r) -M mat (r, θ) selects the SOL plasma where ϕ relaxes to ΛT e /e. Similarly, combination M mat (r, θ) -M wall (r) selects the limiter, which can be either biased ϕ bias ̸ = 0 or grounded ϕ bias = 0. The latter is chosen in the present manuscript. As a side technical remark, it is best that masks be added rather than multiplied so that the amplitude, transition abruptness, etc. of their combination can be readily inferred by the amplitude, transition abruptness, etc. of the individual mask functions. At last, T b.c. e is the cold electron temperature within limiter and wall, chosen as the minimum T e value within the computational domain, Λ ∼ log( m i /m e ) is given by Eq.(3.1) -see Appendix.A for complete details-and coefficient λ (set to unity in the present study) may be used to alter the inertia of the zonal potential.

Penalisation of the gyrokinetic equation

Flux-or gradient-driven dynamics may be considered, switching on and off terms in Eqs.(A.1) and (A. 2) amounts to a localised heat and momentum sink, its shape set by mask M mat (r, θ) which mimics arbitrarily complex limiter and wall geometries. In order to mimic transitions from plasma to material bourdaries as sharp as possible, one would wish for ν → ∞ in Eq.(A.2). In practice, this is source of potential instabilities so infinite penalisation is preferred [Caschera et al. 2018] in lieu of the BGK term -νM mat (r, θ)[ Fs -G cold ] where ν would be chosen as large. This choice is detailed in Appendix D. Penalisation acts so as to relax Fs to a target cold Maxwellian distribution function

G = G cold = n w (2πT w ) -3/2 exp[-(v 2 G∥ + µB)/2T w ],
characterised by low wall thermal energy T w and target density n w . The former is constrained by velocity-space resolution; we typically choose it an order of magnitude lower that temperature at mid radius whilst the target density n w is chosen so as to ensure particle conservation.

Case-2 is flux-driven, forced with the same heat source S in Eq.(A.1) as in Case-1. Here, D = ν L = γ GD = 0 and M mat (r, θ) = M wall (r) = 0. The heat sink is provided by the remaining term in Eq.(A.2). Its strength is uniformly applied in the poloidal and progressively increases in the radial past r/a ≥ 1, until all fluctuations are damped at r = r max . In Case-2, field lines are all closed. The presence of a SOL-like solution past r/a ≥ 1 is felt through the choice of the target distribution G, chosen as the equivalent Maxwellian computed from exponentially-decaying density and temperature profiles, smoothly continued at r/a = 1 from the initial profiles. The system adapts in the vicinity of r/a ∼ 1, realising a compromise between unaltered flux-driven dynamics in the core and edge and an increasingly constrained SOL-like solution past r/a ≥ 1.

Case-3 is also shown in Fig. 3.1 and mimics gradient-driven-like dynamics when

S = D = M mat (r, θ) = M wall (r) = M SOL (r) = ν S = ν L = 0 and γ GD ̸ = 0.
In that case, the target distribution function F F-D is chosen as the statistical distribution at equilibrium from flux-driven Case-2. The relaxation rate γ GD = 5.43 10 -5 ∼ γ lin /10 is constant and chosen an order of magnitude smaller than the local linear turbulence growth rate γ lin at r/a = 0.7, thereby realising a compromise between sustaining the drive (the background gradients) and affecting as little as possible the linear growth physics. This last term in Eq.(A.2) is also built such as to prevent nonlinearly overdamping zonal modes [McMillan et al. 2008]. Gradient-driven like dynamics resulting from the use of this operator differ from their flux-driven counterpart. This is discussed below for the edge and illustrated for the core plasma in Ref. [Dif-Pradalier et al. 2017b].

Case-3 is the gradient-driven twin of Case-2. Let us note on a minor difference however in the resolution of the quasineutrality equation: Case-2 is usually run with non vanishing mask M SOL (r) in Eq.(??), whilst this is not the case for Case-3. The evolution of the SOL electric potential is thus slightly different; these are modest differences in view of the imposed relaxation in Case-3 towards F F-D from Case-2, everywhere, including in the SOL, which dominates the dynamics.

New light on old conundrums

"If the fool would persist in his folly he would become wise"

William Blake

We now focus on a comprehensive discussion of the mechanisms whereby turbulence organises in the plasma edge when core, edge and SOL dynamically interplay. Given the complexity of understanding early stages of edge pedestal build-up, there is therefore considerable merit to search for a "bare-bones first-principles" approach as the model appropriate to provide resolution of the edge transport shortfall problem. The present work is based on flux-driven electrostatic gyrokinetics. It takes a complementary approach to the common flux-tube, local, gradient-driven approximations. Massive increase in computational power in recent years has been instrumental to permit the present investigation. This study connects to a larger tradition of works that have shown pervasive though usually highly non-trivial impact of boundary conditions and of forcing on nonlinear evolution of complex systems. Recent examples in fluids have experimentally demonstrated that conjugate forcings (either energy-constrained or temperature-constrained) produce non-unique statistical steady states as well as different bifurcations between regimes (eg. Ref. [START_REF] Saint-Michel | [END_REF]). We highlight similar conclusions here for magnetic fusion, with substantial implications for modelling. Notably, we find that the elusive questions of stability of NMsL and the sought-after onset of an edge transport barrier are intimately connected. Key results are:

(i) boundary-mediated interactions trigger spontaneous symmetry breaking and localised destabilisation of the outer plasma edge. Constant and localised turbulence intensity production follows, (ii) which spreads, contaminating stable NMsL on global scales and providing a possible resolution of the shortfall problem. Concomitant to turbulence production, (iii) a persistent transport barrier emerges, at nominal experimental parameters. Diamagnetic currents are found to be surprisingly influential to edge transport barrier onset, advecting vorticity and electric field inhomogeneities.

(iv) At last, removal of the plasma-boundary interface or postulation of scale separations hinders turbulence spreading, resulting in a laminar edge (a shortfall), which prevents access to edge transport barrier onset.

Free energy injection at the magnetic separatrix

Interestingly, a marked difference in edge turbulence (Fig. 3.1) is found upon modification of the outer plasma boundary conditions from limiter (poloidally asymmetric with closed/open field line transition) to radially progressive and poloidally symmetric (closed field lines only).

Case-1 is our reference computation, which mimics tokamak parameters of the Tore Supra shot #45511, using the GYSELA framework [Gra16]. The coupled gyrokinetic and Poisson equations are time stepped in five dimensions, from core (ρ = r/a = 0) to wall (ρ = 1.3).

Cases 2 and 3 are companion computations which only differ from Case-1 either by forcing (flux-driven versus gradient-driven) or by boundary conditions (limited versus poloidally symmetric).

Comparing Case-1 with Cases-2 and 3 under similar plasma parameters probes the role of poloidal asymmetry and of closed/open field line transition in edge barrier onset. Gradient inhomogeneities in the poloidal plane are pervasive in experiments, though often elusive to diagnose and may stem from the presence of a limiter, of a radiating X-point or of high-Z impurities. Their determination is a central matter for transport [Angioni & Helander 2014, Donnel et al. 2018]. Interaction with the limiter in Case-1 triggers a spontaneous symmetry breaking: it onsets and sustains an over-dense shelf in the poloidal vicinity of the limiter and (e)], a poloidal distribution of temperature gradients [Fig. 2-(d) and (f)] and also induces a radial electric field well (Fig. 3.3-b), a signature of improved confinement. All three features are absent in poloidally homogeneous Cases 2 and 3.

Furthermore, precise comparison between flux-driven and gradient-driven forcing indeed addresses the role of turbulence spreading for scale separations in gradient-driven frameworks hinders the back-reaction of fluctuations on the driving mean gradients necessary for the self-advection of turbulence. Quantification of the actual importance of spreading from the primitive equations is largely uncharted. Comparison of Case-3 to Cases 1 and 2 therefore provides a comprehensive view of the actual role of turbulence spreading [START_REF] Garbet | Radial propagation of turbulence in tokamaks[END_REF], Hahm et al. 2004, Heinonen & Diamond 2019] in edge self-organisation.

Linear stability of the L-mode edge

To clarify the nature of the edge free energy in all three cases, extensive linear stability analysis of the three Cases has been performed at various radial-poloidal (ρ j , θ k ) edge Figure 3.3: Spontaneous symmetry breaking from plasma-boundary interplay leads to new free energy sources. Density (c) and temperature (d) profiles are colourcoded with poloidal proximity to the limiter (a). Hot plasma particles rapidly stream towards the limiter along field lines as well as through the action of vertical B × ∇B magnetic drift (here pointing towards the limiter). This dynamics leads in Case-1 to the onset and sustainment of an over-dense poloidal density shelf (Fig. 2-c ande) near the limiter as well as localised temperature gradients (Fig. 2-d andf) and a radial electric field well (Fig. 2-b, averaged poloidally), absent in poloidally homogeneous Cases 2 and 3. This poloidal anisotropy is especially significant within 5% of the separatrix, as seen through radial averaging ⟨•⟩ ∆ρ of equilibrium mean gradients (Fig. 2-e andf). Gradient anisotropy decreases in magnitude in the near SOL and poloidally farther from the limiter. Large and anisotropic equilibrium mean gradients in the limited configuration result in locallyenhanced free energy sources for the turbulence in the outermost 5% layer of the confined edge. Inversely, Cases-2 and 3 display uniform and moderate mean gradients, shown in thick yellow and dotted blue. locations (circles in Fig. 3.1), using the initial value framework of the Gyrokinetic Workshop (Gkw) code [Peeters et al. 2009]. It leads to the following picture (see Appendix F.0.2): (i) edge profiles in Cases-2 and 3 are linearly stable to drift-wave, interchange (DWI) and parallel shear flow instabilities past ρ ≥ 0.9 whilst (ii) a localised region of instability, dominantly of interchange character appears in Case-1 immediately inside the separatrix due to the presence of the limiter, on both sides of it. Stability of this large portion ρ ≥ 0.9 of the edge echoes the weak estimated convective drive in "No Man's Land" (NMsL), bearing witness to its puzzling stability.

Electron-scale instabilities or electromagnetic modes, absent in the present modelling may also be active in actual tokamaks and contribute to edge transport. The present work focuses on robust shear flow and ion-scale modes only and thus has the powerful appeal of relative simplicity, whilst being at the fringes of current-day flux-driven supercomputing capabilities. It allows to define a minimally-relevant set of ingredients that produces transport in the edge consistent with experimental levels and naturally triggers and sustains an edge transport barrier. As such it arguably goes some way towards disentangling the relevant hierarchy of basic mechanisms at play: plasma-boundary interaction builds the local dynamics of the "tail" which helps confinement through edge transport barrier onset and wags the "dog" through scale interplay.

Importantly, mean profiles at the low field side midplane (θ = 0) where experimental measurements are usually performed are very similar between all three cases. On the basis of local analysis of available free energy, undistinguishable nonlinear evolution would be predicted there between all three cases. The limiter modifies this picture. As a cold spot, it tends to create radial-poloidal pressure anisotropies in its vicinity, in a manner akin to thermodiffusion described by off-diagonal (density and momentum pinch) terms in the transport matrix. Such a mechanism provides robust free energy injection nearby a localized heat or momentum sink from combined onset of a radial and poloidal pressure gradient.

The main findings are as follow: plasma-boundary interaction provides a pathway to a novel source of free energy that locally builds ion-scale, electrostatic turbulence in the vicinity of the cold sink. Then turbulence self-advection ("spreading"), which results from the possibility of flux-driven micro-to meso-scale interplay, redistributes fluctuations globally despite underlying convective stability. We note that this convective stability is robust and independent of the details of the electron response, Boltzmann or not (the edge r/a>0.84 is linearly stable, with or without this approximation, as shown by GKW). This is a possible indication that the nature of the instability is probably not as critical to edge dynamics as the combination of plasma-solid interactions and turbulence spreading. Lastly, we find that curing the edge shortfall problem through combination of the above mechanisms naturally opens the route for pedestal build-up, thus highlighting intimate connection between both problems. Vorticity advection appears as the important causal link to reconcile the somewhat counter-intuitive idea that additional turbulence in the edge may lead to pedestal build-up. Self-advection (spreading) of patches of turbulence intensity I may be quantitatively followed through wave-energy budget [Mattor & Diamond 1994]. The relevant conserved quantity is the negative of the entropy density nI, which involves the ambient density profile n and naturally leads to amplification of fluctuations I, as n decreases in the plasma edge. Wave-energy density conservation indeed implies that the fewer the oscillators, the larger the oscillations. The related conservation equation [START_REF] Gürcan | Transport of radial heat flux and second sound in fusion plasmas[END_REF]:

Turbulence spreading: instrumental to edge turbulence

D Dt (nI) + ∇ • Γ s = S (3.2) Γ s (r, θ, t) = 1 2 d 3 v (v E • ∇r) (δf ) 2 f M φ (3.3)
features the energy flux Γ s as kinetic proxy for spatial turbulence spreading. Here D Dt denotes the turbulent convective derivative, S free energy injection and dissipation mechanisms, (v E • ∇r) is the radial E × B velocity, (δf ) 2 = ( Fs -F M ) 2 the departure of the full ion distribution function Fs from an ensemble averaged Maxwellian F M reconstructed from evolving local density and temperature profiles and ⟨•⟩ φ denotes toroidal averaging.

The "beach effect" model [Mattor & Diamond 1994] does not satisfactorily account for our observations: core waves propagating on large radial distances and amplifying via convective instability in the edge is not the main mechanism found here. The maximum of relative fluctuations dn/n [Fig. 3.5-(c)] is also around the top of the high gradient region. This model however provides an interesting mental framework as it emphasises the role of poloidal asymmetry of propagating waves in order to obtain in toroidal geometry a non-zero radial energy flux, i.e. spreading, which we now show to be an essential ingredient.

Self-advection of turbulent patches is best apprehended as ∆S increments. Fig. 3.4 displays times series of poloidal cross-sections of spreading increments ∆S = Γ S (r, θ, t) -Γ s (r, θ, t ref ) for Case-1. Three phases appear: a NMsL devoid of fluctuations (Phase I, subplots (a) through (d)) is clearly visible in the early stages of Case-1 for all poloidal angles and for 0.85 ≤ r/a ≤ 1, echoing aforementioned linear stability of underlying profiles. Strong and persistent inward advection of turbulence intensity originates at the edge-limiter boundary (around θ = -π/2), propagates (white arrows) radially inwards and poloidally anticlockwise all the way to the top (θ = +π/2), in about 0.1ms [16, 700 Ω -1 ci ]. This turbulence intensity front propagates inwards until about r/a = 0.82, which amounts to a radial penetration depth of about 60 local Larmor radii, i.e. 10 to 12 local turbulence correlation lengths ℓ c .

From there [Phase II, subplots (e)-(f)], outwards radial spreading (black arrows) accompanied by clockwise poloidal motion from the plasma high field side and top regions to the low field side midplane (θ = 0) fills-in NMsL with turbulence in about 0.5ms. This to and fro redistribution of turbulence intensity bridges the region of free energy injection near the limiter with the upstream confined core. As both regions connect, core turbulence spills over, further enhancing edge fluctuation levels. This complex radial-poloidal dynamics finds an echo in the synthetic reconstructions [Fig. 3.4-(h) and (i)] of the radial profiles of turbulent fluctuations δn/n around θ = 0, as would be measured in Tore Supra. These synthetic profiles are plotted against typical actual measurements using fast-swept reflectometry (shaded grey) [Clairet et al. 2010]. A dynamic equilibrium (Phase III, subplot (g)), characterised by quasi-periodic relaxations of the edge turbulence is reached at later times (from 1.55ms onward). Bursts of clockwise rotating outgoing patches of turbulence (black arrows in Fig. 3.4-(g) and (i)) equilibrate incoming anticlockwise limiter-borne fluxes of turbulence intensity.

Comparatively, Case-2 [Fig. 3.5-(a)] and ] equilibrate on faster timescales, with spreading patterns dominantly outwards and reaching equilibrium in about 0.6ms. A clear shortfall is observed , top] past r/a ≥ 0.90 for Case-2 and past r/a ≥ 0.85 for Case-3, again expressing linear stability of underlying gradients, uncom- pensated by fluxes of turbulence activity from distant active (core or boundary) regions. Spreading in gradient-driven Case-3 is modest, at most about a few local turbulence correlation lengths, as in earlier works [Hahm et al. 2004], and insufficient to explain the required fluctuation levels in the edge.

Normalized contributions ∆ sp [j] = δn/n[j] -δn/n[j + 1] further quantify the weight of spreading , bottom] to overall statistical equilibrium, j denoting the Case index. Quantity ∆ sp [1] assesses outside→in spreading. Clearly, near-separatrix→core contamination of limiter-borne turbulence activity accounts almost in full for the fluctuation dynamics in the outer radial 10% of the confined plasma. Similarly, ∆ sp [2] quantifies the weight of inside→out spreading of core turbulence intensity towards the edge, amplified by the beach effect. Mean gradients in Cases 2 and 3 being the same by construction, the reversing sign of ∆ sp [2] illustrates the natural tendency of the flux-driven system with respect to the gradient-driven framework to radially redistribute patches of turbulence intensity, both outward (0.8 ≤ r/a ≤ 0.92) and inward (r/a ≤ 0.4). Extra spreading in flux-driven approaches, illustrated in subplot (a) by the black contour about θ 1 and lack thereof in Case-3 [same contour, subplot (b)] explains the gradient-driven underprediction of fluctuation levels in the deep core r/a ≤ 0.45 and outer edge regions 0.85 ≤ r/a ≤ 0.90 as well as its overprediction in the intermediate linearly-unstable region 0.55 ≤ r/a ≤ 0.75.

From turbulence in the edge to edge transport barrier onset

"Any fool can know. The point is to understand"

Albert Einstein

We have established that the combination of free energy injection at the plasmaboundary interface and of turbulence spreading provides a robust pathway to a turbulent edge. The last part of the triptych is to substantiate how turbulent eddies organise so as to unlock access to improved confinement. A fecund way to address this problem is to examine causal relationships and discuss the chain of events that presides over edge transport barrier onset. To this end, we apply a causality-seeking algorithm, the Transfer Entropy method to the gyrokinetic vorticity balance, which elucidates the dominant terms that may lead to flow onset. The "Transfer Entropy" (T E) method is a simple nonlinear extension of the Granger causality [Granger 1980], introduced by Schreiber [Schreiber 2000] and investigated in magnetised plasmas by Van Milligen et al. [van Milligen et al. 2014]. The idea behind T E is simple: let's consider a time series (x i ) of realisations of observable X, with 0 ≤ i ≤ n. If one can better predict its next realisation x n+1 using additional data from another time series (y j ) of observable Y with 0 ≤ j ≤ n, then "Y transfers information (i.e. causes) X", or more accurately as "Y forecasts X", which constitutes the definition of causality here. This idea is quantified measuring deviation of transition probabilities from independence, i.e. from a stationary Markov process. In its simplest expression, if processes X and Y are independent, then the following generalised Markov property holds for all 0

≤ k ≤ n: p(x n+1 |x n-k , y n-k ) = p(x n+1 |x n-k ).
The standard notation for conditional probabilities is used here: p(a|b) is the probability of a knowing b. If now processes X and Y are not independent, the ratio of these two transition probabilities provides a measure of how much information they may share. In other words, how much knowing values within Y in addition to past values in X may help to better evaluate next-step x n+1 . This idea leads to the following definition of the Transfer Entropy (T E) from process Y to process X:

TE Y →X (k) = p(x n+1 , x n-k , y n-k ) log p(x n+1 |x n-k , y n-k ) p(x n+1 |x n-k ) (3.4)
where k is thus a time lag and represents the k-past of times series X and Y . The summation process is detailed below, in Eq.(3.6). T E can equivalently be recast as a conditional mutual information and represents the additional amount of information that must be added to adequately represent the studied process p(x n+1 |x n-k , y n-k ) with respect to its reference Markov process p(x n+1 |x n-k ). In the absence of information flow from Y to X, the logarithm vanishes as state Y has no influence on the transition probabilities of X. It also follows that T E is directional, i.e. T E Y →X ̸ = T E X→Y , effectively allowing to infer causality between processes X and Y . T E displays interesting properties: it is independent of the relative magnitudes of signals X and Y ; it may apply to either linear and nonlinear regimes; it is easy to evaluate directly in real space rather than in Fourier space and is typically less demanding in terms of statistics than bispectral techniques. Practically, T E is evaluated expressing the conditional probabilities as joint probabilities

p(x n+1 |x n-k , y n-k ) = p(x n+1 , x n-k , y n-k )/p(x n-k , y n-k
) and computing the 4 multidimensional probability density functions (pdfs):

TE Y →X,α (k) = p α (x n+1 , x n-k , y n-k ) log α p(x n+1 , x n-k , y n-k ) p(x n-k ) p(x n+1 , x n-k ) p(x n-k , y n-k ) (3.5)
as a function of time delay k and normalised such that 0 ≤ T E ≤ 1. The 4 pdfs in Eq.(3.5) result from a binning process of times series X and Y , such that Eq.(3.5) is estimated in practice as:

TE Y →X,α (k) = β i=1 β j=1 β l=1 p 3d (i, j, l) α log α p 3d (i, j, l) p 1d (j) p 2d xx (i, j) p 2d xy (j, l) (3.6)
where p 3d , p 2d xx , p 2d xy and p 1d are the discretised versions of respectively p(x n+1 , x n-k , y n-k ), p(x n+1 , x n-k ), p(x n-k , y n-k ) and p(x n-k ). In order to have sufficient statistics, a bin size β = 2 or β = 3 is typically chosen, depending on the available length of the time series (the longer the times series, the larger β can be). We introduced here the additional exponent α ≥ 1, which effectively represents a nonlinear threshold: low probabilities will be further reduced and higher ones amplified. In a complex setting, information may flow both ways, from Y to X and inversely. It is thus especially useful to define the net transfer entropy

∆ X,Y (T E)[k] = T E Y →X [k]-T E X→Y [k]
, which provides the net flow of information between processes X and Y , at timelag k. In the manuscript, pdfs in Eq.(3.5) are discretised using β d = 2 d bins, with d the dimensionality of the pdf. The nonlinear threshold exponent α is set to unity and X and Y are discretised at the same rate and enter the T E calculation with zero temporal mean. Further details may be found in Ref. [Dif-Pradalier et al. 2021] and appendix G.

Vorticity balance

"Ce qui est simple est faux. Ce qui est compliqué est inutilisable"

Paul Valéry

We wish to apply the method above to radial electric field dynamics. To this end, we derive a vorticity balance [Abiteboul et al. 2011, Ajay et al. 2020, Sarazin et al. 2021] from the primitive gyrokinetic-quasineutrality equations, Eqs.(A.1)-(A.3). The following vorticity equation can be inferred from the primitive gyrokinetic equations including E × B drift and finite Larmor radius at leading order:

∂ t ⟨Ω r ⟩ + 1 r ∂ r r ⟨v Er Ω r ⟩ + ⟨v ⋆r Ω r ⟩ -v ⋆θ 1 r ∂ θ E r = r.h.s (3.7) r.h.s ≈ -∂ t ⟨Ω θ ⟩ - 1 r ∂ θ (v Eθ + v ⋆θ )Ω θ -∂ r 1 2r ∂ θ v 2 Eθ + 1 2r 3 ∂ θ ∂ r r 2 v 2 Er - 1 r ∂ θ v ⋆r 1 r ∂ θ v Eθ (3.8) Ω r = ∂ r (r∂ r ϕ) /r & Ω θ = ∂ 2 θ ϕ/r 2 (3.9) v Er = -∂ θ ϕ/r & v Eθ = ∂ r ϕ = -E r (3.10) v ⋆r = -∂ θ p ⊥ /r & v ⋆θ = ∂ r p ⊥ (3.11)
where ⟨•⟩ denotes an average over toroidal angle φ. All terms in Eqs.(3.7) and (3.8) are radially (0.95 ≤ r/a ≤ 1.0), poloidally (0 ≤ θ ≤ 2π), toroidally (0 ≤ φ ≤ π/5) and temporally (0 ≤ t Ω -1 ci ≤ 50, 000) averaged. Their algebraic sum is displayed as "error" in Fig. 3.6-(a) (the value -7 corresponds to a 0.9% error). The terms in Eq.(3.8) are those that are individually small, with r.h.s denoting their sum. Vorticity balance is precisely satisfied (0.9%) during these early stages of radial electric field build-up, which allows to write the leading order vorticity equation:

∂ t ⟨Ω r ⟩ + 1 r ∂ r r ⟨v Er Ω r ⟩ + ⟨v ⋆r Ω r ⟩ -v ⋆θ 1 r ∂ θ E r = r.h.s ≈ 0.
(3.12)

Causality in transport barrier onset

"Le mode-H [La justice] c'est comme la Sainte Vierge. Si on la voit pas de temps en temps, le doute s'installe"

Michel Audiard "Politicians use statistics in the same way that a drunk uses lamp-posts -for support rather than illumination"

Andrew Lang

Radial vorticity ⟨Ω r ⟩ (flow shear) in Eq.(3.12) dominantly evolves through the combined influence of 3 fluxes: the usual E × B (radial) advection of vorticity ⟨v Er Ω r ⟩ [hereafter denoted "Reynolds force"] and two seldom discussed mechanisms: diamagnetic (radial) advection of vorticity ⟨v ⋆r Ω r ⟩ and diamagnetic (poloidal) advection of poloidal inhomogeneities of the radial electric field v ⋆θ 1 r ∂ θ E r [hereafter, "field advection"]. Unexpectedly, in the early stages of radial electric field build-up both latter contributions display magnitudes larger than that of the Reynolds force ]. Dynamical significance of each term to vorticity build-up however does not straightforwardly follow.

To elucidate this point, we systematically apply the T E algorithm Eq.(3.6) to actual time series X and Y from the flux-driven Case-1 computation, including limiter boundary conditions in the last 5% inside the separatrix where the spontaneous onset of a persistent transport barrier is observed. The sharp radial electric field well at the plasma edge ] is clear signature of it. All pairwise possible permutations of quantities in the following set are considered:

(X, Y ) ∈ ⟨Ω r ⟩ , ⟨v Er Ω r ⟩ , ⟨v ⋆r Ω r ⟩ , -v ⋆θ 1 r ∂ θ E r (3.13)
which corresponds to the leading-order terms in Eq.(3.12). Application of T E on electric field data from the primitive flux-driven equations opens captivating possibilities to confront with basic mechanisms. Relevant interactions are portrayed in Fig. 3.6-(c), (d) and (e). We focus on the outer corona (0.95 ≤ r/a ≤ 1.0) and divide the plasma volume in six areas, shown in Fig. 3.6 subplots (b) through (e). We especially focus on areas ❺ (-2π/3 ≤ θ ≤ -π/3), centred about the limiter where the radial electric field well originates and area ❻ (-π/3 ≤ θ ≤ 0) into which turbulence later spreads. Dominant information transfer from systematic pairwise T E computations is diagrammatically represented, for two time intervals. During ∆t 1 (0.04ms to 0.17ms), turbulence is confined to area ❺ [Fig. 5-(d background is responsible for persistent diamagnetic currents, whose importance is not commonly stressed. Interestingly here, diamagnetic (radial) advection of vorticity ⟨v ⋆r Ω r ⟩ proves to be at barrier inception the dominant causal agent ] which directly generates vorticity ⟨Ω r ⟩ and forecasts the Reynolds force ⟨v Er Ω r ⟩.

(ii) This central role of ⟨v ⋆r Ω r ⟩ endures at later times ] even when the Reynolds force (towards which ⟨v ⋆r Ω r ⟩ keeps transferring information) starts contributing more to vorticity production.

(iii) The oft-expected dominant transfer mechanism: Reynolds force causing vorticity production becomes major only at later times ∆t 2 and in area ❻ [Fig. 3.6-(c)], naturally stable without limiter. It has required turbulent fluctuations to have spread there, downstream of production region ❺.

(iv) T E analysis also highlights another (new to our knowledge) mechanism for both vorticity production and shear dissemination: (poloidal) diamagnetic currents contribute to vorticity build-up whilst poloidally propagating radial electric field poloidal inhomogeneities v ⋆θ and (e)]. Farther downstream [Fig. 3.6-(c)] vortical structures ⟨Ω r ⟩ transfer information to field advection. This mechanism contributes to expanding the new-sprung E r well, initially a localized feature of area ❺ and making it a poloidally global feature spanning regions ❶ through ❻.

1 r ∂ θ E r [
Causal analysis of transport barrier build-up emphasizes the central and somewhat unexpected role of diamagnetic flows (v ⋆θ = ∂ r p ⊥ and v ⋆r = -∂ θ p ⊥ /r), shedding light on new or low-keyed mechanisms. This stresses the role of pressure inhomogeneities and finite Larmor radius (FLR) effects in barrier build-up. The importance of the latter could have interesting side effects, foremost on I-mode or H-mode accessibility as different isotopes or different classes of particles (electrons, main ions, energetic particles or impurities) may thus differently contribute to vorticity (shear) production, affecting access to bifurcated states of enhanced confinement.

Some words of caution

"The only statistics you can trust are those you falsified yourself"

Sir Winston Churchill

The bivariate T E method employed above can achieve high sensitivity for short time series (larger e.g. than bivariate mutual information) whilst it is unable to control false positives as the amount of data increases. In that limit, one should target multivariate methods. This is currently a matter of investigation. For this reason, we limited ourselves to the short times series up to about 30, 000 to 50, 000 Ω -1 ci , which amounts to about 0.2 to 0.3ms. This is a short time duration, accessible in modelling but certainly challenging to diagnose from an experimental point of view. Longer times series indeed start to show false positives and we have thus discarded them. The sampling time resolves the evolution of the vortical structures. Within this 0.2 to 0.3ms temporal window, we have further smoothed the complex turbulent dynamics by (i) averaging over the toroidal angle (and have checked that for a sufficiently large averaging, results from the TE algorithm are unchanged) and (ii) repeatedly run the T E algorithm independently at various radial locations or when averaging radially over a typical eddy scale. What is shown here are the robust features from all these analyses.

To the difficult question: "is there a hidden variable that may be driving the whole dynamics of both ⟨v ⋆r Ω r ⟩ and ⟨v Er Ω r ⟩", the short answer is "I do not know". However, as the pressure gradient is not simply a passive scalar, a strong feeling, hard to unambiguously confirm or refute is that the dynamics of density is key, results from a primitive poloidal anisotropy and is largely responsible for the vorticity dynamics that lead to transport barrier onset. This is a personal view at onset of the barrier; the causality in steady state is not as clear. On the one hand the same free energy source keeps pumping density fronts near the cold pulse. On the other hand, in steady-state, through conservation of potential vorticity, dynamics of density fronts and dynamics of vorticity are of course intimately linked.

Towards a fully gyro-kinetic edge & SOL

"Quand les hommes de 100 kg parlent, ceux de 60 kg les écoutent" Michel Audiard, contradicted by sheath physics Sheath physics builds up due to the difference of mobility between heavy ions and light electrons. Fast electron loss to the wall drives the potential drop at the plasma-material boundary interface. At leading order, the boundary electric field which develops at the Debye scale is in the parallel, provided normal magnetic field incidence on the material boundaries. In the above, a Boltzmann response for the electrons has been assumed in the core, edge and SOL. This simplified response allows to penalise the gyrokinetic-quasineutral system, effectively circumventing the requirement to resolve sheath physics. The electric potential can indeed be made to relax in the SOL towards its theoretical fluid prediction. As discussed in section 3.2.1.1, generalisation to a kinetic response of trapped electrons in the confined core is straightforward, provided passing electrons, electrons in the SOL and electrons intercepting a material boundary remain treated adiabatically. This is detailed in Appendix E.

Consistent treatment of electron dynamics however requires a fully kinetic framework, both in the edge and the SOL. In the perspective of a fully consistent gyro-kinetic treatment of the SOL, typical sheath scales are subgrid. Resolved scales are indeed commensurable to the (typically ion) Larmor radius, roughly 2 orders of magnitude larger than the Debye scale. Fine grid resolution at the Debye scale, on global tokamak scales is currently hardly accessible to supercomputing. Furthermore, the gyrokinetic Poisson equation involves the transverse Laplacian of the electric potential, as illustrated e.g. in simplified Eq.(A.3). This poses a non-trivial twofold problem: (i) one wishes to incorporate subgrid sheath physics in gyrokinetics without resolving the relevant Debye scale. Additionally, (ii) sheath physics constrains at leading order the parallel electric field, leaving the potential in the transverse plane, relevant to the gyrokinetic framework, largely unconstrained. Strategies to manage this problem are currently under active investigation. Our approach uses the reduced (1d, 1v) Voice code [Munschy et al. 2023b], a fully kinetic version of Gysela tailored to address sheath physics, at the Debye scale.

Assuming that the above problems can be resolved, one needs to further model particle sources. This is a requirement in order to access genuine steady-states and because interactions between plasma and neutrals in the vicinity of limiters or divertors give rise to recombination and ionisation processes. The reaction rates (magnitude) and anisotropy (radial-poloidal distribution) of such particle sources are known to critically impact heat and particle decay lengths in the SOL as well as deposition patterns on the material boundaries. Various approaches have been studied in the literature -mostly within fluid approximations [Frerichs et al. 2012, Bufferand et al. 2014, Ciraolo et al. 2019, Bufferand et al. 2021, Coroado & Ricci 2022]. They range from coupling to a first-principles Monte Carlo framework such as Eirene [Reiter et al. 2005] or to advanced approaches such as Degas2 [START_REF] Stotler | Neutral Gas Transport Modeling with DEGAS 2[END_REF], Gtneut [Mandrekas 2004] or Neut2d [Shimizu et al. 2003] to more heuristic (and fast) frameworks. These are as well important questions.

3.5 Build with confidence: confront modelling plurality "Debugging is twice as hard as writing the code in the first place.Therefore, if you write the code as cleverly as possible, you are, by definition, not smart enough to debug it" Brian Kernighan -C/Unix guru The problem of transport in the edge region of fusion plasmas is among the difficult and critical issues for the operation of present and coming devices. Advancing along the lines now briefly laid down should help make way with some degree of confidence:

• a large spectrum of models is key. Relevant models span from fully kinetic approaches ({3; 2; 1}d, {3; 2; 1}v) to reduced fluid or gyro-fluid ({3; 2; 1}d) and more heuristic models (e.g. κ -ε, machine-learnt or scaling laws). Reduced models allow to address a diversity of relevant problems: scan parameter space, assess robustness, couple plasma response with atomic physics, run in almost real time in the control room, etc. Physically-motivated (if not systematic) comparisons between models would certainly lead to cross-fertilisation, to the benefit of both basic understanding and practical operation. In particular, one could wish for:

systematic comparisons between different plasma-surface boundary conditions: logical or conducting sheaths, penalisation techniques, neutral models, etc.

the statistically-relevant long energy and particle confinement timescales are key to assess particle and heat loads on material boundaries, to test theoretical ideas (e.g. [Norscini et al. 2022]) or to build large databases to train reduced models (e.g. through Artificial Intelligence techniques [Cartier-Michaud et al. 2016, Heinonen & Diamond 2020, Asahi et al. 2021]). Long timescales and large parameter scans are hardly accessible to the more complex models. A hierarchy of inter-connected models is thus evidently key.

• standardised frameworks can prove important for code verification and to confront with experiments.

It is clear that the so-called Cyclone Base Case [Dimits et al. 2000] has helped structure the gyrokinetic and gyrofluid approaches in the plasma core. Illustration of the critical role of zonal flows near marginal stability, following the work by [START_REF] Rosenbluth | [END_REF]] has been an important landmark.

-In this spirit, a common framework of computations for the edge and SOL has been proposed by Philippe Ghendrih, dubbed the Mistral Base Case for obvious reasons [Dif-Pradalier et al. 2011b]. In particular, assessing the specific (or not) role of kinetic effects in the physics of the edge, SOL and divertor is important.

-The role of open field line geometry (diverted, limited, snowflake, etc.), especially when coupled to atomic physics is another important topic. Reduced approaches are important there for often the only ones accessible with current computing capabilities. Divertor configuration simulations are demanding: the cold plasma that builds-up in the divertor leads to very small Larmor radii. The characteristic scales of turbulent eddies are governed by the Larmor radius, k ⊥ ρ i ∼ 0.3. Very large meshes are thus required, following typically a ρ -3 ⋆ scaling. Tore Supra experiments, such as the framework described in section 3.3 provide an interesting alternative. The circular geometry (i) is simpler and allows one to simplify the magnetic geometry requirements. The (ii) hotter SOL plasmas in a limiter configuration are less demanding with respect to the more extreme ρ ⋆ values of diverted configurations. This directly relaxes requirements on mesh size and computational length. Regardless, (iii) there is a wealth of interesting observations and quality experimental measurements that relate turbulence, transport and the choice of the limiter configuration with very small changes of the magnetic equilibrium. Last (iv) Tore Supra is characterised by significant ripple that can be varied with a proper positioning of the plasma within the chamber. This topic is under active investigation, see e.g. the works of Varennes et al. [START_REF] Varennes | Synergy of Turbulent Momentum Drive and Magnetic Braking[END_REF]. It offers an interesting opportunity to investigate plasma rotation and validate the modelling effort, especially given the critical importance that rotation shear plays in barrier formation at the edge of the confined plasma.

Summary and outlook

"Trop nombreuses les clartés séparées aggravent l'obscurité générale"

Olivier Rey

In the early days of fusion research, where plasma and material boundaries would meet was not the focus of specific attention. This situation has evolved. In part because of the requirement to maintain fluxes of heat or mass incident on the material boundaries at manageable levels and with controlled deposition patterns. In part as well because of the experimental discovery of transport bifurcations, chiefly dependant on the organisation of the peripheral turbulent plasma. This latter question is the focus of our present work. What we establish is threefold:

(i) Plasma-boundary interaction deeply modifies convective stability next to the magnetic separatrix;

(ii) Resulting locally-borne eddies spread out and destabilise distant regions of the edge and core. A globally organised state emerges, 'nonlocally' [Dif-Pradalier et al. 2010, Ida et al. 2015, Hahm & Diamond 2018] controlled by fluxes of turbulence activity;

(iii) Flow shear builds as eddies (vorticity) are advected, primarily through pressure inhomogeneities. The expanding interface organises into a stable peripheral transport barrier, i.e. into improved confinement.

Mechanism (i) above is generic and likely applies beyond limiter configurations to physical cases with radiation near a magnetic X point (a magnetic saddle point). Direct interplay between confined plasma and material boundaries is reduced in X point configurations. Neutral particles thus certainly play an essential role there as they provide localised dissipation in the vicinity of the X point. An additional source of free energy, linked to the penetration depth of neutrals and akin to the free energy source described here for a limiter can thus be expected in X point configurations. Interestingly, where sharp gradients and a radial electric field well are expected from experiments, state-of-the-art gyrokinetic modelling of ASDEX Upgrade plasmas in X point geometry tend to show without plasma-neutrals interplay smoothly varying profiles of density, temperature and radial electric field across the separatrix, at nominal parameters [Michels et al. 2022]. This discrepancy may further hint at the central importance of poloidally-localised free energy injection near the separatrix to initiate edge turbulent dynamics, later made global through spreading (ii). The present argument however does not alone provide an explanation for edge turbulence growth and subsequent confinement improvement in the case where the X-point is not dissipative and the radiation front located further out in the SOL, near the strike points. This case should be carefully considered.

Self-advection (spreading) of turbulent fluctuations (ii) is shown to play a central role in the global equilibration of edge turbulence. The connection between spreading and confinement is complex. On the one hand, spreading contributes to disseminating turbulence activity, spatially, so it is usually thought of as detrimental to confinement. On the other hand, the same dynamics that propagates pressure inhomogeneities is shown to contribute (iii) to vorticity fluxes, thus leading to shear production and being beneficial to confinement. Having both turbulence spreading and barrier formation is less mutually exclusive than it may seem.

This possibility has important implications for modelling. It certainly nudges towards flux-driven frameworks. It also provides a different view on transport barrier formation from a bath of active turbulence. Much focus has indeed been given to understanding the origin of the fluctuations (from which instability do they stem: electrostatic, electromagnetic, at ion or electron scales, at which exact location, etc.). The rationale being that the nature of the instability carries over nonlinearly and may critically influence the route towards improved confinement.

Our results suggest that a complementary perspective could prove fecund. Provided that a source of instability exists, notably near a cold dissipative region, eddies are spawned and turbulent activity spreads. These properties are robust and a priori independent of the nature of the underlying instability. Key questions become: given a distribution of sources (external heating, recycling,. . . ) and sinks (collisional dissipation, interaction with the boundaries) under which conditions is free energy production sustained? What is the branching ratio between flow shear reinforcement and dissemination of turbulence activity? First steps have been provided here and should be expanded as a comprehensive understanding of these questions would likely reshape our understanding of confinement bifurcations, yet a central question for fusion research.

Chapter 4

Conclusion and perspectives

"Whenever you find yourself on the side of the majority, it is time to pause and reflect"

Mark Twain

Microscale turbulence is well known to nonlinearly excite secondary mesoscopic structures such as fronts, avalanches or zonal flows. Such structures, which may appear as mutually exclusive, are yet broadly observed, especially near marginal stability (akin here to the so-called 'Dimits shift' region). Such regimes are likely to be especially relevant in the core plasma of future large, hot and moderately externally-driven experiments such as Iter. These organised structures are sensitive to assumptions of scale separation and are best understood in a flux-driven framework. Self-organisation of these competing mesoscopic structures and their impact on transport raise acute questions. Near marginality, (i) distinct turbulent states and flux levels can be obtained from almost identical underlying mean plasma gradients. This fact underlies a multivalued (bistable mostly) [Diamond et al. 1997], non-monotonic flux-gradient landscape, with competition between mean (equilibrium) flow shear and zonal (fluctuating) flow shear [Dif-Pradalier et al. 2017b]. In addition, (ii) turbulence spreading [START_REF] Garbet | Radial propagation of turbulence in tokamaks[END_REF]] is ubiquitous in flux-driven regimes, its more visible manifestations being the observation of fronts or avalanches. This fact underlies departure from a local or Fickian transport paradigm [Dif-Pradalier et al. 2010]. Connection between spreading and departure from local, single-valued flux-gradient landscapes is still debated. Besides breakdown of Fickian transport at microscales, flux-driven frameworks report (iii) mesoscale organisation near marginal stability with emergence of a tertiary structure, the "ExB staircase" [Dif-Pradalier et al. 2010]. Staircases are a sequence of zonal mean flow layers co-located with profile corrugations which statistically enclose regions of avalanching (strong transport). The zonal mean flow inter-spacing scale is about 40 to 50 ion Larmor radii, with a fat tail [Dif-Pradalier et al. 2017b, Hornung et al. 2017, Milovanov et al. 2021a]. This emergent mesoscale near marginality (iv) is consistent with a transition from Bohmlike (avalanche-dominated) transport for small plasma sizes to gyro-Bohm like (staircaseregulated avalanching) transport for large plasma sizes [Dif-Pradalier et al. 2017b]. Above marginal stability, as staircase organisation fades, avalanching activity is less regulated and Bohm-like confinement can be observed [START_REF] Nakata | [END_REF]. Interplay between layering and spreading, permeability of shear layers (staircase-related or otherwise) to extended transport events are seldom explored though could prove important for a deeper understanding of transport. Regarding robustness, (v) staircase organisation is numerically observed with adiabatic or kinetic [Rath et al. 2021, Qi et al. 2021] electron responses, in electrostatic or electromagnetic [Rath & Peeters 2022] frameworks. Gradient-driven approaches hinder emergence of such structures. Experimentally, density fluctuations measured with fast swept reflectometry on Tore Supra [Dif-Pradalier et al. 2015, Hornung et al. 2017] have yielded the first identifications of regularly-spaced shear flow layers consistent with the staircase phenomenology. Many more machines have since confirmed this trend and shown additional features such as avalanche regulation by shear flows [Qi et al. 2021]. More experiments are required to fully document staircase regimes. Several (vi) models have been proposed to account for staircase formation and have been discussed [Kosuga et al. 2014, Ashourvan & Diamond 2016, Garbet et al. 2021, Milovanov et al. 2021a]. There is still need for a quantitative assessment of the relative merits of each approach with respect to both experimental measurements and flux-driven gyrokinetics. This step is important to gain understanding of these phenomena and to build reduced models with confidence. Last but not least, (vii) preliminary observations have yielded significant discrepancies in heat transport predictions between flux-driven approaches and gradient-driven or quasilinear approaches, in near marginal regimes with observed staircase patterning [Gillot et al. 2022]. These observations are drawn whilst performing a comprehensive comparison of three leading reduced models used by the community -nonlinear local gyrokinetics, local quasilinear and flux-driven quasilinear approaches-to reference flux-driven nonlinear gyrokinetics. Above marginality (these are accessible regimes in current small tokamaks) all 4 models compare favourably; they do not near marginal stability. As experiments increase in size, the plasma response becomes stiffer and (nonlinear) marginal stability becomes practically relevant. In this regime, ion-scale turbulence interplays with equilibrium mean scales and turbulence spreading (especially through avalanching processes) becomes important. In this regime all three reduced models display large and systematic under-predictions of fluxes, flows and transport. Whereby this shortfall arises is comprehensively discussed. A first guess would be to question the assumed assumption of 'quasilinearity' of fluctuations. We show that somewhat surprisingly it is actually robust, despite exceeding the quasilinear domain of validity. Current choices for turbulence saturation rules on the other hand are found to be inadequate near marginal stability. This is good news for model reduction as it allows to provide constructive recommendations as to how reduced models must be improved near marginality. These observations (which e.g. may push towards multivalued flux-gradient landscapes) are still in their infancy and require further investigation, especially whilst relaxing the assumption of an adiabatic response for the electrons. They may yet have important implications for reduced modelling, especially when extrapolating to Iter and to future devices. We all know that boundary conditions or forcing matter, yet as thorny issues their precise influence often remains unquantified. For over three decades, the observation in magnetised plasmas of rapid core confinement improvement upon favourable modifications of edge operating conditions has been a nagging source of puzzlement. Exploring mechanisms whereby bifurcations to improved confinement may occur is important for the safe achievement of fusion energy production on Earth. Transitions to improved confinement are believed to result from spontaneous self-regulation of turbulence in the plasma edge, where confined plasma and material boundaries interplay. Notably in experiments, turbulent plasmas robustly self-organise to higher energy states upon application of additional free energy sources or modification of edge operating conditions. Mechanisms whereby such bifurcations occur have been actively debated for decades. Enhanced confinement occurs in the plasma edge, where a shortfall of predicted turbulence intensity has been puzzling scientists for decades. We have carried out large-scale numerical analyses based on the primitive flux-driven gyrokinetic equation, based on actual experimental parameters from Tore Supra. The computations span core, edge and Scrape-Off Layer (SOL) self-consistently -in particular they include coupling between confined and unconfined plasma-whilst somewhat crudely taking into account plasma-solid interactions. Our results shed new light on 2 long-standing problems [START_REF] Dif-Pradalier | [END_REF]: the shortfall of turbulence power in the plasma edge and the puzzling relative inability from first princi-ples to describe spontaneous onset of improved confinement. We show that both problems are intimately connected: interplay between confined plasma and material boundaries results in a novel source of free energy which triggers separatrix-borne turbulence. Inward self-advection (spreading) of this turbulence fills-in the plasma edge, curing the "shortfall" problem. The turbulent edge in turn organises such as to trigger and sustain the onset of a transport barrier, locking-in an improved confinement state. A comprehensive discussion of the underlying mechanisms is proposed. Further work is needed, especially to include description of the convective nature of the turbulence in the SOL and edge. These results are a first step to addressing highly relevant questions in the quest for magnetic fusion. Interestingly also, they are generic enough to possibly inform problems in fluids and plasmas where turbulence self-advection (spreading/entrainment) competes with transport barrier onset and sustainment -e.g. from the stability of the solar tachocline to the onset and stability of jets in astrophysical or geophysical fluids.

Appendices

Appendix A

Model equations with adiabatic electrons

Low-frequency microturbulence in weakly collisional magnetised plasmas is appropriately described within the gyrokinetic framework [Brizard & Hahm 2007]. The Gysela code [Grandgirard et al. 2016] solves the governing coupled gyrokinetic:

B * ∥s ∂ Fs ∂t + ∇ ∇ ∇ • B * ∥s dx G dt Fs + ∂ ∂v G∥ B * ∥s dv G∥ dt Fs = B * ∥s (rhs) (A.1) rhs = C + S + D -νM mat (r, θ)[ Fs -G cold ] -γ K Fs -F F-D 1 + ⟨ Fs -F F-D ⟩ ⟨F F-D ⟩ (A.2)
and quasi-neutrality: equations for the guiding-center distribution function Fs of ion species s, evolved with no separation between equilibrium and perturbation in five dimensional guiding-centre space (x G , v G∥ , µ) and time. In Eq.(A.2), ⟨•⟩ = J v J x • dv ∥ dµdθdφ, with J v and J x being the velocity and space Jacobians. The charge density of guiding-centers ρ is computed as:

∆n e = ρ + 1 n e0 s Z s ∇ ⊥ • n s0 BΩ s ∇ ⊥ ϕ (A.
ρ(x, t) = 1 n e0 s Z s dµJ µ . J v dv G∥ ( Fs (x, v, t) -Fs,eq (r, θ, v G∥ )) (A.5)
with J µ . the gyro-average operator. Notations are those of Ref. [Grandgirard et al. 2016]. The computational domain extends from inner core (r/a = 0) to the material boundaries (r/a = 1.3). Flux-or gradient-driven dynamics may be considered. For flux-driven evolution, γ K = 0 and the distribution function evolves according to volumetric sources S [Sarazin et al. 2010] and penalised heat and momentum sinks M mat (r, θ), M SOL (r) and M wall (r) that can mimic from poloidally-uniform boundary conditions (Case-2) to the more complex limiter and wall geometries (Case-1). The latter case allows description of the closed to open field lines transition in the Scrape-Off Layer (SOL). Gradient-drivenlike dynamics may also be considered whilst imposing S = M mat (r, θ) = M wall (r) = 0.

In Case-3, the target distribution function F F-D is chosen as the statistical distribution at equilibrium from flux-driven Case-2 and the relaxation rate γ K = 5.43 10 -5 ∼ γ lin /10 is an order of magnitude smaller than the local linear turbulence growth rate γ lin at r/a = 0.7. Imposing M SOL (r) as in Case-2 or cancelling this mask does not alter the dynamics which is dominated by the BGK operator [last term of Eq.(A.2)], specifically The magnetic surfaces in the MFE system are to a first approximation concentric torii with the B = Bz field direction at a small angle with the toroidal direction φ. The so-called safety factor q, ratio of the toroidal to the poloidal magnetic flux is a measure of how the magnetic field lines wind up on the flux surfaces, weaving the magnetic geometry. It is compactly expressed as the pitch of the field q ∼ dφ/ dθ for non-rational q values, θ being the poloidal angle. The safety factor has a fixed value on each toroidal surface so is essentially a function of radius. Rational values of q(r) correspond to surfaces with closed field lines and are especially vulnerable to electromagnetic pertubations. Electrostatically too q rationals play an important role in linear theory as turbulent modes resonate at these locations, a reason that often leads to promoting a picture of "radially pinned turbulence". Low-order rationals have sometimes been shown experimentally to favour the onset of transport barriers. On the basis of these linear or quasi-linear pictures, low order q rationals are often believed to play a critical role in turbulence organisation. Asking about the connection between these rational surfaces and the localisation of the staircase steps is certainly natural and addressed in Ref. [Dif-Pradalier et al. 2017b]. The connection between low order q rationals and staircase steps is not ubiquitous -this point is further addressed in section ix. Bottomline is that the plasma staircase in MFE appears as a genuinely nonlinear phenomenon with memory of the linear proclivities lost or largely subdominant.

This complex magnetic geometry is embedded in the gyrokinetic results of sections A and on; for the illustrative purpose of the current section however the magnetic field is uniform so that the toroidal φ and B directions coincide. The poloidal plane (r, θ) in HM is thus the analogue of the horizontal β-plane, the radial direction r being the MFE analogue of the longitude y and the latitudinal (zonal) direction x corresponding to the plasma poloidal direction θ. For both systems the horizontal momentum equation reads:

dv dt + R × v = -∇U (B.1)
where formally the Larmor motion R ≡ eB/m in the MFE case corresponds to the Coriolis term R ≡ f z in GFD. Similarly the electrostatic potential ϕ that enters U ≡ eϕ/m is the analogue of the shallow water layer depth H for the GFD system where U ≡ gH, g being the constant of gravity. An equivalent description in the GFD case in terms of the streamfunction Ψ = gH/f 0 would similarly lead to considering the fluid streamfunction as the analogue of the plasma electrostatic potential. Quasi-geostrophic equilibrium is assumed in the GFD system meaning that (i) the Coriolis term f z × v in Eq.(B.1) dominates the advective terms, the ratio of the two being the (small) Rossby number Ro = U/f 0 L, where U and L are the typical horizontal velocity and length scales. The smallness of the Rossby number [Ro ∼ 0.1 for the Earth atmosphere and ∼ 0.01 in the ocean] is the translation of the fact that planetary rotation is important, which is equivalently stated as

Ω f0 ∼ 1 f0 ∂ ∂t = o(1)
. It is further assumed that (ii) variations in the Coriolis parameter are small, that is |βL| ≪ f 0 . Furthermore, the variations in fluid depth are small compared to its total depth so that an expansion H = H 0 + h with h/H 0 = o(1) is possible, H 0 being the mean layer depth and h the depth fluctuations. More fundamentally it means that the scale of the motion is of the order of the Rossby radius of deformation L d = √ gH 0 /f 0 which is the typical scale at which rotational effects become as important as buoyancy or gravity wave effects in the evolution of the flow.

The ion Larmor radius ρ i in MFE will play a comparable role to the Rossby radius L d : the typical turbulent motions will take place at scales commensurable with ρ i and will be slow compared to the ion Larmor gyration so that

Ωv,i ωc,i ∼ 1 ωc,i ∂ ∂t = o(1), where Ω v,i = ∇ × v is
∇ • v = - d dt log X (B.2)
where X ≡ H for the shallow water and X ≡ n in the HM framework, the momentum equation (B.1) can be recast in the GFD case as the Charney equation:

∂ ∂t -(∇h × z) • ∇ h -∇ 2 h + log H 0 f planetary PV = 0 (B.3)
and in the MFE case as the HM equation:

∂ ∂t -(∇ϕ × z) • ∇ ϕ -∇ 2 ϕ + log n 0 ω c,i plasma PV = 0 (B.4)
Both equations (B.3) and (B.4) can be compactly written as dQ/dt which expresses the material invariance of the quantity Q -i.e. the conservation of Q on each fluid elementknown in both cases as the potential vorticity (PV). Several expressions can be encountered for the quasi-geostrophic PV in rotating shallow water equations on the β-plane. If focusing on the streamfunction Ψ, the planetary PV is often recast as βy + ∇ 2 Ψ -L -2 d Ψ with L d the Rossby radius of deformation. Similarly in MFE, whilst going beyond the simple HM generic attractor for such classes of systems.

We also show that the possibility for the system to self-adjust is key to adequately describing this aforementioned patterning: when interfering too strongly with the possibility for the system to self-adjust, especially whilst over-constraining the system through a tight forcing or enforcing ad-how scale separations the staircase pattern is deeply altered -see Fig. 2.6 and Appendix C. This may have significant consequences on our modeling practices, see e.g. section 2.3.3.

B.1.2 Constitution of a staircase

The patterning of a turbulent fluid into a staircase structure requires three major ingredients. One is a variation on the key idea first introduced by Phillips in the context of ocean circulation [Phillips 1972]: the so-called "Phillips effect" describes the spontaneous layering of a stably stratified fluid (a stable vertical distribution of the fluid density), in which vertically homogeneous stirring produces vertically inhomogeneous mixing of the background buoyancy gradient. In other words, starting with a constant density gradient and supposing that there is in an ocean a turbulent transport mechanism that produces a vertical flux of buoyancy initially constant with depth, the question asked is what may happen to a local disturbance of the background density gradient. The key timescales are that of the density gradient evolution and the response time for the transport of buoyancy. An instability may occur if the propensity of buoyancy to be transported [its 'diffusivity' if the transport of the buoyancy field obeys a Fick's law] is reduced in the steeper gradient region faster than the density gradient increases. In that case the local decrease in flux must be compensated by a tendency to further increase the gradient. Inversely, increased mixing occurs in the region of flattened gradient so that the flux increases causing a convergence of buoyancy near the interfaces and a step-like density profile to form.

The above mechanism, leading to the spontaneous formation of layers in a mixed fluid has been lent strong support by classic experiments in a controlled environment [Ruddick et al. 1989]. In the planetary case, a Phillips-like effect has been advocated [START_REF] Dritschel | Multiple Jets as PV Staircases: The Phillips Effect and the Resilience of Eddy-Transport Barriers[END_REF] to explain differential PV mixing along stratification surfaces, a first necessary step to explain spontaneous jet formation in the atmosphere and the constitution of what is often referred to as "eddy-transport barriers", clearly documented through the observation of chemicals in the upper atmosphere.

In a plasma we would argue that the radial direction of quasi-constant heat flux flowing from core to edge is the analogue of the vertical direction of constant buoyancy flux. The Phillips mechanism translates into the fact that a homogeneous turbulence acting on a similarly homogeneous base state [a smooth background temperature gradient for instance] may naturally lead to the development of inhomogeneities in this mean temperature gradient. In plasmas, any variation on this mechanism may robustly lead to spontaneously-generated inhomogeneities in the mean background gradients themselves.

The so-called "flux-driven" modeling approach in fusion plasmas acknowledges the facts that the system (i) should be forced by a constant flux and that (ii) the profiles themselves should freely evolve. Given these premisses, strong inhomogeneities in the background plasma gradients are repetitively observed in our calculations. These inhomogeneities have been reported as responsible for a departure from the oft-invoked neoclassical poloidal rotation [START_REF] Dif-Pradalier | [END_REF]] -they could contribute to the onset of a transport barrier [Sarazin et al. 2006]-or for the 'containment' of nonlocal heat transport to mesoscales [Dif-Pradalier et al. 2010] (e.g. a genuine macroscale avalanching is found to be rare) -the connection between these inhomogeneities and the organisation of turbulent transport in tokamaks is further discussed especially in section 2.3.2.2. In the following, we shall refer to these inhomogeneities on the mean [or "background"] plasma gradients as "corrugations". A quasi-regular patterning of these corrugations is the most visible manifestation of the plasma staircase.

The possibility of corrugations happening is one thing; them robustly surviving on temporal scales much longer than the local eddy timescale in a strongly turbulent medium is another. To this end some measure of positive feedback, mediated by shear, between these profile inhomogeneities and flows and stresses is necessarily required. Momentum budgets are key to understanding this point, a point well emphasised in the atmospheric context. Two additional ingredients are necessary to understand jet self-sharpening, eventually leading to staircase formation. The first one is the so-called "PV invertibility principle" [Hoskins et al. 1985].

From Ertel's theorem dQ/dt = 0, PV is dynamically conserved and represents the prognostic quantity from which all other dynamical observables may be deduced. In other words this quantity contains everything there is to know about the fluid advective nonlinearity. The possibility of its "inversion" is the second ingredient necessary for staircase formation. It means that there is a quantity (PV) and an inversion operator such that at every moment in time all the other fields -including in the planetary context the velocity, as in Fig. B.2 and buoyancy fields-can be derived diagnostically from the PV field at that instant. In a plasma, the full [kinetic or gyrokinetic] distribution function is a material invariant and contains prognostically all the information regarding the other plasma fields and as such is the analogue of the planetary PV field. In particular, the velocity and pressure fields are diagnostically obtained as fluid moments of it. Given this fact, let us ponder the implications of a Phillips mechanism on an initially smooth planetary or plasma PV profile.

Inhomogeneities in the mixing of PV are naturally created. Through the inversion process above it means that many key observables are also going to be differentially mixed, including the momentum-related variables. To illustrate this let us for now focus on the GFD case and classically consider the small-amplitude expansion of the velocity and PV fields:

v = U 0 (y) x + u 1 x + v 1 y + o(ε) (B.5) Q = Q 0 (y) + q 1 + o(ε) (B.6)
where Q 0 and U 0 are respectively the planetary O(1) large-scale PV and flow profiles expected from the Earth rotation. The primed quantities are first-order in the small parameter ε and second-order terms are unspecified. The material invariance of Q can be recast as:

D t q 1 + v 1 ∂ y Q 0 = 0 (B.7) with D t = ∂ t + U 0 ∂ x .
This expresses the fact that the q 1 perturbation through advection in the meridional y direction is not conserved along the leading-order trajectories in the presence of a background vorticity gradient ∂ y Q 0 . In the case where U 0 is constant on the β-plane, the planetary vorticity gradient simply is

∂ y Q 0 = β ≥ 0. Plane-wave solu- tions Ψ 1 = Ψ exp [ı (kx + ly -ωt)],
where Ψ 1 which satisfies ∇ 2 Ψ 1 = q 1 is the perturbed streamfunction naturally yield to finding the Rossby waves whose dispersion relation reads:

ω = U 0 k - βk k 2 + l 2 (B.8)
Rossby Waves thus owe their origin to the gradient of the tangential speed of the planetary rotation (the planetary vorticity) i.e. when β ̸ = 0. In particular zonally-travelling Rossby waves (U 0 = l = 0) are backward waves, i.e. ones in which the phase (momentum) velocity -β/k 2 is opposite to the group (energy) velocity +β/k 2 . Consequences of this fact are discussed below. At each order in the small-amplitude expansion Eqs.(B.5)-(B.6) the quantities can be further decomposed into a zonally averaged (average along the latitudinal x direction) mean part noted with a bar and a disturbance noted with a prime. As the planetary vorticity

[resp. flow] has no disturbance Q 0 = Q 0 [resp. U 0 = U 0 ] and inversely q 1 = q ′ 1 [resp. u 1 = u ′ 1 and v 1 = v ′ 1 ]
. Noting ℓ the elementary particle displacement in the y-direction which satisfies D t ℓ = v 1 , Eq.(B.7) can be recast as: D t (q ′ 1 + ℓ∂ y Q 0 ) = 0 where ℓ is first-order and ℓ ≡ ℓ ′ 1 . If initialised as ℓ = -q ′ 1 /∂ y Q 0 , this relation holds at later times so that the zonal averaging of Eq.(B.7) gives: D t p = v ′ 1 q ′ 1 where an important new quantity p is defined:

p = -q ′2 1 2∂ y Q 0 (B.9)
which is dimensionally a velocity, is called the zonal pseudomomentum per unit mass and is always negative for Rossby waves. The plasma analogue of this quantity straightforwardly exists in the MFE context [Diamond et al. 2008]. This quantity is second-order in ε and is usually called a 'wave property' in GFD as it is completely determined by the first-order ℓ and q ′ 1 quantities. It is defined as the negative of the wave activity density which being quadratic in the wave amplitude is a measure of the Rossby wave intensity. Physically the pseudomomentum thus represents the momentum per unit mass carried by the Rossby [drift in the plasma case] waves, which as shown below directly controls the rate of change (the acceleration) of the zonal (i.e. latitudinal in GFD and poloidal in MFE) mean flow. Remarking that D t p = ∂ t p we get the important pseudomomentum evolution equation:

∂ t p = v ′ 1 q ′ 1 (B.10)
The flux of potential vorticity v ′ 1 q ′ 1 is dimensionally a force; this is readily seen through the Taylor identity:

-v ′ 1 q ′ 1 = ∂ y v ′ 1 u ′ 1 (B.11)
the right-hand side of Eq.(B.11) being the Reynolds force.

Starting from the material invariance of PV and allowing inhomogeneities in the mixing of PV to develop through a Phillips mechanism, one naturally gets the development of Rossby [drift] waves and the generation of stresses mediated by these waves. The last useful piece of information to close this loop is to ask about the mean flow response to the waves. Substituting the exact decomposition v = (u + u ′ )x + (v + v ′ )y into the zonally-averaged continuity equation (B.2) reads ∂ x u + ∂ y v = 0 which injected into the momentum equation leads to the seeked (exact) evolution equation for the mean latitudinal (x direction) flow u:

∂ t u = -∂ y v ′ 1 u ′ 1 (B.12)
which expresses the fact that the meridional (i.e. y-direction) flux of zonal momentum drives the zonal (i.e. latitudinal in the x-direction) mean flow rate of change. Equations (B.10), (B.11) and (B.12) are the key results and have far-reaching implications. One readily sees that up to second-order terms in ε, the growth of the Rossby wave pseudomomentum equates the zonal (latitudinal) mean flow acceleration:

∂ t u = ∂ t p + o(ε 2 ) ⇒ u = p + U 0 + o(ε 2 ) (B.13)
At a Rossby wavefront p is non-vanishing, negative and given by Eq.(B.9) whereas elsewhere p = 0. The mean flow is thus decelerated as a Rossby wavefront arrives and accelerated back to its mean initial value U 0 as the wavefront leaves. This is the usual so-called "nonacceleration theorem" of GFD [Vallis 2006]: in the absence of dissipative body forces, (i) only transient waves can accelerate a mean flow and (ii) this effect is reversible so that a long-lasting influence of the waves on the flow is out of scope. Jet self-sharpening and staircase formation require long-lasting irreversible effects on the mean flows. The "non-acceleration" picture radically changes when accounting for some degree of dissipation. Dissipation or wave breaking in the celebrated "surf zones" [McIntyre & Palmer 1983] of the atmospheric context is thus the third required ingredient. Including forcing and dissipation is usually done whilst introducing a body force F in the right-hand side of the momentum equation Eq.(B.1) so that the material invariance of Q Eq.(B.7) now reads:

D t q 1 + v 1 ∂ y Q 0 = F 1 (B.14)
where F 1 is first order in ε. The pseudomomentum is still defined as in Eq.(B.9), is negative and its evolution equation is now recast as:

∂ t p = -∂ y v ′ 1 u ′ 1 - q 1 F 1 ∂ y Q 0 (B.15)
As ∂ y Q 0 ≥ 0 regions where q 1 F 1 ≥ 0 are regions of wave forcing whereas regions where q 1 F 1 ≤ 0 correspond to regions of wave dissipation. Eq.(B.12) remains valid provided the body force does not add mean zonal momentum to the system so that the mean flow acceleration now reads:

∂ t u = ∂ t p + q 1 F 1 ∂ y Q 0 + o(ε 2 ) (B.16)
Let us now assume for clarity's sake that a steady train of Rossby wave is emitted in the system so that ∂ t p = 0. We now get the interesting result that in a dissipative region ∂ t u ≤ 0 whilst ∂ t u ≥ 0 in a forcing region.

A zonal (latitudinal) mean flow may thus be forced (i) as in the non-dissipative case through wave transience ∂ t p or (ii) may be accelerated in the Rossby wave forcing region due to wave breaking even though the body force has been assumed not to input any zonal momentum. Interestingly also if q 1 F 1 is weakly dependent on time, Eq.(B.16) admits as approximate solution:

u ≈ p + U 0 + q 1 F 1 β -1 t (B.17)
a superposition of the non-dissipative "non-acceleration" mean flow and of an additional term which gives secular growth of the mean flow in the forcing region where q 1 F 1 > 0. At long times this solution is no longer valid of course yet Eq.(B.17) points towards a result of broad significance that can be proven to hold in a more rigorous perturbation analysis [Bühler 2009]: even small-amplitude waves can contribute in the course of time O(1) mean flow changes. This result is central to obtaining perennial turbulence-generated mean flows [i.e. jet self-sharpening] and staircase patterning. Interestingly also the introduction of any amount of dissipation in the GFD or MFE system alike (eddy-eddy friction, collisional dissipation, irreversibility through boundary conditions, etc.) leads to the irreversibility of the mean flow changes. Where previously the effect of a wavefront propagating through a shear flow was reversible, in the presence of dissipation even a single wavefront will lead to irreversible changes in the rate of change of the mean flow.

The above discussion is yet another way to emphasise the backward nature of the Rossby [and drift] waves previously mentioned: even with a body force that does not input zonal momentum [i.e. the net change in total momentum is zero] the interaction of turbulence [e.g. at mid-latitudes] and wave energy propagation [northward/southward, i.e. away from the region of stirring] in the presence of the large-scale background PV gradient imparted through the planetary rotation leads to: (i) the convergence of momentum into the region of stirring [the region of Rossby wave generation] producing there a net eastward flow in the case of the Earth complemented by (ii) a westward flow in the dissipation region in the flanks of this mean eastward flow that ensures momentum conservation -see Fig. B.2. The reader could interestingly compare the proximity of this result with the radial profile of poloidal velocity displayed in Fig. 2.3 (third from top) next to one of our computed plasma staircase steps.

B.1.3 Summary and further discussion

There are plasma analogues -see Table .B.1.1 to all latter introduced quantities, at each step of the discussion: drift and Rossby waves share the same structure for the dispersion relation, a "plasma Phillips effect" is robustly observed in flux-driven modeling as the turbulence spontaneously tends to organise dynamically into regions of differential mixing, creating corrugations and nucleating mean E × B flows [Dif-Pradalier et al. 2010], material invariants exist in both systems, a plasma pseudomomentum p can be worked out [Diamond et al. 2008], the zonal (latitudinal) mean flow u is the flux-surface averaged E × B plasma mean flow and a Taylor relation Eq.(B.11) has been proven for gyrokinetics [McDevitt et al. 2010]. All these facts tend to emphasise the above discussion as in fact quite generic provided the correct analogies are made.

An overall picture emerges where starting from a latitude-dependent monotonic PV distribution (it monotonically increases on Earth from the South pole to the North pole) (i) a Phillips-like effect starts generating a corrugated PV profile. The pair of a locally flattened (mixed) PV region next to a sharper (unmixed) PV gradient (ii) through PV inversion leads to the nucleation of a jet at the sharper PV gradient location. Regions of substantial cross-stream PV gradients as e.g. on Earth in the vicinity of the main circumpolar vortex [McIntyre & Palmer 1983] are regions able to sustain strong Rossby wave propagation. These large PV gradients are mostly found on the side of the mixing regions (i.e. the flanks of the zonal jets). Such jets (iii) act as waveguides for Rossby waves which as they propagate dynamically organise the fluctuating fields in the sense of introducing systematic correlations between them, shaped by the polarisation of the waves [START_REF] Dritschel | Multiple Jets as PV Staircases: The Phillips Effect and the Resilience of Eddy-Transport Barriers[END_REF]. This dynamic organisation is accompanied by systematic mean effects as Rossby waves propagate on long distances, generating long-range momentum transport [long-range radiation stresses]. These stresses mediate the wave-mean-flow interaction and redistribute the wave momentum with a constant sign as Rossby [drift] waves are backward waves. Increasing with the strength of these stresses is the resilience of the material contours to sideway undulations; in other words a positive feedback on the jet is envisioned where as the local PV gradient increases, further supporting Rossby wave propagation the material contours of the jet are strengthened through increased wave elasticity, further strengthening it. Conversely, strongly-mixed PV regions will support weaker PV gradients; these regions, having weakened wave elasticity will become easier to mix. Based on these arguments, a spatial inhomogeneity where PV is turbulently mixed (i.e. a pair of a locally flatter next to a locally sharper than average PV gradient) on an otherwise large-scale uniform PV field ? andB.2 which may be regarded as a natural outgrowth of the resting planetary [plasma] PV profile, unstable in the presence of Rossby [drift] waves and instabilities. The associated mean zonal wind profile reacts to this and organises the momentum fluxes so as to maintain the staircase structure. This property emphasises the tight link between PV mixing and jet self-sharpening: mixing in the jet flanks suffices to get jet self-sharpening which is key to a perennial staircase. Basic to much of this is the catalysis of potential-vorticity mixing by a Rossby-wave radiation stress together with the one-wayness of Rossby waves.

The essential result is that a mean flow can be maintained by a turbulent stirring [that imparts no net momentum to the medium] in the presence of a gradient of a material invariant [a differential mixing in PV, i.e. a flat mixed region next to a steep PV gradient]. This gradient is provided by the differential rotation in GFD and exists by construction in MFE. Importantly some amount of dissipation is required to irreversibly lock-in the momentum changes in the mean flow due to wave momentum deposition. At last, dissipation and forcing regions spontaneously separate as Rossby [drift] waves, generated in the forcing region propagate meridionally [radially] before dissipating, thus broadening the dissipation distribution. Should they not do so forcing and dissipation would occur at the same location and the "non-acceleration" theorem would hold. These waves are intrinsic to both systems, are efficient carriers of momentum fluxes and mediate a nonlocal wave-induced momentum transfer between the forcing and the dissipation regions.

To finally close the loop an explanation is required as to how turbulence may arise, allowing for mixing and homogenisation of the material invariant -the PV field. The key insight in the Earth context was provided by the observation that as do waves approaching a beach, Rossby waves also break [McIntyre & Palmer 1983] providing for an irreversible deformation of the material flow contours. These are the celebrated "surf zones", the analogues of the "treads" of a staircase (i.e. the flat PV regions in Fig. B.2) which physical reality and their tendency to mix PV has been verified in a vast number of observational and modeling studies, including studies of the stratospheric ozone layer. The steep sections between the treads of the staircase -the "risers"-are where the jets [the plasma mean flows] are located.

The mechanisms above are very generic and therefore attractive. It leads us to a twofold expectation:

• if jets exist initially in the system, the positive feedback on the flows that results from the competition between wavelike and turbulence dynamics can make these jets endure and reinforce;

• if jets do not exist initially and the initial system is homogeneous, even when stirred homogeneously may one expect the spontaneous emergence of long-lived, selfreinforcing jets that will assuredly be associated to large-scale shear straining patterns as these would amplify the corrugation's resilience and amplify the positive feedback on the flow.

Given the close analogies between the plasma and the planetary cases even though the precise mechanisms may differ, our observation of the E × B staircase in tokamak plasmas should appear as less of a surprise. This fact also points towards the interesting new idea that long-lived, self-reinforcing jets may well constitute one aspect of a robust class of attractors for plasmas, a point that has not so far attracted the attention it may deserve, especially with respect to transport barrier formation. In practice, in order to compare flux-and gradient-driven approaches, the target distribution F FD in Eq.(C.2) is constructed from a steady flux-driven distribution function. The following procedure is applied that allows Gysela to mimic usual computations at fixed gradient:

(i) a flux-driven computation is run until flux equilibrium. Mean density n FD , temperature T FD and flow profiles are computed;

(ii) corrugations are removed from the above flux-driven mean profiles resulting in smooth equivalent mean profiles;

(iii) an equivalent Maxwellian distribution F FD is built from these smooth profiles and imposed in Eq.(C.2); (iv) linear stability is computed for the base state F FD . Unusable global modes grow with a maximum growth rate γ lin ;

(v) the "spring constant" γ K in Eq.(C.2) is chosen such that it has no notable effect on the linear instability growth of the ITG and nonlinearly such that the root mean square (RMS) of the electric potential fluctuations (δΦ 2 ) 1/2 , a measure of the turbulence intensity is comparable in both flux-driven and gradient-driven approaches. This is shown in (vi) the gradient-driven calculation is compared to the reference flux-driven computation at quasi steady-state for the flows and the fluctuations.

Through this procedure, flux-and gradient-driven computations have near-identical ambient mean gradients. The curious reader may also have a look at Appendix D where infinite penalization is described as compared to a situation where a strong restoring force of the Krook type is used (i.e. where γ K in Eq.(C.2) would be large). This technique is used to implement a limiter in Gysela, as detailed in section 3.2.1.2.

Note that the gradient-driven approach of Eq.(C.2) contains additional physics as compared to usual gradient-driven approaches. In the latter a scale separation is usually performed in the equations so that formally the total f = F 0 +f and H = H 0 +h are divided into a background distribution F 0 and equilibrium Hamiltonian H 0 , fixed in time and fluctuating parts f and h. The mean evolution ∂ t F 0 is discarded and so are nonlinear advection terms that contribute to its evolution, formally written ⟨[h, f]⟩, with ⟨•⟩ representing an ensemble averaging. Only f dynamically evolves and the consistency of the performed separation of scales f/F 0 = O(ε) has to be checked throughout the computation. In Gysela, the total f and H enter Eq.(C.1) and evolve so that the "mean distribution" ⟨f ⟩ [regardless of whether such a "mean" is correctly approximated by an analytic (Maxwellian) F 0 ] tends to change and corrugate due to the turbulent dynamics-essentially due to the back-reaction of the small turbulence scales ⟨[h, f]⟩ on the larger, slower mean ones. Then the Krook-type operator is applied, partly restoring the mean profiles. The gradient-driven approach of Gysela described by Eq.(C.2) thus contains features of turbulence self-organisation absent when a scale separation is performed. It is thus intermediate between flux-driven and more usual δf global gradient-driven approaches. Actual differences between flux-driven (necessarily full-f) and δf GD approaches, especially regarding structure formation are thus expected to be larger that the differences reported below.

C.2 Mesoscale organisation vs details of forcing models

As expected, the discrepancies highlighted above and e.g. in Fig. 2.6 between flux-driven and gradient-driven systems regarding pattern formation, shear flow generation and transport levels decrease with decreasing strength of γ K in Eq.(C.2) and tend to (transiently) vanish in the limit γ K → 0. Mean profile corrugations, mean shear patterns, more coherent avalanche-type events and staircase organisation tend to appear in gradient-driven systems, as transients, at very low values of the restoring force e.g. γ K ∼ γ lin /100. However at such low values of γ K the mean profiles relax and turbulence decays as its drive is no longer sustained. This remark however tends to emphasise E × B staircase organisation as a natural tendency for the self-organising core plasma.

It is also worth noticing that both decaying turbulence computations as well as computations run with thermal baths at both radial boundaries (this represents a cruder and earlier version of the heat flux-driven computations shown in the paper and illustrated in Fig. C.1, as the amount of heat actually injected in the system is not controlled) also exhibit staircase organisation, though transiently in the decaying case. The robustness of this form of organisation with respect to forcing also tends to emphasise staircase-avalanche interplay as a natural tendency near marginal (nonlinear) stability for such systems, as soon as consistent feedback between flows, shear and mean profiles is allowed.

The second term in the LHS accounts for the adiabatic response of electrons. In particular, β is defined as β (ψ, t) = e ⟨ n adia. e represents the density of electrons that are supposed to have an adiabatic response.

The term γ (ψ) = 1 -M SOL (ψ) corresponds to a modification of the adiabatic response in presence of a scrape-off layer, see section E.3. Finally, the RHS of Eq.(E.1) is expressed as ρ (ψ, θ, φ, t) = i Z i δn Gi (ψ, θ, φ, t) -δn kin. [The paragraphs below are not currently implemented in Gysela. This might be useful to consider for long simulations when the equilibrium profiles have time to evolve, i.e. typically on the confinement time.] To improve the numerical accuracy during the computation of the RHS, one can remove centered Maxwellian distribution functions for both ions and electrons with the constraint i Z i nGi (ψ, θ, t) = ne (ψ, θ, t). The procedure is the following. First compute nGi (ψ, θ, t) (E.9)

= µmax 0 dµ v ∥,max -v ∥,max dv ∥ J v J 0,i • Fi ψ, θ, v ∥ , µ i , t (
Finally, one has Fe = i Z i nGi /η kin. e fe . To limit the numerical cost of this all procedure, one can use the densities and temperatures of the previous time step for the construction of Fi and Fe , as these distribution functions do not need to be exact approximation of the actual distribution function. They just need to be close enough.

E.2 Padé approximation of the quasi-neutrality equation

The Padé approximation of the quasi-neutrality reads [Lanti et al. 2016]:

-∇ ∇ ∇ ⊥ • (α∇ ∇ ∇ ⊥ ϕ) + β (ϕ -γ⟨ ϕ ⟩ FS ) -∇ ∇ ∇ ⊥ • {κ∇ ∇ ∇ ⊥ [β (ϕ -γ⟨ ϕ ⟩ FS )]} = ρ -∇ ∇ ∇ ⊥ • [κ∇ ∇ ∇ ⊥ ρ] (E.10)
where κ = i ρ 2 i = i miTi q 2 i B 2 .

E.3 The hybrid electron model

In the previous section, the general expression of the quasi-neutrality has been derived. In this section, we detail the hybrid electron model where trapped and passing electrons are treated differently, leading to dedicated expressions of n kin. . Moreover, a specific treatment is required in the closed field line region for axi-symmetric modes n = 0, including -but not exclusively -the zonal (i.e. flux surface averaged) mode [Lanti et al. 2020]. Additional ideas on this variation are also detailed in section E.3.3.

E.3.1 Definition of the kinetic and adiabatic electron populations

In the confined region of the hybrid electron model, trapped electrons are considered as kinetic and passing electrons are adiabatic. In the scrape-off layer two different treatments can be adopted: where B LIM (r, θ) corresponds to the amplitude of the magnetic field on grid points corresponding to the limiter. The trapped fraction is then computed with regular expressions. We see from Eq.(E.12) that the position of the limiter will have a direct impact on the quasi-neutrality which is novelty as compared with the adiabatic electron model. The RHS of Eq.E.4 is modified to take into account the response of the adiabatic electron population as ρ (ψ, θ, φ, t) = i Z i δn Gi (ψ, θ, φ, t) -δn kin. In the perspective of the implementation of sheath-like boundary condition in the limiter region, the hybrid electron model is modified to ensure that in the Scrape-Off layer region all electrons are treated as kinetic. There is thus no adiabatic electron response in this part of the plasma. The quasi neutrality equation is identical to Eq.(E.1), and the guiding center charge density ρ expressed as in Eq.(E.4). The kinetic treatment of the Scrape-Off layer is directly embedded in the trapping condition that is now expressed as

v ± ∥ (ψ, θ, µ) = ± 2µ [B max (ψ) -B (ψ, θ)] (1 -M SOL ) ± v ∥,max M SOL , (E.16)
where v ∥,max is the largest parallel velocity resolved in the simulation. In the core plasma the hybrid electron model is unchanged, with the trapped particles treated as kinetic and passing particles as adiabatic. On the contrary in the Scrape-Off layer region the adiabatic density of electrons vanishes while the kinetic density of electron is equal to the total density.

Independently of the chosen electron model to describe the Scrape-Off layer region, the density of kinetic electrons is defined as where L φ is the length of the tokamak in the toroidal direction, J v = 2πB * ∥s /m s is the jacobian in velocity space and where F e is the electron distribution function. The density of adiabatic electrons is deduced from the the kinetic and the total electron density:

n adia.

E.3.2 Treatment of the zonal mode

One of the main difficulties of the hybrid electron model for GYSELA is the peculiar treatment of the zonal mode {n = 0, m ⋆ = 0} which is described with a full kinetic electron response (in the core) as described in [Lanti 2018]. For GYSELA which uses the geometric poloidal angle θ instead of the intrinsic angle θ ⋆ , this can lead to a difficulty. The simplest way is to solve the quasi-neutrality twice, one with the adiabatic and kinetic densities described above ϕ T KE , and another one with the full kinetic response ϕ F KE . The total potential is then One therefore needs to express dθ⋆ dθ which depends on the geometry. In the case of adhoc circular geometry, the relation between θ and θ ⋆ is the following

θ ⋆ = 2 atan √ 1 -ε √ 1 + ε tan θ 2 (E.22)
with ε = r R0 . In this specific case, one has

dθ ⋆ dθ = √ 1 -ε 2 1 + ε cos θ (E.23)
In the case of Culham equilibrium, one has (cf note of Kévin on the Culham equilibrium): Y dφ/L φ , while brackets ⟨Y ⟩ F S refer to flux surface averages. Notice that ⟨Y ⟩ F S = Y m⋆=0,n=0 , where (m ⋆ , n) stand for the Fourier wave numbers related to the (intrinsic) poloidal θ ⋆ and toroidal φ angle, respectively.

θ ⋆ = θ -∆ ′ + r R 0 sin (θ) - 1 2 E ′ - E r sin (
The polarization charge density can be given two different expressions, depending on whether the long wave length (LWA) or the Padé1 approximation is considered:

ρ pol,LW A (ψ, θ, φ, t) = ∇ ∇ ∇ ⊥ • (α∇ ∇ ∇ ⊥ ϕ) (E.30) ρ pol,P ade (ψ, θ, φ, t) = ∇ ∇ ∇ ⊥ • (α∇ ∇ ∇ ⊥ ϕ) 1 -∇ ∇ ∇ ⊥ • (κ∇ ∇ ∇ ⊥ ) (E.31)
with the α and κ functions defined as follows:

α(ψ, θ, t) = The hybrid model consists in solving the Vlasov equation for the entire electron distribution function, while retaining the kinetic contribution of the trapped ones in the QN equation, the passing ones being essentially given an adiabatic response. Such an adiabatic response is actually expected for passing electrons provided they have a low inertia. The immediate advantage of the hybrid approach is to allow one to consider (unrealistic) not-too-small electron to ion mass ratios, that potentially enable relatively large numerical time steps. Indeed, Trapped Electron Mode (TEM) turbulence develops due to wave-particle resonance at the precession frequency, which does not depend on the particle mass. Therefore, TEMs are expected to still be well described in this case. However, as originally highlighted by Idomura [?] and further investigated in [Lanti et al. 2020], special attention should be paid to axi-symmetric modes in general, and zonal modes in particular. These issues are discussed in the following.

The case of non-zonal modes (i.e. non flux-surface averaged, which we shall denote with a tilde . . .) is the most tricky. The response of passing electrons is adiabatic, namely: Remark 1: the adiabatic -or Boltzmann -response of electrons breaks down for flux-surface averaged fluctuations. Therefore, the electron density should not depend on the electric potential for these modes. This is indeed the case in the proposed expression eq.(E.34) because one only retains the flux-surface averages of n e,pass and T e (⟨n e,pass /T e ⟩ F S would fit as well, although less straightforward to justify). Conversely, if one replaces ⟨n e,pass ⟩ F S /⟨T e ⟩ F S by e.g. their toroidal averages ne,pass / Te , then one would find ⟨n e,pass ⟩ F S = ⟨n e,pass ⟩ F S + e ⟨n e,pass ϕ/ Te ⟩ F S -⟨n e,pass / Te ⟩ F S ⟨ϕ⟩ F S which is clearly inconsistent since the right hand side parenthesis is nonvanishing a priori.

Remark 2: notice that this constraint regarding adiabaticity would not prevent one from choosing an adiabatic response of the form: n e,pass = ne,pass + e⟨n e,pass ⟩ F S (ϕ n̸ =0 -⟨ϕ⟩ F S ) /⟨T e ⟩ F S . But, in this case, one may face spurious so-called ω H = (m i /m e ) 1/2 (k ∥ /k ⊥ )ω ci modes 2 that are high frequency and therefore require prohibitive small numerical time steps. In turn, unless addressing electromagnetic effects, the kinetic response of passing electrons to finite k ∥ modes should not be considered. In other words, the passing electron response to such modes has to be adiabatic. Remark 3: as discussed below, zonal modes are actually calculated from another equation than eq.(E.36). Therefore, one might be tempted to drop the zonal electric potential from the expression eq.(E.35). This is however incorrect and would lead to an erroneous estimate of non-zonal axisymmetric modes n = 0. Indeed, because of the poloidal dependency of α eq.(E.2), these modes are coupled to the zonal mode.

Let us now consider flux-surface averaged modes. In this case, the kinetic response of all electrons can be retained without any risk to excite ω H modes (since k ∥ = 0). Therefore, the QN eq.(E.26) reads for these modes:

-⟨ρ pol ⟩ F S = ⟨δρ Gi ⟩ F S -⟨δn e ⟩ F S (E.37)

Whatever the expression retained to compute ρ pol , eq.(E.30) or (E.31), equation (E.37) provides as an output the expression of the electric potential for all axisymmetric modes n = 0, not solely the one of ⟨ϕ⟩ F S . The reason is that, since α eq.(E.2) depends on the poloidal angle, it couples all poloidal mode numbers together, so that the determination of ⟨ϕ⟩ F S requires the knowledge of all axisymmetric modes ϕ m⋆̸ =0,n=0 . However, since the kinetic response of passing electrons should not be included to calculate these modes -apart from the zonal one for which k ∥ = 0 -to avoid ω H modes, these axisymmetric non-zonal modes will not be retained from equation (E.37). Equation (E.37) will therefore only provide ⟨ϕ⟩ F S , while all non-zonal modes will be given by eq.(E.36).

E.3.3.2 Numerical treatment

Solving eqs.(E.36-E.37) can be done in four steps:

1. Solve equation (E.37) -that involves the fully kinetic electron response (FKE) -for all axisymmetric modes, namely:

-ρpol = δ ρGi -δn e (E.38)

The right hand side is only 2D in space (ψ, θ). This gives access to all axisymmetric modes of the electric potential φF KE .

2. Retain only zonal modes ⟨ φF KE ⟩ F S . Indeed, since they have a finite parallel wave vector k ∥ , the other axisymmetric modes may be subject to high frequency ω H modes.

3. Solve eq.(E.36) with ⟨ϕ⟩ F S replaced by its value computed at the previous step: ⟨ φF KE ⟩ F S . This leads to:

e⟨n e,pass ⟩ F S ⟨T e ⟩ F S ϕ -ρ pol = δρ Gi -δn e,trap + n init e,pass -⟨n e,pass ⟩ F S 1 + e⟨ φF KE ⟩ F S ⟨T e ⟩ F S (E.39) with δ ρGi and δñ e,trap given by eqs.(E.27) and (E.28), respectively. This equation gives access to all the modes, both axisymmetric -including zonal ones -and nonaxisymmetric. They are computed assuming a kinetic response of trapped electrons only, hence its label ϕ T KE . 4. The final step consists in excluding zonal modes from this last equation. The final electric potential is then given as follows:

ϕ = ϕ T KE -⟨ϕ T KE ⟩ F S + ⟨ φF KE ⟩ F S (E.40)

Note that ⟨ϕ T KE ⟩ F S and ⟨ φF KE ⟩ F S are different whenever α is 2D and/or when the metric depends on θ. Indeed, in these cases, the flux surface averaged potential is coupled to non-zonal axisymmetric modes, which are computed with an adiabatic or a kinetic response of the passing electrons, respectively.

E.3.4 Kinetic treatment of the Scrape-Off layer

This section is under development, guided by experience from Voice and tailored to the specifics of Gysela. There are special problems (i) due to the gyrokinetic framework: the Poisson equation is written as ∼ ∇ 2 ⊥ ϕ, without knowledge of ∼ ∇ 2 ∥ ϕ at Debye scales which is the natural description for the sheath. Furthermore, (ii) the specifics of the Gysela numerical scheme (semi-Lagrangian and backward) provide additional constraints through global (spline) interpolation due to the importance of correctly dealing with what happens to the foot of characteristics that happen to fall within the limiter/wall/penalised volume.

E.4 Expression of the differential operators in an arbitrary geometry

While the general expression of the quasi-neutrality equation, Eq.(E.1), is unchanged when passing from a circular geometry to an arbitrary geometry, the differential operators (assuming e ∥ = e φ ) are now defined as: where g is the metric tensor for the coordinates of the poloidal cross-section (see the note on the Culham geometry for detailed expressions of the metric tensor). Using these expressions, it can be shown (note "magnetic Poisson" par Emily) that:

∇ ⊥ ϕ =
∇ ⊥ • (F ∇ ⊥ G) =

E.5 Numerical solver

For its main computation, the spline-based solver uses a linear algebra library implemented in Selalib (in ext/selalib-lite) which can solve 2D equations equations of the form

-∇ ⊥ • (F 1 ∇ ⊥ ϕ) -∇ ⊥ • (V ϕ) + F 2 ϕ = ρ (E.45)
where the F i are smooth scalar fields and V is a vector field on the poloidal plane.

The linear solver uses a finite element-like approach to solve the equation in weak form using a conjugate gradient method. Notice that the construction of the matrix in the library assumes that the domain includes the pole, requiring the use of rhomin=0 and skiphole=.true. in Gysela or the implementation of additional interpolations to add new points near the pole when converting from the Gysela grid to the one used by the solver.

E.5.1 Long wavelength approximation with an adiabatic response, fix point method

The quasi-neutrality equation in presence of adiabatic electrons contains an extra term compared with the selalib solver Eq. (E.45). To circumvent this issue while also making use of the linear solver, a fixed-point scheme combining 2 equations is used to compute φ = ϕdφ from which γ⟨ ϕ ⟩ FS = γ⟨ φ ⟩ FS is deduced and moved to the RHS. The equation for φ is obtained by computing the toroidal average of eq. (E.1)

-∇ ⊥ • α∇ ⊥ φ + β φ -γ⟨ φ ⟩ FS = ρ (E.46)
The condition of stop is ∥⟨ ϕ new ⟩ FS -⟨ ϕ old ⟩ FS ∥ ≤ tol • ∥⟨ ϕ new ⟩ FS ∥, where ϕ new and ϕ old correspond respectively to the latest solution and the previous solutions of solver B eq. (E.51).

Once the solver converged to φconv , one can compute -∇ ⊥ • κβγ∇ ⊥ ⟨ φconv ⟩ FS -∇ ⊥ • κ⟨ φconv ⟩ FS ∇ ⊥ (βγ) + βγ⟨ φconv ⟩ FS and therefore compute the RHS of Eq.(E.49) which is then used to compute the 3D electric potential.

E.6 Convergence of the current fixed-point method for the computation of the flux surface average

There is no formal proof for the convergence of either fixed-point methods in the general case. Their formulation is based on considerations in a simplified case, where convergence can be proven, and it is assumed that the spectrum of the differential operator remains close enough that convergence is not lost. In what follows, we will assume:

• circular configuration,

• α = 1 (which is the same as assuming it constant and renormalising),

• β is constant,

• ε ≪ 1, which turns the flux-surface average into a mere averaging over θ and φ.

E.6.1 Diagonalisation of the differential operator

Let's recall that ∇ 2 ⊥ = ∂ rr + 1 r ∂ r + 1 r 2 ∂ θθ . Furthermore, -∇ 2 ⊥ ϕ m = λϕ m , with:

ϕ ′′ m + 1 r ϕ ′ m + λ -
A c = -∇ ⊥ • (α∇ ⊥ ) + c φ -ξ φ + A -1 c β φ -c φ -βγ⟨ φ ⟩ FS -A -1 c ρ = φ (E.55)
where convergence is guaranteed if the application on left-hand side is a contraction i.e. that the eigenvalues of the linear part (without ρ) are smaller than 1 in absolute value. The equation above does not require the assumptions that we used before, however, making them allows calculating the -point-spectrum of that left-hand side: α log α p 3d (i, j, l) p 1d (j) p 2d xx (i, j) p 2d xy (j, l)

σ p = 1 -ξ 1 + β -c -βγδ m
(G.3)
where p 3d , p 2d xx , p 2d xy and p 1d are the discretised versions of respectively p(x n+1 , x n-k , y n-k ), p(x n+1 , x n-k ), p(x n-k , y n-k ) and p(x n-k ). In order to have sufficient statistics, a bin size β = 2 or β = 3 is typically chosen, depending on the available length of the time series (the longer the times series, the larger β can be). We introduced here the additional exponent α ≥ 1, which effectively represents a nonlinear threshold: low probabilities will be further reduced and higher ones amplified. In a complex setting, information may flow both ways, from Y to X and inversely. It is thus especially useful to define the net transfer entropy ∆ X,Y (T E)[k] = T E Y →X [k] -T E X→Y [k], which provides the net flow of information between processes X and Y , at timelag k. In the manuscript, pdfs in Eq.(G.2) are discretised using β d = 2 d bins, with d the dimensionality of the pdf. The nonlinear threshold exponent α is set to unity and X and Y are discretised at the same rate and enter the T E calculation with zero temporal mean. Further details may be found in (42).

We systematically apply the T E algorithm to actual time series from the flux-driven Case-1 computation with limiter boundary conditions in the last 5% inside the separatrix, where the spontaneous onset of a persistent transport barrier is observed. The following vorticity equation can be inferred from the primitive gyrokinetic equations including E × B drift and finite Larmor radius at leading order: where ⟨•⟩ denotes an average over toroidal angle φ. The T E algorithm is applied to many of the possible permutations of quantities in Eq.(G.4) and especially here to the following set:

∂ t ⟨Ω r ⟩ + 1 r ∂ r r ⟨v Er Ω r ⟩ + ⟨v ⋆r Ω r ⟩ -v ⋆θ 1 r ∂ θ E r = r
(X, Y ) ∈ ⟨Ω r ⟩ , ⟨v Er Ω r ⟩ , ⟨v ⋆r Ω r ⟩ , -v ⋆θ 1 r ∂ θ E r (G.9)

  Planetary Lorentz: ω c,i z × v Coriolis: f z × v Eq.(B.1): R ≡ ω c,i z = eB m Eq.(B.1): R ≡ f z = 2Ω sin θ z ion cyclo. freq. ω c,i Coriolis par. f = f 0 + βy, β-plane ion density n = n 0 + ñ atmosphere depth H = H 0 + h quasi-neutral: ñ n0 = eϕ T shallow-water: h H0 = o(1) Eq.(B.1): U ≡ eϕ m Eq.(B.1): U ≡ gH ′ f ′ ≡ elect. potential ϕ ′ f ′ ≡ depth h or streamfunction Ψ B field planetary rotation Ω small param: ρ ⋆ = ρ i /a Rossby # Ro = U/f 0 L Hasegawa-Mima (B.4) Charney (B.3) drift wave (westward) Rossby wave 2D plane ⊥ B 2D plane ⊥ (f z) advection: E × B advection: geostrophic vel. v E×B = -∇ϕ×z |B| v g = -g∇h×z f0 = -∇Ψ × z PV as in Eq.(B.4) ; in general, moments of PV as in Eq.(B.3) distribution function F α [6D] or (J 0 F ) α [5D]

Figure 2

 2 Figure 2.1: The E × B staircase, schematic view.

  2.2 as coherent flux-surface averaged E × B shear γ E×B = r∂ r (E r /rB). The flux-surface averaged poloidal ⟨v r v θ ⟩ and toroidal ⟨v r v φ ⟩ Reynolds stresses, the parallel momentum flux M ∥ and the poloidal v θ and toroidal flows v φ similarly display long-lived structuring close to the mean profile corrugations. Regarding poloidal rotation, the turbulence-borne temperature corrugations at the staircase steps neoclassically drive[START_REF] Helander | Collisional transport in magnetized plasmas[END_REF] localised poloidal rotation ∝ ∇T there, an example of synergistic interplay between collisions and turbulence [Dif-Pradalier et al. 2009]; Nota Bene: In Fig.2.3 are detailed the flow and stress patterns at the location of a temperature corrugation. This is obtained from a zoom of the data displayed in Fig.2.2 in the vicinity of the central staircase between ρ = 0.48 and ρ = 0.61. The profiles are averaged between t = 1860 a/c s and t = 2480 a/c s . As in the planetary staircase, the region of sharp temperature gradient [sharp PV gradient] coincides with the location ρ flow where a strong poloidal shear flow has nucleated and endures. Being sheared, this flow [the jet in the planetary case] is associated with a bipolar shear structure for γ E×B : the flanks of the flow are differentially rotating from its core. The poloidal flow profile shown in Fig.2.3 clearly reminds of the discussion presented in section B.1.2: strongly peaked with a definite sign (here positive) and confined to the region of steep temperature gradient [the analogous of the eastward mean zonal jet in the Earth case] it is flanked on either side by contra-rotating flows in what must be the dissipative regions. The downstream flank of this "plasma jet": ρ > ρ flow -the radial orientation is dictated by the central location of the heat source-corresponds to a region of intense mixing and intense Reynolds stress activity ⟨v r v θ ⟩. This is again reminiscent of the strong Rossby wave-mediated stresses in the flanks of the planetary jets that are key to the onset and sustainment of the planetary staircase.

Figure 2 . 3 :

 23 Figure 2.3: Detail of the flow-mean profile-stress interplay next to a staircase step [a corrugation]; averages are between 1860 a/c s and 2480 a/c s .

Figure 2 . 4 :

 24 Figure 2.4: Autocorrelation of the electrostatic potential. Two slopes are found: the local autocorrelation length ℓ c ∼ 6 ρ s is microscale, the mesoscale is indicative of avalanches and therefore noted ∆ (see section 2.3.2.2).

Figure 2 . 5 :

 25 Figure 2.5: Comparing zonal disturbance (black) and zonal mean shear (red) averaged radially from ρ = 0.4 to 0.8 in flux-and gradient-driven frameworks. Mean is defined through sliding time averages of 0.3 ms. For each case the steady-state zonal shear magnitude (disturbance plus mean, thick blue) is indicated [time averaged between 12 and 16 ms (left); between 17.9 and 18.5 ms (right)]. Corresponding levels of turbulence (δΦ 2 ) 1/2 are shown in Fig.C.3. The average gradient-driven zonal shear magnitude (resp. turbulence intensity) is 40% (resp. 15 -20%) larger than its corresponding flux-driven counterpart.

Figure 2 . 6 :

 26 Figure 2.6: The mean E × B shear, averaged over 0.3 ms (1 collision time at ρ = 0.5) is shown in both flux-and gradient-driven computations that display comparable turbulence intensities (see Fig.C.3). No perennial structures genuinely endure in the gradient-driven approach, mean components are vanishing and the staircase pattern is lost.

Figure 2 . 7 :

 27 Figure 2.7: Autocorrelation of the turbulent heat flux in two separate radial regions of the plasma, at the same time, during 100 a/c s : case (a) in-between two staircase steps distant of ∆ = 58 ρ s and case (b) distant of ∆ = 153 ρ s [shown as white interval]. In both cases the full width at half maximum (black contour), a measure of the radial coherence of heat transport, tracks well the local staircase step width.

Figure 2 . 8 :

 28 Figure 2.8: Staircase microbarrier spacing ∆ normalised to the local ion Larmor radius as a function (a) and (b) of gradient drive (distance to marginality) and of (c) plasma size ρ -1⋆ . The statistics of step sizes ∆ is shown in (d). The most probable microbarrier spacing is ∆ stat = 40 ± 2 ρ s . ∆ however can be large, displaying statistics of heavy-tailed distributions: a Fréchet distribution approximates it well, which is symptomatic of extreme value statistics.

Figure 2 . 9 :

 29 Figure 2.9: The Gysela computed probability density (the noiselike curve-shown in red (light gray) color) versus the Fréchet distribution (smooth humped curve) with κ ≃ 0.67 for which the optimum fit was obtained. The steplike curve (gray color) is coarse experimental distribution and interpolates the discrete distribution as of Fig. 4. The right panel summarises the normalised root-mean-square error in a percentile to the maximum error at κ ≃ 0.85, showing that the normalised error is minimised for κ ≃ 0.67.

Figure 2 .

 2 Figure 2.11: Radial profiles of normalised temperature gradients (blue and red, left axis) and zonal flow shear (green, right axis) for the 'above threshold' (top) and 'near marginal' (bottom) cases. Shaded areas represent temporal standard deviation; black hourglass symbols ( ) linear instability thresholds R/L lin T at vanishing E × B shear; red squares (■) nonlinear thresholds R/L T c estimated by Gkw in the local limit, including the Gysela E × B shear.

Figure 2 .

 2 Figure 2.12: Kubo numbers for the principal nonlinear dynamics in the problem. Plain and plusses: turbulent radial and poloidal E × B velocity effect during a turbulent autocorrelation time. Circles: trapping of particles due to turbulent vorticity during transverse crossing of the turbulent filament.

  3. With analogy to incompressible fluids, particles trapped in turbulent convective cells explore eddies in a typical turn-over time given by the local vorticity τ trap int ∼ B/⟨|∇ 2 ⊥ φ| 2 ⟩ 1/2 . Whilst they undergo this vortical motion they also drift in about τ ⋆ jump ∼ L θ (eB/∇T ) from one turbulent structure to the next at the typical speed of the local diamagnetic velocity. Particle trapping K trap transverse drifts τ ⋆ jump eddy turn-over τ trap int Random walks K {r,θ} diff Lagrang. correlation time τ corr jump E × B random walk τ diff,{r,θ} intTable2.3: Five typical wave-particle and turbulent times lead to three Kubo number combinations, shown in Fig.2.12.

Figure 2 .

 2 Figure 2.13: (a) and (b): gradient-driven QuaLiKiz and Gkw heat fluxes confronted to reference flux-driven Gysela flux levels, based on the Gysela profiles of Fig.2.11 and expressed in gyro-Bohm unit. Gray shaded areas represent a standard deviation of Gysela heat fluxes during the considered time interval, as profiles fluctuate.Red and blue shaded areas represent the sensitivity of QuaLiKiz to these profile variations; the hourglass symbols that of Gkw to a increased input temperature gradients. The reversed approach is followed in panel (c): heat fluxes in QuaLiKiz-Jetto and Gene-Tango are made to match the Gysela reference fluxes; the unknowns are thus the QuaLiKiz-Jetto and Gene-Tango profiles. The remarkable agreement above threshold and large over-prediction of the temperature gradient near marginality are consistent with results in panels (a) and (b). This emphasises model reduction adequacy in 'above threshold' regimes and missing physical ingredients in 'near marginal' regimes.

Figure 2 .

 2 Figure 2.14: Testing spectral distribution of key fundamental quantities: (a) linear crossphases, (c) heat flux, (d) saturation rule and (b) their ratio, proxy to the dispersion relation from Gysela (---), Gkw (×) and QuaLiKiz (•), averaged over radial interval 0.3 ⩽ r/a ⩽ 0.7.

  /2 + P * ⊥,n dθ, where parallel and transverse components of the pressure and the radial component of the E × B drift are computed from 3D output data. The n-Fourier components of the heat flux |Q n | and of the squared potential |ϕ n | 2 , i.e. the saturation rule are respectively shown in panels (c) and (d). Panel (b) displays their ratio, a proxy for the dispersion relation, sometimes called 'QL flux integrals'. Clearly, linear cross-phases display reasonable agreement, within factors of 2 and across all regimes. They are not responsible for the flux discrepancies. Factors of disagreement -and avenues of improvement for QL modeling-are essentially twofold: [i] 90% of the heat flux is carried by modes n ∈ [[5, 50]] in the above threshold regime; all approaches provide similar conclusions. Near marginality however, flux is carried in Gysela and Gkw through n ∈ [[5, 70]], twice the amount of active modes with respect to QuaLiKiz. Furthermore, [ii] saturation rules are clearly responsible for much of the observed flux discrepancies. In the above threshold regime, flux spectra agree well [panel (c)] -though this results from a surprising compensation: the potential, or saturation rule [panel (d)] is under-estimated, the dispersion relation [panel (b)] is over-estimated yet heat flux spectra in the above threshold regime are in reasonable agreement. No such compensation occurs near marginality; there severe under-estimation of the potential spectrum is clearly responsible for the under-prediction of fluxes.

Figure 2

 2 Figure 2.15: Staircase: experimental evidence. Successive staircase shear layers appear as quasi-regularly spaced local contractions of the coherence (a) spectrally averaged over turbulent frequencies, (b) displayed as a radial profile (gray: 95% confidence interval). Signatures of shear flows at the staircase steps are shown through (c) the evolution of the asymmetry of the cross-correlation functions along the radial profile and (d)-(e) differential tilting on each side of the staircase step.

Figure 2 .

 2 Figure 2.16: Localisation versus external heating. (a) pressure profile modification during additional heating; (b) corresponding inward drifts of local minima of the mean coherence length.

Figure 2 .

 2 Figure 2.18: Staircase and LOC-SOC transition. (a) Energy confinement time against mean electron density: staircases observed in SOC only. LOC (b) and SOC (c) frequency power spectra at r/a ∼ 0.66 and their respective coherence lengths.

Figure 2 .

 2 Figure 2.19: Time evolution of the driving temperature gradient and of the responding heat flux near the source region at ρ = 0.3. The flux responds with a time delay τ ≈ 30 a/c.

Figure 2 .

 2 Figure 2.20: Illustrating (a) stable/unstable (dotted line, as in Ref.[Posmentier 1977]) or stable/unstable/stable (solid line, as in Ref.[Balmforth et al. 1998]) flux-gradient relation, as per Eq.(2.10). On the right panel (b), an illustration from Gysela of the instantaneous heat flux-temperature gradient relation, in the near marginal regime where staircasing is observed (yet unpublished figure).

Figure 3

 3 Figure 3.2: An attempt to address the shortfall problem, from Ref.[Görler et al. 2014].

  2):Case-1 is shown in Fig.3.1 and refers to flux-driven evolution with limiter. It allows description of the transition from closed to open field lines. In this case, D = ν S = γ GD = 0. The distribution function evolves according to volumetric sources S [Sarazin et al. 2010]. The remaining BGK operator in Eq.(A.

"

  Le vacarme intermittent du petit coin [de limiteur] me rassure." Paul Valéry, variations on Blaise Pascal ("Le silence éternel des espaces infinis m'effraie.")

Figure 3 . 4 :

 34 Figure 3.4: Spatiotemporal redistribution of induced free energy injection at the plasma-boundary interface (flux-driven Case-1). Subplots (a)-(g): times series of poloidal cross-sections of spreading increments ∆S. Positive (red) increments represent radially-outward fluxes of turbulence intensity; inward fluxes of turbulence intensity are pictured as negative (blue) increments. The limiter is responsible for vigorous mixing in its immediate poloidal (between θ 3 and θ 5 ) and radial (0.95 ≤ r/a ≤ 1.0) vicinity. The drive endures throughout linear -subplot (a)-and nonlinear evolution [(b) through (g)]. The choice of an early reference time t ref = 30, 000 Ω -1 ci allows to follow the full unfolding of the spreading sequence. With a later choice for t ref the early spreading sequence would not appear as clearly; later nonlinear dynamics (Phases II and III) on the other hand would be qualitatively unchanged: cyclic equilibration on the outboard midplane between incoming and outgoing fluxes of turbulence intensity that originate from distinct poloidal regions. Inward spreading (Phase I) is responsible for edge increase of turbulence activity at the outboard midplane [subplot (h)]. As NMsL fills-in, outgoing turbulence activity spreads to the edge [Phase II, subplot (i)], further enhancing relative fluctuation levels in the outer edge. At equilibrium (Phase III), the fluctuation profile hovers about its 1.84ms value, equilibrating incoming and outgoing spreading increments.

Figure 3 . 5 :

 35 Figure 3.5: How forcing and interplay with boundaries impact redistribution of turbulence activity and globally affect transport. Spreading increments for Case-2 (a) and Case-3 (b) at time t ref = 50, 000 Ω -1 ci when both Cases enter nonlinear regime in the edge. A clear shortfall in NMsL is visible at statistical equilibrium (subplot (c)-top) for Cases 2 and 3 with respect to experimental or to Case-1 fluctuation profiles δn/n. Differences ∆ sp [1] and ∆ sp [2] respectively quantify (subplot (c)-bottom) the importance for turbulent fluctuations of separatrix→core and core→edge spreading.

  [START_REF] Hlavackova-Schindler | Causality detection based on information-theoretic approaches in time series analysis[END_REF]].

  Figure 3.6: Causal interactions discriminated with proximity to boundaries in the early stages of barrier build-up (initial 0.3ms, or 50, 000 Ω -1 ci ), between 0.95 ≤ r/a ≤ 1.0 where limiter-induced fluctuations are borne.

  3) T e e ∆n e = ϕ -λ 1 -M SOL (r) ⟨ ϕ ⟩ FS -M mat (r, θ) -M wall (r) ϕ bias -λΛ M SOL (r) -M mat (r, θ) (T e -T b.c.

Figure B. 1 :

 1 Figure B.1: (Colour online) The GFD and MFE geometries and their system of coordinates.

  the vorticity of the ion fluid. Quasineutrality is assumed so that the electron distribution obeys a Boltzmann distribution: n i ≈ n e = n 0 (r) exp(eϕ/T e ) where n i [n e ] is the ion [electron] density, ϕ the electrostatic potential such that eϕ/T e = o(1) and T e the electron temperature. As in the GFD case, the ion density can be expressed perturbatively: log n i ≈ log n 0 + eϕ/T e . At leading order, whilst neglecting the advective terms in Eq.(B.1), the dominant velocities are respectively the geostrophic v g = -g∇h×z f0 and the E × B velocities v E×B = -∇ϕ×z |B| : horizontal advection in GFD is accomplished by the geostrophic winds; it is accomplished by the E × B drift in MFE. Coupled to the continuity equation:

Figure

  Figure B.2: (Colour online) The staircase formation mechanism in the GFD case; the flow and PV profiles are remarkably close to our analogue plasma quantities, see Fig.2.3.

Figure C. 1 :

 1 Figure C.1: (Colour online) Typical heat flux-driven setup.

Figure C. 2 :

 2 Figure C.2: (Colour online) Evolution, as the temperature fluctuates, of the adaptive gradient-driven source/sink Eq.(C.2) from one time step t 0 to a next t 0 + ∆t. The gradientdriven forced profile (solid line) is computed from the flux-driven mean T FD profile in Fig.C.1 and coincides with the flux-driven 'ambient mean profile', as defined on page. For clarity, only a zoom between ρ = 0.37 and 0.57 is shown.

Figure C. 3 :

 3 Figure C.3: (Colour online) The "spring constant" γ K in Eq.(C.2) [here γ K = 5.43 10 -5 ω c,i = γ lin /10] is chosen such that it agrees with usual prescriptions for gradientdriven models γ K ∼ γ lin /10. With this prescription, computed turbulence intensities (δΦ 2 ) 1/2 with Gysela in flux-and gradient-driven approaches agree within 15 -20%.

  Fig.C.3. γ K has been varied between γ lin /3 and γ lin /15. The results in the paper are qualitatively independent of the precise choice of γ K . Unless otherwise stated: γ K = γ lin /10, consistently with [McMillan et al. 2008];

e

  (ψ, θ, φ, t) -δ⟨ n adia. e ⟩ FS (ψ, t) (E.4)For any quantity X, δX = X(t) -X(t = 0). Note that the partition between adiabatic and kinetic electrons can be different between the core and the Scrape-Off layer depending on the model chosen to describe this region of the plasma, see section E.3. The guiding center density is defined asn Gi (ψ, θ, φ, t) = µmax 0 dµ v ∥,max -v ∥,max dv ∥ J v J 0,i • F i ψ, θ, φ, v ∥ , µ, t (E.5)and where n kin. e is the density of kinetic electrons. The different masks are represented in fig.E.1 in the case of a circular cross section.

Figure E. 1 :

 1 Figure E.1: Improved limiter masks used whilst solving Poisson with SOL and Flat or Rounded Limiter.

e

  and n adia.

  e

E

  .3.1.0.1 Adiabatic limiter The first possible treatment constists in treating electrons whose trajectories intercept the limiter as adiabatic. The other electrons that do not intercept the limiter are treated kinetically. The implementation of this condition is straightforward and consists in changing the maximum magnetic field B max in the trapping condition:v ± ∥ (ψ, θ, µ) = ± 2µ [B LIM max (ψ) -B (ψ, θ)] (E.11) with B LIM max (ψ) = 1 -M SOL (r) max θ [B(ψ, θ)] + M SOL (r) min θ B LIM (ψ, θ) (E.12)

e

  (ψ, θ, φ, t) -δ⟨ n adia. e ⟩ FS (ψ, t) (E.13)+β (ψ, t) Λ M SOL -M mat ⟨ T e ⟩ FS -T b.c. e + M mat -M wall ϕ bias (1.0.2 Kinetic SOLThe other solution reads as follows.

  v F e ψ, θ, φ, v ∥s , µ s , t (E.17)

  e

  dv ∥ J v F e ψ, θ, φ, v ∥s , µ s , t (E.19)

ϕ

  tot. = ϕ TKE + -ϕ TKE (m⋆=0,n=0) + ϕ FKE (m⋆=0,n=0) (E.20)Note that for the fully kinetic electron part, only the axisymmetric part of the potential φ needs to be computed using Eq.(E.46). The last difficulty is to extract the component {n = 0, m ⋆ = 0} of a 3D field g (ψ, θ, φ

  r) is the Shafranov shift, E (r) is related to the elongation κ (r) ≃ (r + E (r)) / (r -E (r)) and T (r) is related to the triangularity δ (r) ≃ 4T (r) /r. With the Culham equilibrium, one hasdθ ⋆ dθ = 1 -∆ ′ + r R 0 cos (θ) -E ′hybrid version -April 21, 2023We start from the general expression of the QN equation:δρ Gi + ρ pol -δn e = 0 (E.26)with the following definitions for the ion and electron (δn e = δn e,trap + δn e,pass ) charge densities:δρ Gi (ψ, θ, φ, t) = i Z i d 3 v J 0,i .[F i (ψ, θ, φ, v ∥ , µ, t) -F init i (ψ, θ, v ∥ , µ)] (E.27) δn e,trap (ψ, θ, φ, t) = trap d 3 v [F e (ψ, θ, φ, v ∥ , µ, t) -F init e (ψ, θ, v ∥ , µ)] (E.28) δn e,pass (ψ, θ, φ, t) = pass d 3 v [F e (ψ, θ, φ, v ∥ , µ, t) -F init e (ψ, θ, v ∥ , µ)] (E.29)Note that the initial distribution functions have to fulfill the QN condition, namelyi Z i d 3 vJ 0i .F init i = d 3 vF init e .The velocity space element is given by d 3 v = J v dµ dv ∥ with J v = 2πB * ∥s /m s the jacobian in the velocity space. The trapped and passing integrals refer to the trapped and passing domains defined by trap , toroidally averaged quantities will be denoted with a bar: Ȳ = Lφ 0

  of passing/trapped electrons and of zonal/non-zonal modes

  ñe,pass = ⟨n e,pass ⟩ F S 1 + e ⟨T e ⟩ F S (ϕ -⟨ϕ⟩ F S ) = n adiab. pass = ⟨n e,pass ⟩ F S -n init e,pass + e⟨n e,pass ⟩ F S ⟨T e ⟩ F S (ϕ -⟨ϕ⟩ F S ) = δn adiab. e,pass (E.35) The only poloidal dependence of axisymmetric modes is carried by the electric potential. The QN equation reads for these modes: δñ e,pass -ρpol = δ ρGi -δñ e,trap (E.36) Beware that the notation ρpol refers here to ∇ ∇ ∇ ⊥ • (α∇ ∇ ∇ ⊥ ϕ) (in the long wave length approximation), and not to ∇ ∇ ∇ ⊥ • (α∇ ∇ ∇ ⊥ φ).

  m is the p-th zero of J m .E.6.2 Eigenvalues in the fixed-point formφ -ξ -∇ ⊥ • α∇ ⊥ φ + β φ -γ⟨ φ ⟩ FS -ρ = φ (E.54)

  Figure F.2: Data from flux-driven poloidally-symmetric Case-2, at equilibrium (tΩ ci = 200, 000). (a) Poloidal cross-section snapshot of the fluctuating electric potential. Specific locations are marked, combination of three radial locations: {r 1 , r 2 , r 3 }/a = {0.90 (circles), 0.96 (stars), 1.02 (squares)} and four poloidal locations θ k = {9 • (magenta), 126 • (cyan), -118 • (red), -61 • (yellow)}. Properties at these locations are shown in subplots (b): E × B shear, (c): parallel velocity shear and (d) through (f): maximul linear instability growth rate at vanishing E × B shear, keeping the same symbol-color combination.

Figure F. 3 :

 3 Figure F.3: Same layout as in Fig.F.2. Data is from flux-driven Case-1 with limiter at two different times: subplots (a) through (f) are early in the nonlinear development of turbulence (tΩ ci = 30, 000); subplots (g) through (l) at thermal equilibrium (tΩ ci = 250, 000). An additional location near the limiter at (r/a = 0.96, θ = -75 • ), marked by the large white triangle is shown and corresponds to the region of maximum linear instability growth.

v

  .h.s (G.4) r.h.s ≈ -∂ t ⟨Ω θ ⟩ -1 r ∂ θ ⟨(v Eθ + v ⋆θ )Ω θ ⟩ -Er = -∂ θ ϕ/r & v Eθ = ∂ r ϕ = -E r (G.7) v ⋆r = -∂ θ p ⊥ /r & v ⋆θ = ∂ r p ⊥ (G.8)

  ].

	Relevant co-authorship:
	[Dif-Pradalier et al. 2009]
	[Dif-Pradalier et al. 2010]
	[Kosuga et al. 2014]
	[Ghendrih et al. 2014b]
	[Ghendrih et al. 2014a]
	[Ghendrih et al. 2014c]
	[Cartier-Michaud et al. 2014]
	[Norscini et al. 2014]
	[Dif-Pradalier et al. 2015]
	[Ida et al. 2015]
	[Dif-Pradalier et al. 2017b]
	[Hornung et al. 2017]
	[Milovanov et al. 2021a]
	[Garbet et al. 2021]
	[Gillot et al. 2021b]
	[Sarazin et al. 2021]
	[Yamamoto et al. 2023]
	[

  E.6) where Fi ψ, θ, v ∥,i , µ i , t = ni (ψ, θ, t)

							√	2πv T i	3 exp -	v 2 ∥ + 2µiB(ψ,θ) mi T i 2v 2	(E.7)
	with vT i (ψ, θ, t) =		Ti(ψ,θ,t) mi	. Then compute
	η kin. e	(ψ, θ, t) = 2	0	µmax	dµ	v ∥,max ∥,e (ψ,θ,µ) v +	dv ∥ J v fe ψ, θ, v ∥ , µ, t	(E.8)
	where	fe ψ, θ, v ∥ , µ, t =	√	1 2πv T e	3 exp -	v 2 ∥ + 2µB(ψ,θ) me 2v 2 T e

  ∂ϕ ∂x j g jk e k =

				∂ϕ ∂r	g rr +	∂ϕ ∂θ	g θr e r +	∂ϕ ∂r	g rθ +	∂ϕ ∂θ	g θθ e θ	(E.41)
	and
	∇ • X ⊥ =	1 |det g ab |	∂ ∂x j	|det g ab |X j			(E.42)
	=	1 |det g ab |	∂ ∂r	|det g ab |X r +	1 |det g ab |	∂ ∂θ	|det g ab |X θ	(E.43)

This designation is disputed, on the grounds that transition from gas to plasma occurs continuously and not through a phase transition.

Integrated modelling frameworks, while flux-driven in the sense that a stationary state is reached at power-balance, still implicitly make the scale-separation assumption, since the transport PDEs (macro model) call the reduced turbulence model (micro model) for an instantaneous update following a macro timestep, typically on order of 1ms.

The Padé approximation is given a symbolic expression in eq.(E.31). Details of its numerical implementation will be given in a following section

These modes are understood as the electrostatic limit of kinetic Alfvén waves[Lee 1987, Scott 1997]. They disappear (actually they transform into kinetic Alfvén waves characterized by much smaller frequencies) as long as β ≳ (k ⊥ ρs) 2 (me/m i ), condition that is easily fulfilled in tokamak plasmas.

See https://mathworld.wolfram.com/BesselFunctionZeros.html for the values of the first few zeros of J 0 through J 5

Appendix A. Model equations with adiabatic electrons

Figure A.1: The various masks used for penalisation in the gyrokineticquasineutrality system. The resolved domain spans from the very core r/a = 0 to the outer wall region at r/a = r max = 1.3 (a). The outer wall is circular and within 1.25 ≤ r/a ≤ 1.3 (b). The minimum radius where mask M mat (r, θ) = 1 2 defines the location of the magnetic separatrix at r/a = r separatrix = 1. An appealing aspect of penalisation is the ease with which the shape of masks (i.e. the geometry of the material boundaries) can be altered. Two examples are shown in panels (c) and (d); all limiter computations shown here are performed in geometry (d).

described in Ref. [Dif-Pradalier et al. 2017b] (see Appendix C) and built such as to prevent overdamping zonal modes [McMillan et al. 2008].

Penalisation [Paredes et al. 2014] modifies the equations through introduction of a series of masks M mat (r, θ), M SOL (r) and M wall (r), combinations of hyperbolic tangents, adjustable in location, shape (for M mat (r, θ)) and stiffness. They are illustrated in Fig. A.1. Electrons have a Boltzmann response modified by penalisation such that the electric potential ϕ in the quasi-neutrality equation is relaxed towards its expected presheath condition ΛT e /e in the SOL. For an extension with trapped kinetic electrons in regions where magnetic field lines do not intercept the limiter, see section 3.4 and Appendix E. Additionally, ϕ may be biased to ϕ bias in the limiter (ϕ bias = 0 in the current study) and is freely evolving elsewhere. T b.c. e is the cold electron temperature within limiter and wall, chosen as the minimum T e value within the computational domain, Λ = -1 2 ln 2π me mi 1 + Ti Te and coefficient λ (set to unity in the present study) may be used to alter the inertia of the zonal potential. In the gyrokinetic equation, infinite penalisation [Caschera et al. 2018] relaxes Fs to a target cold Maxwellian distribution function G cold = n w (2πT w ) -3/2 exp[-(v 2 G∥ + µB)/2T w ], characterised by low wall thermal energy T w and target density n w . The former is constrained by velocity-space resolution; we typically choose it an order of magnitude lower that temperature at mid radius whilst the target density n w is chosen so as to ensure particle conservation. Articifial diffusion D can be added but will not be considered here. Collisions C are discussed in Ref. [Donnel et al. 2019].

B.1 Formal analogies between fusion plasmas & planetary fluids

There are clear analogies between aspects of astrophysical fluid dynamics (AFD) or geophysical fluid dynamics (GFD) and magnetised fusion plasmas (MFE). In order to illustrate this point, we briefly discuss two simple yet illustrative systems: one is the celebrated systematic scaling theory for large-scale atmospheric motions by Charney and the derivation of the quasi-geostrophic equations [Vallis 2006], the other is the well-known Hasegawa-Mima [HM] model [START_REF] Hasegawa & Mima | Pseudo-three-dimensional turbulence in magnetized nonuniform plasma[END_REF] of MFE. The more general Hasegawa-Wakatani [HW] framework in MFE or the generalised Charney equations in GFD are often preferred due to the possibility that they offer to describe either drift wave or interchange instabilities (i.e. particle transport) in MFE or their atmosphere-ocean counterpart, the barotropic instability. Anyhow, the HM model (or Charney-Hawegawa-Mima) captures the essential physics of Drift-Wave (or Rossby Wave) Zonal Flow formation while involving the minimal number of scalar fields relevant [Diamond et al. 2008]. Let us venture the following 2 remarks:

• as will be clear below, the common structure of both models in these 2 very different systems underlies deep commonalities between GFD and MFE which stresses that e.g. the physics of jet sustainment and patterning and the secondary establishment of staircases should thus be expected as well in MFE. The mere simplicity of the above models to illustrate some of the basic ideas that robustly manifest in such disparate frameworks is certainly appealing. Their importance in GFD to weather forecasting is a possible indication of the importance of zonal flows and staircases in MFE as well;

• however, beyond clear commonalities, important differences exist between both systems (regarding instability mechanisms, the possibility for fronts of avalanches, the existence of an inertial range, boundary conditions, . . . ) so that one should remain weary of straightforward conclusions based on the sole knowledge of the other system.

PV as in Eq.(B.3)

Table B.1: The planetary-tokamak analogies summarised. framework towards a kinetic description of the magnetised plasma an infinite number of material invariants exist, the more obvious being either the 6D kinetic distribution function F itself or in the gyrokinetic approach the 5D gyroaveraged distribution J 0 F later introduced as f = Fs in Eq.(A.1), J 0 being the gyroaverage operator. The key importance of a material invariant for staircase formation is discussed below.

Focussing on the similar structure of Eqs.(B.3) and (B.4), regardless of the actual variables: either height and the Coriolis parameter or density, electrostatic potential and magnetic field, the equations above are often referred to as the Charney-Hasegawa-Mima (CHM) equation. In the presence of inhomogeneities, the CHM predicts the existence of a linear wave that will admit in both GFD and MFE the same structure for its dispersion relation and thus similar propagating properties. Waves in fluids result from the action of restoring forces on fluid parcels displaced from their equilibrium positions. In the GFD case the dominant restoring force comes from rotation and the associated wave is the well-known westward Rossby wave; in plasmas rotation or gravity are negligible and the dominant effects are either electromagnetic (Alfvèn waves are ignored here) or come from compressibility. The so-called drift waves are akin to sound waves and are the MFE analogue of the Rossby waves.

At last, one of the key points for staircase formation that shares common grounds between both systems concerns the forcing. In the planetary case, as the Sun's energy reaches the Earth unevenly, the heat source budget is that of a hot equator (θ = 0) and a cold pole (θ = ±π/2). The wind pattern and the various gradients freely adjust according to this quasi-steady source distribution and to the planetary rotation. In a tokamak, the core (r = 0) is hot and the edge (r = a) is cold. A quasi-steady centrally located heat source imposed from the outside similarly exists which in combination with the given magnetic geometry leads one to naturally expect self-adjusting flow and stress patterns maybe not so different from those observed in the GFD case as well as patterns in the pressure gradient distribution. One of these patterns is of course the staircase that appears as a possibly

Appendix C

Flux-driven (FD) versus gradient-driven(-like) (GD) forcing

C.1 Same code, two frameworks

The usual flux-gradient paradigm [Dimits et al. 2000] relates local values for heat fluxes to local values of the gradients η-or R/L T in the case of ITG with low to moderate L n , with fixed mean values for the turbulence drive, i.e. fixed mean profiles. This approach is referred to as "fixed-gradient" or "gradient-driven" (GD) as opposed to "flux-driven" (FD) where both fluxes and profiles evolve consistently [Garbet et al. 1999] and are unknown functions of the dynamics. Formally, the flux-driven gyrokinetic equation solved in Gysela (explicit in Appendix A, Eqs.(A.1)-(A.2)) can be rewritten as: Switching between a flux-and a gradient-driven framework is done whilst replacing the source term in Eq.(C.1) by a Krook-type operator that acts as both a source or a sink, depending on the dynamics (the sign of f -F FD ) and driving a restoring force with strength γ K towards the chosen mean profiles F FD :

Infinite penalisation versus large BGK restoring force

This appendix is detailedly explained in [Caschera et al. 2018] and reproduced here for completeness. We detail here the difference in the numerical application of an infinite penalization as described in Section 3.2.1.2 and a very strong restoring force of the Krook type.

The numerical algorithm of Gysela [Grandgirard et al. 2016] includes a time-splitting technique which allows to solve the time advance due to each operator of Vlasov equation separately. Thus the equation for the Krook term corresponds to

To numerically handle ν → ∞ one possible choice is to apply a standard Krook term with a very large ν. However for increasing ν the shape of the coefficient ν(x) = νM (x) is close to a step function, which can lead to numerical instabilities. In this case, one should first define the maximum gradient numerically acceptable for the function ν(x) and next, given the restoring force strength ν, recover a modified mask function M (x) = ν(x)/ν to use in the restoring force term. Instead, the following alternative approach can be pursued: at each time step ∆t, inside the limiter and wall region, the distribution function is fully replaced by its target value g as

This formula differs from the application of a standard Krook term with very large but finite ν only in the transition region, where 0 < M (x) < 1. We show the proof Hereafter.

The equation for the Krook term solved by Gysela corresponds to

Where we substitute νM (x) = ν(x), allowing a spatial dependence of the restoring force strength which does not depend directly on the mask function. The analytical solution of this equation over a time step of amplitude ∆t, corresponding to the time step of one iteration, is

After k iterations the solution becomes

On the other hand, when applying an infinite penalization as in equation D.2, the spatial dependence of the restoring force is given by the mask M (x) only. After k iterations

The results D.5 and D.6 are equivalent if

The right hand side of Equation D.7 does not depend on time step ∆t or restoring force strength ν, as opposed to the left hand side. Thus in the case of infinite penalization as in equation D.2 the effective strength of the restoring force is ν(x) = log(1 -M (x))/∆t and depends on the time step ∆t. Substituting ν(x) = ν M (r, θ) in the relation D.7 we distinguish three different cases: (i) the case M (x) = 0 corresponds to ν = 0 and the restoring force is not active in both the formulations; (ii) in the case M (x) = 1 the relation D.7 exactly gives the limit ν → ∞; (iii) the case 0 < M (x) < 1, and so 0 < ν(x) < ν, is restricted to a narrow region ∼ 5ρ i and corresponds to the only difference between the two methods. Hence, the formulation D.2 can be preferred to the application of a very strong Krook term when one would enforce the limit ν → ∞, taking into account that the effective restoring force strength depends on the time step only in the narrow transition region.

Appendix E

Normalized quasi-neutrality (QN) equations for various (adiabatic, kinetic) electron responses, with and without limiter penalisation "Moi, j'ai appris à lire. Et ben, je souhaite ça à personne !"

Alexandre Astier/Leodagan

This appendix is part of a technical note detailing the various strategies actually implemented or considered for implementation in Gysela. This is a collective and still under construction endeavour; the especially difficult case of a fully kinetic electronic response in the SOL, with a limiter boundary is still under active investigation [Munschy et al. 2023b, Munschy et al. 2023a] as of Sept. 2023. This point is briefly touched upon in section 3.4.

E.1 Long wavelength approximation of the quasineutrality equation

We explicitly write the quasi-neutrality in the most general form to simplify the comparison between models and their implementations. Compared with the previous quasi-neutrality solver in GYSELA, there are three novelties: (i) the geometry can be arbitrary (but axisymmetric), (ii) the spatial dependence of the magnetic field is taken into account in the polarization term, (iii) densities are allowed to evolve in time and to have a poloidal dependence.

Whatever the model chosen for electrons (adiabatic, kinetic or hybrid), the quasineutrality equation, in the long wavelength approximation is of the generic form

The first term on the left hand side corresponds to the ion density polarization. It is in all cases defined as

where all quantities with a bar are toroidally averaged. Here the ψ, θ coordinates are used to explicitly emphasize the fact that one deals with 2D spatial functions. But θ actually refers to any periodic coordinate in the poloidal plane. Also, we use ψ as a flux surface label but any other label is possible.

The two equations that are solved for the fixed-point method are

which is efficient at converging low radial modes, and

which is efficient for high ones.

The selalib solver can be used to solve both eq.(E.47) and eq.(E.48). The procedure is then: Note that all the solvers are seeded with the solution of their last call as initial guess, allowing their conjugate gradient to convergence much faster.

E.5.2 Padé approximation with an adiabatic response, fix point method

The method is similar to the one developed for the LWA case but with more terms. It can be shown that Eq.(E.10) is equivalent to:

The general solver Eq.(E.45) can therefore be used at the condition to be able to compute

To do so, a fixed point approach, similar approach to the one used in the LWA case is applied. The two equations solved are:

The selalib solver can be used to solve Eqs.(E.50) and eq.(E.51). The procedure becomes:

1. Initialise by solving where the circular configuration and ε ≪ 1 also allowed us permuting A c and the fluxsurface average, here represented as a dirac.

It is then useful to distinguish "low radial modes" for which Z m p /a 2 ≪ 1 (which is true for the first few zeros of any m since a ≪ 1 3 ) and "high radial modes" for which Z m p /a 2 ≫ 1 (for l points per ion Larmor radius, the highest modes are roughly Z m p /a 2 ∼ l 2 π 2 ).

E.6.2.1 Getting high radial modes to converge

Asymptotically, the eigenvalues for those modes read

i.e. the modes converge as long as 0 ≤ ξ ≤ 2, with the central value ξ = 1 being optimal (independently of c unless very large values are used).

E.6.2.2 Getting low radial modes to converge

When Z m p /a 2 ≫ 1, it is important to use c ≫ 1/a 2 to prevent the denominator from creating very large eigenvalues as compensating via a very small ξ would prevent any convergence in other cases. However, for purely adiabatic electrons, where γ = 1, the mode m = 0 behaves differently as c ≫ 1/a 2 leads to eigenvalues of the form 1 -ξ Z m p /a 2 /c ∼ 1 and thus very poor convergence. Instead, c = 0 is ideal for this particular case (again, with ξ = 1), despite making the other low modes modes diverge. Although it does not seem possible to separate the modes in the real equation, where the Jacobian also depends on θ (leading to a more complex interaction between the flux-surface averaging and the double gradient), filtering can be applied to φ while iterating the fixed-point to take advantage of this convergence without its side-effects.

E.6.3 Iterating the fixed-point schemes

Combining the previous considerations, we expect the first fixed-point scheme for c = β and ξ = 1: φn+1

to converge to the right solutions, although the low radial modes of m = 0 would do so extremely slowly when γ = 1. Conversely, the second fixed-point scheme for c = 0 and

is expected to recover such modes while causing the divergence of all the other ones. The final procedure is thus based on a call to the latter scheme every few iterations of the former one. The reference Tore Supra shot 45511 had 2MW of Ion Cyclotron Resonance Heating on top of 1MW of Ohmic heating injected in a deuterium plasma of normalised size ρ ⋆ = ρ i /a = 1/500 at mid radius and aspect ratio a/R 0 = 1/3.3, a and R 0 being respectively the minor and major radius. The plasma current is I p = 0.8MA, the magnetic field on axis is B 0 = 2.8T and the density and temperature at mid-radius respectively read: n = 4 10 19 m -3 and T = 0.8keV. In GYSELA, a 3MW volumetric heat source comparable in shape to that in the experiment is injected in the central 40% of a deuterium torus of same aspect ratio. Initial density, electron and ion temperature profiles are the same as in the experiment up to the separatrix. In the core T e /T i > 1 whilst this ratio reverses in the edge and SOL. To slightly offset the numerical cost of the computations, run on Tier-0 Joliot-Curie at GENCI@CEA and MareNostrum at Barcelona Supercomputing Center, we assume a reduced magnetic field on axis: B 0 = 1.7T, which amounts to computing a plasma column of slightly smaller size ρ ⋆ = 1/300 on a 1/4 wedge torus with (r, θ, φ, v ∥ , µ) = (512, 1024, 64, 128, 64) grid.

Core-Edge-SOL appendices

Extensive tests have been performed to assess the robustness of the reported main results: observed gradient anisotropy and magnitudes in the limited configuration are robust whilst varying distribution function initialisation (local or canonical Maxwellians), presheath values in the SOL through Λ scans (from 0 to 5, its nominal value for Deuterium being Λ ∼ 4.1), target penalised temperature T w in the limiter and wall, limiter shape (large and flat, narrow, rounded) and poloidal location (bottom, top). Further scans have also been performed in a flux-driven poloidally symmetric setting akin to Case-2, varying the T e /T i ratio, changing the safety factor q and magnetic shear s and altering the experimental density gradient. The latter two are illustrated in Fig. F.1. Of course, details of edge transport vary with the various parameters scanned. Precise discussion of sensitivity is not within the scope of the current manuscript; the bottomline conclusion is that the reported dynamics of Case-1 (modified outer edge linear stability, nonlinear destabilisation and ensuing transport barrier onset) requires the combined possibility for plasma-boundary interplay and turbulence spreading. Absence of one or both prevents as in Cases-2 and 3 transport barrier onset and leads to transport shortfalls in NMsL of varying severity.

F.0.2 Linear stability of the plasma edge

Linear stability analysis of the poloidally-symmetric and limited GYSELA profiles is performed using the initial value framework of the Gyrokinetic Workshop (Gkw) code (11), Figure F.1: Series of flux-driven computations at ρ ⋆ = 1/316, akin to Case-2 but for the outer boundary conditions. In lieu of the penalised poloidally symmetric SOL of Case-2, an outer diffusive buffer region between 1 ≤ r/a ≤ 1.3 surrounds the confined plasmas through application of the diffusion operator D in Eq.(A.2), detailed as Eq.( 50) in ( 24). The goal is twofold: assess sensitivity (beyond experimental uncertainties) of edge transport dynamics to variations or uncertainties in input profiles and evaluate incidence of surrounding artificial diffusion for edge transport. Density gradients have a stabilising effect on the dominant Ion temperature Gradient (ITG) instability. Subplots (a) and (b) display incremental relaxation of density gradients past r/a ≥ 0.6 and their incidence on fluctuation levels. Similarly, lower magnetic shear and safety factor values past r/a ≥ 0.6 modestly increase fluctuations there. In all cases, edge fluctuation levels are lower than in comparable Case-2. In the absence of outer edge turbulent activity, core to edge spreading appears subdominant. Technically, the choice of an outer diffusive boundary region, though efficient numerically is thus not innocuous (and probably quite poor) on physical grounds. based on the gradient-driven and local (flux-tube) approximations. Unless stated otherwise (see Strategy III below) Boltzmann electrons have been assumed, as in Gysela. Growth rates for the most unstable poloidal wavenumbers k θ ρ i = 0.6 in the poloidally-symmetric (Fig. (i) at a given radial location r 0 , compute the set

values of the Gysela plasma parameters, q being the safety factor, s = (r/q)dq/dr the magnetic shear, ν ⋆ is the collisionality, U ′ the parallel flow shear, γ E denotes the E × B shear, R the tokamak major radius and L -1 X = -d(LogX)/dr with X = {T, n} the local logarithmic gradients of respectively ion temperature and density;

(ii) within Gkw, the local approximation requires mean gradients to be constant over the computational domain, the numerical representation of a flux tube. Coarse-graining of the flux-driven Gysela values is thus performed over a typical radial turbulence length scale ∆r = 10ρ i and in time over the observed linear growth ∆t Ω ci = 25 10 3 of turbulent fluctuations. Furthermore, as the computational domain of GKW winds around the torus parsing both poloidal and toroidal angles, the background state is assumed to be uniform. This implies poloidal homogeneity along the flux tube and requires further averaging the Gysela values on a flux-surface (θ, φ): For all Cases, the dominant instabilities inside r/a ≤ 0.84 are found to be of interchange character. With a Boltzmann electron response, the ion temperature gradient (ITG) is dominant. With a kinetic electron response the instability inside r/a ≤ 0.84 is a combination of ITG and Trapped electron modes (TEM).

• In poloidally-symmetric Cases-2 and 3, Gkw finds the edge to be marginally stable at vanishing E × B shear:

When including E × B shear, strategy I predicts the edge to be nonlinearly unconditionally stable past r/a ≥ 0.9, with γ eff < 0 for all radial-poloidal combinations considered. Strategies II and III confirm that location (r/a = 0.96, θ = 9 • ) is indeed stable. Contributions of

n | ≈ 0, when averaged on a flux surface. Here M ∥ is the parallel Mach number. Local values however of parallel velocity shear about the limiter may locally reach up 3 to 5 times the mean with γ PVG ≳ 0, which does not rule out the possibility for localised yet significant free energy sources from parallel shear flow instability.

It is important to note that even though this procedure provides information on the nature of instabilities at play, two major and oft-made approximations in local approaches should be relaxed to accurately interpret flux-driven dynamics near the separatrix. Indeed, the local E × B shearing rates obtained from limited Gysela with respect to poloidallysymmetric cases and (h) to ] are large enough to significantly decrease linear growth, as currently estimated through Gkw. Strategy I balancing maximum growth with E × B shear only provides partial insight into the nature of active instabilities. Secondly, postulating poloidal (parallel) homogeneity as flux tubes wind around the torus is reasonably accurate in the core but clearly less justifiable to describe the near-separatrix in flux-driven limited configurations. Much of the nonlinear destabilisation of linearly stable edge is consequence of the onset of such poloidal inhomogeneities. Furthermore, additional stabilisation mechanisms are possible, such as profile coupling or poloidal shift of the envelope mode. This linear analysis thus likely provides an upper-bound estimate of the actual flux-driven instability growth. It has the merits however to confirm that the edge under poloidally-symmetric boundary conditions is linearly stable and that it gets destabilised locally in the vicinity of the cold sink, emphasising the central role of turbulence spreading in understanding global equilibration of the turbulence dynamics.

Appendix G

The Transfer Entropy method 4.Causal inference. Causality detection in information theory is actively debated (38). The "Transfer Entropy" (T E) method is a simple nonlinear extension of the Granger causality (39), introduced by Schreiber (40) and investigated in magnetised plasmas by Van Milligen et al. (41). The idea behind T E is simple: let's consider a time series (x i ) of realisations of observable X, with 0 ≤ i ≤ n. If one can better predict its next realisation x n+1 using additional data from another time series (y j ) of observable Y with 0 ≤ j ≤ n, then "Y transfers information (i.e. causes) X", or more accurately as "Y forecasts X", which constitutes the definition of causality here. This idea is quantified measuring deviation of transition probabilities from independence, i.e. from a stationary Markov process. In its simplest expression, if processes X and Y are independent, then the following generalised Markov property holds for all 0 ≤ k ≤ n: p(x n+1 |x n-k , y n-k ) = p(x n+1 |x n-k ). The standard notation for conditional probabilities is used here: p(a|b) is the probability of a knowing b. If now processes X and Y are not independent, the ratio of these two transition probabilities provides a measure of how much information they may share. In other words, how much knowing values within Y in addition to past values in X may help to better evaluate next-step x n+1 . This idea leads to the following definition of the Transfer Entropy (T E) from process Y to process X: where k is thus a time lag and represents the k-past of times series X and Y . The summation process is detailed below, in Eq.(G.3). T E can equivalently be recast as a conditional mutual information and represents the additional amount of information that must be added to adequately represent the studied process p(x n+1 |x n-k , y n-k ) with respect to its reference Markov process p(x n+1 |x n-k ). In the absence of information flow from Y to X, the logarithm vanishes as state Y has no influence on the transition probabilities of X. It also follows that T E is directional, i.e. T E Y →X ̸ = T E X→Y , effectively allowing to infer causality between processes X and Y . T E displays interesting properties: it is independent of the relative magnitudes of signals X and Y ; it may apply to either linear and nonlinear regimes; it is easy to evaluate directly in real space rather than in Fourier space and is typically less demanding in terms of statistics than bispectral techniques. Practically, T E is evaluated expressing the conditional probabilities as joint probabilities p(x n+1 |x n-k , y n-k ) = p(x n+1 , x n-k , y n-k )/p(x n-k , y n-k ) and computing the 4 multidimensional probability density functions (pdfs):