Kinetic study of the bifurcation of resonant magnetic perturbations for edge localized mode suppression in ASDEX Upgrade

The correlation between the bifurcation of resonant magnetic perturbations (RMPs) to the unshielded state and edge localized mode (ELM) suppression in ASDEX Upgrade is studied using a kinetic plasma response model numerically and analytically. For the numerical studies, the linear kinetic Maxwell solver KiLCA for cylindrical geometry and the quasilinear transport code QL-Balance are used in combination with the ideal MHD solver GPEC to account for realistic tokamak geometry. Based on this modelling, a numerical local bifurcation criterion is introduced which estimates the effect of RMP-induced temperature plateau formation in the resonant layer. Its analytical form is derived in constant-psi approximation. The kinetic model reproduces the known gyrocenter resonance, Er=0 , and the electron fluid resonance. In contrast to MHD theory, the latter is located at the zero of the perpendicular electron fluid velocity computed only with half of the electron temperature gradient. The application of the criterion to experimental data shows a correlation between bifurcation and the ELM suppression phase. Moreover, an electron density limit is found resembling the one observed in experiments.

The correlation between the bifurcation of resonant magnetic perturbations (RMPs) to the unshielded state and edge localized mode (ELM) suppression in ASDEX Upgrade is studied using a kinetic plasma response model numerically and analytically. For the numerical studies, the linear kinetic Maxwell solver KiLCA for cylindrical geometry and the quasilinear transport code QL-Balance are used in combination with the ideal MHD solver GPEC to account for realistic tokamak geometry. Based on this modelling, a numerical local bifurcation criterion is introduced which estimates the effect of RMP-induced temperature plateau formation in the resonant layer. Its analytical form is derived in constant-psi approximation. The kinetic model reproduces the known gyrocenter resonance, E r = 0, and the electron fluid resonance. In contrast to MHD theory, the latter is located at the zero of the perpendicular electron fluid velocity computed only with half of the electron temperature gradient. The application of the criterion to experimental data shows a correlation between bifurcation and the ELM suppression phase. Moreover, an electron density limit is found resembling the one observed in experiments. 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.

Introduction
The operation of tokamaks in high confinement mode or H-mode is usually accompanied by edge localized modes (ELMs) [1]. These instabilities are characterized by large transient heat loads to plasma-facing components [2]. Since these heat loads are projected to be unacceptably high in bigger devices, e.g. ITER [3], methods to suppress ELMs are required. ELMs occur if the plasma pressure gradient at the edge exceeds a limiting value that is caused by the formation of a transport barrier [2]. Thus, suppressing ELMs in H-mode plasmas requires lowering the edge pressure gradient.
As experimentally demonstrated in several tokamak devices [4][5][6][7][8] and observed as a density 'pump-out', lowering the plasma-edge pressure gradient can be achieved by applying small (δB/B ≈ 10 −4 ) externally generated resonant magnetic field perturbations (RMPs). RMPs distort the axisymmetry of the equilibrium flux surfaces, and, for high enough amplitude, cause the non-ideal plasma to reconnect magnetic field lines at rational surfaces and thus create magnetic island chains [9]. A rotating plasma, however, responds to RMPs by generating localized shielding currents at respective rational surfaces [10][11][12]. Nevertheless, if the amplitude of an RMP mode exceeds a threshold value, the mode can penetrate and a bifurcation to a 'reconnected' state occurs [13,14]. Bifurcation of an RMP mode resonant at the pedestal top was shown to be involved in ELM suppression [15][16][17]. However, the specific conditions necessary for bifurcation and, by extension, ELM suppression are not fully understood yet [18][19][20]. The purpose of this paper is to advance our understanding of bifurcation and its connection to ELM suppression.
As recently reported in [21], we study bifurcation and ELM suppression with an extended linear and quasilinear kinetic model [22]. This model is composed of the cylindrical Maxwell solver KiLCA [23] and the 1D radial transport code QL-Balance [22]. The former determines the linear kinetic plasma response to RMPs. The latter evolves the plasma parameters in time by considering RMP-induced and anomalous fluxes. The model includes a qualitative estimation of poloidal mode coupling by rescaling the cylindrical plasma response currents determined by KiLCA with the toroidal plasma response current of the linear ideal MHD code GPEC [24]. Within this model, we identify bifurcation with a local criterion [21] that can approximately be written in an analytical form using the constant-ψ approximation. Building the bridge to the experiment, we apply the bifurcation criteria to data from ASDEX Upgrade and assess their correlation to ELM suppression.
The content of this paper is ordered as follows. We give an overview of the model and discuss the linear and quasilinear plasma response as well as the role of fluid resonances in section 2.1. In section 2.2 we explain the used ASDEX Upgrade data and its preprocessing. Then, in section 3, we investigate a specific bifurcating case and the accompanying effects. In section 4, we discuss the approximate local bifurcation criterion that we apply to an ASDEX Upgrade discharge. We investigate the scaling of the bifurcation criterion with electron temperature and density and determine an analytical expression for the criterion in section 4.2. Finally, in section 5, we summarize and discuss the important findings.

Model
The model employed here is based on previous work [11,22,23,25]. In this section, a brief formulation of the RMP penetration model developed earlier is presented together with an improved method to couple our model with realistic tokamak geometry. Besides that, general expressions derived earlier are studied here in more detail, resulting in approximate analytical expressions for the penetration factor and for quasilinear transport coefficients used later in the formulation of the analytical scaling of the RMP penetration threshold.

