A numerical approach to the calculation of the Alfvén continuum in the presence of magnetic islands

Mode activity in the presence of a magnetic island has been experimentally observed in FTU without the presence of fast ions.¹ 

The theoretical investigation of the Alvénic spectrum of magnetic islands began roughly a decade ago with analytic work by Biancalani et al.^2–4 Subsequently, an alternative analytic approach using a WKB technique was developed by Cook and Hegna⁵ resulting in an analytic expression for the Alfvén spectrum of the branches with the same helicity as the island. Both methods consistently delivered an up-shift of the continuous spectrum inside and outside the island. However, the theory in Refs. 2 and 3 applies also to modes with a different helicity than the island. The arising gap can be the location where a global mode with little continuum damping may reside.

As such modes can be driven unstable by interaction with particles (either fast or thermal ones) and may cause transport, they need to be understood. The calculation of the Alfvén continuum is the first step in this direction. An island-induced Alfvén eigenmode (IAE) has indeed been found in a calculation with the SIESTA code⁶ for an equilibrium of the Madison Symmetric Torus (MST) explaining an observed Alfvénic mode activity during neutral beam injection (NBI).⁷ 

There have been experimental observations of the effect of the magnetic island on the beta induced Alfvén Eigenmodes (BAE) agreeing for moderate perturbation levels with the theoretical predictions.⁸ 

In TJ-II, discrete shear-Alfvén eigenmodes, which are presumably induced by a non-rotating magnetic island [magnetic island induced Alfvén eigenmodes (MIAE)],² have been observed experimentally.⁹ Recently, investigating resonant magnetic perturbations in J-TEXT, MIAE have been found,¹⁰ changing their frequency with the island consistently with the predictions by Biancalani et al.² Also for other experiments, such as the Wendelstein 7-X stellarator, magnetic islands have a certain importance, because this experiment uses a 5/5-island chain at its outer boundary for an island divertor. Moreover, due to the small magnetic shear, the islands are relatively large. Also in W7-X, there has been an experimental scan of the rotational transform where islands appeared also inside the plasma.¹¹ Low frequency mode observations and transport phenomena associated with island effects in that device have been discussed in Refs. 12 and 13.

The intention of this work is to provide a numerical solution to the calculation of an Alfvén continuum including islands, which can easily be applied to experimental situations in recent stellarators such as W7-X, LHD, or TJ-II.

In Sec. II, we re-derive the governing equations in Boozer coordinates and describe the implementation of their solution into the CONTI code, which is routinely used for W7-X and TJ-II. In Sec. III, we compare with earlier results to confirm the validity of our approach.

Then, in Sec. IV, we discuss issues concerned with numerics and convergence. Finally, we apply the extended code to a W7-X configuration with an island chain inside the plasma to demonstrate how a realistic island size may alter the Alfvén spectrum. For that purpose, we use the HINT code¹⁴ to find the island size.

We start with the unperturbed magnetic field without islands given in Boozer coordinates and follow the notation used earlier in Ref. 15. In short, $( s , ϑ , φ )$ (with $s ∈ [ 0 , 1 ] , ϑ ∈ [ 0 , 1 ]$⁠, and $φ ∈ [ 0 , 1 ]$ in one field period or $φ ∈ [ 0 , N P ]$ in the full torus, where N_p denotes the number of field periods) form a left-handed coordinate system with the normalized flux coordinate $s = F T / F T ( a )$ with the minor radius a and F_T being the toroidal flux. We introduce the island using a perturbing field $B I$ of the following form:

with α being

and for later use, we introduce $α ¯$ as

Then, the new magnetic field is

Throughout this work, we regard A as a small constant specifying the amplitude of the field perturbation. Due to its smallness, we need not consider corrections to the equilibrium force balance and can follow earlier work.^2,3,5

The unperturbed magnetic field vector is tangent to the unperturbed flux surfaces, i.e.,

and due to the perturbation, the flux surfaces are distorted to new flux surfaces $s *$⁠, such that

and $s * = s * ( s , ϑ , φ )$⁠.

We use the representation of B₀ in Boozer co-ordinates,

where $F P ( s ) ′$ and $F T ( s ) ′$ are the derivatives of the poloidal and the toroidal flux with respect to s, respectively. J(s) and I(s) are the total poloidal and toroidal current, while $β ̃ ( s , ϑ , φ )$ is a small, finite β correction term discussed in Refs. 16 and 17. That of B_I is

