# Probing non-linear MHD stability of the EDA H-mode in ASDEX Upgrade

A. Cathey1, M. Hoelzl1, L. Gil2, M.G. Dunne1, G.F. Harrer3, G.T.A.
Huijsmans4,5, J. Kalis1, K. Lackner1, S.J.P. Pamela6, E. Wolfrum1, S. Günter1,
the JOREK team111see the author list of M. Hoelzl et al 2021 Nucl. Fusion 61
065001, the ASDEX Upgrade Team222See the author list of U. Stroth et al 2022
Nucl. Fusion 62 042006, and the EUROfusion MST1 Team333see the author list of
B. Labit et al. 2019 Nucl. Fusion 59 0860020 1 Max Planck Institute for Plasma
Physics, Boltzmannstr.2, 85748 Garching, Germany  2 Instituto de Plasmas e
Fusão Nuclear, Instituto Superior Técnico, Universidade de Lisboa, 1049-001
Lisboa, Portugal  3 Institute of Applied Physics, TU Wien, 1040 Vienna,
Austria  4 CEA, IRFM, 13108 Saint-Paul-Lez-Durance, France  5 Eindhoven
University of Technology, P.O. Box 513, 5600 MB Eindhoven, The Netherlands  6
CCFE, Culham Science Centre, Abingdon, Oxon, OX14 3DB, United Kingdom
<EMAIL_ADDRESS>
###### Abstract

Regimes of operation in tokamaks that are devoid of large ELMs have to be
better understood to extrapolate their applicability to reactor-relevant
devices. This paper describes non-linear extended MHD simulations that use an
experimental equilibrium from an EDA H-mode in ASDEX Upgrade. Linear ideal MHD
analysis indicates that the operational point lies slightly inside of the
stable region. The non-linear simulations with the visco-resistive extended
MHD code, JOREK, sustain non-axisymmetric perturbations that are linearly most
unstable with toroidal mode numbers of ${n=\\{6\dots 9\\}}$, but non-linearly
higher and lower $n$ become driven and the low-n become dominant. The poloidal
mode velocity during the linear phase is found to correspond to the expected
velocity for resistive ballooning modes. The perturbations that exist in the
simulations have somewhat smaller poloidal wavenumbers (${k_{\theta}\sim
0.1~{}\mathrm{to}~{}0.5~{}\mathrm{cm^{-1}}}$) than the experimental
expectations for the quasi-coherent mode in EDA, and cause non-negligible
transport in both the heat and particle channels. In the transition from
linear to non-linear phase, the mode frequency chirps down from approximately
35 kHz to 13 kHz, which corresponds approximately to the lower end of
frequencies that are typically observed in EDA H-modes in ASDEX Upgrade.

††:

##  1 Introduction

Tokamak operation in the so-called high confinement mode (H-mode) naturally
exhibits the onset of repetitive edge localised modes (ELMs), which can expel
${5}$ to ${15~{}\%}$ of the thermal energy stored in the magnetically confined
plasma [1, 2]. For present-day devices, the energy expelled by ELMs is not a
cause for concern as it does not exceed the material limits of the plasma
facing components. However, for future devices like ITER, which will confine
much larger amounts of thermal energy, it is foreseen that large ELMs will not
be tolerable [3, 4]. Techniques that mitigate or suppress ELMs are actively
investigated, but it is presently not possible to determine operational spaces
of applicability in future machines. The same problem exists for regimes of
operation that are either completely devoid of large ELMs or show only small
ELMs that may be tolerable by the plasma facing components. Such no- and
small-ELM regimes host some transport mechanism that prevents the edge
pressure gradient and current density to grow unconstrained (which constitutes
the reason why ELMs become excited) and flushes unwanted impurities out of the
confined plasma [5, 6, 7].

One specific example of no-ELM regimes is the enhanced D-alpha (EDA) H-mode,
which was first observed in the high-field tokamak, Alcator C-mod, with ion-
cyclotron heating after a fresh boronisation [8, 9]. The EDA constantly
exhibits a quasi-coherent mode (QCM) primarily localised in the low-field side
with a radial width of ${\sim 2-5~{}\mathrm{mm}}$ and frequencies of
${50-150~{}\mathrm{kHz}}$. The QCM which appears to be located in the steep
gradient region of the pedestal (at the radial electric field minimum [10],
but further outside in Ref. [11] which probes, and may influence, the plasma).
The influence of the QCM is observed in fluctuations of the density,
electrostatic potential, and poloidal magnetic field [12, 13, 11]. As such,
the QCM is thought to be an instability with electromagnetic character that
constantly regulates the pedestal below the conditions where type-I ELMs are
excited by causing continuous particle transport.

In addition to C-mod, the EDA has also been observed recently in EAST (closely
related to a small ELM regime) [14], in ASDEX Upgrade with electron-cyclotron
heating (ECRH) [15], and with a mixture of ECRH and neutral beam injection
(NBI), including Argon seeding [16]), and also in DIII-D [17]. For the latter
two machines, which have similar sizes and current/magnetic field and are
larger than C-mod, the QCM frequency
${f_{\mathrm{QCM}}=15-40~{}\mathrm{kHz}}$, and the poloidal wavenumbers range
between ${k_{\theta}=0.5-0.6~{}\mathrm{cm^{-1}}}$, while in C-mod they range
between ${0.8-2.0~{}\mathrm{cm^{-1}}}$ [18]. These differences can be expected
due to the different machine sizes; the toroidal mode numbers in question
correspond to ${n\approx 20}$. In all devices, the QCM moves in the electron
diamagnetic direction in the lab frame444In Alcator C-mod, mode propagation
with respect to the plasma frame of reference (i.e., ${v_{\mathrm{ExB}}}$) has
been separately reported in the ion diamagnetic direction for ICRF-heated
plasmas as measured with gas puff imaging [10] and in the electron diamagnetic
direction for ohmically-heated plasmas as measured with mirror Langmuir probes
[11]. and its frequency starts at larger values (roughly double), chirp down
and remain in the aforementioned ranges for several confinement times.
Magnetohydrodynamic (MHD) analysis of EDA indicates that the operational
points lie slightly inside of the stability boundary for peeling-ballooning
modes (the instabilities responsible for ELMs) with medium/high toroidal mode
numbers ($n=\\{10\dots 30\\}$) and has been hypothesised that the QCM possibly
relates more closely to resistive ballooning modes [19, 20, 21, 22] or to
drift ballooning modes [12].

The theoretical framework that relates to the QCM in EDA has not been yet
fully explained [10]. Approaches from modelling include BOUT simulations for
Alcator C-mod which connect resistive ballooning modes to the QCM in terms of
several distinct features (poloidal wavenumber, frequency, velocity) [20].
BOUT++ has also been applied to study the QCM in C-mod and found good
agreement with experiments; such studies suggest that resistive ballooning
modes and drift-Alfvèn waves are the dominant instabilities during the EDA.
The gyrokinetic code, GENE, has been used to study an EDA in ASDEX Upgrade,
which did not find clear indications of the QCM (possibly due to difficulties
encountered in resolving flux surfaces too close to the separatrix) [23].
Based on the same experimental equilibrium on which the GENE simulations had
been based, first non-linear extended MHD simulations with the JOREK code [24,
25] are presented.

