Formation, propagation and conversion of transport barriers triggered by dynamical critical gradient in tokamak plasmas

In this work, we propose a reduced model with a dynamical critical gradient to study the formation, propagation, and conversion of the transport barriers. In contrast to the commonly adopted static critical gradient, an evolving critical gradient self-consistently softens the profile stiffness, so as to facilitate the generation of transport barriers. This is especially crucial to the internal transport barrier (ITB) formation. Numerically, we show that the inhomogeneity of turbulent and neoclassical diffusivities can induce the global wave front propagation of the transport barrier. When the heating power ramps quickly, the ITB propagates unidirectionally to the edge region and converts into an edge transport barrier. For slow power ramping, the propagating ITB will bifurcate into bidirectional wavefronts and finally convert into a steady double transport barrier state. Our model uncovers the vital role of a dynamical ‘profile-stiffness’ in depicting the global dynamics of the transport barrier.

Multiple-equilibrium-state is a ubiquitous feature of complex systems.Classical examples include spin glass state in magnetic materials [1] and predator-prey phenomena in ecological * Author to whom any correspondence should be addressed.
Original Content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence.Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.systems [2].Correspondingly, the selection rule of which equilibrium state the system will be in is an essential topic in the field of complex system theory.The magnetically confined plasma system is well-known for its multi-confinement states (i.e.equilibrium states) [3].Typical scenarios include the L-mode (L: low confinement), the H-mode (H: high confinement) and the I-mode (I: improved), etc et al [4,5].An 'order parameter' depicting the confinement state is the transport barrier, measured by the shearing rate of the mean E × B flow.Transport barrier is a key structure for maintaining the high confinement state of magnetically confined plasma systems.When the transport barrier lies in the edge region, it is called edge transport barrier (ETB) state or H-mode.When in the core region, it is called internal transport barrier (ITB) state.Transitions between different types of transport barriers (typically, the ITB and the ETB) cause the conversions of different confinement states.The coexistence of the ITB and the ETB forms a confinement state, dubbed as double transport barrier (DTB) state (D: doubled).In the Deuterium-Tritium plasma experiment, it is observed that the DTB state has benign features, such as enhanced core confinement, lower power threshold and greater fusion power compared to the standard H-mode [6,7].The synergism between the global dynamics (i.e.propagation) and the local dynamics (i.e.transition) is critical to understanding the transition mechanisms of different confinement states.While, existing studies usually only focus on one side of the process, either on the transition [8][9][10][11], or on the propagation [12,13].
It is generally known that the plasma profiles tend to stay in a near marginal state with a critical gradient.This selforganized criticality (SOC) property is the origin of the socalled profile stiffness [14,15].The profile stiffness exerts a strong constraint on the transition of confinement states.Particularly, the plasma core region is free of boundary conditions and hence is more attracted to the marginal state.The formation of transport barrier requires weakening the profile stiffness, so as to break the old SOC state and form a new one.Since the dynamics associated with the SOC involve strong multi-scale couplings between the turbulence fluctuations and the plasma profiles, the use of reduced models is indispensable to identifying the key physics that govern the nonlinear dynamics of the transport barrier.In this respect, the most comprehensively employed reduced model is the predator-prey model [8,10,11].For L-H transition, it includes two feedback loops: the turbulence-zonal flow (fast) and the turbulence-mean flow (slow).The fast loop is the trigger of the L-H transition [11], while the slow loop facilitates the final establishment of the ETB.To study the global dynamics of the transport barrier, we need continue the 1D predator-prey model in the edge region [11] to the entire system, i.e. from the core to the edge.An important, yet not sufficiently explored element in the existing reduced models is how to address the profile stiffness, properly.It is conventionally modeled as an adjustable parameter [16], i.e. a critical gradient, exceeding which the turbulence bursts abruptly.However, there has been experimental evidence suggesting that the profile stiffness is correlated with the rotational flow, i.e. by increasing the flow shear, the profile stiffness could be softened [17,18].Thus, the critical gradient should be consistent with the flow evolution and it is a dynamical quantity, other than an adjustable parameter.
In this work, we present a 1D reduced model to study the global dynamics of the transport barriers.By incorporating a dynamical critical gradient (DCG), we demonstrate the spatial dependence of the profile stiffness induced by a mean flow shear, and observe the ITB formation, where the profile stiffness is substantially reduced.A key finding is the spatiotemporal bifurcation of the wavefront propagation of the ITB, with the front being driven by the gradient of the turbulent and the neoclassical diffusivities.Depending on the ramping rate of the heating power, the ITB can either convert into an ETB (fast ramping) or a DTB (slowly ramping).
Dynamical critical gradient.The evolution equation of the turbulence intensity I (e.g. the density fluctuation [8]) has a generic reaction-diffusion structure: ∂ t I = S I + r −1 ∂ r (rχ I ∂ r I), with χ I the turbulent diffusivity.The reaction term S I includes various local interactions, which are in turn classified into a driving effect ∝ I and three nonlinear dissipation effects: selfinteraction ∝ −I 2 , micro-mesoscale interaction ∝ −E ZF I and micro-macroscale interaction ∝ −E V I.For modeling, E ZF and E V are usually set as ⟨V ′ ZF ⟩ the zonal flow shear rate and the mean E × B flow shear rate.The turbulence intensity evolution then follows as [11]: where r and t are normalized by the small radius a and a/c s , respectively, with c s = √ T e0 /m i the ion acoustic velocity.The coupling coefficients (α I , ∆ω, α ZF and α V ) depend on the specific instabilities, and their typical ranges have been comprehensively explored [8,10,11].The turbulent diffusivity has a form of χ I (I) = χ N I β , where β = 1 in weak turbulence approximation [19] and χ N /χ GB = O(1) with χ GB the gyro-Bohm diffusivity.
A standard way to model the SOC property of the system is utilizing the critical gradient hypothesis: where ϵ T = R/L T is the normalized temperature gradient, R is the major radius, H(x) is the Heaviside function and γ I is renormalized by c s /R.So that, as ϵ T approaches the upper boundary of the critical gradient ϵ T,c .The divergence of ∂ ϵT γ I reflects the strong stiffness of the plasma profile [20].The L-H transition corresponds to a Dimits shift, such as ϵ T,c → ϵ T,c + ∆ϵ T,c with ∆ϵ T,c induced by the flow shear.Via ∆ϵ T,c , the system is locked to a new SOC state with a steeper critical gradient.However, an important question (yet not addressed in the critical gradient model) is how does the profile stiffness get weakened?To address this question, a logical step is taking ∆ϵ T,c as a dynamical shift, other than a parameter.∆ϵ T,c is caused by the flow shear (i.e.E ZF and E V ) and the flow shear is coupled to the turbulence intensity through the above mentioned feedback loops, so that ∆ϵ T,c must be a dynamical quantity.To see this explicitly, we consider the feedback of γ I on ϵ T,c , and an infinitely-nested functional relation is then obtained: The stiffness near ϵ + T,c becomes: The rapidly growing of It is rightaway to see: With ϵ ′ ′ T,c > 0, the stiffness is reduced by a factor . The rapidly growing of γ I with respect to ϵ T tends to invalidate the perturbation expansion.As a workable model, according to the feedback loop relation: , which can be understood as a dynamical Dimits shift.The DCG model then follows as: where ϵ T,c ≡ ϵ T,c (0).A stimulating way to look at the DCG effect is reexpressing equation ( 4) as Therefore, the stiffness becomes , where ∂ ϵ T, eff γ I is the stiffness of equation ( 2).Generally ∂ ϵT E V > 0, so the original stiffness is reduced by a factor of (1 − α 0 ϵ T,c ∂ ϵT E V ).The coefficient α 0 in equation ( 4) is not a free parameter, but can be calculated when the feedback loops reach a fixed point and it is (see the supplement for details): where E V,0 is mean E × B shear at the equilibrium state, and . The evolution of E ZF has a general form [10,11]: in which the turbulence decoupling effect by E V is modeled as a factor (1 + ζ 0 E V ) −1 , with ζ 0 represents the inhibition of zonal flow growth by mean flow shear [10,21].γ d is the collisional drag, which is a linear saturation mechanism of the zonal flow.Our numerical study shows that adding a nonlinear saturation term in equation ( 6) (e.g.∝ E 2 ZF ) does not significantly impact the conclusion, so we adopt a simple description.By the radial force relation: (the toroidal ion flow is not included here), E V is coupled to the particle, thermal and poloidal momentum transports through and where ρ s is the Larmor radius and the coefficient α θ can be found in [22,23].The neoclassical viscosity is µ neo = ν ii, * q(r) 2 R 2 , with ν ii, * the ion-ion collision frequency.The inward particle pinch is ) [24,25].The heat source is deposited in the core region ), and the density source in the edge region S n = S n, max exp(−(1 − r) 2 /2L 2 n,dep ), with L p,dep and L n,dep the deposition width.Equations ( 1) and ( 6)-( 10) constitute the reduced model with a DCG.More details on the coefficients can be found in [11].It will be shown immediately that the DCG and the inhomogeneity of the transport coefficients play an essential role in the formation, propagation and conversion of the transport barriers.
Figure 1 shows how the DCG impacts the local stiffness of the plasma profile.For E V ≈ 1.5 × 10 −4 , the stiffness parameter ∂ ϵT γ I tends to be divergent near ϵ T,c .At E V = 0.2, the profile stiffness is substantially reduced, so that it facilitates the formation of the ITB.Figures 2(a)-(c) are the steady profiles for different α 0 .The heating power is below that of the L-H transition (without DCG effect).For α 0 = 0, the system is in L-mode state (black lines).For a finite α 0 , ITBs form both in the density-and the temperature profiles and its spatial extension will expand as α 0 increases.The nonlinear stiffness mitigation provides a noticeable improvement of the temperature in the core region, compared with the standard critical gradient model.The core density only increases mildly because the reduction of the inward particle pinch velocity and the lack of density source in the core region.Figures 2(d)-( f ) show the mean flow profiles and their components across three cases, indicating that ⟨V E×B ⟩ peaks in ITB and ETB regions and its dominant contributor is diamagnetic flow (−ρ s L −1 p ).The nonlocal evolution of the ITB is driven by two processes: (1) a pure diffusion that tends to flatten the transport barrier; (2) a wavefront propagation of the plasma profile induced by the spatial inhomogeneity of the diffusivities [12].First, we study the neo-classical diffusivity induced wavefront effect.Since −∂ r χ neo = −2ν ii R 3/2 ρ 2 s qq ′ , a positive magnetic shear (q ′ > 0) will induce an inward front, while a negative magnetic shear will induce an outward one.As shown in figure 3, if q is uniformly upshifted by ∆q = 0.1 (black solid line), the ITB is weakened due to the increase of χ neo .When increasing q ′ on the right of the inflection point of the pressure profile, an accelerated left-propagation wavefront would steepen the ITB.When adding a reversed magnetic shear, head-on wavefronts on the two sides of the inflection point form (subplot in figure 3(b)) and ITB is further enhanced.The wavefront propagation due to ∂ r χ neo also provides an interpretation of how the reversed magnetic shear configuration facilitates the formation of ITB [26,27].
In light of understanding the global, spreading and conversion of the ITB, we linearly elevate the input power from Q p, max = 5 × 10 −4 to Q p, max = 1 × 10 −2 in the duration (t = 1 × 10 4 , t = 6 × 10 4 ).The spatiotemporal structure of the transport barrier can be identified through the 'valleys' of I or the 'peaks' of E V .Figures 4(a)-(c) shows that for α 0 = 0, a standard L-H transition occurs through limitcycle-oscillations (LCOs) and it belongs to a 'soft' L-H transition.By increasing the ramping rate of the heating power, a stimulated L-H transition can occur, i.e. a 'hard' transition, as has been reported in [11,28].In the DCG model, with an initial steady ITB state, for the case of slow power ramp (figures 4(d)-( f )), LCOs in the I − E ZF feedback loop appear first, but having a short duration.Then an avalanchelike outward propagation of the ITB is triggered and causes a global quench of the turbulence (figure 4( f )).However, the fully quenched state is only transient and the front of the transport barrier immediately bifurcates into an inward-and an outward wavefronts.The outward branch propagates into the edge region and forms an ETB.The inward one propagates back to the core region and eventually restores the ITB, accompanied by continuously splitting off subfronts that propagate into the edge region and strengthen the ETB.Finally, a stable DTB state builds up (figure 4( f )).For quick power ramp, the DTB disappears.Instead the initial ITB propagates outward and finally converted into an ETB, i.e. the H-mode.The corresponding initial and final steady-state profiles are summarized in figure 5, with figures 5(d)-( f ) revealing a stronger mean flow in the ETB region than in the ITB.Additionally, figure 5( f ) exhibits a reduced core mean flow in the steady DTB state relative to the low-power ITB state.
The detailed mechanism of the front propagation of the transport barrier can be glimpsed by inspecting the field configuration of ∂ r χ turb .With equation ( 7), the mean E×B shear is decomposed into , where R p = d(ln p ′ )/dr the curvature of the pressure profile and L f = d(ln f)/dr with f = p, n and ⟨V θ ⟩.Taking the ITB to DTB transition for illustration, figure 6(a) shows that the curvature of the pressure profile dominates the , except for a thin layer in the edge region.Since |∂ r χ neo | ≪ |∂ r χ turb | when the turbulence is excited, we need only consider ∂ r χ turb -driven propagation.The bidirectional propagations of the transport barrier are plotted in figures 6(b)-(g), respectively, which corresponds to t = 24 000 − 26 400 and t = 31 000 − 31 400 in figure 4( f ).Alongside the propagation, it clearly shows that the head-on wavefronts tend to steepen the profile and the back-propagating wavefronts tend to flatten the profile.
Figure 7 illustrates the development from the standard L-H transition to the 'avalanche' scenario as the nonlinear stiffness coefficient α 0 increases.For α 0 ⩽ 5, the DCG effect is weak and ETB is formed through LCO, i.e. the standard L-H transition, and an example is given for α 0 = 5 (figure 7(a)).Figure 7(b) is an intermediate situation where both ITB and ETB are formed via local transport processes.As α 0 increases and exceeds a certain threshold, the ITB 'avalanche' events occur and the DTB states are maintained for both power ramping rates.Until α 0 reaches a relatively large value (e.g.For α 0 = 8.5), the DTB becomes unstable for the quick power ramping and it induces a DTB→ETB transition.
Finally, we explored how collisions affect the transition process, via the zonal flow damping rate γ d , the poloidal flow's neoclassical viscosity µ neo , and neoclassical diffusivities χ neo /D neo .These quantities are proportional to the collisionality ν ii, * .Figure 8 revisits the cases in figures 4(d)-(i) under varying collisionalities.With weaker collisions (figure 8(a)), the system struggles to achieve a steady state, experiencing intermittent ITB 'avalanche' events.This occurs because the DCG effect is more prominent than the collisional effect, akin to how turbulence in fluids becomes more prevalent with increased Reynolds numbers.Figures 8(a   increasing collisionality.Specifically, at ν ii, * = 3.6 × 10 −3 , ITB is difficult to stabilize at elevated input power even with a slow ramping rate (figure 8(d)).Further increases in collisionality ultimately revert the system to the standard local L-H transition (figure 8(e)).
In this work, we mainly focus on the physics of the DCG effect.The toroidal rotation is an important factor in producing the radial electric field through radial force balance relation.To include it self-consistently, one needs to incorporate the toroidal momentum evolution equation, which is beyond the current work, but is definitely doable.This is also a subject we will pursue in future work.It's worth noting that while not accounted for in the current model, MHD stability, bootstrap current evolution (affecting the q-profile and its shear), Shafranov shift, and α-stabilization, along with related parameter β p , are also expected to be crucial for sustaining high plasma performance in high-pressure (gradient) regimes.These elements are essential for a comprehensive evaluation of a feasible fusion scenario and should be incorporated into the future expanded version of the current DCG model.
This project was supported by the National MCF Energy R&D Program (2018YFE0311400) and the Natural Science Foundation of China (NSFC) No.12 075 013.     ) where t in = 10 4 is the particle injection time, τ SBMI = 250 is the duration of fuelling, r dep = 0.97 is the deposition location, ∆r = 0.05 is the deposition width and the strength of particle injection is ∆S SMBI = 3.2 × 10 −3 .The initial state is subcritical and is in the L-mode with input power of Q p,max = 5 × 10 −4 ≈ 0.5Q th .Here, Q th represents the threshold power for L-H transition (α 0 = 0) under the parameters we adopted.The particle injection leads to an abrupt increase in local density, and the mean flow shear is immediately excited.This results in the suppression of edge turbulence and the subsequent continuous attenuation of zonal flow.LCOs are absent in such stimulated transition.Compared to the initial state, the density and temperature remain higher even after the injection.These results are consistent with the published work [28].

