**5. Global numerical analyses of MHD modes: AEGIS code formalism**

In Sec. 2.3 analytical or semi-analytical theories are presented to describe four types of MHD modes in toroidal geometry. Due to the developments of modern numerical method and computer hardware, conventional asymptotic expansion methods for global modes have become outdated and been substituted by direct numerical computation. Several excellent numerical codes have been developed in the past to study linear MHD stability of toroidally confined plasmas, such as such as PEST Grimm et al. (1976) Chance et al. (1978), GATO Bernard et al. (1981), DCON Glasser (1997), AEGIS Zheng & Kotschenreuther (2006), etc. In this section we focus on description of AEGIS code, in view of that AEGIS is an adaptive MHD shooting code capable to study MHD continuum Zheng et al. (2005). Through describing AEGIS formalism, we can further explain the general features of MHD eigen modes in toroidally confined plasmas.

Let us first describe the toroidal system to be investigated. The core part is plasma torus, which is surround by a resistive wall; Between plasma torus and resistive wall there is inter vacuum region and outside the resistive wall there is outer vacuum region, which extends to infinity. For simplicity, it is assumed that the wall is thin. We denote respectively the interfaces between plasma torus and inner vacuum region, inner vacuum region and wall, and wall and outer vacuum region as *<sup>ψ</sup>a*, *<sup>ψ</sup>b*−, and *<sup>ψ</sup>b*+.

22 Will-be-set-by-IN-TECH

*dεdμB*

Here, the superscript ∗ represents complex conjugate. Matching inner (Eq. (89)) and outer

Here, we note that kinetic effects from core plasma should also be taken into account in outer region. As proved in Ref. Zheng & Tessarotto (1994a) this results in the so-called apparent mass effect and leads Ω in the first term of Eq. (90) to become complicated function of actual

Similarly, for KDMs of TAE type, for example EPMs, one need to consider even and odd modes. For even modes the dispersion relation is given by Refs. Zheng et al. (2000) and Tsai

<sup>2</sup>

where *δTf* represents MHD fluid contribution and *δTK* is energetic-particle contribution to the

The dispersion relations in Eqs. (90) and (91) extend respectively MHD ballooning modes in diamagnetic gap and TAEs in Alfvén gap to respective continua. Kinetic drives are the causes

In Sec. 2.3 analytical or semi-analytical theories are presented to describe four types of MHD modes in toroidal geometry. Due to the developments of modern numerical method and computer hardware, conventional asymptotic expansion methods for global modes have become outdated and been substituted by direct numerical computation. Several excellent numerical codes have been developed in the past to study linear MHD stability of toroidally confined plasmas, such as such as PEST Grimm et al. (1976) Chance et al. (1978), GATO Bernard et al. (1981), DCON Glasser (1997), AEGIS Zheng & Kotschenreuther (2006), etc. In this section we focus on description of AEGIS code, in view of that AEGIS is an adaptive MHD shooting code capable to study MHD continuum Zheng et al. (2005). Through describing AEGIS formalism, we can further explain the general features of MHD eigen modes in

Let us first describe the toroidal system to be investigated. The core part is plasma torus, which is surround by a resistive wall; Between plasma torus and resistive wall there is inter vacuum region and outside the resistive wall there is outer vacuum region, which extends to infinity. For simplicity, it is assumed that the wall is thin. We denote respectively the interfaces between plasma torus and inner vacuum region, inner vacuum region and wall, and wall and

<sup>|</sup>*v*�| *<sup>ω</sup>dδgh*.

*<sup>p</sup>* <sup>−</sup> (*<sup>s</sup>* <sup>−</sup> *<sup>α</sup>* cos *<sup>θ</sup>*)<sup>2</sup> *p*2

 |*ζ*| 2 ,

+ *δTf* + *δTK* = 0, (91)

− *i*Ω + *δWf* + *δWK* = 0. (90)

where

mode frequency.

& Chen (1993)

quadratic form in inner region.

toroidally confined plasmas.

outer vacuum region as *<sup>ψ</sup>a*, *<sup>ψ</sup>b*−, and *<sup>ψ</sup>b*+.

to make causality conditions satisfied.

*δWf* =

*δWk* =

 +∞ −∞ *dθ ∂ζ ∂θ* 2 − *α* cos *θ*

 +∞ −∞

> − *i* Ω<sup>2</sup> <sup>−</sup> <sup>−</sup> <sup>Ω</sup><sup>2</sup>

Ω<sup>2</sup> <sup>+</sup> − <sup>Ω</sup><sup>2</sup>

**5. Global numerical analyses of MHD modes: AEGIS code formalism**

*<sup>d</sup>θζ*<sup>∗</sup> <sup>1</sup> *p*1/2

(Eq. (85)) solutions one obtains the dispersion relation Tsai & Chen (1993)

#### **5.1 MHD equations and numerical solution method for plasma region**

In this subsection we describe how to reduce MHD equations for global mode analyses. The starting equation is the single fluid MHD equation (14). This is a vector equation and can be projected onto three directions to get scalar equations. The parallel projection has be derived in Eq. (33). From parallel equation one can solve for ∇ · *ξ*, which is the only unknown needed for two-perpendicular equations to become a complete set of equations. In principle the parallel motion can not be described by MHD model, since particles are not localized along magnetic field line. There are wave-partcle resonance, trapped particle, and so-called small parallel particle speed effects, etc. Nevertheless, from analyses in Sec. 4.3 one can see that in low frequency regime the parallel coupling results only in the so-called apparent mass effect, while in intermediate regime the parallel coupling mainly gives rise to the 2nd TAEs. Note that apparent mass effect can be absorbed by rescaling mode frequency and inclusion of the 2nd TAE effect is straightforward as discussed in Sec. 4.4. For brevity we limit ourselves to treat only two perpendicular components of Eq. (14) with Γ set to zero. AEGIS-K code has been developed to include parallel dynamics using kinetic description Zheng et al. (2010).

Using general flux coordnates in Eq. (22), the magnetic field line displacement is decomposed as follows

$$\mathbf{f} \times \mathbf{B} = \xi\_s \nabla \psi + \xi\_\psi \chi'(\nabla \zeta - q \nabla \theta). \tag{92}$$

Since we deal with linear problem, the Fourier transform can be used to decompose perturbed quantities in poloidal and toroidal directions,