Linear plasma response, non-ideal parallel current.
To evaluate the non-axisymmetric electromagnetic field perturbation induced by the RMP coils, we use the kinetic Maxwell solver KiLCA [23]. In this code, the total electric and magnetic fields are represented as a sum of equilibrium and perturbation field, E = E 0 + δE and B = B 0 + δB. The time dependence of the perturbation field is assumed harmonic, i.e. {δB, δE} = Re({ B, E} e iωt ), with the perturbation frequency ω. The time harmonics of the perturbation fields are then computed from the set of Maxwell equations, where c is the speed of light. The current densities on the right hand side of Ampere's equation are j RMP , the RMP coil current density located outside the plasma, and j, the linear plasma response current density computed in the kinetic approximation. Maxwell's equations are solved for a simplified tokamak geometry in the form of a radially inhomogeneous straight plasma cylinder with identical ends (infinite aspect ratio limit of a circular tokamak). Introducing cylindrical coordinates (r, θ, z), respectively being the minor radius, the poloidal angle and the axial position (toroidal coordinate z = R 0 φ where R 0 is the major radius and φ is the toroidal angle), the cylindrical components of the electromagnetic field and the currents are Fourier-expanded as a = m,n a m (r)e i(mθ+kzz) .
Here, a ∈ { E i , B i , j i }, k z = n/R 0 , and m ≡ (m, n) are poloidal and toroidal mode numbers, respectively. This Fourier expansion reduces the set of three dimensional partial differential equations to a set of one dimensional ordinary differential equations (ODEs) which is solved in a moving reference frame where the perturbation frequency is finite. The moving reference frame is chosen since we consider static RMPs, while the plasma response current is expressed in KiLCA solely via the perturbation electric field. Note that the Fourier expansion in cylindrical geometry neglects the coupling of poloidal modes that is present in real tokamak geometry [26]. This geometrical issue is discussed in section 2.1.4. The plasma response current density j which is generally determined by an integral conductivity operator is modelled in KiLCA in the form of a finite Larmor radius (FLR) expansion. The conductivity kernel is computed from the solution of the linearized kinetic equation, including an energy preserving Fokker-Planck type collision operator (Ornstein-Uhlenbeck operator with energy conserving term, for details see [23]). In this case, the contravariant components of the plasma response current density for a single species can be expressed via the covariant components of the perturbation electric field as [23] In this constitutive relation, the sum goes up to the desired order of the FLR expansion N FLR , σ il kk ′ ,α is the plasma conductivity kernel and α ∈ {e, i} is the species index. In the following, the FLR expansion is used to first order which is sufficient to reproduce the ideal plasma response described by a single second order ODE of Furth et al [27]. This is valid in the majority of the plasma volume, except for the resonant layers localized around rational magnetic flux surfaces (see [11]). Respective 'ideal' response currents, i.e. diamagnetic and Pfirsch-Schlüter current, are contained in the first order FLR expansion terms in (3). Further, the dominant, zero-order contribution to the parallel current density responsible for the shielding of RMPs is, in the laboratory reference frame, well represented by the drift kinetic expression, equation (60) of [22]. This drift kinetic expression is given as Here, e α , m α , n α , T α , v Tα = (T α /m α ) 1/2 , and ν α are α species charge, mass, unperturbed density, temperature, thermal velocity, and the transverse collision frequency [28], respectively. Thermodynamic forces are given by where E 0r = −∂Φ 0 /∂r is the equilibrium radial electric field. Further, are the tangential component of the perpendicular perturbation electric field and the normal component of the perturbation magnetic field, respectively. These are defined with respect to the unperturbed magnetic flux surface, which is labelled with the radius r. Here, h = B 0 /B 0 is the unit vector along the unperturbed magnetic field. The remaining quantities in (4) are complex susceptibility functions I kl ≡ I kl (x 1 , x 2 ) defined by equations (56)-(59) of [22] as where I kl D is given by equations (A1) and (A2). The arguments of I kl = I kl (x 1 , x 2 ), are the normalized distance to the resonant surface and the normalized radial electric field which also plays the role of collisionality parameter, respectively. In these expressions, and is, up to the minus sign, the perturbation frequency in the moving reference frame where the equilibrium electric field is zero. As shown in [22], the zeroth order parallel current is absent if the perturbed electrostatic potential is constant within the perturbed magnetic flux surfaces (see equation (47) there). Using for the susceptibility functions the relations (A4) derived in appendix A, the parallel current density (4) can be expressed via a single component of the perturbed electrostatic field as Here, E MA m⊥ is the Fourier amplitude of the tangential part of the total perpendicular electric field with respect to the perturbed flux surfaces. For infinitesimal perturbations, it is given as m⊥ , Here, r [ψ] = r + r is the perturbed magnetic flux surface label satisfying the magnetic differential equation B · ∇r [ψ] = B 0 · Figure 1. The magnitude of the electron (top) and the ion (bottom) parallel current density in kinetic (equation (3)) and drift kinetic (equation (4) or equation (11)) approximation for the parameters of section 2.2. The position of the resonant surface is indicated by the vertical dashed line.
∇ r + B r = 0. The field E MA m⊥ results from the misalignment of equipotential and magnetic flux surfaces, and is absent in ideal MHD where these surfaces are perfectly aligned. Within linear theory, perfect alignment is only possible with complete RMP shielding at the resonant surface since the part of the perpendicular field resulting from the corrugation of flux surfaces, E [ψ] m⊥ , is singular at this surface if B r m is finite there. This singularity, however, is an artifact of the definition of E MA m⊥ and is absent in the original expression (4).
The parallel current density from the first order FLR expansion (3) is compared to the drift kinetic approximation equation (4) (or equation (11)) in figure 1, where the electromagnetic field is computed by KiLCA using (3). We see that for electrons, the drift kinetic model approximates the kinetic current density well. The ideal MHD contribution, i.e. the Pfirsch-Schlüter current given by the first order terms of the FLR expansion, is seen in the kinetic ion current outside of the resonant layer. This contribution is of the same order for electrons and for ions. However, due to the electron current density being two orders of magnitude higher, the ideal MHD contribution is not visible for the electrons. Thus, ions play no direct role in the shielding of RMPs but contribute only indirectly via their effect on the electrostatic field E ⊥ and on the misalignment field E MA m⊥ . The latter, together with its two components, equation (12), is shown in figure 2. We see that E MA m⊥ rapidly drops to zero in the ideal region outside the resonant layer due to the mutual cancellation of the perturbation electrostatic field E m⊥ and the flux surface corrugation induced field E [ψ] m⊥ (see also figures 6-8 in [22]). The latter actually dominates in most of the resonant layer, especially close to the resonant surface. This means that one can actually ignore E m⊥ there and use, for the drift kinetic current, only the terms in the second line of equation (4) describing the current due to parallel streaming along the perturbed magnetic field. The magnitude of the misalignment field and its components, equation (12), for the same parameters as in figure 1. The resonant surface is indicated by the dashed line.