We insert both into Eq. (6) to finally obtain the following differential equation for $s * = s * ( s , ϑ , φ )$⁠,

This equation can be simplified to

with the shorthand notation for the radial derivative of a quantity g: $g ′ = ∂ g ∂ s$⁠. We solve it by making the following ansatz (with f a function to be determined),

and obtain

Using the relation for the equilibrium pressure,

we can estimate the last term on the right-hand side at the position s_I of the island, given by $ι I ( s I ) = N P F p ′ F T ′ = − n I m I$⁠, to be small (note, that $F P ′$ is the derivative of the poloidal flux per period),

Thus, the term is not only quadratic in the island amplitude, which is assumed to be small compared with the main field, but it is also related to the pressure gradient, which has to be zero throughout the island to avoid a diverging current density. Thus, we are left with

We solve the latter equation integrating over s and obtain for $s *$⁠,

an equation describing the new flux surfaces depending on the unperturbed coordinates $( s , ϑ , φ )$⁠. Figure 1 illustrates the arising topology of the flux surfaces. Topologically, the regions inside and outside the island can be distinguished. For orientation, we remark that $s * = 0$ is the O-point and $s * = ± A$ marks the separatrices.

Obviously, outside the island, $s *$ is a unique function of s and the angle $2 π ( m ϑ + n N p )$ if a cut is made at s = s₀. To illustrate that, we approximate (see  Appendix A) the integral in Eq. (19) with the result,

leading to

The sign ± distinguishes the sub-regions between the plasma center and the separatrix (−) and between the separatrix and the boundary (+). A one-to-one mapping of the “old” flux label coordinates $s , ϑ , φ$ (without the island) and the flux label $s * , ϑ , φ$ (with the island) can only be reached if a solution for each angle $2 π ( m I ϑ + n I N p )$ is possible, which is the case for $| s * | ≥ | A |$⁠. If the equality holds, the minimal deviation from s₀ is 0 and the maximal deviation is

to each side, resulting in the expression of

for the half width of an island measured in units of the normalized toroidal flux. While we have solved the problem of the coordinate mapping, we need to discuss a further mapping step. In this work, we aim at a description of islands that is suitable for comparison with experimental observations. Thus, when representing results, we will map the rather abstract $s *$-contour to the maximum s-value that it assumes for $s < s 0$ and to its minimum s-value if $s > s 0$⁠, as illustrated in Fig. 2(a). This way, the continua can be attributed to an approximate radial position. For the coordinates inside the island, we will provide an analogous mapping as well.

Again, we consider Eqs. (20) and (21), but now for the case that $| s * | ≤ | A |$⁠. The $s *$-contours are closed and we are inside the island. Then, not for all angles $2 π ( m I ϑ + n I N p φ )$⁠, a solution s is possible for a given $s *$ and the ± sign in Eq. (21) refers to the same flux surface. Consequently, we get two values of s for one given coordinate triple of $( s * , ϑ , φ )$ implying that this system of coordinates is not unique and needs to be replaced.

We chose a new angular variable β measuring an angle around the center of the island (see Fig. 3). At first, we define the auxiliary angle ξ,

allowing us to transform

such that we can write Eq. (20) as

with $C 1 = − 1 2 F T ′ N p ι 0 ′ I ( s 0 ) = 2 A w s 2$ being a constant. Solving now for $s − s 0$⁠, we obtain

This motivates the definition of the new angular variable β as

Using the new variable, Eq. (27) can be rewritten as

resembling the parametric representation of a shifted circle. An illustration, of how the angle β can be understood, is shown in Fig. 3. Finally, using Eqs. (28) and (24), we can write down a mapping from the new coordinate β to the old coordinate ϑ,

To complete the new coordinate system, we keep the toroidal angle $φ$ and arrive at the transform $( s , ϑ , φ ) ⇒ ( s * , β , φ )$⁠.

Taking the ϑ- and $φ$-derivatives of Eq. (28), we obtain

and

with k being

If $g *$ is the Jacobian of the new coordinates such that

it is straightforward to show that

if we keep in mind that from Eq. (28): $β = β ( s * , ϑ , φ )$⁠. With Eq. (35), the contravariant β-component of $B →$ can be expressed as

