Kinetic ballooning modes as a constraint on plasma triangularity in commercial spherical tokamaks

Commercial viability of fusion reactors require low capital cost and high plasma performance. In these regards, a particularly promising class of design is the spherical tokamak (ST), which differs from a conventional tokamak by having a smaller aspect ratio ($A\,\equiv\,R_{0}/r_\mathrm{minor}\,\lt$ ∼2, where R₀ is the major radius and $r_\mathrm{minor}$ the minor radius). This reduces the (radial) machine size (hence cost), increases the plasma β (plasma pressure/magnetic pressure) [1] and has been found to improve some plasma stability properties [2, 3]. A recent study predicted an increase in fusion triple product $nT\tau_\mathrm{E}$ (where n, T and $\tau_\mathrm{E}$ are the density, temperature and energy confinement time respectively) by a factor of 3 could be achieved in an ST compared to a conventional tokamak with similar values of fusion power and field strength but different machine size [4]. STs are thus receiving significant attention in the fusion community. Several devices have recently been built or upgraded, such as NSTX-U [5], ST40 [6] and MAST-U [7]. Others are currently being designed or built, such as SMART [8] and STEP [9]. However, reactor-relevant equilibria must have sufficiently low turbulent transport to provide sufficient confinement time to maintain density and temperature profiles with modest heating power. The tuning of ST equilibria to optimise turbulence is an area of active research.

A novel approach to optimisation is 'negative triangularity', whereby the plasma poloidal cross section resembles a backward 'D' (rounded on the inboard, high-field side). Experiments in TCV [10] and DIII-D [11] reported reduced turbulence with negative triangularity plasmas. For DIII-D, which operated negative triangularity in a limiter configuration, stored energy increased by $25\%$ and electron energy confinement time by $26\%$ compared to L-mode positive triangularity shots. However, no L-H transition occurred for the negative triangularity discharge, even at high heating power. This was found to coincide with the (modelled) H-mode pedestal being unstable to (toroidal mode number) $n\,=\,\infty$ ideal ballooning modes (IBMs) [12], an instability driven by normalised plasma pressure gradient $\partial \beta / \partial \psi$ (where ψ is the poloidal flux function). Saarelma et al [12] proposed that L-H transition in negative triangularity is only achievable when the pedestal is stable to IBMs for all pressure gradients; the so-called 'second stability window'. Without second stability access, strongly driven turbulence above a critical $\partial \beta / \partial \psi$ prevents pedestal formation and H-mode is suppressed. IBMs are known to be destabilised by negative triangularity [13], and this raises the question of whether negative triangularity is viable in commercial STs, which typically access second stability across the full radius [14].

Although the above discussion portrays IBMs as the responsible instability, the mechanism for the suppression is spatially small (∼ion Larmor radius) turbulence, for which kinetic effects cannot, in general, be ignored. The relevant instability is therefore the kinetic ballooning mode (KBM) [15, 16], which may be considered the 'kinetic analogue' of the IBM. While KBM stability is often assumed to track the IBM, it may be altered by effects such as diamagnetic stabilisation [17] or destabilisation by trapped particles [18]. A kinetic theory is thus required for turbulence studies.

In this paper, we seek to answer whether KBMs prohibit strongly negative triangularity in commercial ST reactors, and whether KBMs are likely to be problematic in positive triangularity ST equilibria. To do this, we use gyrokinetic simulations to examine the ion-scale instabilities in three strongly shaped hypothetical ST equilibria, constructed with a commercial reactor in mind.

Our equilibria (which are publicly available [19]) were generated by the fixed-boundary equilibrium code SCENE [20], which simultaneously solves the Grad–Shafranov equation and the neoclassical current contributions to generate an equilibrium with self-consistent current profiles. Some important quantities regarding these equilibria are shown in table 1, and radial profiles of electron density $n_\mathrm{e}$, electron temperature $T_\mathrm{e}$ and safety factor q in figure 1. For simplicity we take $T_\mathrm{e}\,=\,T_\mathrm{i}$. Since fusion αs preferentially heat electrons, it is plausible that $T_\mathrm{e}$ would exceed $T_\mathrm{i}$ in practice. However, the focus of our study is ballooning modes, which are typically sensitive to the total pressure gradient rather than to the contribution of individual species.

