Modeling of frequency-sweeping Alfvén modes in the TJ-II stellarator

VMEC is an equilibrium solver built to model three-dimensional fusion devices. It uses a steepest descent method to find the minimum of the plasma potential energy function, which can be written as:

$$\begin{align} W = \int\left(\frac{|\mathbf{B}^2|}{2\mu_0} + \frac{p}{\gamma -1}\right) \textrm{d}^3x \end{align} \tag{ 5 }$$

where B is the magnetic field, p is the pressure, γ is the adiabatic index and the integral is taken over the volume of the plasma. The resulting VMEC equilibrium file serves as the input to both the STELLGAP and AE3D solvers. Additional details about the VMEC code can be found in [21].

STELLGAP [22] solves for the Alfvén continuum spectrum in a three-dimensional toroidal device. The continuum spectrum defines the frequencies, radial locations, and toroidal and poloidal mode numbers of shear Alfvén waves that are expected to be seen in a toroidal device with a particular equilibrium configuration (as defined by the input provided by VMEC).

To calculate the Alfvén continuum spectrum, STELLGAP solves the following equation in straight field line Boozer coordinates [25]:

$$\begin{align} \mu_0 \rho \omega^2 \frac{|\nabla\psi|^2}{B^2} E_{\psi} + \mathbf{B} \cdot \mathbf{\nabla} \left[ \frac{|\nabla \psi|^2}{B^2} \left(\mathbf{B} \cdot \mathbf{\nabla}\right)E_{\psi} \right] = 0 \end{align} \tag{ 6 }$$

where ρ is the plasma mass density, and E_ψ is the covariant ψ component of the electric field.

The resulting Alfvén continuum plots, showing the continuum frequency as a function of radius for the mode harmonics m and n involved, serve as a 'map' indicating the radial locations and frequencies of modes that are most likely to be driven unstable by high-energy particles. Compressibility effects can be optionally included in the calculation of these continuum plots in STELLGAP.

Modes are likely to be destabilized most easily by energetic particles when they are located in the gaps of the continuum, defined by a local minimum in a continuum line located above a local maximum in another continuum line. The frequencies of the modes in the continuum gaps do not cross the continuum lines and therefore avoid strong continuum damping, meaning they are weakly damped and easy to excite with energetic particles. Modes that are observed experimentally are generally expected to be those which are localized in the gaps of the continuum and in the following sections we will exploit this relationship to validate simulations against experiment. Note that in STELLGAP, due to the helicity in TJ-II (the sign of the iota profile, which in inputs is negative), toroidal mode numbers are reported as negative. All profiles shown in the following sections are plotted with iota positive.

AE3D [23] calculates the eigenfrequencies and mode structures of shear Alfvén waves in a three-dimensional device. The governing equations of AE3D are the ideal MHD Ohm's law and the vorticity equation. When combined and simplified, they lead to the following eigenvalue problem:

$$\begin{align} & \omega^2 \nabla \cdot \left( \frac{1}{v^2_A} \nabla \phi \right) + \left(\mathbf{B}\cdot \nabla\right) \left[\frac{1}{B} \nabla^2 \left(\frac{\mathbf{B}}{B} \cdot \nabla \phi \right) \right] \notag\\ & \quad + \nabla \zeta \times \nabla \left( \frac{\mathbf{B}}{B} \cdot \nabla \phi \right) \cdot \nabla \frac{J_{||0}}{B} = 0 \end{align} \tag{ 7 }$$

where $J_{||0}$ is the equilibrium parallel current, φ is the electrostatic potential, and ζ is the toroidal angle coordinate.

The resulting eigenfrequencies and mode structures can be used in conjunction with the continuum plot calculated by STELLGAP to identify which modes are located within gaps in the continuum, and therefore are most likely to be destabilized. AE3D does not include compressibility effects in its calculation of eigenfrequencies and mode structures.

Recall that STELLGAP provides the locations of gaps in the Alfvén continuum for a given equilibrium configuration, and that modes that are most likely to be destabilized are those which are localized inside the gap region. Also note that simulation input parameters are constrained to those resulting from the experimental analysis that was described in section 2.

In this section, we study the localization of the reversed shear-induced gap region for various input plasma ion densities. The goal of this analysis is to calibrate the plasma ion mass density, which is a required input to STELLGAP and AE3D, taking into account the constraints on mode number, frequency and ι_min that were set by experimental analysis. While electron density is measured experimentally in TJ-II, ion density is not; therefore, the experimentally-measured frequencies will be used to calibrate ion density. The relationship between Alfvén frequency and ion density is described by equations (1) and (2).

Input iota profiles for this study are shown in figure 5 and a detailed discussion about the particular functional shape of these iota profiles will be provided in section 4.2.2. For now, they should be taken as an ansatz that is constrained by the values of ι_min that were predicted for the m = 6, n = −9 mode as shown in figure 2.

