MHD simulations of small ELMs at low triangularity in ASDEX Upgrade

The development of small and no-ELM regimes for ITER is a high priority topic due to the risks associated with type-I ELMs. By considering non-linear extended magnetohydrodynamic (MHD) simulations of the ASDEX Upgrade tokamak with the JOREK code, we probe a regime that avoids type-I ELMs completely, provided that the separatrix density is high enough. The dynamics of the pedestal in this regime are observed to be qualitatively similar to the so-called quasi-continuous exhaust regime in several ways. Repetitive type-I ELMs are substituted by roughly constant levels of outward transport, caused by peeling-ballooning modes (with dominant ballooning characteristics) which are localised in the last 5% of the confined region (in normalised poloidal flux). The simulated low triangularity plasma transitions to a type-I ELMy H-mode if the separatrix density is sufficiently reduced or if the input heating power is sufficiently increased. The stabilising factors that play a role in the suppression of the small ELMs are also investigated by analysing the simulations, and the importance of including diamagnetic effects in the simulations is highlighted. By considering a scan in the pedestal resistivity and by comparing the poloidal velocity of the modes to theoretical estimates for ideal and resistive modes, we identify the underlying instabilities as resistive peeling-ballooning modes. Decreasing the resistivity below experimentally-relevant conditions (i.e. going towards ideal MHD), the peeling-ballooning modes that constrain the pedestal below the type-I ELM stability boundary display sharply decreasing growth rates.


Introduction
The thermonuclear experimental fusion reactor ITER is foreseen to operate in high-confinement mode (Hmode), which is characterised by the quasi-periodic excitation of type-I ELMs (edge localised modes) [1].In H-mode there are reduced levels of turbulent transport in the edge of the confined region, thus forming a narrow transport barrier that creates a 'pedestal' in the density and temperature profiles.Type-I ELMs are macroscopic magnetohydrodynamic (MHD) instabilities that are destabilised by high pressure gradient and/or high toroidal current density (and current density gradient).In standard ELMy H-mode, the pedestal rises (thus increasing ∇p which simultaneously increases j tor by the formation of the bootstrap current) until type-I ELMs are excited.The appearance of type-I ELMs causes the pedestal to crash in a timescale of 10 2 ∼ 10 3 µs, and it is followed by a quiet inter-ELM phase that lasts 10 1 ∼ 10 3 ms.The large transient heat loads associated to these ELMs must be avoided in future tokamaks in order to achieve an acceptable divertor lifetime [2,3].
Naturally ELM-free (e.g., QH-mode, EDA Hmode) and ELM-mitigated (e.g., RMP mitigation, pellet pacing, grassy ELMs, small ELMs) operational conditions have been successfully achieved using several methods in different existing tokamaks.However, naturally ELM-free regimes and ELM-mitigated scenarios are only achievable in reduced parameter regimes that differ between different tokamaks and, therefore, it is uncertain which methods and regimes will be accessible in future tokamaks [4].Extrapolating such regimes to ITER becomes increasingly uncertain because ITER parameters cannot be simultaneously accessed with existing tokamaks (e.g., density or collisionality, not both at the same time) [5].
Certain small ELM regimes are completely free of type-I ELMs, maintain other desirable features, and could, if accessible, signify an attractive option for ITER.One such regime is the quasi-continuous exhaust (QCE) regime, which is routinely operated in ASDEX Upgrade (AUG) and TCV and completely avoids type-I ELMs while maintaining good confinement properties [6,7].In order to access the QCE regime, it is necessary to operate at high separatrix density (n e,sep 0.3 × n GW ), with high triangularity and close to double null [7].The physical mechanism that constrains the pedestal beneath the type-I ELM stability boundary remains unclear, but it is thought to be high-n ballooning modes that are located near the separatrix which cause sufficient outward transport.These modes give rise to a quasicontinuous exhaust of heat and particles which impinge onto the divertor.It has been found that the power fall-off length is larger during the QCE regime than expected from the empirical Eich scaling [8].Other small ELM regimes that can be free of type-I ELMs include grassy ELMs (observed in JT-60U [9], JET [10], EAST and it is foreseen as a potential operational regime for CFETR [11]), type-III ELMs (not reactor-relevant because they cause confinement degradation [12]), and a low density small ELM regime in JET [13].
Naturally ELM-free operation includes regimes like QH-mode (in DIII-D [14], AUG and JET with carbon wall [15,16], and JT-60U [17]), I-mode (features a pedestal in the temperature but not in the density profile [4]), and EDA H-mode (found in Alcator C-mod [18]).The latter, i.e., the enhanced D-alpha H-mode, was recently achieved in AUG [19].First experiments of EDA H-mode in AUG were performed at low triangularity with pure electron heating and observed a narrow operational window in terms of the applied heating power.Recently it has been found that the operational window can be extended to higher heating powers by increasing the plasma triangularity; heating above said operational window results in a type-I ELMy H-mode [20].This operational regime always features an edge quasi-coherent mode (in a frequency range between 20 and 80 kHz).
The JOREK non-linear extended MHD code [21,22] has been extensively used to simulate macroscopic edge instabilities in tokamaks plasmas.In particular for AUG, it has been used to produce realistic simulations of type-I ELMs [23,24], RMP-ELM mitigation and suppression [25], and pellet-triggered ELMs [26,27].The present article details JOREK simulations of small ELMs at low triangularity in AUG and discusses their relation to the small ELMs that underlie the QCE regime and to the QCM of EDA H-mode.The simulations presented here follow from the approach for modelling the pedestal build-up described in Ref. [23].In simulations at sufficiently high separatrix density (n e,sep ≈ 0.4 × n GW ), small ELMs appear beneath the type-I ELM stability boundary and feature medium-n resistive peelingballooning modes near the separatrix which cause quasi-continuous heat exhaust.The present article is structured as follows.A brief description of different types of small ELMs in AUG and some features of the EDA H-mode are presented in section 2. The JOREK model used for the present simulations together with the simulation set-up details and the axisymmetric pedestal build-up are presented in section 3.In section 4 the results of the non-axisymmetric simulations are presented and a detailed analysis is provided.Two different paths to leave the small ELM regime and reach a type-I ELMy H-mode are presented in section 5. Finally, conclusions and outlook for future work are discussed in section 6.