In order to reach converged frequencies in the island center, it is crucial to use some analytical cancelation in the leading order terms to occur, resulting in

(for the details, see  Appendix B).

It is worthwhile to investigate with what sort of angle the new angle β can be associated. Evaluating its definition, Eq. (28), at the separatrix leads to

i.e., β is basically an angle related to the helicity of the island. However, in a plane of a constant $φ$⁠, as Fig. 3 shows, it works as a poloidal angle of the island chain. Another choice of $φ$ will shift the island chain but not alter its shape. That means that the island has been transformed into a quasi-two dimensional elongated tube encircling the torus.

An illustration of this feature is provided in Fig. 4. The metric coefficient and the Jacobian $g$ are dominated by the m = 2 harmonic. (The m Fourier index is related to β.) In analogy to a tokamak, it reflects the elliptic shape of the $s * = const .$ surfaces. Specifically close to the separatrix other components get larger because of the more complex structure of the surfaces there.

Another important point is the extension of the domain of $φ$⁠. The ratio $− n I m I$ determines the helicity of the island. The number $gcd ( m I , − n I )$ gives the number of separated island tubes winding around the torus introduced by the perturbation in Eq. (1). [Note, that the function $gcd ( x , y )$ returns the greatest common denominator of both its arguments x and y.] To close on itself, such an island flux tube has to encircle the torus $m I gcd ( m I , − n I )$ times. Note that for stellarators, in virtually all types of numerical linear calculations, the computation is done on a single field period, i.e., $0 ≤ φ ≤ 1$⁠. This effort-saving approach is no longer possible since we now have $0 ≤ φ ≤ N p m gcd ( m I , − n I )$⁠.

The equation for the Alfvén continuum frequency ω in general, geometry is given in Ref. 18 

which is valid on each flux surface $s * = const$⁠. In this equation, ψ is the poloidal flux and ξ^s is the eigenfunction of the Alfvén continuum. Note that we have omitted the sound modes, which would have added a coupled second equation. Using $ψ = ψ ( s * )$⁠, we can write Eq. (41) as

With the coordinate representations derived before, we can straightforwardly write down the equations for the continuum inside and outside the island and we are left with a generalized eigenvalue problem in $ω 2$⁠.

Note that the equations for the full MHD spectrum, i.e., including sound waves are given by Cheng and Chance¹⁸ and many others as well. The sound waves couple to Alfvén waves mediated by the curvature of the magnetic field and the pressure gradient. As both, the effect of the magnetic island and the sound waves affect mainly low frequencies, their inclusion would be desirable. We have omitted the sound waves to limit the computational burden.

As new flux surfaces are given by $s * = s * ( s , ϑ , φ )$⁠, we obtain

The parallel derivatives are calculated straightforwardly as

Introducing a Fourier representation for ξ^s

leads to

with the following definitions of the matrix elements

resembling the equations originally implemented in the CONTI code.¹⁹ 

Formally, the equation for the Alfvén continuum is the same as that outside the separatrix. However, here, the new coordinates $s *$ and β allow us to write $s = s ( s * , β )$ [see Eq. (28)]. Therefore,

The parallel derivative of ξ^s is given by

where the representation of $B β$ is critical for the convergence of the numerical solution and is given by Eq. (39), while

The Fourier-space representation of ξ^s is now

leading to the following Fourier matrices identical to Eqs. (48) and (49) but with θ replaced by β

To solve the generalized eigenvalue problem formulated in Sec. II, we use the CONTI code, which has been applied to Alfvén and sound continuum calculations in three-dimensional devices.^19–21 

The CONTI code uses quantities in Boozer coordinates initially stored pointwise in one- or three-dimensional arrays. Thus, they have to be transformed into coordinate systems $( s * , ϑ , φ )$ or $( s * , β , φ )$⁠, respectively. For that purpose, all Boozer quantities need to be interpolated with one- or three-dimensional cubic splines (using the Princeton Spline and Hermite Cubic Interpolation Routines²²) so that they can be calculated at any position in the new coordinates allowing us to fill in the arrays needed for the computation. The kernels of the integral expressions for the matrix elements in Eqs. (48) and (49) or, likewise, in Eqs. (54) and (55) are initially calculated point-wise and Fourier transformed afterwards using the FFTW²³ library. Note that the grid points for the equilibrium quantities accumulate in the radial coordinate close to the separatrix at either side. For the calculation of the eigenvalues, the matrix elements are interpolated between these grid points.

