2. General remarks on ICRH modelling

Ion cyclotron resonance heating involves processes occurring on very different time scales. The most rapid time scale is that of the RF dynamics. Since the RF frequency is in the MHz range, all variations of the background against which these effects occur can safely be neglected. For that reason the wave equation is typically solved as a driven problem, assuming the densities and temperatures of the background species are known and prescribed. Power transferred from the waves to the particles is further redistributed among the plasma constituents through Coulomb collisional interaction. This process – happening on time scales that are a few orders of magnitude larger than the RF time scale – is described by the Fokker–Planck equation. Energy and particle transport across magnetic surfaces is commonly assumed to occur on time scales that are longer than those on which slowing down processes occur. When this is not the case, the Fokker–Planck equation needs to be solved accounting for temporal evolution of the distribution functions. Solving the wave equation yields the RF perturbed electric field, while solving the Fokker–Planck equation allows us to determine the slowly varying distribution function $F_o$. Provided the wave–particle interaction is modelled self-consistently, the power density as seen by the wave equation is identical to that of the Fokker–Planck equation,

(2.1)\begin{align} P_{\textrm{abs}} &= \frac{\partial}{\partial t}\left[ \int \textrm{d}\boldsymbol{v}\,\textrm{d}\boldsymbol{x}\epsilon F_o\right] = \left[ \int \textrm{d}\boldsymbol{v}\,\textrm{d}\boldsymbol{x}\epsilon \left.\frac{\partial F_o}{\partial t}\right|_{\textrm{RF}} \right] = \frac{1}{2}Re \left[ \int \textrm{d}\boldsymbol{v}\,\textrm{d}\boldsymbol{x}\epsilon Q(F_o) \right] \nonumber\\ &= \frac{1}{2}Re \left[ \int \textrm{d}\boldsymbol{v}\,\textrm{d}\boldsymbol{x} \frac{q}{m} [\boldsymbol{E} + \boldsymbol{v} \times \boldsymbol{B} ]^* . [ \nabla_{\boldsymbol{v}}\epsilon ] f_{\textrm{RF}} \right] = \frac{q}{2}Re \left[ \int \textrm{d}\boldsymbol{v}\,\textrm{d}\boldsymbol{x}\boldsymbol{E}. \boldsymbol{J}_{\textrm{RF}} \right], \end{align}

in which $Q$ is the RF quasilinear diffusion operator needed in the Fokker–Planck equation,