Zoom In Zoom Out Reset image size

Figure 6 shows the results of the study for three ion density profiles with varying on-axis values; the starting point for the shape of these profiles is the experimental electron density as given in [15]. In these results, the ion to proton mass ratio is fixed to 1 for simplicity. The results show that as the central ion density is decreased, the area of the reversed shear-induced gap generally becomes wider and shifts upward in frequency. The three central ion densities, distinguished by the areas shaded in grey, blue, and green, are overlaid against a sampling of points from the sweeping mode observed in TJ-II discharge 27 806. Recall that the perturbed density spectrogram data from which these points are sampled is shown in figure 2.

Zoom In Zoom Out Reset image size

Figure 6 suggests that the true value of central ion density should fall somewhere between $1\times 10^{19}$ (1/m³) and $2\times 10^{19}$ (1/m³), which is roughly in line with the experimentally-observed electron density [15]. The profile with central density $1\times 10^{19}$ (1/m³) is selected as the modeling input. It will be shown in section 4.2.2 that the error on the resulting modeled mode frequencies is relatively insensitive to variations in density.

It should be noted that although the locations of the gaps shown in figure 6 were calculated with the particular iota profile shape from figure 5, the resulting gap location is measured at ρ = 0.65 and is therefore specifically dependent on the iota profile's value at ρ = 0.65 and largely independent of its values elsewhere. Therefore, in calibrating the density, we have made no additional assumptions that would affect the output beyond assuming the results of the analysis described in section 2 and shown in figure 2 (specifically, the values of the x-axis as mapped to the frequency of the observed modes) to be correct.

In the upcoming sections, we will attempt to reproduce the circled frequency-sweeping mode from figure 2 with two different types of ansatz rotational transform profiles: monotonically increasing and non-monotonic.

4.2.1. Monotonic iota profiles.

In this section, we demonstrate that the m = 6, n = −9 mode could not have been observed in a TJ-II discharge with a monotonic iota profile. Figure 7 shows ansatz iota profiles input to STELLGAP and AE3D. The functional shape of the profiles is identical to that of the TJ-II vacuum profile, and the profile is translated upward over five runs of the simulation to mimic the effect of dynamically changing the profile's minimum during a discharge.

Zoom In Zoom Out Reset image size

The continuum plot calculated in STELLGAP for each of the five iota profiles is shown in figure 8. Continuum lines are shown in blue, orange, green, red and purple, while individual eigenmodes with toroidal mode number n = −9 as calculated by AE3D are shown in black.

Zoom In Zoom Out Reset image size

From these plots, we can observe that there are no eigenmodes (black lines) contained within any gap formed by the continuum lines (with the exception of a radially-extended gap mode around 280 kHz that appears when the minimum value of the iota profile is 1.49, but this mode is outside the frequency range of interest). We can conclude, therefore, that all of the modes with n = −9 for this rotational transform profile shape are subject to continuum damping and therefore unlikely to be destabilized.

This is an expected result since the input iota profiles are monotonic (continuously increasing with increasing radius) and therefore have normal shear. In order for reversed-shear Alfvén eigenmodes such as the ones suspected to have been observed in TJ-II discharge 27 804 to be present, the rotational transform profile must be non-monotonic [26]. Therefore, a monotonic iota profile shape can be ruled out.

4.2.2. Non-monotonic iota profiles.

In the previous section, we saw that the TJ-II vacuum iota profile, which is monotonically increasing with ρ, does not induce continuum structures that would accommodate the observation of the m = 6, n = −9 frequency-sweeping mode in discharge 27 804. This behavior is expected, as the modes of interest were observed during discharges with neutral beam injection, which can modify the vacuum iota profile. Neutral beam injection induces a plasma current $I_{\textrm{p}l}$ which, when assumed to be small, can be related to perturbations to the equilibrium in a linear way.

The left plot in figure 3 shows the modification of the vacuum iota profile (shown in black) due to various values of NBI-induced plasma current. These currents are of the same order as those computed by the ASCOT5 Monte Carlo orbit-following code for TJ-II [27].

Co-directional neutral beam injection induces a positive plasma current I_pl. From the result in figure 3, we see that in both cases with positive plasma current (I_pl = $+1, +2$) the vacuum iota profile is modified to demonstrate strong magnetic shear.

Since we are modeling a discharge with positive plasma current, we create an ansatz iota profile with a strong degree of shear.

A second factor to consider when creating the ansatz iota profile for the strong shear non-monotonic case is the prediction of the location of the iota profile's minimum value. The methodology for arriving at the predicted value is described in detail in section 2. The result of the analysis was that the minimum value of the iota profile should occur at the radial location ρ = 0.65.