After this step, the infrastructure of the CONTI code is engaged to solve the generalized eigenvalue problem using either the NAG library or LAPACK routines.

During the more involved post-processing, the mapping of each eigenvalue to its appropriate radial position is done as laid out in Secs. II A and II B (see Fig. 2)—amendments have been made in the code.

Previous analytic work treated the lowest part of the spectrum close to the x-point of the island, i.e., close to the resonant value of ι. If there is no island present, all branches of the spectrum with the same helicity as the island, i.e., with mode numbers m and n fulfilling the relation $ι r = − n I m I = − n m$⁠, go to zero at this point.

Figure 5(a) shows a part of the Alfvén continuum in a cylinder without an island, while Fig. 5(b) shows the spectrum of the same cylinder with an $m I = 3 , n I = − 1$ island present outside the island enclosed by the separatrix. The profile of the rotational transform has a constant shear to match the requirements of analytical theory. The other parameters of the calculation can be found in  Appendix D.

The orange symbols in Fig. 5(b) mark features emerging in the presence of an island. The most prominent one is the up-shifting branches with the same helicity as the island (lower ellipse). For comparison, the results of analytical theory,⁵ including the islands, are provided for the branches with the lowest m and n numbers. Agreement is best for the lowest branch. Close to the separatrix, the deviations are relatively large, but as discussed in Sec. IV B, the solutions converge further out.

While the analytic theory⁵ is limited to the region close to the separatrix and a single helicity, it takes notice of the up-shift but cannot describe it as part as a new spectral gap induced by the island. This—very small—gap is indicated by the orange arrows and is more clearly visible in the high resolution, Fig. 6. The upper ellipse marks the breaking-up of the spectrum at the separatrix in a rather unusual fashion, not forming a spectral gap where the branches bend away and leaving a frequency domain open. Also, the branches on the left-hand side of the separatrix end there with a different frequency than those on the other side. Usually, the difference between left both sides gets smaller with better convergence, i.e., more modes included. However, it is hard to judge if it completely disappears. The most likely explanation is that the $s *$ contours giving rise to those values are calculated on completely different radial positions, as can be seen from Fig. 1. These qualitatively new features related to the appearance of an island will be discussed in more detail below and illustrated in Figs. 7(a) and 7(b).

While, in Fig. 5, the intention was to give an impression of the continuum and to compare with the analytical result, in Fig. 6, the spectrum has been re-calculated using a larger number of Fourier modes (⁠ $− 36 ≤ m ≤ 36$ and $− 12 ≤ n ≤ 12$⁠). A gap appearing in the Alfvén spectrum is now clearly visible. As expected, it is due to the interaction of Alfvén waves with the magnetic island (⁠ $m I = 3 , n I = − 1$⁠) as the difference between the intersecting branches is $Δ m = 3$ and $Δ n = − 1$⁠. The gap does not reach the lowest part of the spectrum, but it can be shown that this is due to the slow convergence with the number of modes. In fact, the $m = 0 , n = 0$ branch is actually the lowest one (ω = 0) such that the $Δ m = 3 , Δ n = − 1$ gap opens when $m = 3 , n = − 1$ is shifted up.

Therefore, the upshift found in earlier analytic work^2–5 turns out to be the lowest part of the island-induced gap. It is worthwhile to repeat that the theory of Refs. 2 and 3 should be able to reproduce these results, while that in Refs. 4 and 5 applies only to modes with the same helicity as the island. As the perturbation causing the island is usually small, also the size of the emerging gap is small compared with gaps due to field inhomogeneities in tokamaks or stellarators. In principle, it is possible to find global modes in these gaps. Therefore, as introduced earlier by Biancalani et al.² in the following, we will call these hypothetical eigenmodes Magnetic Island induced Alfvén Eigenmodes (MIAE), mark them with the helicity of the gap and use it, as a shorthand, also to name the gap itself. Thus, we have a $MIAE 3 , − 1$ gap in Fig. 6.