In this paper, non-linear extended MHD simulations with initial conditions
chosen from an EDA discharge in ASDEX Upgrade are presented. The description
of the experimental discharge and the set-up of the simulations are presented
in section 2, which includes results from linear ideal MHD calculations for
the experimental equilibrium. Thereafter, in section 3 the linearly unstable
modes in the visco-resistive JOREK simulations are described in terms of their
structure, growth-rates, mode velocity and direction of propagation and,
finally, the resulting non-linear coupling that excites the linearly stable
modes. The dominant modes during the non-linear phase (when the non-
axisymmetric perturbations interact with the background plasma), the changes
caused onto the magnetic fields, and the nature of the fluctuations across the
pedestal are detailed in section 4. Finally, section 5 summarises and
discusses the results presented.

##  2 Experimental discharge and simulation set-up

The simulations presented here take their initial conditions from AUG
discharge #36330 in the time range 6.115 to 6.190 s, which is an argon-seeded
EDA H-mode first described in Ref. [16], and that has been simulated with GENE
[23]. This section presents several details on the experimental discharge as
well as newly analysed data from the magnetic pick-up coil signals in section
2.1, and the JOREK simulation set-up together with relevant details from the
model are discussed in section 2.2.

###  2.1 AUG #36330 – Ar seeded EDA H-mode

The discharge used to set-up the simulations presented in this paper featured
a stable EDA H-mode for several seconds. Detailed accounts of this
experimental discharge can be found in Refs [16, 23]; here is only a brief
summary of the characteristics relevant for the present modelling work. The
discharge features a lower single null plasma with elongation ${\kappa=1.6}$
and high triangularity (${\delta_{\mathrm{lower}}=0.49}$ and
${\delta_{\mathrm{upper}}=0.31}$, ${\delta_{\mathrm{avg}}=0.40}$) operating in
favourable ${\nabla B}$ drift direction (where the ion $\bm{B}\times\nabla B$
drift points towards the active X-point). The toroidal magnetic field is
$-2.5~{}\mathrm{T}$, the plasma current ${0.8~{}\mathrm{MA}}$, and the edge
safety factor ${q_{95}=4.86}$. Figure. 1 shows the traces of several
quantities including the line-averaged density for lines-of-sight that probe
the core and the edge (a), the divertor currents (b) where it can be seen that
ELMs are not present, a frequency spectrogram from a helium beam signal (c)
with a quasi-coherent mode (QCM) with fluctuation frequency of $\sim\\{20\dots
40\\}~{}\mathrm{kHz}$, and the heating mix (d). The latter is comprised of
ECRH and NBI with ${2.5~{}\mathrm{MW}}$ each from 5.9 to 6.3 s. The JOREK
simulations are initialised from a time point within this phase.

Figure 1: Traces of the AUG#36330 discharge for the time frame 5.3 to 6.3 s.
(a) shows the core and edge line-integrated density measurements from
interferometry, (b) the inner and outer divertor currents. The frequency
spectrogram from a helium beam signal displays a dominant frequency around
${f_{\mathrm{QCM}}\approx 30~{}\mathrm{kHz}}$ in (c); there are phases without
signal because the helium valve is turned on and off periodically in order to
subtract the background emission. And the heating power mixture together with
the radiated power is shown in (d).

The absence of ELMs observed in fig. 1(b) is likely attributed to the
electromagnetic QCM with ${f_{\mathrm{QCM}}\approx 30~{}\mathrm{kHz}}$ (ranges
between ${20-40~{}\mathrm{kHz}}$), which can be observed in fig. 1(c). The
transport caused by the cross-field transport, including the QCM, causes a
peak heat-flux onto the divertor targets of ${5-10~{}\mathrm{MW/m^{2}}}$
(derived from Langmuir probes) [16].

The high-resolution equilibrium reconstruction with bootstrap constraint is
generated with CLISTE [26] from data in the time range 6.115 to 6.190 s, and
it is used to initialise the JOREK simulations that will be discussed in depth
in the following sections. Linear ideal MHD stability analysis with MISHKA
[27] is performed for this equilibrium, and shown in fig. 2. The linear
stability analysis finds the operational point to be slightly inside of the
(${n\approx 15}$) ballooning boundary. However, the MISHKA simulations neglect
all non-ideal effects, notably the destabilizing resistivity and the
stabilizing ExB and diamagnetic flows. Non-linear [28, 29] and linear [30]
resistive MHD simulations have shown to move the ballooning boundary to lower
pressure gradients, such that the operational point could easily be unstable
to resistive (peeling-)ballooning modes. Indeed, the JOREK simulations that
will be presented in the following sections find that the equilibrium is in
fact unstable to resistive peeling-ballooning modes at realistic plasma
resistivity.

Figure 2: MISHKA ideal linear MHD analysis of the equilibrium reconstruction
at ${6.165~{}\mathrm{s}}$. The operational point at ${6.165~{}\mathrm{s}}$ is
found to be just inside of the ideal ballooning (${n\approx 15}$) boundary.

In addition to the QCM fluctuations that are observed in the helium beam
diagnostic (which can also be observed in spectrograms of reflectometry,
interferometry, electron cyclotron emission, magnetic pick-up coils closest to
the plasma [15]), magnetic fluctuations are measured during the EDA H-mode
phase. By considering several magnetic pick-up coils in the equatorial
midplane, but at different toroidal angles, it is possible to identify the
relevant toroidal mode numbers of ${n=\\{3\dots 7\\}}$. The cross-phaseogram
that allows such analysis is shown in fig. 3, but the precise nature of these
modes, and their relation (if any) to the QCM is unclear as of yet. For this
particular case, the frequencies of the bands with distinct toroidal mode
numbers can be written as ${f_{n}=n\times 32.9~{}\mathrm{kHz}}$. However, it
must be noted that in several EDA discharges which display visible bands in
the magnetic pick-up coils cross-phaseogram the difference between the base
frequency with the QCM frequency can be ${\lesssim 15~{}\mathrm{kHz}}$.

Figure 3: Cross phaseogram of several signals from magnetic pick-up coils in
the equatorial midplane at different toroidal angles. Different colours
indicate distinct toroidal mode numbers. In particular, $n=3\dots 7$ are
observed to be the relevant toroidal harmonics.

###  2.2 JOREK simulation set-up

As mentioned before, the JOREK simulations are initialised from an
experimental reconstruction corresponding to ${6.165~{}\mathrm{s}}$ in AUG
discharge #36330. The corresponding simulations are carried out with the
single-temperature reduced MHD model within JOREK and considers a perfectly
conducting wall as boundary condition. Note that reduced MHD can be applied
here allowing to reduce computational costs, since excellent agreement of
reduced and full MHD models for pedestal applications has been shown [31]. For
all simulations, the diamagnetic drift and bootstrap current source extensions
are considered. The latter means that when profiles change non-linearly, the
bootstrap current fraction changes according to the Sauter analytical model
[32, 33]; if profiles remain stationary, the bootstrap current fraction
remains stationary as well. Using the former means that the radial electric
field well, which is characteristic of experimental $E_{r}$ profiles in H-mode
and proportional to ${\nabla p/n}$, is present in the simulations as has been
shown in previous work [34]. A detailed account of the JOREK code, and of the
single temperature reduced MHD model can be found in Ref. [24].