RMP shielding and fluid resonances.
In order to understand the relation between the parallel electron and ion response currents and to analyse the so-called 'fluid resonances', which are well known from two-fluid MHD theory (see, e.g. [29,30]), let us study equation (11) in more detail. First, the electron fluid resonance, is seen in (11) in absence of an electron temperature gradient, A α 2 = 0. In this case, the parallel current density scales with the combination of thermodynamic forces A α 1 + A α 2 = A α 1 and is proportional to the perpendicular fluid velocity V α⊥ . In general, the latter is given by where V αd is the diamagnetic velocity and D Bα = cT α /(e α B 0 ) is the Bohm diffusion coefficient. Considering electrons and the case A e 1 = A e 2 = 0, for which V e⊥ = 0 at the resonant surface, the dominant electron response current (11) practically vanishes, and thus shielding is lost for the respective RMP mode.
The trend j α m∥ ∝ V α⊥ remains valid close to the resonant surface also in the case of a finite temperature gradient A α 2 ̸ = 0. This can be seen from the behaviour of the susceptibility functions near x 1 = 0 in figure 3. As shown in appendix A, close to the resonant surface the susceptibility functions in question are related as I 31 = 3I 11 (see equations (A15) and (A25)) and, via equation (A4), as I 21 = 3I 10 . These relations indicate that the current density (11) scales linearly with A α 1 + 5A α 2 /2 and, respectively, with V α⊥ , even in the case of a non-vanishing T α gradient. Since in the region near the resonant surface E MA m⊥ → E [ψ] m⊥ we can use there equation (4) with E m⊥ = 0. Further, substituting in equation (4) the limit (A15) and using (13) and (8) we arrive at In fluid theory, the same result follows from the projection parallel to B 0 of the stationary α-species momentum Equation, where the parallel friction force density is considered in the form R α∥ = −m α n α ν α V α∥ and omitting the diamagnetic term in the convective derivative of the parallel fluid velocity perturbation. The latter term vanishes anyway if the Braginskii [31] viscosity η 4 is taken into account ('gyro-viscous cancellation', see e.g. [32]), which is naturally the case in kinetic theory. For the electrons, equation (14) reduces at high collisionality to j e m∥ = −σ ∥ V e⊥ B r m /c where σ ∥ = e 2 e n e /(m e ν e ) is the parallel conductivity (see equation (2.11) of [33]), i.e. j m∥ ≈ j e m∥ results from the dynamo effect of the electron fluid moving in the perturbed magnetic field.
Note that equation (14) follows from the kinetic model of KiLCA for both, electrons and ions, because the collision operator in the model was purposely chosen such that the energy is conserved but not the momentum. This is actually a proper approximation for electrons which collide with slow (practically immobile) and heavy ions while for the ions a rigorous treatment of the straight cylinder (or slab) magnetic field geometry requires momentum conservation. Such a conservation would lead to the vanishing of ν i in the ion species equation (14) and, respectively, to the increased role of (anomalous) shear viscosity at small ω E . However, the ion momentum conservation has been ignored in KiLCA on intention (it can be taken into account as an option, see [34]), in order to model the neoclassical parallel momentum loss rate in the real tokamak geometry in banana and plateau regimes typical for the present day devices and future reactors. (This effect leading to 'poloidal friction' is accounted for e.g. in the reduced-MHD JOREK fluid model [35] using a heuristic approach.) Although such a collision model exaggerates this loss rate by an aspect ratio factor, this has little effect on the shielding as ions are practically excluded from j m∥ . This is due to ions being mostly collisionless (ν i ≪ ω E ), while near the point where ω E = 0 ions are collisional, the aspect ratio still cannot compete with the (m i /m e ) 1/2 ratio between currents (14). In fact, the factor 200 > (m i /m e ) 1/2 between electron and ion current densities we see in figure 1 is actually due to the collisionless ions.
The scaling of the parallel current with V e⊥ does not mean that V e⊥ = 0 corresponds to the maximum penetration point, since this point is essentially determined by the integral parallel current (see section 2.1.4) defined as Here, r m is the resonant surface radius, and δ m ≪ ∆r ≪ L r , where δ m and L r are the typical resonant layer width and the radial scale length of plasma parameter profiles. Contributions of various species to this current are well approximated by the drift kinetic expression (11), where we ignored the variation of plasma parameters and B 0 within the resonant layer and replaced the radial integration with an integration over x 1 using (8) and, approximating k ∥ by a linear function of radius, Here, q and s are the safety factor and the magnetic shear parameter, respectively. For an estimation of integral (16) we again ignore the electrostatic field E m⊥ and use for B r m the so-called 'constant-ψ' approximation, i.e. we assume it constant within the layer. Thus, we arrive at where r Dα = 4π n α e 2 α /T α −1/2 is the Debye radius of species α, I kl int are integrated susceptibility functions defined in equation (A13) and is the temperature gradient contribution to its diamagnetic velocity. Using equations (A24) and (8) we can cast (18) into In this relation, we identify two types of resonances where the response current vanishes. The first one is the so-called 'fluid resonance', where is the resonance condition. The origin of this resonance is the same as for the vanishing of all other equilibrium currents (perpendicular, Pfirsch-Schlüter, bootstrap). All of them can be presented in the form of a linear combination of thermodynamic forces A α k (equation (5)) allowing for a zero. Very close to the resonant surface, the parallel response current density scales with A α k similarly with the perpendicular current and, respectively, with the perpendicular fluid velocity, i.e. they are proportional to each other via equation (14). This, however, does not need to be the case for the whole resonant layer and, respectively, for the total current (20), whose zero is generally shifted from the actual fluid resonance point V α⊥ = 0 in the presence of a temperature gradient [22]. Nevertheless, we will use the well established term 'electron fluid resonance' for this zero-point of the electron response current since it results in the loss of shielding of the resonant field B r m due to the overall dominance of electrons.
The other resonance we find in equation (20) is the electric or gyrocenter resonance ω E = 0. For this resonance, both the electron and the ion current vanish simultaneously. Remarkably, the magnitude of the current density (14) is even slightly increasing when approaching this point. However, the integral current (20) is decreasing because of a decreasing width of the resonant layer. The gyrocenter resonance becomes rather narrow at low collisionalities since it does not exist in the collisionless limit ω E /ν α → ∞, where the decreasing width of the resonant layer δ m ∼ ω E (k ′ ∥ v Tα ) −1 (see equations (A12) and (17)) is balanced by an increasing current density (14) leading to the unchanged integral current (20). In the collisional limit, on the other hand, the resonant layer width scales as The ratio of the radial magnetic field component at the resonant surface computed by KiLCA in plasma, B tot r = |B r m (r m )|, to the value computed in vacuum, B vac r , is shown in figure 4 as a function of the electric drift velocity V E×B (r m ), equation (10). There, the location and the magnitude of the resonances can be seen. In these scans of a perturbation mode m = (6, 2), the profile of the equilibrium radial electric field E 0r was varied by scaling the toroidal rotation velocity profile, which determines E 0r via the ion radial force balance (26). Here, we use parabolic model profiles of the form y(r) = y(r m )F(r), where y(r m ) is the value the profile has at the resonant surface, F(r) = (r 2 − a 2 )/(r 2 m − a 2 ) and a is the plasma radius. The values at the resonant surface for the density and electron temperature were chosen as n e (r m ) = 5 · 10 12 cm −3 and T e (r m ) = 15 keV. There are two cases shown in figure 4. First, the case of constant temperature, T e (r) = T e (r m ) and n e (r) = n e (r m )F(r), and second, the case of constant density, n e (r) = n e (r m ) and T e (r) = T e (r m )F(r). In both cases, the profile of the electron fluid velocity V e⊥ is the same. However, the fluid velocity V e⊥ is equal to the 'resonant' velocity V e res , equation (21), only in the first case (figure 4(a)) since in this case V e dT = 0 and, thus, we see the fluid resonance peak in the penetration factor at V E×B = −V ed . Note that V ed is taken here in the negative direction by convention. Due to the finite electron temperature gradient in the second case (figure 4(b)), we see the fluid resonance peak (or 'kinetic' resonance peak since the peak is shifted by kinetic effects) at V E×B = − 1 2 V ed , which is in agreement with the resonant velocity V e res . Here, the different collisionalities were obtained by formally scaling the Coulomb logarithm leaving the plasma parameters unchanged.
In both cases, we see that the electron fluid resonance is located precisely at zero of V e res at low collisionality while at high collisionality it is shifted towards higher V E×B values. The latter results from a violation of the constant-ψ approximation which is more significant in the high collisionality case where the layer width δ m is ten times larger than at low collisionality. Due to the increase of B r m (r) with radius, the mid-point of the resonant response current (11) shifts from the resonant surface outwards which requires a respective outward shift of the V e res = 0 point for minimization of this current.
Note that the loss of shielding is incomplete in the presence of a finite electron temperature gradient resulting from a perturbation electric field and non-constant-ψ effects ignored in the derivation of equation (20). The latter has a significant effect on the bifurcation criterion discussed below. Besides that, the trend of the gyrocenter resonance to vanish at low collisionality is also observed.

2.1.
3. Quasilinear transport model. As discussed above, the penetration of RMPs essentially depends on the equilibrium plasma parameter profiles. Those, however, are affected by the RMPs. Within quasilinear theory, the slow evolution of the unperturbed plasma parameters, namely, the plasma particle density n e = n i , the toroidal rotation frequency V φ , the electron and the ion temperature profiles T e and T i , respectively, is described by a set of 1D transport equations (see equations (63)-(66) of [22]), ∂ ∂t

∂ ∂t
where ⟨. . . ⟩ denotes the flux surface average and S the flux surface area. The flux surface label r is fixed by the condition ⟨|∇r|⟩ = 1. Further, √ g is the metric determinant of the flux coordinates (r, ϑ, φ), B ϑ 0 and B 0φ are the contra-variant poloidal and covariant toroidal components of the equilibrium magnetic field, and g φφ = R 2 0 is the covariant toroidal component of the metric tensor (the same in any coordinates with rotational symmetry). The quantities S 0 n and S 0 w,α denote the equilibrium sources of particles and energy, respectively, whereas T 0 φ is the equilibrium toroidal torque density. These three equilibrium quantities are determined in a way that, without the perturbation field, the transport equations result in a steady state for the experimental profiles of density, rotation frequency and temperatures. The equilibrium radial electric field E 0r is linked with the toroidal rotation V φ via the ion radial force balance equation (see, e.g. [22,36]), where the coefficient k is computed with the code NEO-2 [36,37]. Note that only the solid body rotation is included in equation (23) because the modification of the ion temperature profile by the RMPs is weak. The effect of anomalous transport is modelled with the anomalous particle Γ A α and heat fluxes Q A α , given by [11,22] and the anomalous shear viscosity µ A ⊥ is approximated as µ A ⊥ ≈ D a . The anomalous diffusion coefficient D a accounts for the net-velocity-independent perpendicular diffusion due to turbulence in a reference frame where the equilibrium radial electric field is zero [22]. In this simplified model, anomalous transport is ambipolar and thus does not produce any torque that has to be included in the momentum equation (23).
In equation (25), we included the main effect of neoclassical physics, which is a contribution to the ion heat flux given by [38] where ω ci is the ion cyclotron frequency and τ i is the ion collision time containing the Coulomb logarithm Λ. The transport resulting from the electromagnetic perturbations induced by the RMPs is described here within quasilinear theory [22]. The particle and heat flux density, Γ EM α and Q EM α , respectively, are linked to the thermodynamic forces (5) by the Onsager-symmetric matrix of transport coefficients, D ql α,ij , as Within quasilinear theory, these transport coefficients have a quadratic dependence on the electromagnetic field induced by the RMP coils and are approximated by the straight cylinder geometry expressions, given by equations (49) of [22] as positive definite quadratic forms in E m⊥ and B r m , Equation (6). Similar to the parallel current density (11), the coefficients can be simplified by using the relations (A6) of appendix A to Thus, we find that the perpendicular electric field perturbation E MA m⊥ caused by the misalignment of equipotential and magnetic flux surfaces is solely responsible not only for the shielding current but also for the quasilinear transport. Considering the quasilinear diffusion coefficients for a single RMP mode m = (m, n) in the limit x 1 = 0, i.e. close to the respective resonant surface, we can simplify equation (30) by using the limiting case expressions (A26) to which holds for arbitrary collisionality. Here, we considered only the electric field perturbation generated from the corrugation of flux surfaces, E [ψ] m⊥ , since it provides the main contribution.