With both of these two factors accounted for, namely the profile's shear based on the experimental plasma current and the location of the profile's minimum based on the results of analysis from [15], we constructed the ansatz iota profiles shown in figure 5. These are the same profiles that were used as input for the density study discussed in section 4.1. Recall that in experiments in TJ-II, iota profile minimum can be varied linearly over a single shot via variation of the current in the device's magnetic coils. Therefore, each iota profile in the sequence can be thought of as a single snapshot in time of the evolution of the Alfvénic activity during the discharge.

Each of the five iota profiles shown in figure 5 were run through VMEC, STELLGAP and AE3D as described in section 3. In figure 9, results from both STELLGAP and AE3D are shown for each profile.

Zoom In Zoom Out Reset image size

The colored lines on each plot represent the Alfvénic continuum structure for the n = −1 mode family to which the n = −9 mode belongs (cf table 1). Recall that the Alfvén continuum structure indicates the expected frequencies and radial locations of Alfvén waves for a discharge with the respective iota profile. The color of each line indicates the mode's toroidal mode number (n).

Table 1. Table showing input toroidal (n) and poloidal (m) mode numbers for each mode family.

Family 1	Family 2	Family 3	Family 4
n = −1; m = 1–5	n = −2; m = 1–5	n = −3; m = 1–5	n = −4; m = 1–5
n = −5; m = 3–8	n = −6; m = 3–8	n = −7; m = 3–8	n = −8; m = 3–8
n = −9; m = 5–15	n = −10; m = 5–15	n = −11; m = 5–15	n = −12; m = 5-15
n = −13; m = 7–17	n = −14; m = 7–17	n = −15; m = 7–17	n = −16; m = 7-17
n = −17; m = 9–27	n = −18; m = 9–27	n = −19; m = 9–27	n = −20; m = 9-27

On each of these plots, there is a gap in the continuum structure localized approximately at ρ = 0.65, formed by the local minimum of the n = −13 continuum line located above the local maximum of the n = −9 continuum line. Recall that the minimum value of each iota profile, or the location of the reversal of the profile's magnetic shear, is also at ρ = 0.65. Based on these coincident radial locations, we conclude that this gap is a reversed-shear induced gap in the Alfvén continuum.

The black horizontal line on each plot represents the full width at half maximum for a mode that falls inside the aforementioned gap as calculated by AE3D, and the full mode structures for each time slice are shown in figure 10. Recall from section 3.2 that modes inside the gaps of the Alfvén continuum are the most likely to be destabilized, and therefore the ones most likely to be observed in experiments.

Zoom In Zoom Out Reset image size

Each plot in figure 9 also includes the modified n = −9 continuum line when compressibility effects, or coupling to sound waves, are included in the simulations [28]. The effect of including compressibility is to shift the continuum line upward in frequency, as demonstrated by the comparison between the green and pink continuum lines. The difference between the maximum frequency of the compressible and incompressible continuum lines increases with increasing value of ι_min, reaching a maximum of 30 kHz of separation in the ι_min = 1.490 case. In the incompressible case, the mode frequency calculated by AE3D moves downward in frequency along with the peak of the n = −9 continuum line (green). It is inferred that the same would be true of the relationship between the n = −9 compressible continuum line (pink) and the sound-coupled mode frequency, resulting in an overall upward shift in frequency of the mode at each value of ι_min. Given that AE3D does not include compressibility effects, it cannot be used to confirm the above hypothesis. For such an assessment, a code such as CAS3D could be used instead [29].

As the minimum value of the iota profile sweeps linearly upward, the frequency of the mode sweeps linearly downward, indicating that these modes can, in fact, be identified as reversed-shear Alfvén eigenmodes. A plot of the frequency of the mode versus the minimum value of the iota profile for which it was observed is shown in figure 11. Error bars are included, showing the calculated error in frequency when a relatively small error of ±0.005 in iota and ±10% in density are assumed. The error of ±0.005 on iota was chosen from the full measured error range of 0.004–0.037 in iota as shown in figure 11 of [15]. Each of these two errors are considered independent from one another, and are propagated through equations (1)–(3) via addition in quadrature. The vast majority ($\gt$90%) of the resulting error on the frequency value is due to the uncertainty in iota.

Zoom In Zoom Out Reset image size

We see that simulations in which the non-monotonic profile minimum is varied as predicted by analysis give decent agreement with experimental observations. For the first four simulated points, the measured mode falls comfortably within the simulation's error bars. The fifth point, at ι_min = 1.490, captures the second displayed measured mode within its error bar range when an error of ±0.01 in iota is considered.

We expect that the comparison between the experimental and computed mode frequencies would further improve by taking into account the compressibility effects that are not included in the AE3D results in figure 11. We expect that taking into account these effects would result in an upward shift of the calculated mode frequency with the magnitude of shift increasing with the minimum iota, as suggested by the STELLGAP results for the n = −9 continuum line shown in figure 9.

