# Analysis of energetic particle driven modes and interchange modes in LHD plasma using the gyro-fluid code FAR3d

J. Varela, K. Nagaoka, S. Ohdachi, K. Watanabe, Y. Takemura, Y.

Narushima, T. Tokuzawa, K. Tanaka, K. Ida, Y. Todo, et al.

► To cite this version:

J. Varela, K. Nagaoka, S. Ohdachi, K. Watanabe, Y. Takemura, et al., Analysis of energetic particle driven modes and interchange modes in LHD plasma using the gyro-fluid code FAR3d. Journal of Fusion Energy, 2025, 10.21203/rs.3.rs-6718076/v1. hal-05468262

HAL Id: hal-05468262

https://hal.science/hal-05468262v1

Submitted on 20 Jan 2026

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
Characterization, analysis and optimization of
energetic particle driven modes and interchange
modes in LHD plasma using the gyro-fluid code
FAR3d

J. Varela

E-mail: jacobo.rodriguez@austin.utexas.edu
Institute for Fusion Studies, Department of Physics, University of Texas at
Austin, Austin, Texas 78712, USA

K. Nagaoka

National Institute for Fusion Science, National Institute of Natural Science,
Toki, 509-5292, Japan

S. Ohdachi

National Institute for Fusion Science, National Institute of Natural Science,
Toki, 509-5292, Japan

K. Y. Watanabe

National Institute for Fusion Science, National Institute of Natural Science,
Toki, 509-5292, Japan

Y. Takemura

National Institute for Fusion Science, National Institute of Natural Science,
Toki, 509-5292, Japan

Y. Narushima

National Institute for Fusion Science, National Institute of Natural Science,
Toki, 509-5292, Japan

T. Tokuzawa

National Institute for Fusion Science, National Institute of Natural Science,
Toki, 509-5292, Japan

K. Tanaka

National Institute for Fusion Science, National Institute of Natural Science,
Toki, 509-5292, Japan

K. Ida
National Institute for Fusion Science, National Institute of Natural Science,
Toki, 509-5292, Japan

Y. Todo

National Institute for Fusion Science, National Institute of Natural Science,
Toki, 509-5292, Japan

M. Yoshinuma

National Institute for Fusion Science, National Institute of Natural Science,
Toki, 509-5292, Japan

R. Seki

SOKENDAI, Department of Fusion Science, Toki/Gifu and National Institute
for Fusion Science, National Institute of Natural Science, Toki, 509-5292, Japan

K. Nagasaki

Institute of Advanced Energy, Kyoto University, Uji, Japan

S. Kobayashi

Institute of Advanced Energy, Kyoto University, Uji, Kyoto, Japan

P. Adulsiriswad

National Institute for Fusion Science, National Institute of Natural Science,
Toki, 509-5292, Japan

S. Yamamoto

National Institutes for Quantum and Radiological Science and Technology,
Naka, Ibaraki 311-0193, Japan

D. A. Spong

Oak Ridge National Laboratory, Oak Ridge, Tennessee 37831-8071, USA

L. Garcia

Universidad Carlos III de Madrid, 28911 Leganes, Madrid, Spain

A. Cappa

Laboratorio Nacional de Fusion CIEMAT, Madrid, Spain

Y. Ghai

Oak Ridge National Laboratory, Oak Ridge, Tennessee 37831-8071, USA

J. Ortiz

Universidad Carlos III de Madrid, 28911 Leganes, Madrid, Spain
W. A. Cooper

Ecole Polytechnique de Lausanne (EPFL), Centre de Recherches en Physique des Plasma (CRPP), CH-1015 Lausanne, Switzerland

X. Du

General Atomics, P.O. Box 85608, San Diego, California 92186-5608, USA

S. Sharapov

UKAEA, Culham Science Centre, Abingdon, Oxfordshire OX14 3DB, UK

F. L. Waelbroeck

Institute for Fusion Studies, Department of Physics, University of Texas at Austin, Austin, Texas 78712, USA

B. Breizman

Institute for Fusion Studies, Department of Physics, University of Texas at Austin, Austin, Texas 78712, USA

D. Zarzoso

Aix-Marseille Université, CNRS, Centrale Marseille, M2P2 UMR 7340, Marseille, France

C. Hidalgo

Laboratorio Nacional de Fusion—CIEMAT, Madrid, Spain

A. Azegami

National Institute for Fusion Science, National Institute of Natural Science, Toki, 509-5292, Japan

**Abstract.** Energetic particles (EP) generated by the neutral beam injectors (NBI) in Large Helical Device (LHD) destabilize Alfvén Eigenmodes (AE) and energetic particle modes (EPM), leading to a reduction of the device performance. In particular, AE/EPM induce EP losses before thermalization that causes a lower plasma heating efficiency and damages to plasma facing components. Pressure gradient driven modes (PGDM) as interchange and ballooning modes also hamper the capability of LHD to confine the thermal plasma, limiting the maximum thermal plasma $\beta$ of the discharge. The present study summarizes AE/EPM and PGDM characterization and optimization analysis performed in LHD plasma using the gyro-fluid code FAR3d. The linear and saturation phase of AE/EPM and PGDM leading to the destabilization of MHD burst, energetic-ion-driven resistive interchange mode (EIC) burst, internal collapse and saw-tooth like events are discussed. Optimization trends with respect to the NBI operation regime, magnetic field configuration, thermal plasma properties, external actuators and multiple EP populations effects are explored, comparing experiment and simulation data. Improved operation scenarios are identified and confirmed in dedicated experiments thanks to the minimization or avoidance of AE/EPM and PGDM along the discharge.
PACS numbers: 52.35.Py, 52.55.Hc, 52.55.Tn, 52.65.Kj

Keywords: Stellarator, LHD, EIC, MHD, AE, energetic particles, interchange modes
# 1. Introduction

The performance of fusion devices can be improved by avoiding or mitigating the instabilities triggered in the plasma. For the case of Stellarators, and particularly in the Large Helical Device (LHD), the pressure gradient driven modes (PGDM) hamper the device performance by a degradation of the device confinement leading to larger thermal plasma fluxes, limiting the maximum thermal $\beta$ of the operation scenarios [1, 2, 3, 4, 5]. A typical example is the resistive interchange modes (RIC) destabilized in the magnetic hill region linked to the rational surface 1/1 resonating at the plasma periphery [6, 7, 8]. In addition, energetic particle (EP) driven instabilities reduce the LHD plasma heating efficiency [9, 10, 11, 12, 13, 14, 15, 16, 17] causing a lower device performance [18, 19, 20, 21, 22, 23, 24, 25, 26].

EP driven instabilities increase the turbulent transport of EPs generated by heating mechanism as neutral beams and ion cyclotron resonance heating (ICRH) in LHD plasma [27, 28, 29, 30], a phenomena observed in multiple Stellarator devices as W7-AS, TJ-II and Heliotron J as well as Tokamak such as JET and DIII-D [31, 32, 33, 34, 35, 36, 37, 38, 39]. The EP losses are observed if there is a resonance between the frequency of an unstable mode and the EP drift, bounce or transit frequency.

Alfvén Eigenmodes (AE) are driven in the spectral gaps of the shear Alfvén continua [40, 41]. Each AE family is linked to periodic variations of the Alfvén speed, for example the toroidicity induced Alfvén Eigenmodes (TAE) couple $m$ with $m+1$ modes [42, 43, 44] ($m$ is the poloidal mode number), ellipticity induced AE (EAE) coupling $m$ with $m+2$ modes [45, 46] and noncircularity induced AE (NAE) coupling $m$ with $m+3$ or higher [47, 48]. In addition, beta induced AEs (BAE) are driven by compressibility effects [49], Reversed-shear AE (RSAE) due to local minimum in the safety factor $q$ profile [50], Global AE (GAE) are observed nearby the minimum of the Alfvén continua [51, 52] and helical AE (HAE) if different toroidal modes and helicities ($m/n$) are coupled [11] (with $n$ the toroidal mode family). Energetic particle modes (EPM) can be triggered inside the shear Alfvén continua if the continuum damping is not strong enough and/or the drive can keep them above marginal stability thresholds [53, 54, 55, 56, 57].

The present review is dedicated to discuss the

main contributions of the FAR3d gyro-fluid code in the analysis of AE / EPM and interchange modes in LHD plasma. Several studies were dedicated to characterize the linear and saturation phase of instabilities that trigger dangerous plasma relaxations linked to an important reduction of the device performance, specially the MHD burst [58], EIC burst [59], internal collapse [60, 61, 62, 63] and saw-tooth like events [64, 65].

LHD plasma is heated by three tangential NBI lines almost parallel to the magnetic axis injected at the plasma core with a beam energy of 180 keV, as well as two NBIs perpendicular to the magnetic axis injected at the plasma periphery with a beam energy of 32 keV [66, 67, 68]. There are multiple observations of AEs destabilized by the passing EPs injected by the tangential NBI [69] and by the helically trapped EP generated by the perpendicular NBI [70]. LHD can be operated modifying the location of the vacuum magnetic axis ($R_{ax}$), namely inward shifted configurations if the magnetic axis is located at $R_{ax} < 3.65$ m and outward shifted configurations if $R_{ax} > 3.75$ m [71]. Inward shifted configurations shows adverse MHD stability because the magnetic well only exists in the plasma core for high $\beta$ operations (due to the effect of the Shafranov shift) still showing a broad magnetic hill at the plasma periphery [72]. Nevertheless, stable plasma in inward shifted configurations are observed for thermal plasma $\beta > 3\%$ [73, 74, 75]. On the other hand, dangerous plasma relaxations reducing the device performance as internal collapses and saw-tooth like oscillations are observed. Outward shifted configurations show stable interchange mode because a magnetic well covers main part of the plasma, although fast growing ballooning modes are destabilized at the plasma periphery for regions with local bad curvature [76], leading to the onset of core density collapses [77, 78, 79]. It should be noted that not only the plasma stability is affected by the vacuum magnetic axis configuration, inward shifted configurations show a favorable confinement of the thermal plasma while EP transport in weakly affected [80, 81].

The numerical analysis discussed in this review are performed using the gyro-fluid code FAR3d that solves the reduced nonlinear resistive MHD equations coupled with the moment of the gyro-kinetic equation, particularly the EP density and parallel velocity [82, 83, 84, 85, 86, 87, 88]. The Landau damping/growth
is included by Landau closure relations, reproducing the linear wave-particle resonance effects on six field variables that evolve from a three dimensional equilibria generated by the VMEC code [89]. The main advantage of FAR3d is the computational efficiency, the code is designed to minimize the kinetic effects included in the model while retaining the essential elements to reproduce the experimental observations. FAR3d is applied to perform parametric and profiles scans required for optimization and design studies, the identification of AE stability trends as well as complementing more complex numerical models by providing a first estimate of AE / EPM linear and non-linear stability narrowing down the parameter selection.