Account of realistic geometry.
The Fourier expansion given in (2) treats toroidal and poloidal modes separately since they are not coupled in cylindrical geometry. However, poloidal mode coupling in realistic, toroidal tokamak geometry may result in a considerable change of amplitudes of the magnetic field perturbation compared to the poloidal spectrum of the external field [39,40]. In the previous model [22], realistic tokamak geometry was accounted for by rescaling the Fourier amplitude of each electromagnetic field harmonic with the ratio of the perturbation field vector potentials from the two vacuum models, C m = A tor θm /A cyl θ,m , where A tor θm is the Fourier amplitude of the poloidal component of RMP coil vector potential computed in straight field line flux coordinates in realistic tokamak geometry and A cyl θ,m is the respective amplitude of the KiLCA vacuum solution in cylindrical geometry. Accordingly, squares of these factors were used for rescaling the contributions of different modes computed by KiLCA to quasilinear diffusion coefficients (30). However, poloidal mode coupling in plasma generally differs from coupling in vacuum leading, in particular, to a significant phase shift of the resonant field components upon phasing the RMP coil currents [4]. To account for this effect, we upgrade the model by redefining the rescaling coefficients as where I ∥,m denotes the integral of the Fourier amplitude of the resonant plasma response current density for the mode m = (m, n). The integral is performed over the poloidal crosssection, φ = const. For the denominator, I cyl ∥m , we take the integral parallel current (15) determined in the cylinder. In particular, we integrate j m∥ computed by KiLCA in a thin annulus around the resonant surface estimating the width of this resonant layer via fitting a Gaussian function to the parallel current density and taking five times the standard deviation. The numerator, I tor ∥m , takes the same quantity but calculated in flux coordinates for the realistic tokamak geometry. For this purpose, we use the generalized perturbed equilibrium code GPEC [24,39] that calculates the ideal MHD plasma response in toroidal geometry.
The rescaling with the coefficient (32) follows from the similarity of the solutions from ideal MHD and kinetic theory in the case of strong RMP shielding. The ideal MHD solution of GPEC assumes the intact (in the leading order, see below) embedding of flux surfaces in the presence of RMPs, which requires that all resonant components of the perturbation magnetic field, (B r /B φ 0 ) m , computed in symmetry flux (PEST) coordinates, are zero at respective rational magnetic surfaces. This is achieved by the presence of current sheets at each resonant surface. In turn, the current sheets result in jumps of the components of respective Fourier amplitudes of the perturbation field that are tangential to the flux surface. The integrated currents of these sheets, I tor ∥m , are directly related to the jumps of the tangential field (and respective jumps of radial derivatives of (B r /B φ 0 ) m ), which are required to annihilate the resonant components. Since the jumps are fully determined by the 'outer' solution of the ideal MHD equations in the plasma volume excluding rational surfaces the currents I tor ∥m are so as well. In particular, in the straight cylinder geometry, jumps of the physical poloidal RMP compon- where k 2 ⊥ = m 2 /r 2 m + k 2 z follows from (9) at the resonant surface r = r m and h z = B 0z /B 0 . For an estimate ignoring the plasma response in the 'outer' solution and assuming a large aspect ratio (m ≫ k z r m and h z ≈ 1), one can relate the jump in the radial derivative to the vacuum field at the rational surface by and estimate the integral current required to shield the vacuum perturbation as In the non-ideal case, i.e. for visco-resistive MHD or kinetic models, where the width and the amplitude of the shielding current density are finite, the integral current I ∥m and the allow us to perform quantitative estimations for toroidal geometry even though our model is formulated in a straight cylinder.
It should be noted that in the case of finite RMP amplitudes, the embedding of flux surfaces is destroyed in higher order of the perturbation amplitude by the formation of residual islands even in ideal MHD with current sheets fully eliminating the normal magnetic field component at the resonant surface [41,42]. The radial size of these islands, however, scales linearly with RMP amplitude and, therefore, is negligibly small compared to the unshielded island width which scales as a square root of the amplitude. Another issue occurring at finite RMP amplitudes is a 'breakdown' of linear theory [43]. Formally, the breakdown occurs as an overlap of perturbed flux surfaces r [ψ] = const computed within linear theory (see the text below equation (12)). A criterion for this overlap is given by |∂ r/∂r| > 1, i.e. if the modulus of the radial derivative of the radial displacement r exceeds one elsewhere in the plasma volume an overlap is signified. However, such an event does not necessarily mean that linear MHD results cannot be used as an order of magnitude estimate. Note that the absence of such events also does not mean that linear MHD is strictly valid. In cases of [43] such a 'breakdown' happens in a limited part of the volume or does not happen at all which indicates that nonlinear terms which prevent the overlap of perturbed flux surfaces and which are ignored in linear MHD model are of order one, and thus the results are still valid as an estimate. In our case such a derivative is much smaller than one, i.e. linear approximation of GPEC is well applicable.

ASDEX Upgrade data
The experimental data used for the analysis within this paper was taken from the ASDEX Upgrade shot 33353 which aimed for the optimization of ELM suppression. In this shot, the toroidal field was B t = −1.8 T, the plasma current I p = 0.9 MA with a corresponding edge safety factor q 95 = 3.65, normalized beta β N = 1.9, elongation κ = 1.65, and upper and lower triangulartiy δ u = 0.23 and δ l = 0.42. The equilibrium data is of type EQH computed by the code CLISTE at ASDEX Upgrade [44]. Note that this equilibrium is not constrained by kinetic profiles. From the equilibrium data, we calculate the toroidal magnetic field B 0z and the safety factor q. We fitted the plasma profiles via the tool augped by free-knot-spline fitting, with the poloidal flux coordinate ρ pol = √ ψ n as the radial variable, where ψ n is the normalized poloidal flux. The effective radius r is then defined via the toroidal flux ψ tor and the toroidal field at the magnetic axis B ref by where the assumption of a circular toroidal cross-section area was made. The profiles include the particle density n e , toroidal rotation frequency V φ , electron temperature T e and ion temperature T i . As an example, profiles for a single time point can be seen in figure 6. The data for the electron temperature and density profiles is from the core and edge Thomson diagnostics [45]. The data for the T i and the V φ profiles is from the charge exchange recombination spectroscopy diagnostics [46]. To fit a profile, we used data from a symmetric time interval (100 ms) centered around the time point (time slice) of interest.
To use the plasma profiles in KiLCA, they need to be extended beyond the separatrix (r max ≈ 62 cm, varying with equilibrium) and the RMP coil (r coil = 67 cm) up to the position of the ideally conducting wall (r W = 70 cm). This is done by applying a smoothed step function, where T orig e is the original profile with polynomial extrapolation to r > r max and T ∞ e = const is the value at r → ∞. We chose parameters d = 0.1 and r c = r max + 0.2 cm resulting in the profiles shown in figure 6. We also extend the safety factor profile beyond the separatrix, though this case is more subtle, since the safety factor diverges at the separatrix. Substituting in the relation between the physical equilibrium field components, qR 0B0ϑ = rB 0z , the safety factor profile from the 2D equilibrium and using the radial force balance, the magnetic field strength profile B 0 (r) is obtained from the ODE where p is the plasma pressure. The solution to (39) fully determines the components of the equilibrium field and, via Ampere's law, the components of the equilibrium current density. The later components are extended outwards using the function (38) and are, in turn, integrated within Ampere's law to get the new magnetic field. This is then used to calculate the modified safety factor profile which coincides with the original profile up to q ≈ 4.5 so that all relevant resonances are unaltered, as shown in figure 7.
The rational mode numbers that we consider are all for toroidal mode number n = 2. This is the one dominating the coil spectrum in the ASDEX Upgrade shot we study. The poloidal mode numbers that we examine are m = 5, 6 and 7 which correspond to rational surfaces further inside the plasma, close to the pedestal top and inside the gradient region, respectively.
Finally, profiles of the poloidal rotation velocity V pol and of the radial electric field E 0r follow from (26). In the straight cylinder geometry with √ g = rR 0 , B 0z = B 0φ /R 0 , V pol = rV ϑ , Figure 7. Plot of the reconstructed profiles of the poloidal rotation velocity V pol , the radial electric field E 0r and the safety factor q for shot 33353 at time t = 2.9 s. The poloidal rotation and the radial electric field were reconstructed with the code NEO-2 [36,37] using the extended profiles. Compared to the original safety factor profile from CLISTE, the modified safety factor profile does not diverge at the separatrix (dash-dotted line) and has a discontinuous derivative due to the extension of the profiles beyond the separatrix.
An example of such reconstructed profiles is given in figure 7.
Since V pol and E 0r are calculated after the extension of the required profiles they are already defined up to the wall (r W = 70 cm). It should be noted that, due to the rescaling procedure of section 2.1.4, the details of the profile extension have no significant effect on the final results as long as the resonant layers for the Fourier modes of interest are well separated from the profile extension region. It should also be noted that the negative electric field at the separatrix and beyond results from using the neoclassical reconstruction in the whole radial range of the main plasma volume and the subsequent extension outside the separatrix. Such a reconstruction loses its validity in some vicinity of the separatrix (at the pedestal foot) where the shear viscosity and the ergodization of the magnetic field play a significant role in the formation of the E 0r profile, and where the experimentally measured electric field tends to zero instead. This vicinity, however, is well separated from the resonant layers of interest in this study.