The input profiles of the density, temperature, radial electric field, and
toroidal current density at the outboard midplane (flux-surface average for
${j_{\varphi}}$) are shown in fig. 4 in full black lines. Black crosses
represent profiles after ${10~{}\mathrm{ms}}$ of axisymmetric evolution (i.e.,
in the absence of perturbations) are shown. The radial electric field profile
in full green line (which only spans until 1.00) comes from charge exchange
recombination spectroscopy measurements for the impurity toroidal and poloidal
velocity of and the force balance equation ${E_{r}=\nabla
p_{\alpha}/(eZ_{\alpha}n_{\alpha})-\bm{v}_{\alpha}\times\bm{B}}$, where
$\alpha$ is the ${\mathrm{N}^{7+}}$ impurity species. The uncertainty of the
$E_{r}$ minimum is roughly ${\sim\pm 4~{}\mathrm{kV/m}}$ [23]. It is evident
that the axisymmetric profile does not match closely to the measured profile,
which is partly because the present simulations do not include a source for
the toroidal rotation caused by the NBI torque and partly because the single
temperature treatment cannot account for lower ion pressure gradients than
electron and drives $E_{r}$ to be more negative in the pedestal middle.
However, when the non-axisymmetric simulation (${n=0\dots 13}$) develops
(dashed purple line) the minimum value of $E_{r}$ matches with the
experimental measurements due to modifications to the profiles caused by the
perturbations which will be described in detail in section 4.

Figure 4: The input profiles of density, temperature, radial electric field,
and toroidal current density. The latter is comprised of ohmic and bootstrap
contributions and it is a flux-surface averaged profile, while the rest are
profiles at the outboard midplane. Full black lines are the initial profiles
and black crosses represent profiles after $10~{}\mathrm{ms}$ of axisymmetric
simulation. A profile fitted from CXRS measurements for the $E_{r}$ is shown
in full green line from ${\psi_{\mathrm{norm}}=0.85}$ to 1.00, and a
toroidally-averaged midplane profile from a non-linear simulation (${n=0\dots
13}$) in a dashed purple line.

The axisymmetric evolution of the profiles is determined by diffusion
coefficients and sources of particles and heat, which were designed to
maintain the initial temperature, density, and current density profiles
constant, as depicted by the fact that the full black lines and the black
crosses in fig. 4 overlap. The physical meaning of the diffusion coefficients
corresponds to neoclassical and turbulent transport that cannot be simulated
with JOREK, but is included with these ad-hoc profiles.

The grid resolution required to resolve the simulations comprises 214 in the
radial direction (i.e., from the axis to the perfectly conducting wall) which
are accumulated primarily in the pedestal region such that the radial extent
of one grid element is ${1.28~{}\mathrm{mm}}$, 354 in the poloidal direction,
and 64 toroidal planes, which are needed to include $14$ toroidal Fourier
harmonics in the simulation: ${n=0\dots 13}$. These parameters make it
possible to include a realistic parallel-to-perpendicular heat diffusion
anisotropy of the order ${10^{9}}$. The temporal resolution in JOREK is not
constrained by the CFL criterion since an implicit time stepping scheme is
used; the time step used for the simulations corresponds to
${\sqrt{\mu_{0}\rho_{0}}=0.6687~{}\mathrm{\mu s}}$ (except in the non-linear
phase when it is sometimes necessary to decrease it by half to improve
convergence of the iterative solver). Since the MHD model is valid for
perturbations of low frequency (relative to the ion cyclotron frequency,
${\Omega_{i}=eB/m\sim 10^{8}~{}\mathrm{1/s}}$), perturbations that are faster
cannot be adequately resolved with JOREK even if the time step were set to
${\delta t\lesssim 1/\Omega_{i}\sim 0.01\sqrt{\mu_{0}\rho_{0}}}$.

Having described the experimental discharge, the simulation set-up and initial
conditions, the focus is turned to describing the results during the linear
phase.

##  3 Linear growth phase and dependencies

The present section describes the simulation results during the linear growth
phase including linear growth rates, and the dynamics of the non-axisymmetric
perturbations during the early non-linear phase, where three-wave interactions
allow linearly stable modes to become non-linearly excited. The location and
velocity of the modes during the linear phase is further discussed and their
relation to resistive ballooning modes is highlighted.

###  3.1 Growth rates and mode localisation

After ${1~{}\mathrm{ms}}$ of axisymmetric simulation (during this time the
parallel and poloidal flows establish and saturate), non-axisymmetric
perturbations with toroidal Fourier harmonics ${n=1\dots 13}$ are initialised
at noise-amplitude. Figure 5 shows the non-axisymmetric perturbation magnetic
energies (top) and their growth rates during the linear and early non-linear
phases (bottom). It can be distinguished that there are linearly unstable
modes (${n=\\{5\dots 11\\}}$) and linearly stable modes (${n=\\{1\dots
4\\},~{}12}$, and 13) that become non-linearly excited during the early non-
linear phase through three-wave interactions. Namely, a mode $n_{3}$ is driven
by $n_{1}$ and $n_{2}$ provided that ${n_{3}=n_{1}\pm n_{2}}$ and the
resulting growth-rate corresponds to
${\gamma_{n_{3}}=\gamma_{n_{1}}+\gamma_{n_{2}}}$ [35]. In this case several
pairs of modes can contribute to driving a linearly stable mode. For example,
${n_{1}=6}$ and ${n_{2}=6}$ drive ${n_{3}=12}$, which is also driven by
${n_{1}=5}$ and ${n_{2}=7}$ and any other combination that satisfy
${n_{3}=n_{1}\pm n_{2}}$. The resulting growth-rate of $n=12$ is given by
contributions from all pairs of modes that satisfy ${n_{1}\pm n_{2}=12}$, but
it is mostly determined by the dominant mode pair that drive it.

Figure 5: Top: magnetic energies of the non-axisymmetric perturbations during
the linear phase (${\lesssim 1.5~{}\mathrm{ms}}$) and early non-linear phase
(${\gtrsim 1.5~{}\mathrm{ms}}$). Bottom: growth-rates for the different
toroidal mode numbers during the linear phase (black squares) and during the
early non-linear phase (purple circles). Linearly stable modes become non-
linearly destabilised due to three-wave interactions.

###  3.2 Mode location and velocity