$$\left\{\tilde{\xi}\exp\{-in\zeta\}\right\}=\sum\_{m=-\infty}^{\infty}\tilde{\xi}\_{m}\frac{1}{\sqrt{2\pi}}\exp\{i\left(m\theta-n\zeta\right)\},\tag{93}$$

with *<sup>ξ</sup><sup>m</sup>* = *<sup>π</sup>* <sup>−</sup>*<sup>π</sup> <sup>d</sup>θξ* exp{−*imθ*}/ <sup>√</sup>2*π*. With the toroidal symmetry assumed, only a single toroidal Fourier component needs to be considered. As usual, equilibrium quantities can be decomposed as matrices in poloidal Fourier space, for example

$$\mathcal{J}\_{mm'} = \frac{1}{2\pi} \int\_{-\pi}^{\pi} d\theta J(\theta) e^{i(m'-m)\theta}.$$

In the poloidal Fourier decomposition, the Fourier components are cut off both from lower and upper sides respectively by *m*min and *m*max. Therefore, the total Fourier component under consideration is *M* = *m*max − *m*min + 1. We also use bold face (or alternatively [[··· ]]) to represent Fourier space vectors, and calligraphic capital letters (or alternatively �···�) to represent the corresponding equilibrium matrices (*e.g.*, J for *J*) in poloidal Fourier space.