This paper is organized as follows. There is a
basic introduction of the numerical model in section ??.
Studies dedicated to identify and analyze instabilities
are discussed in section 3. Optimization analysis are
introduced in section 4. Multiple EP effects studies
are presented in section 5. The saturation phase of
AE/EPM and PGDM is explored in section 6. Next,
preliminary results of ongoing studies are introduced
in section 7. Finally, the conclusions are presented in
section 8.

## 2. Numerical model

FAR3d code solves a set of reduced MHD equations
retaining the toroidal angle dependency in a three-
dimensional VMEC equilibrium [89, 90]. The EP
perturbation is included by moments of the gyro-
kinetic equation: the EP density (n_f) and the EP
velocity parallel to the magnetic field lines (v_||_f).
The truncation of the gyro-kinetic equation moments
hierarchy is done using the Landau closure coefficients,
obtained from the answer function of gyro-kinetic
simulations matching the analytic TAE growth rates of
the two-pole approximation of the plasma dispersion
function, leading to a Lorentzian energy distribution
function for the EP. The two-moment gyrofluid model
used here is based on a Lorentzian distribution function
that can be matched to a Maxwellian or to a slowing-
down distribution by choosing an equivalent average
energy. The averaged Maxwellian EP distribution
function used in the simulations has the same second
moment, the effective EP temperature, as that of
the equivalent slowing down distribution. It must
be mentioned that a single EP Maxwellian distribution
cannot reproduce the same resonance as a slowing
down distribution, because the gradient of the phase
space distribution determines the drive of the AE
modes, although parametric studies with respect to the
EP energy and β performed using a set of Maxwellian
can be used to approximate the resonances triggered

by a slowing down distribution function. More details
of the model equations and numerical scheme are found
in the reference [91].

Previous benchmarking studies validated FAR3d code results with respect to gyro-kinetic and hybrid codes [92]. FAR3d code was applied to model the MHD activity in DIII-D [93, 94, 95, 96, 97], TJ-II [98, 99, 100], Heliotron J [101, 102, 103] and EAST [104] plasma, showing a reasonable agreement with the experimental observations. Also, the nonlinear version of the code reproduced the AE saturation in DIII-D plasma [105, 106]. In addition, the FAR3d code was applied to analyze the AE stability in plasma with multiple EP populations for DIII-D, ITER and CFETR [107, 108, 109] and to predict the MHD stability of devices as JT60SA, CFQS and QPS [110, 111, 112].

FAR3d can provide timely support/interpretation
during and in between experimental campaigns
because linear simulations can be performed in seconds
or minutes depending on the number of modes
included in the model. Thus, the favorable code
performance helps to reduce the time gap between
theoretical analysis and experimental results, linked to
the large computational demand of more sophisticated
models and the accuracy of the input data (profiles,
equilibrium, realistic EP distributions, etc). Reduced
models as FAR3d can help to reduce this time
lag by providing a first order approach for the
characterization of the AE / EPM activity measured
in the experiments as well as information of the AE
/ EPM stability trends with respect to the main
operation parameters of the device.

The first step of FAR3d analysis consist in
the identification of the EP resonance that causes
the destabilization of the AE / EPM measured in
the experiment. On that aim, the code output
is compared with the data of several diagnostics,
considering successful the instability identification if
the dominant modes, frequency range, radial location
and eigenfunction structure obtained in the linear
simulations is consistent with the experimental data.
On that aim, the simulation output is compared with
the magnetic fluctuations measured by Mirnov coils
or Langmuir probes (B̃), electron density (ñe) and
temperature fluctuations (T̃e) from electron cyclotron
emission (ECE) plus Thomson scattering diagnostics,
heavy ion beam probe diagnostic measurements of the
plasma electrostatic potential (Φ) and tomographic
reconstructions of the instability using soft X-ray
diagnostics. Examples are provided in section
3. Next, if the linear simulations can provide a
reasonable description of the experiment, a set of
parametric studied are performed with respect to
simulation variables directly connected to experimental
parameters, identifying potential optimized scenarios with reduced AE / EPM activity. Theoretical optimized scenarios are confirmed by performing dedicated experiments. These kind of analysis are introduced in section 4. The last step of the study is analyzing the saturation phase of the AE / EPM, performing nonlinear simulations that provides information of the turbulent EP transport, the energy exchange between different toroidal / helical families and towards the thermal plasma as well as the generation of zonal structures. Examples are shown in the section 6.

3. Mode identification

The first task of the gyro-fluid version of FAR3d code
was the identification and characterization of the AE
/ EPM activity measured in LHD discharges. In
particular, the analysis were dedicated to study the
stability of the TAEs and 1/1 EIC.

3.1. Toroidal Alfven Eigenmodes

LHD operation scenarios with low magnetic field and
density show a rather large AE activity if the plasma is
strongly heated by the tangential NBIs, indicating an
efficient resonance between EPs and AEs. LHD inward
shifted configurations with magnetic fields intensities
of ≈ 0.5 T and high thermal β show bursting MHD
activity measured by the Mirnov coils in the frequency
range of 40 – 80 kHz, as well as an enhancement of the
EP fluxes detected by the tangential neutral particle
analyzer [58, 113]. FAR3d simulations were used to
reproduce the instabilities measured in the experiment,
identifying unstable n = 1 and 2 TAEs in the frequency
range of 50 to 80 kHz [88]. Figure 1, panel a,
shows wide TAE gaps covering all the minor radius
(continuum gaps are calculated using the Stellgap
code including the sound wave coupling [114]). The
n = 1 and 2 TAEs calculated in the simulations are
located in the inner-middle plasma region, consistent
with the observations. Panels b and c show the
growth rate and frequency of the n = 1 and n = 2
perturbations with respect to the EP β, identifying
a critical EP β = 0.012 for n = 1 TAE and < 0.01
for the n = 2 TAE destabilization. Similar results
are obtained in simulations performed by the MEGA
code [115, 116]. Panels d and e (f and g) shows
the evolution of the n = 1 (n = 2) perturbation
in a 2D poloidal cut for an EP β below and above
the critical value required for the TAE destabilization.
Plots on the left show interchange modes and TAEs
on the right. TAEs destabilization lead to a wider
perturbation, showing up-down asymmetric patterns
and the poloidal rotation of the mode. The study
also include the analysis of the TAE stability with

respect to the energy of the EP population and the
radial location of the EP density profile, linked to
the NBI voltage and deposition region, respectively.
This analysis can be considered one of the first studies
dedicated to explore optimized operation scenarios in
fusion devices with respect to the AE stability.

3.2. Energetic-ion-driven resistive interchange mode

EICs can be destabilized in low density and high
ion temperature discharges strongly heated by the
perpendicular and tangential NBIs. EIC are identified
by the Mirnov coils as sharp magnetic fluctuations
similar to the fishbone oscillations [59, 117]. EICs are
destabilized at the plasma periphery (around r/a ≈
0.88) above a given critical value of the helically
trapped EP β injected by the perpendicular NBIs,
showing a bursting phase once the mode frequency
up sweeps from 4 to 9 kHz. The EIC is an
energetic particle mode seeded by the 1/1 resistive
interchange mode [4, 5]. The RIC resonates with the
precessional motion of the helically trapped EP around
the frequency range of 10 kHz, destabilizing the EIC
that enhances the EP transport, reducing the plasma
heating efficiency of the perpendicular NBI limiting
the maximum thermal β of the discharge [118, 119].
The FAR3d code was used to analyze the resonance
of the helically trapped EP and reproduce the EIC
destabilization [120]. Figure 2, panel a and b, show
a 3D view of the plasma column and a 2D poloidal
cut of the 1/1 EIC density perturbation calculated in
the simulation (2/3 interchange mode and 3/4 EPM
are unstable in outer plasma region too). Panels c
and e show the critical β of the helically trapped EP
(β<sub>f,⊥</sub>) required to destabilize the 1/1 EIC is 0.0028
and an increment of the mode frequency from 3 to 5
kHz as the EP β rises. Below the EP β threshold a
RIC is unstable, panel d, evolving to a tearing like
eigenfunction structure once the EIC is destabilized,
panel f. Consequently, the linear simulations can
approximate the phase I and early phase II of the
EIC observed in the experiment (please see fig. 7 of
[59]). The analysis also suggests the possible role of
the helical couplings in the complex mode structure
and mode chirping, indicating the critical EP β and
frequency range of the EIC induced by the helical
family n = 1, 9, 11 are higher compared to the n = 1
simulations.

4. Optimization

Several techniques are tested to ameliorate or stabilize
AE / EPM instabilities, for example by modifying
the thermal plasma parameters, the magnetic field
configuration, the heating patterns and the presence
of different EP populations in the plasma. External
<figure><img src="image_2.png" /><figcaption>**Figure 1.** (a) Alfven continuum gaps for *n* = 1 to 5 toroidal mode families. Dashed horizontal lines indicate the frequency range and radial location of *n* = 1 (cyan) and 2 (orange) TAEs. (b) Growth rate and (c) frequency of the *n* = 1 and 2 TAEs with respect to the EP *β* value. 2D plots of the electrostatic potential perturbation for the *n* = 1 perturbation (d) below and (e) above the critical EP *β* to destabilize the TAE. Same plots for the *n* = 2 perturbation (f) below and (g) above the critical EP *β*. Reproduced courtesy of IAEA, IOPScience, Nuclear Fusion journal. Figure adapted from [88] and cite170. Copyright (2017 and 2022) IAEA.</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>**Figure 2.** (a) Isosurface and poloidal contours of the EP density perturbation. The yellow rectangle indicates the location of the poloidal contour. (b) Poloidal contour of the EP density perturbation. Instability growth rate (c) and frequency (e) for different values of the trapped EP *β* (fixed passing EP *β* = 0.002). Eigenfunction of the electrostatic potential of (d) resistive interchange mode and (f) 1/1 EIC. Reproduced courtesy of IAEA, IOPScience, Nuclear Fusion journal. Figure adapted from [120] and cite171. Copyright (2019 and 2021) IAEA.</figcaption></figure>

actuators as the injection of electron cyclotron waves (ECW) generate non inductive currents that modifies the magnetic field topology affecting the stability

of the AE / EPM and the Alfven gaps structure [121, 122]. Several experiments demonstrated the electron cyclotron current drive (ECCD) [101, 123, 124,
125, 126] can improve the stability of pressure/current gradient driven modes and AEs in tokamaks [127, 128, 129, 130, 131, 132] and stellarators [133, 134, 135, 136, 137]. Another example is the neutral beam current drive (NBCD) [138, 139] that causes broad modification of the magnetic field [140, 141, 142, 143] also affecting the AE / EPM and PGDM stability [144, 145]. Regarding the thermal plasma parameters, altering the thermal ion density is an effective way to change the AE / EPM stability because the plasma Alfvén velocity is modified ($V_A = B/\sqrt{N\mu_0m_pn_i}$ with $B$ the magnetic field intensity, $\mu_0$ the vacuum magnetic permeability, $n_i$ the thermal ion number density, $m_p$ the proton mass and $N$ the ion species mass/proton ratio). The EP resonance can be described by the ratio of the EP velocity ($V_f = \sqrt{T_{f,i}/Nm_pn_i}$ and $V_A$ with $T_{f,i}$ the energy of a given EP population) and the Alfvén velocity, leading to a resonance enhancement if $V_f/V_A$ is close to unity (or a fraction of) that is traduced to the destabilization of AE / EPM for a lower EP $\beta$. On top of that, the EP $\beta$ of the plasma decrease as the thermal ion density increase because the NBI deposition rate and EP slowing down time are altered. The thermal plasma temperature has a role in the AE / EPM stability too because the EP slowing down time ($\tau_{ep}$) is proportional to $T_e^{3/2}/n_e$ with $T_e$ ($n_e$) the thermal electron temperature (density) [146]. Indeed, it was demonstrated the capability of an off-axis electron cyclotron heating (ECH) to stabilize EICs by increasing the temperature of the plasma periphery [147]. Other techniques are dedicated to improve the efficiency of the plasma heating by modifying the operational regime of the heating mechanism inducing a weaker EP resonance, for example by balancing the NBI voltage and power while maintaining the same EP $\beta$ [148, 149] or by modifying the NBI deposition region [150, 151]. In addition, resonant magnetic perturbations (RMPs) can improve the stability of AEs [152].