Zoom In Zoom Out Reset image size

Table 1. Some key quantities for the equilibria used in this study.

	'−ve tri'	'high $q0$'	'low $q0$'
Aspect ratio	1.67	1.67	1.67
R0	2.50	2.50	2.50
$q0$	2.58	2.71	1.38
$q_\mathrm{LCFS}$	4.50	8.97	9.16
$\delta_\mathrm{LCFS}$	−0.300	0.543	0.543
$\kappa_\mathrm{LCFS}$	2.80	2.80	2.80
β	$18.1\%$	$18.6\%$	$18.5\%$
βN	4.21	5.47	5.50
Toroidal current (MA)	16.5	16.5	16.5
Toroidal B (magnetic axis) (T)	1.93	1.83	1.83
Plasma volume (m3)	310	287	288
Heating power (MW)	60	60	60
Fusion power (MW)	514	808	839

The equilibria were constructed with identical major radius and elongation on the last closed flux surface (LCFS), total toroidal current and auxiliary heating/current drive power. The density and temperature profiles are reasonable models, but not the results of transport simulations. These equilibria operate in H-mode, though the conclusions regarding core microstability should be broadly applicable to commercially relevant L-mode plasmas. Another benefit of using H-mode equilibria is that it may give indication as to whether L-H transition is feasible (similar to Saarelma et al [12]), although that is not the focus of this work.

Two equilibria (which we label 'high $q0$' and 'low $q0$') were given strongly positive LCFS triangularity (0.543), the main difference being their on-axis safety factor $q0$. Both are unstable to external kink modes, though the 'high $q0$' case is stabilised by a wall at $r_\mathrm{wall}\,=\,1.4 r_\mathrm{minor,LCFS}$. The 'low $q0$' case remains unstable to $n\,=\,1$ external kink even at $r_\mathrm{wall}\,=\,r_\mathrm{minor,LCFS}$, and so is not a robust equilibrium from a magnetohydrodynamic (MHD) perspective. The third ('−ve tri') has LCFS triangularity of −0.3. Although this leads to a greater plasma volume, the density pedestal gradient is lower, resulting in a reduced fusion power of $40\%$ compared to the positive triangularity cases. MHD stability has not been examined for this equilibrium.

For this work, we use the freely available local gyrokinetic software GS2 [21] to (a) calculate the ideal ballooning stability for our equilibria and (b) perform local $\delta f$ gyrokinetic simulations. These are described in the following sections.

3.1. Ideal ballooning stability analysis

We calculate the $n\,=\,\infty$ ballooning stability using the GS2 module ideal_ball , which has previously been used for IBM studies in other devices [22].

An example of IBM stability for a flux surface as a function of normalised magnetic shear $\hat{s}\,\equiv\,\partial q / \partial \psi_N$ and normalised pressure gradient $\alpha\,\equiv\,\beta \cdot (1/p) \cdot \partial p / \partial \psi_N$ ($\psi_N\,\equiv\,\psi/\psi_\mathrm{LCFS}$ will be used as our radial coordinate in this work) is shown in figure 2, where we consider the cyclone base case [23]. For $\hat{s}\,\gt\,\hat{s}_\mathrm{min}$, the plasma is ballooning unstable for $(\alpha_\mathrm{crit} (\hat{s})\,\lt\,\alpha\,\lt\,\alpha_\mathrm{2nd} (\hat{s}))$; $\alpha_\mathrm{crit}$ is the boundary separating the first stable region from the unstable region, and $\alpha_\mathrm{2nd}$ separates the unstable and second stable regions. $\hat{s}\,\lt\,\hat{s}_\mathrm{min}$ represents the second stability window, in which the flux surface is stable for all α. The window size (i.e. the value of $\hat{s}_\mathrm{min}$) is particularly important for ST equilibria (e.g. [14]), since a high core β requires $\alpha(\psi_N)$ to be large, and typically $\alpha\,\gt\,\alpha_\mathrm{crit}$ is required.