In Figs. 7(a) and 7(b), the newly emerging local gap-like features are shown again [cf. Fig. 5(b)], but this time with the larger mode window at different frequencies. More precisely, it includes $ω / ω A = 1 / 3$ [Fig. 7(a)] and $ω / ω A = 1$ [Fig. 7(b)]. Additionally, at the position of the resonant surface with $ι = ι I$⁠, it is

for the mode numbers of two arbitrary intersecting branches. Consequently, without an island, more than two branches may intersect at this point. In the presence of an island, we observe that these intersection points split up and can apparently even reconnect [Fig. 7(b)]. In neither of the figures, an additional global gap seems to form. Interestingly, the emerging patterns are different from gaps, as here the sum and not the difference between the mode numbers of the branches in, e.g., Fig. 7(b) yield those of the island helicity. (The pairs are blue–green and yellow–red.) Note that the branches are reshaped close to the separatrix only.

Inside the separatrix, the structure of the spectrum is entirely different. In earlier work, this inner-island-spectrum has been described as well, using analytical theory.^2–5 The spectrum is primarily characterized by the smallness of the rotational transform within the island. (As the island tube is already following the global resonant $ι I$⁠, only a residual part due to the much smaller flux of the island field leads to a rotational transform inside). The branches were named with Fourier indices or other “quantum” numbers and with “odd” and “even” attributes referring to their eigenfunctions.

The approach used here reveals additional details. Note that the meaning of m, n used to describe the spectral branches has changed with respect to the region outside the separatrix. As laid out in Sec. II B, m is the Fourier index corresponding to the angle β, while n, as before, is the Fourier index for the toroidal angle $φ$⁠.

In Fig. 8(a), it is shown that there are indeed odd and even branches if the spectrum is restricted to n = 0 branches only. Here, the terms odd and even refer to the properties of the dominant m component of the eigenmodes as shown in Fig. 9, which is here its real part.

However, as is visible from Fig. 8(b), the odd/even feature disappears if also n = −1 contributions to the spectrum are considered in the calculations. Then, the usual naming and attribution of branches to pairs of (m, n) numbers are again in place. As in normal fusion devices, also the gaps can now be named. Also in earlier work,^2,3 the even gaps and especially $Δ m = 2$ where found to dominate. Looking at the properties of the associated angle β (see Sec. II B), also the name EAE gap (Ellipticity induced Alfvén Eigenmodes) for $Δ m = 2$ seems appropriate. However, as we will see later, also the much smaller odd components, will induce gaps in the spectrum.

For comparison, the analytical result from Ref. 5 is indicated with dashed black lines. Agreement is quite good apart from the fact that the analytical result misses some of the branches. It seems that, in the analytical theory, the finite extent of the regions between the gaps had not been accounted for. This might be due to the applied approximations (WKB). Note that the theory devised by Biancalani et al.^2–4 qualitatively agrees much better than the WKB approach as it is able to reproduce the splitting (already in the version with a single helicity only⁴) and the other continuum structure (⁠ $n ≠ 0$⁠) as well.^2,3

Close to the separatrix, again the convergence of the numerical results is poor. (See the next section for a detailed discussion.)

If more n harmonics (i.e., additionally to n = 0) are taken into account for the calculation, the space between the gaps fills in more and more. Figure 10 illustrates the change in the spectrum and the crossing of the new branches. Enlarged parts of those spectra are shown in Fig. 11. The bands split again and also gaps with odd m appear, which had not been found before. However, these gaps are quite small, closed toward the edge and, because of the small value of the rotational transform within the island, they are associated with very high mode numbers (⁠ $m ≈ 70$⁠). Having in mind that the m numbers inside the island are multiples of the island harmonics, would bring the mode numbers with respect to the full device to $m ≈ 200$⁠. This range is certainly outside the validity of ideal MHD.

Nevertheless, it has been shown that the Alfvén spectrum could be calculated for the whole cylinder with full radial resolution. It shows good agreement with analytical theory and even reveals additional details such that an application to realistic devices is possible.

Already in the cylindrical computations advantages and drawbacks of the chosen approach are visible and will be discussed in this paragraph.

As outlined in Sec. IV B, the new coordinates can be calculated with good accuracy out of a set of initial Boozer coordinates using cubic splines. This is reflected in the good convergence properties of the code away from the separatrix.

