The effect of time-varying flow-shear on the nonlinear stability of the boundary of magnetized toroidal plasmas

We propose a phenomenological yet very general model in a form of generalized complex Ginzburg-Landau equation to understand the dynamics of the quasi-periodic fluid instabilities (called edge-localized modes) in the boundary of toroidal magnetized high-temperature plasmas. The model reproduces key dynamical features of the boundary instabilities observed in the high-confinement state plasmas on the KSTAR tokamak, including quasi-steady states characterized by field-aligned filamentary eigenmodes, transitions between different eigenmodes, and rapid transition to non-modal filamentary structure prior to the relaxation. It is found that the inclusion of time-varying perpendicular sheared flow is crucial for reproducing the observed dynamical features.


I. INTRODUCTION
Relaxation phenomena in magnetized plasmas are widespread in nature. 1,2 A notable example is the explosive flares on the surface of the Sun. Another example is the semi-periodic explosive bursts appearing at the boundary of toroidally-confined high-temperature plasmas (e.g., tokamak). In toroidal magnetic confinement devices, sufficient heating of the plasma can lead to a transition from low-confinement state (L-mode) to high-confinement state (H-mode) if the heating power exceeds a threshold. During the transition, a transport barrier (called pedestal) spontaneously appears at the edge of plasma via strong E × B flow shear which reduces heat and particle transports. However, this barrier is quite unstable and prone to a class of fluid instabilities called edge-localized modes (ELMs) driven by the large gradient of density, temperature, current density, and flow. [3][4][5][6][7][8] It is believed that these instabilities are responsible for the relaxation (or crash) of the transport barrier, i.e., rapid expulsion of heat and particles. The expulsion events are commonly called ELM crash. The H-mode plasmas are characterized by semi-periodic cycles between slow transport barrier buildup and its fast relaxation.
The ELM crash must be controlled because the natural or uncontrolled crashes induce significant heat and particle fluxes which can damage the plasma-facing walls of the confinement device.
Magnetic perturbations have been used successfully to mitigate or suppress the crash 9-11 but the underlying mechanisms of mitigation and suppression are still unclear. Accordingly, it is crucial to understand the dynamics of ELM for more reliable and robust methods to avoid the crash. For this reason, a nonlinear mathematical analysis is required beyond linear stability analyses. 12 For the purpose of studying the nonlinear behavior, a nonlinear model for the perturbed pressure was derived in a form of complex Ginzburg-Landau equation based on a 1D reduced MHD model. 13 The numerical solutions to the model equation showed nonlinear relaxation oscillations with the characteristics of type-III ELM. Inspired by Ref. 13, we mathematically studied the model equation to understand the effect of perpendicular flow shear on the nonlinear behavior of the perturbed pressure during the ELM cycle. 14 More precisely, it was shown that there exists a linearly stable symmetric steady state for small shear and the first eigenvalues of unstable states for the case of zero shear are bounded below by a positive constant. In the case of large shear, a theoretical clue was found for the long-time behavior of the solutions: 1) nonlinear oscillation; 2) convergence to 0. The theoretical results were supported by numerical verifications.
However, in Ref. 14, the shear strength was set constant in time, which was insufficient to explore clues for the various phenomena observed in experiments on the Korea Superconducting Tokamak Advanced Research (KSTAR) device such as quasi-steady state with a single eigenmodelike structure 15 and fast transitions between the quasi-steady states. 16 In this paper, the effect of time-varying flow shear is analyzed as the key for accessing different dynamical states. The remaining of the article is organized as follows. In section II, we present the analysis of the model for the case of a single-mode. In section III, we extend the model to treat the case of two coupled modes. In section IV, we discuss the results and give a conclusion with a prospect for comparison with future numerical simulation including the effect of time-varying flow shear.