Anomalous diffusion.
In the previous version of the model [22], the anomalous diffusion coefficient, due to micro turbulence, was estimated either by a constant value D a 0 = 10 4 cm 2 s −1 or by a simple function D a = D a 0 (1 − 0.8(r/a) 3 ), where a is the minor radius of the plasma. In this work, we improve the previous estimations. We calculate the anomalous diffusion coefficient for the time slice t = 2.9 s with the transport code ASTRA [47], as this case is used in the analysis of scaling relations in section 4.2. For the other time slices, we estimate the anomalous diffusion coefficient based on heat fluxes in the plasma. In ASDEX Upgrade, the plasma is heated by neutral beam injection and electron cyclotron resonance heating. Subtracting from the heating power the radiation losses and assuming an equipartition of input power between electrons and ions, results in the total input power into the electrons of P e,inp = (P heat − P rad ) /2. For P rad we take the radiated power in the main plasma from the BPD diagnostic, which is the total radiated power of the plasma above the X point derived from the bolometric line integrals of radiation [48]. Respectively, the electron heat flux density is given by a steady-state value Q e = P e,inp /S, where the flux surface area S is calculated from the flux tube volume as S = dV/dr. Assuming the heat transport to be purely conductive, Q e ≈ −n e D a ∂T e /∂r, we arrive at an approximate expression for the anomalous diffusion coefficient A comparison of the anomalous diffusion coefficient profiles calculated with ASTRA and the flux estimation procedure is shown in figure 8. In this example, the results agree well in the region we are interested in (50-60 cm). As the flux surface area goes to zero when approaching the center, the flux estimation of the anomalous diffusion diverges. For time slices without results from ASTRA, equation (41) improves the estimation of the anomalous diffusion significantly, compared to the previously chosen constant value of D a = 10 4 cm −2 s −1 (dashed horizontal line in figure 8).
During the quasilinear time evolution, the anomalous diffusion would change due to the effects the RMPs have on the profiles, see e.g. [19]. In particular, since the transport barrier is eroded by the RMP-induced effects the turbulence becomes more important. This is for example prominent in the influence the RMPs have on the LH transition [49]. However, our model does not consider the evolution of turbulence as it is focused on the RMPs.

Identifying bifurcation in the quasilinear evolution
To identify the bifurcation of a mode for a given plasma configuration, i.e. profiles and equilibrium field, we linearly ramp up the RMP coil current during the quasilinear time evolution, while keeping the equilibrium magnetic field and shear fixed. We do it slowly, over a time scale of about two seconds which is much larger than the typical profile relaxation time at the edge. The quantity providing the information about the penetration of the RMP mode is the penetration factor B tot r /B vac r , which is evaluated at the resonant surface in question. Here, B tot r (B vac r ) is the radial magnetic field perturbation with (without) plasma response. A low value of this ratio indicates strong shielding. However, we are not interested in the value per se but in its behaviour during the ramp-up of the coil current. If the linear ramp-up of the current results in an abrupt non-linear incline in the penetration factor, we define this to be an indication for bifurcation.
In the following, we consider the ASDEX Upgrade discharge 33353 at time 2.9 s which is some time after the RMP coil current was turned on (∼1.7 s) and shortly after the transition to ELM suppression (∼2.77 s). The quasilinear evolution of the penetration factor in figure 9 shows a non-linear behaviour for mode m = (6, 2), which is indicating bifurcation of this mode. Moreover, the mode bifurcates approximately when the experimental MP coil current value is reached (I MP = 1.3 kA), meaning that for this plasma configuration, the mode might bifurcate. For mode m = (7, 2), there is no increase in the penetration factor visible at all. The mode m = (5, 2), on the other hand, seems to start bifurcating at the very edge of the ramp-up. However, the trend is not abrupt enough yet to be deemed bifurcating. If the RMP coil current value would be increased further, and, thus, significantly beyond the experimentally possible value, we might find this mode bifurcating.

Cause for bifurcation
Now, since we observe bifurcation of mode m = (6, 2) of the time slice 2.9 s we are interested in the specific processes involved in the bifurcation. To investigate the involved mechanisms, we scan the penetration factor at different stages of quasilinear evolution over the range of V E×B values evaluated at the resonant radius. We perform this scan by scaling the toroidal rotation velocity profile V tor and keeping the other profiles fixed. Figure 10 shows scans for different evolution stages during the ramp-up of the RMP coil current. We find that the quasilinearly evolving V E×B (green line) and the electron fluid resonance (orange line) move towards each other and eventually meet. There are two different mechanisms leading to the bifurcation.
First, the change of the quasilinearly evolved V E×B originates in the RMP-induced non-ambipolar particle transport that leads to a torque changing the toroidal rotation and, respectively, the E 0r profile globally. For an estimation of this trend, we apply a standard analysis [50] using the steady state momentum balance equation (23) alone. The only momentum source for the change in toroidal rotation frequency due to RMPs, This torque density is strongly localized around the resonant surface (see appendix B). Thus, using a zero boundary condition for the momentum flux at the axis the momentum equation can be integrated once over the radius, Here, T tot φm is the integral torque, Θ(x) is the Heaviside step function, and we ignored the finite width of the resonant layer. With a zero boundary condition for ⟨∆V φ i ⟩ at the separatrix we integrate (42) once more. Assuming the density and the temperatures unchanged, we express the resulting value of ⟨∆V φ i ⟩ at the resonant surface via the change in the V E×B velocity using equations (26) and (10) as Using the estimate (B4), we see that ∆V E×B = −|const|/V e res tends to reduce the distance to the fluid resonance V e res = 0, see equation (21). This is the main mechanism responsible for the bifurcation of the modes in the core plasma [13,14]. For the quantitative comparison below, it is convenient to estimate the integral torque T tot φ m assuming for simplicity no electron temperature gradient, where we used (29) and the explicit flux-force relation in (23). Using equations (21) and (13) for the resonant velocity V e res = V e⊥ and replacing all functions of radius in the sub-integrand of (43) with their values at the resonant surface, we can roughly estimate the relative change of the resonant velocity ∆V e res = ∆V E×B driven by the torque as Here, we used ⟨g φφ ⟩ ≈ R 2 0 and √ gB ϑ 0 /B 0 ≈ r m /q(r m ), and defined the distance to the separatrix, which is of the order of pedestal width, as L ped = a − r m . Further, the poloidal sound gyroradius is ρ pol = qR 0 v s /(r m ω ci ) with v s = (T e /m i ) 1/2 . Expressing the resonant layer width using (9) and (17) as  (44) as Estimating L ped ∼ 7 cm from figure 7 and T e ∼ 1 keV from figure 6 we get v s ∼ 3 · 10 7 cm s −1 and ρ pol ∼ 3 cm, which, together with ν e ∼ 2.5 · 10 6 s −1 , qR 0 ∼ 500 cm and k ⊥ ∼ 0.1 cm −1 results in C t/p ≈ 0.5 ∼ 1. The second mechanism, which moves the fluid resonance towards the gyrocenter resonance (black dash-dotted line in figure 10), is the result of the formation of a plateau in the electron temperature and density profiles by the RMP-induced transport within the resonant layer [22,51,52] (with a finite E 0r in the layer, this effect may lead even to positive density gradients [51]). Considering that a steadystate particle flux at the resonant layer is determined by the sources in the core plasma and, therefore, stays mostly unchanged with the appearance of the RMP-induced transport, Similar to (44), this change tends to reduce V e res via the reduction of A e 1 (see equations (21) and (13)). Due to C t/p ∼ 1 and µ A ⊥ ∼ D a , both effects are of the same order, as seen in the evolved profiles in figure 11. Note that the decrease in the shift of the resonance (purple arrow) also demonstrates the reduction of the electron temperature gradient. This has a similar effect on ∆V e res as the plateau formation and can be estimated by replacing in (47) D ql e,11 with D ql e, 22 being, in the end, of the same order.
It should be noted, that the formation of the plateau in quasilinear theory is not the result of magnetic islands opening. The width of the plateau structure is determined in this theory by the linear width of the resonant layer assuming that the size of the islands is small. At the same time, when the field fully penetrates the islands can become comparable in size or even larger than the linear resonance layer width. In the second case, the quasilinear theory loses its validity.

