Fast ion generation by combined RF-NBI heating in W7-X

Abstract

Good fast ion confinement in quasi-isodynamic stellarators still has to be demonstrated experimentally. In the absence of fusion alphas, auxiliary heating systems can be used to generate fast ions. For this purpose, W7-X has neutral beam injection and its newly installed ion cyclotron heating system. In practice, both systems can run simultaneously, the synergetic running of both is the focus of the current modelling development. With the upgrades to the SCENIC code described here, we find that, compared with pure neutral beam injection (NBI), the fast ion tail is significantly enhanced when radio frequency (RF) heating is enabled. But the effectiveness of heating and fast ion tail generation strongly depends on the plasma parameters and antenna frequency, which calls for further study, both numerical and experimental.

1. Introduction

Good alpha particle confinement is a requirement for a future stellarator reactor. This favourable property has yet to be demonstrated experimentally in W7-X. Due to the lack of fusion alphas, it is foreseen to generate fast ions using neutral beams, ion cyclotron resonance heating (ICRH) or both. It is these fast ions that function as ‘mock-up’ alpha particles for confinement tests. At an energy of 60 keV their Larmor radius divided by the minor radius is comparable to that of fusion alphas in a HELIAS reactor (assuming a scaled-up version of W7-X with 5 T magnetic field) (Drevlak et al. 2014). However, a perfect comparison will not be possible because particles in this energy range are still influenced by the radial electric field, while fusion alphas are not (Kolesnichenko et al. 2006; Faustin et al. 2016). NBI provides protons in this energy range, but only about a third are at full energy. Therefore, heating by RF waves may aid an exploration in this regime. As the name ICRH suggests, it relies on a wave–particle resonance

(1.1)\begin{equation} \omega = n\omega_c +k_\parallel v_\parallel, \end{equation}

where $\omega$ is the antenna frequency, $\omega _c$ the cyclotron frequency, $n$ the harmonic number, $k_\parallel$ the parallel wavenumber and $v_\parallel$ the parallel velocity. The cyclotron frequency depends on the particle position, so for a given $\omega, n, k_\parallel, v_\parallel$, the positions where (1.1) holds describe a surface (in three dimensions) at constant magnetic field strength. The wave–particle interaction is strongest on this surface. Therefore, a common approximation is to assume that diffusion in velocity space only occurs where (1.1) holds, also known as resonant diffusion. In a Monte Carlo orbit following code this is done by applying a kick in velocity, only when a test particle crosses the resonant surface. As some particles reside in a sweet spot, they continue to receive kicks from the wave for an extended time, thereby reaching energies considerably higher than the bulk temperature. Therefore, a side effect of the bulk heating by ICRH is the generation of a fast ion population, if the conditions are suitable. The diffusion in velocity space is constrained by

(1.2)\begin{equation} \frac{{\rm d}v_\parallel}{{\rm d}v_\perp} = \frac{v_\perp}{v_\parallel}\frac{\omega - n \omega_c}{n\omega_c}, \end{equation}

where $v_\perp$ is the perpendicular velocity (Kennel & Engelmann 1966). In pure RF heating schemes (so no NBI), the ion population that interacts with the wave typically starts with the same temperature as the background ions. So the numerator in (1.2) is nearly zero, and the kicks mainly increase the perpendicular velocity. Consequently, the high energy population that forms will have a lot of trapped ions, which are rapidly lost in the magnetic field of W7-X. One way to improve this is by using a so called three-ion scheme (Kazakov et al. 2015). In that case the power per ion is significantly enhanced, so higher energies can be obtained before the particles are lost. Another approach is to use neutral beams to seed a fast ion population, which is then further accelerated using RF heating. This is attractive for several reasons: the beams provide an extra 6.8 MW of power, in addition to a projected 1.5 MW of coupled RF power. The starting population is already above thermal energy, so it experiences less friction due to collisions as the ions are accelerated. And lastly, because the beam ions start with a considerable parallel velocity (and thus Doppler shift), the kicks are more isotropic, see (1.2). This results in a smaller trapped particle fraction compared with pure RF schemes.

Successful fast ion generation with combined RF-NBI heating was already demonstrated in stellarator experiments on Heliotron-J (Okada et al. 2016). Where a considerable enhancement of the fast ion tail was observed when ICRH was enabled. In order to aid experiments and extrapolate to future scenarios, extensive modelling of fast ions is a crucially important tool. In the present paper the code package SCENIC (Jucker 2010) is used to study such synergetic heating scenarios. It aims to self-consistently compute the magnetic equilibrium, the RF wave field and the fast ion distribution function. For this purpose the following codes are used respectively: VMEC (Hirshman & Whitson 1983) (or its anisotropic pressure version ANIMEC Cooper et al. 2009), LEMan (Popovich, Cooper & Villard 2006) and VENUS-LEVIS (Pfefferlé et al. 2014). These three are iterated until a converged self-consistent state is reached. However, modelling of combined RF-NBI poses additional challenges over pure RF and pure NBI scenarios, which will be discussed in detail in § 2. In short, the thermal ions cannot simply be ignored because they provide a source of fast ions due to RF heating. A pragmatic way to deal with this is to simulate both fast and thermal markers (test particles) of the minority species. Transport of thermal ions in many experiments far exceeds the predicted values by neo-classical transport. This ‘anomalous’ transport is thought to be caused by turbulence (Freidberg 2007, § 6.7). The true nature of this transport is very complex, but it is often well approximated as a diffusive process. The anomalous diffusion coefficient is typically of the order of 1 m$^2$ s$^{-1}$ (Wesson 2004, § 4.16). This type of diffusion is implemented as a Monte Carlo operator in VENUS-LEVIS, which has recently been made compatible with stellarator geometry (§ 2). Using this diffusion model, a steady state sets in after approximately 200 ms for pure NBI (§ 3). Enabling RF leads to a collisionless power deposition which strongly depends on the chosen parameters (§ 4). These parameters also have a major effect on the fast ion tail formation. If chosen correctly, the number of fast ions is enhanced considerably with RF-NBI (§ 5). The results are discussed in § 6 and summarised in § 7. Note that the goal of this study is not to perform an extensive parameter scan, but instead the focus is on identifying the key aspects that make a scheme suitable for fast ion generation. Hence, this paper describes the new feature (the radial diffusion model) in detail, and exploits it for an application where it is needed, in particular the RF-NBI scheme of W7-X.

