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

HAL
open science

# Island stability in phase space

Alexandra V Dudkovskaia, Xavier Garbet, Maxime Lesur, Howard R Wilson

► To cite this version:

Alexandra V Dudkovskaia, Xavier Garbet, Maxime Lesur, Howard R Wilson. Island stability in phase space. Journal of Physics: Conference Series, 2018, 1125, pp.012009. 10.1088/1742-6596/1125/1/012009. [https://doi.org/10.1088/1742-6596/1125/1/012009]

HAL Id: hal-02489986

https://hal.science/hal-02489986v2

Submitted on 4 Sep 2023

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_2.png" /></figure>

HAL
open science

# Island Stability in Phase Space

Alexandra V Dudkovskaia, Xavier Garbet, Maxime Lesur, Howard R Wilson

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

► To cite this version:

Alexandra V Dudkovskaia, Xavier Garbet, Maxime Lesur, Howard R Wilson. Island Stability in Phase Space. Journal of Physics: Conference Series, 2018, 1125, pp.012009. 10.1088/1742-6596/1125/1/012009. hal-02489986

HAL Id: hal-02489986

https://hal.science/hal-02489986

Submitted on 24 Feb 2020

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.
PAPER • OPEN ACCESS

# Island Stability in Phase Space

To cite this article: A V Dudkovskaja et al 2018 J. Phys.: Conf. Ser. **1125** 012009

View the [article online](#) for updates and enhancements.
# Island Stability in Phase Space

A V Dudkovskaia¹, X Garbet², M Lesur³, H R Wilson¹,⁴

¹ York Plasma Institute, Dept. of Physics, University of York, Heslington, York YO10 5DD, UK

² CEA, IRFM, F13108 St. Paul-lez-Durance cedex, France

³ IJL, UMR 7198 CNRS, University de Lorraine, Nancy, France

⁴ CCFE, Culham Science Centre, Abingdon Oxon OX14 3DB, UK

E-mail: avd512@york.ac.uk

**Abstract.** Starting with a conventional bump on tail problem, which is equivalent to finding a solution of the Vlasov/Fokker-Planck equation in the presence of the phase space island, we obtain a primary equilibrium state. The stability of this state is investigated as a function of the effective velocity-space drag and diffusion, as well as the width of these phase space islands. The secondary instabilities have been found in a certain range of plasma parameters and wave numbers. Solving the full Vlasov/Fokker-Planck – Poisson system, we obtain the dispersion function, which provides information about the secondary mode onset and allows an estimate of the secondary mode growth rate for different input plasma parameters.

## 1. Introduction

Energetic particles (EPs), generated by Resonance Heating or Neutral Beam Injection, as well as fusion alpha particles can excite Alfvén modes, resonating with plasma waves in a tokamak [1]. These instabilities, in turn, might cause EP losses [2] and degrade confinement [3]. Since fusion produced alpha particles are expected to be the main heating source in ITER and future power plants, it is crucial to predict the particle-wave interaction and control/mitigate its consequences.

In the simplest case, this problem is referred to as the bump on tail problem. The main thermal distribution provides a Maxwellian background, and the EP fraction appears as an additional peak in the tail of this Maxwellian, centred about the beam velocity, $V_b$, in velocity space [4]. The EP distribution function has a positive slope close to this beam velocity, which makes the mode unstable, provided $V_b$ is sufficiently large. Originally, this description had been applied to Langmuir waves [5], i.e. a purely electrostatic problem, but is also valid for toroidal Alfvén modes in magnetically confined plasmas [6,7]. It is well known that in the absence of sinks and sources on the right hand side of the kinetic equation, the distribution function has a "plateau" in the vicinity of the particle-wave resonance. Adding collisions always drives the particle distribution towards a Maxwellian, and hence its final shape is formed by a competition of these two processes.

To address both, purely electrostatic and magnetic cases, we apply the Hamiltonian formalism, introducing the Hamiltonian of the considered system as $H_0 = p^2/2 - \omega_b^2 \cos\xi$ with {$p, \xi$} being a pair of the generalised momentum and coordinate and $\omega_b$ denoting the bounce frequency of deeply trapped particles. Contours of constant $H_0$ in the ($p, \xi$) plane represent an island in phase space. $H_0 \ge \omega_b^2$ corresponds to the region outside the phase space island, i.e. to the passing particles, while $-\omega_b^2 \le H_0 \le \omega_b^2$ describes the contours inside the island

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

Content from this work may be used under the terms of the Creative Commons Attribution 3.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.
region, i.e. to represent particles, trapped in phase space (see Figure 1). When collisions are neglected, the EP distribution function is flattened inside the island, which decreases the drive for secondary modes. However, the distribution function gradient is steep in the vicinity of the island separatrix, which enhances the drive.

<figure><img src="image_5.png" /><figcaption>**Figure 1.** Sketch of $H_0$ vs. $\xi$ at $p = 0$. $\xi$ varies from $-\pi$ to $\pi$ for passing particles and between the bounce points, $\xi_{b1,2}$, given by $H_0 = -\omega_b^2 \cos \xi_{b1,2}$ for trapped particles.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 2.** The primary equilibrium distribution function for the case of pure diffusion (triangle down markers) and diffusion vs. drag (triangle up and square markers). Dashed line represents the initial equilibrium, i.e. in the absence of the island. $\hat{\omega}_b = 1$.</figcaption></figure>

Adding the collision operator to the problem will modify the EP distribution. In this work we focus on a joint consideration of the pitch angle scattering, velocity diffusion and slowing down effects on the form of the EP distribution function. Working in the wave reference frame, we first seek the time independent solution of the initial Fokker-Planck equation, localised to the island vicinity. Once this distribution function is found, we address the Vlasov/Fokker-Planck – Poisson system to explore the stability of this new obtained equilibrium.

We have to highlight that there exist several asymptotic solutions of the Vlasov/Fokker-Planck equation we consider for the integrable Hamiltonian system. One of them [8] was found in the framework of the Landau damping for plasmas in the full absence of collisions. It is valid in nonlinear regimes, where particle trapping is important, and includes a full calculation of the distribution function and electric field time evolution. Later it was extended to the bump on tail problem, where it predicts a saturated nonlinear state [9]. The difference is in the distribution function behaviour close to the island separatrix. The O'Neil solution assumes that as the distribution function is flat within the island due to phase mixing, it then joins smoothly the unperturbed distribution function in the island vicinity, again using the phase mixing ansatz. In contrast, Berk's prescription corresponds to the discontinuous distribution function, which is flattened inside the island, and joins the unperturbed distribution exactly at the separatrix. Other types of the functional behaviour close to the island edge are also allowed and, in general, depend on a contribution of barely passing particles, i.e. particles in the vicinity of the separatrix from the outer island region. Another solution is the Zakharov and Karpman solution [10], found in the presence of strong velocity diffusion also in steady state, and then exploited by Berk et al in [11]. The other branch of solutions is less dissipative than the Zakharov-Karpman limit and is characterised by a formation of a clump-hole pair. First, it had been found by Berk and Breizman [12,13] and was further investigated by Lilley and Nyqvist
[14], and Eriksson [15]. To consider the transition from the Zakharov-Karpman case to the Berk-Breizman solution, it is convenient to Taylor expand the initial equilibrium distribution function (i.e. distribution function in the absence of the island) in the vicinity of the resonant surface, i.e. $f_{eq} = f_{eq}|_{res} + \frac{\partial f_{eq}}{\partial p}|_{res} p$. Therefore, the Berk and Breizman solution is expected to prevail when $D_p/\omega_b^3 \sim v_{f,p}/\omega_b^2 p/\omega_b^2$ (see Eqs. (1,2) in Section 2 for more detail). Here $D_p$ and $v_{f,p}$ are the velocity diffusion and drag rates in $p$ space, and a dot denotes the partial time derivative). After some algebra, an equivalent condition can be obtained: $D_p/\gamma_L^3$ ($\gamma_L/\delta\omega$) ($\gamma_d/\gamma_L$), where $\delta\omega$ is the frequency shift in $p$ space that corresponds to a position of the resonant surface. $\gamma_L$ and $\gamma_d$ are the linear (in the absence of dissipation processes) and the ad-hoc damping rates, defined as shown in [13]. In addition, we note that the detachment condition imposes $\omega_b' \gamma_L \delta\omega$, i.e. which requires the system to be weakly unstable. The above condition on $D_p$ is roughly in agreement with the phase diagram, found by Lesur and Idomura [16], and the pioneering work by Vann et al [17], in the sense that clump/hole pairs are always found for a weak dissipation in the Vlasov/Fokker-Planck equation.

Our solution for the distribution function, presented in Section 2, is dissipative and is obtained in the wave frame. In this sense, it can be approached as an extension of the Zakharov and Karpman solution, but includes a detailed treatment of the separatrix layer, which is crucial for the secondary mode analysis. Investigating the stability of a single island in phase space, we find the primary EP distribution function, localised to the island region. From the physical point of view, we keep only one primary mode, as for instance only one toroidal Alfvén mode would appear first in the context of EP driven modes.

## 2. Plasma Response to Phase Space Island

We consider a plasma of three species: main electrons, main ions and the EP component. Assuming a Maxwellian background, we start with a calculation of the EP distribution, $f_{0,j}$, associated with the presence of the island in phase space ($j$ denotes the fast particle species: energetic electrons/ions). This distribution function corresponds to the new primary equilibrium state, whose stability we will explore in the next section. It describes the EP population of lower density, compared to the bulk plasma, and hence is to be derived as a solution of the Fokker-Planck equation:

$$ \frac{\partial f_{0,j}}{\partial t} - \{H_0, f_{0,j}\} = C(f_{0,j} - f_{eq,j}) + S. \quad (1) $$

Here curly brackets represent the conventional Poisson's bracket. $f_{eq,j}$ is the initial equilibrium distribution in the absence of the wave, shown in Figure 2 by the dashed line. Solving the kinetic equation of a form, given by Eq. (1), allows us to find $f_{0,j}$ in the Langmuir wave problem, as well as the EP distribution, associated with the Alfvén modes in tokamak plasmas. Working in the slab geometry, we introduce a single wave of the form $\Phi(x, t) = \Phi_0 \cos(k_0x - \omega_0 t)$ for the primary mode, where $\Phi$ is the electrostatic potential, and $k_0$ and $\omega_0$ are the wave number and frequency, respectively. The generalised coordinate and momentum are $\xi = k_0x - \omega_0 t$ and $p = k_0V - \omega_0 t$ with $x, V$ and $t$ being position, velocity and time. We note that $p$ is related to the particle velocity in the frame of reference of the primary wave, if $V$ here is the particle velocity in the laboratory frame. The energy in the laboratory frame, $H = \frac{m_j v^2}{2} + eZ \Phi$, is not an invariant of motion, as $\Phi$ is time dependent. On the other hand, the primary Hamiltonian function, $H_0$, introduced above, is a motion invariant, which can be easily verified. The bounce frequency of deeply trapped particles is then $\omega_b^2 = k_0^2 \frac{eZ}{2} \Phi_0$, where $eZ$ and $m_j$ are the particle charge and mass. Provided $\rho_{\theta e, i} r$, for a burning tokamak plasma the pair $\{p, \xi\}$ can be interpreted as $V/qR$ ($q$ is the safety factor and $R$ is the major radius of a tokamak, see [18] for more detail) and $n_0\phi - m_0\theta - \omega_0 t$ with $m_0/n_0$ being the poloidal/toroidal primary mode number, and $\phi$ and $\theta$ denoting the toroidal and poloidal angles, respectively. $\rho_{\theta e, i}$ is the electron/ion
poloidal Larmor radius and *r* is the minor radius of a tokamak. The small inverse aspect ratio, ε 1, circular cross section, toroidal geometry [19] has been assumed here. On the right hand side of Eq. (1), C and S are the collision operator and the source term. We take their sum as shown in Eq. (2.74) of [20]. It includes the following operators in velocity space: pitch angle scattering, diffusion and dynamical friction. Following the Berk and Breizman approach [21] and projecting the right hand side operator of Eq. (1) on the resonant surface, we reduce its dimension from 2D to 1D to arrive at a simple combination of diffusion and drag in *p*-space:

$$C + S = D_p \frac{\partial^2}{\partial p^2} + v_{f,p} \frac{\partial}{\partial p'} \quad (2)$$

where $D_p$ and $v_{f,p}$ are the corresponding diffusion and drag rates. To solve Eqs. (1,2), we impose a linear behaviour of the distribution function far from the island, where it stays unperturbed by the energetic fraction. Solving Eq. (1) in the wave reference frame for $g_{0,j} = f_{\theta,j} f_{eq,j}$ with the right hand side, Eq. (2), we obtain the full distribution function, valid inside/outside the phase space island, as well as in the separatrix layer, where the collisional effects are comparable to the free streaming contribution, ~ $p\partial/\partial\xi$. $g_{0,j}$ is treated as $g_{0,j}(\xi, H_0(p, \xi); \sigma_p) = g_{0,j}(p, \xi)$, where $\sigma_p$ is the sign of *p*. Reconstructing $f_{0,j} = f_{eq,j} + g_{0,j}$, we obtain a "perturbed" equilibrium that fully accounts for the presence of an island in phase space. $f_{0,j}$ vs. *p* is plotted in Figure 2 for different values of $D_p$ and $v_{f,p}$ (the dotted vertical lines represent the island separatrix). Hats are

used to denote the normalised quantities: $\hat{g}_{0,j}, \hat{f}_{0,j} = g_{0,j}, f_{0,j} \cdot \frac{\partial f_{eq}}{\partial p} \cdot \text{res}^{-1}$, $\hat{p} = p/(\gamma_L - \gamma_d)$, $\hat{H}_0 = H_0/(\gamma_L - \gamma_d)^2$, $\hat{D}_p = D_p/(\gamma_L - \gamma_d)^3$, $\hat{v}_{f,p} = v_{f,p}/(\gamma_L - \gamma_d)^2$ and $\hat{\omega}_b = \omega_b/(\gamma_L - \gamma_d)$. Here $\gamma_L$ is the EP contribution to the growth rate of the wave, while $\gamma_d$ is the wave damping rate due to dissipation processes. The described normalisation has been chosen to establish a connection with [22]. In a pure diffusion case, the EP distribution function remains flattened in the island region for any chosen value of $D_p$. Adding the drag term modifies its form significantly, creating a hole inside the island (see Figure 2). A similar destabilising dynamical friction effect was shown by Lilley in [22] in a pure electrostatic case in the slab geometry. Estimations, made in [22], demonstrate that the slowing down effects may be dominant over the collisional diffusion close to the resonance region. As we will see in Section 3, the form of the secondary mode dispersion function depends significantly on the inclusion of drag. $f_{0,j}$, we have found, describes the new primary equilibrium state, the stability of which will be the subject of next section.

### 3. Seeking Secondary Instability

To produce the stability analysis of the obtained equilibrium, we address the Vlasov/Fokker-Planck – Poisson system, which we write as

$$-i\delta\omega f_{j\omega} - \{H_0, f_{j\omega}\} = \{h_\omega, f_{0,j}\} \quad (3)$$

together with

$$\frac{\partial L}{\partial h_{\omega}^{*}} = 0. \quad (4)$$

Here we have assumed secondary waves of the form $\Phi(x, t) = \Phi_{k\omega}e^{ikx - i\omega t} + \text{c.c.}$ in the laboratory frame, which becomes $\Phi(\xi, t) = \Phi_{k\omega}e^{i\xi - i\delta\omega t} + \text{c.c.}$ in the primary wave frame, $l = k/k_0$ and $\delta\omega = \omega - l\omega_0$ ($k_0$ and $\omega_0$ are the primary wave number and frequency). For the full Hamiltonian and the full EP distribution function we impose $H(\xi, p) = H_0(\xi, p) + \delta H$ and $f_j(\xi, p) = f_{0,j}(\xi, p) + \delta f_j$, respectively. $H_0$ and $f_{0,j}$ are "perturbed" equilibrium quantities, defined in the previous section, while $\delta H$ and $\delta f_j$ are secondary perturbations, associated with the secondary modes and written as $\delta H = h_\omega(\xi, p)e^{-i\delta\omega t} + \text{c.c.}$ and $\delta f_j = f_{j\omega}(\xi, p)e^{-i\delta\omega t} + \text{c.c.}$.
$h_\omega$ is defined as $eZ_j\Phi_\omega$. These definitions are valid for Langmuir waves and also can be applied to consider toroidal Alfvén modes. Note: to extend the above slab formulation to a tokamak case, we have to replace $kx$ above by $n\phi - m\theta$.

The Lagrangian $L$ in Eq. (4) is given by

$$L(\omega, l) = -l^2 |h_{k\omega}|^2 + \sum_j L_j(\omega) \qquad (5)$$

with

$$L_j(\omega) = \omega_{pj}^2 \int_{-\pi}^{\pi} \frac{d\xi}{2\pi} \int_{\mathbb{Z}} f_{j\omega}(\xi, p) h_{\omega}^{\dagger}(\xi, p) dp. \qquad (6)$$

Poisson's equation requires the extremum of the functional $L$ for any given $\Phi_\omega^*$, which in the bump on tail problem is equivalent to the extremum in $h_\omega^*$. The first term of Eq. (5) is the field contribution, while the rest comes from each particle species, i.e. main ions, main electrons and the EP component, whose distribution function we have determined in Section 2. $\omega_{pj}$ is the plasma frequency of each species, and $f_j$ is normalised to the density of the considered species.

To solve Eq. (3), we split $f_{j\omega}$ into the adiabatic contribution $\frac{\partial f_{j\omega}}{\partial H_0} h_\omega$ and the resonant term $g_{j\omega}$ (note that $g_{j\omega}$ here differs from $g_{0,j}$, introduced in Section 2):

$$f_{j\omega} = \frac{\partial f_{0,j}}{\partial H_0} h_{k\omega} e^{il\xi} + g_{j\omega}. \qquad (7)$$

Thus, $g_{j\omega}$ needs to satisfy

$$-i\delta\omega g_{j\omega} + p \frac{\partial g_{j\omega}}{\partial \xi} = i\delta\omega \frac{\partial f_{0,j}}{\partial H_0} h_{k\omega} e^{il\xi}. \qquad (8)$$

Its analytic solution is

$$g_{\bar{j}\omega} = i\delta\omega \frac{\partial f_{0,j}}{\partial H_0} h_{k\omega} e^{i\delta\omega Q} \int_{\xi_{b1}}^{\xi} \frac{d\xi^j}{p^j} e^{i(\xi - \delta\omega Q^j)} + C(\sigma_p), \qquad (9)$$

where $p^j$ and $Q^j$ denote $Q^j = Q(\xi^j, H_0; \sigma_p)$ and $p^j = (\xi^j, H_0; \sigma_p)$, respectively. $Q$ by definition is

$$Q(\xi, H_0; \sigma_p) = \int_0^\xi \frac{d\xi^j}{p^j(\xi^j, H_0; \sigma_p)}, \qquad (10)$$

which has an equivalent representation through the incomplete elliptic integral of the first kind,

$$2\sigma \int_0^\xi \frac{d\xi^j}{p^j} + \omega_{pj}^2 \int_{H_0+\omega_{pj}^2}^{H_0} \frac{d\xi^j}{p^j} . C(\sigma)$$

$C(\sigma)$ is a constant of integration and is to be determined for both, passing and trapped branches. For passing particles, we simply require the distribution function to be the same at $\xi = \xi_{b1}$ and $\xi = \xi_{b2}$ for each sign of $p$. While for trapped particles, we need to ensure that $g_{j\omega}$ matches at $\xi = \xi_{b2}$ after half a bounce on $[\xi_{b1}, \xi_{b2}]$ and at $\xi = \xi_{b1}$ on the way back, i.e. $[\xi_{b2}, \xi_{b1}]$. $\sigma_p$ is different in each case. Thus, to simplify the algebra below, we introduce an angle variable, $\alpha$ [23]:

$$\alpha = \Omega_b \int_0^\xi \frac{d\xi^j}{p^j}, \quad p > 0 \qquad (11)$$

and

$$\alpha = \pi + \Omega_b \int_0^\xi \frac{d\xi^j}{p^j}, \quad p < 0 \qquad (12)$$
with $\Omega_b(H) = \int_{\xi_{b,1}}^{\xi_{b,2}} \frac{d\xi}{\pi |p|}^{-1}$. $\Omega$ is the bounce frequency, $\omega_b$ is its limit at the deeply trapped end. For passing particles, $\alpha$ is given by

$$ \alpha = \Omega_b \int_0^\xi \frac{d\xi'}{p^\jmath} \qquad (13) $$

with $\Omega_b(H_0) = \sigma_p \int_{-\pi}^{\pi} \frac{d\xi}{2\pi|p|}^{-1}$. Rewriting Eq. (9) in terms of $\alpha$ and imposing the continuity condition as $g_{j\omega}(H_0, \alpha = -\pi/2) = g_{j\omega}(H_0, \alpha = 3\pi/2)$ [23], we obtain

$$ C(\sigma_p) = \frac{-\pi \Omega_b \frac{\delta\omega}{e^{2\pi i \Omega_b} - 1}}{\int_{\pi}^{\pi} d\alpha e^{i \frac{\alpha - \frac{\delta\omega}{\Omega_b}}{\Omega_b}} \alpha} \qquad (14) $$

Substituting the obtained solution for the perturbed distribution function, Eq. (7) with Eqs. (9,14), into Eq. (5) provides

$$ D(\delta\omega, l) = -l^2 + l^2 \sum_{j=e,i} \left( \frac{\omega_{pj}^2}{\omega^2} + 2i \frac{\omega\gamma}{\omega^2 l} + \omega_{pj}^2 \right) \sum_{\sigma_p} \int_{-\infty}^{+\infty} \frac{dH}{0} \frac{\partial f}{\partial H_0} + \sum_{\sigma_p} \int_{-\infty}^{+\infty} \frac{dH_0}{\Omega_b} \frac{\partial f_{0,j}}{\partial H_0} + 2\pi i \omega^2 \sum_{\sigma_p} \int_{-\infty}^{+\infty} \frac{dH_0}{\Omega_b} \frac{\delta\omega}{\partial H_0} \int_{-\pi}^{\pi} \frac{d\alpha}{\omega \Omega_b^\alpha} C(\sigma_p) + \int_{-\infty}^{+\infty} \frac{d\alpha^j}{\Omega_b} e^{i \frac{\alpha^j - \frac{\delta\omega}{\Omega_b}}{\Omega_b}} \cdot \Theta(\sigma - \alpha^j), \qquad (15) $$

$D$ is the dispersion function, defined as $L(\delta\omega, l) = D(\delta\omega, l) h_{\delta\omega}^2$. $\Phi$ denotes the Heaviside step function. The second term of Eq. (15) is introduced for the thermal particle background.

$\gamma_j$ has been defined as $\gamma_j = \frac{1}{2} \frac{\pi}{2} \frac{1}{2} \omega_{pj} \frac{\omega_{pj}^3}{\omega^2 l_j} e^{-\frac{\omega^2 l_j}{2\omega}}$ with $\omega_{pj} = kV_{Tj}$ being the transit frequency, $V_{Tj} = (T_j/m_j)^{1/2}$ is the thermal velocity. This is a common expression of the Landau damping rate of Langmuir waves in Maxwellian plasmas. The last two contributions come from the adiabatic and resonant parts of the EP distribution function. It can be shown that $D(\delta\omega, l)$ has an equivalent form with an explicit resonance [23]:

$$ D(\delta\omega, l) = -l^2 + l^2 \sum_{j=e,i} \left( \frac{\omega_{pj}^2}{\omega^2} + 2i \frac{\omega\gamma_j}{\omega_{pj}^2} - \sum_{\sigma_p} \int_{-\infty}^{+\infty} \frac{dH_0}{n\Omega_b} \frac{\partial f_{0,j}}{\partial H_0} \cdot \frac{1}{2} \right) + \sum_{\sigma_p} \int_{-\infty}^{+\infty} \frac{dH_0}{\Omega_b} \frac{\delta\omega}{\partial H_0} \cdot \frac{1}{2} \qquad (16) $$

with

$$ \bar{h}_{n\omega} = \int_{-\pi}^{\pi} \frac{d\alpha}{2\pi} e^{i(\frac{\alpha}{\Omega_b} - n\alpha)}, \qquad (17) $$

that can be obtained by a direct discrete Fourier method in $\alpha$ space, applied to Eq. (8) [23]. $\delta\omega$ is complex here and takes the form $\delta\omega + i\gamma$. $\gamma$ corresponds to the secondary mode growth/decay rate. Eq. (15) and Eqs. (16,17) describe a full dispersion function of the secondary mode in non-resonant and resonant forms, respectively. $D(\delta\omega, \gamma) = 0$ provides the secondary mode dispersion relation. To analyse the stability of this mode, we plot constant $|D(\delta\omega, \gamma)|$ contours in the $(\delta\omega, \gamma)$ plane. Any root of $|D(\delta\omega, \gamma)|$ will appear as a pole of $|D(\delta\omega, \gamma)|^{-1}$. For simplicity, we
focus on the fast electron component, dropping the background ion contribution in the dispersion function, as $\omega_{pi}$ and $\omega_{pe}$, provided the plasma quasineutrality is maintained. The fraction of EPs is small by assumption. In Figure 3 we plot the secondary mode growth/decay rate against $\hat{D}_p$, $\hat{v}_{f,p}$ and $l$ (for convenience $\delta\omega$ and $\gamma$ are normalised to $\omega_{pe}$, $\delta\hat{\omega} = \delta\omega/\omega_{pe}$ and $\hat{\gamma} = \gamma/\omega_{pe}$, respectively). As can be seen from the figure, $\hat{\gamma}$ varies monotonically with velocity diffusion and dynamical friction rates for given $\hat{\omega}_b$ and $l$. Keeping $\hat{\omega}_b$ fixed and fixing the velocity diffusion and slowing down rates, we vary $l$ in the secondary mode dispersion relation to obtain the $l$ dependence of $\hat{\gamma}$, shown in Figure 3b. $\gamma$ as a function of $l$ is not monotonic and has two roots, $l_{crit}$ and $l_{sat}$. Between these roots $\hat{\gamma}$ is positive, which corresponds to the unstable region of secondary modes. The critical $l$ value, below which modes are stable, $l_{crit}$, against the slowing down rate is plotted in Figure 3c.

<figure><img src="image_7.png" /><figcaption>**Figure 3.** The normalised secondary mode growth/decay rate, $\hat{\gamma}$, as a function of (a) the normalised velocity diffusion growth rate, $\hat{D}_p$, (red square markers) at fixed $\hat{v}_{f,p} = 50$. and $l = 0.35$; the normalised dynamical friction rate, $\hat{v}_{f,p}$, (green round markers) at fixed $\hat{D}_p = 50$. and $l = 0.35$, and (b) as a function of $\ln l$. The bounce frequency at deeply trapped end is fixed, $\hat{\omega}_b = 1.2$. Arrows indicate roots of $\hat{\gamma} = \hat{\gamma}(l)$. $l_{crit}$ and $l_{sat}$ are the first and the second roots, respectively. $l_{crit}$ is a critical $l$ value, below which secondary modes are stable; $l_{sat}$ corresponds to the saturation level, above which $\hat{\gamma}$ becomes negative again. Solid lines are used to approximate the dependencies in the regions just before and after the main maximum. (c) $l_{crit}$ vs. $\hat{v}_{f,p}$.</figcaption></figure>

Provided $\omega_0/k_0$ is the primary island resonant velocity and $\omega/k$ is the secondary mode resonant velocity, we can estimate the $l$ value that corresponds to the maximum growth rate of the secondary modes from $\omega/k \approx \omega_0/k_0 \pm \omega_b/k_0$. Here we expect to obtain the maximum growth rate when the secondary mode resonant velocity is close to the boundary of the primary island $(\pm 2\omega_b/k_0)$ / the primary island O-point $(\pm \omega_b/k_0)$. The first case is explained by steepening of the distribution function in the vicinity of the island separatrix due to its flattening inside the island in a pure diffusion model. The second case is due to a hole in the perturbed equilibrium when drag is added. As $\omega \approx \omega_0 \approx \omega_{pe}$ typically, the latter condition roughly
<figure><img src="image_8.png" /><figcaption>**Figure 4.** Contours of constant $\log_{10} |D(\delta\hat{\omega}, \hat{\nu})|^{-1}$ (a.u.) in the $(\delta\hat{\omega}, \hat{\nu})$ plane for different values of the bounce frequency at deeply trapped end: (a) $\hat{\omega}_b = 0.6$, (b) $\hat{\omega}_b = 1.0$, (c) $\hat{\omega}_b = 1.2$, (d) $\hat{\omega}_b = 1.4$, (e) $\hat{\omega}_b = 1.6$, (f) $\hat{\omega}_b = 1.8$. In each case $\hat{D}_p = 50$, $\hat{v}_{f,p} = 40$. and $l = 0.1$.</figcaption></figure>

becomes $1 \pm \omega_b/\omega_{pe} \approx k_0/k = 1/l$, which, in turn, gives $l = 0.5$. The other important parameter is the width of the island in phase space, $2\hat{\omega}_b$. In Figure 4 we plot $\log_{10} |D|^{-1}$ in the $(\delta\hat{\omega}, \hat{\nu})$ plane for different widths of the phase space island, characterised by the bounce frequency of deeply trapped particles. $\hat{\nu}$ increases monotonically with $\hat{\omega}_b$, reaches maximum and then decreases. This functional behaviour allows us to introduce a marginal island width, below which $\hat{\nu} < 0$ and hence the mode is stable, as well as a saturation level, where $\hat{\nu}$ as a function of $\hat{\omega}_b$ has a second root. The analysis of this marginal stability is produced in [23] in more detail.

## 4. Conclusions

To summarise, in the absence of an island in phase space, created by trapping of EPs in a plasma wave, the equilibrium distribution function is simply Maxwellian, i.e. linear in the vicinity of $V_b$. Adding EPs creates the island and hence modifies the total particle distribution. Its form in the island region depends on the ratio of the diffusion and dynamical friction rates in $p$ space. Once this new equilibrium state is obtained, we investigate its stability by solving the Vlasov/Fokker-Planck-Poisson system. We have found secondary modes, which correspond to $\gamma > 0$ for a certain range of plasma parameters. $\gamma$ has been calculated as a function of the diffusion and the slowing down rates, as well as $l = k/k_0$, based on a full dispersion function, Eqs. (15,16,17). Two stable and one unstable region of $l$ have been found in range of the slowing down rates. This work was motivated by a paper by Lilley [14], where the secondary instabilities are referred to as a generic mechanism for triggering holes and clumps in phase space, and therefore explain mode frequency chirping in a general context of Alfvén waves, driven by EPs. The results, presented here, confirm that the secondary modes are possible, but under certain conditions on the primary mode amplitude and for a certain range of wave numbers of these secondary instabilities. Although the secondary mode onset is demonstrated, it still does not explain the erosion of the island separatrix. Berk and Breizman suggested the asymptotic dynamics. However, the actual physics requires more rigorous calculations and is far from being clear. We emphasise that the analysis, produced here, can be applied to both, Langmuir waves and toroidal Alfvén modes as the primary wave, and depends on a definition of the Hamiltonian function. Although we do not specify a class of the secondary modes in this work, we admit that, in principle, the above results may be applied to the EP-driven geodesic acoustic modes (EGAMs).

Investigating the stability of a single island in phase space, we took as an assumption the
existence of only one primary mode. It is worth mentioning that in a general situation, a number of harmonics can be resonant. If they resonate on the same resonant surface, the island configuration is maintained but will be deformed. If several harmonics are resonant on different surfaces, several islands will be generated and can overlap. In the region of their overlapping, stochasticity will arise, which potentially could flatten the distribution function everywhere between the two adjacent resonant surfaces, and hence prevent secondary modes. Some theories suggest that this can be relevant to the situation in ITER advanced scenarios, where several toroidal Alfvén eigenmodes (TAEs) appear simultaneously. We leave this situation for future investigation.

In this work we have focused on the role of the dynamical friction in formation of the phase space island and the secondary mode growth/decay rate behaviour. The next step will be to estimate the threshold phase island width, above which modes saturate, and hence to determine the region of the marginal stability of secondary modes.

## Acknowledgments

The authors would like to acknowledge the 9th Festival de Théorie in Aix-en-Provence, France, where this work was initiated.

## References

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

[2] S E Sharapov et al, Energetic Particle Instabilities in Fusion Plasmas, 24th IAEA Fusion Energy Conference 2012, San Diego, USA.

[3] M Ishikawa, *Nucl. Fusion* 47 (2007) 849.

[4] N J Balmforth, *Commun. Nonlinear Sci. Numer. Simulat.* 17 (2012) 1989.

[5] V E Zakharov and V I Karpman, *Sov. Phys. JETP* 16 (1963) 351.

[6] H L Berk and B N Breizman, *Phys. Fluids B* 2 (1990) 2226.

[7] M Lesur et al, *Phys. Plasmas* 17 (2012) 122311.

[8] T O'Neil, *Phys. Fluids* 8 (1965) 2255.

[9] H L Berk et al, *Phys. Plasmas* 3 (1996) 1827.

[10] V E Zakharov and V I Karpman, *Zh. Eksp. Teor Fiz.* 43 (1962) 490, and *Sov. Phys. JETP* 16 (1963) 351.

[11] H L Berk, B N Breizman and H Ye, *Phys Rev Lett* 68 (1992) 3563.

[12] H L Berk, B N Breizman and N V Petviashvili, *Phys. Lett. A* 234 (1997) 213.

[13] H L Berk et al, *Phys. Plasmas* 6 (1999) 3102.

[14] M K Lilley and R M Nyqvist, *Phys Rev Lett* 112 (2014) 155002.

[15] F Eriksson, R M Nyqvist and M K Lilley, *Phys. of Plasmas* 22 (2015) 092126.

[16] M Lesur and Y Idomura, *Nucl. Fusion* 52 (2012) 094004.

[17] R G L Vann et al, *Phys. Plasmas* 10 (2003) 623.

[18] X Garbet et al, *Journal of Physics: Conference Series* 775 (2016) 012004.

[19] R L Miller et al, *Phys. of Plasmas* 5 (1998) 973.

[20] M Lesur, The Berk-Breizman Model as a Paradigm for Energetic Particle-driven Alfvén Eigenmodes, a *dissertation submitted for the Degree of Doctor of Philosophy* (2010) 117.

[21] H L Berk, B N Breizman and M S Pekker, *Plasma Phys. Rep.* 23 (1997) 778.

[22] M K Lilley, B N Breizman and S E Sharapov, *Phys. Rev. Lett.* 102 (2009) 195003.

[23] A V Dudkovskaia, X Garbet, M Lesur and H R Wilson, in preparation (2018).