Plateau formation in experiment
As discussed above, quasilinear evolution with RMP-induced transport shows the formation of a local plateau in the electron temperature and density profiles. Due to the small radial extent of the plateau, it is not directly observable in experimental measurements. Thus, instead of identifying the plateau directly, we establish an estimation that may hint towards its existence. Figure 12. Assessment of the existence of the local electron temperature and density plateaus during the ELM suppression phase of discharge 33353 in ASDEX Upgrade. The existence of the plateaus is indicated by the rational surface m = (6, 2) lying between the modified fluid resonance with perfect electron temperature plateau (purple ×) and the fluid resonance without any plateau (blue triangle). For comparison, the location of a modified fluid resonance (red circle), with V evol e⊥ being the fluid velocity determined with evolved profiles in a multi-mode RMP run (m = 5, 6) evaluated at rm = r (6,2) , is indicated.
Note that the toroidal rotation velocity V tor measured in the experiment is that of boron impurities seeded to determine the rotation, though so far, we assumed it to be the one of deuterium. This assumption has an impact on the calculation of the poloidal rotation V pol and thus on the calculation of the radial electric field profile, E 0r . Therefore, to increase accuracy, for the plateau formation analysis we determined the radial electric field with the multi-species neoclassical transport code NEO-2 [36,37] assuming the toroidal rotation to be the one of boron. This is in contrast to the approach of calculating E 0r described in section 2.2 where we attributed toroidal rotation to the deuterium ions and computed the rotation in the cylindrical geometry limit (40).
We crudely estimate the existence of the plateau in the experiment in the following way. During the course of the ELM suppression phase of the discharge, where the plateau should already have formed, we determine various versions of the electron fluid resonance for different time slices in figure 12. Here, we assume that the alignment of the fluid resonance with the resonant surface m = (6, 2) causes bifurcation. Considering the original profiles, i.e. without any plateau, the electron fluid velocity crosses zero at some radial position marking the location of the (MHD) fluid resonance (blue triangle). If a perfect plateau is established in both electron temperature and density, i.e. the gradients are identical to zero, the diamagnetic part of the fluid velocity vanishes and the fluid velocity equals V E×B . Thus, the fluid and gyrocenter resonance coincide (green triangle). Additionally, the case of only a perfect electron temperature plateau but no density plateau is shown (purple ×). This case lies in between the gyrocenter resonance and the fluid resonance without plateaus.
To improve this crude estimation, we use our model to evolve the plasma parameters for each time slice while ramping up RMP modes m = (5, 2) and m = (6, 2) until a steady-state is reached. Note that we use here multiple modes since a bifurcation of mode m = (5, 2) might shift the fluid resonance towards the core. In the corresponding evolved profiles, plateaus have formed at the respective rational surfaces. To take this into account, we shift the original electron fluid velocity profile by the 'height' of the resulting local change in the evolved fluid velocity at the surface m = (6, 2), i.e. by (V evol e⊥ − V e⊥ ) r (6,2) . Here, V evol e⊥ is the fluid resonance calculated with the evolved density and electron temperature profiles. The result is shown in figure 12 by the red circles. The closeness of this 'most realistic' case to the relevant rational surface (m = (6, 2)) suggests the existence of the local plateau there.

Local bifurcation criterion
Based on the estimates of section 3.2 we define the approximate local bifurcation condition [21] of the RMP mode m via quasilinear, D ql e, 22 , and anomalous, D a , electron heat diffusion coefficients evaluated at the relevant resonant surface r = r m as D ql This heuristic condition is reasoned on the observation that an RMP mode bifurcates during the ramp-up when this condition is approximately fulfilled. The onset of plateau formation caused by D ql e, 22 reduces V e res which leads to the increased field penetration and respective increase in D ql e, 22 (see equations (20) and (31)) which further pulls the electron fluid resonance to the resonant surface, resulting in a positive feedback loop. If a very strongly shielded mode fulfills the criterion (48) this bifurcation can be rather abrupt (not the case in figure 11). Note that the unity on the right-hand side of the condition (48) is based Figure 13. Local bifurcation criterion (48) for shot 33353 at t = 2.9 s. For the rational surfaces m = (6, 2) the diffusion ratio is nearly above the threshold, in accordance with the exact criterion (see figure 9). To estimate the location of the pedestal top, the total pressure ptot is plotted. on experience and, thus, it has to be seen with error margins in mind.
In application to RMPs at the edge, the condition (48) appears more practical than the condition on the rotation [13,14] since it is local and is not affected by the interaction with other modes which may pull the rotation in opposite directions.
Considering the local criterion at time t = 2.9 s in figure 13, we find that the ratio of diffusion coefficients is nearly above the threshold for m = (6, 2) indicating that this mode might be close to bifurcation or might start to bifurcate. This agrees with the exact criterion, cf figure 9.

Local criterion during experiment
To check the correlation of the local bifurcation criterion with ELM suppression in the experiment, we determine the diffusion ratio for various time slices during shot 33353. We find that the local criterion for the RMP mode m = (6, 2) exceeds the threshold shortly before ELM suppression starts (t ∼ 2.7 s), as shown in figure 14. This indicates the bifurcation and the penetration of the respective mode. Additionally, this correlation supports the assumption that the bifurcation of the mode on top of the pedestal is responsible for ELM suppression.
Note that mode m = (7, 2) is already in the steep gradient region which is characterized by high plasma rotation, hence, this mode is strongly shielded. This fact is indicated by the low value of the diffusion ratio compared to the threshold value. Also, despite the diffusion ratio of mode m = (5, 2) being above the threshold quite early after the onset of the RMP coil current, there is no correlation to ELM suppression. This is as expected since the corresponding rational surface is already too far inside the plasma to be responsible for suppressing ELMs.
Moreover, the density evaluated at the resonant surfaces is plotted in figure 14. We see that ELM suppression starts shortly after the pedestal density (m = (6, 2)) is below the empirical threshold value [4] n e,ped = 3.3 · 10 19 m −3 . After some time in the ELM suppression phase, a single ELM event occurs (t ∼ 3.2 s), despite the density staying below the threshold value. This event is forecasted by the decline in the local bifurcation criterion of mode m = (6, 2) after ELM suppression starts. Nevertheless, the exact reason for this event is unknown.