In the following, several applications of the FAR3d code dedicated to explore optimized LHD operation scenarios with respect to the AE / EPM and PGDM stability are discussed. The techniques included in the analysis are the injection of ECCD and the generation of NBCD to change the magnetic field configuration, the modification of the thermal plasma parameters or the NBI operational regime.

### 4.1. External actuators

The effect of the ECCD injection on the AE stability was explored in several LHD configurations, applying the FAR3d code to calculate the AE growth rate considering the distortion induced in the magnetic configuration [153]. Figure 3, panels d and g, compare discharges without ECCD and with ctr-

ECCD, respectively (thermal plasma density and NBI heating pattern are the same in both discharges). The AE / EPM activity in the discharge with ctr-ECCD is significantly lower, thus the ECCD injection stabilizes all the modes at frequency ranges above 30 kHz. The spectrogram of the discharge without ECCD indicates unstable AE / EPM in the frequency range between 40 to 180 kHz identified as 1/2 EPM, 1/2 GAE and 2/3 – 2/4 TAE in the simulations (panels a to c). Low frequency modes calculated in the discharge with ctr-ECCD are identified as 1/1 EPM or PGDM and (f) 2/6 BAE or PGDM. Panels h and i compare the continuum gap structure in discharges injecting co- or ctr-ECCD, respectively. The ctr-ECCD case shows slender gaps with respect to co-ECCD configuration because the ctr-ECCD injection modifies the iota profile enhancing the magnetic shear. That means, AEs in ctr-ECCD configurations are confined in gaps of slender frequency bands showing a narrow eigenfunction width, easily damped by continuum and FLR effects, leading to a higher EP $\beta$ threshold for the AEs destabilization. Consequently, the ECCD is an useful tool to improve the AE stability in LHD discharges by reinforcing the magnetic shear and closing the Alfvén gaps in the radial location where the largest EP density gradients are located, minimizing the free energy required to trigger the AEs and increasing the destabilization EP $\beta$ threshold.

The stability of AEs and PGDM was also analyzed with respect to the net plasma current generated by the tangential NBIs in LHD inwards and outwards shifted configurations [154, 155]. Figure 4, panels a to c, show the variation of the AE activity in the outward shifted discharge 147372 due to an unbalanced NBI injection as the net plasma current increases. During the first part of the discharge a co-NBCD is generated followed by a second phase of ctr-NBCD that caused a down sweeping of the AE frequency. Panels d and e indicate the modification of the Alfvén gaps induced by different intensities of the co-NBCD that may explain the evolution of the AE stability along the discharge. FAR3d simulations indicate a decrease of the AE growth rate as the co-NBCD intensity increase, panels f and g, consistent with the weaker Mirnov coil signal observed in the frequency range of 200 kHz as the current intensity increases along the co-NBCD phase of the experiment. Figure j shows the iota profile distortion induced by the co-NBCD for different current intensities, leading to an up-shift of the profile and a modification of the magnetic shear, particularly enhanced in the middle plasma region. Consequently, the NBCD is an useful actuator to modify the AE and PGDM stability. Nevertheless, the analysis does not identify any configuration showing an improved AE and PGDM stability at the same time. A larger
<figure><img src="image_4.png" /><figcaption>**Figure 3.** Eigenfunction of the electrostatic potential perturbation of (a) 1/2 EPM, (b) 1/2 GAE and 2/3 – 2/4 TAE. (d) Spectrogram (magnetic data) of the discharge 138675 (no ECCD). The colored star symbols indicate the modes calculated by FAR3d. Eigenfunction of the electrostatic potential perturbation of (e) 1/1 EPM and (f) 2/6 BAE. (g) Spectrogram (magnetic data) of the discharge 138677 (ctr-ECCD). The colored star symbols indicate the modes calculated by FAR3d. Continuum gaps in configurations with co-ECCD (h) and ctr-ECCD (i). Reproduced courtesy of IAEA, IOPScience, Nuclear Fusion journal. Figure adapted from [153]. Copyright (2020) IAEA.</figcaption></figure>

co-NBCD in inward shifted configurations leads to a weaker AE activity although interchange modes are further destabilized. On the other hand, a more intense co-NBCD in outward shifted configurations barely modifies the PGDM growth rate and the AEs are less unstable. Figure 5 shows examples of dedicated experiments performed to improve the AE and PGDM stability thanks to the NBCD. Panels a to c indicate the 1/2 interchange mode is unstable during the ctr-NBCD phase of the discharge 167800 although stable in the co-NBCD phase. 1/2 interchange mode stabilization is caused by the modification of the magnetic configuration induced by the NBCD because the thermal $\beta$ is almost constant and even slightly increases once the 1/2 mode signal vanished. Panel d shows the evolution of the iota profile at different discharge times indicated in panels a to c by dashed vertical lines, measured using the Motional Stark Effect (MSE) diagnostic at $t = 5.5$ (red line) and 5.8 (blue line) s although using a VMEC equilibrium reconstruction at $t = 5.0$ s (black line) because no MSE data is available in the ctr-NBCD phase [156, 157]. The MSE data indicates the iota profile increases along the co-NBCD phase above 0.5, thus 1/2 rational surface

shows a transition from resonant to non resonant that may explain the 1/2 interchange mode stabilization. Panel e shows the growth rate of the $n = 1$ PGDMs linked to the most important rational surfaces of the equilibria, approximately reproducing the stability trends of 1/2 and 1/1 interchange modes along the discharge, in particular the 1/2 mode stabilization at $t = 5.8$ s once the 1/2 rational surface goes non resonant. Panels f to h indicate the stabilization of an AE in the frequency range of 70 kHz along the ctr-NBCD phase of the discharge 167805. No MSE data is available in the ctr-NBCD phase although measurements at the beginning of the co-NBCD phase indicate the iota profile at the inner plasma may decrease below 0.3 during the ctr-NBCD phase, panel i. Panel j shows the evolution of the continuum gaps as the current intensity in the ctr-NBCD phase increases and the iota profile down shifts. Configurations with down shifted iota (cases with $t_{min} = 0.1$ and 0.3) show slender gaps in frequency range and radial width as well as a reduction of the TAE and EAE gaps frequency range, leading to AEs with a narrower eigenfunction that may suffer a stronger stabilizing effect by continuum and FLR dampings. This result
may explain the AE stabilization above a given threshold of the current intensity along the ctr-NBCD discharge phase.

## 4.2. Thermal plasma

The 1/1 EIC stability was analyzed for different LHD operation scenarios with respect to the thermal plasma parameters, operational regime of the perpendicular NBI, magnetic field topology and intensity as well as thermal plasma and perpendicular NBI species [158]. Figure 6, panels a and b, show the 1/1 EIC is unstable below a given threshold of the thermal ion density ($n_i$) and thermal electron temperature ($T_e$). In particular, the 1/1 EIC is stable if the thermal plasma $\beta$ is above 0.25% (dotted purple lines). On the other hand, EICs are unstable if the thermal plasma $\beta$ is 0.2% with $n_i < 0.15 \cdot 10^{20} \text{ m}^{-3}$ and $T_e > 1.5$ keV or $n_i > 0.25 \cdot 10^{20} \text{ m}^{-3}$ and $T_e < 0.5$ keV. From the simulations data generated in the parametric study the stabilization threshold of the 1/1 EIC is calculated considering the non linear fit $n_i = AT_e^\alpha$ with $A = 0.23$ and $\alpha = -1.09$. The theoretical 1/1 EIC stabilization threshold shows a reasonable agreement with the experiments showing unstable 1/1 EICs, indicated in the figure by star symbols, because the main part of the stars are located below the theoretical threshold for the 1/1 EIC stabilization. Panels c and d indicate the decrease of the 1/1 EIC growth rate as the thermal ion density and thermal electron temperature increase while the growth rate of the RIC increases. Once the growth rate of the RIC is dominant, the 1/1 EIC is stabilized. Panels e analyzes the dependency of the 1/1 EIC with the thermal plasma species, showing a reduction of the growth rate as the atomic mass of the thermal plasma species increases, pointing out the resonance of the helically trapped EP is weaker because the 1/1 EIC stabilization is observed for a lower thermal ion density. Likewise, panel f shows the EP $\beta$ threshold to destabilize the 1/1 EIC is higher as the atomic number of the thermal plasma species increases. Consequently, 1/1 EIC stability improves in operation scenarios with higher thermal plasma atomic number species. It should be noted that several of the optimization trends identified in the analysis are consistent with the experimental observations [18, 147].

The effect of the thermal plasma density on the AE / EPM and PGDM stability was explored in different operation scenarios performing a set of FAR3d simulations and dedicated experiments [155]. Figure 7, panels a to d, show the MHD activity measured in discharges with averaged thermal electron densities ($<n_e>$) ranging from 0.2 to $3.2 \cdot 10^{19} \text{ m}^{-3}$ (thermal plasma $\beta$ from 0.21 to 1.46). Low density discharges, $<n_e> \le 0.8 \cdot 10^{19} \text{ m}^{-3}$, show robust AE activity in the frequency range of 40 to 70 kHz although stable