II. ANALYSIS OF SINGLE-MODE
We consider the following single-mode equation for the perturbed pressure P(t, x, y) in cylindrical magnetized plasma assuming local slab geometry with the magnetic field direction z, the local radial direction x, and the perpendicular direction y: where W K (x) = tanh (Kx) is the prescribed shear flow with x ∈ [−1, 1], K > 0 is the inverse of the shear layer width, and A ≥ 0 is the shear flow strength. Eq. (1) may be considered an extension of Ginzburg-Landau equation (GLE) with constant complex coefficients. Note that P represents the complex-valued amplitude of a Fourier mode, i.e. δP(x, y, t) = P(x, t)e iky + c.c. Here, γ N , γ L and η are constant coefficients for the nonlinear, the linear growth and the dissipative terms respectively. It was observed that the behavior of a solution to Eq. (1) is completely different with the presence of the flow-shear for both the Dirichlet and Neumann boundary conditions. 13,14 Since it is unclear which boundary condition is reasonable in real experiments, both types of boundary conditions are considered here to understand the long-time behavior of a solution P (t, x) to Eq. (1): Inspired by Ref. 14, we will consider two subjects for the model Eq. (1). The first subject is to characterize the long-time behavior of a solution P(t, x) for the fixed large shear strength A so that we can distinguish the regions of either convergence to 0 or nonlinear oscillations in the γ L -η parameter space. The second subject is to characterize the long-time behavior of a solution P(t, x) between nonlinear oscillations and convergence to nontrivial steady states in the A-K parameter space under suitable fixed parameters γ L and η such that non-trivial solutions are guaranteed. We find a threshold A K > 0 for each K such that solutions converge to a nonzero steady state of Eq. (1) for A < A K and nonlinearly oscillate for A > A K . Combining these results, we propose that the salient features of the ELM dynamics observed in the KSTAR H-mode plasmas can be explained based on time-varying perpendicular shear flow.