The location of the linearly unstable modes is around
${\psi_{\mathrm{norm}}\approx 0.98}$, which corresponds to the maximum
temperature gradient and is consistent with experimental measurements from
Alcator C-mod [10], and with recent observations from the EDA in ASDEX Upgrade
with Helium beam [36]. To clearly illustrate this, fig, 6 (top) shows the
derivatives of the electron temperature, pressure, and density with full black
line, dashed purple line, and green line with marks, respectively. And fig. 6
(bottom) shows the absolute value of the dominant complex ${m/n}$ Fourier
coefficients in blue for ${n=6}$, red for $7$, and black for $8$. With thicker
lines, the resonant surface of ${q=6}$ is shown for each harmonic with the
respective colours. Unsurprisingly, the dominant poloidal mode number during
the linear phase is ${m=\\{35\dots 56\\}}$. The location of the mode maximum
is shown with a vertical black line that spans both plots and corresponds to
${\psi_{\mathrm{norm}}=0.983}$, which sits almost exactly at the
${\nabla_{R}(T_{e})}$ maximum.

Figure 6: Top: derivatives (with respect to the major radius) at the outer
midplane for $T_{e}$, $p_{e}$, and $n_{e}$ profiles; the maximum gradient is
located at ${\psi_{\mathrm{norm}}=0.980}$, 0.977, and 0.995, respectively.
Bottom: absolute value of various $m/n$ Fourier coefficients for ${n=6}$ in
blue, ${7}$ in red, and $8$ in black. The respective thick lines show the
${m/n=36/6}$, $42/6$, and $48/6$; namely at the $q=6$ rational surface.

Figure 7 shows the poloidal velocity of the perturbations along the flux-
surface where the modes are located (${\psi_{\mathrm{norm}}\approx 0.980}$)
during the linear phase (${t=1.60~{}\mathrm{ms}}$) with green circles and
during the non-linear phase (${t=2.70~{}\mathrm{ms}}$) with black squares.
Additionally, the plasma rotation velocity
(${v_{\mathrm{ExB}}+\bm{v}_{\parallel}\cdot\bm{b}_{\theta}}$) is also shown
for either phase with a green dashed line and a black full line, respectively.
The reduction in mode velocity between linear and non-linear phases will be
discussed further in section 4. The movement of the non-axisymmetric
perturbation is tracked along individual flux-surfaces and the distance
travelled is used to obtain the poloidal mode velocity. From
${\psi_{\mathrm{norm}}=0.929}$ until ${0.987}$ (which is the outermost flux-
surface with a reliable measurement of the velocity), the mode velocity at the
outboard midplane is roughly ${-15~{}\mathrm{km/s}}$, where the negative sign
indicates movement in the electron diamagnetic drift direction. In the
laboratory frame, the modes travel in the electron diamagnetic direction
(negative velocities in fig. 7) at all radial locations. In the plasma frame
(i.e., the difference between the green circles and dashed line) the modes
move in the ion diamagnetic direction at the outboard midplane with
${v_{\mathrm{mode,pl.}}\approx 1~{}\mathrm{km/s}}$ only close to the location
where the mode amplitude is maximised (from ${\psi_{\mathrm{norm}}\approx
0.970}$ until $0.982$). Inward and outward of these locations, the modes move
in the electron diamagnetic direction in the plasma frame. The mode velocity
in the plasma frame is faster in radial locations further away from
${\psi_{\mathrm{norm}}\approx 0.980}$ since the modes experience rigid body
rotation. For instance,
${v_{\mathrm{mode,pl.}}(\psi_{\mathrm{norm}}=0.934)=-9~{}\mathrm{km/s}}$ and
${-4.5~{}\mathrm{km/s}}$ at ${0.997}$, but where the mode amplitude is
maximised at ${0.980}$, as mentioned before, it is approximately
${+1~{}\mathrm{km/s}}$.

Figure 7: Poloidal velocity of the non-axisymmetric perturbations during the
linear phase (at ${t=1.60~{}\mathrm{ms}}$) and non-linear phase
(${t=2.70~{}\mathrm{ms}}$) together with the poloidal velocity of the plasma
along the ${\psi_{\mathrm{norm}}\approx 0.98}$ flux-surface.

Experimental measurements from Alcator C-mod have measured the mode velocity
in the laboratory frame to be in the electron diamagnetic direction always. In
ohmic heated-only EDA H-modes, it has been reported that the QCM moves in the
electron diamagnetic direction in the plasma frame (with measurements with
mirror Langmuir probes) [11] and in ICRF-heated EDA in the ion diamagnetic
direction in the plasma frame (with measurements from gas puff imaging) [10].
The latter showed that the mode amplitude is maximised in the maximum electron
temperature gradient and in that location ${v_{\mathrm{mode,pl.}}}$ is
measured to be in the ion diamagnetic direction, but at larger velocities than
observed in the simulations in the range ${\psi_{\mathrm{norm}}=\\{0.95\dots
0.99\\})}$, the modes move with ${\\{+15\dots+7\\}~{}\mathrm{km/s}}$ (in these
experiments, the mode velocity is compared with the ExB velocity, but the
poloidal projection of the parallel velocity is neglected). At this stage it
is worth mentioning that the simplified scrape-off layer model used in the
single-temperature JOREK simulation results in a radial electric field that is
unrealistic in the SOL. Future work may then involve a more advance SOL model
[37] and a separation between ion and electron temperatures to understand the
influence onto the mode dynamics, but this goes beyond the scope of the
present work.

From a ballooning mode dispersion relation, the velocity of ideal ballooning
modes has been reported in Ref. [38] to be
${v_{\mathrm{ExB}}+\bm{v}_{\parallel}\cdot\bm{b}_{\theta}+v^{*}_{i,\theta}/2}$,
where ${v_{\mathrm{ExB}}}$ lies solely in the poloidal direction (this is one
of the assumptions for the reduced MHD model in JOREK [34]), and for resistive
ballooning modes to be
${v_{\mathrm{ExB}}+\bm{v}_{\parallel}\cdot\bm{b}_{\theta}}$. The mode velocity
extracted from the simulations is then compared to the local
${v_{\theta}=\left(\bm{v}_{\mathrm{ExB}}+\bm{v}_{\parallel}\right)\cdot\bm{b}_{\theta}}$,
which varies across and along flux-surfaces. The mode velocity is found to
match best with the local ${v_{\theta}}$ where the mode amplitude is largest
(i.e., ${\psi_{\mathrm{norm}}\approx 0.980}$), as mentioned before. This is an
indication that the relevant modes are resistive peeling-ballooning modes.

The nature of the $n=7$-dominated perturbation along a toroidal line at
${(R,Z)=(2.130,0.042)~{}\mathrm{m}}$, which is roughly at the maximum
temperature gradient location (${\psi_{\mathrm{norm}}=0.98}$), is shown in
fig. 8. It shows variations of the density, temperature, and electrostatic
potential during the linear phase (top) and during the non-linear phase
(middle, and further explained in section 4). The cross-phase in the toroidal
angle, $\varphi$, for ${\langle n_{e},T_{e}\rangle}$ in dark-red and for
${\langle n_{e},\Phi\rangle}$ in dark-green are shown for the linear phase
with squares and for the non-linear phase with circles (bottom). During the
linear phase, the perturbations of density and temperature are ‘perfectly’ in-
phase with each-other, but have a cross-phase of ${\sim 10^{\circ}}$ with the
electrostatic potential.