An example for the convergence outside the separatrix is given in Fig. 12. Here, only modes with the same helicity as the magnetic island (m, n) have been considered, i.e., their mode numbers are $m = k · m I$ and $n = k · n I$ with k being an integer. Comparison is made with analytical theory⁵ again. The more modes are included in the computation, the better the convergence. For a very high number of modes, however, the accuracy of the computation breaks down near the separatrix and even negative eigenvalues are obtained. This is a sign that rounding errors have taken over as the accuracy of the matrix elements does not allow further improvement.

It has been found that this behavior cannot be fought with an increased number of grid points for the equilibrium quantities. Usually, 200 × 200 points are taken in the angular direction, while in the radial direction, 400 points are being used, which are quadratically accumulated around the separatrix.

As is visible in Fig. 6, convergence issues also exist at rational surfaces where the branches with resonant mode numbers normally go to zero, which is not the case here. For ordinary equilibria, without islands, CONTI usually converges well also for these points, since straight-field line (Boozer) coordinates are being used. When adapting to other, non-field-aligned coordinates, similar convergence problems close to points where $k | | = 0$⁠, have been observed.

Consequently, for improved numerical behavior, a transition to straight-field-line coordinates would be beneficial.

Similar convergence behavior can be observed inside the island. To avoid diverging expressions and to achieve a reasonable accuracy here, the partial cancellation of terms was used as outlined in  Appendix B. As this was not always possible, furthermore, the grid for the angular variable β is slightly shifted to avoid that the zeros of the angular functions involving β are exactly met. The results are found to be largely independent of that shift.

Figure 13 shows that the n = 0 branches converge well except very close to the separatrix while the n = −1 branches require a very high number of modes to converge. This is due to the fact that the rotational transform inside the island is very small. Additionally, it approaches zero at the separatrix, which makes the problem even bigger. The non-converged branches point upwards like an animal tail. They fall into the foreseen area for the “bands” between the gaps the higher the number of modes is. For a very high mode number, again the rounding errors take over and even negative eigenvalues may appear as the accuracy of the matrix elements does not suffice.

But, in spite of these issues, the numerical approach demonstrated convergence, revealed new features of the spectra, and allows an estimate of gap sizes and upshift, which is accurate enough to be useful for experimental applications.

In the latest operational phase of W7-X, a scan with respect to different values of the rotational transform was made¹¹ where marked MHD activity due to islands depending on their sizes was observed. We have chosen one of the equilibria (FQM001) of this scan with a rather large island for a stellarator application of the CONTI code. We did not aim at explaining experimental observations, but rather demonstrate that the code is ready to provide valuable information about the Alfvén spectrum in a realistic equilibrium.

Usually, the calculations of equilibria in W7-X and other stellarator experiments are done using VMEC.²⁴ (A new implementation of the theoretical ideas with an improved radial discretization is made in the GVEC code^25,26) These codes assume the existence of nested flux surfaces, which is not necessarily satisfied in reality, as magnetic islands may appear at rational flux surfaces, due to the interaction with external and internal perturbations. They distort the flux surfaces in the vicinity of those locations and change the character of the equilibrium. The island size is not known from observations but has to be determined by a numerical calculation using an appropriate tool. In this work, we employ the HINT code¹⁴ for this purpose. In Fig. 14(a), the result of the calculation is shown using a plasma β of 1%. The island size of the inner 5/5 island is measured to be roughly $11 cm$⁠. The outer 5/4 island chain (not considered in the calculations of this paper) is also visible. To measure the island size with respect to the effective minor radius which corresponds to a flux surface label of the initial calculation with VMEC, the VMEC and the HINT calculation have been overlayed [Fig. 14(b)]. From the difference in the outer boundaries, the $r eff$-size of the island is determined to be $8.7 cm$⁠. From this value, the size in the normalized toroidal flux can be calculated, and, finally, the necessary size of the magnetic perturbation used in the island model of this paper. The results are given in  Appendix D.

Looking back at the cylinder calculation for the spectrum outside the island, similarly, we would expect that an island induced gap (MIAE gap) with the helicity $Δ m = 5$ and $Δ n = − 5$ should appear. Following earlier notation,^19,21,27 we single out the fivefold periodicity and speak of the $HAE 5 , − 1$ gap (i.e., a gap in which Helical Alfvén Eigenmodes with this $Δ m$ and $Δ n$ may reside). Equivalently, the gap may be named as $MIAE 5 , − 1$⁠.