Zoom In Zoom Out Reset image size

3.2. Gyrokinetic analysis

In the '−ve tri' equilibrium, we find that the second stability is blocked for the entire plasma ($\hat{s}_\mathrm{min}\,\lt\,0$) and $\alpha(\psi)\,\gt\,\alpha_\mathrm{crit}(\psi)$ (figure 3). As a result, the equilibrium is ideal MHD ballooning unstable over the core. In the pedestal, α strays into the second stable region ($\alpha\,\gt\,\alpha_\mathrm{2nd}$), coinciding with a sharp drop in the normalised gyrokinetic growth rate $\gamma(\psi_N)$ (lower plot). This supports the picture (validated in section 6.1) that, where the plasma is IBM-unstable, long-wavelength KBMs are the dominant instability. It should be noted that the gyrokinetic growth rate γ is normalised to the thermal velocity associated with the flux surface, so the absolute growth rates have a different (but qualitatively similar) trend.

Zoom In Zoom Out Reset image size

At $\psi_N\,\gt$ ∼0.97, there are strongly growing instabilities, as the equilibrium again becomes ballooning-unstable. However, our estimated $\vert \omega_\mathrm{HB} \vert$ is high in this region, so flow shear (as well as collisions) would likely need be included in simulations to give a realistic picture of stability.

Figure 4 shows $\gamma(k_y\rho_\mathrm{r}, \theta_0)$ for several values of ψ_N . For $\psi_N\,=\,0.5,0.8$, the instability spans a wide a range of low $k_y \rho_\mathrm{r}$, suggestive of an 'ideal' KBM, largely governed by ideal MHD physics [17, 18]. These KBMs are also wide in θ₀, and thus less easily stabilised by flow shear. $\psi_N\,=\,0.9$ shows the KBM being less unstable for modes with $k_y\rho_\mathrm{r}\,\lt$ ∼0.1, before the onset of a tearing parity mode at $k_y\rho_\mathrm{r}\,\lt$ ∼0.01.

Zoom In Zoom Out Reset image size

Although nonlinear gyrokinetic simulations would be needed to calculate the turbulent fluxes of particles and energy, our linear results indicate that those fluxes would be large. It should be noted that reliable nonlinear results may be difficult to obtain, since electromagnetic nonlinear simulations commonly fail to saturate for equilibria beyond marginal stability [30].

4.1. Second stability access in the negative triangularity regime

Section 4 supports a simple picture: KBM growth rates tracking IBM stability, and likely unacceptably high for steady-state operation. This prompts the question: is the closing of second stability a general feature of negative triangularity STs, or could a configuration exist with second stability access?

To test this we parametrise the $\psi_N\,=\,0.5$ surface using the model developed by Miller et al [31], which provides an analytic description of the flux surface in terms of nine dimensionless parameters ('Miller parameters') (a description of our parametrisation process is given in appendix C). We then test the effect of varying the triangularity (δ) and elongation (κ) in isolation, and find (in agreement with [13]) that this strongly affects stability. The role of reducing δ is to reduce $\hat{s}_\mathrm{min}$ and reduce $\alpha_\mathrm{crit}$ i.e. to close the second stability window and shrink the first stability region. For negative δ, the effect of increasing κ is to increase $\alpha_\mathrm{crit}$ but reduce $\hat{s}_\mathrm{min}$. This is illustrated in figure 5.

Zoom In Zoom Out Reset image size

These results indicate that, where ballooning modes are concerned, $\delta\,\lt\,0$ is not an attractive option, unless β is sufficiently low that the equilibrium remain in the first stable region. For $\psi_N\,=\,0.5$ for example, $\alpha_\mathrm{crit}\,=\,0.16$, compared with the equilibrium value of $\alpha\,=\,0.42$. Fixing the plasma size, a low-β device would require a large magnetic field. While high temperature superconductors may provide a possible pathway, they present engineering challenges.