Figure 8: Perturbations of density, temperature, and electrostatic potential
during the linear phase (top), during the non-linear phase (middle) along a
toroidal line with constant
$(R,Z_{\mathrm{axis}})=(2.130,0.042)~{}\mathrm{m}$. This location is at the
outer midplane roughly where the temperature gradient is maximised. The cross-
phases ${\langle n_{e},T_{e}\rangle}$ in dark-red and ${\langle
n_{e},\Phi\rangle}$ in dark-green (bottom) are shown with squares for the
linear phase and circles for the non-linear phase.

During the linear and early non-linear phase, where the non-axisymmetric
perturbations do not cause changes to the axisymmetric background, the
perturbations are well described with a fluctuation frequency of ${\sim
35~{}\mathrm{kHz}}$. The frequency of the quasi-coherent mode in the
experimental discharge can be observed in fig. 1 (d), and it is clear that the
fluctuation present in the simulation during the linear phase corresponds well
with the experimental $f_{\mathrm{QCM}}$. It is worth reiterating at this
stage that during the linear phase the $E_{r}$ minimum in the simulations
reaches ${-29~{}\mathrm{kV/m}}$, while in the experiments it is calculated to
be ${-20\pm 4~{}\mathrm{kv/m}}$.

##  4 Non-linear simulation results

The previous section described the linear phase in terms of growth-rates, mode
structure (dominant toroidal and poloidal mode numbers), location, and
velocity. The present section details the temporal dynamics at play when the
non-axisymmetric perturbations become large enough to interact with the
axisymmetric background plasma. The non-linear phase observes a shift of
dominant mode numbers from higher (${n=\\{6\dots 9\\}}$) to lower
(${n=\\{2\dots 5\\}}$). The mode activity causes perturbations to the
background magnetic field strong enough to create an ergodised layer in the
pedestal (outward of ${\psi_{\mathrm{norm}}\approx 0.965}$) and enhance
parallel heat transport, which results in a depletion of the temperature
pedestal. In addition, convective cells are formed and particle transport is
generated which also causes a reduction of the pedestal top density.

###  4.1 Shifting dominant toroidal mode numbers

The non-linear evolution of the non-axisymmetric perturbations (their magnetic
energies) is shown in fig. 9 from 2 to ${6~{}\mathrm{ms}}$ in linear (top) and
logarithmic (bottom) scales, and three different temporally-averaged spectra
(considering sample sizes of 1, 2, and ${3~{}\mathrm{ms}}$) centered around
${4~{}\mathrm{ms}}$ (bottom). From the first two figures it is clear that the
dominant mode numbers shift from higher-to-lower in time. There are two
reasons for this effect 1) non-linear mode coupling gives energy from the
higher-n perturbations that are linearly unstable (${n=\\{5\dots 11\\}}$, but
in particular ${n=6,}$ 7, and 8) to the lower-n that are linearly stable
(${n=\\{1\dots 4\\}}$) and 2) once the linearly most-unstable modes start to
have an impact on the axisymmetric background they cause the pedestal density
and temperature to become depleted, which in itself always has an effect on
the linear spectrum.

Figure 9: Linear (top) and logarithmic (middle) magnetic energies of the non-
axisymmetric perturbations during the non-linear phase (from 2 to
${6~{}\mathrm{ms}}$). Bottom: three different temporally-averaged spectra
centered around ${4~{}\mathrm{ms}}$ with sample sizes of 1, 2, and
${3~{}\mathrm{ms}}$.

Figure 9 (bottom) clearly shows that once the non-linear phase is underway the
dominant modes become ${n=2}$ and 3. The ${n=3}$-dominated structure can be
distinguished in fig. 8 (middle), where also several (11) local maxima can be
observed. The presence of this $n=11$ structure is a clear indication that
high toroidal mode numbers are involved in the non-linear dynamics as well.
During the non-linear phase, density, temperature, and electrostatic potential
are roughly in-phase (circles in the bottom part of fig. 8). The effect of the
non-axisymmetric perturbations onto the background plasma is detailed in the
following.

###  4.2 Interaction with the background plasma

As the non-linear phase gets underway, the interaction between background
plasma and non-axisymmetric perturbations causes heat and particle transport
in a continuous manner. This is in opposition to the transport caused by edge
localised modes in JOREK simulations which, when studied at realistic
resistivity and including diamagnetic effects and bootstrap current, feature
transient transport [39]. The rate of heat and particle losses resulting from
the non-axisymmetric perturbations correspond to ${\sim 2.5~{}\mathrm{MJ/s}}$
and ${\sim 6.9\times 10^{21}~{}\mathrm{ions/s}}$, respectively. The evolution
of edge density and temperature profiles (toroidally-averaged at the outboard
midplane) during the non-linear phase is shown in fig. 10. This causes a
steady depletion of the pedestal density, which may be partly responsible for
the reduction in the dominant toroidal mode numbers during the non-linear
phase. The temperature profile also decrease, but it is less affected because
the heating power remains constant and the number of particles in the pedestal
is constantly decreasing. The heat and particles that escape the magnetic
confinement due to the mode activity are ultimately deposited into the
simplified divertor targets.

Figure 10: Toroidally-averaged edge density and temperature profiles at the
outboard midplane during different times of the non-linear phase. The steady
depletion of the pedestal due to the edge non-axisymmetric perturbations can
be clearly observed (particularly in the density profiles).

###  4.3 Fluctuations across the pedestal

Locally, the temperature and density fluctuations in the pedestal are very
large, of around ${\delta T_{e}/\langle T_{e}\rangle\sim 50~{}\%}$ and
${\delta n_{e}/\langle n_{e}\rangle\sim 40~{}\%}$, which is consistent with
observations of the QCM for both C-mod [11] and AUG [15]. This, together with
a spectrogram of the temperature fluctuations signal is shown in fig. 11. In
the spectrogram it is possible to determine that in the early non-linear phase
(up to ${\sim 2.2~{}\mathrm{ms}}$) the fluctuations cover the frequency range
${\\{30\dots 45\\}~{}\mathrm{kHz}}$ and their frequency chirps down when the
non-linear phase gets underway. This is due to the fact that the radial
electric field well becomes less deep due to the interaction between non-
axisymmetric perturbations and background plasma. During the non-linear phase,
the frequency remains relatively constant in a frequency range of ${\\{8\dots
18\\}~{}\mathrm{kHz}}$. On average, during this phase, the $E_{r}$ minimum is
${\sim-15~{}\mathrm{kV/m}}$, which is shallower than the experimental
measurement (${-20\pm 4~{}\mathrm{kV/m}}$). From the experimental measurements
through helium beam (fig. 1), it is clear that the QCM frequency sits in the
frequency range of ${f_{\mathrm{QCM}}\approx\\{20\dots 40\\}~{}\mathrm{kHz}}$.
This frequency range overlaps with what is obtained in the simulations during
the linear phase (where the profiles have not yet experienced the depletion
due to non-linear interactions with the perturbations). This is an indication
that the resistive peeling-ballooning modes in the simulation could be related
to the QCM.