In Fig. 15, it is seen that such a gap already exists somewhat away from the rational surface even if there is no island present (green symbols). At the resonance, all branches with the same helicity as the island go to zero as regular flux surfaces are assumed. For the continuum drawn with black symbols in Fig. 15, an island was added to the VMEC equilibrium using the model described above and taking the parameters emerging from the HINT calculation. Because of the considerable island size and the more complex geometry, the calculation resembles not only features known from the cylinder looked at in Sec. IV A. An up-shift of approximately $10.8 kHz$/ $13.6 kHz$ appears at the island position, which again is the lowest part of an island induced gap. This time, however, the $HAE 5 , − 1$ gap already existed beforehand (green symbols) and appears enhanced by the island effects. Unlike the cylinder, there is a difference in the up-shift between the plasma to the left and to the right of the island. Therefore, two values of the up-shift had been given. At a first glance, one might have argued that both values should be equal as they address the same point and, in fact, in the simple cylindrical equilibrium, they are equal for the lowest and almost converged eigenvalues. However, in a more complex equilibrium, the new flux surfaces left and right of the separatrix obviously “feel” the radial dependence of many equilibrium quantities which may lead to different results. Of course, it cannot entirely be ruled out that the convergence issue discussed above is responsible for this deviation.

For the case, we consider here, the islands are “natural,” meaning that they comply with the fivefold symmetry of the W7-X. Therefore, also the concept of the mode families discussed earlier^19,21,28 remains intact. While Fig. 15 was given for the N = 1 mode family meaning that only n with $n ≡ 1 mod 5$ can couple, the other mode families give a very similar gap structure, which is not shown here. In this respect, it is again seen that the up-shift of the spectrum at the resonance position is actually the lowest part of an opening gap, not a pure shift: For the N = 0 mode family, the lowest branch is, as in the cylinder, $m = 0 , n = 0$⁠, the up-shifted one is $m = 5 , n = − 5$⁠. For the N = 1 mode family, the lowest branch is the $m = − 1 , n = 1$⁠. It is not up-shifted. The lowest up-shifted branch is $m = 4 , n = − 4$⁠. For the N = 2 mode family, the lowest branch is $m = − 2 , n = 2$ while the up-shifted branch is $m = 3 , n = − 3$⁠. Always, the difference is $Δ m = 5 , Δ n = − 5$⁠, matching the island helicity.

In addition to this lowest order gap with respect to the island helicity, there are also higher-order gaps, like $MIAE 9 , − 2$ and $MIAE 11 , − 2$⁠. Interestingly, only their n mode numbers are integer multiples of the island $n I$⁠. The m numbers deviate from $2 n I$ by ±1. Note that these higher order gaps extend to higher frequencies as shown in Fig. 16. This figure is the continuation of Fig. 15 to frequencies higher than the toroidal Alfvén eigenmode (TAE) gap. Again the green symbols mark the continuum without an island present while the black symbols stand for the continuum with island effects. The dashed red line is at the position of the resonant rotational transform. Apart from the now visible third order gaps $MIAE 14 , − 3$ and $MIAE 16 , − 3$ which are rather small, the $MIAE 9 , − 2$ gap reaches a considerable size at small s values. Also the third order gaps deviate in their m number from $3 m I$⁠, but this time by ±2.

Figure 17 shows that also the inner part of the Alfvén spectrum is qualitatively similar to the cylindrical result. In Fig. 17(a), one can see that also in a more complex equilibrium the naming of the pure n = 0 continuum using the parity does not apply. (The naming with respect to the dominating m number does also not apply in this case what is not shown.) With also an n = 1 component present [Fig. 17(b)], the modes can be named as before and show a similar structure with alternating gaps and bands as before in a cylinder. Also, the convergence improvement with a higher number of modes is demonstrated. The $Δ m = 2$ gap varies between $13 − 14 kHz$ and approximately corresponds to the up-shift of $10.8 kHz$/ $13.6 kHz$ observed in the spectrum outside the island.

A numerical approach to the calculation of the Alfvén continuum with magnetic islands has been developed. It is incorporated into the CONTI code and can be applied to all magnetic equilibria, which can be calculated with VMEC/GVEC. One island can be introduced as a perturbation to the main field.