2. Upgrade to VENUS-LEVIS

The fast ion population can be accurately modelled by sampling it with markers, and following their trajectories. In the case of pure NBI heating, the markers are born at high energy, and slow down on the background due to collisions. Since the main interest is in the fast ions, markers can be removed from the simulation as soon as they thermalise.¹ A balance is reached when the inflow of markers due to the beam equals the outflow due to losses and thermalisation.

Contrary to pure NBI, where the markers are born fast, and are mainly slowing down, RF provides a mechanism to speed up (thermal) markers. For this reason the thermalised markers are not removed from the simulation. In addition, in the absence of the real fuelling source (NBI), the particle balance has to be achieved in an ad hoc way. To ensure this balance in VENUS-LEVIS, the markers are ‘re-injected’ by sampling from the initial distribution function whenever they are lost, see figure 1.

Figure 1. Schematic overview of modelling in VENUS-LEVIS, inflow at the top, outflow at the bottom. The population in the grey box is tracked by VENUS-LEVIS. (a) pure NBI, (b) pure RF, (c) RF-NBI.

When modelling both heating methods simultaneously, the real particle source is known. But because RF is involved the thermal markers need to be included.² It has been found that this approach leads to an unphysical accumulation of thermalised beam ions, because the inflow far exceeds the losses in the simulation (Patten 2019).

This imbalance is expected as VENUS-LEVIS is designed for fast ion modelling. The guiding centre equations it solves are valid for both fast and thermal particles, however, the collision operator is incomplete for thermal ions. The transport experienced by these thermal ions, which includes turbulence, is complicated. In order to avoid modelling this complex physics, a radial diffusion model has been incorporated.

Note also that VENUS-LEVIS does not include self-collisions. Meaning, markers only experience a friction due to collisions with plasma species other than their own. This is justified as long as the marker species is a minority, where ‘species’ in this context refers to distinct populations which can be of the same particle type. For example, a hydrogen plasma with hydrogen beams will have electrons, thermal hydrogen and a high energy tail population which can be modelled as its own hydrogen ‘species’. VENUS-LEVIS would then be employed to only resolve this fast H population, which is indeed a small fraction of the total number of plasma particles, so neglecting self-collisions is justified. If thermal markers are not removed from the simulation this leads to an underestimation of the collisional drag. However, as will be shown later on, the particle fraction of markers is still only a few per cent in the examined cases.

2.1. The RF-NBI model, step by step

Typically, when simulating NBI using a Monte Carlo orbit tracer such as VENUS-LEVIS, all markers are spawned in the plasma at the start of the simulation, contrary to a continuous injection. This is allowed since every marker moves independently, and the procedure provides sufficient information to reconstruct the steady-state beam profiles.³ They are then traced until they are either lost, or thermalised. Along the way the state of the markers is saved a number of times (e.g. every millisecond). From these data the beam profiles can be reconstructed.

To make the plasma quasi-neutral the beam species density profile is then subtracted from that of the background ion species’. These profiles are then passed on to LEMan. One limitation of the model is that the steady-state concentration of the beam is unknown at the start of the simulation, so there is no way to ensure the final beam concentration matches that of an experiment. Measurements provide the plasma profiles of all the species, and the injection rate of the beam is fixed too. Which leaves just one free parameter, the radial diffusion rate. This coefficient can be varied to find a best fit of the beam species concentration to that of the experiment (such a scan is not performed here).

The problem with the single injection (all markers spawned at once) is that it ends when all markers have been lost, thermalised or a given time is exceeded. The last state of the markers is therefore empty or nearly empty. For this reason a different approach is used when RF is also involved: markers are continuously injected and lost, and the simulation is run until particle balance sets in. Once this state is reached the profiles can be computed, LEMan is then run and VENUS-LEVIS can resume its run using the last marker state. This contains both freshly injected and slowed down beam ions.