It should be noted that the other Miller parameters, such as q, also affect $\hat{s}_\mathrm{min}$. A general study of the parametric dependencies of the IBM is beyond the scope of this work, but we note that relatively large changes to the Miller parameters of this surface are required to open second stability at fixed $(\delta, \kappa)$. As an example, achieving $\hat{s}_\mathrm{min}\,\gt\,0$ for our negative triangularity scenario requires increasing q from 3.0 to 5.4 with all other parameters kept fixed.

We also explored adjusting the magnitude of the local Shafranov shift parameter, $\vert \partial R_0/\partial \psi_N\vert$. This is generally higher in the positive triangularity equilibria than for negative triangularity (e.g. 'high $q0$' has 10% greater shift than '−ve tri' at $\psi_N\,=\,0.5$), probably due to the reduced core pressure in the latter, and one may wonder whether this modifies stability. We found that $\hat{s}_\mathrm{min}$ decreased (i.e. became more negative) as the magnitude of shift increased from its original Miller-fitted value. Whilst this is not exactly the same as recalculating the global equilibrium with an increased core β_N , this result suggests that attempting to operate at higher β is unlikely to alleviate the problem of second stability access.

We now consider the positive triangularity equilibria. Consistent with remarks made in section 4.1, we find these have better IBM stability, with $\hat{s}_\mathrm{min}\,\gt\,0$ everywhere as shown in figure 6. Despite this, 'low $q0$' is unstable near the magnetic axis ($\psi_N\,\lt$ ∼0.5), illustrating the importance of current distribution (via its effect on q) on ballooning stability. A 'high $q0$' is achieved with a hollow current profile, and this is IBM-stable across the plasma.

Zoom In Zoom Out Reset image size

The gyrokinetic $\gamma(\psi_N)$ (shown in figure 7) show similar behaviour between 'high $q0$' and 'low $q0$' for $\psi_N\,\gt$ ∼0.5 (consistent with the similarity of $q, \kappa, \delta$). In this region, instabilities with shorter wavelengths ($k_y \rho_\mathrm{r}\,\sim\,0.5$) are dominant for $\theta_0\,=\,0$, peaking near the pedestal top. However, as shown in figure 8, $\gamma (\theta_0)$ is very narrow, and hence highly susceptible to $\boldsymbol{E}\,\times\,\boldsymbol{B}$ shear stabilisation.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

However, $\gamma(\psi_N)$ differs in the core ($\psi_N\,\lt$ ∼0.5), coinciding with differences in IBM stability. 'low $q0$' exhibits a long-wavelength 'ideal' KBM ($\gamma(k_y \rho_\mathrm{r})$ shown in figure 8), consistent with the plasma being on or over the ideal ballooning boundary. 'high $q0$' has small but finite γ at long wavelength. Identification of these dominant instabilities is discussed in the following section.

For both equilibria, $\gamma(\psi_N)$ is lower than for '−ve tri', unstable over a smaller range of $k_y \rho_\mathrm{r}$ (figure 8) and of similar magnitude to $\vert \omega_\mathrm{HB} \vert$ over much of the plasma. This represents a dramatic improvement in microstability, as a direct result of making the triangularity positive.

6.1. Negative triangularity

Here we seek to identify the core and pedestal instabilities in the '−ve tri' equilibrium, by considering $(k_y \rho_\mathrm{r}\,=\,0.05,$ $\theta_0\,=\,0)$ for the $\psi_N\,=\,0.5$ (core) and $\psi_N\,=\,0.9$ (pedestal) flux surfaces.

We first verify that the dominant instability tracks the IBM stability boundary in $\hat{s}$–α space (as shown in figure 9). For both surfaces, the complex frequency Ω tracks the IBM boundary fairly well; γ is high across the unstable region and the frequency is smooth and in the ion diamagnetic direction ¹ .

Zoom In Zoom Out Reset image size