The code finds that the up-shift of the Alfvén continuum outside the island close to the island position found in earlier work is the actually lowest part of an island induced gap.

The agreement with the analytical results is generally satisfying, even excellent if the radial distance to the x-point is not too small. Otherwise, slow convergence and rounding errors are an issue. Although the structure of the spectrum is completely different inside and outside the island, the up-shift at the separatrix should be the same. This is approximately the case. While the outside value approximately approaches the analytical value, it does only converge very slowly converge inside the island since the island rotational transform is very small.

Nevertheless, the code converges well at the o-point of the island and over most of the island radius. The splitting with $Δ m = 2 k$ (k being an integer) gaps is in accordance with the results found by Biancalani et al.,^2–4 while the small gaps with an odd $Δ m$ found in this work had not been observed before. However, due to their small size and the very high mode numbers, it is unlikely that a prominent global mode may reside within the odd gaps.

The numerical convergence of the results is generally satisfactory, but it is relatively poor close to $k | | = 0$ positions and close to the x-point of the island. While the first problem could likely be tackled by using straight-field-line coordinates, the convergence close to the x-points would be certainly much harder to improve as here many modes are coupling. At higher frequencies, close to points where many branches cross, a splitting of the crossing point has been observed in the presence of an island. The local re-connection of the branches seems not to be linked to the formation of a gap.

For W7-X, islands for a realistic configuration (FQM001) have been calculated with the HINT code. Their width had been determined and was used to calculate the Alfvén spectrum, which turned out to be quite similar to that in a cylinder. It could be shown that the island can enhance already existing gaps outside the separatrix (somewhat away from the island position) and creates new (higher order) gaps at higher frequencies. The width of the gaps varies and may at some radii be comparable to already existing gaps in the equilibrium. It was shown that the up-shift of the spectrum at the island position itself, similar to the cylinder, is just the lowest part of the island-induced gap. Although the metric looks more complex than that of the cylinder, the three-dimensionality of the configuration seems not to add more features to the island's appearance in the Alfvén spectrum.

In its current state, the CONTI code can be applied to island-related frequency observations in all magnetic configurations that can be described with equilibria calculated with VMEC/GVEC and, thus, be used to attribute observed mode activity to the presence of a magnetic island.

The CONTI code can only calculate the continuum and, therefore, only provide a frequency range for island-induced gaps where a global mode may reside. For this reason, the convergence issues of the code appear to be tolerable.

The authors would like to thank A. Biancalani for his interest and his very useful comments, P. Helander and T. Andreeva for valuable discussions, and M. Borchardt and H. Leyh for their support with computer issues.

This work has been carried out within the framework of the EUROfusion Consortium, funded by the European Union via the Euratom Research and Training Program (Grant Agreement No. 101052200—EUROfusion). Views and opinions expressed are however those of the author(s) only and do not necessarily reflect those of the European Union or the European Commission. Neither the European Union nor the European Commission can be held responsible for them.

One of us, Jinjia Cao, was supported by a grant from the China Scholarship Council (CSC), File No. 201908430060, the key scientific research program of the Education Department of Hunan Province (Grant No. 20A417), the Hunan Nuclear Fusion International Science and Technology Innovation Cooperation Base (Grant No. 2018WK4009), the ITER Project of Ministry of Science and Technology (Grant No. 2022YFE03080002) and wishes to thank P. Helander and the Stellarator theory division of IPP for the kind hospitality.

The authors have no conflicts to disclose.

Axel Konies: Conceptualization (equal); Investigation (lead); Software (lead); Writing – original draft (lead). Jinjia Cao: Investigation (supporting); Software (supporting). Ralf Kleiber: Conceptualization (equal); Investigation (supporting). Joachim Geiger: Investigation (supporting); Software (supporting).

The data that support the findings of this study are available from the corresponding author upon reasonable request.

We approximate the flux function $s *$ as in Refs. 2 and 5. In Boozer coordinates, ι is defined as

using a local approximation in the vicinity of the island, and the poloidal flux $F P ′ ( s )$ can be approximated as

Note that $F T ′$ is a constant, and the safety factor $ι 0 = ι ( s 0 )$ at the island position and its derivative $ι 0 ′ = ι ′ ( s 0 )$ are taken to be constant. Furthermore, we set

Considering that the island helicity defines ι₀,

inserting into Eq. (19) gives