As compared with the monotonic profiles in the previous section, non-monotonic profiles can produce frequency-sweeping Alfvén cascades, which supports theoretical predictions.

As an extension of the study beyond the m = 6, n = −9 mode in the ι_min range of 1.47–1.49, we have carried out simulations for a wider ι_min range of 1.42–1.55 which is similar to the range explored experimentally (cf figure 2). Our objective in this extended study is to investigate whether additional cascades beyond the m = 6, n = −9 mode could occur according to modeling and, if so, with what mode numbers and frequencies.

First, STELLGAP simulations were carried out where the non-monotonic iota profiles considered in section  4.2.2 were shifted up or down to match the predicted ι_min values from the analysis results shown in figure  2. In these simulations, four mode families with input toroidal and poloidal mode numbers as shown in table 1 were considered. Moreover, the ι_min range of 1.42–1.55 was divided in four parts and in each part, the family of input toroidal mode numbers was selected to coincide with the family of the toroidal mode number n given by the post-experiment analysis in that ι_min range (deduced mode numbers m n⁻¹ are indicated by the labels of the red arrows in figure 2). Consequently, the mode families 3, 2, 1 and 4 were considered in the ι_min range of 1.445–1.465, 1.46–1.475, 1.47–1.49 and 1.54–1.55 respectively, while all other variables including the plasma ion density and the location of the minimum of the iota profile (ρ = 0.65), were held constant.

Note that this approach results in a simplification of the experimental situation for which the post-experiment analysis [15] suggested that the minimum value of the iota profile evolved radially outward from a minimum of ρ = 0.2 to a maximum of ρ = 0.8 during the time interval shown in figure 2. The radial evolution of the minimum value of the iota profile is neglected in this work to limit the number of independent variables in each simulation, though it could be considered in future work.

The resulting locations of gaps in the Alfvén continuum spectrum as a function of the minimum of the iota profile are shown in figures 12–15 for the first sweeping mode for each mode number family.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In their frequency and ι_min values, the gap locations could roughly accommodate some of the experimentally-observed modes shown in figure 2 (recall that the locations of destabilized cascade modes correspond to the location of gaps in the Alfvén continuum).

In order to corroborate the STELLGAP results shown in figures 12–15, a simple analytical model was constructed. Combining equations (1)–(3), we find that the Alfvén frequency for a mode with a toroidal mode number n and poloidal mode number m can be written as

$$\begin{align} f_\textrm{A} = C |m\iota - n| \end{align} \tag{ 8 }$$

where C is a constant and ι is the minimum value of the iota profile. The constant C = 625 kHz was determined using STELLGAP results with $f_\textrm{A}$ = 75 kHz for the m = 6, n = 9 mode at ι = 1.48 (cf figure 12). The allowed frequencies of Alfvén modes were then calculated with this simple formula for seven toroidal mode numbers in each family shown in table 1. The resulting modes with frequency below 250 kHz in the iota range of 1.42–1.55 are shown in figures 12–15 together with the STELLGAP results.

Comparing the gap boundaries calculated by STELLGAP for the first sweeping mode for each mode number family with the lines calculated by the simple analytical model shown in figures 12–15, we see that the analytical model reproduces STELLGAP points very well. The points marked with stars in figures 12–15 indicate locations where a continuum line n and another continuum line n + 4 coincide in frequency. These points represent locations in the frequency-ι_min space where a downward-sweeping cascade mode could appear. The interpretation of these points as the locations where cascades could be instantiated is supported by the results from STELLGAP and AE3D. According to STELLGAP and AE3D, it is in the gap of the maximum and minimum of the n = 9 and n = 13 continuum lines, respectively, that the observed n = 9, m = 6 mode appears. Note also that the n = 9, m = 6 mode follows closely the maximum of the n = 9 continuum line as the minimum value of iota increases (cf figure 9).

When there is a clear observation that the measured frequency reaches a minimum, as in figure 2, it can be used to determine the corresponding rational iota value. Figure 12 demonstrates that the n = 9, m = 6 mode is expected to reach its minimum value in frequency at the rational iota point of 1.5. The measured mode shown in figure 2 appears to reach its minimum at ι_min = 1.51. This suggests that the true error on the experimental measurement of iota is likely 0.01, which is included in the aforementioned full measured error range of 0.004–0.037 as shown in figure 11 of [15].

The simple analytical model correctly predicts the points where each of the four modes identified in figure 2 are instantiated. It also reproduces a trend seen in experimental results, namely that the frequency at which the downward-sweeping mode is instantiated increases with increasing ι_min. It also suggests additional modes which could explain the experimentally-observed modes in figure 2 which were not identified by the post-experiment analysis reported in [15], though it should be noted that the existence of a gap in the Alfvén continuum spectrum does not directly imply the existence of a cascade in that frequency range.