On the possibility of limit-cycle-state of peeling mode near stability boundary in the quiescent H-mode

A model is proposed for the edge harmonic oscillation, in which the stationary coherent mode is sustained in the almost linear phase as has been observed in JT-60U. We study the coupled dynamics of the peeling mode amplitude and edge pressure gradient. The limit cycle oscillation is predicted. The peeling mode (which is almost in the linear phase) is in a dynamical stationary state with amplitude modulation. In this model, the time scales for the change of parameters that specify magnetic structures (such as magnetic shear and edge plasma current) are assumed to be much slower, so that are decoupled from the limit cycle dynamics. The condition that the limit cycle state appears is shown. The oscillation frequency of the modulation is given by the hybrid mean of the typical growth rate of the peeling mode and the additional loss rate of pressure gradient by the peeling mode.


Introduction
The high confinement mode (H-mode) in magneticallyconfined toroidal plasmas [1] has brought about a large impact on the controlled thermonuclear fusion research, by demonstrating that the energy confinement time can be 7 Deceased on 18 July 2019. * Dedication Prof. Sanae-I. Itoh has passed away before completing this article. The rest of authors dedicate this article to the memory of her.
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. enhanced substantially. At the same time, the H-mode has shown a problem, i.e. the edge-localized mode (ELM) appears, which is a repetitive pulse of strong energy and particle losses from the H-mode plasma [2,3]. Such intermittent strong loss is a large obstacle in designing the fusion experimental reactor like ITER [4]. Many kinds of ELMs (from phenomenological classification) have been reported (see a review, e.g. [5],). Because of the many variations of ELMs, several kinds of theoretical models have been proposed: some are focusing on magnetohydrodynamic (MHD) stabilities [6][7][8][9], or on limit cycles combined with transport bifurcation [10][11][12][13] (see a review, e.g. [14].). Nevertheless, comprehends understanding of ELMs, in general, is still unsatisfactory. Among various ELMs, the Type-I ELM [2] is considered to be most serious, and the importance of MHD peeling mode has attracted attentions.
As a typical example of the H-mode plasmas free from ELMs, the quiescent H-mode (QH-mode) has been observed in various experimental devices [15][16][17][18]. In the QH-mode, a coherent MHD mode that is localized near the transport barrier appears, and is called edge harmonic oscillation (EHO) [19,20]. It has been conjectured that the additional plasma loss owing to the EHO influences the transport barrier so that the ELMs do not appear in the QH-mode. Plausible candidates of EHO are the peeling mode or peeling-ballooning mode [9,[21][22][23]. From experimental point of view, the parameter region where the QH-mode appears is close to the instability boundaries of peeling mode or peeling-ballooning mode [24,25], so that the EHO is considered to be related with the peeling mode or peeling-ballooning mode. Theories of the peeling mode and peeling-ballooning mode have been elaborated by taking into account of the effects of E × B flow [26][27][28][29][30][31] and of ion diamagnetic drift [32][33][34].
In the experiments on DIII-D, ASDEX-U and JET, the EHO appears as a nonlinear wave, which is accompanied by many higher harmonics. This suggests that the nonlinear saturation mechanisms by higher harmonics keep the EHO at almost constant amplitude. One hypothesis in the nonlinear saturation mechanism of highly nonlinear wave is the phaseslip model, in which possibility of the phase slip between the pressure and velocity perturbations are analyzed [35]. A test of this model in DIII-D experiments has also been reported [36].
Motivated by these important roles of the peeling mode, an identification of the peeling mode has been reported by JT-60U experiments [37]. Details of the data analysis are reported in a separate article [38]. The distinct features of the QH-mode in JT-60U are that (a) the Q-H mode appears in the cases of co-and counter injection of neutral beam injection (implying a relatively weak dependence on plasma rotation or its shear), and that (b) the EHO is not associated with higher harmonics (at most with an weak second harmonics). In other words, the peeling mode keeps stationary amplitude almost in a quasilinear state, without strong nonlinear MHD interactions. It is important to understand the mechanism why the EHO is sustained stationary almost in its linear eigenfunction.
In this article, we study the possibility that EHO is sustained in the almost linear phase of peeling mode, and propose a model to explain the observed phenomena in JT-60U. We study the coupled dynamics of the peeling mode amplitude and edge pressure gradient, by introducing a time-scale separation. That is, the time scales for the change of parameters that specify magnetic structures (such as magnetic shear and edge plasma current driven by the steep pressure gradient in the transport barrier) are assumed to be much slower than the dynamical time of the pressure gradient and peeling mode amplitude. As a result of the assumption of time scale separation, we show that the pair of edge pressure gradient and peeling mode amplitude shows a limit cycle oscillation. The peeling mode (which is almost in the linear phase) is in a dynamical stationary state with amplitude modulation. The condition that the limit cycle state appears is shown, and the period of limit cycle oscillation is predicted. The oscillation frequency of the modulation is given by the hybrid mean of the typical growth rate of the peeling mode and the additional loss rate of pressure gradient by the peeling mode. Issues for the future experimental tests are also identified. At the end of this article, the effects of inhomogeneous E × B flow and of ion diamagnetic drift are discussed.

Time scale separation
In this article, the dynamics of the peeling mode near the stability boundary in the H-mode plasma is investigated. The geometry of the plasma and magnetic structures (minor radius, a, major radius, R 0 , safety factor, q, and other elements, etc) are considered as prescribed. Under such conditions, the growth and decay of the peeling mode are controlled by the parameters near the plasma boundary (s, α, J edge ): They represent the magnetic shear, normalized pressure gradient and the edge current, which is driven by the Bootstrap current associated with the steep pressure gradient in the transport barrier, respectively. The normalized pressure gradient is defined as α = −2q 2 R 0 µ 0 B −2 dp/dr (q: safety factor, R 0 : major radius, B: magnetic field density, p: plasma pressure, r: minor radius). As a fundamental assumption, we introduce the time scale separation. That is, the time scale for the change of equilibrium magnetic field is much longer than those of the amplitude of the peeling mode and the gradient of plasma pressure. Along this line of thought, we here separate the parameter α out among the controlling parameters (s, α, J edge ), and assume that s and J edge are constant in time. The effects of inhomogeneous E × B flow and of ion diamagnetic drift are neglected in this article for the transparency of the arguments. Influences of these additional effects are discussed at the end of this article; it is shown that the conclusion in this article is qualitatively unchanged.
Based on the assumptions of the time scale separation, we study the coupled dynamics of peeling mode amplitude and the pressure gradient near edge, α, and other parameters are treated as given constants.

Model of the growth rate of the peeling mode
Motivated by the observation on the JT-60U, in which the amplitude of the peeling mode is small enough and the parameters are close to the stability boundary, the evolution of the peeling mode is modeled by the linear growth. We choose the situation, which is illustrated in figure 1, i.e. for fixed values of s and J edge , the growth rate of the peeling mode is larger for larger pressure gradient near the stability boundary following [34]. Based on the assumption of time scale separation, the dependence of the growth rate of the peeling mode, γ, on the pressure gradient is kept as where I is the amplitude of the peeling mode. Here we employ a model form, by keeping the lowest order term in Taylor expansion, as where α c is the critical pressure gradient of the linear instability, and γ 0 is a typical value of linear growth rate of unstable peeling mode. (Note that in ideal MHD calculations, growth rate can have singular dependence on the controlling parameter at the stability boundary, e.g. Nevertheless, if one takes into finite dissipations and resonances, etc, the leading term in the Taylor expansion of the growth rate becomes like equation (2) [39][40][41]. We employ equation (2) for the transparency of the arguments.) This simplified model is used to characterize the dependence on the pressure gradient and to introduce the shortest time scale, γ 0 −1 .

Evolution of the pressure gradient
For the evolution equation of the pressure gradient, we must employ the model of the heat flux. We here introduce a very crude model that the heat flux is given by the diffusive model. Noting the limitations of diffusive model of transport in tokamaks, we employ this simple one owing to the following reason. If we see the amplitude of the peeling mode in JT-60U, which is of our interest, the radial excursion length of this weak mode is of the order of few mm [38]. This excursion length is much smaller than the width of the transport barrier. Therefore, we take a model that the contribution by peeling mode on the mean radial heat flux is modeled by the enhancement of diffusion coefficient, as the first step of the analysis. The radial heat flux q r is modeled in this article as where K 0 is the modified thermal conductivity in the absence of the peeling mode and K peeling is the increment owing to the peeling mode, respectively. The relation between the standard thermal conductivity and assumptions to employ this model is explained in the appendix A. The power balance equation is rewritten, by taking the radial derivative and appropriate normalization, as where P norm = − 2q 2 R 0 µ 0 B −2 dP/dr. The approximation to deduce equation (6) from equation (5) is also explained in the appendix A. The radial derivative is approximated by the derivative in the Cartesian coordinate, x, because we are studying the thin region near the plasma surface.

Linearization
The deviation of the pressure gradient from the critical value is denoted by the parameter g, and is taken as a smallness parameter here. The modification of the transport coefficient by the peeling mode, K peeling , in equation (4) is also considered as a smallness parameter, in which the lowest order correction with respect to the peeling mode perturbation is kept. Equation (6) is deformed, by neglecting higher order terms of g and K peeling , as dg/dt = S + ∆κ 0 g + ∆κ peeling α c (8) where the first term in the RHS of equation (8), S, is given as and represents the external operation condition that slowly (i.e. in the transport time scale) drives the plasma into the peeling unstable domain. In studying the coupling dynamics of g and the amplitude of the peeling mode, the term S is treated as the prescribed constant. We are studying the peeling mode in the regime of the quasi-linear response, so that the additional thermal conductivity is proportional to the squared intensity of perturbation A as where K peel,0 is the value of K peeling at the peak value of the peeling mode, I 0 , and the squared amplitude is normalized as A = I 2 /I 0 2 . Thus, one has We next approximate the partial differential equation by the ordinary differential equation. For this purpose, the spatial derivative operator in equations (8) and (11) are evaluated by a scalar coefficient. Following the work [37,38], one sees that the radial eigenfunction of the peeling mode in the condition of our interest is localized near the plasma edge. So are the mean pressure gradient and perturbation amplitude. Therefore, the quantity K peel,0 α c A is considered to have a convex shape. That is, the term ∆ K peel,0 α c A takes a negative value. Using a positive coefficient τ peel , it is evaluated as This indicates that the peeling mode peels the plasmas near edge, so that the plasma pressure gradient becomes smaller by the effects of the peeling mode growth.
At this moment of the analysis, we assume that the second term in the RHS of equation (8) is smaller than the third term. This is based on the ansatz that the unperturbed transport coefficient K 0 is suppressed in the transport barrier. In addition, the amplitude of perturbation g is much smaller than that of A, as is shown later in section 3.1. Thus, in the following analysis, the second term in the RHS is neglected (in comparison with the third term). The effect of the second term in the RHS of equation (8) is discussed in the appendix B.
Summarizing preceding arguments, a set of coupling equations for the peeling mode intensity, A, and pressure gradient perturbation, g, are summarized as: and

Integral of equations and amplitude of oscillation
The set of coupled equations (13) and (14) has an integral of motion. Rewriting equation (14) as and substituting it into the last term of RHS of equation (13), one has Rewriting equation (13) as 2(γ 0 / α c ) g = A −1 dA/dt, and substituting this form into the RHS of equation (16), one has This equation can be integrated as where C is a constant of motion. The integral equation (18) shows the limit cycle dynamics of the peeling mode near the marginal stability condition. We choose the constant of integral as C = 1, because the variable A is the squared amplitude normalized to the peak value, i.e. g = 0 at A = 1. One finally has the relation between the perturbed pressure gradient and peeling mode amplitude as When the condition is satisfied, there exists a lower boundary A low, (which is the solution of the equation 1 − A + τ peel S ln(A) = 0 and satisfies the condition 0 < A low < 1), the RHS of equation (19) becomes positive in the region This shows that the limit cycle solution appears. (The case, where equation (20) is not satisfied, is discussed later in section 3.4). Figure 2 illustrates a typical example of the limit cycle solution. The trajectory follows counter-clockwise direction. For instance, if the initial point is chosen at g = 0 and the minimum of A, the gradient increases and the growth of mode amplitude follows. When the mode amplitude becomes large enough, the gradient starts to decrease, and, as the point g = 0 is crossed, the mode amplitude begins to decrease. This leads to the enhancement of g later, reproducing the initial point of trajectory.
The variation of the pressure gradient associated with this limit cycle oscillation can be estimated by equation (19). Under the condition of equation (20), the RHS of equation (19) is O (1). Therefore, the amplitude of the pressure gradient modulation is estimated as The parameter is a square root of the ratio between shorter time (γ 0 −1 ) and longer time, τ peel , so that the perturbation of the pressure gradient remains small. Therefore, this limit cycle allows the dynamical steady state in which the plasma state stays near the stability boundary, even under the condition that the external operation condition drives continuously the plasma into the peeling unstable domain, S > 0, (in the transport time scale).

Period of limit cycle oscillation
From equations (13) and (19), the time evolution can be discussed. Rewriting equation (19) as and substituting this expression into the RHS of equation (13), one has

Equation (23b) shows that temporal evolution A(t) is given as an implicit function; For instance, in the phase where A decreases, one has
where t 0 is chosen as A(t 0 ) = 1. This shows that the amplitude A is given by a nonlinear oscillating function. Even though this is not solved by use of elementary functions, the time scale of oscillation is given by because the integral of RHS of equation (23b) over the one period of oscillation is O (1). It is noted that this time scale of limit cycle oscillation is given by the hybrid time scale, (τ peel / γ 0 ) 1/2 , that is, the hybrid mean of very short time scale (MHD mode growth time), 1/γ 0 , and the intermediate time scale for the peeling of surface plasma by the enhanced transport by the peeling mode, τ peel .

Case study of example
In order to illustrate the qualitative feature of the abovementioned limit cycle, we choose the example of parameters similar to those in JT-60U [37,38] as follows. The amplitude of the peeling mode, when excited, is such that to induce the radial excursion of the plasma as δ ∼ 2 mm. The correlation time is employed based on the full-width half-maximum as ∆f ∼ 0.5 kHz. Combination of these choices leads to the order of magnitude estimate for the thermal conductivity that is induced additionally by the peeling mode [42] as, Taking an example that the width, where the peeling mode has large amplitude, d, is about 2 cm, the rate of the additional change of gradient by the peeling mode, τ peel , is estimated as A characteristic value for the MHD mode may be evaluated as Substituting equations (26) and (27) into equation (24), the characteristic time for the limit cycle oscillation is given as This characteristic time scale is much shorter than the magnetic diffusion time near plasma edge to change the edge current and magnetic shear. Therefore, it is consistent with the assumption of the time-scale separation, by which parameters s and J edge are kept constant in the dynamical change of (g, A).
The modulation amplitude of the perturbation of the pressure gradient is also evaluated by substituting equations (27) and (28) into equation (22) as As is illustrated by equation (29), the pressure gradient is bound near the stability boundary.
For these parameters, the condition equation (20) is rewritten as That is, the global parameters (associated with the pressure gradient) that drive the plasma state into unstable region, develops slower than 50 ms. (20) Before closing this section, we briefly discuss the case where equation (20) is not satisfied. One can show that limit cycle solution exists, even if equation (20) is not satisfied.

On condition equation
When equation (20) is not satisfied, i.e.
S τ peel > 1 (31) holds, the value g 2 , which is given by the integral equation (19), is positive definite in the region where A high is the solution of the equation 1 − A + τ peel S ln(A) = 0 and satisfies the condition This shows that the limit cycle solution appears. In figure  3, the trajectory (α, A) is shown in the case of τ peel S > 1. In this limit-cycle solution, the variable A is re-defined such that A is normalized by the minimum value of the amplitude. All coefficients, e.g. τ peel , etc, are defined by use of the minimum value of the peeling mode amplitude. There is no qualitative difference between the cases of equations (20) and (31).
Considering the results for equations (20) and (31), one can conclude that the dynamical solution of variables surface pressure gradient and peeling mode amplitude can show the limit cycle solution if holds. Therefore, so long as the change of global parameters (associated with the pressure gradient) drives the plasma into unstable region of figure 1, and the time scale separation is allowed, the limit cycle can appear.

Summary and discussion
In this article, a coupled dynamics of the edge pressure gradient and the peeling mode amplitude is studied, by assuming the time scale separation, i.e. the magnetic shear s and the edge current J edge are kept constant in the dynamical change of (g, A). The limit cycle oscillation is possible to occur, in which the amplitude of peeling mode (which is in a quasilinear state) and the pressure gradient oscillate slowly in time.
The dynamical stationary state of the peeling mode is possible to occur, even under the condition that the external operation condition drives the plasma condition into the unstable domain of peeling mode instability. In this dynamical stationary state, the peeling mode repeats the growth and damping periodically. In order to keep the amplitude of the peeling mode at quasistationary value, the nonlinear stabilization effects by higher harmonics of the peeling mode is not necessary. Based on this model the modulation frequency of the peeling mode amplitude is derived. The amplitude modulation of the pressure gradient is also obtained.
Necessary condition for the limit cycle is shown as equations (2) and (34), that is (a) the peeling mode becomes unstable in the higher-pressure-gradient region across the stability boundary condition, and that (b) 0 < τ peel , i.e. the growth of the peeling mode causes the reduction of the pressure gradient near edge.
The condition (a) indicates the main contrast to the result of the pioneering work by Connor [9], in which the peeling mode becomes unstable in the lower-pressure-gradient region across the stability boundary condition (for fixed J edge ). While the possibility of explosive growth of the peeling mode was indicated in [9], the dynamical stationary state of quasilinear mode is demonstrated here. There is a third possibility that the stability boundary is along the constant J edge surface in the parameter space. In this case, the change of α does not modify the growth rate of peeling mode. Then, the limit cycle of quasilinear mode and pressure gradient does not appear: the stationary state, if it is possible, is realized by the nonlinear stabilization by higher harmonics. This classification may explain the contrast of EHO in JT-60U to those in DIII-D where strong higher harmonics were observed.
As an initial application of this model to the observation of QH-mode on JT-60U, examples from [37,38] are compared to the present analysis. A tentative calculation of the eigenfunction seems to support the hypothesis 0 < τ peel [38]. By using the observed amplitude of the peeling mode and the estimated value of the relaxation time, (a) the frequency of the amplitude modulation of peeling mode is estimated as about 200 Hz, and (b) the relative modulation amplitude of α is O(1%). It seems that this hypothesis is promising in understanding the QH-mode sate of JT-60U, so that further detailed test is encouraged. It is noted that the plasma equilibrium is subject to a slow change in studying the QH-mode on JT-60U [37,38]. However, the rate of this change (e.g. being evaluated by the rate (dR 0 /dt)/R 0 , where R 0 is the major radius) is much smaller than the oscillation frequency of the limit cycle. Thus the simplification to neglect the plasma motion is not a bad approximation.
Before closing this article, we add some discussion on the validity of this simple modeling. As is mentioned in the introduction, the roles of inhomogeneous E × B velocity, V ′ E×B , and of ion diamagnetic drift frequency, ω * i , have attracted attentions in calculating stability of the peeling mode: they are neglected in the model equation (1), in which the influence of pressure gradient is kept, in order to keep the transparency of the argument. The extension of this model to include these neglected effects is briefly sketched here. This extension does not change the conclusion qualitatively. When the pressure gradient is perturbed by the amount of δα, the gradient of E × B velocity V ′ E×B and ion diamagnetic frequency ω * i are modified by the mount of δV ′ E×B and δω * i , respectively. As a result of these modifications, the growth rate γ is modified by the mount δγ, as We here assume that the gradients of density, electron and ion temperatures vary coherently with fixed phase in the peeling mode. (In this sense this assumption is an opposite limit of phase-slip model [35]. We note that in the experiment on JT-60U, which is of our interests, the phase-slip was not observed.) Under this circumstance, the variation of perturbation of E × B velocity is governed by the change of pressure gradient, because other sources of E × B velocity (like ion orbit loss, turbulence driven Reynolds stress, etc) are quenched in the H-mode [43]. For small perturbation, one has δα.
A similar argument applies to the modification of ion diamagnetic drift frequency, Substituting equations (36) and (37) into equation (35), one has δα. (38) Near the stability boundary (for the small perturbation δα), equation (38) may be rewritten as Equation (39a) has the same structure as equation (2). Thus, the qualitative feature of the limit cycle oscillation does not change even if one takes into accounts of inhomogeneous E × B velocity, V ′ E×B , and of ion diamagnetic drift frequency, ω * i .
Here we explain assumptions and reasoning to employ unconventional expressions for the energy evolution equations, equations (3)-(6). The energy evolution equation is usually written as where p is a product of number density and temperature, nT, Γ is the particle flux, and P is the heating power density. We neglect the convective heat flux, for the simplicity of argument, because the subject in this article is to study the coupled dynamics between two-variables (pressure gradient and peeling mode amplitude), and we have.
where x is the local Cartesian coordinate to denote the radial direction, because we are studying the thin region near the plasma surface. The radial heat flux is given, with the approximation of diffusion model of heat flux, where χ is the thermal conductivity. The pressure gradient is chosen as a main variable of the analysis, so that one has a relation between heat flux and pressure gradient with modified thermal conductivity K as The parameter K is modified with a numerical coefficient of the order unity with respect to χ. As is explained later, the plasma is almost frozen with the perturbed magnetic field line, so that the coefficient treated here as constant under the modulation of peeling mode amplitude. (In this sense, this model is an opposite limit of the [35], in which the change of relative phase between density and temperature is a key.) Taking into account of the enhancement of thermal diffusivity χ by the peeling mode as the coefficient K is written as where K 0 is the value in the absence of the peeling mode and K peeling is the increment owing to the peeling mode, respectively. Combining equations (A2) and (A4a), and rewriting p in terms of normalized value α = −2q 2 R 0 µ 0 B −2 dp/dr, one has where P norm = − 2q 2 R 0 µ 0 B −2 dP/dr. In changing the variable from p to α in the second term in the RHS of equation (A6), the spatial derivative of normalizing coefficient q 2 R 0 µ 0 B −2 is neglected, based on the assumption that the changes of K and α associated with the modulation change of the peeling mode amplitude is much strongly localized near the surface than other equilibrium parameters. With these assumptions and bases of modeling, the basic equations in the main text, equations (4) and (6) are employed in this article.

Appendix B. Effect of diffusion of pressure perturbation
In equation (8), the diffusion of pressure perturbation (the second term in the RHS) is considered to be much smaller than the effect of steepening by the peeling mode (third term in the RHS), and is neglected in the following analysis. In this appendix, the (small but finite) effect of diffusion of pressure perturbation is discussed.
We take the assumption that the second term in equation (8), ∆ K 0 g, is much smaller than the effect of peeling mode, this pressure-diffusion effect is denoted in this point model as where the minus sign indicates that this effect reduces perturbation of gradient, and the ε denotes a proportionality coefficient parameter. Equation (14) is modified as Combining equations (13) and (B2), and following the same procedure to obtain equation (17) in section 3, one has Comparing equation (B3) to equation (17), one sees that there appears a damping term in the RHS of equation (B3).
Equation (B3) is solved perturbation way. An unperturbed solution is equation (19), If one chooses the initial condition (A = 1, g = 0) at t = 0, one has where abbreviation G 2 = (τ peel γ 0 / α c ) g 2 is used. The additional term (4th term in the RHS) that is proportional to ε is negative definite. Owing to this finite damping term, the oscillation amplitude of pressure-gradient modulation, G 2 , gradually becomes smaller. Because the integrand of the 4th term in the RHS is O (1), and other terms is also O(1), the magnitude of G 2 substantially reduces and approaches to zero, when the elapse time is t ∼ O ( ε −1 τ peel ) . Therefore, the characteristic time, t decay , that the amplitude of limit cycle oscillation decays, is estimated as Oscillation period of limit cycle is estimated as equation (24), τ lc ∼ (τ peel /γ 0 ) 1/2 . Therefore, during the time period of equation (B5), number of oscillation is estimated as