To get scalar equations, we project Eq. (14) respectively onto two directions *<sup>J</sup>*2∇*<sup>θ</sup>* × ∇*<sup>ζ</sup>* · **<sup>B</sup>** <sup>×</sup> [··· × **<sup>B</sup>**]/*B*<sup>2</sup> and (1/*qχ*� )*J*2∇*<sup>ζ</sup>* × ∇*<sup>ψ</sup>* · **<sup>B</sup>** <sup>×</sup> [··· × **<sup>B</sup>**]/*B*2, and then introduce the Fourier transformation in Eq. (93) to two projected equations. These procedures lead to the following set of differential equations in matrices

$$\left(\mathcal{B}^{\dagger}\mathfrak{F}\_{s} + \mathcal{D}\mathfrak{F}\_{\psi}^{\prime} + \mathcal{E}\mathfrak{F}\_{\psi}\right)^{\prime} - \left(\mathcal{E}^{\dagger}\mathfrak{F}\_{s} + \mathcal{E}^{\dagger}\mathfrak{F}\_{\psi}^{\prime} + \mathcal{H}\mathfrak{F}\_{\psi}\right) = 0,\tag{94}$$

$$
\mathcal{A}\mathfrak{F}\_s + \mathcal{B}\mathfrak{F}\_\psi' + \mathcal{C}\mathfrak{F}\_\psi = 0. \tag{95}
$$

in Toroidal Plasma Confinement 25

Overview of Magnetohydrodynamics Theory in Toroidal Plasma Confinement 25

where <sup>F</sup> <sup>=</sup> D−B†A−1B, <sup>K</sup> <sup>=</sup> E−B†A−1C, and <sup>G</sup> <sup>=</sup> H−C†A−1C. These matrices can be

*<sup>N</sup>*A*<sup>i</sup>* + Q(*n*G<sup>23</sup> + G33M)

I <sup>−</sup> *<sup>i</sup>γ*<sup>2</sup> *N*C*i* 

 .

**u**� = L**u**, (99)

 <sup>A</sup>−1<sup>C</sup>

G31Q − *g*�

, where **u**<sup>2</sup> = F*ξ*� + K*ξ*, Eq. (96) is

**c***p*, (100)

, (97)

. (98)

*N*A*i* − *γ*2

*<sup>N</sup>*A*<sup>i</sup>* + (*n*G<sup>23</sup> + MG33)Q

*<sup>N</sup>*A*<sup>i</sup>* + Q(*n*G<sup>23</sup> + G33M)

G<sup>33</sup> − *iχ*�

 −F−1K F−<sup>1</sup> G−K†F−1K K†F−<sup>1</sup>

We note that *ξ* and **u**<sup>2</sup> in plasma region are related to the magnetic field and pressure as follows

− [[J (**B** · *δ***B** − *ξ* · ∇*P*)]] = **u**2. The set of eigen mode equations in Eq. (99) can be solved numerically by independent solution method together with multiple region matching technique as described in Ref. Zheng & Kotschenreuther (2006). With *M* boundary conditions imposed at magnetic axis, there remain

**u**1

where the superscripts are used to label independent solutions. We use the cylinder limit to describe the boundary condition at magnetic axis, *i.e., ξψ*,*<sup>m</sup>* ∝ *rm*. The general solution can be

where **c***<sup>p</sup>* is a constant vector with *M* elements. Without loss of generality (by defining **c***<sup>p</sup>* =

[[J ∇*ψ* · *δ***B**]] = *i*Q*ξ*,

*<sup>ξ</sup>*1, ··· , *<sup>ξ</sup><sup>M</sup>*

<sup>2</sup>, ··· , **<sup>u</sup>***<sup>M</sup>* 2

 ,

*<sup>p</sup>* ), we can set Ξ*<sup>p</sup>* to be unity I at plasma edge. Therefore, at

− [[J (**B** · *δ***B** − *ξ* · ∇*P*)]] = *i*W*p***c***p*. (102)

[[J ∇*ψ* · *δ***B**]] = −Q**c***p*, (101)

 *ξ* **u**2 

)�

further simplified as follows Glasser (1997)

<sup>F</sup> <sup>=</sup> *<sup>χ</sup>*�<sup>2</sup> *n*2 

<sup>K</sup> <sup>=</sup> *<sup>χ</sup>*� *n i γ*2

×A−<sup>1</sup> *γ*2

−Q

Introducing the expanded 2*M* unknowns **u** =

reduced to the set of 2*M* first order equations

where 2*M* × 2*M* matrix

only *M* independent solutions:

*<sup>p</sup>* and <sup>W</sup>*new*

*<sup>p</sup>* <sup>=</sup> <sup>W</sup>*p*Ξ−<sup>1</sup>

plasma-vacuum interface *ψ<sup>a</sup>* we have

Ξ−<sup>1</sup> *<sup>p</sup>* **c***new* QG33<sup>Q</sup> <sup>+</sup> *<sup>γ</sup>*<sup>2</sup>

*χ*��G<sup>23</sup> + (*qχ*�

L =

 Ξ*<sup>p</sup>* W2

then obtained as a combination of the *M* independent solutions,

 ≡ 

 *ξ* **u**2 = *i* Ξ*<sup>p</sup>* W*p*

Here, the equilibrium matrices contain two parts: plasma/field force and inertia contributions, e.g., A = A*<sup>p</sup>* + *γN*A*i*, where

A*<sup>p</sup>* = *n*(*n*G<sup>22</sup> + G23M) + M(*n*G<sup>23</sup> + G33M), B*<sup>p</sup>* = −*iχ*� [*n*(G<sup>22</sup> + *q*G23) + M(G<sup>23</sup> + *q*G33)] , C*<sup>p</sup>* = −*i χ*��(*n*G<sup>22</sup> + MG23)+(*qχ*� )� (*n*G<sup>23</sup> + MG33) −*χ*� (*n*G<sup>12</sup> + MG31)Q + *i*(*g*� Q − *μ*0*nP*� J /*χ*� ), <sup>D</sup>*<sup>p</sup>* <sup>=</sup> *<sup>χ</sup>*�<sup>2</sup> [(G<sup>22</sup> <sup>+</sup> *<sup>q</sup>*G23) + *<sup>q</sup>*(G<sup>23</sup> <sup>+</sup> *<sup>q</sup>*G33)] , <sup>E</sup>*<sup>p</sup>* <sup>=</sup> *<sup>χ</sup>*� *χ*��(G<sup>22</sup> + *q*G23)+(*qχ*� )� (G<sup>23</sup> + *q*G33) <sup>−</sup> *<sup>i</sup>χ*�2(G<sup>12</sup> <sup>+</sup> *<sup>q</sup>*G31)<sup>Q</sup> <sup>+</sup> *<sup>μ</sup>*0*P*� J , <sup>H</sup>*<sup>p</sup>* <sup>=</sup> *<sup>χ</sup>*�� *χ*��G<sup>22</sup> + (*qχ*� )� G23 + (*qχ*� )� *χ*��G<sup>23</sup> + (*qχ*� )� G33 +*iχ*� *χ*��(MG<sup>12</sup> − G12M)+(*qχ*� )� (MG<sup>31</sup> − G31M) <sup>+</sup>*χ*�2QG11<sup>Q</sup> <sup>+</sup> *<sup>μ</sup>*0*P*� *χ*��J /*χ*� + *μ*0*P*� J � − *g*� *q*� *χ*� I, <sup>A</sup>*<sup>i</sup>* <sup>=</sup> *<sup>B</sup>*<sup>2</sup> 0 *X*2 0 *q*2 0 <sup>J</sup> *<sup>ρ</sup><sup>N</sup> <sup>B</sup>*<sup>2</sup> |∇*ψ*<sup>|</sup> 2 , <sup>C</sup>*<sup>i</sup>* <sup>=</sup> *<sup>B</sup>*<sup>2</sup> 0 *X*2 0 *q*2 0 *χ*� <sup>J</sup> *<sup>ρ</sup><sup>N</sup> <sup>B</sup>*<sup>2</sup> (∇*<sup>ψ</sup>* · ∇*<sup>ζ</sup>* <sup>−</sup> *<sup>q</sup>*∇*<sup>ψ</sup>* · ∇*θ*) , <sup>H</sup>*<sup>i</sup>* <sup>=</sup> *<sup>B</sup>*<sup>2</sup> 0 *X*2 0 *q*2 0 *<sup>χ</sup>*�2<sup>J</sup> *<sup>ρ</sup><sup>N</sup> <sup>B</sup>*<sup>2</sup> (|∇*ζ*<sup>|</sup> <sup>2</sup> <sup>+</sup> *<sup>q</sup>*2|∇*θ*<sup>|</sup> <sup>2</sup> <sup>−</sup> <sup>2</sup>∇*<sup>θ</sup>* · ∇*ζ*) ,

B*<sup>i</sup>* = D*<sup>i</sup>* = E*<sup>i</sup>* = 0, M*mm*� = *m*I*mm*� , Q*mm*� = (*m* − *nq*)I*mm*� , *γ<sup>N</sup>* denotes the dimensionless growth rate normalized by the Alfvén frequency at magnetic axis, *ρ<sup>N</sup>* is the dimensionless mass density normalized by the mass density at magnetic axis, subscript "0" refers to quantities at magnetic axis, and

$$\begin{aligned} \mathcal{G}\_{11} &= \langle f(\nabla \theta \times \nabla \xi) \cdot (\nabla \theta \times \nabla \xi) \rangle, \\ \mathcal{G}\_{22} &= \langle f(\nabla \xi \times \nabla \psi) \cdot (\nabla \xi \times \nabla \psi) \rangle, \\ \mathcal{G}\_{33} &= \langle f(\nabla \psi \times \nabla \theta) \cdot (\nabla \psi \times \nabla \theta) \rangle, \\ \mathcal{G}\_{12} &= \langle f(\nabla \theta \times \nabla \xi) \cdot (\nabla \zeta \times \nabla \psi) \rangle, \\ \mathcal{G}\_{31} &= \langle f(\nabla \psi \times \nabla \theta) \cdot (\nabla \theta \times \nabla \zeta) \rangle, \\ \mathcal{G}\_{23} &= \langle f(\nabla \xi \times \nabla \psi) \cdot (\nabla \psi \times \nabla \theta) \rangle. \end{aligned}$$

We can reduce the set of equations (94) and (95) into a set of first order differential equations as in the DCON formalism Glasser (1997). By solving Eq. (95), one obtains

$$\mathfrak{f}\_s = -\mathcal{A}^{-1}\mathcal{B}\mathfrak{f}'\_{\psi} - \mathcal{A}^{-1}\mathcal{C}\mathfrak{f}\_{\psi}.$$

Inserting this solution into Eq. (94), we get

$$\frac{d}{d\psi} \left( \mathcal{F}\mathfrak{F}' + \mathcal{K}\mathfrak{F} \right) - \left( \mathcal{K}^\dagger \mathfrak{F}' + \mathcal{G}\mathfrak{F} \right) = 0,\tag{96}$$

where <sup>F</sup> <sup>=</sup> D−B†A−1B, <sup>K</sup> <sup>=</sup> E−B†A−1C, and <sup>G</sup> <sup>=</sup> H−C†A−1C. These matrices can be further simplified as follows Glasser (1997)

$$\mathcal{F} = \frac{\chi'^2}{n^2} \left\{ \mathcal{Q}\mathcal{G}\_{33}\mathcal{Q} + \gamma\_N^2 \mathcal{A}\_i - \left[ \gamma\_N^2 \mathcal{A}\_i + \mathcal{Q}(n\mathcal{G}\_{23} + \mathcal{G}\_{33}\mathcal{M}) \right] \right\}$$

$$\times \mathcal{A}^{-1} \left[ \gamma\_N^2 \mathcal{A}\_i + (n\mathcal{G}\_{23} + \mathcal{M}\mathcal{G}\_{33})\mathcal{Q} \right] \Big\}, \tag{97}$$

$$\mathcal{K} = \frac{\chi'}{n} \left\{ i \left[ \gamma\_N^2 \mathcal{A}\_i + \mathcal{Q}(n\mathcal{G}\_{23} + \mathcal{G}\_{33}\mathcal{M}) \right] \mathcal{A}^{-1} \mathcal{C} $$

$$-\mathcal{Q} \left[ \chi'' \mathcal{G}\_{23} + (q\chi')' \mathcal{G}\_{33} - i\chi' \mathcal{G}\_{31}\mathcal{Q} - g'\mathcal{T} \right] - i\gamma\_N^2 \mathcal{C}\_i \right\}. \tag{98}$$

Introducing the expanded 2*M* unknowns **u** = *ξ* **u**2 , where **u**<sup>2</sup> = F*ξ*� + K*ξ*, Eq. (96) is reduced to the set of 2*M* first order equations

$$\mathbf{u}' = \mathcal{L}\mathbf{u},\tag{99}$$

where 2*M* × 2*M* matrix

24 Will-be-set-by-IN-TECH

Here, the equilibrium matrices contain two parts: plasma/field force and inertia

)�

)�

*χ*��J /*χ*� + *μ*0*P*�

*<sup>B</sup>*<sup>2</sup> (∇*<sup>ψ</sup>* · ∇*<sup>ζ</sup>* <sup>−</sup> *<sup>q</sup>*∇*<sup>ψ</sup>* · ∇*θ*)

as in the DCON formalism Glasser (1997). By solving Eq. (95), one obtains

F*ξ*� + K*ξ*

*<sup>ξ</sup><sup>s</sup>* <sup>=</sup> −A−1B*ξ*�

 − 

<sup>2</sup> <sup>+</sup> *<sup>q</sup>*2|∇*θ*<sup>|</sup>

Q − *μ*0*nP*�

(G<sup>23</sup> + *q*G33)

)�

B*<sup>i</sup>* = D*<sup>i</sup>* = E*<sup>i</sup>* = 0, M*mm*� = *m*I*mm*� , Q*mm*� = (*m* − *nq*)I*mm*� , *γ<sup>N</sup>* denotes the dimensionless growth rate normalized by the Alfvén frequency at magnetic axis, *ρ<sup>N</sup>* is the dimensionless mass density normalized by the mass density at magnetic axis, subscript "0" refers to

G<sup>11</sup> = �*J*(∇*θ* × ∇*ζ*) · (∇*θ* × ∇*ζ*)�, G<sup>22</sup> = �*J*(∇*ζ* × ∇*ψ*) · (∇*ζ* × ∇*ψ*)�, G<sup>33</sup> = �*J*(∇*ψ* × ∇*θ*) · (∇*ψ* × ∇*θ*)�, G<sup>12</sup> = �*J*(∇*θ* × ∇*ζ*) · (∇*ζ* × ∇*ψ*)�, G<sup>31</sup> = �*J*(∇*ψ* × ∇*θ*) · (∇*θ* × ∇*ζ*)�, G<sup>23</sup> = �*J*(∇*ζ* × ∇*ψ*) · (∇*ψ* × ∇*θ*)�. We can reduce the set of equations (94) and (95) into a set of first order differential equations

)�

(*n*G<sup>23</sup> + MG33)

J /*χ*� ),

*χ*��G<sup>23</sup> + (*qχ*�

J � − *g*� *q*� *χ*� I,

> ,

<sup>2</sup> <sup>−</sup> <sup>2</sup>∇*<sup>θ</sup>* · ∇*ζ*)

*<sup>ψ</sup>* − A−1C*ξψ*.

<sup>K</sup>†*ξ*� <sup>+</sup> <sup>G</sup>*<sup>ξ</sup>*

= 0, (96)

)� G33 

 ,

(MG<sup>31</sup> − G31M)

<sup>−</sup> *<sup>i</sup>χ*�2(G<sup>12</sup> <sup>+</sup> *<sup>q</sup>*G31)<sup>Q</sup> <sup>+</sup> *<sup>μ</sup>*0*P*�

J ,

contributions, e.g., A = A*<sup>p</sup>* + *γN*A*i*, where

C*<sup>p</sup>* = −*i*

<sup>E</sup>*<sup>p</sup>* <sup>=</sup> *<sup>χ</sup>*�

<sup>H</sup>*<sup>p</sup>* <sup>=</sup> *<sup>χ</sup>*��

<sup>A</sup>*<sup>i</sup>* <sup>=</sup> *<sup>B</sup>*<sup>2</sup> 0 *X*2 0 *q*2 0 <sup>J</sup> *<sup>ρ</sup><sup>N</sup> <sup>B</sup>*<sup>2</sup> |∇*ψ*<sup>|</sup> 2 ,

<sup>C</sup>*<sup>i</sup>* <sup>=</sup> *<sup>B</sup>*<sup>2</sup> 0 *X*2 0 *q*2 0 *χ*� <sup>J</sup> *<sup>ρ</sup><sup>N</sup>*

<sup>H</sup>*<sup>i</sup>* <sup>=</sup> *<sup>B</sup>*<sup>2</sup> 0 *X*2 0 *q*2 0 

quantities at magnetic axis, and

Inserting this solution into Eq. (94), we get

*d dψ* 

−*χ*�

+*iχ*�

A*<sup>p</sup>* = *n*(*n*G<sup>22</sup> + G23M) + M(*n*G<sup>23</sup> + G33M), B*<sup>p</sup>* = −*iχ*� [*n*(G<sup>22</sup> + *q*G23) + M(G<sup>23</sup> + *q*G33)] ,

*χ*��(*n*G<sup>22</sup> + MG23)+(*qχ*�

(*n*G<sup>12</sup> + MG31)Q + *i*(*g*�

<sup>D</sup>*<sup>p</sup>* <sup>=</sup> *<sup>χ</sup>*�<sup>2</sup> [(G<sup>22</sup> <sup>+</sup> *<sup>q</sup>*G23) + *<sup>q</sup>*(G<sup>23</sup> <sup>+</sup> *<sup>q</sup>*G33)] ,

*χ*��(G<sup>22</sup> + *q*G23)+(*qχ*�

)� G23 + (*qχ*�

*χ*��(MG<sup>12</sup> − G12M)+(*qχ*�

*χ*��G<sup>22</sup> + (*qχ*�

<sup>+</sup>*χ*�2QG11<sup>Q</sup> <sup>+</sup> *<sup>μ</sup>*0*P*�

*<sup>χ</sup>*�2<sup>J</sup> *<sup>ρ</sup><sup>N</sup>*

*<sup>B</sup>*<sup>2</sup> (|∇*ζ*<sup>|</sup>

$$
\mathcal{L} = \begin{pmatrix} -\mathcal{F}^{-1}\mathcal{K} & \mathcal{F}^{-1} \\ \mathcal{G} - \mathcal{K}^{\dagger}\mathcal{F}^{-1}\mathcal{K}\ \mathcal{K}^{\dagger}\mathcal{F}^{-1} \end{pmatrix}.
$$

We note that *ξ* and **u**<sup>2</sup> in plasma region are related to the magnetic field and pressure as follows

$$\begin{aligned} [[\mathcal{I}\nabla\psi \cdot \delta \mathbf{B}]] &= i\mathcal{Q}\mathbf{f}, \\ -[[\mathcal{J}\left(\mathbf{B}\cdot\delta \mathbf{B} - \mathbf{f}\cdot\nabla P\right)]] &= \mathbf{u}\_2. \end{aligned}$$

The set of eigen mode equations in Eq. (99) can be solved numerically by independent solution method together with multiple region matching technique as described in Ref. Zheng & Kotschenreuther (2006). With *M* boundary conditions imposed at magnetic axis, there remain only *M* independent solutions:

$$
\begin{pmatrix} \Xi\_p \\ \mathcal{W}\_2 \end{pmatrix} \equiv \begin{pmatrix} \mathfrak{f}^1 \, \cdot \, \cdot \, \cdot \, \mathfrak{f}^M \\ \mathfrak{u}^1\_{2'} \, \cdot \, \cdot \, \cdot \, \mathfrak{u}^M\_2 \end{pmatrix} \prime
$$

where the superscripts are used to label independent solutions. We use the cylinder limit to describe the boundary condition at magnetic axis, *i.e., ξψ*,*<sup>m</sup>* ∝ *rm*. The general solution can be then obtained as a combination of the *M* independent solutions,

$$\begin{pmatrix} \mathfrak{F} \\ \mathfrak{u}\_2 \end{pmatrix} = i \begin{pmatrix} \Xi\_p \\ \mathcal{W}\_p \end{pmatrix} \mathfrak{c}\_{p\prime} \tag{100}$$

where **c***<sup>p</sup>* is a constant vector with *M* elements. Without loss of generality (by defining **c***<sup>p</sup>* = Ξ−<sup>1</sup> *<sup>p</sup>* **c***new <sup>p</sup>* and <sup>W</sup>*new <sup>p</sup>* <sup>=</sup> <sup>W</sup>*p*Ξ−<sup>1</sup> *<sup>p</sup>* ), we can set Ξ*<sup>p</sup>* to be unity I at plasma edge. Therefore, at plasma-vacuum interface *ψ<sup>a</sup>* we have

$$[[\mathcal{J}\nabla\psi \cdot \delta \mathbf{B}]] = -\mathcal{Q}\mathbf{c}\_{p\prime} \tag{101}$$

$$-\left[\left[\mathcal{J}\left(\mathbf{B}\cdot\delta\mathbf{B}-\mathfrak{F}\cdot\nabla P\right)\right]\right] = i\mathcal{W}\_p\mathbf{c}\_p.\tag{102}$$

in Toroidal Plasma Confinement 27

Overview of Magnetohydrodynamics Theory in Toroidal Plasma Confinement 27

In the inner vacuum region, the independent solutions can be constructed, for example, with the use of an inward numerical shooting Zheng & Kotschenreuther (2006), with the following

> = I O

= T I

where O is *M* × *M* zero matrix. Since the boundary conditions in Eq. (105) give *δ***B** · ∇*ψ* = 0 at wall, these conditions correspond to a set of solutions that corresponds to the perfectly conducting wall type. On the other hand, since the boundary conditions in Eq. (106) guarantee that the independent solutions to be continuous with outer vacuum solutions, these conditions correspond to a set of solutions that corresponds to the no-wall type. Using the general expression for the solutions in Eq. (104), we can express the normal and parallel

The solutions in the plasma and vacuum regions described in the last two subsections can be used to construct the eigen value problem Zheng & Kotschenreuther (2006). The normal magnetic field component and the combined magnetic and thermal pressures are required to be continuous at the plasma-vacuum interface. Matching plasma [Eqs. (101) and (102)] and

[[J ∇*ψ* · *δ***B**]] = −V1**c***v*<sup>1</sup> − V2**d***v*1, (107)

−[[J **B** · *δ***B**]] = *i*Q (U1**c***v*<sup>1</sup> + U2**d***v*1). (108)

<sup>1</sup> *<sup>δ</sup>*W*bδ*W−<sup>1</sup> <sup>∞</sup> <sup>F</sup>2**c***v*1, (109)

<sup>U</sup><sup>2</sup> − U1V−<sup>1</sup>

1 V2 *ψa* ,

<sup>1</sup> ]*ψa*Q, F<sup>1</sup> = Q

. Note that *δ*W<sup>∞</sup> and *δ*W*<sup>b</sup>* correspond to the energy matrices

∇×∇× *δ***B** = −*γμ*0*σδ***B**, (110)

, (105)

, (106)

 U1 V1

 U2 V2

magnetic fields at the plasma-vacuum interface as follows:

vacuum [Eqs. (107) and (108)] solutions at the interface *ψ<sup>a</sup>* gives

Maxwell equation ∇ · *δ***B** = 0 and the thin wall assumption lead to

**<sup>d</sup>***v*<sup>1</sup> <sup>=</sup> <sup>F</sup>−<sup>1</sup>

<sup>2</sup> ]*ψa*Q, *<sup>δ</sup>*W*<sup>b</sup>* <sup>=</sup> <sup>W</sup>*<sup>p</sup>* − Q[U1V−<sup>1</sup>

without a wall and with a perfectly conducting wall at *ψb*, respectively, as can be seen from

We now consider the matching across the thin resistive wall. For the radial magnetic field, the

**v**|*ψb*<sup>−</sup> = **v**|*ψb*<sup>+</sup> = **d***v*1. The current in the resistive wall causes a jump in the scalar magnetic potential. This can be

*ψb*<sup>−</sup>

*ψb*<sup>−</sup>

boundary conditions imposed at *<sup>ψ</sup>b*−:

**5.3 Eigenvalue problem**

where *<sup>δ</sup>*W<sup>∞</sup> <sup>=</sup> <sup>W</sup>*<sup>p</sup>* − Q[U2V−<sup>1</sup>

obtained from the Ampére law

<sup>U</sup><sup>1</sup> − U2V−<sup>1</sup>

2 V1 *ψa*

the boundary conditions in Eqs. (105) and (106).

and F<sup>2</sup> = Q

#### **5.2 The solution of vacuum region**

For completeness, in this subsection we briefly review the vacuum solutions in Ref. Zheng & Kotschenreuther (2006). The vacuum regions are described by the Laplace equation

$$
\nabla^2 \mu = 0,\tag{103}
$$

where *u* is the magnetic scalar potential and is related to the perturbed magnetic field by *δ***B** = −∇*u*. Here, we note that this representation of vacuum magnetic field, although being simple, excludes the consideration of *n* = 0 modes. To study *n* = 0 modes, one more scalar is needed to represent the vacuum magnetic field. For the sake of conciseness, we outline the general solutions for inner and outer vacuum regions simultaneously.

As in the plasma region, Fourier decompositions are introduced for both poloidal and toroidal directions to solve Eq. (103). Then Eq. (103) becomes a set of second-order differential equations of number *M* for **u**. This set of second-order differential equations can be transformed into a set of first-order differential equations of number 2*M*, by introducing an additional field **v** = −[[J ∇*ψ* · *δ***B**]], which is related to the magnetic scalar potential in Fourier space as follows:

$$\mathbf{v} = \left< \mathcal{J} |\nabla \psi|^2 \right> \frac{\partial \mathbf{u}}{\partial \boldsymbol{\psi}} + \left< i \mathcal{J} \nabla \psi \cdot \nabla \theta \right> \mathcal{M} \mathbf{u}.$$

There are 2*M* independent solutions for Eq. (103), which can be used to construct the following independent solution matrices:

$$
\begin{aligned}
\begin{pmatrix}
\mathcal{U}\_1 \\
\mathcal{V}\_1
\end{pmatrix} &\equiv \begin{pmatrix}
\mathbf{u}^1 & \cdots & \mathbf{u}^M \\
\mathbf{v}^1 & \cdots & \mathbf{v}^M
\end{pmatrix}, \\
\begin{pmatrix}
\mathcal{U}\_2 \\
\mathcal{V}\_2
\end{pmatrix} &\equiv \begin{pmatrix}
\mathbf{u}^{M+1} & \cdots & \mathbf{u}^{2M} \\
\mathbf{v}^{M+1} & \cdots & \mathbf{v}^{2M}
\end{pmatrix}.
\end{aligned}
$$

The general solutions in the vacuum regions can be expressed as a linear combination of the independent solutions:

$$
\begin{pmatrix} \mathbf{u} \\ \mathbf{v} \end{pmatrix} = \begin{pmatrix} \mathcal{U}\_1 \\ \mathcal{V}\_1 \end{pmatrix} \mathbf{c}\_{\upsilon} + \begin{pmatrix} \mathcal{U}\_2 \\ \mathcal{V}\_2 \end{pmatrix} \mathbf{d}\_{\upsilon} \tag{104}
$$

where **c***<sup>v</sup>* and **d***<sup>v</sup>* are constant vectors in the independent solution space. To distinguish the inner and outer vacuum solutions, we let **c***v*<sup>1</sup> and **d***v*<sup>1</sup> denote the constants for inner vacuum region and **c***v*<sup>2</sup> and **d***v*<sup>2</sup> for outer vacuum region.

In the outer vacuum region, the scalar potential **u** is subjected to *M* boundary conditions at infinite *ψ*. With these M boundary conditions imposed, there are only *M* independent solutions left. Without loss of generality, we can set **c***v*<sup>2</sup> to be zero in this case. Consequently, eliminating **d***v*<sup>2</sup> in Eq. (104), we obtain

$$\mathbf{u}|\_{\psi\_{\mathbb{b}+}} = \mathcal{T}\mathbf{v}|\_{\psi\_{\mathbb{b}+}\prime}$$

where the *<sup>M</sup>* <sup>×</sup> *<sup>M</sup>* matrix <sup>T</sup> is given by <sup>T</sup> <sup>=</sup> <sup>U</sup>2V−<sup>1</sup> <sup>2</sup> |*ψb*<sup>+</sup> . The matrix T can be computed by means of the Green function method Chance (1997).

In the inner vacuum region, the independent solutions can be constructed, for example, with the use of an inward numerical shooting Zheng & Kotschenreuther (2006), with the following boundary conditions imposed at *<sup>ψ</sup>b*−:

$$
\begin{pmatrix} \mathcal{U}\_1 \\ \mathcal{V}\_1 \end{pmatrix}\_{\psi\_{\mathbb{P}\_-}} = \begin{pmatrix} \mathcal{I} \\ \mathcal{O} \end{pmatrix}\_{\prime} \tag{105}
$$

$$
\begin{pmatrix} \mathcal{U}\_2 \\ \mathcal{V}\_2 \end{pmatrix}\_{\psi\_{\mathbb{P}\_-}} = \begin{pmatrix} \mathcal{T} \\ \mathcal{I} \end{pmatrix}\_{\prime} \tag{106}
$$

where O is *M* × *M* zero matrix. Since the boundary conditions in Eq. (105) give *δ***B** · ∇*ψ* = 0 at wall, these conditions correspond to a set of solutions that corresponds to the perfectly conducting wall type. On the other hand, since the boundary conditions in Eq. (106) guarantee that the independent solutions to be continuous with outer vacuum solutions, these conditions correspond to a set of solutions that corresponds to the no-wall type. Using the general expression for the solutions in Eq. (104), we can express the normal and parallel magnetic fields at the plasma-vacuum interface as follows:

$$[[\mathcal{J}\nabla\Psi \cdot \delta \mathbf{B}]] = -\mathcal{V}\_1 \mathbf{c}\_{v1} - \mathcal{V}\_2 \mathbf{d}\_{v1\prime} \tag{107}$$

$$-[[\mathcal{J}\mathbf{B}\cdot\delta\mathbf{B}]] = i\mathcal{Q}\left(\mathcal{U}\_1\mathbf{c}\_{v1} + \mathcal{U}\_2\mathbf{d}\_{v1}\right).\tag{108}$$

#### **5.3 Eigenvalue problem**

26 Will-be-set-by-IN-TECH

For completeness, in this subsection we briefly review the vacuum solutions in Ref. Zheng &

where *u* is the magnetic scalar potential and is related to the perturbed magnetic field by *δ***B** = −∇*u*. Here, we note that this representation of vacuum magnetic field, although being simple, excludes the consideration of *n* = 0 modes. To study *n* = 0 modes, one more scalar is needed to represent the vacuum magnetic field. For the sake of conciseness, we outline the

As in the plasma region, Fourier decompositions are introduced for both poloidal and toroidal directions to solve Eq. (103). Then Eq. (103) becomes a set of second-order differential equations of number *M* for **u**. This set of second-order differential equations can be transformed into a set of first-order differential equations of number 2*M*, by introducing an additional field **v** = −[[J ∇*ψ* · *δ***B**]], which is related to the magnetic scalar potential in Fourier

There are 2*M* independent solutions for Eq. (103), which can be used to construct the following

The general solutions in the vacuum regions can be expressed as a linear combination of the

where **c***<sup>v</sup>* and **d***<sup>v</sup>* are constant vectors in the independent solution space. To distinguish the inner and outer vacuum solutions, we let **c***v*<sup>1</sup> and **d***v*<sup>1</sup> denote the constants for inner vacuum

In the outer vacuum region, the scalar potential **u** is subjected to *M* boundary conditions at infinite *ψ*. With these M boundary conditions imposed, there are only *M* independent solutions left. Without loss of generality, we can set **c***v*<sup>2</sup> to be zero in this case. Consequently,

**u**|*ψb*<sup>+</sup> = T **v**|*ψb*<sup>+</sup> ,

**<sup>u</sup>**1, ··· , **<sup>u</sup>***<sup>M</sup>* **<sup>v</sup>**1, ··· , **<sup>v</sup>***<sup>M</sup>*

**<sup>u</sup>***M*<sup>+</sup>1, ··· , **<sup>u</sup>**2*<sup>M</sup>* **<sup>v</sup>***M*<sup>+</sup>1, ··· , **<sup>v</sup>**2*<sup>M</sup>*

*∂ψ* <sup>+</sup> �*i*J ∇*<sup>ψ</sup>* · ∇*θ*�M**u**.

 ,

> .

**d***v*, (104)

<sup>2</sup> |*ψb*<sup>+</sup> . The matrix T can be computed by

<sup>∇</sup>2*<sup>u</sup>* <sup>=</sup> 0, (103)

Kotschenreuther (2006). The vacuum regions are described by the Laplace equation

general solutions for inner and outer vacuum regions simultaneously.

**v** = J |∇*ψ*| 2 *∂***u**

> U1 V1 ≡

> U2 V2 ≡

 **u v** = U1 V1 **c***<sup>v</sup>* + U2 V2 

region and **c***v*<sup>2</sup> and **d***v*<sup>2</sup> for outer vacuum region.

where the *<sup>M</sup>* <sup>×</sup> *<sup>M</sup>* matrix <sup>T</sup> is given by <sup>T</sup> <sup>=</sup> <sup>U</sup>2V−<sup>1</sup>

means of the Green function method Chance (1997).

eliminating **d***v*<sup>2</sup> in Eq. (104), we obtain

**5.2 The solution of vacuum region**

space as follows:

independent solution matrices:

independent solutions:

The solutions in the plasma and vacuum regions described in the last two subsections can be used to construct the eigen value problem Zheng & Kotschenreuther (2006). The normal magnetic field component and the combined magnetic and thermal pressures are required to be continuous at the plasma-vacuum interface. Matching plasma [Eqs. (101) and (102)] and vacuum [Eqs. (107) and (108)] solutions at the interface *ψ<sup>a</sup>* gives

$$\mathbf{d}\_{v1} = \mathcal{F}\_1^{-1} \delta \mathcal{W}\_b \delta \mathcal{W}\_{\infty}^{-1} \mathcal{F}\_2 \mathbf{c}\_{v1} \tag{109}$$

$$\text{where } \delta \mathcal{W}\_{\infty} = \mathcal{W}\_{\mathcal{P}} - \mathcal{Q}[\mathcal{U}\mathcal{Q}\_2\mathcal{V}\_2^{-1}]\_{\mathfrak{F}\_1} \mathcal{Q}\_\ast \delta \mathcal{W}\_{\mathfrak{b}} = \mathcal{W}\_{\mathcal{P}} - \mathcal{Q}[\mathcal{U}\_1\mathcal{V}\_1^{-1}]\_{\mathfrak{F}\_1} \mathcal{Q}\_\ast \mathcal{F}\_1 = \mathcal{Q}\left[\mathcal{U}\_2 - \mathcal{U}\_1\mathcal{V}\_1^{-1}\mathcal{V}\_2\right]\_{\mathfrak{F}\_1}.$$

and F<sup>2</sup> = Q <sup>U</sup><sup>1</sup> − U2V−<sup>1</sup> 2 V1 *ψa* . Note that *δ*W<sup>∞</sup> and *δ*W*<sup>b</sup>* correspond to the energy matrices without a wall and with a perfectly conducting wall at *ψb*, respectively, as can be seen from the boundary conditions in Eqs. (105) and (106).

We now consider the matching across the thin resistive wall. For the radial magnetic field, the Maxwell equation ∇ · *δ***B** = 0 and the thin wall assumption lead to

$$\mathbf{v}|\_{\varnothing\_{b-}} = \mathbf{v}|\_{\varnothing\_{b+}} = \mathbf{d}\_{v1}.$$

The current in the resistive wall causes a jump in the scalar magnetic potential. This can be obtained from the Ampére law

$$
\nabla \times \nabla \times \delta \mathbf{B} = -\gamma \mu\_0 \sigma \delta \mathbf{B} \,, \tag{110}
$$

in Toroidal Plasma Confinement 29

Overview of Magnetohydrodynamics Theory in Toroidal Plasma Confinement 29

In this chapter we have given an overview of MHD theory in toroidal confinement of fusion plasmas. Four types of fundamental MHD modes in toroidal geometry: interchange, ballooning, TAEs, and KDMs, are discussed. In describing these modes we detail some fundamental analytical treatments of MHD modes in toroidal geometry, such as the average technique for singular layer modes, ballooning representation, mode coupling treatment in TAEs/KDMs theories. Note that analytical approach is often limited for toroidal plasma physics. Global numerical treatment of MHD modes is also reviewed in this chapter, especially the AEGIS code formalism. These theories are reviewed in ideal MHD framework. Here, we briefly discuss kinetic and resistive modifications to ideal MHD, as well as the

Let us first discuss kinetic effects. Since strong magnetic field is used to contain plasmas in magnetically confined fusion experiments, MHD theory can be rather good to describe fusion plasmas in the direction perpendicular to magnetic field. This is because strong magnetic field can hold plasmas together in perpendicular movement. Therefore, MHD is a very good model to describe perpendicular physics, if FLR effects are insignificant. However, in parallel direction the Lorentz force vanishes and particle collisions are insufficient to keep particles to move collectively. Consequently, kinetic description in parallel direction is generally necessary. Kinetic effect is especially important when wave-particle resonance effect prevails in the comparable frequency regime *ω* ∼ *ωsi* Zheng & Tessarotto (1994a). In the low frequency regime *ω* � *ωsi*, wave-particle resonances can be so small that kinetic description results only in an enhancement of apparent mass effect. Kinetic effect in this case can be included by introducing enhanced apparent mass. Another non-resonance case is the intermediate frequency regime *ωsi* � *ω* � *ωse*. In this regime kinetic description results in a modification of ratio of special heats. By introducing proper Γ MHD can still be a good approximation. Recovery of perpendicular MHD from gyrokinetics has been studied in details in Ref. Zheng

Fig. 2. CITM physics picture. The dot-dashed line represents mode rational surface. Perturbed current at rational surface due to interchange modes leads to field line

reconnection and formation of magnetic islands.

**6. Summary and discussion**

et al. (2007).

connection of MHD instabilities to transport.

where *σ* is the wall conductivity. Equation (110) can be reduced to

$$\mathcal{V}\left(\mathbf{u}|\_{\psi\_{\mathbb{H}+}} - \left.\mathbf{u}\right|\_{\psi\_{\mathbb{H}-}}\right) = \tau\_{\mathcal{U}}\gamma\_{N}\mathbf{d}\_{\upsilon1\prime} \tag{111}$$

where *τ<sup>w</sup>* = *μ*0*σdb*/*τA*, *d* is the wall thickness, *b* is the average wall minor radius, and

