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.Footnote ³ 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.Footnote ⁴

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 Reference Boozer and Kuo-Petravic1981; Eriksson & Helander Reference Eriksson and Helander1994), 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 Reference Hirvijoki and Kurki-Suonio2012). 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.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 Reference Jucker2010), 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. Reference Dendy, Hastie, McClements and Martin1995). 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 Reference Patten2019). 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.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 Reference Patten2019) and (Machielsen, Graves & Cooper Reference Machielsen, Graves and Cooper2021). 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.Footnote ⁵ 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 (Reference Faustin2017, § 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.