We also tested whether the mode was sensitive to α, but not specific contributions ($\partial T/\partial \psi_N$, $\partial n/\partial \psi_N$, β). This was done by scanning the quantity $f\,\equiv\,\beta \cdot(1/T) \cdot (\partial T/\partial \psi_N)$, in three cases: (1) scaling $(1/T) \cdot (\partial T/\partial \psi_N)$ at fixed β, adjusting $(1/n) \cdot (\partial n/\partial \psi_N)$ to keep α constant (2) scaling $(1/T) \cdot (\partial T/\partial \psi_N)$ at fixed β but keeping $\eta\,\equiv\,(n \cdot \partial T/\partial \psi_N)/$ $(T \cdot \partial n / \partial \psi_N)$ fixed (such that f changes α, with $f\,=\,0$ corresponding to $\alpha\,=\,0$) (3) keeping the gradients constant but scaling β. In all cases, we kept the magnetic geometry constant i.e. not changed to be consistent with α.

The result of this scan (shown in figure 10) shows Ω and the parallel electric field max $(E_\parallel(\theta) )\,\equiv\,$ max $(- \partial \phi / \partial \theta\,-\,\partial A_\parallel / \partial t)$ (normalised to the electrostatic potential φ) are constant with f when α is fixed (case 1). Cases (2) and (3) demonstrate that the kinetic gradients and β are approximately interchangeable in driving the mode (provided $\beta\,\gt\,0$). Moreover, the normalised mode structures ($\vert \hat{\phi} (\theta)\vert$, $\vert \hat{A}_\parallel (\theta) \vert$, $\vert \hat{B}_\parallel (\theta) \vert$) (figure 11) are virtually identical when f is scanned at fixed α for $\psi_N\,=\,0.5$ and very similar for $\psi_N\,=\,0.9$. This confirms the mode is indeed driven by α.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We conclude that the instability observed across the core is an 'ideal MHD'-like KBM, with the following features: (a) γ tracks the ideal MHD boundary well (b) the frequency is in the ion diamagnetic direction (c) the mode has a low parallel electric field $E_\parallel(\theta)$, indicative of MHD-like modes [18] (d) the mode amplitude is greatest in the bad curvature region and (e) has twisting parity (f) the mode preferentially occurs at long wavelength (g) the mode structure and complex frequency are sensitive to α but not to kinetic gradients or β individually. The pedestal, though IBM-stable, has all of the above properties, although γ is more sensitive to $(k_y \rho_\mathrm{r})$ at long wavelength (figure 4). We speculate this is a 'non-ideal' KBM, destabilised by some kinetic effect(s), allowing it to slightly exceed the ideal boundary.

6.2. Positive triangularity

We repeated this analysis for $(\psi_N\,=\,0.5, k_y\rho_\mathrm{r}\,=\,0.05, \theta_0\,=\,0)$ for the 'high $q0$' and 'low $q0$' equilibria. Departures from the ideal MHD stability boundary are seen in $\hat{s}$–α space (figure 12), with the KBM smoothly extending across the IBM boundary into the second stability window. To confirm this smooth transition, we performed a fine scan in $\hat{s}$, keeping α at its equilibrium value (shown in figure 13).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Like the KBMs discussed in section 6.1, the mode has a frequency in the ion diamagnetic direction and twisting parity. $E_\parallel$ is low, but rises smoothly as $\hat{s}$ leaves the IBM unstable region. Scanning $f\,\equiv\,\beta \cdot(1/T) \cdot (\partial T/\partial \psi_N)$ for 'high $q0$' (figure 14), we again find the mode is α-driven. These features were also found for other flux surfaces sampled across the core, and for $k_y\rho_\mathrm{r}\,=\,0.5$ (f scan for $(\psi_N\,=\,0.8, k_y\rho_\mathrm{r}\,=\,0.5)$ shown in figure 14). These observations support the conclusion that, although ideal ballooning stable, the dominant long-wavelength instability in these equilibria is a 'non-ideal' KBM.

Zoom In Zoom Out Reset image size