A. Long-time behavior of P(t, x) on γ L and η
Notice that the Dirichlet boundary condition does not allow nonzero uniform steady states of Eq. (1) even without the shear in contrast with the Neumann boundary condition. Nevertheless, we obtained similar results for both boundary conditions. Fig. 1 represents the long-time behaviors of a solution P(t, x) on γ L and η for a fixed large A = 50 in both boundary conditions. The blue regions in Fig. 1(a)-(b) display that P(t, x) converges to 0 as t →∞. Conversely, red regions in Fig. 1(a)-(b) display that P (t, x) oscillates nonlinearly in time. These results show a certain relation between η and γ L which determines the long-time behavior of P (t, x). Inspecting Fig. 1, it is clear that nonlinear oscillations are guaranteed only if the ratio γ L /η is sufficiently large. Otherwise, P(t, x) converges to 0. Note that the parameters in Eq. (1) are related to heat flux Q as (see Ref. 13), where Q c is the threshold heat flux related to the critical pressure gradient for linear instability, p 0 is the reference pressure, and a denotes the radius of the cylinder (see Ref. 13 for detail). Therefore, even if the heat flux Q exceeds the linear threshold Q c , nonlinear oscillations may not occur if 0 < Q Q c 1 such that γ L 1 and (γ L /η) 1. This is consistent with experiment observations since it is known that ELM crash does not immediately occur after Q exceeds Q c (see Fig. 1 in Ref. 17).
It is also possible to interpret the case of Q Q c < 0 (γ L < 0) as L-mode. γ L < 0 guarantees the long time behavior of P (t, x) such that lim t→∞ |P (t, x)| → 0. Therefore, Eq. (1) provides a reasonable explanation of the overall ELM dynamics. We need to discuss the effect of γ N . Our expectation is that the stability of the zero solution is crucial to determine the long-time dynamics of P(t, x) for a fixed A 1. In consideration of the analysis result in Ref. 14, it is natural to think that P(t, x) will oscillate nonlinearly if the zero solution is unstable, but converge to 0 if the zero solution is stable. For this prediction, we linearized Eq. (1) around the zero solution P = 0 and proved that the stability of the zero solution is independent of γ N , as expected: Accordingly, it is reasonable to expect that γ N cannot affect the long-time behavior of the zero solution for large A > 0. Conversely, γ N is expected to affect the long-time behavior of the non-zero solution for large A > 0. Under this prediction, we confirmed numerically that γ N does not affect the qualitative long-time behavior of the solutions illustrated in Fig. 1(a)-(b). Instead, γ N can affect the amplitude of nonlinear oscillations. The change of the amplitude |P (t, 0)| in our model is strongly  B. Long-time behavior of P (t, x) on A and K Fig. 2 suggests that there exists a threshold flow shear amplitude A K for given K for both boundary conditions. If 0 < A < A K (blue regions), the solution P(t, x) converges to a nonconstant steady state P s (x) for any given initial condition. On the other hand, the qualitative long-time behavior of P(t, x) abruptly changes if A > A K (red regions). P(t, x) oscillates nonlinearly and never converges to any steady state in the red regions. These numerical results show that there is a certain stability/instability criterion A K of A for each K > 0 for both boundary conditions. According to Fig. 2, we can also predict that ELM crash only occurs under sufficiently strong flow shear. We can also observe that as approaching the threshold line in Fig. 2, the amplitude of nonlinear oscillations (in the red regions) increases and the central value |P(0)| for a nonzero steady state P(x) (in the blue regions) decreases but remain finite (i.e. nonzero). Besides, it is also observed that A K and K are inversely correlated for small K for both boundary conditions, but A K barely changes for large K.
Mathematical clues for the two different dynamic behaviors illustrated in Figs. 1 and 2 can be explained in the case of the Neumann boundary condition. Let where θ = ∂ x θ. In Eq. (3), the shear term AW K (x) affects the amplitude R only indirectly via the phase-gradient θ . Without flow-shear (A = 0), the steady-state P = (γ L /γ N ) 1/2 is the only stable equilibrium. 18 Hence, without flow-shear, the phase-gradient θ converges to 0. However, for finite flow-shear, the term ηRθ 2 in Eq. (3) is nonzero and causes R to decay in time. If the shear is large, the term ηRθ 2 dominates the linear growth term γ L R in a neighborhood of x = 0, so R (t, 0) decays due to the phase-gradient θ until a critical phase-gradient θ = θ c is reached. After decaying, however, the term γ N R 3 is weak close to 0 and the term η∂ 2 x R grows so large that R(t, 0) tends to return to its original state with the help of the linear drive γ L R. This interaction between decay and growth terms makes the nonlinear oscillation. However, if γ L is too small, i.e., the mode is linearly stable, the term η∂ 2 x R is insufficient to fully dominate the term ηRθ 2 . Accordingly, it is impossible to return to the initial state and R(t, x) converges to 0 instead. Similar explanations for the behavior of nonlinear  oscillations were introduced in Refs. 13 and 14. In addition, it can be proved that K is not an important parameter in Fig. 2

C. The effect of time-varying A
Nevertheless, we could not observe non-oscillating quasi-steady state for the prescribed shear flow AW K (x) for both boundary conditions when A > A K . The existence of a quasi-steady state is important for the validation of our model because the ELM dynamics observed on the KSTAR consists of distinctive stages including quasi-steady states, transition phase, and crash phase. 15 We believe that it is impossible to obtain a quasi-steady state for time-independent coefficients. Indeed, if |∂P L /∂t| 1, then a solution should be close to a steady state. However, there is no reasonable  steady state P s A,K (such that ∂ x P s A,K (x) ≤ 0 in 0 ≤ x ≤ 1 and ∂ x P s A,K (x) ≥ 0 in 1 ≤ x ≤ 0) for a sufficiently large fixed shear, 14 so we cannot expect a quasi-steady state. In real experiment, it is natural to think that shear flow evolves, i.e., A and K vary in time. Thus, it makes sense that in the quasi-steady state phase, the parameters A and K are initially located in a region where solutions converge to a steady state (blue regions in Fig. 2), but as time flows, a shear flow gradually increases, and A and K gradually change. As a consequence, as A exceeds the critical point A K , i.e., A moves from the blue regions to the red regions in Fig. 2, the quasi-steady state can no longer exist, which may amount to the sudden crash observed in each ELM cycle. The existence of quasi-steady states with time-varying A(t) is numerically illustrated in Figures 3  and 4 for both boundary conditions. These numerical examples suggest that the change of A induces different stages in the ELM dynamics. Based on these results, we expect that magnetic perturbations can reduce the shear flow strength A such that quasi-steady ELMs can persist without crash, which would correspond to the suppression (absence) of ELM crashes.