interchange modes, because a thermal $\beta = 0.61$ is below the PGDM destabilization threshold. On the other hand, discharges with $<n_e> \ge 2.1 \cdot 10^{19} \text{ m}^{-3}$ show stable AEs although unstable interchange modes with frequencies below 10 kHz, because a thermal $\beta = 1.39$ is above the PGDM destabilization threshold. Than means, optimization trends of AEs and PGDM with respect to the thermal plasma density are the opposite. That means, high thermal $\beta$ operations with rather large thermal ion density may lead to configuration with stable AEs and an efficient plasma heating although strongly unstable PGDMs that cause an unfavorable plasma confinement. Panel e compares the gap structure in configurations with low (blue line) and high (pink line) thermal plasma density, showing slender gaps in configurations with high thermal plasma density. Consequently, AEs stability improves in high density configurations because continuum and FLR damping effects are reinforced increasing the EP $\beta$ threshold. Panels f and g show the growth rate and frequency of the $n=1$ perturbation (black dots) and 1/2 interchange mode (red dots) with respect to the thermal ion density. Modifying the simulation thermal ion density changes the properties of the EP resonance given by the $V_f/V_A$ ratio, fixed $V_f$. The simulations indicate the stabilization of the $n=1$ AE / EPM if $n_i \ge 1.6 \cdot 10^{19} \text{ m}^{-3}$, consistent with the experiment observations. Calculating a regression of the $n=1$ AE / EPM growth rate indicates a non linear dependency with thermal ion density, $\gamma\tau_{A0} = 0.10 - 0.03n^{2.5}$ because the EP resonance is proportional to $\sqrt{n_i}$. The 1/2 interchange mode growth rate increases linearly with the thermal ion density, $\gamma\tau_{A0} = 0.08 + 0.03n^{1.07}$, explained by the linear dependency with the thermal plasma $\beta$. The green stars show the 1/2 interchange mode growth rate and frequency if the simulations include the decrease of the EP $\beta$ as the thermal ion density increases, that is to say, the simulations consider the decrease of the EP slowing down time and the reduction of the EP population with the thermal plasma density. The simulations indicate a decrease of the 1/2 interchange mode frequency as the thermal ion density increases, consistent with the experiment observations.

## 5. Multiple EP effects

In the burning plasma of future fusion reactors multiple EP populations will coexist, including fusion born alpha particles and EPs generated by different heating mechanism as NBI and ion cyclotron resonance heating (ICRH) [159, 160, 161, 162]. Multiple EP populations may have an important impact on the AE / EPM and thermal plasma stability. Such effects were originally observed in TFTR device DT experiments, plasma
<figure><img src="image_5.png" /><figcaption>Figure 4. AE stability with respect to the NBCD in the shot 147372. (a) Averaged $\beta$, (b) toroidal current with the NBI injection pattern and (c) raw signal of the magnetic probe. Continuum gaps in the configuration with (d) $I_p = 0$ kA/T and (e) 30 kA/T. Growth rate and frequency of the unstable modes calculated by FAR3d in the configurations with (f) $I_p = 0$ kA/T and (g) 30 kA/T. Instability mode number measured by the magnetic probe arrays in the (h) co-NBCD and (i) ctr-NBCD phases. (j) Iota profile for different $I_p$ values. Reproduced courtesy of IAEA, IOPScience, Nuclear Fusion journal. Figure adapted from [154]. Copyright (2020) IAEA.</figcaption></figure>

with a significant amount of fusion born alpha particles and NBI EPs. AEs destabilized by alpha particle were only observed once the NBI stops, indicating the lower energy NBI EPs have a stabilizing effect through Landau damping on the perturbation induced by the alpha particles [163, 164, 165]. Similar conclusions were obtained in the last JET D-T campaigns [166, 167]. This is explained by the EP mode stability dependency on the combined contribution of all the EP species, that is to say, the gradient sign and amplitude of each EP distribution function on the wave-particle resonance [168].

Multiple EP effects are explored in LHD plasma, particularly the AE / EPM stability resulting in the combination of tangential and perpendicular NBIs or tangential NBI lines with different operation regimes [120, 158, 169]. Figure 8, panel a and b, indicate a decrease of the 1/1 EIC growth rate and frequency if the power injection of the tangential NBI increases for a fixed injection intensity of the perpendicular NBI, that is to say, the passing EPs injected by the tangential NBI have a stabilizing effect on the perturbation induced by the helically trapped EPs generated by the

perpendicular NBI. Panel c indicates the time interval between EIC burst ($\Delta t$) measured in discharges with different NBI heating patterns fixing the perpendicular NBI injection power to $P_{P,NBI} = 13.5 \pm 0.5$ MW (blue diamonds) or $P_{P,NBI} = 17.5 \pm 0.5$ MW (red circles) and changing the tangential NBI injection intensity. Experimental evidence suggest $\Delta t$ increases as the tangential NBI power enhances. Linear regressions provide a first crude approximation of this trend: $\Delta t = 0.001 \cdot P_{T,NBI}$ (17.5 MW case) and $\Delta t = 0.0008 \cdot P_{T,NBI}$ (13.5 MW case). Panel d indicates different EP density profiles calculated using the analytic expression:

$$n_b(r) = \frac{(0.5(1 + \tanh(\delta_{flat} \cdot (r_{peak} - r))) + 0.02)}{(0.5(1 + \tanh(\delta_{flat} \cdot r_{peak}))) + 0.02)}$$

The location of the profile gradient is controlled by the parameter ($r_{peak}$) and the flatness by ($\delta_{flat}$). These set of EP density profiles are used in the analysis of the $n = 1$ and 2 AE stability for multiple EP configurations generated by the combination of two tangential NBIs with different operational regimes, panel e to h. In particular, the deposition region of the tangential NBI B is modified leading to different
<figure><img src="image_6.png" /><figcaption>**Figure 5.** Time evolution of the (a) averaged thermal $\beta$, (b) NBI heating pattern and toroidal plasma current and (c) normalized magnetic field perturbation induced by the most unstable modes in the discharge 167800. The black number in panel (b) indicates the local maximum of the toroidal current. The dashed vertical lines indicate the discharge times selected to analyze the evolution of the iota profile: $t = 5.0$ s (black), 5.5 s (red) and 5.8 s (blue). (d) Iota profile evolution for the selected discharge times. The green dashed line indicates the radial location of the 1/2 rational surfaces. (e) Growth rate of the unstable PGDMs calculated by FAR3d in the selected discharge times. Time evolution of the (f) averaged thermal $\beta$ and (g) NBI heating pattern and toroidal plasma current in the discharge 167805. The dashed vertical pink lines indicate the discharge period while the AE is unstable. The dashed vertical lines at $t = 5.5$ s (black) and 5.8 s (red) show the discharge times with iota profile measurement by the MSE diagnostic. (h) Magnetic probe signal. (i) Iota profiles measured by MSE diagnostic. (j) Alfvén gaps evolution as the iota profile decrease along the discharge ctr-NBCD phase. Reproduced courtesy of AIP Publishing, Physics of Plasma journal. Figure adapted from [155]. Copyright (2024) AIP Publishing.</figcaption></figure>

EP density profiles while the configuration of the tangential NBI A is fixed. In the study, the EP energy and $\beta$ are fixed to 48 keV and 0.0212 for both NBI lines. Here, the growth rate and frequency of the dominant $n = 1$ (black) and 2 (red) AEs are indicated by a solid line for the multiple EP simulations, a dashed line for the fixed tangential NBI A and a dotted line for the variable tangential NBI B. Multiple EP damping effects lead to the stabilization of the $n = 1$ AE if the EP density profile generated by the tangential NBI B has $r_{peak} > 0.5$ or $\delta_{flat} < 0.5$. On the other hand, $n = 2$ AE is further destabilized in multiple EP simulations except if $r_{peak} = 0.4$, configuration leading to a weak stabilizing effect but not the mode stabilization. Consequently, the analysis indicates the on-axis deposition of both NBI lines leads to the most unstable configuration, the combination of on-axis and off-axis NBI deposition is required to observe a stabilizing effect between beams. The study also indicates the most unstable AEs are observed in

configurations with the same EP energy for both NBI lines (same NBI voltage), trend further enhanced as the NBI power increases. Indeed, no optimizations trends are identified linked to the variation of the EP energy or $\beta$, only with respect of the EP density profile.

## 6. Saturation of AE/EPM and interchange modes

Non linear simulations are performed to study the saturation phase of an instability. In particular, the study of the AE / EPM saturation provides information of the induced turbulent EP transport, energy transfer towards the thermal plasma and redistribution between different toroidal or helical mode families, generation of zonal structures as shear flows and zonal current, nonlinear interaction between the electromagnetic fields linked to different EP populations in addition to other effects. Likewise, studying the saturation of PGDM is also useful to
<figure><img src="image_7.png" /><figcaption>**Figure 6.** (a) Growth rate and (b) frequency of the $n = 1$ perturbation in simulations with different values of the thermal plasma density and temperature linked to the $t = 1$ rational surface. The dashed purple lines indicate the iso-lines of the simulations with the same thermal $\beta$. The solid white line shows the transition between scenarios with dominant 1/1 EIC and RIC fitted by the non-linear curve $n = aT^b$. The stars show the EIC destabilized in LHD discharges with respect to the thermal plasma density and temperature at the $t = 1$ rational surface. The pink diamonds indicate the reference case of the theoretical analysis. Growth rate of the 1/1 perturbation in parametric scans with respect to the (c) thermal ion density and (d) electron temperature. The black line and black dots represent the 1/1 EIC and the blue line with blue diamonds the RIC. The pink rectangle indicates the reference case. (e) Growth rate of the 1/1 EIC for different thermal plasma densities (EP $\beta = 0.01$) in Hydrogen plasma with Hydrogen NBI (red line and dots), Helium + Hydrogen plasma with Hydrogen NBI (orange line and circles), Deuterium plasma with Deuterium NBI (green line and stars) and Deuterium + Helium plasma with Deuterium NBI (dark green line and stars). (f) Growth rate of the 1/1 EIC for different EP $\beta$ (fixed $n_i = 0.25 \cdot 10^{20} m^{-3}$) and thermal plasma species. Reproduced courtesy of IAEA, IOPScience, Nuclear Fusion journal. Figure adapted from [158]. Copyright (2020) IAEA.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 7.** Magnetic probe spectrogram in the discharges (a) 167793, (b) 167779, (c) 167799 and (d) 167790. The yellow dashed horizontal line indicates the upper frequency limit of the BAE gap if the electron temperature is 2 keV. The averaged density values are normalized to $10^{19} m^{-3}$. (e) Continuum gaps for the models with $n_i = 0.2$ (blue line) and $2 \cdot 10^{19} m^{-3}$ (pink line). Growth rate (f) and frequency (g) of $n = 1$ perturbation (black dots), 1/2 interchange mode (red dots) and 1/2 interchange mode including the decrease of the EP $\beta$ with the thermal ion density (green stars) for the discharge 167805. Reproduced courtesy of AIP Publishing, Physics of Plasma journal. Figure adapted from [155]. Copyright (2024) AIP Publishing.</figcaption></figure>

understand the energy redistribution between different modes and the generation of zonal structures in the experiment. Consequently, the saturation phase of

AE / EPM and PGDM should be analyzed to gain a improved comprehensions of LHD plasma heating efficiency and confinement.
<figure><img src="image_9.png" /></figure>