In this work, we have assessed the ballooning stability of three hypothetical ST equilibria, constructed with a commercial power-producing reactor in mind. In all three cases, the normalised pressure gradient α is sufficiently high to exceed the first stable region for ballooning modes; the equilibrium is either in the unstable region, or in the second stability window. In the negative triangularity case, the second stability window does not exist at positive magnetic shear, meaning the plasma is ideal ballooning unstable across the core; this coincides with strongly growing KBMs and indicating that this is unlikely to be a feasible equilibrium in the context of transport. We have considered the dependency of second stability access on shaping parameters, and found that second stability may be obtained either by making triangularity more positive, or reducing the elongation. These point to KBMs prohibiting negative triangularity as an option in commercial STs, unless they can be operated at very high field (which then introduces engineering challenges).

The equilibria with positive triangularity show better stability properties. In these, second stability access exists across the plasma and, provided the on-axis safety factor is not too low, the equilibrium occupies the second stability window across its full radius. Microinstability growth rates are correspondingly low. However, there remains a weakly growing instability with KBM-like properties, namely: pressure-driven, twisting parity, smoothly connecting to the 'ideal MHD' KBM, low parallel electric field and, in some cases, the modes are pervasive at long wavelength. We speculate that the KBM is destabilised by kinetic effects. Being weakly growing, these could feasibly be stabilised by flow shear or may impose a soft β limit on STs. Identifying these effects would be an interesting area of future work, and may be necessary if one wished to build a predictive model of ST H-mode pedestals in a similar manner to EPED [32].

The authors are grateful for productive discussions with M Anastopoulos-Tzanis, B Patel, C M Roach, B F McMillan and S Biggs-Fox. This work was supported by the Engineering and Physical Sciences Research Council (EP/L01663X/1) and (EP/R034737/1). This project was undertaken on the Viking Cluster, which is a high performance compute facility provided by the University of York. We are grateful for computational support from the University of York High Performance Computing service, Viking and the Research Computing team. This research used the GS2 version latest commit 142c78781c1bc4c06d234b32fee70eaf100d68e4 (8.1).

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

We use an approach similar to that used by Applegate et al [27]. We take the shearing rate as

$$\begin{equation} \omega_\mathrm{HB}\,=\,\frac{(\partial \psi/\partial r_\mathrm{minor})^2}{B} \frac{\partial^2 \Phi}{\partial \psi^2}, \end{equation} \tag{ A.1 }$$

where B is the magnetic field strength and Φ is the equilibrium electrostatic potential. We estimate Φ by taking the equilibrium force balance

$$\begin{equation} \frac{\mathbf{\nabla} p}{n_\mathrm{i} e}\,=\,\boldsymbol{V}\,\times\,\boldsymbol{B}\,+\,\boldsymbol{E}, \end{equation} \tag{ A.2 }$$

(where p is the ion pressure, $\eta_\mathrm{i}\,=\,(T_\mathrm{i}^{^{\prime}}n_\mathrm{i})/(T_\mathrm{i}n_\mathrm{i}^{^{\prime}})$, ^ʹ denoting a derivative with respect to ψ, V the ion flow velocity and E the equilibrium electric field), and setting $\boldsymbol{V}\,=\,0$. A little manipulation yields

$$\begin{align} \omega_\mathrm{HB}\,=\,\frac{(\partial \psi / \partial r_\mathrm{minor})^2}{B \cdot n_\mathrm{i} e} \left(\frac{\partial^2 p}{\partial \psi^2}\,+\,\frac{1}{p \cdot (1\,+\,\eta_\mathrm{i})}\left(\frac{\partial p}{\partial \psi}\right)^2 \right). \end{align} \tag{ A.3 }$$

We take B as the field strength on the surface at the location of the magnetic axis.

Since the data used for these simulations is publicly available, we only describe noteworthy parameters here. Their values were chosen to minimise computational cost without compromising accuracy of the simulations; we justify our choice by comparing a selection of $(\psi_N, k_y\rho_\mathrm{r}, \theta_0)$ for the three equilibria to higher-fidelity simulations. In particular, the difference in growth rate $\Delta \gamma\,=\,\gamma_{0}\,-\,\gamma_\mathrm{hf}$ (where γ₀ is the growth rate found with our parameters used and $\gamma_\mathrm{hf}$ the growth rate for a higher-fidelity simulation) is used as a figure of merit. In some cases $\gamma_{0,\mathrm{hf}}$ is low, resulting in large fractional errors, so here we quote the absolute error (in normalised units of $v_\mathrm{th,r}/r_\mathrm{LCFS}$).