However, continuous injection has some challenges as well: the final number of confined markers (in steady state) is not yet known at the start of the simulation. Also the necessary simulated time needed to reach equilibrium is unknown. This is incompatible with VENUS-LEVIS because it is designed to model a fixed number of markers for a given time. A workaround is to introduce the concept of ‘cycles’. A cycle represents one run of VENUS-LEVIS, for a given fixed time (e.g. 100 ms). Instead of generating a marker list for just one cycle, it is made to include all future cycles as well. If after completing cycle $j$ the simulation is not yet deemed converged, cycle $j+1$ can be launched by simply shifting the starting index in the marker list, which determines the freshly injected markers during cycle $j+1$. After some time the simulation contains markers from cycles $j+1,j,j-1,\ldots$, but most of the old ones will be lost (inactive). This leads to a larger, but still easily manageable memory footprint. In addition, the computational overhead associated with keeping track of the extra inactive markers is negligible.⁴

For simplicity, in this work an isotropic pressure profile is used to compute the magnetic equilibrium, so using VMEC, not ANIMEC. In addition, it is run only once, instead of including it in the iterative loop. In summary, the modelling procedure is as follows:

1. (i) compute magnetic equilibrium using VMEC;

2. (ii) compute NBI deposition (i.e. where the beam particles ionise) with VENUS-LEVIS;

3. (iii) run VENUS-LEVIS (for multiple cycles) until steady state is reached;

4. (iv) generate plasma profiles of the beam species and adjust the background density to make it quasi-neutral;

5. (v) run LEMan to generate the wave field;

6. (vi) repeat step 3, but now with ICRH enabled; and

7. (vii) repeat steps 4-6 until converged.

2.2. Diffusion operator

In order to include the effect of collisions on markers, VENUS-LEVIS uses Monte Carlo operators. These are assumed to be independent. So the effects of pitch-angle scattering, energy kicks, interaction with the wave field, etc. can be treated separately. For related work see (Boozer & Kuo-Petravic 1981; Eriksson & Helander 1994), however, note that VENUS-LEVIS is not orbit averaged. VENUS-LEVIS also has the option to enable the effect of radial diffusion on markers, which too is implemented as a Monte Carlo operator. Now adopting the same diffusion operator as used by ASCOT, which is valid for any flux surface shape, including stellarators (Hirvijoki & Kurki-Suonio 2012). Its derivation can be summarised as follows:

the term responsible for radial diffusion can be written as

(2.1)\begin{equation} \left( \frac{\partial f}{\partial t}\right)_{d} = \boldsymbol{\nabla} \boldsymbol{\cdot} \left( D \boldsymbol{\nabla} f \right) . \end{equation}

With $f$ the distribution function of the marker species and $D$ the diffusion coefficient. The subscript $d$ indicates that this is just the contribution of diffusion. Defining a generic moment as follows:

(2.2)\begin{equation} \left\langle h\right\rangle = \int_{\varOmega} f h \,{\rm d}V, \end{equation}

where the integration is over the entire plasma volume $\varOmega$. The time derivative, due to diffusion alone, is

(2.3)\begin{equation} \left(\frac{{\rm d}\left\langle h\right\rangle}{{\rm d}t}\right)_{d} = \int_{\varOmega} \left(\frac{\partial f}{\partial t}\right)_{d} h \,{\rm d}V = \int_{\varOmega} h \boldsymbol{\nabla} \boldsymbol{\cdot} \left( D \boldsymbol{\nabla} f \right) \,{\rm d}V . \end{equation}

Integrating by parts (twice) gives

(2.4)\begin{equation} \left(\frac{{\rm d}\left\langle h\right\rangle}{{\rm d}t}\right)_{d} = \int_{\varOmega} f \boldsymbol{\nabla}\boldsymbol{\cdot} (D\boldsymbol{\nabla} h)\,{\rm d}V + \int_{\partial \varOmega} D\left(h\boldsymbol{\nabla} f - f\boldsymbol{\nabla} h\right) \boldsymbol{\cdot} \boldsymbol{n} \,{\rm d}S = \left\langle\boldsymbol{\nabla} \boldsymbol{\cdot} (D\boldsymbol{\nabla} h)\right\rangle . \end{equation}

The surface integral appears due to the divergence theorem, but $f$ is assumed to vanish at the plasma boundary, so this surface term disappears. Note that this assumption will be easily satisfied later on as a delta function will be used for $f$. Now the following coordinate system $(u^1,u^2,u^3)$ with $u^1$ the radial flux label is used. Substitution of $h=u^1$ or $h=(u^1)^2$ in (2.4) yields, respectively, the following in curvilinear coordinates:

(2.5)\begin{gather} \left(\frac{{\rm d}\left\langle u^1\right\rangle}{{\rm d}t}\right)_{d} = \left\langle\frac{1}{\sqrt{g}}\frac{\partial}{\partial u^i}\left(\sqrt{g}D g^{1i}\right)\right\rangle, \end{gather}

(2.6)\begin{gather}\left(\frac{{\rm d}\left\langle (u^1)^2\right\rangle}{{\rm d}t}\right)_{d} = 2\left\langle\frac{1}{\sqrt{g}}\frac{\partial}{\partial u^i}\left(\sqrt{g}D u^1g^{1i}\right)\right\rangle , \end{gather}

where the Einstein summation convention is applied, $g^{ij}$ is the metric tensor and $\sqrt {g}$ the Jacobian determinant. The time evolution of the variance $\sigma ^2$ can be written as