Figure 2 .
Figure 2. Steady-state profiles of (a) density n, (b) temperature T and (c) pressure p for different α 0 .The dashed lines are their gradients.The colored shadow regions label the sampling regions in the figure 1. Subplots (d)-( f ) show the corresponding mean flow profiles and their components.

Figure 3 .
Figure 3. (a) Temperature gradient variation ∆ϵ T = ϵ T − ϵ T,q ref for different q-profiles with q ref = 1 + 2r 2 (dashed line).(b) The corresponding q-profiles.The solid black line is an up-shift of q ref .
figure 4( f ).Alongside the propagation, it clearly shows that the head-on wavefronts tend to steepen the profile and the back-propagating wavefronts tend to flatten the profile.Figure7illustrates the development from the standard L-H transition to the 'avalanche' scenario as the nonlinear stiffness coefficient α 0 increases.For α 0 ⩽ 5, the DCG effect is weak and ETB is formed through LCO, i.e. the standard L-H transition, and an example is given for α 0 = 5 (figure7(a)).Figure7(b) is an intermediate situation where both ITB and ETB are formed via local transport processes.As α 0 increases and exceeds a certain threshold, the ITB 'avalanche' events occur and the DTB states are maintained for both power ramping rates.Until α 0 reaches a relatively large value (e.g.For α 0 = 8.5), the DTB becomes unstable for the quick power ramping and it induces a DTB→ETB transition.Finally, we explored how collisions affect the transition process, via the zonal flow damping rate γ d , the poloidal flow's neoclassical viscosity µ neo , and neoclassical diffusivities χ neo /D neo .These quantities are proportional to the collisionality ν ii, * .Figure8revisits the cases in figures 4(d)-(i) under varying collisionalities.With weaker collisions (figure8(a)), the system struggles to achieve a steady state, experiencing intermittent ITB 'avalanche' events.This occurs because the DCG effect is more prominent than the collisional effect, akin to how turbulence in fluids becomes more prevalent with increased Reynolds numbers.Figures8(a)-(e) clearly demonstrate a gradual weakening of the DCG effect with