Figure 11: Temperature (top) and density (middle) fluctuations relative to the
toroidally-averaged value in the location where the mode amplitude is
maximised (${\psi_{\mathrm{norm}}\approx 0.98}$). Bottom: spectrogram of the
temperature fluctuations. At the end of the early non-linear phase
(${t\lesssim 2.2~{}\mathrm{ms}}$) the pedestal fluctuations correspond to a
frequency of roughly $35~{}\mathrm{kHz}$, but they chirp down to around
${13~{}\mathrm{kHz}}$ in the non-linear phase.

The poloidal structure of the non-axisymmetric perturbations (along a single
flux-surface) features a predominantly ballooning structure. In the outer
midplane during the non-linear phase, the range of poloidal wavenumbers is
${k_{\theta}\sim\\{0.1\dots 0.5\\}~{}\mathrm{cm^{-1}}}$. Conversely, the
reported poloidal wavenumbers for the QCM in C-mod are ${\sim
1.5~{}\mathrm{cm^{-1}}}$. For AUG, based on measurement with helium beam the
range of poloidal wavenumers for the QCM is ${\sim 0.6}$ to
${0.9~{}\mathrm{cm^{-1}}}$, and for fluctuations during the quasi-continuous
exhaust regime ${k_{\theta}\sim\\{0.4\dots 0.9\\}~{}\mathrm{cm^{-1}}}$. As
such, the poloidal wavenumbers hold a closer resemblance to the fluctuations
measured during the QCE regime than the quasi-coherent mode during the EDA.

##  5 Discussion and conclusions

In this paper, first non-linear extended MHD simulations that consider initial
conditions from an experimental discharge in ASDEX Upgrade that corresponds to
the enhanced D-alpha H-mode (EDA) were presented. The equilibrium
reconstruction was carried out with the CLISTE code [26] and linear ideal MHD
stability analysis with MISHKA [27], which found the experimental point inside
the stable region near the high-n ballooning boundary. The experimental
discharge under consideration features the EDA with a quasi-coherent mode
(QCM) in a frequency range of 20 to ${40~{}\mathrm{kHz}}$ and is completely
devoid of large edge localised modes. The cross-field transport that prevents
the pedestal from building up is thus thought to be caused by the QCM. In
Alcator C-mod, where the EDA was first observed, and routinely obtained,
magnetic fluctuations were not observed in magnetic pick-up coils, but in
ASDEX Upgrade, magnetic pick-up coils in the outboard midplane indeed observe
fluctuations with toroidal mode numbers around ${n=\\{3\dots 7\\}}$. The
fundamental (i.e., $n=1$) frequency of these fluctuations in the discharge
under consideration is ${32.9~{}\mathrm{kHz}}$. The relationship between these
magnetic fluctuations and the QCM is not yet clear; the fundamental frequency
in some discharges matches that of the QCM, while in others it displays a
difference to ${f_{\mathrm{QCM}}}$ by up to $\sim 15~{}\mathrm{kHz}$.

The simulations presented in this paper are carried out with the JOREK code
[24, 25]. The linearly most unstable modes are found to be resistive peeling-
ballooning modes with toroidal mode numbers in the range ${n=\\{6\dots 9\\}}$.
The identification of such instabilities as resistive peeling-ballooning modes
comes from their spatial structure (predominantly on the low-field side) and
from their poloidal rotation velocity, which corresponds to that of resistive
ballooning modes: ${v_{\mathrm{ExB}}+\bm{v}_{\parallel}\cdot\bm{b}_{\theta}}$
[38]. The location of the maximum mode amplitude corresponds to the maximum of
the electron temperature gradient, and at this location the poloidal mode
velocity is in the ion diamagnetic direction in the plasma frame; both such
observations are in qualitative agreement with experimental measurements from
gas puff imaging in Alcator C-mod [10]. As the linearly most unstable modes
(notably $n=7$) grow in amplitude, they are able to non-linearly drive the
linearly stable modes through three-wave interactions. As such, perturbations
with lower and higher toroidal mode numbers (than ${n=4}$ and $9$,
respectively) become excited. This early non-linear phase gives way to the
fully non-linear phase when the non-axisymmetric perturbations grow large
enough to interact with the background plasma. As the non-linear phase gets
underway, the dominant toroidal mode number decreases from ${n=7}$ to ${n=3}$
but intermediate and high toroidal mode numbers remain relevant. During this
phase, particles and heat are quasi-continuously expelled by the cross-field
transport generated by the non-axisymmetric perturbations. As such, the
pedestal density starts to become depleted.

The density and temperature fluctuations during the non-linear phase have
large amplitudes, which is qualitatively consistent with measurements from EDA
in ASDEX Upgrade [15]. Before the early non-linear phase comes to an end,
these fluctuations lie on a frequency range of ${\\{30\dots
40\\}~{}\mathrm{kHz}}$, which is consistent with the observations for the
present discharge (before the density and temperature profiles are affected by
the non-linear interaction between the non-axisymmetric perturbations and the
background plasma). During the fully non-linear phase, when the density
pedestal is depleting and the radial electric field well is becoming
shallower, the fluctuation frequencies decrease to ${\\{8\dots
18\\}~{}\mathrm{kHz}}$ (which still sits at the lower end of the observed
values of ${f_{\mathrm{QCM}}}$ in AUG [15]).

In summary, from the analysis performed of the JOREK simulations of an EDA
H-mode, several quantitative and qualitative similarities are found with
respect to experimental features and measurements of the quasi-coherent mode
and magnetic fluctuations captured by the pick-up coils. Further comparisons
between non-linear extended MHD simulations and experimental measurements from
the EDA H-mode are left for future work. In particular, a focus on more
quantitative comparisons of the EDA H-mode itself and of its boundaries,
including its relation to the quasi-continuous exhaust (QCE) regime is
foreseen. To this end, taking into account separate electron and ion
temperatures, resistive walls instead of the ideal wall boundary condition
employed here, advanced SOL models, toroidal rotation sources, and potentially
higher order finite Larmor radius corrections are under consideration.
Additional dependencies on the viscosity values, heat and particle sources,
and impurity radiation can additionally be taken into account.

## Acknowledgements

This work has been carried out within the framework of the EUROfusion
Consortium, funded by the European Union via the Euratom Research and Training
Programme (Grant Agreement No 101052200 — EUROfusion). Views and opinions
expressed are however those of the author(s) only and do not necessarily
reflect those of the European Union or the European Commission. Neither the
European Union nor the European Commission can be held responsible for them.
In particular, A. Cathey received funding from a EUROfusion Researcher Grant
and work package Tokamak Exploitation, and are acknowledged. The non-linear
JOREK simulations were performed using the Marconi-Fusion supercomputer within
the FUA36_MHD project.

## References