(2.7)\begin{equation} \left(\frac{{\rm d}\sigma^2}{{\rm d}t}\right)_{d} = \left(\frac{{\rm d}}{{\rm d}t} \left[ \left\langle (u^1)^2\right\rangle- \left\langle u^1\right\rangle^2 \right]\right)_{d} = \left(\frac{{\rm d}\left\langle (u^1)^2\right\rangle}{{\rm d}t}\right)_{d} -2\left\langle u^1 \right\rangle\left(\frac{{\rm d}\left\langle u^1\right\rangle}{{\rm d}t}\right)_{d} . \end{equation}

Initially, a particle is located in just one point $u^1(0)$, so its distribution function is $\delta (u^1-u^1(0))$. Note that this Dirac delta is weighted such that $\left \langle 1\right \rangle = 1$. Thus, for $t=0$ the volume in (2.6) can be eliminated

(2.8)\begin{equation} \left. \begin{aligned} \left.\left(\frac{{\rm d}\left\langle u^1\right\rangle}{{\rm d}t}\right)_{d}\right|_{t=0} & = \frac{1}{\sqrt{g}}\frac{\partial}{\partial u^i}\left(\sqrt{g}D g^{1i}\right),\\ \left.\left(\frac{{\rm d}\left\langle (u^1)^2\right\rangle}{{\rm d}t}\right)_{d}\right|_{t=0} & = \frac{2}{\sqrt{g}}\left( u^1(0) \frac{\partial}{\partial u^i}\left(\sqrt{g}D g^{1i}\right) +\sqrt{g}Dg^{11}\right),\\ \left.\left(\frac{{\rm d}\sigma^2}{{\rm d}t}\right)_{d} \right|_{t=0} & = 2Dg^{11} , \end{aligned} \right\} \end{equation}

where the right-hand side is also evaluated at $t=0$. From this the stochastic radial kick is constructed

(2.9)\begin{align} \Delta u^1 & = \Delta t \left.\left(\frac{{\rm d}\left\langle u^1\right\rangle}{{\rm d}t}\right)_{d}\right|_{t=0} + R \sqrt{\Delta t \left.\left(\frac{{\rm d}\sigma^2}{{\rm d}t}\right)_{d} \right|_{t=0} }\nonumber\\ & = \frac{1}{\sqrt{g}}\frac{\partial}{\partial u^i}\left(\sqrt{g}D g^{1i}\right) \Delta t + R \sqrt{2Dg^{11}\Delta t}, \end{align}

with $R$ a randomly selected plus/minus sign. For this derivation no assumption has been made on the magnetic equilibrium, nor on the spatial dependence of $D$. But for simplicity, a homogeneous diffusion coefficient has been chosen in this work. Along the same lines as laid out above, a Monte Carlo operator for diffusion in the angular directions $u^2,u^3$ can also be derived. However, this type of diffusive transport in the angular variables is unimportant compared with the normal transport channels in those directions, hence it is omitted in VENUS-LEVIS.

The impact of diffusion is studied in a pure NBI scenario, see figure 2. Here, thermal markers are defined as having an energy less than 1.5 times the half-radius electron temperature, i.e. $E<1.5 T_e(s=\rho ^2=0.5)$, and ‘fast’ is anything above that. The diffusion mostly effects low energy markers, but even in the range 20–50 keV a small but noticeable effect is seen. Above 50 keV the effects are negligible, as seen by the convergence of the lines around this energy in figure 2. To explain this effect, note the following: in the absence of RF there is no mechanism for markers to speed up (apart from the random collisions that lead to some markers increasing their energy). Thus the lower energy markers have been in the plasma for longer, during which diffusion has had a chance to kick out of some of them. For the freshly injected markers on the other hand, diffusion has not had the time. Hence the losses accumulate at the low energy range.

Figure 2. Comparison of the energy distribution for three cases: no diffusion, diffusion on thermal markers only or diffusion on all markers. The data represent the NBI slowing down distribution using base case parameters after 200 ms, normalised such that the integral over energy gives the total number of H ions in the plasma.

Besides the direct effect of diffusion on fast ions, there is also an indirect contribution, namely, the thermal concentration affects the dielectric properties of the plasma. And so in scenarios involving RF heating, the wave propagation and polarisation are affected. This in turn modifies the RF power per ion. As was already touched on in § 2, VENUS-LEVIS is designed for modelling fast ions, not thermal ions. This deficiency is remedied by adding a diffusion operator. However, to be consistent with the assumption that collisions are accurately modelled for fast ions, the diffusion is only applied to the thermal ions in the rest of this work. To be explicit, the following expression is adopted for the diffusion coefficient:

(2.10)\begin{equation} D(E) = \begin{cases} D_0, & E<1.5T_e(s=0.5),\\ 0, & E \ge 1.5T_e(s=0.5) . \end{cases} \end{equation}

3. Neutral beam injection

Using the newly upgraded version of the code, the capabilities of fast ion generation in W7-X are assessed for several NBI and RF-NBI test cases, and their parameters are listed in table 1. The plasma consists of a deuterium background, with H beams. Balanced injection is used (injectors $3,4,7,8$), which results in a total inflow of approximately $1.2 \times 10^{21}$ protons per second. This amounts to a total injected power of approximately 6.84 MW. See Rust et al. (2011) for more information on the W7-X NBI system. Note that D can be swapped with He-4 with little change to the wave physics: in cold plasma two species with the same charge to mass ratio, and the same plasma frequency, make equal contributions to the dielectric tensor. Because of quasi-neutrality, the density of He-4 would be half that of D, therefore their plasma frequency is the same. However, it should be mentioned that the contributions are not exactly the same in hot plasma due to the difference in thermal velocity.