Figure 5 .
Figure 5. (a) Density, (b) temperature, and (c) pressure initial/final steady-state profiles corresponding to the cases shown in figure 4. Subplots (d)-( f ) compare the profiles of the components of the mean flow at steady state under the high input power (Qp,max = 10 −2 ).

Figure 6 .
Figure 6.(a) The three components (ρsL −1p R −1 p , −ρsL −1 p L −1 n , and − ⟨V θ ⟩ c −1 s L −1 V θ ) of the mean E × B shear.Subplots (b)-(d)and (e)-(g) display the inward and the outward fronts of transport barrier, respectively.The black solid line is the pressure profile, and the red dashed line is its curvature.Blue arrows indicate the amplitude and the direction of wavefront speed.

Figure 7 .
Figure 7.The spatiotemporal evolution of E V with different α 0 .Other parameters are same as the cases in figures 4(d)-(i).Slow ramp of input power is applied in (A1)-(E1), and quick one is applied in (A2)-(E2).The white lines label the duration of power ramping.

Figure 8 .
Figure 8.The spatiotemporal evolution of E V with different collisionalities.Other parameters are same as the cases in figures 4(d)-(i).Slow ramp of input power is applied in (A1)-(E1), and quick one is applied in (A2)-(E2).The white lines label the duration of power ramping.

Figure 10 .
Figure 10.Lissajou curve between zonal flow and mean flow from t = 45 250 to 49 000 at the radial position r = 0.973.The arrows label the direction of temporal evolution.

Figure 11 .
Figure 11.The sptiotemporal evolution of (a) I, (b) E ZF and (c) E V for L-H transition induced by SMBI.The subplots (d)-( f ) exhibit the evolution of n, p, I, E ZF and E V at r = 0.962.