## References

  * [1] Hartmut Zohm.  Edge localized modes (elms).  Plasma Physics and Controlled Fusion, 38(2):105, 1996.
  * [2] Anthony W Leonard.  Edge-localized-modes in tokamaks.  Physics of Plasmas, 21(9):090501, 2014.
  * [3] T Eich, B Sieglin, AJ Thornton, M Faitsch, A Kirk, A Herrmann, W Suttrop, et al.  Elm divertor peak energy fluence scaling to iter with data from jet, mast and asdex upgrade.  Nuclear Materials and Energy, 12:84–90, 2017.
  * [4] J.P. Gunn, S. Carpentier-Chouchana, F. Escourbiac, T. Hirai, S. Panayotis, R.A. Pitts, Y. Corre, R. Dejarnac, M. Firdaouss, M. Kočan, M. Komm, A. Kukushkin, P. Languille, M. Missirlian, W. Zhao, and G. Zhong.  Surface heat loads on the ITER divertor vertical targets.  Nuclear Fusion, 57(4):046025, mar 2017.
  * [5] N Oyama, P Gohil, L D Horton, A E Hubbard, J W Hughes, Y Kamada, K Kamiya, A W Leonard, A Loarte, R Maingi, G Saibene, R Sartori, J K Stober, W Suttrop, H Urano, W P West, and the ITPA Pedestal Topical Group.  Pedestal conditions for small ELM regimes in tokamaks.  Plasma Physics and Controlled Fusion, 48(5A):A171–A181, apr 2006\.
  * [6] E. Viezzer.  Access and sustainment of naturally ELM-free and small-ELM regimes.  Nuclear Fusion, 58(11):115002, sep 2018.
  * [7] E. Viezzer, M.E. Austin, M. Bernert, K.H. Burrell, P. Cano-Megias, X. Chen, D.J. Cruz-Zabala, S. Coda, M. Faitsch, O. Fevrier, L. Gil, C. Giroud, T. Happel, G.F. Harrer, A.E. Hubbard, J.W. Hughes, A. Kallenbach, B. Labit, A. Merle, H. Meyer, C. Paz-Soldan, P. Oyola, O. Sauter, M. Siccinio, D. Silvagni, and E.R. Solano.  Prospects of core–edge integrated no-elm and small-elm scenarios for future fusion devices.  Nuclear Materials and Energy, page 101308, 2022.
  * [8] Y. Takase, R. L. Boivin, F. Bombarda, P. T. Bonoli, C. Christensen, C. Fiore, D. Garnier, J. A. Goetz, S. N. Golovato, R. Granetz, M. Greenwald, S. F. Horne, A. Hubbard, I. H. Hutchinson, J. Irby, B. LaBombard, B. Lipschultz, E. Marmar, M. May, A. Mazurenko, G. McCracken, P. O’Shea, M. Porkolab, J. Reardon, J. Rice, C. Rost, J. Schachter, J. A. Snipes, P. Stek, J. Terry, R. Watterson, B. Welch, and S. Wolfe.  Radiofrequency-heated enhanced confinement modes in the alcator c-mod tokamak.  Physics of Plasmas, 4(5):1647–1653, 1997.
  * [9] M. Greenwald, R. Boivin, P. Bonoli, R. Budny, C. Fiore, J. Goetz, R. Granetz, A. Hubbard, I. Hutchinson, J. Irby, B. LaBombard, Y. Lin, B. Lipschultz, E. Marmar, A. Mazurenko, D. Mossessian, T. Sunn Pedersen, C. S. Pitcher, M. Porkolab, J. Rice, W. Rowan, J. Snipes, G. Schilling, Y. Takase, J. Terry, S. Wolfe, J. Weaver, B. Welch, and S. Wukitch.  Characterization of enhanced d-alpha high-confinement modes in alcator c-mod.  Physics of Plasmas, 6(5):1943–1949, 1999.
  * [10] C Theiler, J L Terry, E Edlund, I Cziegler, R M Churchill, J W Hughes, B LaBombard, T Golfinopoulos, and the Alcator C-Mod Team.  Radial localization of edge modes in alcator c-mod pedestals using optical diagnostics.  Plasma Physics and Controlled Fusion, 59(2):025016, jan 2017.
  * [11] B. LaBombard, T. Golfinopoulos, J. L. Terry, D. Brunner, E. Davis, M. Greenwald, and J. W. Hughes.  New insights on boundary plasma turbulence and the quasi-coherent mode in alcator c-mod using a mirror langmuir probe.  Physics of Plasmas, 21(5):056108, 2014.
  * [12] A. E. Hubbard, R. L. Boivin, R. S. Granetz, M. Greenwald, J. W. Hughes, I. H. Hutchinson, J. Irby, B. LaBombard, Y. Lin, E. S. Marmar, A. Mazurenko, D. Mossessian, E. Nelson-Melby, M. Porkolab, J. A. Snipes, J. Terry, S. Wolfe, S. Wukitch, B. A. Carreras, V. Klein, and T. Sunn Pedersen.  Pedestal profiles and fluctuations in c-mod enhanced d-alpha h-modes.  Physics of Plasmas, 8(5):2033–2040, 2001.
  * [13] J A Snipes, B LaBombard, M Greenwald, I H Hutchinson, J Irby, Y Lin, A Mazurenko, and M Porkolab.  The quasi-coherent signature of enhanced d-alpha h-mode in alcator c-mod.  Plasma Physics and Controlled Fusion, 43(4):L23, apr 2001.
  * [14] P. J. Sun, Y. D. Li, Y. Ren, X. D. Zhang, G. J. Wu, Y. M. Wang, T. H. Shi, B. Lyu, Y. M. Duan, J. Z. Zhang, F. D. Wang, H. Q. Liu, S. T. Mao, L. Q. Xu, Q. Li, P. Li, J. Bi, L. Q. Hu, and J. G. Li.  Experimental study of quasi-coherent mode using east tangential co2 laser collective scattering diagnostic in far-forward mode.  Physics of Plasmas, 26(1):012304, 2019.
  * [15] L. Gil, C. Silva, T. Happel, G. Birkenmeier, G.D. Conway, L. Guimarais, A. Kallenbach, T. Pütterich, J. Santos, P.A. Schneider, M. Schubert, E. Seliunin, A. Silva, J. Stober, U. Stroth, E. Trier, E. Wolfrum, and and.  Stationary ELM-free h-mode in ASDEX upgrade.  Nuclear Fusion, 60(5):054003, apr 2020.
  * [16] A. Kallenbach, M. Bernert, P. David, M. G. Dunne, R. Dux, E. Fable, R. Fischer, L. Gil, T. Görler, F. Janky, R. M. McDermott, W. Suttrop, G. Tardini, ASDEX Upgrade Team M. Wischmeier, and EUROfusion MST1 team.  Developments towards an elm-free pedestal radiative cooling scenario using noble gas seeding in asdex upgrade.  Nuclear Fusion, 61(1):016002, nov 2020.
  * [17] C Paz-Soldan and the DIII-D Team.  Plasma performance and operational space without elms in diii-d.  Plasma Physics and Controlled Fusion, 63(8):083001, jun 2021.
  * [18] J.L. Terry, N.P. Basse, I. Cziegler, M. Greenwald, O. Grulke, B. LaBombard, S.J. Zweben, E.M. Edlund, J.W. Hughes, L. Lin, Y. Lin, M. Porkolab, M. Sampsell, B. Veto, and S.J. Wukitch.  Transport phenomena in the edge of alcator c-mod plasmas.  Nuclear Fusion, 45(11):1321, oct 2005.
  * [19] D A Mossessian, P B Snyder, M Greenwald, J W Hughes, Y Lin, A Mazurenko, S Medvedev, H R Wilson, and S Wolfe.  H-mode pedestal characteristics and mhd stability of the edge plasma in alcator c-mod.  Plasma Physics and Controlled Fusion, 44(4):423, mar 2002.
  * [20] A. Mazurenko, M. Porkolab, D. Mossessian, J. A. Snipes, X. Q. Xu, and W. M. Nevins.  Experimental and theoretical study of quasicoherent fluctuations in enhanced ${D}_{\alpha}$ plasmas in the alcator c-mod tokamak.  Phys. Rev. Lett., 89:225004, Nov 2002.
  * [21] J. R. Myra and D. A. D’Ippolito.  Edge instability regimes with applications to blob transport and the quasicoherent mode.  Physics of Plasmas, 12(9):092511, 2005.
  * [22] M. Porkolab, J.C. Rost, N. Basse, J. Dorris, E. Edlund, Liang Lin, Y. Lin, and S. Wukitch.  Phase contrast imaging of waves and instabilities in high temperature magnetized fusion plasmas.  IEEE Transactions on Plasma Science, 34(2):229–234, 2006.
  * [23] Karl Stimmel, L. Gil, T. Görler, M. Cavedon, P. David, M. Dunne, R. Dux, R. Fischer, F. Jenko, A. Kallenbach, and et al.  Gyrokinetic analysis of an argon-seeded eda h-mode in asdex upgrade.  Journal of Plasma Physics, 88(3):905880315, 2022.
  * [24] M. Hoelzl, G.T.A. Huijsmans, S.J.P. Pamela, M. Bécoulet, E. Nardon, F.J. Artola, B. Nkonga, C.V. Atanasiu, V. Bandaru, A. Bhole, D. Bonfiglio, A. Cathey, O. Czarny, A. Dvornova, T. Fehér, A. Fil, E. Franck, S. Futatani, M. Gruca, H. Guillard, J.W. Haverkort, I. Holod, D. Hu, S.K. Kim, S.Q. Korving, L. Kos, I. Krebs, L. Kripner, G. Latu, F. Liu, P. Merkel, D. Meshcheriakov, V. Mitterauer, S. Mochalskyy, J.A. Morales, R. Nies, N. Nikulsin, F. Orain, J. Pratt, R. Ramasamy, P. Ramet, C. Reux, K. Särkimäki, N. Schwarz, P. Singh Verma, S.F. Smith, C. Sommariva, E. Strumberger, D.C. van Vugt, M. Verbeek, E. Westerhof, F. Wieschollek, and J. Zielinski.  The JOREK non-linear extended MHD code and applications to large-scale instabilities and their control in magnetically confined fusion plasmas.  Nuclear Fusion, 61(6):065001, may 2021.
  * [25] GTA Huysmans and O Czarny.  Mhd stability in x-point geometry: simulation of elms.  Nuclear fusion, 47(7):659, 2007.
  * [26] PJ Mc Carthy.  Analytical solutions to the grad–shafranov equation for tokamak equilibrium with dissimilar source functions.  Physics of Plasmas, 6(9):3554–3560, 1999.
  * [27] A.B Mikhailovskii, G.T.A. Huysmans, W.O.K. Kerner, and S.E. Sharapov.  Optimization of computational mhd normal-mode analysis for tokamaks.  ”Plasma Physics Reports”, 23(10):844–857, 1997.
  * [28] A Cathey, M Hoelzl, G Harrer, M G Dunne, G T A Huijsmans, K Lackner, S J P Pamela, E Wolfrum, S Günter, the JOREK Team, the ASDEX Upgrade Team, and the EUROfusion MST1 Team.  Mhd simulations of small elms at low triangularity in asdex upgrade.  Plasma Physics and Controlled Fusion, 64(5):054011, apr 2022.
  * [29] A. Kleiner, N.M. Ferraro, A. Diallo, and G.P. Canal.  Importance of resistivity on edge-localized mode onset in spherical tokamaks.  Nuclear Fusion, 61(6):064002, apr 2021.
  * [30] H. Nyström, L. Frassinetti, S. Saarelma, G.T.A. Huijsmans, C. Perez von Thun, C.F. Maggi, J.C. Hillesheim, and JET contributors.  Effect of resistivity on the pedestal mhd stability in jet.  Nuclear Fusion, 62(12):126045, oct 2022.
  * [31] S. J. P. Pamela, A. Bhole, G. T. A. Huijsmans, B. Nkonga, M. Hoelzl, I. Krebs, and E. Strumberger.  Extended full-mhd simulation of non-linear instabilities in tokamak plasmas.  Physics of Plasmas, 27(10):102510, 2020.
  * [32] Olivier Sauter, Clemente Angioni, and YR Lin-Liu.  Neoclassical conductivity and bootstrap current formulas for general axisymmetric equilibria and arbitrary collisionality regime.  Physics of Plasmas, 6(7):2834–2839, 1999.
  * [33] Olivier Sauter and Clemente Angioni.  Erratum: “neoclassical conductivity and bootstrap current formulas for general axisymmetric equilibria and arbitrary collisionality regime” [phys. plasmas 6, 2834 (1999)].  Physics of Plasmas, 9(12):5140–2839, 2002.
  * [34] A Cathey, M Hoelzl, S Futatani, P T Lang, K Lackner, G T A Huijsmans, S J P Pamela, S Günter, the JOREK Team, the ASDEX Upgrade Team, and the EUROfusion MST1 Team.  Comparing spontaneous and pellet-triggered elms via non-linear extended mhd simulations.  Plasma Physics and Controlled Fusion, 63(7):075016, jun 2021.
  * [35] Isabel Krebs, Matthias Hoelzl, Karl Lackner, and Sibylle Günter.  Nonlinear excitation of low-n harmonics in reduced magnetohydrodynamic simulations of edge-localized modes.  Physics of plasmas, 20(8):082506, 2013.
  * [36] J Kalis.  personal communication.
  * [37] S.Q. Korving, G.T.A. Huijsmans, J.-S. Park, A. Loarte, and the JOREK team.  Development of the neutral model in the nonlinear mhd code jorek: Application to exb drifts in iter pfpo-1 plasmas.  Physics of Plasmas (submitted), 2023.
  * [38] JA Morales, Marina Becoulet, X Garbet, F Orain, G Dif-Pradalier, M Hoelzl, S Pamela, GTA Huijsmans, P Cahyna, A Fil, et al.  Edge localized mode rotation and the nonlinear dynamics of filaments.  Physics of Plasmas, 23(4):042513, 2016.
  * [39] A. Cathey, M. Hoelzl, K. Lackner, G.T.A. Huijsmans, M.G. Dunne, E. Wolfrum, S.J.P. Pamela, F. Orain, S. Günter, the JOREK team, the ASDEX Upgrade team, and the EUROfusion MST1 team.  Non-linear extended mhd simulations of type-i edge localised mode cycles in asdex upgrade and their underlying triggering mechanism.  Nuclear Fusion, 60(12):124007, nov 2020.