Table 1. Description of the various test cases used. The density $N$ and temperature $T$ are on-axis values. For the profiles used see Appendix A. The magnetic equilibrium is generated using base case parameters.

3.1. Steady-state base case

After two cycles of 100 ms a steady state is reached, meaning loss rate $\approx$ birth rate, see figure 3. Note that the birth rate is less than the injection rate because of shine through. After one cycle of 100 ms, steady state is almost reached, but the amount of new markers still exceeds the number of lost ones. In the second cycle, 1 990 656 are born, and 1 918 915 are lost, demonstrating that a balance is reached. Note that perfect balance is not required, as multiple iterations of VENUS-LEVIS will be run in more complete simulations. After reaching this balance, moments of the slowing-down distribution function can be computed, such as the density as shown in figure 4(a), and the anisotropic pressure. As described in (Jucker 2010), the parallel and perpendicular temperatures can then be derived from this pressure, see figure 4(b). The central hydrogen concentration is approximately 7.3 %. To make the plasma quasi-neutral, this profile is subtracted from that of deuterium.

Figure 3. Loss rate of beam particles vs time. Using base case parameters, table 1.

Figure 4. Hydrogen plasma profiles vs radial flux label $\rho$. The coordinate $\rho$ is defined as the square root of the normalised toroidal flux. Base case parameters used. (a) Density, (b) temperature.

In order to use the marker population in the computation of the wave field, it is first split into a thermal background (isotropic Maxwellian) and a fast part (modified bi-Maxwellian Dendy et al. 1995). This splitting is unrelated to that used for the diffusion; it is chosen to find a reasonable fit to the marker population using a combination of the two aforementioned distribution functions. In this particular case ‘fast’ was defined as having $E>6T_e(s)$, so using the local electron temperature. This criterion has also been used to generate figure 4. The distribution of particle pitch $\lambda = v_\parallel /v$ is seen in figure 5.

Figure 5. Particle pitch distribution for the slowing down (SD) distribution after 2 cycles (200 ms), normalised such that the integral over pitch gives the total number of H ions in the plasma. In red, with corresponding $y$-axis on the right, the deposition of particles for one whole cycle. The symmetry is a result of using balanced beams.

Injectors 4,8 are more tangential, creating peaks at $\lambda = \pm 0.46$, while 3,7 are more radial, resulting in peaks at $\lambda = \pm 0.30$ (figure 5). The exact pitch of the peaks in figure 5 depends on the magnetic field. After some time has passed, these peaks have become smoothed due to pitch-angle scattering, which is most effective at low energy. This explains why the ‘low’ energy ($<20$ keV) part of the distribution is nearly isotropic in pitch. Eventually the distribution reaches a steady state which is well described by a Maxwellian plus a high energy population. This low density, high temperature population will never become isotropic no matter the waiting time, because the particle source itself is anisotropic.

Note that the initial NBI population consists overwhelmingly of passing particles, see figure 6, which shows the distribution of fresh beam ions in terms of their reflection field $B_{ref}=E/\mu$. Where $E$ is the kinetic energy, and $\mu = m v_\perp ^2/(2B)$ the magnetic moment. A particle is classified as passing when $B_{ref} > \max (B)$, and trapped otherwise, where the maximum is to be taken along its full trajectory. In other words, reflection would occur at $B=B_{ref}$, but if $B$ is always less than $B_{ref}$ this bounce point is never encountered, hence it is passing. Figure 4(b) shows that the perpendicular temperature of the fast population exceeds its parallel temperature. This does not contradict the finding that most particles are passing; indeed $v_\perp >|v_\parallel |$ on average, which results in the anisotropy, yet the ions are mostly passing because of the injection sites. Which are at high magnetic field strength, meaning small $\mu$, and thus large $B_{ref}$.

Figure 6. Distribution of $B_{ref}=E/\mu$ values for the freshly injected beam ions, normalised so that the integral over $B_{ref}$ gives the total ionisation rate. For comparison, the value of the magnetic field strength on axis is approximately $2.43$ T in front of the antenna.

Note that in all of the runs in this work the same high mirror (KJM) configuration has been used. The reason for this choice is because the high mirror scenario results in better core NBI confinement than in low/standard mirror (Patten 2019). To clarify, if the empirical goal is simply RF bulk heating, one may be better off using the low mirror because the magnetic field strength varies less in toroidal direction. Suppose on-axis heating is used in front of the antenna, then ions will still resonate near the axis at all toroidal angles. So the share of thermal ions in the core which can interact with the wave is maximised. In the high mirror, on the other hand, the resonant surface even leaves the plasma in the triangular cross-section. However, remember that bulk heating is not the aim of this work. The goal is to maximise fast ion generation for advanced studies, such as testing confinement, exciting Alfvén waves and transport suppression, which can only be achieved if the NBI ions are well confined, and if these ions can be further accelerated by RF.

4. The ICRF wave

4.1. Doppler shift