Scaling relations
In the following, we want to determine scaling relations for the RMP coil current corresponding to the threshold D ql e, 22 = D a which follows from the local criterion (48) indicating bifurcation. To separate the main effect of the plasma parameters on the resonant layer from their effect on the 'outer' ideal MHD solution (the latter effect is of the order one or smaller), we assume that the shielding current in the resonant layer scales only with the RMP coil current but not with plasma parameters, I ∥m ∝ I RMP . Besides that, we fix also the anomalous coefficient D a . Thus, due to the quadratic scaling of the quasilinear diffusion coefficients with the shielding current (36), the scaling of the coil current threshold is determined solely by the scaling of the quasilinear electron heat diffusion coefficient

Numerical scaling.
We determine the scaling relations of D ql e, 22 given in equation (30) by scaling the original electron density and electron temperature profiles with constant factors, and keeping the ion temperature and toroidal rotation velocity profiles unchanged. Note that in this case the radial electric field profile, E 0r , is unchanged. Moreover, one should note that rescaling the density with a constant factor does not change the diamagnetic velocity which keeps the fluid resonance velocity (21) unchanged, in contrast to a temperature scaling.
We are interested in the density scaling, By scanning the quasilinear diffusion D ql e,22 in density space for multiple electron temperature values ( figure 15), we find a value range for the scaling for all three considered modes as Comparing the temperature averaged valuesᾱ m , we find that for mode m = (5, 2) the scaling is steepest withᾱ 5 = 1.92, whereas for modes m = (6, 2) and m = (7, 2) it is slightly smaller withᾱ 6,7 ≈ 1.4. Overall, this result agrees well with two-fluid nonlinear MHD simulations done by Nazikian et al [53] for DIII-D, which show a density scaling for the penetration threshold of the radial magnetic field for a mode resonant on top of the pedestal as B r ∝ n 0.7 e . Note that the scaling of the quasilinear heat diffusion coefficient with density in combination with the local criterion implies a density limit for bifurcation. This limit qualitatively resembles the density threshold for ELM suppression found in AUG [4] and DIII-D [54]. The AUG threshold value [4] is indicated in figure 15 by a vertical dashed line.  Density scan of the local bifurcation criterion for multiple electron temperatures for shot 33353 at t = 2.9 s. The empirical density threshold limit 3.3 · 10 19 m −3 determined in [4] is marked (vertical dashed line). The original pedestal temperature value is T e,ped = 1.14 keV.

Analytical scaling.
To get an analytical approximation of the trends of the local bifurcation criterion for given plasma parameters, we consider the expression for the quasilinear diffusion coefficient close to the resonant surface (31).
In constant-ψ approximation [13], we can use the kinetic integral plasma response current I e m∥ from equation (20) to substitute the radial magnetic field perturbation in (31) and we can write the local bifurcation criterion analytically as Here, the occurring equilibrium plasma parameters are evaluated at the resonant surface r = r m and V e res is defined in (21). A corresponding vacuum field required for bifurcation can be roughly estimated by replacing here I e m∥ with I cyl ∥m given by equation (35).
It should be noted that equation (51) is valid in the case of strong shielding, for which B r m given by equation (20) is much smaller than its vacuum value, |B r m | ≪ I ∥m /(cr m ), see section 2.1.4. If the saturation of |B r m | around its vacuum value is taken into account [13,50] in the case of a loss of shielding, equation (51) would not be singular at the resonances V e res = 0 and ω E = 0.
To assess the analytical version of our local bifurcation criterion (51), we compare it to the numerical one where D ql e,22 is calculated with equation (30) using the electromagnetic field provided by KiLCA. As shown in figure 16, we find that the analytical version of the bifurcation criterion gives the same qualitative picture as the numerical one, which is discussed in section 4.1. Thus, the analytical bifurcation criterion is advantageous when processing a large amount of cases.
Further, with equation (51), we check the density scaling of the local bifurcation condition for different collisionality limits. Recall that we assume D a to be constant. In case of high collisionality, ν e ≫ ω E , we get Whereas in the case of low collisionality, ν e ≪ ω E , we have Thus, the analytical approximation for both cases agrees with the numerically obtained scaling (50). Note that equation (51) as well as the numerical criterion (48) are valid for the case of intermediate or high collisionality, |x 2 | = |ω e |/ν e ≲ 1, which are typical for AUG and ITER parameters at the pedestal top. In case of low collisionality, |x 2 | ≫ 1, the maxima of the shielding current density and quasilinear diffusion coefficients (which are driven by Landau damping in this case) are split at the resonant surface. There, in turn, these quantities tend to zero due to the absence of resonant particles (see figure 3 and equation (A21)). In this limit, one should use in equation (48) with C ∞ ≈ 4.97. The same result, up to the coefficient of the order one, can be obtained formally by replacing ν e → max (ν e , |ω E |) in equation (51). Note that the analytical bifurcation criterion (51) and (54) as well as the torque estimate (B4) are based on the integral current (20) computed while ignoring the perturbation of the electric field E m⊥ . This field annihilates the corrugation field E [ψ] m⊥ in the ideal regions outside the resonant layers. In some regimes, this may reduce the effective resonant layer width compared to δ m of equation (45) which is suggested by the susceptibility functions, as seen, e.g. in figure 2. In fact, this effect is the main reason for the very low contribution of the ions to the shielding current since their δ m is larger than the  (51)) bifurcation criterion. Indicated are the original parameters (star), the threshold of the bifurcation criterion (thick dashed line), the empirically determined ELM suppression limits [4] (straight dashed lines), and the two-fluid MHD threshold determined with TM-1 by Nazikian et al [53] (red line). electron layer width in the typical case of collisionless ions. In the case that the electric perturbation field reduces the resonant layer width, equation (20) overestimates the shieling, and, therefore, the analytical criterion (51) should be viewed as a sufficient rather than a necessary condition for bifurcation, as it yields lower values of D ql e,22 compared to the numerical criterion (see, e.g. figure 16 for mode m = 7 at t < 2.7 s).

Density and electron temperature parameter space.
Similar to the operational boundaries in pedestal electron density and electron temperature space of figure  15 in [4], we want to establish regions of bifurcation in the two-dimensional density and electron temperature parameter space. To do so, we scale the respective profiles as described in section 4.2.1 and check D ql e, 22 .
Mode m = (7, 2) lies in the steep gradient region of the considered profiles, i.e. the high rotation causes strong shielding. In figure 17(c), we find that the threshold of bifurcation for this mode lies at very low pedestal densities. Though since the bifurcation of a mode in the steep gradient region may deteriorate the edge transport barrier completely [53] it is not desirable for this mode to bifurcate anyway. Mode m = (5, 2) shown in figure 17(a), which lies already too far inside the plasma to be responsible for ELM suppression, seems to be close to an extended region in which bifurcation occurs that can be identified as the fluid resonance.
Finally shown in figure 17(b) is mode m = (6, 2), were we find that the threshold of the bifurcation criterion qualitatively agrees with the structure of figure 15 in [4], as well as a comparable scaling established in [53] for DIII-D. Our results show a 'knee' structure close to which the experimental parameters (star) lie. This structure is absent in the fluid analysis (red line) and is caused by the fluid resonance V e res = 0. Hence, the fluid resonance provides a more or less extended operational region for which bifurcation occurs. Further, comparing the fluid resonance structure with the empirically observed thresholds [4] n emp e,ped = 3.3 · 10 19 m −3 and T emp e,ped = 1.0 keV, we find qualitative agreement, as the knee resembles boundaries in both density and electron temperature. Additionally, a limit similar to the lower collisionality limit (ν * e = 0.25) shown in figure 15 (dotted line) in [4] occurs in figure 17(b) at the bottom side of the knee. Still, due to the qualitative nature of our approach, quantitatively the limits do not coincide.
Further, scans done with the analytical criterion (51) are shown in figure 17(d)-(f ). Quantitatively, the constant-ψ approximation introduces differences. In particular, the fluid resonance 'knee' of mode m = (6, 2) is squeezed to lower temperature values. However, qualitatively, the analytical result provides a similar picture as the numerical one.