III. ANALYSIS OF COUPLED MODES
In this section, we consider two coupled modes with the Neumann boundary condition to study the mode transitions during the quasi-steady observed on the KSTAR. 16 Let W (x) be a prescribed shear flow profile and the pressure P be written as where P = P(t, x) is the slowly time-varying equilibrium pressure and P = P(t, x, y) is the pressure perturbation: with |k 1 | |k 2 |. Extending the single mode model in Ref. 13, we consider the following model: where η > 0, A > 0, b > 0, c > 0, d > 0, C 1 ≥ 0, and C 2 ≥ 0 are constants. With the help of the slaving approximation ∂P ∂t ≈ 0 , 13 we can obtain from Eq. (8) for a constant e ∈ R using ∫ 1 0 | P| 2 dy = |P 1 | 2 + |P 2 | 2 . Therefore, substituting Eq. (9) into Eqs. (6)-(7) yields Denoting γ N bc d , we can rewrite Eqs. (10)-(11) as Let P 1 = R 1 exp (iθ 1 ) and P 1 = R 2 exp (iθ 2 ) . Then Eqs. (12)-(13) can be written aṡ Here, we can interpret γ N , γ L 1 , γ L 2 and η as constant coefficients for the nonlinear term, the linear growth terms for P 1 and P 2 , and the dissipative term respectively. In this paper, we only consider positive values of γ L 1 , γ L 2 , η and γ N . The only difference from Eq. (1) to Eqs. (12)-(13) is the presence of the coupling terms γ N P 1 |P 2 | 2 and γ N P 2 |P 1 | 2 in the equations for P 1 and P 2 respectively, which can account for the mode transition observed in Ref. 16.
In both cases, the mode with higher γ L becomes dominant eventually as expected. However, there is a subtle difference in the time scale between Fig. 5 and Fig. 6. We can explain the difference as follows. For the case of Fig. 5, γ L 1 > γ L 2 and |k 1 | < |k 2 | mean that P 1 has stronger linear growth and, at the same time, less suppression due to the shear compared to P 2 so that P 1 will quickly become dominant. However, in the case of Fig. 6, although γ L 1 < γ L 2 , it takes longer for P 2 to become dominant because P 1 is less suppressed than P 2 by the shear.