The wave–particle resonance is given by (1.1). In this work the antenna frequency is chosen such that the hydrogen resonates at the fundamental cyclotron frequency, i.e. $n=1$. In addition electrons will resonate at $n=0$ and deuterium at $n=2$. The latter being a hot plasma effect, which LEMan can now simulate due to its recent upgrades, detailed in (Patten 2019) and (Machielsen, Graves & Cooper 2021). Including hot plasma effects is important in order to get an accurate estimate of the power fractions of the various species.

The thermal ions will resonate near the cold resonance, $\omega =\omega _c$, but freshly injected beam ions experience a considerable Doppler shift, as illustrated in figure 7. As the wave enters the plasma from the low-field side it will first encounter a resonance with the beam, then the cold resonance and then another resonance with the beam. This Doppler shift itself cannot be controlled in any feasible way. On the other hand, an operator can adjust the antenna frequency up or down (i.e. ‘shifting’ it) so that one of the beam resonances goes through the core.⁵ In the simulation the antenna frequency is set indirectly by choosing a resonant magnetic field $B_{\textrm {res}} = m \omega /q$.

Figure 7. Simple cartoon to illustrate where the resonance is. The thermal ions resonate near the cold resonance $(\omega =\omega _c)$, indicated in blue. The red dashed lines indicate where the beam ions resonate, and the $\times$ highlights the magnetic axis. In panel (a) the antenna frequency is shifted up, in (b) it is shifted down.

Choosing the optimal antenna frequency for fast ion generation is not obvious. Ideally, the beam ions should resonate in the core, where the confinement is best. However, the magnetic field strength on axis depends on toroidal angle. Secondly, there is a spread of parallel velocities figure 5, and $k_\parallel$ is also not uniform in space. The optimal antenna frequency should follow from an extensive simulation (or experimental) scan, but a rough estimate can be obtained by assuming the following values: $\lambda =\pm 0.4$ and full beam energy (55 keV), hence $v=3.2\times 10^6\ \textrm {m}\ \textrm {s}^{-1}$. The parallel wavenumber is computed by LEMan, the on-axis value near the antenna is found to be $k_\parallel =14\ \textrm {m}^{-1}$. This results in a frequency shift of approximately $2.9$ MHz (up or down), corresponding to a shift in resonant magnetic field of $0.19$ T. For the chosen magnetic configuration the field strength on axis in front of the antenna is approximately $2.43$ T, meaning $B_{\textrm {res}} = (2.43 \pm 0.19)$ T, and the antenna frequency should be $\omega /(2{\rm \pi} ) = (37.1 \pm 2.9 )$ MHz. As a starting point, a value of $B_{\textrm {res}}=2.5$ T is used unless specified otherwise. In order to evaluate the dielectric tensor, LEMan assumes that locally in each point a single parallel wavenumber dominates (or set of wavenumbers). To replicate dipole phasing, both $k_\parallel, -k_\parallel$ are used. See Faustin (2017, § 4.5.1) for a description of the antenna model. The LEMan simulations in this work used a radial resolution of 200 elements and $30\times 300$ (poloidal $\times$ toroidal) Fourier modes.

4.2. Wave field

Several different plasma scenarios are compared, and their parameters summarised in table 1. With respect to the base case, an increased density leads to a lower power fraction of H. This is because the injection rate is kept the same, while the background density has increased, i.e. there is a lower H density fraction. It should be mentioned that the lower minority concentration leads to an improved wave polarisation at the cyclotron resonance. However, the H concentration in the base case was not yet high enough to significantly deteriorate the wave polarisation, therefore lowering the concentration results in less minority heating. An increased diffusion rate leads to a lower thermal concentration, which results in a somewhat increased fast H power fraction. If instead the temperature is raised, the damping on electrons and D is more effective. This is not surprising, because D only absorbs via second harmonic heating, which is a finite Larmor radius (FLR) effect. In addition, the increased temperature results in a broader resonance layer (Doppler broadening due to larger $T_\parallel$), which could also explain part of the increased power fraction of D. Electrons also absorb more power, in part because transit time magnetic pumping plays a bigger role as temperature is increased.

The power fractions in figure 8 only indicate the direct wave heating. Subsequently, the power is redistributed over the plasma species via collisions. This collisional power is not computed by LEMan, but is obtained by the follow-up VENUS-LEVIS run, see e.g. figure 15.

Figure 8. Fraction of the RF power directly deposited on the various plasma species. A distinction is made between the thermal (Maxwellian) hydrogen and the fast hydrogen.

It can be seen in figure 9 that the RF power deposition is not peaked on axis. There are multiple reasons for this: the ion cyclotron resonance happens near a surface where $\omega = \omega _c$. This surface gets close to the magnetic axis in front of the antenna, but departs from the magnetic axis at other toroidal angles, meaning that, at those other angles, the power is mainly deposited off axis. Secondly, the aforementioned surface cuts through the plasma, and resonance therefore does not just occur on axis but also at other radial locations (all the way to the plasma edge). And lastly, as can be seen from figure 10, figure 11, the wave amplitude is not at all homogeneous. It would be a coincidence if a peak in $|E^+|$ was exactly aligned with the magnetic axis.

Figure 9. Flux surface averaged power density vs radial position. (a) Base case, (b) density $\times 3$.

Figure 10. A plot of the amplitudes of the left handed (a) and right handed (b) polarised component of the electric field of the wave in front of the antenna (base case). Note, $E^+$ is mainly responsible for heating H. (arbitrary units used)