$$\begin{split} \mathcal{V} &= \mathcal{M} \left< \mathcal{J} |\nabla \psi| |\nabla \theta| - \mathcal{J} |\nabla \psi \cdot \nabla \theta|^2 / (|\nabla \psi| |\nabla \theta|) \right> \mathcal{M} \\ &+ n^2 \left< \mathcal{J} |\nabla \phi|^2 |\nabla \psi| / |\nabla \theta| \right>. \end{split}$$

Since **c***v*<sup>2</sup> = 0, we find that Eqs. (104) - (106) yield

$$\|\mathbf{u}\|\_{\psi\_{\mathbb{H}+}} - \mathbf{u}\|\_{\psi\_{\mathbb{H}-}} = -\mathbf{c}\_{\upsilon 1}.\tag{112}$$

From Eqs. (109), (111), and (112) we find the eigen mode equations

$$\mathcal{D}\_0(\gamma\_N) \mathbf{d}\_{v1} \equiv \tau\_w \gamma\_N \mathbf{d}\_{v1} + \mathcal{V} \mathcal{F}\_2^{-1} \delta \mathcal{W}\_\infty \delta \mathcal{W}\_b^{-1} \mathcal{F}\_1 \mathbf{d}\_{v1} = \mathcal{O} \dots$$