Summary and discussion
In this work, we advanced an earlier kinetic model [11,22,23,25] by coupling the cylindrical Maxwell solver KiLCA [23] to the ideal MHD code GPEC [24] to incorporate toroidal effects. We used this extended model to qualitatively study the bifurcation of RMPs, which are used to suppress ELMs. Based on observations in the quasilinear time evolution, we defined an approximate local bifurcation criterion (48) that suggests bifurcation if the ratio of the quasilinear electron heat diffusivity is of the order of the anomalous diffusivity. A more approximate, analytical form of the local criterion (51), which can be viewed as a sufficient condition for the bifurcation of a given RMP mode, was derived in constantψ approximation, i.e. assuming a constant radial magnetic field perturbation over the resonant layer, and ignoring the electric field perturbation. The criteria were used to analyse multiple time slices of an ASDEX Upgrade discharge. Both criteria suggest bifurcation of the RMP mode m = (6, 2) just before the transition to ELM suppression (figure 16). Further, using the criteria to investigate electron density and electron temperature scalings of bifurcation we find a limit in the electron density that resembles an empirically observed limit [4]. The limit we find varies with the electron temperature and is strongly influenced by the electron fluid resonance.
At this point, it is important to note two caveats. First, the mode we find to be correlated with ELM suppression, m = (6, 2), differs from the one suggested by experiments. In the latter case, the expectation is that mode m = (7, 2) or m = (8, 2) are connected to ELM suppression [4]. So far, we do not know the reason for this mismatch. Second, experimental measurements are inevitably accompanied by uncertainties, however, so far our modelling lacks an uncertainty analysis. This will be remedied in future works.
In the following, we summarize and discuss the five most important points that we observe within our model.
(i) Both RMP shielding and RMP-driven quasilinear transport are determined solely by the perpendicular electrostatic field arising due to the misalignment of perturbed equipotential surfaces and perturbed magnetic flux surfaces. This field is absent in ideal MHD. (ii) The shielding of an RMP mode breaks down if the electron fluid resonance is close to the resonant surface. However, we find in our model that the electron fluid resonance condition, V e res ≡ V e⊥ − 1 2 V e dT = 0, which must be realized at the resonant surface r = r m , differs from the MHD predicted value of this resonance at zero electron fluid velocity, V e⊥ = 0. Therefore, the electron fluid resonance point is normally located closer to the gyrocenter resonance E 0r = 0 than the zero of the electron fluid velocity. This can be one of possible reasons of experimentally observed ELM suppression in case that the point V e⊥ = 0 moves further inside from the pedestal [54]. This change in the fluid resonance position in the presence of an electron temperature gradient is the case for arbitrary plasma collisionality. For the collisionless regime, this has been pointed out earlier [22]. Besides the above shift, nonconstant-ψ effects, i.e. a radially increasing radial magnetic field perturbation, tend to shift the fluid resonance position V e res (r pen ) = 0, required for maximum field penetration, radially outwards from the resonant surface, r pen > r m . This typically results in an opposite trend as compared to the effect of a finite electron temperature gradient. (iii) In contrast to the result of the simplified analytical constant-ψ approximation, a full penetration at the fluid resonance is not observed in the presence of an electron temperature gradient at low and mild collisionalities ( figure 4). Although the reduction of the maximum penetration factor is of the order one, it effectively increases the bifurcation threshold for the RMP coil field. (iv) The analytical result of our kinetic plasma response model can be compared to the drift-MHD response regimes of Cole and Fitzpatrick [14] via the measure for the strength of the current sheet ∆. In our model, the measure is determined by equating the currents in equations (20) and (33), where ω 2 pe = 4π n e e 2 /m e is the electron plasma frequency. At high collisionality, ν e ≫ |ω E |, and in the absence of an electron temperature gradient, V e res = V e⊥ , this formula agrees with the second semi-collisional regime (SCii, table I) of Cole and Fitzpatrick [14] to within 4%, which is due to different numerical prefactors. Note that at low collisionalities, where the shielding is fully determined by Landau damping, ∆ is independent of the mode numbers. In this case, no corresponding regime is found in drift-MHD theory. Note also that the ratio of the perturbation field in plasma to its value in vacuum can be approximated by equation (55) in a view of equation (34) as B r m /B r Vm ≈ 2m/(r m ∆). (v) The bifurcation criterion (51) strongly depends on the electron density and temperature but not on the collision frequency. This means that the effective ion charge number Z eff is predicted to have a weak or even no effect on the penetration. Hence, bifurcation is predicted to occur under almost the same conditions for plasmas with different ion species.
Note that RMP shielding and quasilinear electron transport induced by the misalignment electric field in the resonant layer is only weakly affected by toroidal effects because the whole electron population (essentially passing electrons) is involved in these processes. An account of the misalignmentinduced transport of trapped particles in general toroidal geometry [55] gives a correction in the resonance layer of the order (|ω E |/ν e )(a/R 0 ) 3/2 ≪ 1. This transport, however, may be significant outside resonant layers. In the case that RMPs couple with an electrostatic (drift) mode propagating further in the pedestal region, this could play an important role.
Suppressing ELMs requires bringing the plasma-edge pressure below a stability limit. In that regard, the experimentally observed density 'pump-out' plays a key role and needs to be explained. Concluding this article, it appears that the bifurcation of RMPs is a necessary but not a sufficient condition for ELM suppression. There could be more than one mechanism involved. Locally, there is the formation of vacuumsized island chains after the bifurcation which prevents pedestal propagation to the core [16,56]. A further mechanism may be due to the strongly increasing RMP-driven torque [57][58][59]. Such a torque may modify the edge rotation globally, in particular in the pedestal, leading to the change in the turbulence spectrum and, consequently, in the anomalous transport [60]. Moreover, as already mentioned, coupling of RMPs to electrostatic drift modes might produce relevant transport across the pedestal. To fully resolve drift modes, an upgrade of the FLR expanded conductivity operator to a full integral one is the necessary next step in the improvement of our model. With the integral model available, analysing drift modes will shed light on their role in density 'pump-out' and, subsequently, on ELM suppression.
where only the first term in (A19) responsible for Landau damping within drift-orbit resonance between parallel motion and E × B rotation contributes while the contribution of adiabatic particles integrates to identical zero. An expression which reproduces both, collisional (A14) and collisionless (A23) limits is the following, So far, the authors were unable to check analytically the observation that equation (A24) perfectly fits the result of numerical evaluation for any x 2 (see figure 18). In close vicinity of the resonant surface where z α ≫ 1 using the respective asymptotic of Kramp function (A19) functions (A22) approximately are This expression would be the same with (A15) used formally in the limit x 2 ≫ 1. Actually, it can be checked using equation (A10), which is valid for small x 1 but not necessarily small x 2 , that equation (A15) is valid in the limit x 1 → 0 in the general collisionality case.
here since it does not produce torque on average [50]. (Since the Lorentz force density equals the pressure gradient in ideal MHD, it cannot produce any integral torque.) Assuming a single perturbation mode m and integrating (B2) over the volume in constant-ψ approximation, i.e. replacing all functions of radius except for j α m∥ with their values at the resonant surface r = r m , we get the integral torque from a single mode as T tot φm = − π σ ϑ r m R 2 0 2c q 2 (r m )R 2 0 + r 2 m Re I ∥m B r * m (r m ) .
Here, I ∥m is the integral parallel current (15), σ ϑ = sign √ gB ϑ , and we have set √ g = rR 0 and B ϑ 0 /B 0 = q 2 R 2 0 + r 2 −1/2 for the straight cylinder. Substituting here a vacuum estimate for the shielding current, equation (35) agrees with equation (67) of [50]. Assuming strong shielding by electrons, i.e. using equation (20) for B r * m with I ∥m ≈ I e ∥m , we can express the integral torque (B3) via the shielding current as where V e res is the electron fluid resonance velocity (21). By similar arguments as in the discussion of equation (51), in the case that shielding is lost in the vicinity of a fluid resonance, the torque has no singularity there. In fact, it smoothly goes through zero at V e res = 0 which is a stable point of the toroidal momentum balance, e.g. see figure 9 in [25]).
It should be noted that the above estimates of the (dominant) electron torque ignoring the perturbation electric field E m⊥ should not be used in close vicinity of the gyrocenter resonance ω E = 0, where the effect of E m⊥ is significant. Especially, for ions, whose torque exceeds the electron torque there, the ω E = 0 resonance is a stable point as well (see figure  11 of [22])