Figure 11. Plot of the left and right handed polarised components of the wave field (density x3).

As can be seen from the electric field in figure 12, figure 13 the wave propagates less far in the high density case (stronger absorption). The heating of the plasma is therefore also more localised near the antenna.

Figure 12. Cross-section of the plasma, top-down view, colour indicates $|E^+|$ in arbitrary units. The hot spot is located just in front of the antenna.

Figure 13. Same, but for density x3 case.

5. Fast ion distribution

In order to compute the effect of the ICRF wave on the markers, VENUS-LEVIS needs the correct wave amplitude. To this end VENUS-LEVIS is run first for a number of time steps to measure the absorbed wave power. The field is then rescaled such that this power matches the expected coupled power (assumed to be 1.5 MW here), multiplied by the power fraction of the H species. The simulation is restarted after this step. Later in the simulation the field is rescaled at regular time intervals,⁶ but since these are typically minor adjustments the simulation is not reset to the state of the last time interval.

The addition of RF provides a clear improvement in the fast ion tail over pure NBI, see figure 14 (and appendix A for all cases). However, the effectiveness strongly relies on the plasma parameters. In the high density case the collisionality is too high, preventing a significant tail from forming. The opposite is true for the case with increased temperature (reduced collisionality).

Figure 14. Comparison of the H particle energy distributions, normalised such that the integral (i.e. summing over bin values times bin widths) gives the total number of H particles in the whole plasma. The vertical lines indicate the injection energy (full, half, third).

The binning in figure 14 has been truncated at 150 keV, where all distributions disappear into noise. This can be improved in the future by using more markers, which reduces the noise and higher energies will be resolved. However, this would require excessive marker counts. Alternatively, the effect of noise can be mitigated by using markers with adaptive weights (Schuster et al. 2021). The basic idea is to continuously split and delete markers to increase the resolution in the high energy range. Another way to improve the energy resolution is by replacing the thermal population by an analytical particle source term, so only fast ions need to be modelled. In this method the thermal beam ion population is assumed to be Maxwellian, at some given concentration and temperature known from experiments. Using this distribution and the wave field, an expression for the fast ion particle source can be computed in principle. In other words, the left most scheme in figure 1 would be used, instead of that on the right. Except that it will be supplemented by an influx of markers computed based on the RF heating of the bulk population. As the number of thermal markers vastly outnumbers the fast ones, this approach would potentially lead to the biggest improvement. However, deriving such a source term is far from trivial. To be clear, these three approaches are not mutually exclusive, a combination of them can be used.

VENUS-LEVIS also computes the collisional power on the background species (electrons and deuterons), figure 15. This collisional power is dominated by the beam, meaning that adding RF does not significantly affect this profile. Also observe the spike in electron power near the plasma edge. This is because the beam temperature is much larger than that of the background plasma. A more detailed explanation is found in Appendix B. This spike is an artefact emerging from the choice of plasma profiles, however, it constitutes only a minor fraction of the total power as seen from figure 15(b). Hence its effect on the core physics is negligible.

Figure 15. (a) Collisional power density vs radial position (base case). For comparison, the NBI-only case (second cycle) has been plotted too. (b) Total (cumulative) collisional power within given flux surface.

So far the antenna frequency has been set to 38.1 MHz for simplicity, however, it can be varied in the range of 25–38 MHz (Ongena et al. 2020). The antenna frequency dictates $B_{\textrm {res}}$, which strongly affects which ions resonate and where, see (1.1). The more $B_{\textrm {res}}$ deviates from the on-axis value (2.43 T), the more the resonant surfaces are displaced. Which means that only particles with larger $|v_\parallel |$ can still resonate in the core where the density is high. As a result, the larger the shift, the more power goes to the fast H vs the thermal H, figure 16. This is a trade-off, the antenna frequency should be shifted to optimally accelerate fast ions present in the core even further. But shifting it too far means that the slower ions do not get a chance to reach this optimal level. This was also investigated in Patten (2019, § 5.2.3). Note that $B_{\textrm {res}}=2.7$ T requires an antenna frequency of 41.2 MHz, which is outside of the accessible range, so this is a more academic test case. In a real experiment the coil currents can be reduced, while maintaining the same frequency shift, e.g. $B_{\textrm {res}}=2.6$ T, while having $B=2.33$ T on axis in front of the antenna. Another option would be to shift the antenna frequency down instead of up. This is yet another variation that should be investigated. Note that the resonant surface at $\omega =\omega _c$ is surrounded by two Doppler shifted surfaces, as illustrated in figure 7. Since the injection is from the outboard side, the wave first resonates with the beam ions, then with the thermal H and then with the beam again. By shifting the frequency down, the beam ions in the core encounter the wave last, which may lead to less acceleration of beam ions. Also note that, by shifting $\omega$ too much (up or down), only two of the three resonances remain in the plasma (near the antenna).

Figure 16. H particle energy distributions for the base case, but repeated with different resonant magnetic field strengths.

The particle pitch is also investigated, and found to be anisotropic, see figure 17. It is mostly similar to the pitch in figure 5, except with enhanced energy. Particles with low pitch do exist, but they are lost rapidly, see figure 18. The losses can also be binned in terms of toroidal Boozer angle, see figure 19. It is nearly, but not fully, fivefold periodic.