**Figure 8.** Growth rate (a) and frequency (b) of the 1/1 EIC for different $\beta$ values of the passing EP fixed the EP $\beta$ of the helically trapped EP to 0.015. (c) Time interval between EIC events for different tangential NBI injection intensities for a fixed value of the perpendicular NBI injection intensity: blue diamonds show the discharges with $P_{P,NBI} = 13.5 \pm 0.5$ MW and the red circles $P_{P,NBI} = 17.5 \pm 0.5$. The plot includes the linear regression $\Delta t = bP_{P,NBI}$ for a perpendicular NBI injection of $13.5 \pm 0.5$ MW (dashed blue line) and $17.5 \pm 0.5$ MW (dashed red line). (d) EP density profiles defined with respect to the parameters $r_{peak}$ (radial location of density profile gradient) and $\delta_{flat}$ (slope of density profile gradient). (e) Growth rate and (g) frequency of the dominant AEs in multiple EP configuration with 2 tangential NBI line if $r_{peak}$ of the NBI line B EP density profile is modified. (f) Growth rate and (h) frequency if $\delta_{flat}$ for the NBI line B EP density profile is modified. The solid lines show the result of multiple NBI simulations, the dotted lines the single variable NBI B simulations and the dashed lines the single fixed NBI A simulations. Reproduced courtesy of IAEA, IOPScience, Nuclear Fusion journal. Figure, panels a to c, adapted from [158]. Copyright (2020) IAEA. Reproduced courtesy of AIP Publishing, Physics of Plasma journal. Figure, panels d to h, adapted from [169]. Copyright (2019) AIP Publishing.

This section is dedicated to analyze the saturation phase of TAEs and 1/1 EIC linked to the destabilization of MHD and EIC bursts, respectively, leading to important EP losses and a lower tangential and perpendicular NBIs heating efficiency. In addition, the saturation of interchange modes in inward shifted configurations is explored, studying the destabilization of internal collapse and saw-tooth like oscillations, events that hamper the LHD performance by reducing the thermal plasma confinement and limiting the maximum thermal $\beta$ of the operation scenario.

## 6.1. MHD burst

TAEs saturation is studied to identify the driver of the MHD bursts [170]. Figure 9 panel a shows the spectrogram calculated using the poloidal magnetic

field perturbation along the non linear simulation, indicating TAE activity in the same frequency range than the experiment during the MHD burst, multiple modes in the frequency range of 40 to 80 kHz, as well as unstable interchange modes and zonal structures below 20 kHz. Panels b and c show the energy transfer from the 1/1 – 1/2 TAE to $n = 2$ and 3 toroidal mode families that cause the destabilization of the 2/3 – 2/4 and 3/4 – 3/5 TAEs (initial non linear phase). The overlapping between $n = 1$ to 3 TAEs in the inner-middle plasma region triggers the MHD burst, observed as a local maxima of the TAEs energy in the simulations. In addition, the energy transfer from the TAEs to the thermal plasma lead to the generation of zonal structures, particularly zonal flows. Panels d indicates a flattening of the pressure profile in the plasma region where the TAEs are destabilized,
indicating the TAE saturation causes a deformation of the thermal plasma flux surfaces. Panel e shows how the zonal currents induce excursion of the iota profile in the inner-middle plasma region, that is to say the TAEs saturation modifies locally the magnetic field structure. Panel f shows the partial lost of the EP population during the MHD burst. Panel g indicates the poloidal contour of the EP density and electrostatic potential perturbation, showing a redistribution and partial lost of the EP population as well as the generation of radial electric fields powering the shear flows, respectively. Consequently, the analysis shows LHD performance can be seriously hampered if multiple EP resonance overlap, leading to the destabilization of bursting events and large EP losses, particularly if a broad and intense perturbation channels energy non linearly to other mode families. Nevertheless, MHD burst can be avoided by improving the AE stability of the operation scenario using the techniques discussed in the previous section, preventing the destabilization of multiple TAEs with wide eigenfunction and their overlapping.

## 6.2. EIC burst

Nonlinear simulations were also performed to analyze the different 1/1 EIC saturation phases: unstable RIC followed by the 1/1 EIC destabilization above a given helically trapped EP density threshold (phase I), the further enhancement of the 1/1 EIC leading to a wider mode width and inward propagation (phase II), the bursting 1/1 EIC activity and frequency up sweeping to 9 kHz showing a complex mode structure between $r/a = [0.6, 0.85]$ (phase III) and the 1/1 EIC stabilization after large EP losses before the RIC is unstable again (phase IV) [171]. Figure 10, panel a, shows the evolution of the thermal pressure along the simulation that approaches the electron temperature fluctuation measured in the experiment (figure 3b in the reference [59]), indicating a wide perturbation caused by the 1/1 EIC as well as the destabilization of the 3/4 EPM and 2/3 interchange mode at $r/a = 0.75$ and 0.68, respectively, during the Phase II. Once the 1/1 EIC burst begins, vertical solid yellow line, the simulation shows a partial overlapping between the 1/1 EIC, 3/4 EPM and 2/3 interchange mode perturbations, consistent with the inward propagation of the instability observed in the experiment from $r/a = 0.88$ to 0.65 during the phase III. Panels b and c indicate a wide 1/1 EIC perturbation during the burst, in particular due to the destabilization of the 2/2 and 3/3 overtones, as well as a large energy channeling from the 1/1 EIC towards the thermal, observed in the plots as a large amplitude 0/0 mode. Panels d and e indicate an enhancement of the EP density perturbation at the peripheral plasma region

linked to the 1/1 EIC, 3/4 EPM and 2/3 interchange mode. Panels f to i show a flattening of the pressure and EP density profiles linked to the saturation of the 1/1 EIC and 3/4 EPM, indicating the deformation of the thermal plasma flux surfaces as well as the EP population redistribution and losses. Consequently, the non linear simulation suggests the large EP losses measured in the experiment can be induced once the 1/1 EIC perturbation width is large enough to overlap with 3/4 and 2/3 resonances. Thus, an improved LHD operation should avoid the transition from phases II to III and the destabilization of the 1/1 EIC bursts, identifying the width of the 1/1 EIC as well as the destabilization of the 3/4 and 2/3 modes as precursors of the burst event.

## 6.3. Internal collapse

Interchange modes in LHD inward shift configurations can be destabilized limiting the LHD device performance to confine the plasma [172, 173, 174], although showing a rather weak effect on the energy flows leaving the system [1]. There is experimental evidence of partial plasma collapses caused by the interchange mode 1/1 [63, 175] and 1/2 [61, 176] in LHD plasma as well as in other devices as Heliotron E [177, 178], leading to an important degradation of the plasma confinement. Non linear simulations are performed to analyze the internal collapses induced by the interchange mode 1/2 [179, 180]. Figure 11, panel a to h, indicate the evolution of the magnetic field structure during the internal collapse reproduced in the simulation. Initially, at $\Delta t = 0$ s, the magnetic surfaces are well defined and only 1/2 magnetic islands are observed in the middle plasma (panel a) showing slender magnetic field stochastic regions along the plasma (panel b). Once the internal collapse is triggered at $\Delta t = 0.01$ s, wide magnetic island triggered by different toroidal mode families overlap between the middle and inner plasma (panel c) generating extended stochastic region that strongly hamper the plasma confinement. After the internal collapse the interchange modes saturate and the width of the magnetic islands decrease, no longer overlapped, thus close magnetic surfaces recover reducing the radial extension of the stochastic regions and the plasma confinement improves (panels e and g). Panel i indicates large deformations of the pressure profile in the radial location where the 1/2 interchange mode is destabilized at the simulation time the internal collapse is triggered (blue line). Next, the profile flattening linked to the interchange mode 2/5 indicates the inwards propagation of the perturbation while the 1/2 interchange mode continues the saturation in the middle plasma (green line). Finally, the inner plasma is stabilized and the perturbation is again localized in the middle plasma region where the interchange modes
<figure><img src="image_10.png" /><figcaption>**Figure 9.** (a) Spectrogram of the poloidal component of the magnetic field perturbation along the nonlinear simulation. Temporal evolution of *n* = 1 to 3 TAEs energy during the (b) initial nonlinear and (c) bursting phases. The colored lines indicate the different TAEs: black 1/1–1/2, red 2/3–2/4 and blue 3/4–3/5. Temporal evolution of the (d) normalized pressure, (e) iota profile and (f) normalized EP density. The vertical dashed lines indicate the radial location of the main resonant rational surfaces: red 1/2, light green 3/5, dark yellow 2/3, green 3/4 and blue 1/1. (g) Poloidal contour of the normalized perturbations of the EP density and electrostatic potential during the MHD burst. Reproduced courtesy of IAEA, IOPScience, Nuclear Fusion journal. Figure adapted from [170]. Copyright (2022) IAEA.</figcaption></figure>

1/2 and 2/3 continues the saturation (red line). Panel j shows negligible variations of the iota profile during the internal collapse thus the radial location of the rational surface barely change along the simulation. Panels k to o indicate a strong deformation of the flux surfaces once the internal collapse is triggered, leading to the generation of torn surfaces and the outwards flux of plasmoid linked to a partial lost of the plasma confinement (panel m). Flux surfaces recover after the internal collapse. Consequently, improved operation scenarios in LHD must avoid the destabilization of internal collapses linked to broad PGDM, particular 1/1 or 1/2 interchange modes, as well as the overlapping between instabilities induced by different toroidal mode families. It should be noted that the size of the magnetic islands and stochastic regions are overestimated in the simulations because the plasma resistivity is larger compared to the real plasma. The overlapping between different interchange modes is weaker in real LHD discharges, the observations indicate the destabilization of internal collapse once the width of

a single mode overcomes a given threshold. Nevertheless, islands overlapping could be an issue for the LHD plasma confinement in operations with large thermal $\beta$ because, as the Minov coils measurements show, multiple interchange modes can be unstable along the discharge at the same time.

### 6.4. Saw-tooth like events

High density plasma ($n_e > 10^{20}$ m$^3$) fueled using pellets can show saw-tooth like oscillations linked to an enhancement of the pressure gradient in the plasma core. Tangential viewing soft x-ray camera can measure three islands nearby the magnetic axis as precursor of the saw-tooth oscillations, evolving to a two island pattern located in the middle plasma after the saw-tooth crash [64, 65]. A set of non linear simulations are performed to reproduce 1/3 saw-tooth oscillations in LHD [181, 182]. Figure 12, panels a to f, indicate three islands linked to the 1/3 rational surface nearby the plasma core (panels a and c) and
<figure><img src="image_11.png" /><figcaption>**Figure 10.** (a) Spectrogram of the thermal pressure perturbation along the nonlinear simulation. Temporal evolution of $n = 1$ to 3 TAEs energy during the (b) initial nonlinear and (c) bursting phases. The colored lines indicate the different TAEs: black $1/1-1/2$, red $2/3-2/4$ and blue $3/4-3/5$. Eigenfunction of the electrostatic potential during the MHD burst at (b) $t = 39000\tau_{A0}$ and (c) $t = 40000\tau_{A0}$. Poloidal contour of the EP density perturbation during the MHD burst at (d) $t = 39000\tau_{A0}$ and (e) $t = 40000\tau_{A0}$. Evolution of the normalized (f) EP density and (g) pressure linked to the 1/1 EIC. Evolution of the normalized (f) EP density and (g) pressure linked to the 3/4 EPM. Reproduced courtesy of IAEA, IOPScience, Nuclear Fusion journal. Figure adapted from [171]. Copyright (2021) IAEA.</figcaption></figure>