The parameters governing the parallel coordinate θ are ${\textbf{nperiod}}\,=\,3$, ${\textbf{ntheta}}\,=\,192$ (although it was sometimes necessary to adjust ntheta owing to difficulties encountered in GS2 's gridgen module). The former specifies the extent of the simulation domain in ballooning space ($\theta\,=\,-5\pi$ to $\theta\,=\,5\pi$ in our case). This is insufficient to resolve extended microtearing structures which have been seen in high-β equilibria [22], but reasonable for KBMs. Preliminary research did not reveal any extended structures at long wavelength. Comparing ${\textbf{nperiod}}\,=\,3$ to ${\textbf{nperiod}}\,=\,5$ we found ${\textbf{max}}(\vert \Delta \gamma \vert)\,=\,4\,\times\,10^{-3}$ with rms value $\Delta \gamma_\mathrm{rms}\,\equiv\,\sqrt{\sum_N (\Delta \gamma)_N^2 / N}\,=\,4\,\times\,10^{-4}$. Comparing ${\textbf{ntheta}}\,=\,192$ to ${\textbf{ntheta}}\,=\,394$, we found $\Delta \gamma_\mathrm{rms}\,=\,5\,\times\,10^{-3}$ for $\psi_N\,\leqslant\,0.9$ and $\Delta \gamma_\mathrm{rms}\,=\,0.02$ for $\psi_N\,\gt\,0.9$, reflecting a difficulty in the calculation of geometric quantities in GS2 for strongly-shaped numerically-prescribed surfaces.

The equilibrium data contains density profiles $n(\psi_N)$ for electrons, the main ion species (with a mass of 2.5 times the proton mass and $Z\,=\,1$), helium ash and two impurity species: tungsten (with typical density $n_\mathrm{W}/n_\mathrm{e}\,\sim\,10^{-5}$) and xenon ($n_\mathrm{W}/n_\mathrm{e}\,\sim\,10^{-4}$) (both assumed fully ionised). Lacking simulated temperature profiles, it was assumed that $T(\psi_N)$ is identical for all species (the effect of fast ions may be an interesting area of future research.) Comparing simulations with the five kinetic species to the 'reduced' simulation with two kinetic species (in which $n_\mathrm{i}\,=\,n_\mathrm{e}$ to ensure quasi-neutrality), we found $\Delta \gamma_\mathrm{rms}\,=\,0.02$.

Collisions were ignored in the results presented here. Comparing collisionless to collisional simulations, we found $\Delta \gamma_\mathrm{rms}\,=\,0.01$ for $\psi_N\,\leqslant\,0.9$ and $\Delta \gamma_\mathrm{rms}\,=\,0.09$ for $\psi_N\,\gt\,0.9$ (consistent with the lower temperatures, and hence increased collisionality, near the edge).

The Miller parametrisation [31] describes each surface by nine dimensionless parameters: aspect ratio A, elongation κ, triangularity δ, the radial derivatives of major radius, elongation and triangularity s_κ , s_δ , $\partial R_0/\partial \psi_N$, the safety factor q, normalised magnetic shear $\hat{s}\,\equiv\,\partial q / \partial \psi_N$, and normalised pressure gradient $\alpha\,\equiv\,-\partial \beta / \partial \psi_N$. These specify the shape of the flux surface and its poloidal magnetic field:

$$\begin{align} &R\,=\,R_0\,+\,r\cos[\theta_M\,+\,x \sin \theta_M], \end{align} \tag{ C.1 }$$

$$\begin{align} &Z\,=\,\kappa r \sin \theta_M, \end{align} \tag{ C.2 }$$

$$\begin{align} & B_\mathrm{p}\,=\,\frac{(\partial \psi / \partial r) \kappa^{-1} R^{-1} \left[\sin^2(\theta_M\,+\,x \sin \theta_M)(1\,+\,x \cos \theta_M)^2\,+\,\kappa^2 \cos^2\theta_M \right]^{1/2}}{\cos(x \sin \theta_M)\,+\,(\partial R_0 / \partial r) \cos \theta_M\,+\,[s_\kappa\,-\,s_\delta \cos \theta_M\,+\,(1\,+\,s_\kappa)x \cos \theta_M]\sin \theta_M \sin(\theta_M\,+\,x \sin \theta_M)}, \\ \nonumber \end{align} \tag{ C.3 }$$

where R and Z are cylindrical coordinates (with the magnetic axis located at $Z\,=\,0$), $B_\mathrm{p}$ the poloidal magnetic field and θ_M a poloidal coordinate ranging from 0 to 2π.

We fit Miller parameters to the numerical SCENE flux surface data as follows. Firstly, the numerical values of R, Z and $B_\mathrm{p}$ for the flux surface of interest are extracted from the SCENE data. R and Z are used to construct $\theta_\mathrm{geo}$, a poloidal coordinate defined by $\theta_\mathrm{geo}\,=\,\arctan ((Z\,-\,Z_\mathrm{maxis})/$ $(R\,-\,R_\mathrm{maxis}))$, where $(R_\mathrm{maxis}, Z_\mathrm{maxis})$ is the location of the magnetic axis, and is taken from the SCENE equilibrium. Splines of $R(\theta_\mathrm{geo})$ and $Z(\theta_\mathrm{geo})$ are used to upsample R and Z, and from this calculate the values for A, r, κ, δ and construct Miller's poloidal coordinate by defining

$$\begin{equation} \theta_M\,\equiv\,\arcsin((Z\,-\,Z_\mathrm{maxis})/(\kappa r)) . \end{equation} \tag{ C.4 }$$

$B_\mathrm{p}(\theta_M)$ is then constructed and fitted to equation (C.3), with fitting parameters s_κ , s_δ , $\partial_r R_0$, $\partial \psi / \partial r$ (the latter determining q).

The above process yields an estimate of the Miller parameters describing a particular surface. However, it treats the parameters (A, r, κ, δ) on a different footing to (s_κ , s_δ , $\partial_r R_0$, q); calculating the former from R and Z does not guarantee an optimal fit of $B_\mathrm{p}$, or indeed of (R,Z) (since most of the R, Z points are ignored.) Therefore, we apply an additional step of simultaneously fitting ($R_M(\theta_M)$, $Z_M(\theta_M)$, $B_{\mathrm{p}M} (\theta_M)$) to ($R(\theta_M)$, $Z(\theta_M)$, $B_\mathrm{p} (\theta_M)$) allowing all eight parameters to vary freely, using the previously calculated values as initial guesses. This method ensures that $B_\mathrm{p}$ is treated on the same footing as R and Z, at the expense of allowing (A, r, κ, δ) to deviate from their geometric values.

Finally, the values of $\hat{s}$ and α are selected. Since $(R, Z, B_\mathrm{p})$ do not depend on $(\hat{s}, \alpha)$, these can be scanned independently in GS2 .

A caveat of the scheme presented above is that for reasons of convenience, we use the numerical equilibrium value of q rather than the value derived from $\partial_r R_0$.

Plots of the Miller parametrisation for several flux surfaces for each equilibrium are shown in figure C1. $R,Z,B_\mathrm{p}(\theta_M)$ and ($R,Z$) show good agreement between the SCENE data and Miller fit in the core, with slightly worse agreement towards the edge; in particular, the fitted $B_\mathrm{p}$ for the '−ve tri' has a tendency to artificially peak on the inboard side. Another systematic feature of the fitting is that the elongation tends to be underestimated in the positive triangularity cases and overestimated in the negative triangularity case, with the effect becoming more exaggerated at higher ψ_N .

Zoom In Zoom Out Reset image size