B. Long-time behavior for time-varying A
The analysis shown in Figs. 5-6 still cannot explain the transitions between quasi-stable modes observed in experiments. 16 Now, we consider time-varying A in Eqs. (12)- (13) to understand the mode transitions for the case with γ 2 > γ 1 and k 2 > k 1 . P 2 is dominant for sufficiently small A. If A increases in time, it is expected that P 2 is more suppressed than P 1 because k 2 > k 1 means that P 2 is more sensitive to A than P 1 , so P 1 can become dominant finally. Figs. 7-8 show the behaviors of |P 2 (t, 0)| and |P 1 (t, 0)| with growing A, supporting our prediction. Note that |P 2 (t, 0)| is highly oscillating before convergence to 0 in Fig. 8, but not in Fig. 7. We should mention that the numerical examples presented here capture the importance of time-varying A and offer qualitative explanations for various types of mode transitions observed in experiments. FIG. 7. The time behaviors of |P 1 (t, 0)| and |P 2 (t, 0)| for time-dependent A(t) with η = 1, k 1 = 1, k 2 = 3, γ L 1 = 10, γ L 2 = 11, γ N = 1 and W (x) = tanh(25x). We imposed weak initial conditions P 1 (0, x) and P 2 (0, x) as (γ L 1 /γ N ) 1/2 /1000 and (γ L 2 /γ N ) 1/2 /1000 respectively. A increases linearly, reaching the value 6.4393 at the end of the horizontal x-axis in the figure. First, P 2 (t, 0) is dominant and quasi-steady when the shear is small. As the shear increases beyond a critical value, |P 1 (t, 0)| increases rapidly while |P 2 (t, 0)| vanishes rapidly. After that, P 1 (t, 0) remains in a quasi-steady state until it falls to 0 abruptly.

IV. CONCLUSION
In summary, we considered two cases of ELM dynamics based on the extended complex GLE, Eq. (1). In the case of the single-mode, we studied the long-time behavior of the solution with fixed model coefficients and showed that the linear drive γ L and the time-varying shear flow A(t) determine the long time behavior of the solution. If the linear growth term is sufficiently large, the nonlinear oscillations are guaranteed for large shear flow. Conversely, the solution converges to a nonzero steady state for weak shear flow (Fig. 2). The long-time behavior for the small linear growth term is interesting because a solution converges to 0 for large flow shear (Fig. 1). Combining these results, we conclude that it is insufficient to consider the fixed coefficients on time to realize the quasi-steady states which are observed in experiments. 15 Therefore, by imposing time-varying shear flow, we obtained quasi-steady states numerically (Figs. [3][4]. Note that the plasma rotation (flow) and its shear are well recognized as important factors for the stabilization of ELMs in MHD simulations, in particular for high toroidal mode numbers. [19][20][21][22] In recent nonlinear simulations using the nonlinear resistive MHD code JOREK for KASTAR H-mode plasmas, 22 it was shown that the shear of the self-consistent poloidal rotation profile becomes stronger toward the onset of ELM crash for the case of a single harmonic (eigenmode) simulation with the measured toroidal rotation profile.
Based on our numerical analysis, we expect that the quasi-steady mode can persist if the shear flow is reduced below the critical threshold by application of external magnetic perturbations, which may provide a candidate mechanism for the non-bursting quasi-steady modes in the ELM crash suppression experiment. 11 To study the dynamics of coupled modes P 1 and P 2 , we derived equations (12)-(13). We confirmed that the linear growth terms are crucial to determine the long-time behavior of P 1 and P 2 (Figs. [5][6]. Inspired by these results, we considered the increasing A(t) on time and showed that rapid mode transition occurs (Figs.7-8), reproducing qualitatively the observed mode transitions in experiments. 16 Although we dealt with the equations (12)-(13) for coupled modes, it is also possible to obtain equations for more than two modes and show that each mode solution is successively dominant with suitable time-dependent A(t).
Note that our extended complex GLE model does not resolve the trigger problem for the ELM crash because the model does not directly show the existence of non-modal solitary perturbation and its burst (which initiates the pedestal collapse) observed in the experiments. 15,23,24 Nevertheless, the rapid disappearance of eigenmodes in our model may allude to a condition for the rapid unidirectional transition to the non-modal perturbation. In this regard, our model may be considered complementary to the existing trigger theory 25,26 and the nonlinear multi-harmonic simulations where the energy of the linearly dominant mode is redistributed into a broad spectrum as approaching to the ELM crash. 22,27,28 To conclude, it is critical to consider the time-varying A for explanation of dynamic features in ELM phenomena using the given models (1) and (12)-(13) for single and coupled-modes, respectively.