an extension of the plasma stochastic region from the middle to the inner plasma (panels b and d). After the saw-tooth crash in the simulation, the inner plasma region recovers the close magnetic surfaces although magnetic islands linked to the rational surface 1/2 remain in the middle plasma (panel e) while the width of the stochastic regions decreases (panel f), consistent with the experimental observations. Panel g indicates a deformation of the pressure profile caused by the interchange modes destabilized in the inner plasma, particularly interchange modes linked to the rational surfaces 1/3 and 2/5, as well as the 1/2 interchange mode in the middle plasma. Panel h shows a rather large excursion of the iota profile in the inner plasma induced by the zonal currents generated during the interchange modes saturation, leading to the 1/3 rational surface transition from resonant to non resonance once the saw-tooth crash is triggered in the simulation. Panel i shows a weaker deformation of the flux surfaces during a saw-tooth like event compared to the internal collapse, pointing out the deterioration of the plasma confinement is smaller. In summary, saw-tooth crashes limit the thermal $\beta$ of LHD operations by inducing a recursive decrease of the device confinement performance once the pressure gradient in the plasma

overcomes a given threshold. Avoiding a resonant 1/3 rational surface nearby the plasma core may prevent the most deleterious effects of the saw-tooth like events on the LHD performance. Such optimized scenario requires stable or marginal unstable interchange modes in the inner plasma to avoid the down shift of the iota profile by the zonal current induced during the instabilities saturation.

## 7. Ongoing studies

New research lines dedicated to calculate the linear and nonlinear evolution of AE / EPM and interchange modes in LHD plasma are performed. In particular, shear flows induced by AE / EPM and consequence on the thermal plasma turbulence, the AE stability in discharges with heating patterns combining NBI and ICRH, optimizations studies mixing ECCD and ECH injection with NBCD, bursting activity induced by AE / EPM, characterization of 1/1 instabilities in the plasma periphery and the stabilizing effect of the EPs on PGDM and tearing modes are explored, between other topics.

There is experimental evidence that indicates a stabilizing or destabilizing effects of the EP particles
<figure><img src="image_12.png" /><figcaption>**Figure 11.** Simulation of a 1/2 internal collapse. Evolution of the Poincare plots of the magnetic field including the main resonance surfaces (a, c, e and g) and all the modes in the simulations (b, d, f, h). Evolution of the (i) pressure profile and (j) iota profile. The colored dashed lines indicate the radial location of the main rational surfaces. (k, m, n, o) Evolution of the poloidal section of the pressure. Red circles indicate plasma regions where flux surfaces are torn. Reproduced courtesy of IAEA, IOPScience, Nuclear Fusion journal. Figure adapted from [179, 180]. Copyright (2012 and 2014) IAEA.</figcaption></figure>

on the stability of thermal plasma instabilities, for example by increasing or reducing the waiting time between saw-tooth linked to the destabilization of kinks or tearing modes in JET, ASDEX, EAST and DIII-D plasma [183, 184, 185] or by damping the growth rate of interchange and ballooning modes in LHD plasma [186, 187, 188]. This is explained by the resonance between EPs and plasma waves, active even if the EP population is below the EP $β$ threshold required to destabilize EP-driven modes [189, 190, 191]. Figure 13 shows some preliminary results of the study dedicated to analyze EP effects on the PGDM stability in LHD plasma, particularly in the linear and saturation phases of the RIC. Panel a indicates the growth rate and frequency of the $n = 1$ perturbation, showing an unstable RIC below the EP $β = 0.004$ required to destabilize the 1/1 EIC. Once the 1/1 EIC is destabilized, increasing the EP $β$ leads to a faster growing EIC. On the other hand, below the EIC EP $β$ threshold, increasing the EP population density causes a decrease of the RIC growth rate, pointing out the stabilizing effect of the EPs on the RIC. Panel b and

c indicate the different eigenfunction structure of the RIC (green box) and 1/1 EIC (red box). Panels d and e show the evolution of the poloidal magnetic field perturbation in nonlinear simulations with and without the EPs contribution for configuration with unstable RIC (reason why the EP $β$ is below the threshold to trigger the 1/1 EIC). The simulation without EPs indicates perturbations one order of magnitude higher compared to the simulation with EPs. Panels f and g also show the RIC amplitude is one order of magnitude larger in the simulation without EPs. Consequently, a first exploratory analysis may indicate there is an important stabilization effect induced by the EPs on the RIC, result consistent with the rather weak MHD activity in LHD inward shifted operations for thermal $β$ values up to 5% [74, 75]. Consequently, the presence of an EP population may explain why previous theoretical studies deviates from the observations predicting stronger MHD instabilities in LHD high $β$ inward shifted configurations not observed in the experiments [72, 192, 193]. Further analysis are required to generalize and verify these
<figure><img src="image_13.png" /><figcaption>**Figure 12.** Simulation of a 1/3 saw-tooth like event. Evolution of the Poincare plots of the magnetic field including the main resonance surfaces (a, c, e) and all the modes in the simulations (b, d, f). Evolution of the (g) pressure profile and (h) iota profile. The colored dashed lines indicate the radial location of the main rational surfaces. (i) Evolution of the poloidal section of the pressure. Reproduced courtesy of IAEA, IOPScience, Nuclear Fusion journal. Figure adapted from [181, 182]. Copyright (2012 and 2014) IAEA.</figcaption></figure>

preliminary results.

## 8. Conclusions and discussion

The gyro-fluid code FAR3d is a fast tool dedicated to study the stability of AE / EPM, pressure gradient and current driven modes in nuclear fusion devices, including Tokamaks and Stellarators. Computation efficiency is essential to reduce the gap between numerical analysis and short term studies required in the experimental campaign, reason why reduced model as FAR3d are demanded.

FAR3d code is routinely used to reproduce the AE / EPM activity measured in LHD plasma. The analysis dedicated to study the stability of TAEs and 1/1 EIC provided a reasonable agreement with the experimental observations, identifying the same radial location, dominant mode numbers, frequency range and eigenfunction structure. Further studies were dedicated to the saturation of these instabilities, linked to bursting activity observed in the experiments, namely the MHD and EIC bursts. Nonlinear simulations are able to reproduce the different 1/1 EIC phases and the onset of the MHD burst, identifying the role of the EP resonances overlapping as the possible driver of the bursting activity in LHD plasma. Such information is useful to develop strategies dedicated

to prevent the destabilization of bursting events. Optimized LHD operations must avoid the large EP losses and inefficient plasma heating performance linked to the bursts, limiting the maximum thermal $\beta$ of the discharges.

Fast parametric studies performed by FAR3d code can be used to identify optimized operation regimes with respect to LHD key operation parameters. In particular, improved operation scenarios with reduced AE / EPM and PGDM activity are explored considering the linear and saturation phases of the instabilities. The optimization parameters studied in LHD plasma include the thermal plasma configuration (thermal ion density and electron temperature), the magnetic field configuration (application of ECCD and NBCD), multiple EP population effects (combination of tangential and perpendicular NBIs or tangential beams with different operation regimes) and the NBI operation regime (injection power, voltage and deposition region). Several of the improved operational scenarios identified theoretically were validated in dedicated experiments.

The thermal plasma configuration can be adapted to avoid the destabilization of 1/1 EICs by increasing the thermal ion density, electron temperature and the atomic number of the thermal plasma species. The 1/1 EIC destabilization threshold with respect to the
<figure><img src="image_14.png" /><figcaption>**Figure 13.** (a) Growth rate and frequency of the $n = 1$ perturbation for different EP $\beta$ values. Eigenfunction of the electrostatic potential perturbation of the (b) RIC (green box) and (c) 1/1 EIC (red box). Evolution of the poloidal magnetic field perturbation in nonlinear simulations if (d) EP $\beta = 0$ and (e) 0.001. Eigenfunction of the electrostatic potential perturbation in the nonlinear simulation if (f) EP $\beta = 0$ and $t = 5.5 \cdot 10^4 \tau_{A0}$ and (g) EP $\beta = 0.001$ and $t = 2.1 \cdot 10^6 \tau_{A0}$.</figcaption></figure>

thermal ion density and electron temperature was calculated using a data base of FAR3d simulations, as well as the critical EP $\beta$ required in LHD operations with Hydrogen or Deuterium thermal plasma, showing a fair agreement with the experiments. The analysis of the AE / EPM and PGDM stability with respect to the thermal ion density identified opposite optimization trends. A higher thermal ion density improves the AE / EPM stability by reducing the EP population (smaller slowing down time), weakening the EP resonance (modifying the plasma Alfven velocity) and enhancing the continuum / FLR damping effects (slender gaps in frequency range and radial width). On the other hand, a larger thermal $\beta$ causes an enhanced PGDM activity (higher pressure gradients). Consequently, LHD configurations with optimized plasma heating efficiency (stable AE/EPM) and confinement (stable PGDM) require operating at high thermal ion densities and the application of external actuators to avoid strongly unstable interchange mode (particularly RIC and 1/2 interchange mode) in inward shifted configurations and ballooning modes in outward shifted configurations. The application of RMP, NBCD, ECCD as well as off-axis ECH, ICRH and perpendicular beams to heat the plasma periphery are useful tools to improve the PGDM stability in high thermal ion experiments [147, 152, 175, 194, 195].

The effect of external actuators as the ECCD and NBCD were explored analyzing how the modification of the magnetic field by non induce current can improve the stability of AE / EPM and interchange modes. In

both cases, the simulations indicate the enhancement of the magnetic shear in the radial location where the AE / EPM are unstable can stabilize the modes by modifying the Alfven continuum structure, particularly by reducing the frequency range and radial width of the gaps leading to an enhancement of the continuum and FLR damping effects, traduced to a higher EP $\beta$ threshold for the EP driven modes destabilization. This is the case for the ctr-ECCD and ctr-NBCD configurations studied in LHD plasma. It should be noted that the ECCD injection is a more accurate actuator compared to the NBCD, inducing a localized modification of the iota profile compared to the broad alterations caused by the NBCD. On the other hand, NBCD is an useful tool to improve the stability of PGDM, particularly the 1/2 interchange mode, inducing an up-shift of the iota profile and a transition of the 1/2 rational surfaces from resonant to non-resonant, reflected in a important improvement of the LHD performance.

Far3d code is also applied to analyze the stabilizing effect of the passing particle generated by the tangential NBI on the perturbation induced by the helically trapped EP injected by the perpendicular NBI. The experimental evidence suggest the EIC burst waiting time increases with the tangential NBI power (larger population of passing EP), reason why the 1/1 EIC linear stability was studied with respect to the passing EP $\beta$ for a fixed helically trapped EP $\beta$. The simulations show a decay of the 1/1 EIC growth rate as the passing EP population increases,
consistent with the observations. Nevertheless, nonlinear simulations must be performed to confirm this result in the saturation phase of the instability. In addition, multiple EP damping effects were explored in configurations combining two tangential beam lines with different operation regimes. The analysis indicate that modifying the radial deposition and spreading of different beam lines (proportional to the radial location and stiffness of the EP density profile) may provide a new stabilizing mechanism for the AE/EPM, leading to the stabilization of $n = 1$ AEs for a given range of parameters for the configuration analyzed.