Figure 17. Fast ions binned in pitch and energy. The colour indicates $\log _{10}$ of the total number of particles per bin. Base case used; (a) $B_{\textrm {res}}=2.5$ T, (b) $B_{\textrm {res}}=2.7$ T.

Figure 18. Fast ion losses binned in pitch and energy, for the base case. The colour of a bin indicates $\log _{10}$ of the lost particles per second; (a) $B_{\textrm {res}}=2.5$ T, (b) $B_{\textrm {res}}=2.7$ T.

Figure 19. Particle loss rate vs toroidal Boozer angle $\phi$, normalised such that the integral over $\phi$ gives the total loss rate in particles/second. The black vertical line indicates the position of the antenna, and the green lines the positions where the NBI deposition is strongest.

6. Discussion

The current modelling approach presents a number of notable strengths and weaknesses, summarised below.

Strong points:

1. (i) Sampling particles and tracing their orbits is an effective method to deal with the high-dimensional problem associated with the evolution of the distribution function.

2. (ii) The model finds a self-consistent solution for the fast ion distribution function and RF wave field in the inherently three-dimensional plasma of a stellarator.

3. (iii) VENUS-LEVIS includes finite orbit width effects.

4. (iv) The wave code uses a localised antenna model in three dimensions.

5. (v) RF heating on all species is computed, including FLR effects. Based on the H power fraction, the electric field amplitude is rescaled in VENUS-LEVIS.

Weak points:

1. (i) Geometric limitations: the antenna model is tied to the magnetic equilibrium. This equilibrium is computed by VMEC, so it assumes nested magnetic flux surfaces, and runs up to the LCFS only. Moreover, VMEC is not included in the self-consistent loop, meaning $\boldsymbol {B}_0$ is assumed to be fixed. Lastly, the NBI injection geometry is simplified compared with that in ASCOT/BEAMS3D.

2. (ii) Antenna model is an approximation of dipole phasing.

3. (iii) A steady-state radial electric field is not included.

4. (iv) The diffusion model requires an additional input, namely, the diffusion coefficient. However, this ad hoc parameter can be used to match the minority concentration to that measured in an experiment.

5. (v) Mode conversion is not included in the wave code.

General comments:

VENUS-LEVIS does not use any type of bounce averaging, it solves the full four-dimensional guiding centre equations. This means that the model includes finite orbit width effects, which are important for getting a realistic estimate of the particle losses and collisional power deposition profile. Coulomb collisions are included as well in the form of a Monte Carlo operator.

Regarding the full-wave code, the antenna is implemented with a given radial, poloidal and toroidal extension (‘localised in three dimensions’). Additionally, the dielectric tensor in the hot model makes no approximation on the size of the Larmor radius compared with the wavelength of the electric field.

Compared with the previous RF-NBI work with SCENIC (Patten 2019; Patten et al. 2020), several features have been enhanced/added: the use of cycles to resume NBI runs, and the radial diffusion model in order to reach particle balance. In addition, the plasma profiles have been updated and the aforementioned hot plasma model is enabled. Additionally, the NBI deposition in the new version of SCENIC is more complete, taking into account the beam focussing.

Some further improvements could be made in the future. Most notably, on the subject of the geometric effects. It should be mentioned that this is not a fundamental limitation, in the sense that the code could be modified to take in a more general magnetic geometry that extends to the metal wall. Fast ion hot spots on the wall depend rather sensitively on the magnetic configuration (Äkäslompolo et al. 2018). In order to get good agreement with experiment on this front, these limitations should be addressed within the SCENIC package. Recently, however, strides have already been made to compute the fast ion wall loads for RF minority heating. This was done by taking the markers from VENUS-LEVIS that crossed the LCFS, and importing them into ASCOT. These were then traced to the actual metal wall (Slaby et al. 2022).

A steady-state radial electric field arises from the ambipolarity condition. However, for simplicity it has not been included in these runs. The radial electric field (ion root) typically helps to confine fast ions because of the resulting $\boldsymbol {E}\times \boldsymbol {B}$ drift, see e.g. Kolesnichenko et al. (2006); Faustin et al. (2016). Therefore, by omitting this effect the fast ion population is likely underestimated. Especially for trapped ions in the NBI range of energies, which might otherwise be heated by the RF wave. It should be noted that the radial electric field depends on the plasma profiles, so it will be different for every case in table 1.

It should furthermore be noted that in these simulations the placement of the antenna differs somewhat with respect to the real antenna position. In reality it is offset in toroidal direction with respect to the bean section, by approximately $8^{\circ }$, away from the injection points of the beams. In the model it has been offset $8^{\circ }$ in the other direction, i.e. closer to the beam injection sites.

7. Summary

The radial diffusion operator in VENUS-LEVIS has been generalised to a stellarator geometry. This, together with some other code improvements, enables the marker count to saturate. That is, the particle inflow due to NBI matches the outflow due to losses. Several different scenarios involving synergetic RF-NBI heating have been investigated. The best performance (in terms of fast ion generation) was obtained with a low collisionality plasma, i.e. low background density and high temperature. In addition, for the RF-NBI scheme it is crucial to have the right amount of antenna frequency shift. However, given the vastness of the parameter space, further study is needed to properly optimise this scenario for W7-X.