Transport in fusion plasmas: is the tail wagging the dog?

Turbulent plasmas notably self-organize to higher energy states upon application of additional free energy sources or modiﬁcation of edge operating conditions. Mechanisms whereby such bifurcations occur have been actively debated for decades. Enhanced conﬁnement occurs at the plasma edge, where a shortfall of predicted turbulence intensity has been puzzling scientists for decades. We show, from the primitive kinetic equations that both problems are connected and that interplay of conﬁned plasma turbulence with its material boundaries is essential to curing the shortfall of predicted turbulence and to triggering spontaneous transport barrier onset at the plasma edge. Both problems determine access to improved conﬁnement and are central to fusion research. A comprehensive discussion of the underlying mechanisms is proposed. These results, highly relevant to the quest for magnetic fusion may also be generic to many problems in ﬂuids and plasmas where turbulence self-advection is active.

(Dated: June 10, 2022) Turbulent plasmas notably self-organize 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 at the plasma edge, where a shortfall of predicted turbulence intensity has been puzzling scientists for decades. We show, from the primitive kinetic equations that both problems are connected and that interplay of confined plasma turbulence with its material boundaries is essential to curing the shortfall of predicted turbulence and to triggering spontaneous transport barrier onset at the plasma edge. Both problems determine access to improved confinement and are central to fusion research. A comprehensive discussion of the underlying mechanisms is proposed. These results, highly relevant to the quest for magnetic fusion may also be generic to many problems in fluids and plasmas where turbulence self-advection is active.
In the quest for magnetic confinement fusion, field geometry plays an important role. Magnetic field lines in tokamaks or stellarators are built such as to trace out closed toroidal flux surfaces with a high degree of symmetry. Field symmetry is known to bolster confinement, enabling entrapment of the ionised plasma. Symmetry breaking however is common and usually results in net particle, energy or momentum sinks and ultimately in degradation of confinement. In particular, toroidal Established practice oft distinguishes between a hot confined 'core' region, dense and hot, an unconfined peripheral boundary layer (the 'Scrape-Off Layer' or 'SOL') and an in-between 'edge' region, loosely defined, set between core and separatrix. The SOL is cold and rarefied; it starts at the magnetic separatrix and is mapped out by the open field lines which connect magnetically to the material boundaries. Core and SOL have been extensively studied, mostly independently, the edge usually serving in modelling as fixed boundary condition for both, its dynamics difficult to apprehend.
Strict decoupling however between all three regions is increasingly being questioned. Tokamak plasmas are indeed prone to self-organisation and mounting evidence suggests interplay between core, edge and SOL. 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.
The near-separatrix edge region is a cornerstone of fusion research, where spontaneous transitions from low-confinement "L-mode" to high-confinement "H-mode" occur [1]. The H-mode branch of operation is one of several improved confinement states that have been experimentally discovered, revitalizing the fusion program towards ITER. Microturbulence dominates the transport processes in the L-mode edge. Description from first principles of the transport processes there, prior to bifurcations, is largely missing yet key to establishing whereby transport bifurcations occur.
We discuss a generic situation, based on experimental parameters where the following conundrum is found: experimentally, the edge is measured to be turbulent, with fluctuations increasing with proximity to the separatrix [2]. In contrast, local analysis of the profiles predicts convective stability in the edge and increased stability with proximity to the separatrix. The edge region would thus appear as unfit to produce or sustain a turbulent state. This opposes the experimental trend and precludes the possibility for turbulence-induced bifurcations to improved confinement. This problem was first fomulated decades ago and is still often referred to as the "transport shortfall" conundrum.
We show, from the primitive kinetic equations, a possible resolution for this problem. Understanding the origin of turbulence activity in the edge requires considering the interplay between closed and open field lines. Magnetic connection to the material boundaries deeply modifies global convective stability at the separatrix. An additional source of free energy arises there, resulting in the confined plasma being driven unstable. Fluctuations, initially localised in a narrow peripheral area of the plasma at the plasma-material boundary interface expand beyond their region of convective drive and "spread" throughout the stable edge. Surprisingly, this narrow peripheral region controls the dynamics of much larger plasma core and edge through fluxes of turbulence activity, a case of the tail wagging the dog.
Model equations. -Low-frequency microturbulence in weakly collisional magnetised plasmas is appropriately described within the gyrokinetic framework [3]. The GYSELA code [4] solves the governing coupled gyrokinetic: and quasi-neutrality: equations for the guiding-center distribution functionF s 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: with J µ . the gyro-average operator. Notations are those of Ref. [4]. 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 [5] and penalised [6,7] 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 Ref. [8] and built such as to prevent overdamping zonal modes [9].
Penalisation [6] 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 , 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. This setting is easily extrapolable to a kinetic response for electrons inside the separatrix (for r/a ≤ 1), retaining an adiabatic electron response in the SOL. Relaxation of this assumption in the SOL is currently under investigation [10].
Modelling strategies: flux-and gradient-driven approacheswhich has three major implications: (i) as touched upon above, a flux-driven framework is required implying that profiles, flows, stresses and flux patterns are unknowns of the dynamical problem and the main outcomes of the numerical experiment. The plasma self-organises on a fraction of a confinement time and alters its internal profiles in a self-consistent fashion whilst retaining memory of this reorganisation. It opens the possibility to dynamically describe a bifurcation. In addition, (ii) the full volume of the plasma from dense core to rarefied SOL next to the material boundaries is resolved. Each region of the plasma volume is a priori dynamically connected to the rest so that all parts may interplay. At last, (iii) no temporal and spatial scale separation is assumed as many relevant scales commonly associated to 'mean' quantities or to 'fluctuations' do interplay, from the short ion-scale Larmor radius (ρ i ∼ 10 −3 m) of turbulent eddies driven by Ion Temperature Gradient Such properties make the problem computationally intensive yet are instrumental to the results discussed here. The multiple length and time scales are indeed known to interplay in the core, e.g.
through predator-prey-like dynamics between zonal flows and turbulence [11] or through the onset of staircases [8]. As characteristic scales in the edge all become comparable due to the large gradients of mean plasma quantities when approaching the separatrix, separations either in time or in space, though computationally favourable in gradient-driven models become questionable. Furthermore, the modelling perspective on turbulence changes between the plasma core and its periphery. In the confined core, provided known constraints such as the minor radius of the torus, its aspect ratio or its magnetic configuration, the central focus is on predicting the energy, momentum or particle confinement time. Diffusion in that matter is the oft-invoked paradigm for the analysis of transport.
The question of estimating a confinement time is thus often recast as a problem of estimating heat, momentum or particle diffusivities. Through the argument that transport is stiff, mean gradients are thus often taken out of the dynamic feedback loops that determine transport and are assumed to be known and to act as reservoirs of free energy. Gradient-driven approaches thus allow to compute a transport matrix at a reduced numerical cost. They are the current workhorse for estimating transport in the confined core of turbulent fusion plasmas, documented objections to their implicit assumptions of locality and diffusive transport notwithstanding [12,13]. However, when crossing the separatrix the situation is drastically modified. The confinement time is no longer the primary issue in the vicinity of the separatrix and in the scrape-off layer for it is tightly constrained by the magnetic connections to the boundaries. Cross-field transport impacts flux deposition patterns, which is of primary concern for fusion. As sources and sinks are intrinsically not uniformly distributed along the field lines, the determination of the poloidal distribution of gradients is a central question, as well as the radial propagation of turbulence on both sides of the separatrix. The edge-SOL interface thus needs addressing in a flux-driven framework.
Primitive modelling is especially useful if it strives to remain predictive, not relying upon stiffness arguments, adjustable gradients to match experimental observations. In that respect, flux-driven strives to do just this: to remain predictif. Linear stability analysisof the poloidally-symmetric and limited GYSELA profiles is performed using the initial value framework of the Gyrokinetic Workshop (Gkw) code [14], 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.A-3) and limited (Fig.A-4) configurations are estimated by the following procedure.
(i) at a given radial location r 0 , compute the set Σ r 0 = {q, s, ν , T i /T e , U , γ E , R/L T , R/L n } of local 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  the instability drive as if located at its ballooning angle, effectively maximising it; (iii) for a chosen radial location r j , knowing Σ r j ∆t ∆r θ ϕ allows to compute with Gkw 2dimensional maps of instability growth rates γ lin as a function of the logarithmic gradients (iv) we now estimate local values of the Gysela local growth rates γ lin at 13 different radial-poloidal (r j , θ k ) locations (shown in Fig.1 as well), Causal inference. -Causality detection in information theory is actively debated [16]. The "Transfer Entropy" (T E) method is a simple nonlinear extension of the Granger causality [17], introduced by Schreiber [18] and investigated in magnetised plasmas by Van Milligen et al. [19]. 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. (A-8). 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): as a function of time delay k and normalised such that 0 ≤ T E ≤ 1. The 4 pdfs in Eq.(A-7) result from a binning process of times series X and Y , such that Eq.(A-7) is estimated in practice as: 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, , which provides the net flow of information between processes X and Y , at timelag k. In the manuscript, pdfs in Eq.(A-7) 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. [20].
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 [20][21][22] 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.(A-9) and especially here to the following set: (X, Y ) ∈ Ω r , v Er Ω r , v r Ω r , − v θ 1 r ∂ θ E r (A-14) Conclusions. -Methods described here have allowed to establish in Ref. [23] the following: (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.
The expanding interface organises into a stable peripheral transport barrier, i.e. into improved confinement.
This certainly stresses a key role for finite Larmor radius effects, i.e. for perpendicular pressure fluctuations. Though challenging, experimental measurements of such fluctuations and of their associated radially-inward, poloidally-anisotropic advections (spreading) could prove especially important to diagnose to deeply test our understanding of underlying mechanisms in the plasma edge, especially prior or at transport barrier inception. * guilhem.dif-pradalier@cea.fr