Nonlinear simulations of the interchange mode saturation phase were dedicated to analyze the internal collapse and saw-tooth like events in LHD inward shifted configurations. The simulations of the internal collapse were main focused on configurations with a strongly unstable 1/2 interchange mode, leading to wide magnetic island induced in the middle plasma region where the 1/2 rational surface resonates. The internal collapse reproduced in the simulations indicate the large decrease of the LHD confinement performance is linked to the generation of wide magnetic field stochastic region between the inner and middle plasma, caused by the overlapping of the 1/2 magnetic islands with the perturbations induced by other rational surfaces. It should be noted that there is not experimental evidence of internal collapses caused by the overlapping of multiple interchange modes, only if the 1/1 or 1/2 modes width overcomes a given threshold. The discrepancy between modeling and observations is related to the magnetic Lundquist number used in the simulations, reproducing a plasma with a resistivity significantly larger than the real case, leading to an overestimation of the interchange modes growth rates and the triggering of instabilities not observed in the experiment. On top of that, the stabilizing effect of the EPs on PGDM as well as the thermal ion FLR damping effect were not considered in the simulations. Nevertheless, the destabilization of a broad 1/2 interchange mode that evolves to a tearing mode, inducing the re-connection and stochastization of the magnetic field is a reasonable explanation of the experimental observations. The simulations dedicated to study the saw-tooth like oscillations, particularly the events caused by the rational surface 1/3, show a reasonable agreement with the experiment reproducing the 1/3 perturbation measured nearby the magnetic axis by the tangential viewing soft x-ray camera as precursor of the saw-tooth, as well as the 1/2 perturbation observed after the crash. Indeed, a synthetic diagnostic developed to calculate the plasma emissivity using the simulations data approximately reproduces the soft x-ray collapse during the saw-tooth crash in the inner plasma. The simulations indicate the

decrease of the LHD performance during the saw-tooth like oscillations could be caused by the stochastization of the magnetic field in the inner plasma region once the 1/3 rational surface resonates nearby the magnetic axis, perturbation that propagates outwards destabilizing the middle plasma region and the onset of the 1/2 interchange mode after the crash, consistent with the experiment observations.

Several research lines involving FAR3d code and LHD plasma are ongoing, particularly optimization studies dedicated to improve the LHD plasma heating efficiency and confinement with respect to the activity of AE / EPM and interchange modes. In addition, mode characterization analysis in plasma heated by the combination of NBI and ICRH are being performed. Another topic is the effect of the EPs on the stability of PGDM and tearing modes. Preliminary analysis indicate EPs have a stabilizing effect on the interchange modes linear and saturation phases, result consistent with the experimental evidence and numerical studies performed by other authors, a possible explanation of why LHD inward shifted configurations remains MHD stable, or marginal unstable, in operation scenarios with a thermal $\beta$ up to the 5%.

## Acknowledgments

The authors want to thank the LHD technical staff for their contributions in the operation and maintenance of LHD. This work was supported by the US DOE under grant DE-FG02-04ER54742, Comunidad de Madrid under the project 2019-T1/AMB-13648, U.S. Department of Energy Office Contract DE-AC05-00OR22725 with UT-Battelle LLC and NIFS07KLPH004.

[1] Watanabe, K. Y. et al *Nucl. Fusion*, **45**, 1247, (2005).

[2] Sakakibara, S. et al *Nucl. Fusion*, **50**, 124014, (2008).

[3] Sakakibara, S. et al *Fusion Sci. Technol.*, **58**, 176, (2010).

[4] Sakakibara, S. et al *Plasmas Phys. and Control. Fusion*, **50**, 124014, (2008).

[5] Watanabe, K. Y. et al *Phys. Plasmas*, **18**, 056119, (2011).

[6] Toi, K. et al *Nucl. Fusion*, **44**, 217, (2004).

[7] Toi, K. et al *Phys. Plasmas*, **12**, 020701, (2005).

[8] Ueda, R. et al *Phys. Plasmas*, **21**, 052502, (2014).

[9] Yamamoto, S. et al *27th IAEA FEC, EX/1-3Ra* (2018).

[10] Toi, K. et al *Nucl. Fusion*, **40**, 1349, (2000).

[11] Yamamoto, S. et al *Phys. Rev. Lett.*, **91**, 245001, (2003).

[12] Toi, K. et al *Plasma Phys. Control. Fusion*, **46**, S1, (2004).

[13] Yamamoto, S. et al *Nucl. Fusion*, **45**, 326, (2005).

[14] Higaki, H. et al *J. Plasma Fusion Res.*, **1**, 34, (2006).

[15] Toi, K. et al *Fusion Science and Technology*, **58**, 186, (2010).

[16] Nagaoka, K. et al *Nucl. Fusion*, **51**, 083022, (2011).

[17] Chen, W. et al *J. Plasma Fusion Res.*, **7**, 2402079, (2012).

[18] Ohdachi, S. et al *27th IAEA Fusion Energy Conference, EX/1-3Rb*, (2018).

[19] Takemura, Y. et al *J. Nucl. Fusion*, **52**, 102001, (2012).

[20] Ueda, R. et al *Phys. Plasmas*, **21**, 052502, (2014).
[21] Sakakibara, S. et al *J. Nucl. Fusion*, **55**, 083020, (2015).

[22] Takemura, Y. et al *J. Nucl. Fusion*, **59**, 066036, (2019).

[23] Ogawa, K. et al *Nucl. Fusion*, **50**, 084005, (2010).

[24] Ogawa, K. et al *Nucl. Fusion*, **53**, 053012, (2013).

[25] Ogawa, K. et al *Nucl. Fusion*, **59**, 076017, (2019).

[26] Ogawa, K. et al *Nucl. Fusion*, **60**, 112011, (2020).

[27] Osakabe, M. et al *Fusion Sci. Technol.*, **58**, 131, (2010).

[28] Nuga, H. et al *J. Plasma Fusion Res.*, **14**, 3402075, (2019).

[29] Nuga, H. et al *J. Plasma Phys.*, **86**, 815860306, (2020).

[30] Ogawa, K. et al *Nucl. Fusion*, **62**, 112001, (2022).

[31] Wilson, J. R. et al *Bull. Am. Phys. Soc.*, **37**, 1380, (1992).

[32] Sharapov, S. E. et al *Nucl. Fusion*, **39**, 373, (1999).

[33] Heidbrink, W. W. et al *Nucl. Fusion*, **31**, 1635, (1992).

[34] Duong, H. H. et al *Nucl. Fusion*, **33**, 749, (1993).

[35] Kusama, Y. et al *Nucl. Fusion*, **39**, 1837, (1999).

[36] Shinohara, K. et al *Nucl. Fusion*, **42**, 942, (2002).

[37] Castejón, F. et al *Plasma Phys. Control. Fusion*, **58**, 094001, (2016).

[38] Castejón, F. et al *Nucl. Fusion*, **57**, 102022, (2017).

[39] Garcia-Munoz, M. et al *Plasma Phys. Control. Fusion*, **054007**, 942, (2019).

[40] Heidbrink, W. W. *Plasma Phys.*, **15**, 055501, (2008).

[41] Todo, Y. et al *Rev. Mod. Plasma Phys.*, **3**, 1, (2019).

[42] Kieras, C. E. et al *J. Plasma Phys.*, **28**, 395, (1982).

[43] Cheng, C.Z. et al *Phys. Fluids*, **29**, 3695, (1986).

[44] Wang, Z. et al *Phys. Plasmas*, **22**, 022509, (2015).

[45] Betti, R. et al *Phys. Fluids B*, **4**, 1465, (1992).

[46] Kramer, G. J. et al *Phys. Rev. Lett.*, **80**, 2594, (1998).

[47] Betti, R. et al *Phys. Fluids B*, **3**, 1865, (1991).

[48] Bass, E. M. et al *Phys. Plasmas*, **20**, 012508 (2013).

[49] Turnbull, A.D. et al *Phys. Fluids B*, **5**, 2546, (1993).

[50] Deng, W. *Nucl. Fusion*, **52**, 043006, (2012).

[51] Kimura, H. et al *Nucl. Fusion*, **38**, 1303, (1998).

[52] Evans, T. E. et al *Phys. Rev. Lett.*, **53**, 1743, (1984).

[53] Chen, L. et al *Phys. Plasmas*, **1**, 1519, (1994).

[54] Ido, T. et al *Nucl. Fusion*, **51**, 073046, (2011).

[55] Ido, T. et al *Nucl. Fusion*, **55**, 083024, (2015).

[56] Wang, H. et al *Phys. Plasmas*, **22**, 092507, (2015).

[57] Wang, H. et al *Phys. Rev. Lett.*, **120**, 175001, (2018).

[58] Osakabe, M. et al *Nucl. Fusion*, **46**, S911, (2006).

[59] Du, X. D. et al *Phys. Rev. Lett.*, **114**, 155003, (2015).

[60] Weller, A. et al *Fusion Sci. Technol.*, **50**, 158, (2006).

[61] Kobayashi, T. et al *Nucl. Fusion*, **60**, 036017, (2020).

[62] Okamoto, M. et al *Nucl. Fusion*, **61**, 046005, (2021).

[63] Takemura, Y. et al *Phys. Plasmas*, **29**, 092505, (2022).

[64] Ohdachi, S. et al *21 IAEA Fusion Energy Conference*, 16 - 21 October Chengdu, China, EX/P8-15, (2006).

[65] Ohdachi, S. et al *22nd IAEA Fusion Energy Conference*, 13-18 October Geneva, Switzerland, EX/8-1Rb, (2008).

[66] Tsumori, K. et al *Plasma Sci. Technol.*, **8**, 24, (2006).

[67] Takeiri, Y. et al *Fusion Sci. Technol.*, **58**, 482, (2010).

[68] Nagaoka, K. et al *J. Plasma Fusion Research*, **3**, S1013, (2008).

[69] Ogawa, K. et al *Nucl. Fusion*, **52**, 094013, (2012).

[70] Du, X. D. et al *Nucl. Fusion*, **56**, 016002, (2016).

[71] Yamada, H. et al *Plasma Phys. Control. Fusion*, **43**, A55, (2001).

[72] Sato, M. et al *Nucl. Fusion*, **57**, 126023, (2017).

[73] Motojima, O. et al *Nucl. Fusion*, **43**, 1674, (2003).

[74] Komori, K. et al *Nucl. Fusion*, **49**, 104015, (2009).

[75] Yamada, H. et al *Nucl. Fusion*, **51**, 094021, (2011).

[76] Varela, J. et al *Plasma Fusion Res.*, **6**, 1403013, (2011).

[77] Miyazawa, J. et al *Plasma Fusion Res.*, **3**, 1047, (2008).

[78] Ohdachi, S. et al *Contrib. Plasma Phys.*, **50**, 552, (2010).

[79] Ohdachi, S. et al *Nucl. Fusion*, **57**, 066042, (2017).

[80] Watanabe, T. et al *Nucl. Fusion*, **46**, 291, (2006).