Small ELMs and EDA Hmode at ASDEX Upgrade
Small/no ELM scenarios feature different transport mechanisms that cause losses below a few percent of the plasma stored energy.Type-III ELMs are observed in AUG as distinct peaks in the D α signal with a roughly constant frequency.Each event causes an expulsion of 5% of the plasma stored energy and their repetition frequency decreases with increasing heating power (f ELM ∝ 1/P heat ); their repetition frequency is typically larger than type-I ELMs, f type−III ∼ 10 3 Hz.They can be obtained either at low pedestal density close to the L-H power threshold or at higher heating power by increasing the pedestal density [12,28].Type-III ELMs are thought to be resistive instabilities, and they are associated with poor confinement properties.JET-like simulations of repetitive edge instabilities that featured an inverse dependency between repetition frequency and heating power have been achieved with JOREK [29].The term "small ELMs" has been used in AUG as a broad category that includes small amplitude ELMs but excludes type-III ELMs.In particular, this considers type-II ELMs and the ELMs that underlie the QCE regime.The QCE regime is posited to be an attractive scenario for ITER because it completely avoids type-I ELMs while maintaining good confinement properties.Another favourable feature of this regime is that it deposits the expelled energy in a quasi-continuous manner and in a broader area than observed in the inter type-I ELM phase [8].Small ELMs act as a transport mechanism that expels heat and particles such that the pedestal cannot build-up to a point where type-I ELMs are excited; it is presently hypothesised that small ELMs are ballooning modes and/or high-n peeling-ballooning modes that are located at, or very near, the magnetic separatrix.The important ingredients for maintaining the QCE regime are high separatrix density (n e,sep /n GW 0.3) and closeness to double null (together with high triangularity) [6,7,30].Because there is little variation of separatrix temperature in a given device (T e,sep ≈ 100 eV for H-modes in AUG [31]), high n e,sep translates to high separatrix collisionality (ν * e ∝ n e /T 2 e ).In existing tokamaks, high ν * e,sep implies also high pedestal collisionality because the temperature cannot increase arbitrarily due to the excitation of ELMs.On the other hand, ITER can reach higher pedestal top temperatures and is expected to operate with high ν * e,sep and low ν * e,ped .Such conditions cannot be achieved simultaneously in existing tokamaks; therefore, fundamental uncertainties exist on whether ITER could operate in the QCE regime [6].
Depending on the triangularity and the edge safety factor, increasing heating power can either lead to a sustained small ELM regime (with "standard" QCE parameters: high triangularity, high q 95 , and close to double null) or to a transition from small ELMs to type-I ELMs (with "standard" QCE parameters, but lower q 95 ), as shown in fig. 1.For instance, discharge #35572 (−2.0 T, 1.0 MA, δ = 0.397, low safety factor q 95 = 3.73, P heat = 10 MW) heated with ICRH and NBI (P rad ≈ 4.1 MW) shows small ELMs, and transitions to a type-I ELMy H-mode upon increasing the heating power to P heat = 14.6 MW (P rad ≈ 5.2 MW).
Even with a single null configuration and low triangularity, small ELMs can be achieved at sufficiently high separatrix density.Nevertheless, such small ELMs are associated to degrading confinement and even to an H-mode density limit (n e n GW ).In such cases as the separatrix density is increased (by increasing the gas puff rate), filamentary transport also increases [32,33] and can lead to a flattening of the pressure gradient.This, in turn, causes a reduction of the edge radial electric field 1 which is what causes the back transition to L-mode, i.e., the H-mode density limit (HDL) [32].Recently, a correlation has been observed between ballooning stability at the separatrix and the onset of the HDL in AUG [35] and in JET-ILW [36].In this context, the deterioration of the H-mode confinement and, ultimately, the breakdown of the H-mode are explained by an excess of cross-field transport caused by small ELMs located Figure 1: From top to bottom, the heating power (radiated power in dotted lines), plasma stored energy, and ELM monitor of two AUG discharges with small ELMs that feature increasing heating power.A discharge with low edge safety factor (left; #35572, q 95 = 3.73) changes from a small ELM-dominant regime to a type-I ELM-dominant regime upon increasing the heating power.The discharge with high q 95 (right; #39565, q 95 = 5.76), on the other hand, remains in a small ELM-dominant regime throughout the heating power steps.near the separatrix.
The enhanced D α H-mode is a no-ELM operational scenario with high density which is potentially attractive for ITER.Type-I ELMs are not destabilised during EDA H-mode operation because the pedestal is not able to build-up sufficiently.The transport mechanism that allows the pedestal to remain below the type-I ELM stability boundary is thought to be an electromagnetic mode dubbed quasi-coherent mode (QCM) [37].At low triangularity, the EDA H-mode lives in a very narrow operational space in terms of the applied heating power; however, in AUG this operational space has been observed to expand when increasing the plasma triangularity [20].Nonetheless, a stationary EDA H-mode inevitably transitions to a type-I ELMy H-mode upon a sufficient increase of the heating power.In AUG, the QCM moves in the electron diamagnetic direction and has fluctuation frequency in the range f ≈ 20 − 80 kHz [19].Similarly, type-II ELMs are accompanied by broadband fluctuations (in magnetic pick-up coils and ECE signals) with frequencies in the range f ≈ 30 − 50 kHz [38].

JOREK simulation set-up and axisymmetric build-up
The simulations presented in this paper were produced with the reduced MHD version of the 3D non-linear extended MHD code JOREK [21,22].Said model simplifies the visco-resistive MHD equations with two considerations.First, the toroidal magnetic field is constrained to be time-independent (B ϕ = F0 R φ, where F 0 is a constant, R is the major radius, and φ is the toroidal coordinate).Second, the poloidal velocity is considered to be comprised of the ExB velocity, which is further considered to lie in the poloidal plane, i.e., v pol = v ExB ; this assumption allows a potential formulation for the poloidal velocity through the electrostatic potential E = −∇Φ.It is then possible to include diamagnetic effects (as an extension to the reduced MHD model) by considering the poloidal velocity to be where m i and p i are the ion mass and pressure respectively, e is the fundamental electric charge and ρ is the mass density.The resulting extended MHD model is a closed system of five equations for the poloidal magnetic flux (ψ), the electrostatic potential (Φ), the parallel velocity (v ), the mass density (ρ), and the single fluid temperature (T ).The inclusion of diamagnetic effects allows the JOREK simulations to recover realistic radial electric fields in the pedestal region, i.e., E r ∝ n −1 i ∇p i [34], which play a fundamental role in the stability of PB modes (particularly those with high toroidal mode numbers) [39].The stability of PB modes is also partly determined by the edge current density (and its gradient), which is comprised of an Ohmic contribution and a bootstrap current contribution.The former is included in JOREK by assigning a current density source determined by the initial current density profile (which is also comprised by Ohmic and bootstrap current contributions, j 0 = j 0,Ω + j 0,bs ).The time-evolving contribution from the bootstrap current density is considered in JOREK by making use of the Sauter analytical expression [40,41].The total current density source then corresponds to j source (t) = j 0 + [j bs (t) − j 0,bs ] = j 0,Ω + j bs (t), and it is included in the induction equation: 1 is the Poisson bracket, and η is the resistivity.The extension to include the bootstrap current density as a source term in JOREK was introduced in Ref. [42], and it has been used in recent simulations of ELMs [23,26,27].All simulations presented in this paper consider diamagnetic effects and the bootstrap current density source unless specified otherwise.The remainder of the section describes the numerical setup of the simulations and the axisymmetric build-up of the pedestal profiles.

Numerical parameters and simulation set-up
Initial conditions which are stable to ideal peelingballooning modes are considered for this work.These are obtained from a post-ELM equilibrium reconstruction of AUG discharge #33616 at roughly 7 s with the CLISTE code [43].The magnetic field at the magnetic axis is 2.5 T and the plasma pressure is I p = 0.8 MA.It considers a lower single null magnetic configuration with low triangularity δ av = 0.29 and with the ion B × ∇B drift direction pointing towards the active X-point.The reconstructed current density profile is constrained from based on the steepness of the density, temperature, and pressure profiles, i.e., a bootstrap current constraint [44].
The electron density, plasma temperature and pressure, and toroidal current density initial profiles at the outboard midplane are shown in fig. 2. The plasma temperature, T , is the sum of the ion and electron temperatures, which are assumed to be equal.The toroidal current density is comprised of an Ohmic contribution together with a Pfirsch-Schlüter and a bootstrap current contribution2 .The bootstrap current is only a small contribution in the initial profiles due to the small steepness of the initial pressure profile.An ideal MHD stability analysis with the MISHKA code indicates that these post-ELM pedestal profiles are stable (to ideal PB modes) and would first need to steepen in order to excite a type-I ELM.
We use a flux-surface aligned grid that considers the confined region, the scrape-off layer, and the private flux region.The grid is comprised of 138 points in the radial direction (120 points in the confined region and 18 in the scrape-off layer) and 354 in the poloidal direction.A convergence scan in the grid resolution has been performed; we observe that the linear growth rates of instabilities with n ≤ 20 do not change by changing the radial and poloidal resolution (higher mode numbers were not probed because the dominant modes are n ≤ 12 due to diamagnetic stabilisation).For the axisymmetric build-up, only one poloidal plane is considered, and for the non-axisymmetric simula- For a plasma temperature of 1 keV and density of 5 × 10 19 m −3 , i.e., ψ N ≈ 0.8 in the post-ELM equilibrium shown in fig.2, χ ,SH e ≈ 1.27 × 10 9 m 2 /s.For the simulations presented in this paper, the parallel heat diffusion is roughly 15 times lower than the Spitzer-Härm value; as motivated by the heat-flux limit which can account for a reduction of χ SH by a factor of 15 − 130 [46].Fast parallel heat transport acts as a stabilising agent for PB modes [47].

Axisymmetric pedestal build-up
Together with the post-ELM equilibrium reconstruction of AUG discharge #33616, a pre-ELM reconstruction from the same discharge was considered.An ideal MHD stability analysis (with the code MISHKA) shows that the pre-ELM profiles are unstable to ideal PB modes.Starting from the post-ELM profiles, we impose ad-hoc diffusion coefficients (that describe a well in the pedestal region) and sources that drive the pedestal towards the pre-ELM profiles.The diffusion coefficients and sources are constant in time as the simulation progresses.This evidently results in a pedestal build-up at fixed pedestal width, which is a simplification of what is experimentally observed.Time-evolving diffusion coefficients and sources would require including several key physical effects that are beyond the scope of MHD, and which will be investigated in future work (more on this subject in subsection 3.3).
As the pedestal n e and T evolve due to the stationary diffusion and sources, the radial electric field and j become driven by the increasing influence of diamagnetic effects and the bootstrap current density, respectively.Figures 3(a)-(d) show the time evolution of the pedestal in terms of the electron density, plasma temperature, radial electric field, and the flux-surface averaged toroidal current density.The colours of the profiles change gradually from purple to blue with increasing time as shown in the colour bar on top of the figure.The profiles are plotted every 0.2 ms during the first 4 ms of an axisymmetric simulation.In the first millisecond, the profiles change shape quickly (as can be most clearly evidenced in the evolution of the density and radial electric field).Section 4 will show that non-axisymmetric instabilities driven by the steepening profiles prevent the pedestal from building up towards the profiles shown in fig. 3. From linear stability analysis performed with the MISHKA code, we know that the pedestal build-up shown in fig. 3 does not cross the ideal peeling-ballooning boundary in the time shown.

Limitations of the present approach
The pedestal build-up considered for the simulations presented in this work assumes a constant level of turbulent and neoclassical transport.Experimentally, however, turbulent and neoclassical transport in the pedestal is known to evolve in sub-millisecond and millisecond timescales, respectively.These changes in transport levels would determine how exactly the pedestal top and width evolve.In that sense, to produce a fully realistic pedestal build-up, which would include the pedestal widening, it is necessary to run neoclassical and gyro-kinetic (or even kinetic) simulations to determine the dynamical turbulent and neoclassical transport throughout the simulation time.This would represent not only an extremely costly endeavour from the point of view of computation time, but also would require significant efforts in terms of code development, which lie well beyond the scope of the present work.In the future, neural network based reduced models for the turbulent transport coefficients could potentially be used for incorporating such effects in the MHD simulations (said reduced models, however, do not include transport in the pedestal).
The ion and electron species can be approximated to have the same temperatures and densities in the pedestal region, but in reality ions and electrons behave in distinct ways due to the large difference in their respective masses.Such effects are neglected in the present simulations since we use the single fluid version of the JOREK code.However, a two temperature model has been developed in JOREK and it will be used in the future to understand the effect of such temperature separation in ELM physics and advance the level of realism in our simulations.Another important physical effect that is not considered in the present approach is the penetration of neutral particles onto the confined region and their interaction with the plasma (more generally speaking, a more complete SOL/divertor model is missing).The ensuing ionisation of the neutral particles would determine the amount of particle fuelling that should be considered at any given time during a simulation.These fuelling effects, in turn, directly influence the density (and ultimately temperature) profiles in the pedestal.Ongoing efforts are underway that permit JOREK simulations to consider such effects either by a kinetic treatment [48] or a fluid treatment [49,50] of the neutrals.

Non-axisymmetric simulations
The present section details the simulation results by first describing the linear growth phase of the instabilities together with the early non-linear phase (phase during which n = 0 modes interact with each other but not with the n = 0 axisymmetric background), which take place roughly during the first millisecond of simulation time.The poloidal velocities of the lin-early unstable modes are measured during this linear growth phase, and the impact of varying the resistivity onto the linear growth rates is probed.These analyses allow us to identify the underlying instabilities as resistive peeling-ballooning modes located near the separatrix.
During the axisymmetric build-up shown in fig.3, the pedestal does not cross the ideal PB boundary (as confirmed by ideal MHD stability analysis with MISHKA).However, non-ideal instabilities can become excited due to the finite resistivity used in JOREK.A simulation which includes the toroidal mode numbers n = 0, 2, 4, . . ., 12 is performed by introducing perturbations with said mode numbers at noise-level amplitudes after t ≈ 0.1 ms of axisymmetric build-up.After a brief time of stability (∼ 0.2 ms) the steepening pedestal that is not unstable to ideal PB modes begins to excite non-ideal PB instabilities with predominant ballooning features.The magnetic energies of the non-axisymmetric perturbations can be observed in fig.4(a).The time frame between the first vertical black line (t = 0.3 ms) and the vertical purple line (t = 0.5 ms) denotes the linear growth phase where only the linearly unstable modes (n = 6, 8, 10, and 12) grow; the time frame between the purple line and the second vertical black line (t = 0.7 ms) denotes the early non-linear growth phase where non-linear mode coupling excites linearly stable modes to grow (n = 2 and 4) 3 .The structure of the PB modes of fig.4(a) during the linear phase (at t = 0.4 ms) is shown in fig. 5 with the perturbations of density, temperature, and poloidal magnetic flux.The flux surfaces at ψ N = 0.95, 1.00, 1.05, 1.10 are also shown in thin black lines.It is observed that the PB modes are localised very close to the separatrix.The modes rotate in the electron diamagnetic direction (counter-clockwise in fig.5); their poloidal velocity will be discussed later.
The initial growth phase is started by an n = 12 mode in these simulation and closely followed by the growth of the n = 10 mode.The growth rates of these modes are very similar, roughly γ n=10,12 ≈ 5 × 104 /s.A separate simulation including all even toroidal mode numbers until n = 20 has been produced, and the magnetic energies of the n = 0 modes is shown in fig.4(b).In such a way, we confirm that indeed the fastest growing mode is the n = 12.The toroidal harmonics with n > 14 grow at a slower rate; they remain sub-dominant well into the non-linear phase (not shown).Quadratic non-linear mode coupling takes place during the early non-linear phase (t = [0.5, 0.7] ms).The non-linear mode coupling that gives rise to the excitation of linearly stable modes has been described in JOREK simulations  from Ref. [51].The upcoming subsection is devoted to studying the influence of the plasma resistivity onto the growth rates of the non-axisymmetric perturbations.

Rigidly scanning the resistivity
In this subsection, we freeze the axisymmetric profiles at t = 0.5 ms and change the resistivity to understand its influence onto the stability of the nonideal PB modes.For this scan we consider several multiplication factors of the nominal resistivity, 0.5, 1.5, 3.0, 6.0, 12.0, 24.0, 48.0.The response of the PB modes to the changes in resistivity display several noteworthy characteristics.The linearly stable modes (n = 2 and 4) at nominal resistivity become linearly unstable by increasing the resistivity (i.e., said modes no longer require non-linear mode coupling to grow).Similarly, decreasing the resistivity by half leads the n = 6 mode to become linearly stable.This is not surprising given the fact that we know (from MISHKA ideal MHD simulations) that ideal PB modes are not unstable for the considered profiles.Varying the resistivity results in a change in the growth rate of the linearly unstable modes, as shown in fig.6.The x-axis represents the resistivity (in Ω m) and the y-axis corresponds to the growth rate (in 1/s); both axes are plotted in logarithmic scales.Different colours and symbols represent different toroidal mode numbers.

Linear growth phase -mode velocity
Considering how the peak of the modes in fig. 5 where n is the toroidal mode number.Figure 7 shows the ψ N = 0.92 and 0.99 flux surfaces together with the colour coded arc length, s, calculated from the inboard midplane (a) and the ψ n perturbation for the different toroidal mode numbers (b)-(g) at t = 0.5 ms.The n ≥ 8 mode amplitudes dominate in the LFS indicating their ballooning nature while the modes with lower n do not have a coherent structure (the n = 6 already has a coherent structure, but its amplitude is too small to be observed in fig.7(e)).
Taking the distance travelled by the peaks of fig.7(b)-(e) in a small time, it is possible to deter- mine the poloidal velocity of the different modes.The poloidal mode velocity of the linearly unstable modes (n = 6, 8, 10, 12) at the outer midplane is approximately constant from ψ N ≈ 0.92 − 0.99, and corresponds to v mode,pol ≈ 11 km/s.The poloidal velocity of a mode may be used to attempt to identify the nature of the mode.In particular, from Ref. [52], ideal and resistive ballooning mode rotation velocities in the laboratory frame have been identified as having the following velocities, Resistive : Ideal : where v * i = ∇p i /(en e B) is the ion diamagnetic velocity, and v ,pol is the poloidal projection of the parallel velocity.At 0.5 ms of simulation time, eqn.(1) results in a poloidal velocity of approximately 10 km/s (moving in the electron diamagnetic direction).On the other hand, eqn.(2) results in a poloidal velocity of roughly 4 km/s (in the e − diamagnetic direction).Comparing these two, it would appear that the mode velocity of the linearly unstable n = 6, 8, 10, 12 modes is closer to that of resistive ballooning modes than to ideal ballooning modes.Further support for the identification of these modes as resistive modes comes from the fact that their growth rates become larger by increasing the resistivity.Similarly, by reducing the resistivity, the growth rates of the unstable modes decrease, as shown in the previous subsection.Therefore, the unstable high-n modes unstable in the present simulations are characterised as resistive peeling-ballooning modes.

Importance of extended MHD
Simulations without the diamagnetic effects have been performed in order to understand their influence onto the underlying instabilities described in the previous section.This is done for simulations with only one toroidal harmonic present, n = 8, and for different applied heating powers.It is observed that the simulations that include v * i have fundamentally different non-axisymmetric dynamics with respect to the simulations that neglect the diamagnetic effects.Figure 8 shows the evolution of the n = 8 magnetic energy in logarithmic scale of (a) simulations with and (b) without diamagnetic effects at four different values of P heat .
Increasing heating power causes a steepening of the temperature and, therefore, of the pressure at the plasma edge.For the simulations that include diamagnetic effects, the edge E r well at the pedestal becomes deeper with steeper pressure profiles.In said simulations, the PB modes become stabilised by the diamagnetic drift together with the E r (and its shear) [39,53,54].For the simulations without diamagnetic effects the pedestal steepens, but E r does not change.In fig.8(a) and (b), the heating power of the different simulations is changed at t = 0.33 ms to the values shown in the key.The simulations that consider diamagnetic effects observe γ n=8 to decrease as the heating power is increased.On the other hand, the simulations that neglect diamagnetic effects result in an increase of γ n=8 with increasing heating power.The results presented in this section emphasise the importance of including the diamagnetic flows for simulations of PB modes.
As mentioned at the beginning of this section, the pedestal profiles for the nominal heating power at these stages is stable to ideal peeling-ballooning modes.Higher pedestal pressure and/or edge current densities are required in order to reach a type-I ELM unstable scenario.The access to a type-I ELM unstable scenario appears to be closed without the inclusion of the two-fluid diamagnetic effects.Indeed this result was previously reported in Ref. [29] in the context of obtaining repetitive ELM cycle simulations, and was extended as a requirement to simulate type-I ELM cycles in Ref. [23].The following subsections are devoted to describing the fully non-linear phase (during which the n = 0 modes interact with each other and with the n = 0 axisymmetric background) of the simulation with n = 0, 2, 4, . . ., 12 at nominal heating and resistivity in the presence of diamagnetic effects, i.e., fig.4(a).

Non-linear phase
For the simulation that includes n = 0, 2, 4, . . ., 12, i.e., figure 4  The linear growth phase gives way to the early nonlinear growth phase until the amplitude of the perturbations becomes large with respect to the background plasma, at which point the non-linear phase begins.During the latter, a dynamic interplay between the n = 0 modes and the background plasma determines the instantaneous profiles observed in the simulations.Due to the persisting PB modes and the lack of a clear cyclical dynamics, the dynamics observed can be characterised as peeling-ballooning turbulence.

Filamentary transport
The non-axisymmetric time evolution of the ϕ = 0 outer midplane pressure gradient and the inner/outer divertor incident power are shown in fig.10(a) and (b), respectively, for 10 ms of simulation time.The incident power is defined as where q(t, , ϕ) is the heat flux at a given time in the divertor location at the toroidal angle ϕ, and R is the major radius.The colour map indicates that the pressure gradient and, therefore, the pressure profile in the outermost edge of the plasma is rapidly fluctuating.The corresponding fluctuations are governed by the non-axisymmetric modes that regulate the pedestal to fluctuate about a mean value, i.e., the PB turbulence.The incident power onto the divertors does not have characteristic spikes, but rather displays a quasi-continuous heat deposition, which is qualitatively similar to the QCE regime or the EDA H-mode in AUG. Figure 11 shows the ϕ = 0 outer midplane pressure in a colour map with logarithmic scale for a reduced time window of 0.4 ms, which is chosen between 4.8 and 5.2 ms, in order to show the dynamics of filamentary structures travelling outwards from slightly inside the separatrix (roughly 2 cm).The y-axis is the major radius, and the separatrix position is represented with a white line.The plasma blobs that travel outwards result from the resistive PB modes that are aligned to the magnetic fields and are moving in the electron diamagnetic direction.
Resistive peeling-ballooning modes that are destabilised below the ideal PB stability boundary regulate the pressure gradient about −250 kPa/m.This is made clearer with pressure and pressure gradient profiles taken in the representative time frame of As mentioned before, the peeling-ballooning modes do not behave in a cyclical fashion, but cause a quasi-continuous power deposition in the inner and outer divertor targets, as shown in fig.10(b).
The fluctuating profile in the last ∼ 7% of the confined region is clearly visible in fig.12(a).It shows how the resistive PB modes regulate the pedestal in such a way that the steepness of the profiles cannot grow to large values.This is why these simulations feature only small ELMs and not a mixed regime with small ELMs and type-I ELMs.Taking the time-varying temperature fluctuations at a single point in the steep gradient region in the outer midplane, (R, Z) = (2.14, 0.06), a spectrogram is performed.As a result, a dominant frequency in the range of 20 − 40 kHz is found, as can be observed in fig.13.A type-II ELMy H-mode in AUG with high triangularity and close to double null reported in Ref. [38] was described as having an electron pressure gradient oscillating about ∼ 150 kPa/m.The oscillating ∇p e was reportedly caused by MHD modes which were associated with electromagnetic fluctuations observed in a wide radial extent peaking in a frequency range of 30 − 50 kHz.Both observations hint at qualitative similarities to the simulation results described in this section.Nevertheless, it must be noted that the present simulations were performed in a different magnetic configuration, i.e., low triangularity and far from double null.Therefore, dedicated comparisons need to be performed in the future to produce quantitative comparisons between experiments and simulations.In particular, such comparisons will have to include variations of the plasma shape.

Divertor heat deposition
To show the quasi-continuous exhaust caused by the non-ideal peeling-ballooning modes excited near the separatrix, the electron temperature at the inner and outer divertor targets is plotted in fig.14(a   respectively.The target electron temperature is considered to be half of the plasma temperature and it is plotted for 10 ms of simulation time.The inner divertor target has a lower target temperature than the outer divertor.Similarly, the incident power to the inner divertor is lower than to the outer divertor, as seen in fig.10(b).There is a slight increase in the maximum target temperature (particularly visible in the outer target) as time progresses.This is due to the chosen heat source in the confined region which slowly increases the thermal energy content inside the separatrix.At any given time point, the heat deposition does not show significant variations in the toroidal direction.In other words, the heat deposition is roughly axisymmetric.It is important to note that the present simulations used only a simplified SOL transport model and, as such, the obtained heat distribution between targets will not necessarily reflect experimental observations.A more advanced SOL/divertor model is being developed in JOREK [50] and it will be used for future simulations.

Two simple paths to type-I ELMs
Based on the simulations presented in the previous section, the heating power is increased to understand the response of the resistive PB modes.Increasing heating power causes the edge temperature (and its gradient) to increase, which, in turn, causes the local resistivity to decrease (η ∝ T −3/2 ) and the pressure gradient and the diamagnetic drifts to grow larger.The stabilising influence of the diamagnetic effects and of E r (and its shear) onto PB modes becomes stronger and eventually completely stabilises them.At this point, the small ELM regime gives way to a type-I ELMy H-mode.Further details regarding this bifurcation determined by the applied heating power are presented in subsection 5.1.The transition from the small ELM regime to a type-I ELMy H-mode can also take place by sufficiently decreasing the separatrix density.Subsection 5.2 describes how decreasing n e,sep (with respect to the pure small ELM simulations) manages to completely stabilise the resistive PB modes and gives way to a type-I ELMy H-mode.The decreasing separatrix density prompts three important stabilising effects to take place: a smaller edge pressure gradient, faster plasma flows since v * i and v ExB are ∝ 1/n i , and a higher bootstrap current density.

Increasing heating power
The nominal heating power in JOREK units is 6.2 × 10 −6 (equivalent to ≈ 13 MW) and it was applied in the simulation shown in fig.20.The magnetic energies of the non-axisymmetric perturbations for the first 7 ms of said simulation are shown in fig.15(a).In the subsequent sub-figures, the heating power is progressively increased4 in small steps (the excess heating power is included from the beginning of each simulation).In fig.15(b), the non-linear behaviour does not show many differences to the simulation with nominal P heat .However, in the next figure, a transient phase where the n = 8 mode hosts most of the total non-axisymmetric energy, Σ n>0 E mag,n , is present.For this case with P heat ≈ 6.4 × 10 −6 (≈ 13.5 MW), the non-ideal PB modes become stabilised after roughly 10 ms.And in figs.15(d) and (e), Σ n>0 E mag,n is reduced until complete stabilisation.This mode stabilization allows the pedestal to build up and give rise to a type-I ELM crash eventually (discussed and shown in more detail later in this section).
Figure 15: Magnetic energies of the nonaxisymmetric perturbations for five different values of input heating power.The applied heating power increases progressively from (a)(≈ 13.0 MW) to (e)(≈ 13.9 MW).As a result of the increasing heating power, the magnetic energies of the n > 0 perturbations become completely suppressed if sufficient additional heating power is considered.
To further understand what governs the transition from small ELMs to type-I ELMs, the radial electric field at the outboard midplane is averaged between 1.0 − 2.0 ms and it is shown in figs.16(a)-(e).An interesting observation is that the three scenarios where the small ELMs become stabilised (c)-(e) have deeper radial electric field wells than the two cases that sustain the small ELMs (a) and (b).It is worth pointing out that even for the lowest heating power, the instantaneous radial electric field profiles at the outer midplane are often deeper than E r ∼ −15 kV/m, a representative value which has been associated to the L-H transition in AUG [55,56].This can be seen in fig.16(a), which shows in gray the instantaneous profiles used to obtain the time-averaged profile (black).The scans with and without diamagnetic effects shown in section 4.3 together with the observations presented so far in this section indicate that Figure 16: Time-averaged profiles of the outboard midplane radial electric field for five different values of input heating power.The applied heating power increases progressively from (a)(≈ 13.0 MW) to (e)(≈ 13.9 MW).
small ELMs feel a stabilising effect from larger diamagnetic drifts and the resulting deeper radial electric field well.
In order to directly show the bifurcation from small ELMs to type-I ELMs, the heating power of the simulation described in section 4.4 is increased at 5.6 ms of simulation time from ≈ 13.0 to ≈ 13.9 MW.Resulting from the heating power increase, the resistive PB modes start to become smaller in amplitude and their radial extent starts to reduce.This process takes roughly 4 ms to complete and, thereafter, a steeper pedestal is allowed to form.Ultimately, the steepening pedestal crosses the type-I ELM stability boundary and an ELM with dominant toroidal mode numbers n = 2 and 4 is excited.The process described in this paragraph can be evidenced in fig.17(a) and (b), which respectively show the ϕ = 0 outboard midplane pressure gradient and the power incident on the inner and outer divertors.The small jump at t ≈ 12 ms takes place because the parallel heat conductivity had been increased roughly to the Spitzer-Härm values in this simulation, but the onset of the type-I ELM takes place regardless of said change.The incident power that reaches the divertor targets is significantly increased when the type-I ELM crash appears.Comparatively, it is clear that the small ELMs cause much weaker heat fluxes to the divertor targets.To directly show the influence of the increased heating power onto the resistive PB modes that cause small ELMs, the plasma pressure in real space together with the position of the separatrix are plotted for a restricted time frame between t = 5 and 10 ms in fig.18.The expelled filaments after the heating power increase seem to have smaller amplitudes and they travel for shorter distances.They eventually disappear completely.
In section 2, we discussed the response of AUG small ELMs at high separatrix density, low triangularity and high edge safety factor towards increasing heating power.Namely, in such experiments the small ELMs can become suppressed with sufficiently high additional heating power.The situation Figure 18: Time evolution of the outboard midplane pressure at ϕ = 0 (a) upon increasing the heating power at 5.6 ms.The filaments expelled from the confined region become weaker when the heating power is increased.They ultimately completely disappear and the pedestal is able to grow further. is similar for the EDA H-mode.For both cases, the role of plasma shaping in the suppression of small ELMs/QCM seems to be pivotal.In our simulations, which feature high n e,sep , low triangularity, and high q 95 , the transition from a regime dominated by small ELMs towards a type-I ELM is obtained by suddenly increasing the input heating power in small ELM simulations.The small ELMs start to weaken and the filaments formed by the small ELMs are gradually reduced in amplitude until disappearing completely.p ped rises with increasing P heat due to the increase of T ped ; n e,ped remains unchanged in the first few milliseconds after the heating power was increasedit only starts rising when the particle transport by small ELMs becomes significantly reduced.Additionally, due to the larger ∇p, the radial electric field well at the plasma edge deepens, and the bootstrap current density starts to rise.Taking one millisecond time-averages, the outer midplane profiles are tracked during the pure small ELMs phase (from 2.6 ms until 5.6 ms) and during the transition phase where the resistive PB modes start to disappear (from 5.6 ms to 9.6 ms) and are plotted in fig.19.
Three time-averaged profiles correspond to the original small ELM phase and show a roughly con-Figure 19: Evolution of the time-averaged ϕ = 0 outboard midplane pressure (a), radial electric field (b), and toroidal current density (c).The time averaging is done for one millisecond intervals.The three time-averaged profiles before the heating power is increased have a shallow E r , and the four profiles after P heat show systematically higher p ped and deeper E r well (as indicated by the black arrow for the latter).stant p ped , a weak E r at the edge, and low toroidal current density.On the other hand, the four profiles at higher heating power show systematically higher p ped , deeper E r well, and a broader and higher edge current density.The last time-averaged profiles in gray, with the highest p ped , min(E r ) ≈ −28 kV/m and a high toroidal current density, is taken during a phase that mostly has suppressed the resistive PB modes.The destabilising effect of the increased ∇p is tied to the stabilising effects of the deepening of the E r well and of the toroidal current density.The additional stabilising effect of the edge resistivity decreasing as the pedestal top temperature increases is also important at this stage.Therefore, the disappearance of the resistive PB modes appears to be due to the E r deepening, −j ϕ increasing, and the decreasing local η.

Decreasing separatrix density
The previous section detailed the bifurcation from a small ELM-dominant regime to a type-I ELM by means of increasing the heating power.Another path towards type-I ELMs, starting from a small ELM regime, is to decrease the separatrix density.In the experiment, this can be achieved by reducing, or completely removing, the particle source given by a gas puff (replacing the particle flux by means of cryogenic deuterium pellet injection can keep n e,ped approximately unchanged).Indeed, it has been reported that separatrix densities below ∼ 0.35n GW do not host dominant small ELM phases [7].
The simulations presented thus far had a separatrix density of n e,sep ≈ 3 × 10 19 m −3 .To reduce the separatrix density in our simulations, we reduce the edge density source; in order to maintain n e,ped unchanged we increase the core density source.Resulting from the lower n e,sep , the pressure gradient is locally reduced, and it is overall shifted slightly inwards-such response of the pressure profile position is also observed in experiments [57,58].The decrease of separatrix density (at fixed n e,ped ) also causes a deeper E r well (because E r ∝ 1/n e ) and a higher bootstrap current density (because the density gradient increases).
An important caveat must be mentioned: many physical effects related to the pedestal position (particularly related to neutrals penetration) are not included in the JOREK model used for these simulations and, therefore, the qualitative agreement will likely not translate to a quantitative agreement at this stage.
A new simulation is then set-up with lower n e,sep = 2 × 10 19 m −3 .The non-zero toroidal modes are included after 0.1 ms (exactly the same time as the small ELM simulations described in section 4).The magnetic energies of the high and low n e,sep simulations are shown in fig.20(a

Conclusions
H-mode operation without large type-I ELMs is an imperative requirement for ITER in highperformance conditions.To this purpose, naturally ELM-free H-modes and ELM mitigated/suppressed regimes are considered and actively researched.In AUG, several such alternatives are under investigation; two of them are the quasi-continuous exhaust (QCE) regime and the enhanced D-alpha Hmode.Both can be operated completely without type-I ELMs.The pedestal is limited by small ELMs for the QCE regime and by a quasi-coherent mode for the EDA H-mode.These transport mechanisms quasi-continuously expel heat and particles from the confined region.It is presently unclear whether or not it will be possible to operate such regimes in ITER.Simulations performed with JOREK, which show several key features of small ELMs have been presented in this paper.Resistive peeling-ballooning modes near the separatrix are identified as the transport mechanism underlying such small ELMs.
Modelling the pedestal build-up at fixed pedestal width, with stationary diffusion coefficients and sources, resistive peeling-ballooning modes that regulate the pedestal below the ideal PB stability boundary are observed under appropriate conditions.Namely, simultaneously high separatrix density and not too much heating power.The necessary conditions to sustain sufficient outwards transport by small ELMs is primarily determined by the separatrix density and the input heating power.In particular, simulations with high n e,sep and low heating power observe phases (longer than 10 ms) with quasi-continuous outwards transport that prevent the pedestal from reaching a type-I ELM unstable scenario.The resistive nature of such PB modes is determined by the fact that they appear below the ideal PB stability boundary, because their growth rates are largely reduced/enhanced by decreasing/increasing resistivity, and because their poloidal mode velocities are measured to be close to that expected for resistive modes.
An important ingredient required in order to properly simulate these resistive PB modes is the inclusion of diamagnetic effects, which (in the simulations) allow the radial electric field well to develop in the pedestal region.In the absence of diamagnetic effects, it is not possible to stabilise the small ELMs by increasing the heating power.In contrast, when diamagnetic effects are included, the small ELMs become completely stabilised and the plasma state moves to a type-I ELMy H-mode by increasing P heat .Similarly, decreasing the separatrix density completely stabilises the small ELMs if diamagnetic effects are included.Another important effect that should be included when modelling these instabilities is the bootstrap current density because it has a stabilising influence onto high-n peelingballooning modes.At the moment, JOREK evolves the bootstrap current density through the Sauter formula [40,41], as explained in section 3.However, the Sauter expression is known to be inaccurate depending on the parameter regime, particularly at high collisionality [59].Therefore, an improvement of the bootstrap current density source in JOREK will be pursued in the future.Finally, the simplified resistivity with only Spitzer temperature dependency used in JOREK should be improved to include the influence of neoclassical effects and effective main ion charge greater than unity.
Commission.In particular, contributions by EUROfusion work packages Enabling Research (EnR) and Medium Size Tokamaks (MST) is acknowledged.The simulations presented in this work were performed on the Marconi-Fusion supercomputer operated by CINECA in Italy.

Figure 2 :
Figure 2: Initial conditions for the electron density and plasma temperature (T = T e + T i ) (a), and plasma pressure (p = p e + p i ) and toroidal current density (b) in the outboard midplane.

Figure 3 :
Figure 3: Outboard midplane pedestal profiles of electron density (a), plasma temperature (b), radial electric field (c), and flux-surface averaged toroidal current density (d) during the imposed pedestal build-up.The first millisecond is defined by a strong steepening of the pedestal, and the pedestal top grows progressively at a fixed gradient.

Figure 4 :
Figure 4: Magnetic energies of the n = 0 modes in logarithmic scale in the first 0.8 ms.(a) shows the simulation with n = 0, 2, 4, . . ., 12 and (b) shows the simulation with n = 0, 2, 4, . . ., 20.This comparison shows that the dominant mode numbers in the linear phase are n = 10 and 12 in both cases; into the nonlinear phase (not shown), the n > 12 modes are subdominant.

Figure 6 :
Figure 6: Growth rates of different toroidal mode numbers for the resistivity scan.Increasing the resistivity leads to the destabilisation of resistive PB modes.
move with time it is possible to determine their poloidal velocity.The modes shown in the figure are the result of perturbations with different toroidal mode numbers; the magnetic flux fluctuation is defined as ψ = nmax n>0

Figure 7 :
Figure 7: (a) Flux surfaces of ψ N = 0.92 and 0.99 with their corresponding colour coded arc lengths.(b)-(g) The variation of the poloidal magnetic flux perturbations along the arc length for the different toroidal mode numbers.The magnetic flux perturbations for n = 2 and 4 do not have any coherent structures as they are linearly stable to the profiles at this point in time (t = 0.5 ms).

Figure 8 :
Figure 8: Evolution of E mag,n=8 for simulations with (a) and without (b) the inclusion of ion pressure gradient-driven diamagnetic flows.Four different input heating powers are considered.Therefore, eight single−n simulations are shown.Increasing heating power in simulations that include diamagnetic flows shows an important stabilisation of the n = 8 PB mode.When neglecting the diamagnetic effects, on the other hand, increasing heating power causes the unstable PB mode to grow even faster due to the steeper pressure profiles.
(a), the magnetic and kinetic energies of the non-axisymmetric modes are shown in fig. 9 in linear scale for 10 ms of simulation time (a) and (c) and in logarithmic scale for the first 2 ms (b) and (d).

Figure 9 :
Figure 9: Magnetic and kinetic energies of the nonaxisymmetric modes in linear scale for 10 ms of simulation time (top) and in logarithmic scale for the first 2 ms of simulation time during which the linear growth phase takes place.

Figure 10 :
Figure 10: Time evolving outer midplane edge pressure gradient at ϕ = 0 in colour scale (a), and inner/outer divertor incident power (b).The varying pressure profile is caused by quasi-continuous outward transport created by non-ideal peelingballooning modes.

Figure 11 :
Figure 11: Time evolving outer midplane edge pressure at ϕ = 0 in logarithmic scale.Plasma blobs travelling outwards are aligned to the magnetic field lines.

Figure 12 :
Figure 12: Pressure (a) and pressure gradient (b) profiles in the time window t = [4.0,6.0] ms together with a time-averaged profile in black.The timeaveraged profile shows a 'staircase' structure with a large pressure gradient in the vicinity of the separatrix.

Figure 13 :
Figure 13: Time evolving (a) and averaged (b) frequency spectrogram of the temperature fluctuations in (R, Z) = (2.14, 0.06).Dominant frequencies in the range 20 − 40 kHz can be observed in both cases.

Figure 14 :
Figure 14: Time evolution of the inner (a) and outer (b) target electron temperature (T e = T /2) caused by the resistive PB modes excited at the very edge of the plasma.The inner target has a lower temperature than the outer target as well as a lower incident power.

Figure 17 :
Figure 17: Time evolution of the outboard midplane pressure gradient at ϕ = 0 (a) upon increasing the heating power at 5.6 ms.And the power incident on the inner and outer divertor targets (b) resulting from small ELMs (t < 10 ms) and from a type-I ELM (t 19 ms).

Figure 20 :
Figure 20: Magnetic energies of the nonaxisymmetric modes in logarithmic scale in the first 4 milliseconds of simulation time for (a) small ELMs at high n sep (∼ 3 × 10 19 m −3 ) and (b) their response to lower separatrix density (∼ 2 × 10 19 m −3 ).
) and (b), respectively.The linear phases are similar between the two cases, with growing n = 8, 10, 12 high-n (peeling-)ballooning modes.But the simulation with low n e,sep (b) deviates as the n < 10 and n = 12 become completely stabilised, and the n = 10 only reaches small amplitudes and does not affect the n = 0 background.The stabilisation of the modes with higher toroidal mode numbers leads to a single toroidal mode number with very small amplitude that does not cause any changes to the background plasma.The simulation at lower separatrix density, fig20(b), sees the pedestal evolve, but it was not continued until a type-I ELM is reached to save computing time.Time-averaged outer midplane pro-

Figure 21 :
Figure 21: Time-averaged outboard midplane profiles of the pressure gradient and E r during the linear growth phase for nominal separatrix density (a) and (c), and for lower separatrix density (b) and (d).