The dispersion relation for this eigen value problem is given by the determinant equation det |D0(*γN*)| = 0. In general the Nyquist diagram can be used to determine the roots of this dispersion relation. For RWMs, however, the growth rate is much smaller than the Alfvén frequency. Therefore, the growth rate dependence of *<sup>δ</sup>*W∞*δ*W−<sup>1</sup> *<sup>b</sup>* can be neglected for determining the stability condition. Consequently, one can use the reduced eigen value problem

$$-\mathcal{V}\mathcal{F}\_2^{-1}\delta\mathcal{W}\_\infty\delta\mathcal{W}\_b^{-1}\mathcal{F}\_1\mathbf{d}\_{v1} = \mathfrak{x}\_\mathcal{w}\gamma\_N\mathbf{d}\_{v1} \tag{113}$$

with the RWM mode growth rate *γ<sup>N</sup>* on the right hand side of this equation used as the eigen value to determine the stability.

#### **5.4 Discussion**

Now let us discuss the connection of current global theory with localized analytical theories described in Sec. 4. The singular layer equation in Eq. (48) is derived by employing mode localization assumption. Only localized mode coupling is considered. The general eigen mode equation Eq. (96) in plasma region contains all side band couplings. Noting that Q ∝ *x*, one can see from Eqs. (97) and (98) that <sup>F</sup> <sup>∝</sup> *<sup>x</sup>*<sup>2</sup> and <sup>K</sup> <sup>∝</sup> *<sup>x</sup>* at marginal stability *<sup>ω</sup>* <sup>=</sup> 0. We can therefore see the root of Eq. (48) in Eq. (96). If ballooning invariance in Eq. (57) is introduced, the set of matrix Eq. (96) can be transformed to a single ballooning equation. The TAE theory in Sec. 4.4 uses just two Fourier components to construct eigen modes. The general Alfvén gap structure can be determined by det |F| = 0. Note that, if an analytical function is given on a line on complex *ω* plane, the function can be determined in whole domain through analytical continuation by using the Cauchy-Riemann condition. Note also that one can avoid MHD continuum by scanning the dispersion relation with real frequency �*e*{*ω*} for a small positive growth rate �*m*{*ω*}. Using the scan by AEGIS one can in principle find damping roots through analytical continuation. Due to its adaptive shooting scheme AEGIS can be used to compute MHD modes with very small growth rate. It has successfully computed Alfvén continuum damping rate by analytical continuation based on AEGIS code Chen et al. (2010).