(2.2)\begin{equation} Q(F_o)=\nabla_{\boldsymbol{v}}. \left[ \frac{q}{m} [\boldsymbol{E} + \boldsymbol{v} \times \boldsymbol{B} ]^* \int_{-\infty}^t \textrm{d} t' \frac{q}{m} [\boldsymbol{E} + \boldsymbol{v} \times \boldsymbol{B} ] . \nabla_{\boldsymbol{v}} \right] F_o, \end{equation}

and $\boldsymbol {J}_{\textrm {RF}}=\int \textrm {d}\boldsymbol {v} \boldsymbol {v}\,f_{\textrm {RF}}$ is the RF perturbed current density needed in the wave equation. Here $f_{\textrm {RF}}$ is the perturbed distribution function,

(2.3)\begin{equation} f_{\textrm{RF}}={-} \int_{-\infty}^t \textrm{d} t' \frac{q}{m} [\boldsymbol{E} + \boldsymbol{v} \times \boldsymbol{B} ] . \nabla_{\boldsymbol{v}} F_o. \end{equation}

In the above, $\epsilon$ is the energy, $q$ the charge and $m$ the mass of the considered particle population; $\boldsymbol {v}$, $\boldsymbol {E}$ and $\boldsymbol {B}$ are the velocity, the RF electric field and the RF magnetic field, respectively. Although the above identity suggests it is straightforward to ensure a self-consistent description, there is not a single suite of RF codes currently in existence that actually solves this problem generally (i.e. simultaneously accounting for all waves supported by the plasma and for the various arbitrary distribution functions of interacting minorities, majorities and beams) while putting the wave and Fokker–Planck equations on the same footing (i.e. using fully self-consistent expressions in the wave's dielectric response and the quasilinear diffusion operator). In practice, the wave and Fokker–Planck equations are solved making somewhat differing approximations, allowing a much more speedy description, while retaining the dominant physics. From the wave equation point of view the main simplification is that the dielectric response assumes the electric field varies slowly during a Larmor gyration – yielding a truncation (commonly to second order) in terms of finite Larmor radius corrections. From the Fokker–Planck point of view all finite Larmor radius effects are typically accounted for, but the electric field is assumed to only have a single perpendicular wave component mode, i.e. the homogeneous plasma response is assumed to capture the main dynamics although the dielectric response is actually a differential or integro–differential operator. A set of codes coming very close to self-consistency is the AORSA/CQL3D pair of solvers (Harvey et al. Reference Harvey and McCoy1992; Jaeger et al. Reference Jaeger, Berry, Ahern, Barrett, Batchelor, Carter, D'azevedo, Moore, Harvey and Myra2006; Petrov & Harvey Reference Petrov and Harvey2016) accounting for all finite Larmor radius corrections but considering a locally homogeneous RF diffusion operator, and very recent work by Lee et al. (Reference Lee, Wright, Bertelli, Smithe, Valeo, Petrov, Jaeger, Berry, Harvey and Bonoli2017) where, oppositely, only up to second-order finite Larmor radius correction terms are retained, but the wave diffusion operator is a differential operator in full agreement with the expansion used in the TORIC code.

All wave and Fokker–Planck codes installed in EU-IM make assumptions, some codes making more drastic simplifications than others. A common assumption is that guiding centres are glued to magnetic surfaces allowing us to study the dynamics of the collisional interaction of charged populations sharing a magnetic surface separate from transport processes radially across surfaces.

Making approximations to cut down on computer time necessarily comes at a price. The fact that the wave and Fokker–Planck equation are not derived making the same assumptions excludes the description of being truly self-consistent. Rather than matching the actual electric fields at the exact point of resonance, the wave and Fokker–Planck descriptions match the local power densities on each of the considered magnetic surfaces. Since the distribution function is a complicated (and for the case of cyclotron harmonic heating a very sensitive) function of the input power, this requires iteratively solving the Fokker–Planck equation on an appropriate (non-uniform) velocity grid covering an energy range, allowing us to capture the formation of high-energy tails as well as the slowing down and thermal ion populations until the proper rescaling factor is found, while freezing the wave polarisation and perpendicular wavenumber. A crude estimate of the maximal energy to be included is provided by locating the first zero of the RF diffusion coefficient (which acts as a barrier beyond which the distribution function drops significantly), typically in the range of a few MeV. If treating the problem fully self-consistently were possible this iteration procedure would be unnecessary and the electric field at the resonance could simply be adopted directly, the various Fokker–Planck equations being solved as an initial value, i.e. time stepping problem.

In the present work the focus is on fast and hence simple methods to model tokamak plasma heating scenarios in which various ion populations are simultaneously heated and the Coulomb collisional interaction of multiple non-Maxwellian (minority, majority and beam) populations is crucial. It is complementary to highly sophisticated but CPU time consuming solvers capable of providing a more detailed description of many aspects. Now AORSA is already routinely coupled with the CQL3D 3-D Fokker–Planck solver (Harvey et al. Reference Harvey and McCoy1992; Jaeger et al. Reference Jaeger, Berry, Ahern, Barrett, Batchelor, Carter, D'azevedo, Moore, Harvey and Myra2006; Petrov & Harvey Reference Petrov and Harvey2016), the EVE code is run together with the Monte Carlo code SPOT (Schneider et al. Reference Schneider, Johnson, Dumont, Eriksson, Eriksson, Giacomelli, Girardo, Hellsten, Khilkevitch and Kiptily2016; Joly et al. Reference Joly, Garcia, Imbeaux, Dumont, Schneider, Johnson and Artaud2019), TORIC is routinely used in conjunction with SSFPQL (Brambilla Reference Brambilla1994; Bilato et al. Reference Bilato, Brambilla, Maj, Horton, Maggi and Stober2011) as well as CQL3D (see for example Lee et al. Reference Lee, Wright, Bertelli, Smithe, Valeo, Petrov, Jaeger, Berry, Harvey and Bonoli2017) and the CYRANO wave code (Lamalle Reference Lamalle1994) was coupled to the BATCH Fokker–Planck code (Van Eester Reference Van Eester2001) when describing RF heating of a D beam at the fundamental cyclotron frequency (Lerche et al. Reference Lerche, Van Eester, Krasilnikov, Ongena and Lamalle2009). The Monte Carlo test-particle approach is increasingly popular since it allows a very intuitive insight by actually tracing a large number of orbits. Aside from the earlier mentioned SPOT, examples are ASCOT (Hirvijoki et al. Reference Hirvijoki, Asunta, Koskela, Kurki-suonio, Miettunen, Sipilä, Snicker and Äkäslompolo2014), SELFO (Hedin et al. Reference Hedin, Hellsten, Eriksson and Johnson2002) and SCENIC (Jucker Reference Jucker2010).

3. A suitable one-dimensional form of the Fokker–Planck equation

When magnetic surfaces are decoupled from each other except for transport terms mocked up by a local loss term (see the discussion later), a suitably simplified steady-state form of the one-dimensional (1-D) Fokker–Planck equation for distributions only depending on the velocity is (Karney Reference Karney1986)

(3.1)\begin{equation} \frac{1}{v^2}\frac{\partial}{\partial v} \left[ v^2 \left(A \frac{\partial F_o}{\partial v} + B F_o \right) \right] + S -\frac{F_o}{\tau_P}=0 \end{equation}

where $A$ has an RF and a collisional contribution ($A=A_{\textrm {RF}}+D_{vv,\textrm {coll}}$) where

(3.2)\begin{equation} A_{\textrm{RF}}=K_{\textrm{RF}}\left|J_{N-1}+\frac{E_-}{E_+}J_{N+1}+J_N\frac{E_{/{/}}}{E_+}\right|^2 , \end{equation}

while $B=B_{\textrm {coll}}+B_{\textrm {loss}}$ has a collisional term ($B_{\textrm {coll}}=-F_{v,\textrm {coll}}$) and a loss term

(3.3)\begin{equation} B_{\textrm{loss}}=\frac{v}{2}\left[\frac{1}{\tau_E}-\frac{1}{\tau_P}\right]; \end{equation}

$\tau _E$ is the energy loss time and $\tau _P$ the particle loss time. The wave polarisations of the RF electric field ($E_+$ is the electric field component rotating in the sense of the ions, $E_-$ rotates in the sense of the electrons and $E_{//}$ is the parallel electric field) are evaluated solving the wave equation.

The reason for only solving the time-independent solution is twofold. (1) It allows a rigorous check on energy conservation. Especially when the interaction of many populations needs to be tracked, energy conservation is a highly desirable asset. (2) It provides the desired solution faster. Solving the time-dependent Fokker–Planck equation showed that the setting up of the final distribution function happens in three characteristic times when all species distributions are allowed to be modified: (i) the direct impact of the RF waves (yielding a much faster response in the case of second harmonic heating than in the case of fundamental heating); (ii) the formation of the distribution under the influence of Coulomb collisions crudely in equilibrium for the ions; (iii) the formation of the distribution under the influence of Coulomb collisions in equilibrium both with the other ions and electrons. To enable simulations of full discharges, the speed of computation is deemed crucial in EU-IM, in spite of the fact that it requires sacrificing certain physics aspect details. Since EU-IM allows for time stepping, time-dependent solving of the Fokker–Planck equation is to be implemented at some future stage, however.

Note that only the relative polarisations appear in the above. Per magnetic surface and per specie, the constant $K_{\textrm {RF}}$ is adjusted to ensure the RF power density corresponding to the distribution obtained solving the Fokker–Planck equation is the externally imposed RF power density. The above Fokker–Planck equation is a simplified version of the full 3-D equation that ought to be solved. It ignores pitch angle effects and does not properly model (neoclassical collision as well as RF induced) radial diffusion across magnetic surfaces. Because of these simplifications, solving the equation is fast, making it a suitable tool for implementation in the ETS. For more detailed analysis of specific aspects related to for example particle trapping, this tool is too crude and more complete models need to be considered. The coefficients of the collision operator are sums on the contributions from each of the background species $b$ influencing the test species $a$, e.g. $F_{v,\textrm {coll}}=\sum _b F_{v}^{a/b}$.

When the background distribution is Maxwellian, the collision operator coefficients are

(3.4)\begin{gather} D_{vv}^{a/b}=\frac{\varGamma^{a/b}}{2v} \left[ \frac{\epsilon (u)-u \epsilon '(u)}{u^2} \right], \end{gather}

(3.5)\begin{gather}F_v^{a/b}={-}\frac{\varGamma^{a/b}}{v^2}\frac{m_a}{m_b} \left[ \epsilon (u)-u \epsilon '(u) \right], \end{gather}

in which $u=v/2^{1/2}v_t$ ($v_t$ is the thermal velocity) and $\epsilon$ is the error function. The constant $\varGamma ^{a/b}$ is

(3.6)\begin{equation} \varGamma^{a/b}=\frac{N_b q_a^2q_b^2ln\varLambda^{a/b}}{4{\rm \pi} \epsilon_o^2m_a^2} \end{equation}

in which $\ln \varLambda$ is the Coulomb logarithm. For a general distribution, these expressions generalise to

(3.7)\begin{gather} D_{vv}^{a/b}= \frac{4 {\rm \pi} \varGamma^{a/b}}{3N_b} \left[ \frac{1}{v^3} \int_0^v \textrm{d} p p^4 F_{o,b} (p) + \int_v^\infty \textrm{d} p p F_{o,b}(p)\right], \end{gather}

(3.8)\begin{gather}F_v^{a/b}={-}\frac{4 {\rm \pi} \varGamma^{a/b}}{N_b} \frac{m_a}{m_b v^2} \int_0^v \textrm{d} p p^2 F_{o,b}(p). \end{gather}

The total $D_{vv}$ and $F_v$ are the sum over the various species $b$ for a test species $a$, the corresponding masses are $m_{\cdots }$ and $N_{\cdots }$ the corresponding densities.

The adopted source is the broadened delta function $\delta (v_o-v)$ where $v_o$ is the birth velocity. The smooth version $\exp [-([v-v_o]/\varDelta _S)^2]/[v^2 {\rm \pi} ^{1/2} \varDelta _S]$ rather than the delta function itself is used to avoid breaks in the distribution's derivatives which may give rise to numerical issues when introducing the obtained distribution function into other EU-IM modules.

As the RF term nor the (non-simplified) collision term create or annihilate particles one readily finds that $\langle S\rangle =N_{a}/\tau _P$ where $\langle M(v)\rangle =4{\rm \pi} \int _0^\infty \textrm {d} v v^2 M(v)$. Similarly $\langle L(F_o) mv^2/2\rangle = -E_{\textrm {tot}}/\tau _E$ in which $E_{\textrm {tot}}$ is the total energy and $L(F_o)=L_E+L_P$ with

(3.9)\begin{gather} L_E(F_o)=\frac{1}{v^2} \frac{\partial}{\partial v} \left[\frac{v^3}{2} \left(\frac{1}{\tau_E}-\frac{1}{\tau_P}\right)F_o\right], \end{gather}

(3.10)\begin{gather}L_P(F_o)={-}\frac{F_o}{\tau_P}. \end{gather}

The incoming wave plus neutral beam power density is balanced by the collisional and non-collisional losses.

Unlike the Stix approach which allows integrating the Fokker–Planck equation for non-beam populations rapidly as a first-order differential equation using a simple semianalytical expression (see later), accounting for beam sources necessitates solving the actual second-order differential equation. Imposing vanishing outgoing flux both at $v=0$ and $v=\infty$ constitutes a logical set of boundary conditions which defines the unique solution compatible with a given source. Since boundary conditions need to be imposed at either end of the integration domain and in view of the fact that the distribution becomes very small at high energy, shooting methods are commonly ill-suited to solving the equation and a more global method is required. The finite element method is robust, but somewhat memory as well as time consuming. Some details of its practical implementation exploiting suitable third-order base functions allowing us to directly obtain $F_o$ and $\textrm {d} F_o/\textrm {d} v$ upon solving the linear system associated with the variational formulation of the differential equation are provided in appendix A.

It goes without saying that a 1-D isotropic distribution function does not capture all details of the actual 5D distributions. Effects such as particle trapping and RF induced neoclassical radial diffusion are not included. Moreover, it is known that RF heating mainly affects the perpendicular energy while leaving the parallel energy roughly untouched. The present tool is suitable for a fast evaluation of sets of interacting distributions and assessing the importance of this effect. The best known example of its use is the identification and qualitative description of synergistic effects: heating a tokamak plasma using beams as well as waves achieves a heating efficiency superior to the efficiency obtained as the sum of the two heating mechanisms apart. The basic reason for this synergy is that beam ‘preheating’ cranks up the temperature of the thermal bulk populations and makes the plasma less collisional. The simple picture of this phenomenon is that steady state of the energy level locally requires the balance of incoming energies and (Coulomb collisionally) lost energies. Reducing the friction allows for reaching higher energies for a given RF power level than is the case in the absence of beam ‘preheating’. In particular for harmonic heating, there is equally a second effect: as harmonic heating is only potent at higher energy, neutral beam injection provides a seed of already energetic particles on which the waves can act more efficiently. There is also an obvious beneficial synergistic effect of RF heating on beam populations: upon Coulomb collisionally slowing down, RF induced high-energy tails transfer a non-negligible fraction of their power to the electrons. As a consequence the electron temperature increases, which increases the critical velocity. This allows more of the beam power to flow to the ions, causing the bulk ion temperature to rise which in turn enhances the fusion rate.

In defence of the here discussed simple model adopted in ETS, it should be remembered that effects such as trapping only become important when a significant variation of the magnetic field strength is felt along particle orbits. When studying auxiliary heating close to the plasma core, trapping can be considered as a secondary effect. Also, the fact that isotropic distributions are considered is less dramatic than it seems at first sight: Stix (Reference Stix1992) showed that the tail energies predicted by a model assuming $F_o$ only depends on the velocity $v$ essentially differ by a factor of $2/3$ from the prediction with a model assuming the parallel velocity is irrelevant and assuming $F_o$ solely depends on the perpendicular velocity component $v_\perp$. Intuitively this makes sense, the incoming energy now needing to be redistributed over only two perpendicular velocities rather than over all three. Comparing two-dimensional (2-D) Fokker–Planck computations with the here adopted 1-D model, one can see that at low RF power densities Coulomb collisions manage to isotropise the distribution so that the approximation $F_o(v)$ is more suitable. At very high RF power densities the 2-D results yield tail energies more closely related to the results obtained assuming $F_o(v_\perp )$. At intermediate power densities the predicted tail energies are in-between the two. So, from a purely pragmatic point of view, one can make a prediction at the nominal power density and repeat it at $3/2$ that density to define the two ‘boundaries’ of the range of tail energies between which 2-D models would predict the actual tail energy. Hence assuming the distribution only depends on a single velocity component delivers information that is meaningful and does so at low computational cost. But evidently a simple model, by definition, should not be used to formulate definitive answers. It merely is a tool to highlight a possible road.

4. The Fokker–Planck equation revisited

In absence of particle sources and losses, the Fokker–Planck equation can be solved immediately: integrating once yields

(4.1)\begin{equation} v^2 \left(A \frac{\partial F_o}{\partial v} + B F_o \right) =ct. \end{equation}

The flux being zero at $v=0$ as well as at $v=\infty$, the integration constant has to be zero. Integrating further yields

(4.2)\begin{equation} F_o\propto \exp\left[-\int_0^v \frac{B}{A}\textrm{d} v\right], \end{equation}

the remaining integration constant being determined by the fact that the density is prescribed. Up to the minor upgrade that an energy loss has been allowed for and that the solution holds for any background and cyclotron harmonic, the obtained solution is the well known solution of Stix (Reference Stix1992); this upgrade was extensively discussed in Van Eester & Lerche (Reference Van Eester and Lerche2011). Rather than solving the Fokker–Planck equation for a single ‘heated’ specie, a set of coupled Fokker–Planck equations modelling the fate of the power transferred from the RF waves to various ions is typically solved by the Fokker–Planck model most standardly used in the ETS. The Fokker–Planck module computing the Coulomb collisional redistribution of RF power for the case when the test population has no source is labelled ‘StixReDist’.

When there is a particle source, a particle loss term has to be included for a steady-state solution to exist. The Stix procedure can now no longer be used but the reasoning in terms of the flux still holds. Integrating the full equation now yields

(4.3)\begin{equation} v^2 \left(A \frac{\partial F_o}{\partial v} + B F_o \right) - \frac{1}{\tau_P}\int_0^v \textrm{d} v v^2 F_0 + \int_0^v \textrm{d} v v^2 S =ct. \end{equation}

As for the sourceless case, the integration constant needs to be zero since the net flux from source and loss terms needs to be zero at infinity (the last two terms on the left-hand side need to balance one another since the first two approach zero, $F_0$ and its derivative vanishing at $v=\infty$). Defining a new dependent variable $Y=\int _0^v \textrm {d} v v^2 F_0$ ($Y$ is the density integral up to the factor $4{\rm \pi}$ arising from the integration over the pitch and gyroangle) and expressing the derivatives of this function in terms of the distribution $F_0$ and its derivative yields the equation

(4.4)\begin{equation} A \frac{\textrm{d}^2Y}{\textrm{d} v^2}+\left[-\frac{2A}{v}+B\right] \frac{\textrm{d} Y}{\textrm{d} v}-\frac{Y}{\tau_P}+\int_0^v \textrm{d} v v^2S=0. \end{equation}

As boundary conditions, we impose that the flux needs to be zero at $v=0$ and $v=\infty$. Filling out the relations $F_0=v^{-2}\,\textrm {d} Y/\textrm {d} v$ and $\textrm {d} F_0/\textrm {d} v=[v\,\textrm {d} ^2Y/\textrm {d} v^2-2\,\textrm {d} Y/\textrm {d} v]/v^3$, this condition can readily be transformed into the intuitive $Y(v=0)=0$, i.e. density integral at $v=0$ is zero, and $Y(v=\infty )=\tau _P \int _0^\infty \textrm {d} v v^2 S=N_a/[4{\rm \pi} ]$, i.e. the density integral up to $v=\infty$ comprises the whole distribution.

Aside from the ‘StixReDist’ module needed to solve the Fokker–Planck equation in the absence of particle sources, both versions of the Fokker–Planck equation have been installed in the ‘FoPla’ Fokker–Planck module of the ETS relevant for beam or fusion $\alpha$ populations and can be used upon properly setting one of the input flags. A number of solvers is also available and can be activated; this is equally done via input flags.

5. Numerical issues when solving coupled Fokker–Planck equations

Both the wave and Fokker–Planck equations are equations that need particular attention to ensure the solution found is accurate. From the wave equation point of view, the difficulty lies in the fact that the plasma admits wave modes of very different wavelengths and that both propagative and evanescent regions occur. For typical parameters in machines such as JET, Bernstein waves are millimetric while the fast wave has a wavelength of a few 10s of centimetres. Both need to be traced, forcing grid densities to be high. Tackling the difficulties of the wave equation falls beyond the scope the present paper. The difficulty of the Fokker–Planck equation lies in the fact that wave heating creates high energy tails, the particle density of which is extremely modest (the distribution functions are often six to eight orders of magnitude smaller in the high energy compared with the thermal region) while they can carry a substantial part of the wave energy. Inaccurately treating the high-energy part of the distribution yields results that at best bear weak resemblance to the actual solution. Back-of-the-envelope computations show that distributions decay exponentially as a function of energy under the influence of Coulomb collisions, but auxiliary and in particular RF heating subtly but importantly modifies this behaviour. Fundamental cyclotron heating is dominantly active in the thermal region but allows minority tails to reach energies in the MeV range. Because the waves most strongly interact with particles with modest energies in regions where the particle density is high, solving the Fokker–Planck equation for this type of heating usually does not require specific attention on v-gridding unless the power densities are high or the concentration of the population low (minority heating). Harmonic heating, on the other hand, is a more ‘touchy’ heating scheme since it does not affect the thermal population (except indirectly by fast particles slowing down) but rather accelerates a very small number of already relatively fast particles to even higher energies. Because it affects a very small subpopulation, the response of the distribution function to the power density is highly nonlinear. A modest change of the power density can cause a non-negligible change of the energy density at high energies without dramatically changing the particle distribution itself. Finding suitable forms of the Fokker–Planck equation that allow fast yet accurate integration is therefore a point of attention. For significant RF power densities, the tail commonly stretches up to the first zero of the quasilinear diffusion function containing a combination of Bessel functions (see (3.2)). Beyond this velocity, the distribution drops by several orders of magnitude. Subsequent zeros of the quasilinear diffusion operator are extra barriers through which only an increasingly small number of particles can leak, forming a distribution looking like a staircase when plotting it in log scale. Because of the drop in density when crossing the zeros, the solution is usually well described by choosing $v_{\max }$ somewhat larger than the velocity corresponding to the first zero. Truncating the integration grid at too small $v_{\max }$ and forcing the distribution or the flux to be zero there by imposing the boundary condition that no particle can escape to higher $v$ deceivingly gives the impression of a converged solution, in particular for harmonic heating cases where only a small number of already energetic particles can be accelerated to even higher velocities by the waves.

Solvers guaranteeing a speedier integration of the equations are of great benefit, especially when they are used in a suite of codes such as in the ETS. Many scenarios involve the interaction of various species, all of which undergo Coulomb collisions and therefore interact with each other. Because the set of coupled Fokker–Planck equations needs to be solved on each magnetic surface, and because of the fact that the quasilinear diffusion operator contains the factor $K_{\textrm {RF}}$ that needs to be adjusted to ensure the power density corresponding to the $F_o$s found when solving the Fokker–Planck equation equals the power density found by the wave equation on each magnetic surface, the Fokker–Planck equation has to be solved many thousands of times during a single run. A too slow Fokker–Planck solver makes the tool useless in practice. In the absence of a beam population the Fokker–Planck equation can be solved swiftly and accurately by a brute force semianalytical integration, as was demonstrated by Stix. This advantage is lost when a particle source term is introduced. Because boundary conditions need to be imposed at either end of the velocity interval in which the distribution exists, a global method of integration is needed. Classical methods such as the finite element or finite difference method allow for the recasting of the equation into a system of coupled linear equations that needs to be solved. Inverting the corresponding system matrix is much more time consuming than the interval-by-interval integration required in Stix's Fokker–Planck resolution method. Moreover, storing the system matrix requires significant memory, especially when not exploiting the fact that the system matrix is sparse. Iterative solvers accounting for the sparseness of the system matrix potentially allow us to push down the required computation time. As a check, the conjugate gradient method was explored alongside the classical Gaussian elimination suitable for full matrices. In spite of requiring a large number ($\approx$1000) of iterations involving matrix–vector multiplications – in which the sparseness is accounted for by storing the matrix in compact form, only keeping the non-zero elements – sparse-matrix solvers outperform the full-matrix solver. It was found, though, that iterative methods are useful to get a first good guess of the solution but that their use is limited in practice since the magnitude of the distribution function differs by many orders of magnitude, making the global error estimates by which ‘convergence’ is assessed by sparse solvers a poor representative of actual local convergence. A method that allows avoiding the need to invert a big matrix is the influence matrix method. Because this is an interesting while relatively unknown technique, it is described in some detail in appendix B.

6. Coupled Fokker–Planck equations, collisions and losses

The dynamics of the interaction of populations of charged particles in the presence of RF heating sources has many aspects. It is customary to treat these aspects by different equations assuming they appear on different time scales. The fastest time scale is that of the wave oscillation at the frequency of the generator (MHz range; modelled by the wave equation prescribing the driver frequency $f$, i.e. assuming all perturbed quantities vary as $\exp (-\textrm {i}\omega t)$ with $\omega =2{\rm \pi} f$, and allowing us to average over an oscillation cycle when defining the quasilinear diffusion coefficients needed in the Fokker–Planck equation). The slowest time scales describes the slow evolution of macroscopic quantities such as the density and temperature (Hz range; modelled by the transport equations and typically time dependent). In-between these time scales lies the time scale of interest for the Fokker–Planck equation (kHz range for bounce orbit dynamics, allowing a smoothing/averaging over many bounce cycles when composing the input sources for the transport equations). Some of the phenomena described by the Fokker–Planck equation are, however, occurring on time scales not well separated from transport phenomena. To keep the computational time reasonable, the Fokker–Planck equation is typically solved on a set of magnetic surfaces, the implicit reasoning being that orbits do not deviate much from an average radial position. This approach is deemed suitable in the present paper. More sophisticated Fokker–Planck models add another dimension to the description and account for the cross-talk between magnetic surfaces, i.e. the Fokker–Planck equation is solved for the whole plasma volume.

The simplest and most routinely used Fokker–Planck solvers look at the problem of Coulomb collisional slowing down of a single minority population that is heated by waves but has reached a steady state. This is, for example, the approach adopted in the TRANSP code (Hawryluk et al. Reference Hawryluk1980). In that case energy losses are commonly not even considered and the net power balance is straightforward: the incoming wave power is collisionally passed on to the various bulk plasma constituents. Because the minority density is small, the self-collisions are ignored and the linearised version of the collision operator can be used. Moreover, because the minority is the only heated population, the other populations must be close to Maxwellian so that simplified expressions for the collision operator can be exploited. This is qualitatively depicted in figure 2(a). The resulting collisional redistribution profile is a source term for solving the energy transport equation for the bulk species. The case of beam heating requires a mild upgrade: when fresh neutral beam source particles are continually injected, a steady state can only be reached if the particle loss equals the particle input. In case the density is assumed to remain the same and if the losses are represented by a simple constant particle loss time, this loss time is not a free parameter: it is prescribed by the local density and the beam source strength. But, there still is no need for an energy loss term.

Figure 2.

Philosophy classically used for minority RF heating (a), generalisation to all species (b) and inclusion of loss terms accounting for losses from local magnetic surfaces through radial diffusion (c).

When multiple populations are being heated and/or are interacting with each other, a set of coupled Fokker–Planck equations rather than a single Fokker–Planck equation needs to be solved. When including majority species in the set of populations described, self-collisions can no longer be neglected and the nonlinear collision operator for arbitrary distributions needs to be used. Since Coulomb collisions do not create nor annihilate particles or energy, the sums of all powers collisionally transferred from one specie to another is zero: the nonlinear collision operator $C(F_{o,A},F_{o,B})$ describing collisional interaction between populations $A$ and $B$ satisfies the conservation relation $\langle C(F_{o,A},F_{o,B})+C(F_{o,B},F_{o,A})\rangle =0$ where $\langle \cdots \rangle$ denotes the energy integral over velocity space; $C(F_{o,A},F_{o,B})$ is the collision term appearing in the Fokker–Planck equation solving for $F_{o,A}$ assuming $F_{o,B}$ is known, while $C(F_{o,B},F_{o,A})$ is the collision term appearing in the Fokker–Planck equation solving for $F_{o,B}$. Adopting the generalised version of the logic sketched for a heated minority, this poses a fundamental problem (see figure 2b): whereas injected beam or wave power was flowing away to ‘bulk’ populations in the case of a single heated minority, hereby allowing a steady state – guaranteeing energy conservation – to be reached, there is now a net local energy influx but there is no net local energy loss. As a consequence it is impossible to reach a steady-state solution. So when describing the Coulomb collisional interaction between all species, an energy loss term needs to be introduced for a steady state to be reached, even when no beams are involved (see the bottom subplot of figure 2). The physical origin of this loss term is readily evident when looking at the 3-D Fokker–Planck equation: radial transport evacuates the local heat to neighbouring magnetic surfaces and the only global physical loss is power leaking out of the confined plasma across the last closed flux surface. That losses are an inherent part of magnetic confinement discharges is well known experimentally: to reach a steady state, the plasma needs to be replenished both by energy and particles.

In the light of the realisation that the Fokker–Planck equation needs to include energy loss channels, we can now return to the isolated minority being heated and for which the introduction of a loss term was not deemed necessary. When a steady state is reached, the power collisionally transferred from the minority to the bulk ions and electrons has to be found back – up to a minus sign – in the power balance corresponding to the Fokker–Planck equations of these bulk species. For example, the bulk ions – which are assumed not to be directly heated by the waves when considering a single heated minority – now receive power from the minority. This power needs to be evacuated for a steady-state solution to exist. Note, however, that there is a shift in the reasoning: when solving the isolated Fokker–Planck equation for the minority, the steady state was achieved when the net effect of the RF waves was balanced by the net effect of the collisions. A higher density would for example give rise to more efficient evacuation of power and to lower tail energies. Also, at a prescribed density, a higher RF power yields higher energy tails. Solving the coupled Fokker–Planck equations for the minority and the bulk ions, the net effect of the collisions is cancelled out and the input RF power gain is ultimately compensated by a transport loss. So the local minority tail energy is partly determined by the amount of power locally deposited by the waves, as before, but also by the efficiency with which bulk energy is transported across magnetic surfaces. The proper choice of a suitable loss time $\tau _E$ for the Fokker–Planck equation is therefore crucial. Nor is the role of the collisions lost: the relative efficiency of the collisional transfer to the bulk ions and electrons determines which fraction of the power will be evacuated to each of the bulk populations. When multiple ion species are involved, the ultimate steady state depends on the relative strengths of the various collisional transfers as well as the loss channels and the strength of the power source. Adding losses to the Fokker–Planck equation – even when no beams are involved – is a deviation from what is traditionally done. It is a necessity when accounting for the Coulomb interactions of all species.

One possible way of defining a tau would be by connecting it with the local energy radial diffusion coefficient and local scale length, i.e. identifying the transport term $\partial /\partial \rho [\chi _{\rho \rho } \partial /\partial \rho \cdots ]$ with the energy loss term $\cdots /\tau _E$ in the Fokker–Planck power balance equation. Transport is, however, not described by a single coefficient. At the very least some connection should be made with the conduction as well as the convection coefficient, both of which are typically finite. Moreover, since the transport and Fokker–Planck equations do not model the same populations (by construction, the transport equation strictly only models the Maxwellian subpopulations) a $\tau _E$ that is relevant for the thermal population is not necessarily suitable for the high-energy subpopulation. Rather than imposing a particular energy confinement time externally and as the focus of the present paper is on describing the impact of ICRH power on a discharge, the way the characteristic time $\tau _E$ is defined here is by ensuring that the loss term in the Fokker–Planck equation is consistent with the reigning transport. This is done by requiring that the initial (experimental) temperature is reproduced when solving the Fokker–Planck equation omitting the added RF power. Only in case no experimental data is available to permit deriving $\tau _E$, it is imposed as an external parameter based on a suitable scaling law. It is implicitly assumed that $\tau _E$ does not vary too much when a finite amount of (supplementary) wave power is launched. Moreover, and not necessarily realistic, a single loss time has been adopted for the whole population whereas radial transport is expected to be different for thermal and high-energy subpopulations of a given population. The power locally lost is simply added to the self-collision's channel when composing the energy source terms fed into the transport solver. Since the losses represent radial transport and hence represent a loss of energy evacuated away from a magnetic surface to its neighbours, this step implicitly assumes the transport on neighbouring surfaces is similar so that the energy lost from a given surface is closely matched by the energy influx into it from the neighbouring ones.

7. Computation of the fusion rate

A key variable of interest for a heating scenario is the fusion power it allows to be produced. In this section it is described how the fusion rate for interacting arbitrary 1-D distributions $F_o(v)$ can be computed. The total fusion reactivity per unit volume is computed by evaluating

(7.1)\begin{equation} \langle\sigma(V) V\rangle=\int \textrm{d}\boldsymbol{v}_1 \int \textrm{d}\boldsymbol{v}_2 F_{0,1} (\boldsymbol{v}_1)F_{0,2}(\boldsymbol{v}_2) \sigma(V) V \end{equation}

in which $V=|\boldsymbol {v}_2-\boldsymbol {v}_1|$ is the relative velocity. Approximate expressions exist for the reactivity. Adopting the parameterisation proposed by Bosch and Hale (Bosch & Hale Reference Bosch and Hale1992) allows us to write the reaction rate in the form

(7.2)\begin{equation} \sigma=\frac{S(E)}{E \exp(B_G/E^{1/2})} \end{equation}

in which

(7.3)\begin{equation} S(E)=\frac{A_1+E(A_2+E(A_3+E(A_4+EA_5)))}{1+E(B_1+E(B_2+E(B_3+EB_4)))}. \end{equation}

Here $B_G$ is the relevant Gamov constant, $E$ is expressed in $keV$ and the resulting reactivity $\sigma$ is given in $mb$ (millibarn; $1\,\textit {barn}=10^{-28}\ \textrm {m}^2$). For relevant fusion reactions, the values of the various constants can be found in Bosch & Hale (Reference Bosch and Hale1992). The energy is defined via the reduced mass and the relative velocity $E=m_r V^2/2$ where $m_r=m_1m_2/[m_1+m_2]$. If one population is thermal while the other is fast, then (assuming species 2 is the slower population) $V\approx v_1$ and $E=[m_2/(m_1+m_2)] E_1$. So while the maximal reactivity occurs at $E=65\ \textrm {keV}$, this maximum falls at $E_D=108\ \textrm {keV}$ for the deuterons when the tritons are thermal, and at $E_T=162\ \textrm {keV}$ for the tritons when the deuterons are thermal.

The cross-section only depends on the magnitude of the relative velocity and hence contains no information on the angular dependence. When assuming that the distribution functions themselves also only depend on the magnitude of the velocities of the two types of species, the computation of the integrated reaction rate can be simplified. Writing $\langle \sigma V\rangle$ explicitly in terms if spherical coordinates yields

(7.4)\begin{equation} \int_0^\infty \textrm{d} v_1 v_1^2 \int_0^{\rm \pi} \textrm{d}\theta_1 \sin\theta_1 \int_0^{2{\rm \pi}} \textrm{d}\phi_1 f(\boldsymbol{v}_1) \int_0^\infty \textrm{d} v_2 v_2^2 \int_0^{\rm \pi} \textrm{d}\theta_2 \sin\theta_2 \int_0^{2{\rm \pi}} \textrm{d}\phi_2 f(\boldsymbol{v}_2) \sigma(V) V. \end{equation}

Choosing the first axis of the second spherical coordinate system to be aligned with the first vector $\boldsymbol {v}_1$ (see figure 3) allows us to write $\boldsymbol {v}_2$ as

(7.5)\begin{equation} \boldsymbol{v}_2=v_2[\cos \theta_2\boldsymbol{e}_{/{/}}+\sin\theta_2 \cos\phi_2\boldsymbol{e}_{{\perp},1}+\sin\theta_2 \sin\phi_2\boldsymbol{e}_{{\perp},2}] \end{equation}

where $\boldsymbol {e}_{//}=\boldsymbol {e}_{\boldsymbol {v}_1}=\boldsymbol {v}_1/|\boldsymbol {v_1}|=\boldsymbol {v}_1/v_1$. We can rewrite

(7.6)\begin{align} V &= |\boldsymbol{v}_2-\boldsymbol{v}_1|=[(v_2\cos\theta_2-v_1)^2+(v_2\sin\theta_2\cos\phi_2)^2+ (v_2\sin\theta_2\sin\phi_2)^2]^{1/2} \nonumber\\ &= [v_1^2+v_2^2-2v_1v_2\cos\theta_2]^{1/2} \end{align}

from which it follows that $2V\,\textrm {d} V=2v_1v_2\sin \theta _2\,\textrm {d}\theta _2$. When the distributions only depend on the magnitudes of the velocities then the only dependence on $\theta _1$ and $\phi _1$ is in the Jacobian. The integral then can be reduced to

(7.7)\begin{equation} \langle \sigma(V)V\rangle=\int_0^\infty \textrm{d} v_1 f_1(v_1) \int_0^\infty \textrm{d} v_2 f_2(v_2) I(v_1,v_2) \end{equation}

in which

(7.8)\begin{equation} I(v_1,v_2)=8{\rm \pi}^2 v_1 v_2 \int_{|v_1-v_2|}^{|v_1+v_2|}\textrm{d} V\sigma(V) V^2. \end{equation}

The latter integral can be tabulated once and for all on a grid on $v$-values, and a cubic interpolation allows us to approximate the values of the integral away from the grid points. Then using the trapezium rule to evaluate the integrals, the integral can be written as the vector–matrix–vector product

(7.9)\begin{equation} \langle\sigma(V)V\rangle=\boldsymbol{F}_{o,1}.\bar{\bar{A}}.\boldsymbol{F}_{o,2} \end{equation}

in which $F_{o,1,i}=F_{o,1}(v_{i})$, $F_{o,2,j}=F_{o,2}(v_{j})$ and

(7.10)\begin{equation} A_{i,j}=8 {\rm \pi}^2 \left[v_{i} {\rm \Delta} v_{i} \right] \left[v_{j} {\rm \Delta} v_{j} \right] \left[ I_o(|v_i+v_j|)- I_o(|v_i-v_j|) \right], \end{equation}

where

(7.11)\begin{equation} I_o(V)= \int_{0}^{V}\textrm{d} V\sigma(V) V^2. \end{equation}

The total fusion yield can be computed fast by performing the above given integrals. As inputs, the distribution functions of the various interacting populations obtained by solving the coupled Fokker–Planck equations are used. Aside from the D–T, D–$^3\textrm {He}$ and both the D–D cross-sections available in Bosch and Hale (Bosch & Hale Reference Bosch and Hale1992), the T–T cross-section has been included (Coster et al. Reference Coster, Basiuk, Pereverzev, Kalupin, Zagorski, Stankiewicz, Huynh and Imbeaux2010).

Figure 3.

Unit vectors for spherical coordinates for which $\boldsymbol {e}_x$ is aligned with $\boldsymbol {v}_1$.

8. Example: ‘baseline’ $({}^3\textrm {He})-\textrm {D}-\textrm {T}-\textrm {D}_{\textrm {NBI}}-\textrm {T}_{\textrm {NBI}}$ scenario

So far, this paper has been a rather technical report on fast tools available in EU-IM allowing the study of heating scenarios involving minorities as well as majorities and beams, all simultaneously – directly or indirectly – being heated by RF waves, energetic (neutral) particle injection and/or collisional interaction. In this section a single example is provided. The accent is on illustrating the type of output generated rather than making a deep physics study, the reason being that the accent of ETS and EU-IM publications is typically on the latter while assuming knowledge about the modules behind has been advertised elsewhere. The interested reader desiring to have more insight in the EU-IM framework on the one hand or on actual physics application on the other is referred to papers listed at the outset.

Now, JET is preparing for a second D–T experimental campaign (DTE2), see for example Garzotti et al. (Reference Garzotti, Challis, Dumont, Frigione, Graves, Lerche, Mailloux, Mantsinen, Rimini and Casson2019); Joffrin et al. (Reference Joffrin, Abduallev, Abhangi, Abreu, Afanasev, Afzal, Aggarwal, Ahlgren, Aho-mantila and Aiba2019). The first (DTE1) campaign, executed in the1990s, showed that fusion could be realised abundantly in JET (see for example Jacquinot et al. Reference Jacquinot, Bhatnagar, Cordey, Horton, Start, Barnsley, Breger, Christiansen, Clement and Davies1999; Keilhacker et al. Reference Keilhacker, Gibson, Gormezano, Lomas, Thomas, Watkins, Andrew, Balet, Borba and Challis1999; Start et al. Reference Start, Jacquinot, Bergeaud, Bhatnagar, Conroy, Cottrell, Clement, Ericsson, Eriksson and Fasoli1999). One major difference between the JET machine back then compared with now, is the presence of the ITER-like wall (Philipps et al. Reference Philipps, Mertens, Matthews and Maier2010; Joffrin et al. Reference Joffrin, Baruzzo, Beurskens, Bourdelle, Brezinsek, Bucalossi, Buratti, Calabro, Challis and Clever2014). A key motivation for moving away from a carbon wall was for better control of the in-vessel tritium (T) inventory, crucial in ITER to meet regulatory requirements, plasma-facing components having a very high hydrogen fuel retention. One of the first key results from JET–ILW was the confirmation of low retention of hydrogen fuel compared with JET-C, see for example Brezinsek et al. (Reference Brezinsek2013). The beryllium wall and tungsten divertor have significantly modified the boundary conditions of experiments. One issue to be tackled is the potential of accumulation of high Z material in the core. Core concentrations of tungsten larger than a few $10^{-4}$ cause the core temperature to collapse when not enough auxiliary heating power reaches the centre to compensate the high radiation. Significant radiation is also observed close to the low field side edge of the plasma, attributed to tungsten accumulation resulting from interaction with the W divertor or W-covered surfaces (e.g. the NBI dump plates to protect the machine against shine-through) in the machine. The temperature collapse in the core can result in a global thermal quench provoking a disruption. Ion cyclotron resonance heating has shown to ensure high Z particles are expulsed from the core (Goniche et al. Reference Goniche, Lerche, Jacquet, Van Eester, Bobkov, Brezinsek, Colas, Czarnecka, Drewelow and Dumont2014; Lerche et al. Reference Lerche, Goniche, Jacquet, Van Eester, Bobkov, Colas, Monakhov, Rimini, Czarnecka and Crombé2015; Van Eester et al. Reference Van Eester, Lerche, Kazakov, Jacquet, Bobkov, Crombé, Czarnecka, Dumont, Eriksson and Giacomelli2016b) although a demonstration in D–T relevant circumstances still needs to be done. Apart from providing input for ITER, one of the goals of the second D–T campaign is to demonstrate that the results obtained in the carbon wall machine can still be achieved in presence of a metal wall. Progress in understanding has, however, led to the hope that the then set standards can be bypassed: a key target of DTE2 is the demonstration of sustained DT fusion power of $10$ to 15 MW for 5 s (Litaudon et al. Reference Litaudon, Abduallev, Abhangi, Abreu, Afzal, Aggarwal, Ahlgren, Ahn, Aho-mantila and Aiba2017). A ‘key performance indicator’ (known as KPI) used to monitor the progress of plasma preparation in deuterium was the sustained D–D neutron rate, with a target of $5\times 10^{16}\ \textrm {n}\ \textrm {s}^{-1}$ (equivalent D–D neutron performance of a plasma capable of 15 MW DT power, according to $0$–D calculations with assumption on $H_{98(y,2)}$).

Neutral beam injection is JET's dominant heating scheme: in the recent campaigns it was shown that 32 MW of beam power can be launched into the JET plasma making use of beam energies of approximately 120 keV. Switching off the beam to terminate the shot has shown to hold a risk: the particle and the energy confinement time not being the same, the core temperature typically starts dropping faster than the core density. To counteract – similar to what was discussed above for the main heating phase – the collapse of the temperature profile, a heating well localised in the core is needed after the main heating has been switched off. While beam deposition profiles are commonly very wide (the beam ionisation starts taking place as soon as the neutral beam penetrates a plasma of sufficient density), ICRH can deposit power in a much more localised way. For that reason, ICRH has become of the tool of reference to chase high Z impurities from the core and to ‘land’ JET discharges (Goniche et al. Reference Goniche, Lerche, Jacquet, Van Eester, Bobkov, Brezinsek, Colas, Czarnecka, Drewelow and Dumont2014; Lerche et al. Reference Lerche, Lennholm, Carvalho, Dumortier, Durodie, Van Eester, Graves, Jacquet and Murari2017; Litaudon et al. Reference Litaudon, Abduallev, Abhangi, Abreu, Afzal, Aggarwal, Ahlgren, Ahn, Aho-mantila and Aiba2017). Understanding the interplay between NBI and ICRH is crucial for optimising high performance JET shots.

The preparations of a D–T campaign have highlighted two main candidates for reaching the goal set: the ‘baseline’ scheme (see for example Kim et al. Reference Kim, Sips, Romanelli, Challis, Rimini, Garzotti, Lerche, Buchanan, Yuan and Kaye2018) operates at high magnetic field $B_o$ and current $I_p$ while maintaining a high plasma density $N_e$. As high $B_o$ and $I_p$ require a high auxiliary heating power, this scheme tends to operate at a modest normalised $\beta$. It holds the promise for a long stationary state but suffers from a somewhat high disruption rate, and tends to have MHD activity (fishbones and NTMs, the latter more often than not triggered by violent sawtooth crashes) that degrade its performance. The ‘hybrid’ scheme (see for example Joffrin et al. Reference Joffrin, Sips, Artaud, Becoulet, Bertalot, Budny, Buratti, Belo, Challis and Crisanti2005) relies on current profile tailoring, and typically operates at higher normalised $\beta$. Because of its higher safety factor this scheme does not have sawtooth crashes and suffers less from disruptions.

The RF heating scheme currently mostly used during the ongoing preparations has been the RF ‘work horse’ since the machine's inauguration in 1983: core fundamental cyclotron heating of an H minority in a D plasma. The neutral beam injectors are typically injecting deuterium, and hence the latter as well as the main ion species absorb power in the plasma centre close to their second cyclotron harmonic layer. In the actual D–T campaign, half of the majority D will be substituted for T. It was shown in DTE1 that adding a small amount of ${}^3\textrm {He}$ to the plasma mix and heating it with ICRH boosts the neutron rate significantly. A ${}^3\textrm {He}$ minority can be equally be heated at its fundamental cyclotron layer in JET, this forming with T a pair such as (H)–D as the second harmonic T cyclotron layer coincides with the fundamental cyclotron layer of ${}^3\textrm {He}$. Both the baseline and the hybrid scenarios operate in H-mode.

In what follows the various simulations will be loosely based on the JET experimental baseline shot 92436 but rather than the plasma composition of that shot (H minority in a D majority plasma in which a D beam is fired) we will consider the case of a ${}^3\textrm {He}$ minority in a balanced D–T majority in which balanced D and T beams are injected. The minority concentration is taken as in the reference experiment (${\approx }1.6\,\%$ in the core) while the beam concentration for each of the beam boxes is of the order of 2–3 %. The adopted density and temperature profiles for the bulk species are depicted in figure 4. The central electron density is $N_{e,0}\approx 7\text {--}8\times 10^{19}\ \textrm {m}^{-3}$ while the density at the pedestal is $N_{e,0}\approx 4.5\times 10^{19}\ \textrm {m}^{-3}$, the central electron temperature is $T_{e,0}\approx 8\ \textrm {keV}$. Since there is no experimental ion temperature data in the actual core (the closest $T_i$ data point is approximately 20 cm from the centre), the ion temperature profiles adopted for the present example (with core ion temperatures reaching $T_{i,0}\approx 11\ \textrm {keV}$) have been computed inside ETS when solving for the ion heat equation using power density input profiles provided by the tools discussed in the present paper. When modelling the original JET baseline scenario, the computed ion temperature matches the experimentally found values in regions where experimental data is available (Huynh et al. Reference Huynh, Lerche, Van Eester, Garcia, Johnson, Ferreira, Kirov, Yadykin and Strand2020). Hence – although the plasma composition is vastly different – there is some confidence that the obtained profiles are reasonable, one reason being that the wave polarisation dynamics associated with the $({}^3\textrm {He})$–T pair (sharing the same cyclotron damping region, the former for minority fundamental and the latter for majority second harmonic heating) is similar to that of (H)–D in the reference shot. More details on this computation as well as a deeper study of the DT scenarios choosing either H or ${}^3\textrm {He}$ as minorities can be found in Huynh's paper. The toroidal magnetic field strength is $B_o=3.28T$ (adapted to ${}^3\textrm {He}$ core heating and hence differing from the actual $B_o$ in shot 92436) while the plasma current is $I_p=3MA$. The effective beam power (after subtracting for example charge exchange and shine-through losses) of the reference shot was approximately $P_{\textrm {NBI}}=26\ \textrm {MW}$ while the RF power was $P_{\textrm {RF}}\approx 5\ \textrm {MW}$. In the present example, the D beam is substituted by a D and a T beam. The energy of each D injector is taken to be $125keV$ and that of each T injector is 118 keV, the latter lower than the maximum possible voltage in view of engineering restrictions for safe operation. Each beam has a full-energy, a half-energy and a third-energy component. For launching the RF waves, dipole phasing was exploited and the driver frequency was chosen to be 33 MHz (aiming to heat the ${}^3\textrm {He}$ minority and the T majority) to ensure core ion heating for the adopted magnetic field.

Figure 4.

Density and temperature profiles of the bulk species.

To compute the beam source profiles, ASCOT (Hirvijoki et al. Reference Hirvijoki, Asunta, Koskela, Kurki-suonio, Miettunen, Sipilä, Snicker and Äkäslompolo2014) was used, while CYRANO (Lamalle Reference Lamalle1994) was exploited for describing the wave propagation and collisionless damping on various populations. Next, StixReDist (Van Eester & Lerche Reference Van Eester and Lerche2011) and FoPla (described earlier) were used for computing the respective distributions of the various interacting populations, the former for test particle populations that do not have a beam source, the latter for the ones that do.

Figure 5 depicts the ICRH and NBI power deposition profiles. Panels (a, b) show the RF power flowing to various of the species for a total incoming RF power of 5 MW. The power density peaks slightly away from the core. The ${}^3\textrm {He}$ minority is the most efficient absorber ($47\,\%$). Sharing its second harmonic layer with the fundamental cyclotron layer of the minority is the majority T, it absorbs $24\,\%$ of the incoming RF power. Also the highest energy T beam absorbs some of the power ($5\,\%$). The bulk D has its fundamental cyclotron layer on the high field side yielding the D absorption maximum at $\rho \approx 0.65\ \textrm {m}$. These species absorb $10\,\%$ of the incoming waves power. The other ions absorb very little wave power. The electrons, finally, have a wide deposition profile but have their highest power density near the plasma core; $11\,\%$ of the power flows from the waves to the electrons.

Figure 5.

The ICRH and NBI local and integrated power deposition profiles for a total integrated beam power of 26 MW and RF power of 5 MW.

Figure 5(c, d) shows the deposition profiles of the various beam populations. There are six populations in total, three of which correspond to the D beam and three to the T beam energy components. The beams are composed of a dominant beam fraction launched with an energy of $125\ \textrm {keV}$ (D) or $118\ \textrm {keV}$ (T), a fraction that is launched at half that energy and a small fraction, finally, that is launched at 1/3 of the energy. The highest energy beams (both D and T) penetrate deepest into the plasma and have a deposition profile that peaks close to the centre. In contrast, the lowest energy beams ionise close to the plasma edge and hence have a density that peaks there. The deposition profile of the half-energy beams is roughly flat throughout the cross-section. Most of the beam energy is carried by the high-energy beams; the full energy component also has the highest density. Especially the third-energy beam is a mere correction to the deposition of the full-energy beam component.

Although the RF power level is much lower than that of the beams, the RF power density is of the same order of magnitude as that of the beams in the core. This is a consequence of the fact that RF heating is well localised and that core wave heating was envisaged. Because the beam deposition profile is determined by ionisation of the incoming neutral beam, and since the plasma density is high (and fairly flat), the beams start depositing power as soon as they penetrate the plasma, yielding a very wide deposition profile. Since power lost closer to the edge is smeared out over a large volume, the corresponding power density is low in spite of the beam power being a factor of six higher than that of the RF waves. As a result of the relative magnitude of the contributed power densities it is fair to say that the width of the temperature profile is mainly set up by the beams while core localised RF power allows to have a significant impact in the highest density, highest temperature region of the plasma.

Figure 6 shows how the beam power is collisionally redistributed in ‘absence’ of RF power. Since the power densities remaining in the beam populations is very modest, these contributions have not been retained in the plot. Because the critical energy (at which the population that is slowing down passes on half of its energy to the electrons and half to ions; $E_{\textrm {crit},j}=14.8 kT_e(eV) A_{j} (\sum _i X_iZ_i^2/A_i)^{2/3}$ where $X_i$ is the concentration and the index $i$ loops on the various ion populations while the index $j$ refers to the test species) scales with the electron temperature, electrons tend to receive more of the beam power closer to the plasma edge. Higher energy beams tend to pass on a bigger fraction of their power to the electrons, their energy exceeding the critical energy more. The combination of the two effects results in high-energy beams heating both ions and electrons in the core while the lower energy beams favour ion heating. When including all beam energies, ion heating is dominant in the core. Near the edge, the high-energy beams mainly crank up the electron temperature. For the lower energy beams, ions and electrons play a similar role towards the outer layers of the plasma. Here, electron heating is the dominant effect. Note that the deuterium is systematically favoured over the tritons for heating via Coulomb collisional power transfer. The reason for this behaviour is related to the relative mass: because T is more massive, its thermal velocity for the same temperature is smaller. Hence, fast populations slowing down will first collisionally interact with D ions and only afterwards with T ions. This is true independent of which is the fast subpopulation.

Figure 6.

Coulomb collisional power distribution of the various beam species considered in absence of ICRH power.

The fate of energy carried to the core of the plasma by waves is critically dependent on the RF power density. Lower power densities prevent high-energy tails being formed and hence favour RF power flowing to the ions. Higher levels of wave power launched allow the transfer of more power to the electrons. This is clearly illustrated in figure 7, depicting the total RF+NBI Coulomb collisional power redistribution for four different RF power levels: $0.1$, $3$, $5$ and $7$ MW at constant $P_{\textrm {NBI}}=26$ MW beam power. Again, the modest power flowing to the beam populations has been omitted in the plot. At low power each of the core bulk ions receives double the power flowing to the electrons, in agreement with the pure NBI heating case shown in figure 6. When the RF power level is $3MW$, the core power density is significantly higher than that in absence of ICRH. Also, the power density collisionally passed on to the bulk ions shows a distinct peak. When 5 MW of RF power is coupled to the plasma, the electron power deposition also peaks. Compared with what was the case in absence of RF power, the electrons now receive almost three times more power in the core. This effect becomes gradually more pronounced if even higher RF power levels can be achieved. The extra power is now preferentially flowing to the electrons. The combined heating of electrons and ions in the core allows the boosting of the core temperatures, the $\textrm {RF}+\textrm {NBI}$ power density exceeding the NBI-only power density by a factor two or more in spite of the most amount of RF power compared with that of the beams. As noticed before, the deuterons are receiving more power than the tritons. The increased electron heating has an indirect but very important effect: increasing the electron temperature increases the critical energy and allows more of the beam power to flow to the ions, allowing a further boost of the ion temperature.

Figure 7.

Collisional power redistribution for four different ICRH powers launched: $0.1$, $3$, $5$ and $7$ MW. The beam power is the same for all cases: 26 MW.

Figure 8 zooms in on the redistribution of the two ions receiving most of the incoming wave power: the ${}^3\textrm {He}$ minority and the T majority. There is a distinct difference between the two behaviours. While the minority passes on similar amounts of power to the electrons and each of the bulk ion species (so that the bulk ion heating dominates the electron heating), the tritons heat the ions much more than they do the electrons. Since T is heated at its second harmonic, this comes as a surprise at first: second harmonic heating creates highly energetic tails and leaves the thermal region untouched. The reason is twofold. (i) Since various populations absorb a part of the power, the power densities are relatively modest and the tails are not particularly energetic, allowing to contribute to ion rather than electron heating. (ii) As mentioned before, the beam is the dominant heater of the plasma. Slowing down beam particles feed energy to the thermal population and increase its temperature. Because the results shown are not the result of solving a standalone Fokker–Planck equation for the beam but of solving the coupled set of equations while accounting fully for collisional interaction of the various species (including self-collisions of the majority population), the ion thermal temperature is being pushed up to higher values. Note that opposite to ‘pure’ slowing down which was more of benefit to D rather than to T ions, it is now the tritons that receive most of the redistributed power, partly through the loss term (radial diffusion was assumed to vary little between neighbouring magnetic surfaces, allowing us to identify the radial diffusion loss to neighbouring surfaces with the radial diffusion gain of power from these neighbouring surfaces) and since the tritons absorb wave power themselves. This underlines that single specie reasoning and multi-specie reasoning are different.

Figure 8.

Coulomb collisional power distribution of the core RF heated non-beam species T majority and ${}^3\textrm {He}$ minority for a total integrated RF power of 5 MW.

What sets the ICRH modelling presented here in the EU-IM framework apart from usual ICRH modelling is best illustrated by figure 9: whereas heating scenarios commonly involve various populations exchanging energy via Coulomb collisions, this aspect is commonly unaccounted for. Presenting a method allowing the study of the simultaneous heating of minorities, majorities and beams while accounting for their collisional interaction is a key point of attention in the present paper. For the baseline scenario at hand, nine ion populations interact, none of which is untouched by either the beam or the RF heating or both. The T ions are heated through ICRH and receive power from ${}^3\textrm {He}$ and, to a lesser extent, beam D and T ions slowing down. In the core, D ions are purely heated by collisions. The D heating on the far high-field side occurs in a region where drawing high-energy tails is inefficient. In the present paper the electron distribution function has not been evolved. There are two reasons for this: (i) expressed in terms of the velocity (the strength of Coulomb collisions is determined by the relative velocity) the electron distribution is very wide. Certainly in high density plasmas, ICRH is unable to deviate the electron distribution much from its Maxwellian character and it is sufficient to feed the power flowing to the electrons to the transport equation which describes how the macroscopic temperature changes under the influence of radial diffusion, local sources and losses. (ii) Because there is little overlap between the electron and the ion distributions, it takes many more iterations to reach a steady state when evolving the electron distribution function together with that of the ions (this was discussed and illustrated in Van Eester & Lerche (Reference Van Eester and Lerche2011)). Moreover, electrons easily reach speeds for which relativistic effects need to be accounted for in high temperature plasmas. The numerical grid required to solve the relevant Fokker–Planck equation for the electrons is very different from that required for the ions. One subtlety to keep in mind when omitting the changes in the electron distribution is that one cannot blindly use the electron contribution to the expressions for the collision operator on a grid suitable for the ion distributions: because most of the electron distribution is cut off, an extra contribution

(8.1)\begin{equation} D_{vv}^{a/e,\textrm{extra}}= \frac{4 {\rm \pi} \varGamma^{a/e}}{3N_e} \int_{v_{\max}}^\infty \textrm{d} pp F_{o,e}(p)=\frac{4 {\rm \pi} \varGamma^{a/e}}{3N_e} v_{\textrm{th},2}^2 F_{o,e}(v_{\max}) \end{equation}

has to be added to the diffusion coefficient of the ion collision operator when computing the operator numerically; $v_{\max }$ is the highest velocity grid point considered for computing the ion distributions. This term dominates the contribution added at lower velocities and hence omitting it represents an underestimation of the diffusion ions undergo. Whereas the convection/friction term in the collision operator forces a slowing down, diffusion allows small fractions of particles to gain energy via Coulomb interaction.

Figure 9.

Distribution functions of all ion species at $\rho =146$ mm.

Having obtained the distribution functions of all ion species on all magnetic surfaces $F_{o,i}(v,\rho )$, the local fusion power densities can be computed using the expressions for the total fusion rates accounting for the Bosch & Hale cross-sections for the various nuclear fusion reactions adopting the method described in § 7. Figure 10 depicts an example in which D–D, D–T as well as T–T fusion reactions have been accounted for, considering the actual distribution functions of the interacting species on each magnetic surface. For a given distribution, no distinction is made between the thermal part and the suprathermal part. The contribution from the bulk ions is dominant in baseline discharges with moderate $\beta$. The contributions due to the interaction between the bulk species and the beams are adding a substantial contribution as well. Note that the contribution due to the D beam interacting with the T bulk dominates that of the T beam interacting with the D bulk, suggesting that using pure D beams is more beneficial than a mix of D and T beams even when D is not directly heated. One of the reasons is that it is likely that $N=2$ heating forms a high-energy RF heated beam tail which is counterproductive for the fusion yield since it is at too high energy to add a contribution. Recall for example that the optimal equivalent Maxwellian temperature for T to contribute to the D–T reaction is 162 keV while it is 108 keV for D ions and 65 keV when suprathermal populations of both exist. In contrast, the D beam is not heated directly by the RF waves but contains a thermalised slowed down contribution with a temperature higher than that of the corresponding T. Also, because of its heavier mass, the thermalised subpopulation of the bulk T ions is shifted to lower velocities than that of the bulk D so that any accelerated population slowing down on the bulk ions is more likely to interact with the D than with the T particles. Aside from the $\textrm {D}_{\textrm {bulk}}$–$\textrm {T}_{\textrm {bulk}}$, $\textrm {D}_{\textrm {bulk}}$–$\textrm {T}_{\textrm {NBI}}$ and $\textrm {T}_{\textrm {bulk}}$–$\textrm {D}_{\textrm {NBI}}$, all other contributions to the neutron rate are negligibly small.

Figure 10.

Fusion power density profile for the case when 5 MW of RF power and 26 MW of beam power is launched into the JET plasma.

The net effect of the ICRH–NBI synergy can be seen in table 1 where the fusion yield for four different ICRH powers has been summarised. The first point to notice is that the dominant fusion yield indirectly (i.e. via Coulomb collisions) comes from the beam: when the wave power is negligible compared with the beam power, the fusion yield is already ${\approx }11\ \textrm {MW}$. Adding RF power, this number is growing for growing applied RF power through the individual increases from the various contributions (bulk D–bulk D, bulk D–beam D, $\ldots$) but the dominant contribution comes from the fusion reactions between the bulk species. The latter fusion rate is modified for two reasons: (i) the temperature of the thermal subpopulation has increased and (ii) suprathermal ion tails have been formed. For a fusion reaction with a distinct maximum near the temperature naturally reached by freshly injected beam particles, the second harmonic RF heating of the bulk as well as beam species creates tails that have energies too high to be of benefit for producing fusion reactions. The main mechanism through which the fusion reactions are triggered is the medium high and high-energy particles slowing down and boosting the temperatures of the thermal subpopulations. Through the use of ${}^3\textrm {He}$, the bulk ion temperature can be increased (${}^3\textrm {He}$ being more massive than H, it tends to slow down more on ions than H does). Also electron heating is beneficial for boosting the ion temperature: increasing the critical energy allows more of the beam power to be transferred to the ions.

Table 1.

Partial (per source pair) and total (summed over all pairs) neutron rate in $10^{17}\ \textrm {n}\ \textrm {s}^{-1}$ and corresponding total generated fusion power in MW for four different RF power levels: $0.1$, $3$, $5$ and $7$ MW.

It is important to point out that the present computations only zoom in on ICRH and NBI aspects. In practice, MHD activity, changes in transport and radiation play a key role. The goal of the EU-IM platform is exactly to incorporate all such effects as best as possible when studying scenarios. As the present paper only focussed on the description of the auxiliary heating and ended at providing source terms for the transport equations (power deposition profiles after Coulomb collisional redistribution of incoming wave and beam power) rather than actually integrating them, the numbers given are inherently connected to the starting assumptions made and hence are only partial. Aside from dynamics fully omitted here, also the parameters related to the ICRH and beam sources themselves are of importance. In the present paper a single concentration and a single minority specie was considered. Exploiting ETS and assuming $P_{\textrm {RF}}=5\ \textrm {MW}$, Huynh finds the fusion power drops from 13–14 MW when the minority concentrations are very low ($N_{\textrm {minority}}/N_e<1\,\%$) to 8–10 MW when they are increased beyond the values usually used in the experiments ($N_{\textrm {minority}}/N_e \approx 5\,\%$), and this – although details differ – both when using H or ${}^3\textrm {He}$ (but adopting the appropriate RF frequency to ensure core heating) (Huynh et al. Reference Huynh, Lerche, Van Eester, Garcia, Johnson, Ferreira, Kirov, Yadykin and Strand2020), suggesting the formation of high-energy tails is beneficial but mainly pointing to the fact that fuel dilution is counterproductive. Finally, since the bulk ion temperature is the most important parameter in determining the fusion yield, any means of modifying the bulk temperature considered will impact on the fusion yield obtained. Remembering that core ion temperatures are not available experimentally in JET, correctly predicting the bulk ion temperature is therefore crucial. The inital temperatures considered here are provided by ETS runs. Repeating the above exercise starting from a plasma that starts from a core electron temperature of 7 keV and an ion temperature of 8.5 keV – the temperatures reached only exploiting the beams – reduces the fusion achieved power from the values in the above table to $7.69$, $8.58$, $9.20$ and 9.70 MW for an added RF power of $0.1$, $3$, $5$ and 7 MW, respectively. The D–T contribution from the thermal bulk populations is roughly half of the fusion yield but is gradually increasing with the RF power level. For the adopted power levels it is $44\,\%$, $48\,\%$, $51\,\%$ and $53\,\%$, respectively. This very mild increase – even the maximum RF power considered is modest compared with that of the beams – possibly reflects the fact that the applied wave power helps raise the ion temperature, partly since more power is directly absorbed by the ${}^3\textrm {He}$ minority, which dumps a fraction of its power onto the bulk ions. More generally, wave power can collisionally end up in the bulk ion channel, the total RF power density being distributed over various species and remaining fairly modest, and a fraction of the beam energy being channelled to the ions on their slowing down (see figures 7 and 8). Studying fusion enhancement brought about by ICRH in a DT plasma but using a much more sophisticated Fokker–Planck model than the one adopted here but only considering one (RF) heated non-Maxwellian population at the time while prescribing the density and temperature, Harvey equally found a very mild increase of the enhancement when reaching higher power densities for JET-relevant densities and temperatures when accounting for the observed mild temperature increase (see figure 5 in Harvey et al. (Reference Harvey, McCoy, Kerbel and Chiu1986)). At very high power densities (not reached at JET by any of the heated species) he observes a rollover caused by the optimal operation point being passed.