[81] Tanaka, K. et al *Fusion Sci. Technol.*, **51**, 97, (2007).

[82] Charlton, L. A. et al *Journal of Comp. Physics*, **63**, 107, (1986).

[83] Charlton, L. A. et al *Journal of Comp. Physics*, **86**, 270, (1990).

[84] Hedrick, C. L. et al *Phys. Fluids B*, **4**, 3869, (1992).

[85] Spong, D. A. et al *Phys. Plasmas*, **1**, 1503, (1994).

[86] L. Garcia *Proceedings of the 25th EPS International Conference, Prague, 1998*, VOL. 22A, Part II, p. 1757.

[87] Spong, D. A. *Plasma Fusion Res.*, **9**, 3403077, (2014).

[88] Varela, J. *Nucl. Fusion*, **57**, 126019, (2017).

[89] Hirshman, S. P. et al *Phys. Fluids*, **26**, 3553, (1983).

[90] Garcia, L. et al *Phys. Fluids B*, **2**, 2162, (1990).

[91] Varela, J. et al *Frontiers of Physics*, under review, (2024).

[92] Taimourzadeh, S. et al *Nucl. Fusion*, **59**, 066006, (2019).

[93] Pace, D. C. et al *Phys. Plasmas*, **25**, 056109, (2018).

[94] Varela, J. et al *Nucl. Fusion*, **58**, 076017, (2018).

[95] Varela, J. et al *Nucl. Fusion*, **59**, 046017, (2019).

[96] Huang, J. et al *Nucl. Fusion*, **60**, 126007, (2020).

[97] Ghai, Y. et al *Nucl. Fusion*, **61**, 126020, (2021).

[98] Varela, J. et al *Nucl. Fusion*, **57**, 126019, (2017).

[99] Cappa, A. et al *Nucl. Fusion*, **61**, 066019, (2021).

[100] Eliseev, E. et al *Physics of Plasma*, **28**, 072510, (2021).

[101] Yamamoto, S. et al *Nucl. Fusion*, **60**, 066018, (2020).

[102] Varela, J. et al *Nucl. Fusion*, **60**, 112015, (2020).

[103] Varela, J. et al *Nucl. Fusion*, **63**, 026009, (2023).

[104] Wang, X. H. et al *Nucl. Fusion*, **64**, 016028, (2024).

[105] Spong, D. A. et al *Nucl. Fusion*, **61**, 116061, (2021).

[106] Varela, J. et al *Plasma Phys. Control. Fusion*, **65**, 125004, (2023).

[107] Varela, J. et al *Phys. Plasmas*, **26**, 062502, (2019).

[108] Varela, J. et al *Nucl. Fusion*, **59**, 076036, (2019).

[109] Varela, J. et al *Nucl. Fusion*, **62**, 036005, (2022).

[110] Varela, J. et al *Nucl. Fusion*, **60**, 096009, (2020).

[111] Varela, J. et al *Nucl. Fusion*, **61**, 026023, (2021).

[112] Ortiz, J. et al *Nucl. Fusion*, **63**, 056010, (2023).

[113] Kamio, S. et al *Nucl. Fusion*, **60**, 112002, (2020).

[114] Spong, D. A. et al *Plasma Phys.*, **10**, 3217, (2003).

[115] Todo, Y. et al *Plasma Phys.*, **24**, 081203, (2017).

[116] Seki, R. et al *Nucl. Fusion*, **59**, 096018, (2009).

[117] Nave, M. F. F. et al *Nucl. Fusion*, **31**, 697, (1991).

[118] Ogawa, K. et al *Nucl. Fusion*, **58**, 044001, (2018).

[119] Ogawa, K. et al *Plasma Phys. Control. Fusion*, **60**, 044005, (2018).

[120] Varela, J. et al *Nucl. Fusion*, **59**, 046008, (2019).

[121] Silva Rosa, P. R. et al *Plasma Phys. Control Fusion*, **42**, 755, (2000).

[122] Prater, R. et al *Phys. of Plasma*, **11**, 2349, (2004).

[123] Sauter, O. et al *Phys. Rev. Lett.*, **84**, 3322, (2000).

[124] Suzuki, T. et al *Nucl. Fusion*, **44**, 699, (2004).

[125] Motojima, G. et al *Nucl. Fusion*, **47**, 1045, (2007).

[126] Fernandez, A. et al *Fusion Sci. Technol.*, **53**, 254, (2008).

[127] Nagasaki, K. et al *Nucl. Fusion*, **45**, 1608, (2005).

[128] Prater, R. et al *Nucl. Fusion*, **47**, 371, (2007).

[129] Fasoli, A. et al *Plasma Phys. Control. Fusion*, **39**, B287, (1997).

[130] Felici, F. et al *Nucl. Fusion*, **52**, 074001, (2015).

[131] Petty, C. C. et al *Nucl. Fusion*, **57**, 116057, (2017).

[132] Kim, J. et al *Nucl. Fusion*, **62**, 026029, (2022).

[133] Sallander, E. et al *Nucl. Fusion*, **40**, 1499, (2000).

[134] Castejon, F. et al *Plasma Phys. Control Fusion*, **47**, B53, (2005).

[135] Nagaoka, K. et al *Nucl. Fusion*, **53**, 072004, (2013).

[136] Melnikov, A. V. et al *Nucl. Fusion*, **56**, 112019, (2016).

[137] Garcia, L. et al *Phys. of Plasma*, **23**, 062319, (2016).

[138] Murakami, M. et al *Nucl. Fusion*, **51**, 103006, (2011).

[139] Pengxiang, C. et al *Plasma Sci. Technol.*, **18**, 1064, (2016).

[140] Chapman, I. T. et al *Nucl. Fusion*, **52**, 063006, (2012).

[141] Murakami, M. et al *Nucl. Fusion*, **49**, 065031, (2009).

[142] Isayama, A. et al *Nucl. Fusion*, **49**, 055006, (2009).

[143] Oikawa, T. et al *Nucl. Fusion*, **40**, 435, (2000).

[144] Van Zeeland, M. A. et al *Nucl. Fusion*, **56**, 112007, (2016).

[145] Sharapov, S. E. et al *Plasma Phys. Control Fusion*, **60**, 014026, (2018).
[146] Stix, T. H. *Physics of Fluids*, **1**(4), 308, (1958).

[147] Du, X. D. et al *Phys. Rev. Lett.*, **118**, 125001, (2017).

[148] Pace, D. C. et al *Phys. Plasmas*, **25**, 056109, (2018).

[149] Collins, C. et al *Nucl. Fusion*, **57**, 086005, (2017).

[150] Van Zeeland, M. A. et al *Phys. Plasmas*, **18**, 056114, (2011).

[151] Fenstermacher, M. E. et al *Nucl. Fusion*, **62**, 042024, (2022).

[152] Osakabe, M. et al *Nucl. Fusion*, **62**, 042019, (2022).

[153] Varela, J. et al *Nucl. Fusion*, **60**, 112015, (2020).

[154] Varela, J. et al *Nucl. Fusion*, **60**, 026016, (2020).

[155] Varela, J. et al *Phys. Plasmas*, accepted, (2024).

[156] Ida, K. et al *Rev. Sci. Instrum.*, **76**, 053505, (2005).

[157] Ida, K. et al *Fusion Sci. Technol.*, **58**, 383, (2010).

[158] Varela, J. et al *Nucl. Fusion*, **60**, 046013, (2020).

[159] Chen, Y. et al *Plasma Phys.*, **17**, 102504, (2010).

[160] Pinches, S. D. et al *Phys. Plasmas*, **22**, 021807, (2015).

[161] Rodrigues, P. et al *Nucl. Fusion*, **56**, 112006, (2016).

[162] Figueiredo, A. C. A. et al *Nucl. Fusion*, **56**, 076007, (2016).

[163] Fu, G. Y. et al *Phys. Rev. Lett.*, **75**, 2336, (1995).

[164] Fu, G. Y. et al *Plasma Phys.*, **3**, 4036, (1996).

[165] Nazikian, R. et al *Phys. Rev. Lett*, **78**, 2976, (1997).

[166] Fitzgerald, M. et al *Nucl. Fusion*, **62**, 106001, (2022).

[167] Fitzgerald, M. et al *Nucl. Fusion*, **63**, 112006, (2023).

[168] Porcelli, F. et al *Phys. Plasmas*, **1**, 470, (1994).

[169] Varela, J. et al *Plasma Phys.*, **26**, 062502, (2019).

[170] Varela, J. et al *Nucl. Fusion*, **62**, 126020, (2022).

[171] Varela, J. et al *Nucl. Fusion*, **61**, 126016, (2021).

[172] Ichiguchi, K. et al *Nucl. Fusion*, **29**, 2093, (1989).

[173] Ichiguchi, K. et al *Nucl. Fusion*, **41**, 181, (2001).

[174] Ichiguchi, K. et al *Nucl. Fusion*, **43**, 1101, (2003).

[175] Takemura, Y. et al *Nucl. Fusion*, **61**, 026011, (2021).

[176] Sakakibara, S. et al *Plasmas Phys. and Control. Fusion*, **44**, A217, (2002).

[177] Harris, J. H. et al *Phys. Rev. Lett.*, **53**, 2242, (1984).

[178] Carreras, B. A. et al *Plasma Phys.*, **5**, 3700, (1998).

[179] Varela, J. et al *Plasma Phys.*, **19**, 082501, (2012).

[180] Varela, J. et al *Plasma Phys.*, **21**, 092505, (2014).

[181] Varela, J. et al *Plasma Phys.*, **19**, 082512, (2012).

[182] Varela, J. et al *Plasma Phys.*, **21**, 032501, (2014).

[183] Chapman, I. T. et al *Nucl. Fusion*, **49**, 035006, (2009).

[184] Chao, Y. et al *AIP Advances*, **9**, 015226, (2019).

[185] Liu, D. et al *Nucl. Fusion*, **62**, 112009, (2022).

[186] Sato, M. et al *Nucl. Fusion*, **59**, 094003, (2019).

[187] Sato, M. et al *J. Plasma Phys.*, **86**, 815860305, (2020).

[188] Pinon, J. et al *J. Plasma Phys. Control. Fusion*, **60**, 075007, (2018).

[189] Rosenbluth, M. N. et al *Phys. Rev. Lett.*, **51**, 1967, (1983).

[190] White, R. B. et al *Phys. Rev. Lett.*, **62**, 539, (1989).

[191] Porcelli, F. et al *Plasma Phys. Control. Fusion*, **33**, 1601, (1991).

[192] Nakajima, N. et al *J. Plasma Fusion Res. Ser.*, **6**, 45, (2004).

[193] Miura, H. et al *Nucl. Fusion*, **50**, 054006, (2010).

[194] H. Yamada et al *Contrib. Plasma Phys.*, **50**, 480, (2010).

[195] S. Ohdachi et al *28th IAEA Fusion Energy Con. (FEC 2020)*, *Suppression of the energetic particle driven interchange mode in the large helical device*, Nice 10 – 15 May 2021, (2020) (https://iaea.org/events/fec-2020).