RABBIT: Real-time simulation of the NBI fast-ion distribution

Knowledge of the fast-ion distribution arising from neutral beam injection (NBI) is important for transport analysis and magnetic equilibrium reconstruction. For sophisticated real-time plasma control, which will be essential for the success of future fusion devices, this distribution function needs to be known already during the discharge. Then, the relevant quantities (e.g. heating profiles, current-drive etc) can be coupled to real-time transport and equilibrium codes like RAPTOR [1–3], which has already been implemented and tested in the discharge control systems of present-day machines like TCV and ASDEX Upgrade. Beyond real-time applications, such fast models are essential for optimization problems, e.g. reactor design studies or discharge planning [2].

The NBI fast-ion distribution f can be calculated by solving the kinetic equation:

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\del}{\partial} \displaystyle \label{eq:kineq} \frac{\del f}{\del t} + (\vec {v} \cdot \nabla_{\vec {x}}\, f + \vec {a} \cdot \nabla_{\vec {v}} f) = \hat{C}(\,f) + \sigma. \nonumber \end{align} \tag{ 1 }$$

Here, the terms correspond, from right to left, to the source term (e.g. due to NBI), the collision operator (e.g. slowing down and pitch angle scattering), the orbit effects (where $\vec {a}=\frac{Z e}{m}(\vec {E} + \frac{\vec {v}}{c}\times\vec {B})$ is the acceleration due to the Lorentz force on ions with charge Ze, mass m, velocity vector $\vec {v}$ in an electric and magnetic field $\vec {E}$ and $\vec {B}$) and the partial time derivative. Several sophisticated models exist, that can calculate this beam ion distribution in good agreement with experimental data [4–6], such as the Monte-Carlo codes NUBEAM [7, 8], ASCOT [9] or NEMO/SPOT [10]. The high accuracy of these codes has, however, to be paid with relatively intensive numerical efforts, which compromises their use in real-time applications. In addition, there are reduced models like FIFPC [11], SINBAD-SSFPQL [12–14], RISK [15], CQL3D [16–18]. However, also these models are not optimized for real-time applications. For example, RISK needs  ≈30 s per time-step according to [15]). Other reduced-model codes like PENCIL [19–22] and the NBI block in ASTRA [23, 24] are faster, but still not fast enough for real-time applications. In addition, they usually make quite strong approximations, e.g. PENCIL neglects finite-orbit-width effects entirely and calculates only steady-state solutions.

In this paper, we present the novel code RABBIT (Rapid Analytical Based Beam Injection Tool). RABBIT currently takes  ≈25 ms per time step, which is roughly a factor of 1000 faster than the NUBEAM code. This execution time is comparable with the energy confinement time (e.g.  ≈100 ms on ASDEX Upgrade), which makes real-time applications possible. The approximations needed to arrive at this goal are discussed in this paper, following the terms of the kinetic equation (1) from the right to the left: a strongly simplified beam geometry is used for calculating the beam attenuation (section 2). The collisions are treated by an analytic solution of the steady-state Fokker–Planck equation in each radial cell, locally in the uniform-plasma limit (section 3). A correction for finite orbit width effects is included by an orbit average of the source term over the first fast-ion orbit (section 4). Our analytic treatment of the time-dependence is discussed in section 5. We continue our discussion with a comprehensive benchmark with NUBEAM (sections 6 and 7), and discussions of further details of the model such as the treatment of plasma rotation (section 8) and speed diffusion correction (section 9).

The calculation of the source term—the fast-ion birth rate—is equivalent to calculating the beam attenuation: each neutral atom, which is lost from the beam through ionization, becomes a newly born fast ion. For the calculation of the beam attenuation, we use BESFM (beam emission forward-model), which was developed for the analysis of beam emission charge exchange spectra [25]. The code solves a collisional radiative model, including excitation and ionization reactions with electrons, main and impurity ions and charge exchange reactions with main and impurity ions. The code calculates the neutral density evolution, resolved for different atomic states, on straight beamlets. The geometry and number of beam-lets representing one source can be modified. For our purposes, we choose the simplest geometry, which is the thin center-line of the source. In addition, we restrict the amount of simulated atomic states to $n\leqslant 2$. The restriction to such a low amount of atomic states greatly shortens the computation time, since the matrices have sizes proportional to n². To compensate for this low $\max(n)$, we multiply the losses into higher states by a factor f_last (see figure 1) [25]. It should be noted that (unlike NUBEAM) we do not treat the reactions of beam neutrals with fast-ions separately. In principle, they behave differently than the reactions with thermal bulk ions, because the relative velocities are different. We neglect this and model the fast-ions as if they were thermal. In a comparison to NUBEAM this approach still gave a good agreement of the shine-through power, even for unrealistically high central fast-ion concentrations of 80%.

Zoom In Zoom Out Reset image size

As a result of BESFM, we get the flux Γ of neutrals along the thin center-line (coordinate l) of the NBI source (in (1/s)). The shine-through power is given by the neutral flux at the end of the calculation grid (i.e. when the beam has left the plasma and will the first wall) multiplied by the energy of the neutrals $P_{\rm shine} = \sum_i \Gamma_i(\infty) E_i$ (summed over all beams and energy components).

The fast-ion birth (or deposition) rate can be calculated from the derivative along the line. To get the fast-ion birth-rate (in (1/s)) between two grid points separated by $\Delta l$, we simply take the difference $\tilde{S}_{\rm line} (l+\Delta l/2) = \Gamma(l+\Delta l) - \Gamma(l)$. We can get a birth profile S (in (m⁻³ s⁻¹)) by rebinning the grid points along the line into a radial grid (e.g. in terms of $\rho_{\rm tor}$), and by dividing with the volume of each radial grid cell ($S(\rho_i) = \tilde{S}(\rho_i)/\Delta V_i$). If we would continue to consider only the thin center-line, this would lead to bad agreement with NUBEAM calculations. This is illustrated by the dashed line in figure 2: in reality, a neutral beam has a quite significant broadening in the poloidal plane, perpendicular to the center-line. In the outer plasma, this is not so crucial for the birth profile, because the flux surfaces are tangential to this broadening. In the plasma center, however, bad agreement is observed with realistic NUBEAM calculations. E.g. in this case, the 'thin beam' model shows an almost divergent behavior, because the considered NBI source (AUG Q3) lies very close to the magnetic axis.

Zoom In Zoom Out Reset image size

In order to improve the birth-profile calculation, we need to consider the poloidal spreading of the beam. The geometry of our model is sketched in figure 3: we assume a Gaussian spreading locally perpendicular to the beam-centerline. Its standard deviation $\sigma=\sigma(l)$ is allowed to be a function of the coordinate along the center-line l. We discuss now our beam-width model at a given point B (with given R,z coordinates R_b, z_b) on the center-line. In order to be fast, we assume locally circular flux surfaces: we calculate the geometric minor radius r_b and the radial flux coordinate $\rho_b$ corresponding to $(R_{\rm b}, z_{\rm b})$ based on the realistic (numeric) equilibrium. For all other points in the poloidal plane (away from the beam center-line), we assume a linear relation:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle r(\rho) = \rho \cdot \frac{r_{\rm b}}{\rho_{\rm b}}. \nonumber \end{align} \tag{ 2 }$$

Zoom In Zoom Out Reset image size

In that sense, it is important to use a flux coordinate ρ which behaves very similar to a radius, and we have chosen $\rho = \rho_{\rm tor}$ (the square root of normalized toroidal flux) for this reason. If we want to calculate the contribution of the birth rate in point B $\tilde{S}_{\rm line}(l_{\rm b})$ into a given radial cell $\rho_i$, we need to consider the crossing-points between the orange line (perp. to the beam) and the flux surface boundaries of that radial cell: h_1l, h_2l, h_1u and h_2u. Due to our assumption of circular flux surfaces, these crossing-points can be calculated analytically by elementary geometrical considerations. The contribution into the ith radial cell $\rho_i$ is then given by:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \tilde{S} (\rho_i) = \tilde{S}_{\rm line}(l_{\rm b})\cdot w(\rho_i)/\mathfrak{s} \label{eq:sline} \nonumber \end{align} \tag{ 3 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle w(\rho_i)=\frac{1}{2}\left({\rm erf}\, \frac{h_{2u}}{\sqrt{2}\sigma} - {\rm erf}\, \frac{h_{1u}}{\sqrt{2}\sigma} \right) + \frac{1}{2}\left({\rm erf}\, \frac{h_{2l}}{\sqrt{2}\sigma} - {\rm erf}\, \frac{h_{1l}}{\sqrt{2}\sigma} \right) \nonumber \end{align} \tag{ 4 }$$

since the integral of our assumed Gaussian broadening is given by the error function. For the parts of the beam that are well inside the plasma, these contributions add up to $\mathfrak{s}=\sum_i w(\rho_i) = 1$ (at least within the numerical accuracy). When the beam-center line is close to the plasma boundary, a significant fraction of the beam-width can be outside of the plasma, and hence $\mathfrak{s}$ will be smaller than 1. In such a situation, the attenuation calculation along the center-line will over-estimate the attenuation, since in reality, a part of the beam is outside the plasma and will fly through the given grid-cell unattenuated. Further outside, a fraction $\Delta\mathfrak{f}$ of the beam may intersect with the wall (at $\rho_{\rm lim}$). We use $\mathfrak{s}$ and $\Delta\mathfrak{f}$ to correct the beam-centerline calculation for these finite-beam-width effects, and modify attenuation and shine-through accordingly:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \frac{\Gamma(l+\Delta l)}{\Gamma(l)}& \rightarrow \mathfrak{s} \frac{\Gamma(l+\Delta l)}{\Gamma(l)} + (1-\mathfrak{s}) - \Delta\mathfrak{f}~ \nonumber \\ & P_{\rm shine} \rightarrow P_{\rm shine} + \Delta\mathfrak{f}\, \Gamma(l) E \nonumber \end{align} \tag{ 5 }$$

where $\Delta l$ is the distance between two center-line grid points. In consistence with that, we use the renormalized contribution $w(\rho_i)/\mathfrak{s}$ in equation (3).

As discussed above, this beam-width correction is most important close to the plasma center. Here, the flux surfaces are also in realistic equilibria very close to a circular shape, which justifies our assumptions. Nevertheless, we do a correction for elongation of the flux surfaces. Elongation effectively means, that the length scales in figure 3 are different along the vertical and horizontal axes. This means also that the length scales along the blue and orange line are different. We can take this into account by rescaling the beam standard deviation to an 'effective' value (on each point on the grid along the beam). We use the following rescaling:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \sigma \rightarrow \sigma \cdot \sqrt \frac {(a\cos\theta)^2 + (b\sin\theta)^2 } {(a\cos\beta)^2 + (b\sin\beta)^2 }. \nonumber \end{align} \tag{ 6 }$$

This is motivated by the equation of an ellipse with the two axis a and b, defined on a given flux surface by:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle a = \left\{\begin{array}{@{}ll@{}} \max (R) - R_0 &{\rm for: }~ R> R_0 \nonumber \\ R_0 - \min (R) &{\rm for: }~ R< R_0 \end{array}\right.;\quad \quad b = \frac{\max (z)-\min (z)}{2} \nonumber \end{align} \tag{ 7 }$$

where R₀ is the major radius of the magnetic axis.

Figure 2 shows the results of this beam width correction for the birth profile. We get a very good agreement to NUBEAM simulations, which do take a fully realistic NBI geometry into account with a computationally much more expensive Monte-Carlo approach. In addition, figure 4 shows a comparison of the shine-through calculated by a 'thin beam' approach, RABBIT (i.e. BESFM including above correction for finite-beam-width), FIDASIM [26, 27] and NUBEAM in a low-density case at ASDEX Upgrade for several beam geometries (comp. figure 17). The comparison between FIDASIM and BESFM/RABBIT is interesting, because these codes have a similar collisional radiative model and use the same cross-sections. This allows to compare the center-line attenuation (BESFM/RABBIT) with a Monte-Carlo-based calculation (FIDASIM) including the fully realistic NBI geometry. We find a good agreement with FIDASIM, which indicates that taking only the center-line is an acceptable approximation for the attenuation calculation. Largest (but still acceptable) deviations are found for the two off-axis and most tangential sources Q6 and Q7. These beams enter the plasma in the most oblique angle and have a turning point in the poloidal plane close to the separatrix (compared to figure 17), where the density gradients are large. This explains the deviations, which are lowered by our above described correction for finite-beam-width.

Zoom In Zoom Out Reset image size

In this section we describe and review the analytic solution of the Fokker–Planck (FP) equation, which we solve in RABBIT. In order to make the kinetic equation (1) analytically solvable, we will neglect the orbit term in this section. A correction for orbit effects will be discussed later in section 4. Locally we assume a uniform plasma for the solution of the FP equation, namely the FP equation is solved on each radial cell with the local density and temperature, and independently from the other radial cells. The distribution function f depends then only on two velocity-space dimensions, the speed $v$ and the pitch $\xi=v_{\parallel}/v$. The third dimension would be given by the gyro-angle, but f is supposed to be symmetric with respect to that. If we write out the collision term in the limit of (thermal electron velocity) $\gg$ (fast-ion velocity) $\gg$ (thermal ion velocity), we get the Fokker–Planck equation [7]:

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\del}{\partial} \newcommand{\pitch}{\xi} \displaystyle \label{eq:FokkerPlanck} \frac{1}{\tau_{\rm s} v^2} &\frac {\del}{\del v} [(v^3+v_{\rm c}^3)\,f]+ \frac{\beta}{\tau_{\rm s}} \frac {v_{\rm c}^3} {v^3} \frac {\del} {\del \pitch} (1 - \pitch^2) \frac{\del f}{\del \pitch}\nonumber \\ & + \frac{1}{\tau_{\rm s}v^2} \frac {\del}{\del v} \left [ \left( \frac{T_{\rm e}}{m_{\rm fi}} v^2 + \frac{T_{\rm i}}{m_{\rm fi}} \frac{v_{\rm c}^3}{v}\right) \frac{\del f}{\del v} \right] = \frac{\del f}{\del t} - \sigma \nonumber \end{align} \tag{ 8 }$$

with the Spitzer time:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \tau_{\rm s} = 6.32 \cdot 10^8 \cdot \frac {A_{\rm fi} } {Z_{\rm fi}^2 \ln \Lambda_{\rm e} }\cdot \frac {(T_{\rm e}\, {\rm (eV)})^{3/2}} {n_{\rm e}\, {\rm ({cm}^{-3})}} \,{\rm s} \nonumber \end{align} \tag{ 9 }$$

the critical velocity:

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\ZZAI}{\left\langle \frac{Z_i^2}{A_i}\right\rangle} \displaystyle v_{\rm c}=5.33 \cdot 10^{4} \cdot\sqrt{T_{\rm e}\, {\rm (eV)}}\cdot {\ZZAI}^{1/3}~{\rm m~s^{-1}} \nonumber \end{align} \tag{ 10 }$$

and:

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\ZAV}{\left\langle {Z_i^2} \right\rangle} \newcommand{\ZZAI}{\left\langle \frac{Z_i^2}{A_i}\right\rangle} \displaystyle \beta =& \frac{\ZAV}{2 \ZZAI A_{\rm fi}}, \quad \ZZAI = \frac{\sum_i n_i (Z_i^2/ A_i) \ln \Lambda_i }{n_{\rm e}\ln\Lambda_{\rm e}}, \quad \nonumber \\ & \ZAV = \frac{\sum_i n_i Z_i^2 \ln \Lambda_i }{n_{\rm e}\ln\Lambda_{\rm e}} \nonumber \end{align} \tag{ 11 }$$

where the sums with index i go over the ion species. The subscript fi stands for the fast-ion (i.e. beam) species. $\ln\Lambda_{\rm e}$ and $\ln\Lambda_{i}$ denote the Coulomb logarithms for fast-ion collisions with the electrons and the ${i}$th ion species, respectively. We calculate them with the same formulas that are also implemented in NUBEAM [28–30]. For a fast ion with energy E in a magnetic field B we calculate the Coulomb logarithms with respect to a thermal bulk species α (this can be now both electrons and ions) by:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:collogs} \omega_\alpha^2=1.74 \cdot Z_\alpha^2/A_\alpha \cdot n_\alpha +9.18\cdot 10^{15} \cdot Z_\alpha^2/A_\alpha^2 \cdot B^2 \nonumber \end{align} \tag{ 12 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle v_{{\rm rel,max,}\alpha}^2=9.58\cdot 10^{10} \cdot (T_\alpha/A_\alpha + 2 E/A_{\rm fi}) \nonumber \end{align} \tag{ 13 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle r_{\rm max}=\sqrt{1/\sum_\alpha(\omega_\alpha^2/v_{{\rm rel,max,}\alpha}^2)} \nonumber \end{align} \tag{ 14 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle v_{{\rm rel,min,}\alpha}^2=9.58\cdot 10^{10} \cdot (3 T_\alpha/A_\alpha + 2 E/A_{\rm fi}) \nonumber \end{align} \tag{ 15 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle r_{{\rm min,cl,}\alpha}=0.137\,93 \cdot |Z_\alpha \cdot Z_{\rm fi}| \cdot (A_\alpha+A_{\rm fi})/(A_\alpha A_{\rm fi} v_{{\rm rel,min,}\alpha}^2) \nonumber \end{align} \tag{ 16 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle r_{{\rm min,qu,}\alpha}=1.9121\cdot 10^{-8} \cdot (A_\alpha+A_{\rm fi})/(A_\alpha A_{\rm fi} v_{{\rm rel,min,}\alpha}) \nonumber \end{align} \tag{ 17 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle r_{{\rm min,}\alpha} = \max(r_{{\rm min,cl,}\alpha}, r_{{\rm min,qu,}\alpha}) \nonumber \end{align} \tag{ 18 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \ln\Lambda_\alpha=\ln(r_{\rm max}/r_{{\rm min,}\alpha}). \nonumber \end{align} \tag{ 19 }$$

In these formulas, $T_\alpha$ and E are in keV and $n_\alpha$ are in m⁻³. $Z_{\rm e}=-1$ and $A_{\rm e}=\frac{1}{1836.1}$ are supposed to be used for the electrons. Since these formulas depend (weakly) on the fast-ion energy E, the coefficients in the FP equation (e.g. $\tau_{\rm s}$ or v_c) become weak functions of $v$ themselves. Collisions between fast ions themselves are neglected in RABBIT and NUBEAM. The effect of a finite fast-ion density is treated instead indirectly by a dilution of the bulk ions, which decreases $n_{\rm i}/n_{\rm e}$ and the critical energy as a consequence. NUBEAM uses the same approach, which is, according to [7], a good approximation even up to relatively high fast-ion concentrations of 50%. In principle, this introduces a non-linearity into the model, which we resolve by taking the fast-ion density from the previous time-step as input for the calculation of the following time-step.

This Fokker–Planck equation describes the physics of fast ions in velocity space, i.e. the collisions with thermal bulk ions and electrons. In general, the term 'Fokker–Planck' equation is used for an equation describing the evolution of a distribution function f due to drifts and diffusion. The drift terms are recognized by a single derivative of f. In our case, the first term in equation (8) is a drift term with respect to $v$: it is the slowing down term (often called also drag or friction), which leads to a drift of the distribution function towards lower velocities. The second and third term have double derivatives and therefore describe diffusion—in the pitch direction (pitch angle scattering) and the velocity direction (speed diffusion). Pitch angle scattering leads to a broadening of the pitch distribution, which is, at the birth energy, defined by the properties of the NBI source, and becomes increasingly more isotropic towards lower energies (as illustrated in figure 5). The scattering is weak for $v\gg v_{\rm c}$ (collisions with electrons dominate) and strong for $v\ll v_{\rm c}$ (collisions with ions dominate).

Zoom In Zoom Out Reset image size

In general, the terms containing v_c correspond to collisions with thermal ions and the remaining terms to collisions with thermal electrons. Pitch-angle scattering is thus a pure effect of collisions with ions, because in collisions with electrons the momentum transfer is negligible. The slowing down and speed diffusion parts of the collision operator have contributions from both types of collisions. The effect of the speed diffusion is very weak for velocities below the birth velocity, since the slowing down term dominates here. Therefore, we neglect the speed-diffusion in this section. Above the birth energy, speed diffusion leads to a high energy tail, which we discuss in section 9.

Finally, σ is the NBI source term in the FP equation. We parametrize it as

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\del}{\partial} \newcommand{\pitch}{\xi} \displaystyle \sigma = \frac {S } {2 \pi v^2} \delta(v-v_0) K(\pitch) \label{eq:NBIsource} \nonumber \end{align} \tag{ 20 }$$

which assumes a mono-energetic injection velocity $v_0$. The typically three energy components of a beam can be modeled by summation, then. A broad pitch distribution $\newcommand{\pitch}{\xi} K(\pitch)$ (e.g. due to the broadening of the beam) can be taken into account. The normalization is defined such, that $\newcommand{\pitch}{\xi} \int K(\pitch) {\rm d} \pitch = 1$ and thus $\int\int\int \sigma {\rm d}^3 \vec {v} = S$. Hence S is the fast-ion birth (or deposition) per volume and time (e.g. in units of (m⁻³ s⁻¹)).

The steady-state solution to the FP equation for such a source term has been found in [31, 32]. Since Legendre polynomials $P_l(\xi)$ are eigenfunctions of the pitch-angle scattering operator, the solution is given as a series of Legendre polynomials:

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\pitch}{\xi} \displaystyle f(v,\pitch) = \frac{1}{2 \pi}\frac {S \cdot \tau_{\rm s}}{v^3+v_{\rm c}^3} \cdot \sum_{l=0}^\infty (l+\tfrac 1 2) u^{l(l+1)} P_l(\pitch) K_l \cdot H(v_0-v) \label{fpsolution} \nonumber \end{align} \tag{ 21 }$$

with:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle K_l = \int K(\xi) P_l(\xi) \,{\rm d} \xi \label{eq:Kn} \nonumber \end{align} \tag{ 22 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle u = \left(\frac {v_0^3 + v_{\rm c}^3 }{v^3+v_{\rm c}^3} \frac {v^3}{v_0^3}\right)^{\beta/3} \nonumber \end{align} \tag{ 23 }$$

and H being the Heaviside step function. For the derivation of this solution it has been assumed that the coefficients $\tau_{\rm s}$, v_c and β are constants that do not depend on $v$, which is contradicted by our Coulomb logarithm formulas (12)–(19). We treat this by evaluating the Coulomb logarithms with respect to a reference energy E_ref: if the distribution function is to be evaluated in the interval $E_1<E<E_0$, we use the weighted average

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:Eref} E_{\rm ref} = 0.6 E_1 + 0.4 E_0 \nonumber \end{align} \tag{ 24 }$$

which is motivated by the fact that the energy dependence of $\ln\Lambda_i$ is stronger towards lower energies.

Figure 6(a) shows an example of this solution in the plasma core. We have calculated the solutions for the full, half and third energy component of the NBI independently of each other, and added them together. This is possible because in our model (and also in NUBEAM), collisions between fast ions are neglected. Figure 6(b) shows the result of a NUBEAM calculation as comparison and the shape of the solutions agree very well.

Zoom In Zoom Out Reset image size

In the end we are interested in integrals of the distribution function (e.g. fast-ion current, pressure, heating profiles etc). Many of these integrals can be solved analytically, which is ideal for real-time calculations. For example, the fast-ion density is given by:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle n_{\rm fi} = \int_{0}^{v_0} \int_{-1}^{+1} f \cdot 2\pi v^2 {\rm d} v {\rm d} \xi = \frac{S \tau_{\rm s}} {3} \ln\left({\frac{v_0^3+v_{\rm c}^3}{v_{\rm c}^3}}\right). \nonumber \end{align} \tag{ 25 }$$

In these integrals, the orthogonality relation of the Legendre polynomials is exploited:

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\del}{\partial} \displaystyle \int\limits_{-1}^{1} P_n(\xi) P_m(\xi)\, \mathrm{d}\xi = \frac{2}{2n+1} \delta_{nm}. \nonumber \end{align} \tag{ 26 }$$

Since $P_0(\xi)=1$, we can interprete the above integration over ξ as $\int {\rm d}\xi = \int P_0(\xi) {\rm d}\xi$, which lets only survive the l  =  0 term of the infinite sum. A subsequent integration over $v$ leads then to the final result.

The first few Legendre polynomials are given by:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle P_0(\xi)=1, \ P_1(\xi) = \xi, \ P_2(\xi) = \frac{1}{2} (3\xi^2 - 1), \ \dots . \nonumber \end{align} \tag{ 27 }$$

In general, $P_l(\xi)$ is a polynomial of order l. Since the integrals of interest contain usually only low powers of ξ (e.g. the fast-ion current $\propto \int\int v \xi f \cdot 2\pi v^2$), we need in fact only very few summands of the infinite sum in equation (21)—e.g. RABBIT considers currently only $l=0, 1$.

For a first comparison with NUBEAM, we consider the fast-ion density and heating to electron and ions (the latter is calculated according to equations (43) and (44)). The results are shown by the dashed lines in figure 12. While the overall agreement is reasonable, the profile shapes do not yet agree very well. All three profiles show an under-estimation in the plasma core and an over-estimation at the plasma edge. These deviations stem from the orbit effects, which we have neglected so far and which we will take into account in the next section.

Solving the full kinetic equation including the orbit term is a much more difficult task, and analytically not possible. This is the reason, why many codes such as NUBEAM rely on Monte-Carlo methods. Here, the source term is represented by a set of Monte-Carlo markers. For each marker, the orbit trajectory is calculated and during each orbit step, a Monte-Carlo collision operator is applied (changing the velocity vector). In a real-time application, such a complete treatment is not possible. We therefore have to consider an ad hoc correction for these orbit effects.

A way to do this is to average the NBI source term over the first (unperturbed) orbit. This is illustrated in figure 7 for a trapped fast-ion. Up to now, we have only considered the birth position and assumed that the ion stays on that initial flux surface. In reality, the ion travels towards other flux surfaces, and thereby also changes its pitch. In the concrete example, which is typical for co-current NBI, the birth position (marked with a circle) is almost on the furthest outside part of the orbit, and the ion moves further inside on the orbit. This leads effectively to a broadened birth distribution of this individual ion, which is broadened towards the plasma center with respect to its initial position. We can take this into account by mapping the orbit into a 2D deposition function $\tilde{S}(\rho, \xi)$, where the weight for each bin is given by the fraction of time the ion spends in the given bin (relative to the duration of a full poloidal orbit turn). If an orbit is unconfined (i.e. promptly lost), $\tilde{S}(\rho, \xi)$ is set to zero. In order to use $\tilde{S}(\rho, \xi)$ in our analytic solution of the FP equation, we further decompose it into Legendre polynomials according to equation (22) (again we need only very few moments, currently $l=0, 1$).

Zoom In Zoom Out Reset image size

The brute force approach to calculate such an orbit average for the whole NBI source would be then to take a MC representation of the source, calculate such a 2D deposition function (i.e. an orbit) for every marker, and sum over all markers [33]. To get reasonable statistics, ≈5000 markers are needed. Using a conventional 4th order Runge–Kutta integrator, roughly 1 s would be needed to calculate these orbits, which is clearly too slow for real-time applications.

To overcome this one could, in principle, either use approximation formulas for the orbits or try to reduce the number of necessary orbits. We initially tried the first approach, which however gave not fully satisfactory results. Although there are quite precise approximation formulas for e.g. the orbit widths, it is still difficult to estimate the 2D deposition function, i.e. the radial profile and the pitch evolution. In addition, a MC approach does not work well with our approach for calculating the beam attenuation along a single line (the beam center-line).

Ultimately, we decided to implement the second approach, by calculating only very few orbits and doing an appropriate interpolation in between. Figure 8 shows an example of all (19) guiding center orbits calculated for the highest energy component. We calculate only orbits along the beam attenuation grid (i.e. the beam centerline), and even here we calculate an orbit only for every nth grid point. The guiding center starting position is slightly elevated from the beam-centerline because we take the Larmor radius into account. In doing so, we need to calculate $3\cdot19=57$ orbits per injector. This is possible with a conventional 4th order Runge–Kutta integrator in  ≈10 ms (without parallelization) and hence fast enough for real-time applications. In order to get the deposition function also for the grid-points along the center-line in between (e.g. at $l_{\rm b}$), we need to do an interpolation, as illustrated in figure 9. Assuming that we have calculated the orbit deposition functions $\tilde{S}_1$ and $\tilde{S}_2$ for $l_1$ and $l_2$ with $l_1<l_{\rm b} < l_2$ we use:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:shiftedinterpol} \tilde{S}_b(\rho)\! =\! \tilde{S}_1(\rho - (\rho_b-\rho_1))\cdot \frac{l_2-l_b}{l_2-l_1} + \tilde{S}_2(\rho - (\rho_b-\rho_2)) \cdot \frac{l_b-l_1}{l_2-l_1}. \nonumber \end{align} \tag{ 28 }$$

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

This corresponds to a shift of $\tilde{S}_1$ and $\tilde{S}_2$ to the new location $\rho_{\rm b}$ and a subsequent linear interpolation. The shift is done with respect to ρ where as the interpolation is done with respect to the beam center-line coordinate l. The shift ensures, that the width of $\tilde{S}_{\rm b}$ is equal to an interpolation of the widths of $\tilde{S}_1$ and $\tilde{S}_2$, whereas a pure interpolation without shifting would lead to an over-estimation of the width of $\tilde{S}_{\rm b}$.

With this method, we are able to calculate the orbit averaging for ions born on the beam center-line. Now, we have to make this compatible with our beam-width correction (introduced in section 2). This is illustrated in figure 10. We start on the center-line at point b, which is in the radial cell $\rho_{\rm ref}$. The beam-width correction is applied on the orange straight line perpendicular to the beam center-line and we need now to extrapolate the orbit distribution function from $\rho_{\rm ref}$ to other radial cells $\rho_i$ along the orange line. We do this in a similar fashion as for the interpolation technique: we shift the distribution function from $\rho_{\rm ref}$ to the new position $\rho_i$:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:orbextra} \tilde{S}_i(\rho) = \tilde{S}_b(\rho - (\rho_i-\rho_{\rm ref})). \nonumber \end{align} \tag{ 29 }$$

Zoom In Zoom Out Reset image size

This is based on the assumption, that the orbits at $\rho_i$ and $\rho_{\rm ref}$ are similar: same width (in terms of ρ), same topology, same pitches and same starting position (e.g. in terms of poloidal angle θ). Most of these assumptions should be fulfilled quite well. Topology changes may occur at the trapped-passing boundary, and with this method we effectively 'smooth' over this boundary, which is done in reality by collisions anyway. The last assumption is well fulfilled in the outer plasma, since the width of the beam imposes only small changes of θ here. In the plasma center, this assumption may be not fulfilled well, for example if the beam center-line approaches closely the magnetic axis as in the example which we have considered earlier in figure 2. Since the orbit widths are usually small close to the magnetic axis, this should however be tolerable.

In order to test the accuracy of the RABBIT-orbit averaging technique, we compare it in figure 11 to a conventional Monte-Carlo orbit average of a FIDASIM [26, 27] birth distribution as described in section 2 of [33]. We get very good agreement in the birth profile, and also a good agreement in the average pitch (which is the first Legendre component). Also, the RABBIT-orbit average produces smoother profiles (due to the interpolations), while the MC profiles contain Monte-Carlo noise. At the same time, the RABBIT-orbit average is orders of magnitude faster than the Monte-Carlo approach, since it requires only  ≈60 orbit calculations compared to the 5000 orbit calculations in the MC calculation.

Zoom In Zoom Out Reset image size

It should be noted that we have observed some cases, where the inner-most radial point of the birth profile (i.e. closest to the magnetic axis) shows an unphysical jump with respect to the rest of the profile. This could be caused because our assumptions (see above) are not so well fulfilled close to the magnetic axis, or because the volume of the inner-most cell is very small, such that it is covered only by very few orbits and is thus more prone to suffer from numerical noise. To cure this, we have implemented an upper limit of the gradient between the first and second profile point (based on a fraction of the maximal gradient in the rest of the birth profile). If this limit is violated, we smooth the spike in the first point of the birth profile by redistributing into the neighboring cells under conservation of the volume-integrated birth rate. Since the volume element of the neighboring cell grows very fast (roughly $\propto \rho_{\rm tor}$, i.e. the second cell has roughly the triple volume of the first cell), their absolute values are much less affected by this smoothing procedure than the first radial point, where the spike is being removed.

As final step, we plug the orbit-averaged source term in the FP equation and calculate the profiles (figure 12). The orbit averaging leads to a very good agreement with the profiles from NUBEAM. Slight deviations remain in the plasma center, affecting only a small fraction of the plasma volume. In figure 13 we compute the volume-integrated heating power to electrons and ions. It is interesting to note, that the orbit average has also an effect on these volume-integrated quantities and improves the agreement with NUBEAM. This is due to the (mostly) linear dependence of the critical energy on T_e: since we consider co-current NBI in this exampe, the orbits shift the fast-ion distribution further inwards—i.e. away from the cold plasma edge (where the critical energy is lower and electron heating more dominant) towards the hot plasma core (where the critical energy is higher and electron heating is thus weaker).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

It should be noted that NUBEAM stops following the Monte Carlo markers when their energy has reached 3/2 of T_i, while in this paper the integrals of f are calculated all the way down to $v=0$. In NUBEAM, the remaining heating power after stopping the MC markers is stored as 'thermalization power'. To ensure a fair comparison between RABBIT and NUBEAM, we add this thermalization power to the NUBEAM ion-heating (since one can expect that at such low energies one is well below the critical energy such that the slowing down is mainly on the background ions). In general, RABBIT has an option to set the lower integration boundary to an user-defined multiple of T_i.

In the previous sections, we have discussed the steady-state solution of the Fokker–Planck equation, which we will label now f_ss. In order to treat the time-dependence, we follow the approach introduced in [32]. Using a Laplace-transform and neglecting speed diffusion, he derived the following time evolution of the distribution function:

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\pitch}{\xi} \displaystyle f(v,\pitch,t) = \frac {\tau_{\rm s}}{v^3+v_{\rm c}^3} \cdot \sum_{l=0}^\infty \left(l+\frac 1 2\right) \left(\frac{v^3}{v^3+v^3_{\rm c}}\right)^{l(l+1)\beta/3} P_l(\xi) \cdot \nonumber \\ \int_v^{\infty} S_l\left(v', t-\frac{\tau_{\rm s}}{3} \ln\frac{v'^3+v_{\rm c}^3}{v^3+v_{\rm c}^3}\right) \cdot \left(\frac{v'^3+v_{\rm c}^3}{v'^3}\right)^{l(l+1)\beta/3} v'^2 {\rm d} v' .\label{fptime} \nonumber \end{align} \tag{ 30 }$$

The key ingredient here is, that thanks to the neglection of speed diffusion, one can calculate the time τ which is needed to slow a fast ion down from $v'$ to $v$:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \tau(v',v) = \frac{\tau_{\rm s}}{3} \ln\frac{v'^3+v_{\rm c}^3}{v^3+v_{\rm c}^3} . \label{eq:tau} \nonumber \end{align} \tag{ 31 }$$

This function appears in equation (30), meaning that ions that have been injected with an injection velocity $v'$ a time $\tau(v', v)$ ago, will appear now with a lower velocity $v$.

It should be noted, that this time evolution still assumes a constant background plasma—and only a time-dependence of the source $S_l(v, t)$ (i.e. injector power) is taken into account. We will now discretize our time steps, and assume that plasma parameters are constant during a time step $\Delta t=t_1-t_0$. This allows us to use the above equation. Furthermore, we assume that the injected power during that time step is a δ-function pulse at the beginning of the time-step: $\newcommand{\del}{\partial} S_l(v, t) = \frac{S K_l}{2\pi v^2 }\delta(v-v_0)\delta(\frac{t-t_0}{\Delta t})$, and we will average over the time-step later. This has the advantage that we only have to track the slowing down of a single fast-ion pulse per time step (according to equation (31)), rather than a continuous set. The δ-function is normalized such, that it yields a constant birth rate of S if it is averaged over the time-step and integrated over ${\rm d}^3 v$.

Plugging this into the above equation yields the following:

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\del}{\partial} \newcommand{\pitch}{\xi} \displaystyle f(v,\pitch,t) =& ~\frac{1}{2 \pi}\frac {S \cdot \tau_{\rm s}}{v^3+v_{\rm c}^3} \cdot \sum_{l=0}^\infty \left(l+\frac 1 2\right) \left(\frac{v^3}{v^3+v^3_{\rm c}}\right)^{l(l+1)\beta/3} K_l P_l(\xi) \nonumber \\ & \cdot \delta\left(t-t_0-\frac{\tau_{\rm s}}{3} \ln\frac{v_0^3+v_{\rm c}^3}{v^3+v_{\rm c}^3}\right) \cdot \left(\frac{v_0^3+v_{\rm c}^3}{v_0^3}\right)^{l(l+1)\beta/3} \cdot\Delta t \nonumber \\ =&~ f_{\rm ss}(v,\pitch) \cdot \delta\left(t-t_0-\frac{\tau_{\rm s}}{3} \ln\frac{v_0^3+v_{\rm c}^3}{v^3+v_{\rm c}^3}\right) \cdot \Delta t .\label{fptimepulse} \nonumber \end{align} \tag{ 32 }$$

This is equivalent to a δ-function-like fast-ion pulse, which moves in time to lower energies. The pulse stays a δ-function, because we have neglected speed diffusion and thus all fast ions slow down at exactly the same rate. Averaging over the duration of the time-step $\int_{t_0}^{t_1} {\rm d} t / \Delta t$ yields then:

$$\begin{align} \newcommand{\bi}{\boldsymbol} \newcommand{\e}{{\rm e}} \newcommand{\pitch}{\xi} \displaystyle f(v,\pitch)\big|^{t_1}_{t_0} = f_{\rm ss}(v,\pitch) \cdot H\left(t_1-t_0-\frac{\tau_{\rm s}}{3} \ln\frac{v_0^3+v_{\rm c}^3}{v^3+v_{\rm c}^3}\right) \nonumber \end{align} \tag{ 33 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\pitch}{\xi} \displaystyle = f_{\rm ss}(v,\pitch) \cdot H\left(v^3 + v_{\rm c}^3 - (v_{\rm 0}^3 + v_{\rm c}^3) \cdot \exp\frac{-3 \cdot \Delta t}{\tau_{\rm s}} \right) \nonumber \end{align} \tag{ 34 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\pitch}{\xi} \displaystyle \equiv f_{\rm ss}(v,\pitch) \cdot H(v - v_1). \label{fptimesol} \nonumber \end{align} \tag{ 35 }$$

In the last step we have called the velocity $v$, for which the argument of the Heaviside function vanishes, $v_1$. $v_1$ is then the velocity, to which the fast-ions have been slowed down during the time-step $\Delta t = t_1-t_0$, starting from a velocity $v_0$ at time t₀.

This result is illustrated in figure 14. In the given example, the fast ions have slowed down from 59 keV to 48 keV during the 10 ms time-step. The resulting distribution function (averaged over the time-step) is hence the normal steady-state distribution function f_ss, multiplied by a Heaviside box function. If we want to calculate moments of $\newcommand{\pitch}{\xi} \newcommand{\bi}{\boldsymbol}f(v, \pitch)\big|^{t_1}_{t_0}$, we can still use the same analytic integrals as described in section 3: the pitch integral is exactly the same, and velocity integrals evaluate as:

$$\begin{align} \newcommand{\bi}{\boldsymbol} \newcommand{\e}{{\rm e}} \newcommand{\pitch}{\xi} \displaystyle \int_0^{\infty} f(v,\pitch)\big|^{t_1}_{t_0} {\rm d} v = \int_{v_1}^{v_0} f_{\rm ss}(v,\pitch) {\rm d} v = \int_0^{v_0} f_{\rm ss}(v,\pitch) {\rm d} v - \int_0^{v_1} f_{\rm ss}(v,\pitch) {\rm d} v. \nonumber \end{align} \tag{ 36 }$$

Zoom In Zoom Out Reset image size

As it turns out, the above is all we need for a consistent time-dependent solution: at the end of a time step, we can use the its final state $f(v, \xi, t_1) = f_{\rm ss}(v_1, \xi)$ as new, 'virtual' source for the next time-step. The source for time-step (2), going from t₁ to t₂ is then defined by the following values:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle S^{(2)} = S^{(1)} \nonumber \end{align} \tag{ 37 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle v_0^{(2)} = v_1^{(1)} \nonumber \end{align} \tag{ 38 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle K^{(2)}(\xi) = f_{\rm ss}^{(1)}(v_1,\xi) / \int f_{\rm ss}^{(1)}(v_1,\xi') {\rm d} \xi' {\rm, or~equivalently:} \nonumber \end{align} \tag{ 39 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle K_l^{(2)} = u^{(1)\, l(l+1)}|_{v=v_1^{(1)}} K_l^{(1)} = \left(\frac {v_0^{3} + v_{\rm c}^3 }{v_1^3+v_{\rm c}^3} \frac {v_1^3}{v_0^3}\right)^{\frac{\beta}{3}l(l+1)} K_l \nonumber \end{align} \tag{ 40 }$$

where we have omitted the ⁽¹⁾ in the last step for brevity.

Equation (37) means that we keep the source strength constant, i.e. we assume there are no losses in between consecutive time-steps. In the future, a modification of this assumption could allow one to model e.g. charge exchange losses or radial diffusion (=losses into other radial bins). As new birth energy we use the end-point of slowing down in the previous time-step, and as new pitch distribution we use the end-point of the pitch-angle scattering in the previous time-step—i.e. $K^{(2)}(\xi)$ will in general be broader and more isotropic than $K^{(1)}(\xi)$.

This scheme is illustrated in figure 15, where the time-steps correspond to different columns. If the beam is still turned on in time-step (2), then we simply need to add a new 'row', with the (current) nominal birth energy and birth pitch distribution. The full solution at each time-point is then given by a sum over the rows. This scheme continues then, e.g. during each time-step where the beam is turned on, a new row gets added. A row gets terminated, once its fast-ion pulse has reached $v=0$.

Zoom In Zoom Out Reset image size

If the plasma parameters and beam source are constant in time, then we ultimately reach the steady-state solution, such that our time-dependent solution is consistent with that. The solution which we get from this approach is in fact slightly more precise, because we evaluate the analytic solutions only in small velocity-intervals. Thus we can take the weak velocity-dependence of the Coulomb-logarithm into account by evaluating it for each interval individually (using equation (24)).

In addition, we can treat changes of the source and/or plasma parameters consistently. The latter is demonstrated in figure 16, where we have introduced a strong change of n_e and T_e from step 1 to step 2. We see that we have now two different steady-state solutions $f_{\rm ss}^{\rm (2a)} \neq f_{\rm ss}^{\rm (2b)}$, because the fast ions in the upper row (2a) have had a different 'past history'. They have endured e.g. weaker pitch angle scattering in step 1, and the resulting effects on the distribution function are taken into account consistently. In that sense our approach is more general than previous attempts [34], where this effect is neglected. It must be noted that when we add together all partial solutions for one time point (i.e. the rows), we can get discontinuities in the distribution function. Since the goal of RABBIT is to provide mainly integrals of the distribution function, this is not of concern here. For future applications, where e.g. the gradients of the distribution function would be of interest, another numerical method should be implemented.

Zoom In Zoom Out Reset image size

For a first comparison of our time-dependent solution with NUBEAM we will consider the heating power. We calculate the heating density rates (units W m⁻³) by an integral over the slowing-down (sd) part of the collision operator:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle P_{\rm tot} = -\int \frac{m v^2}{2} \hat{C}_{\rm sd}(\,f) {\rm d}^3 \vec {v}. \nonumber \end{align} \tag{ 41 }$$

Here, m is the mass of the fast-ion species. The collision operator can be rewritten in terms of a collisional flux $\hat{C}(\,f) = -\nabla \cdot \vec {\Gamma}_{\rm c}$ [33], with $\vec {\Gamma}_{\rm c, sd} = - \frac {1}{\tau_{\rm s}} \frac{v_{\rm c}^3 +v^3}{v^2} f \hat{e}_v$. This and a subsequent partial integration allows us to rewrite:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle P_{\rm tot} &= \int \frac{m v^2}{2} \nabla \cdot \vec {\Gamma}_{\rm c,sd} {\rm d}^3 \vec {v} = - \int \vec {\Gamma}_{\rm c,sd} \cdot \nabla \left(\frac{m v^2}{2}\right) {\rm d}^3 \vec {v}\nonumber \\ & = \int mv |\vec {\Gamma}_{\rm c,sd}| {\rm d}^3 \vec {v}. \nonumber \end{align} \tag{ 42 }$$

The integration over ξ lets survive only the l  =  0 summand from the Legendre expansion of f. Using our time-dependent solution from equation (35) we get for the total heating density:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:ptot} P_{\rm tot} &= \int_{v_1}^{v_0} mv \cdot \frac{v_{\rm c}^3 +v^3}{\tau_{\rm s} v^2} \cdot \frac{\tau_{\rm s} S}{v_{\rm c}^3 +v^3} \cdot v^2 {\rm d} v\nonumber \\ & = \int_{v_1}^{v_0} Smv {\rm d} v = S \left(\frac{m v_0^2}{2}- \frac{m v_1^2}{2}\right). \nonumber \end{align} \tag{ 43 }$$

This is equal to the kinetic energy that the fast-ion pulse looses in slowing down from $v_0$ to $v_1$ and thus the intuitively expected result. For the heating to the ions we shall consider only the part of the collision operator corresponding to collisions between fast-ions and thermal ions (i.e. the terms containing v_c):

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:pi} P_{\rm i} = \int_{v_1}^{v_0} mv \cdot \frac{v_{\rm c}^3}{\tau_{\rm s} v^2} \cdot \frac{\tau_{\rm s} S}{v_{\rm c}^3 +v^3} \cdot v^2 {\rm d} v = Smv_{\rm c}^3 \int_{v_1}^{v_0} \frac{v}{v_{\rm c}^3 +v^3} {\rm d} v \nonumber \end{align} \tag{ 44 }$$

$$\begin{align} \newcommand{\bi}{\boldsymbol} \newcommand{\e}{{\rm e}} \displaystyle = &\frac{Smv_{\rm c}^2}{6} \left.\bigg (\ln(v_{\rm c}^2-v_{\rm c}v + v^2) - 2\ln(v_{\rm c}+v) \right. \nonumber \\ & \left. +2\sqrt{3}\arctan\frac{2v-v_{\rm c}}{\sqrt{3}v_{\rm c}} \right.\bigg)\bigg|^{v=v_0}_{v=v_1}. \nonumber \end{align} \tag{ 45 }$$

With $v_1=0$, these results are identical to [35]. The electron heating is then given by $P_{\rm e} = P_{\rm tot} - P_{\rm i}$.

To compare with NUBEAM, we have chosen the ASDEX Upgrade discharge #33856 where seven (out of eight) NBI sources are turned on alternately. This allows also a benchmark under a large variation of beam geometries (see figure 17). The resulting timetraces of (volume-integrated) heating to ions and total heating are shown in figure 18. We can see, that we find a very good agreement of the temporal evolution.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

It should be noted that RABBIT currently considers three power loss channels: the first is the shine-through, which has already been benchmarked in section 2 (for a lower n_e case). Secondly, prompt orbit losses are calculated during the orbit-averaging. The user provides the coordinates $(R_{\rm lim}, z_{\rm lim})$ of a point on the limiter, which is calculated into a flux radius $\rho_{\rm lim}$. If a guiding center orbit (+gyro-radius) exceeds $\rho_{\rm lim}$ it is considered as being lost. In this benchmark case, orbit losses are however insignificant (in agreement both in NUBEAM and RABBIT) due to co-current NBI. Finally, the third loss channel is plasma rotation, which will be discussed in section 8.

As next step for comparison, we consider the toroidal current carried by the fast ions. The related moment of the distribution function is the toroidal fast-ion current given by:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle j_{\rm fi} = Z_{\rm fi} e \cdot \int\!\!\int v\xi f(v, \xi) \cdot 2 \pi v^2 {\rm d} v {\rm d} \xi \nonumber \end{align} \tag{ 46 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle = S \tau_{\rm s} K_1 \int \frac {v^3}{v^3+v_{\rm c}^3} \left(\frac {v_0^3 + v_{\rm c}^3 }{v^3+v_{\rm c}^3} \frac {v^3}{v_0^3}\right)^{\frac{2 \beta}{3}} {\rm d} v. \nonumber \end{align} \tag{ 47 }$$

The ξ-integral can be evaluated analytically and only the l  =  1 Legendre components survives from the infinite sum, because $\xi = P_1(\xi)$. The remaining $v$-integral has an analytic solution as well, which however involves the hypergeometric function. Since the evaluation of such special functions is typically more complex than a numeric integration, we solve the $v$-integral numerically using Simpson's rule. Figure 19 shows a comparison with NUBEAM of the timetraces and two profiles for Q2 (intermediate on-axis) and Q7 (off-axis). All curves correspond to the fast-ion current (without electron shielding), except for the two lower time-traces in the top row, which include electron shielding (as explained later in this section). The overall shapes of the profiles agree quite well. In the absolute values of the area-integrated total current, we see some significant deviations. If we calculate the ratio between RABBIT and NUBEAM values, we can see that the RABBIT values typically are within  ±15% of the NUBEAM values, which can still be seen as a reasonable agreement (considering that RABBIT is roughly a factor 1000 faster than NUBEAM).

Zoom In Zoom Out Reset image size

The reason, why the deviations between RABBIT and NUBEAM are much higher for the fast-ion current than for other quantities, is that the signed value of $v_{\parallel}$ enters into the integral. This means that knowledge of the orbits (also during slowing down) is very important. For example, trapped particles contribute only very little to the total current, because positive $v_{\parallel}$ on the outer banana leg and negative $v_{\parallel}$ on the inner banana leg largely cancel each other. In RABBIT, we do calculate the initial orbits of the fast ions (for the orbit average of the source). However, we neglect orbit changes during the slowing down. In the top row of figure 19, we have checked whether a more sophisticated Monte-Carlo orbit average of a realistic 3D FIDASIM NBI birth distribution (as explained in section 2 of [33]) would improve the agreement. It gives, however, very similar results as RABBIT. This demonstrates again (as already shown in figure 11) that the fast orbit-averaging technique implemented in RABBIT is in good agreement with a conventional Monte-Carlo orbit-average and the deposition calculated with the simplified beam geometry in RABBIT is in good agreement with a fully realistic one. Instead, the deviations arise mainly from the concept of orbit-averaging itself.

From figures 19 and 20 we can assess that the deviation between RABBIT and NUBEAM is a function of the beam inclination with respect to the magnetic field. Very radial sources (like Q4) are under-estimated by RABBIT while very tangential sources (like Q6 and Q7) are over-estimated. In the following, we will try to discuss reasons for this dependency: for very tangential sources (like Q6 and Q7), the initial orbits are mostly passing. Some of these orbits will later change their topology (e.g. due to pitch angle scattering), become trapped orbits and can no longer significantly contribute to the fast-ion current. Even if they may become passing again in subsequent collisions, they now can become co- or counter-passing with approximately equal probability. This means that once an ion has moved from passing to trapped, it looses its ability to contribute strongly to the fast-ion current. Since we neglect this transport over the trapped-passing boundary, RABBIT tends to over-estimate the current for very tangential sources. For very radial sources, this effect does not occur because many fast ions are already initially trapped when they are born. The under-estimation of these sources could be explained by the neglected speed-diffusion: it would create a high energy tail above the injection energy, which would increase the current. In that sense, the simultaneous neglection of both these effects seems to be beneficial, because they tend to cancel each other (most efficiently for intermediate beams like Q3 or Q8.

Zoom In Zoom Out Reset image size

In the future, the above findings could allow a further improvement of the model: either a suitable ad hoc modeling of the above physics can be found, or, if this is not feasible, one could also think of a correction factor as function of the beam inclination (based on the observation in figure 20). This would, however, first require an investigation of this dependence under a larger variation of plasma parameters.

7.1. Electron shielding current

To calculate the current driven by the neutral beams (NBCD), we need to take into account that a part of the fast-ion current is shielded by electrons following the fast ions. This is done by a factor $\newcommand{\e}{{\rm e}} \eta_{\rm nbcd}$ such that the driven current is given by $\newcommand{\e}{{\rm e}} j_{\rm nbcd} = \eta_{\rm nbcd} \cdot j_{\rm fi}$. Various models and formulas exist for $\newcommand{\e}{{\rm e}} \eta_{\rm nbcd}$. In RABBIT, we have implemented $\newcommand{\e}{{\rm e}} \eta_{\rm nbcd} = 1- Z_{\rm fi}/Z_{\rm ef\,\!f}\cdot(1-L_{31})$, where L₃₁ comes from neoclassical transport theory [36, 37]. For L₃₁ we use a formula that was obtained by Sauter et al [38, 39] by fitting numerical results from the CQLP code:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle L_{31} = (1 + \frac{1.4} {Z+1}) X - \frac{1.9}{Z+1} X^2 + \frac{0.3}{Z+1}X^3 + \frac{0.2}{Z+1}X^4 \quad ([38]14a) \end{align} \tag{ 48 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle X = \frac{f_{\rm t}} {1+(1-0.1f_{\rm t}) \sqrt{\nu^*_{\rm e}} + 0.5 (1-f_{\rm t}) \nu^*_{\rm e}/Z} \qquad\quad ([38]14b) \end{align} \tag{ 49 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \nu^*_{\rm e} = 6.921\cdot10^{-18} \frac{|q| R n_{\rm e} Z_{\rm ef\,\!f} \ln\Lambda_{\rm e} }{T_{\rm e}^2 \epsilon^{3/2}} \qquad ([38]\ 18b\ {\rm with}\ Z = Z_{\rm eff}) \end{align} \tag{ 50 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle\begin{array}{l} \ln\Lambda_{\rm e} =& \ln\left(1.09\cdot10^{11} \frac{T_{\rm e}}{Z_{\rm ef\,\!f}\sqrt{n_{\rm e}/10^6}}\right) \\ & ([38]\ 18e\ \text{with additional}\ Z_{\rm eff}\ {\rm dependence}). \end{array}\end{align} \tag{ 51 }$$

In these equations, n_e is in m⁻³, temperatures are in eV, R  = (R_max  +  R_min)/2 and $\newcommand{\e}{{\rm e}} \epsilon = \frac{R_{\rm max}-R_{\rm min}}{R_{\rm max}+R_{\rm min}}$ is the (local) inverse aspect ratio. The above equations correspond to the NUBEAM implementation with option nmcurb  =  4. Differences to the original Sauter paper [38] are outlined above. The Sauter formulas have been obtained from one-species calculations, and it is therefore a controversial issue [40], which Z is to be used in equations (48) and (49). We follow the NUBEAM implementation (and the original suggestion from Sauter [39]) and use $Z=Z_{\rm ef\,\!f}$ in (48) and (49).

Above formulas include the trapped electron fraction f_t correction by Lin-Liu [36] (NUBEAM nmcurb  =  3) and take, in addition, the collisionality of the electrons $\nu^*_{\rm e}$ into account (which in general increases the electron shielding and thus decreases the current-drive efficiency).

It should be noted that Hager et al [41] introduce further correction factors $\gamma_{3i}$ for the L_3i. These corrections are however very small for L₃₁ ($|\gamma_{31}-1| \lessapprox 10^{-3}$) such that we stick with the Sauter formulas for the sake of simplicity.

The final question to clarify is how we calculate the trapped electron fraction f_t. For a given flux surface, it is defined as:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle f_{\rm t} = 1 - \frac{3}{4} \left\langle B^2\right\rangle \int_0^{1/B_{\rm max}} \frac{\lambda {\rm d} \lambda}{\left\langle\sqrt{1-\lambda B}\right\rangle} \label{eq:ftdef} \nonumber \end{align} \tag{ 52 }$$

where the integration ${\rm d} \lambda$ and averages $\left\langle \right\rangle$ are to be taken along the flux surface. This is a numerically expensive double integral. A very sophisticated estimate can be calculated from a single integral [42]. For our real-time purposes, we need to simplify further and use:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle f_{\rm t} = 1.462\,4256 \cdot\epsilon^{0.5} - 0.46 \cdot\epsilon^{1.5} \label{eq:ftapp} \nonumber \end{align} \tag{ 53 }$$

which is also used in e.g. [37, 43] (we use however more decimals for the first term). With $\newcommand{\e}{{\rm e}} \epsilon = \frac{R_{\rm max}-R_{\rm min}}{R_{\rm max}+R_{\rm min}} \equiv \epsilon_0$, this formula is a good approximation for elliptic flux surfaces in the zero-pressure limit (i.e. with zero Shafranov shift). In figure 21, we have compared this approximation formula (53) with the exact numeric solution (52) for a range of ASDEX Upgrade equilibria (calculated with the CLISTE code). We find a quite significant over-estimation of f_t which increases with the Shafranov shift. This is essentially a geometric effect: due to the Shafranov shift, the low field side of the plasma becomes smaller and the high field side larger. Thus the trapped particle fraction is reduced, since trapped particles occur dominantly on the low field side. We have found a heuristic correction, to take these finite-β effects into account: if we use $\newcommand{\e}{{\rm e}} \epsilon=\epsilon_{\rm ef\,\!f}$ with:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \epsilon_{\rm ef\,\!f} = \frac{R_{\rm max} - R_0}{R_0} \nonumber \end{align} \tag{ 54 }$$

where R_max is the maximal major radius on the flux surface and R₀ the major radius of the magnetic axis, we obtain a much better agreement to the exact calculation, such that we have implemented this in RABBIT. For zero β, $\newcommand{\e}{{\rm e}} \epsilon_{\rm ef\,\!f}$ is equal to the conventional inverse aspect ratio $\newcommand{\e}{{\rm e}} \epsilon_0$, while for finite-β, $\newcommand{\e}{{\rm e}} \epsilon_{\rm ef\,\!f}$ evaluates the minor radius solely on the low-field-side (i.e. $\newcommand{\e}{{\rm e}} \epsilon_{\rm ef\,\!f}<\epsilon_0$), which is more relevant for the trapped particle fraction and therefore leads to a better agreement.

Zoom In Zoom Out Reset image size

The lower curves in the top row of figure 19 show a comparison of the driven current between RABBIT and NUBEAM. We observe the same behavior as with the fast-ion current—which is very reasonable, since we use the same formulas for the current drive efficiency as NUBEAM.

In this section, we want to describe how toroidal plasma rotation $\omega_{\rm tor}$ is included in the RABBIT model. Plasma rotation affects the fast ions from NBI in various ways. The beam attenuation is slightly changed, because the relative velocities between beam neutrals and thermal plasma ions or electrons are changed. This effect has been easily implemented in our beam-attenuation calculation.

We move now to the rotation effects for slowing down and orbit calculations. The toroidal plasma rotation motivates to introduce two different frames of reference: the lab frame, which is an inertial frame at rest with respect to the torus hall, and a plasma frame, which rotates along with the plasma. Our analytic solution of the FP equation (21) does not include any terms with respect to $\omega_{\rm tor}$—therefore this solution must be evaluated in the plasma frame. In contrast, the NBI injector is located at rest in the lab frame. This makes necessary a transformation of the birth velocity vector (e.g. of the injection energy and pitch):

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle v_{\parallel}^{\rm (lab)} = v_{\parallel}^{\rm (\,plasma)} + R \omega_{\rm tor} {\rm sgn}{I_{\rm p}} \nonumber \end{align} \tag{ 55 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle v_\perp^{\rm (lab)} = v_\perp^{\rm (\,plasma)}. \nonumber \end{align} \tag{ 56 }$$

The additional factor containing the sign of the plasma current I_p arises from the fact that we define the parallel direction with respect to the plasma current, and the direction of $\omega_{\rm tor}$ and I_p with respect to the toroidal angle φ. In typical situations, the toroidal plasma rotation has the same direction as the NBI injection (because it is mainly driven by that). In that case the injection energy is effectively reduced in the plasma frame (e.g. from 60 keV to 50 keV for strong rotation), which leads to weaker plasma heating, current drive etc.

In addition, the orbit calculations (necessary for the orbit-average of the source term) are affected by toroidal plasma rotation. In the lab frame, a given plasma rotation with velocity vector $\vec {v}_{\rm rot}$ creates a radial electric field $\vec {E}_{\rm r}$, which can be calculated by solving the radial force balance equation:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \vec {E}_{\rm r} = \frac{\nabla p_\alpha}{n_\alpha Z_\alpha e} - \vec {v}_{{\rm rot,}\alpha} \times \vec {B}. \nonumber \end{align} \tag{ 57 }$$

For our purposes, we neglect the pressure-gradient term and the possibly different rotation velocities among different plasma species α and the poloidal plasma rotation. This allows us to rewrite E_r:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \vec {E}_{\rm r} = - \vec {v}_{\rm rot} \times \vec {B} = - \frac{\omega_{\rm tor}\nabla\varphi}{(\nabla\varphi)^2} \times (\nabla\varphi \times \nabla \Psi + F \nabla\varphi) = \omega_{\rm tor} \nabla\Psi \nonumber \end{align} \tag{ 58 }$$

with Ψ being the reduced poloidal flux (in Wb/rad). $\omega_{\rm tor}$ is assumed to be a flux surface constant, i.e. a function of Ψ (or $\rho_{\rm tor}$). In guiding-center orbit calculations, $\vec E_\text{r}$ enters through an additional $\vec {E} \times \vec {B}$-drift term. In terms of the constants of motion, the potential $V= -\int \omega_{\rm tor}{\rm d}\Psi$ is of relevance: since $\omega_{\rm tor}$ is typically a radial profile, the plasma rotation creates a radial potential well. Since the total energy (i.e. the sum of kinetic and potential energy) is a constant of motion, the kinetic energy of an ion may change along an orbit due to the potential well.

To study this effect, we consider ASDEX Upgrade discharge # 30809 at 4.48 s, which features a rather strong core rotation frequency of 24.2 kHz. Figure 22 shows the NUBEAM fast-ion distribution calculated for an off-axis position on the low field side. This case has only one active NBI source (AUG Q3 with 58.4 keV). There are three beam injection peaks visible at this plot position: from highest to lowest pitches, they arise from deposition on the high-field-side, deposition on the low-field-side (at the plot position) and deposition into a trapped orbit further outside on the low field side. Despite the fact that the distribution function is given in lab-frame coordinates, only the middle peak is located at the nominal injection energy. The ions deposited at the high-field-side have gained kinetic energy on their orbit to the plot position (because they have been born at a more central flux surface), and the trapped orbits have lost kinetic energy (we see them here at the inner leg of their banana).

Zoom In Zoom Out Reset image size

In RABBIT, we have implemented the orbit-average in the plasma frame. This invokes a Lorentz-transformation of the electric and magnetic fields. In the non-relativistic limit ($\omega_{\rm tor}R \ll c$), only $\vec {E}$ undergoes a transformation:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \vec {E}_{\rm r}^{\rm (\,plasma)} = \vec {E}_{\rm r}^{\rm (lab)} + \vec {v} \times \vec {B} = 0. \nonumber \end{align} \tag{ 59 }$$

Within our approximations, the Lorentz-transformation exactly cancels out the electric field in the plasma frame. This is beneficial for the calculation time, since we do not need an additional $\vec {E} \times \vec {B}$-drift term. However, the fact that the plasma frame is not an inertial frame gives rise to fictitious forces (namely Coriolis, centrifugal and Euler force). Both a back-of-the-envelope estimate and a numerical comparison (figure 23) indicate, however, that the drifts resulting from these fictitious forces are weak compared to the gradient-B and curvature drifts and hence can be neglected.

Zoom In Zoom Out Reset image size

Since we assume that the electric field vanishes in the plasma frame, the potential well vanishes as well such that the kinetic energy remains constant along ion orbits. We only need to do one transformation for the birth energy, from the nominal injection energy in the lab frame to the corresponding energy in the plasma frame at the birth position. This introduces however a conceptual problem: it gives rise to different birth energies on the same flux surface, because ions deposited on high or low-field-side may have different birth energies, and the orbit average introduces a further mixing of different birth locations onto a given radial cell. In our analytic solution of the FP equation we have, in contrast, assumed a mono-energetic source. We solve this by taking only the average birth energy per radial cell into account. In doing that, we assume that the real birth energies show only a small variation in a radial cell, which can be justified because the gradients of the rotation profile are small.

To benchmark our treatment of plasma rotation, we choose the time-point 3.48 s from our previous benchmark case (AUG #33856). We set up a simulation run, where all plasma parameters (including equilibrium) are kept constant to this reference point, except for the toroidal plasma rotation. Here, we use the experimentally measured profile from #30809 at 4.48 s (figure 24(a)) and scale it with a time-dependent factor, which starts from 0 to 2 in several steps (figure 24(b)). Increasing time means thus increasing plasma rotation. Figures 24(c) and (d) show then a comparison of heating powers and fast-ion current. Both reduce significantly with increasing rotation, because the fast-ion birth energy becomes lower in the rotating plasma frame. We find excellent agreement with NUBEAM for the heating powers and reasonable agreement for the fast-ion current.

Zoom In Zoom Out Reset image size

In the RABBIT calculations of the fast-ion current, a distinct step appears at 3 s when using our default 20 orbits per energy component for the orbit-average. This step is not visible in the volume-integrated heating powers. It occurs, because at that time, one of the 20 orbits switches from passing to trapped (because the plasma rotation lowers the parallel velocity in the plasma frame). This gives an indication about the uncertainty introduced by the low number of orbits into the current calculation. In a comparative RABBIT run with 150 orbits the curves are smooth and the step does not occur, the overall trend and results however are the same.

In this section, we discuss the calculation of the fast-ion pressure and energy density and perform a comparison with NUBEAM. The energy density is given by the following distribution function moment:

$$\begin{align*} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:energydens} \epsilon &= \frac{1}{2}m_{\rm fi} \int\!\!\!\int v^2 f \cdot 2\pi v^2 {\rm d} v {\rm d}\xi \\&= \frac{m_\text{fi}}{2}\cdot\frac{S \tau_\text{s}}{6} \left(3v^2 + v_\text{c}^2 \log \frac{(v+v_\text{c})^2}{v^2 -v v_\text{c}+v_\text{c}^2} - 2\sqrt{3}v_\text{c}^2 \arctan\frac{2v-v_\text{c}}{\sqrt{3}v_\text{c}}\right) + \mathit{const} \end{align*} \tag{ 60 }$$

and the pressure is given by $\newcommand{\e}{{\rm e}} p =\frac{2}{3} \epsilon$. The integral is solvable analytically. Figure 27 shows a comparison of this calculation (red dashed lines) with NUBEAM (black). The pressure profile shapes agree in general quite well, but in terms of absolute values, there is a systematic underestimation of about 5% (as seen in the volume-integrated fast-ion energy). This under-estimation arises from the neglection of the speed-diffusion: it creates a high-energy tail above the injection energy, which is increasingly important for higher velocity momenta. Below the injection energy, the effect of the speed diffusion is usually considered as negligible [32]. We will therefore focus to model only the high energy tail above the injection energy.

For this purpose, we can ignore the ξ-component (pitch-angle scattering) and consider only drag A and diffusion D with respect to $v$. To make the equation analytically better treatable, we consider constant drag and diffusion rates A and D, evaluated at the injection energy $v=v_0$. This is justified by the fact that the width of the high energy tail is supposed to be small compared to the birth velocity $v_0$ and A and D hence vary slowly. Then, the Fokker–Planck equation reads as:

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\del}{\partial} \displaystyle A \frac {\del f}{\del v} + D \frac {\del^2 f}{\del v^2} = \frac {\del f}{\del t}~{{\rm with:}} \nonumber \end{align} \tag{ 61 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle A = \frac{v_0^3+v_{\rm c}^3}{\tau_{\rm s} v_0^2} \quad D=\frac{T_{\rm e} v_0^3+ T_{\rm i}v_{\rm c}^3} {\tau_{\rm s} m_{\rm fi} v_0^3}. \nonumber \end{align} \tag{ 62 }$$

In the current version of RABBIT, we neglect the time-dependence of the speed-diffusion correction, and consider only the steady-state solution. This is justified, because the high-energy-tail reaches its steady-state on a faster time-scale than the slowing-down time. This can be infered e.g. from a NUBEAM calculation (figure 25) where we compare the build-up after the NBI-onset of the high-energy-tail (orange curve) with the build-up of the fast-ion density below the injection energy. We find the solution:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:speeddiffsol} f(v) = \frac{1 }{2 \pi} \frac{S \tau_{\rm s} }{v_0^3+v_{\rm c}^3} \exp{\frac{-(v-v_0)}{v_{\rm w}}} \nonumber \end{align} \tag{ 63 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle v_{\rm w} = D/A = \frac{T_{\rm e} v_0^3+ T_{\rm i}v_{\rm c}^3}{v_0^3+v_{\rm c}^3} \frac{1}{m_{\rm fi} v_0} \nonumber \end{align} \tag{ 64 }$$

which has also been obtained by [31, 32].

Zoom In Zoom Out Reset image size

Equivalently, one can perform a coordinate transformation from speed $v$ to energy E ($\newcommand{\del}{\partial} \frac{\del}{\del v} = mv \frac{\del}{\del E}$). The differential equation reads then:

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\del}{\partial} \displaystyle A_E \frac {\del f}{\del E} + D_E \frac {\del^2 f}{\del E^2} = \frac {\del f}{\del t} ~{{\rm with:}} \nonumber \end{align} \tag{ 65 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle A_E = m_{\rm fi} \frac{E_0^{3/2}+E_{\rm c}^{3/2}}{\tau_{\rm s} \sqrt{2 E_0 / m_{\rm fi}}} \quad D_E = m_{\rm fi} \frac{T_{\rm e} E_0^{3/2}+ T_{\rm i}E_{\rm c}^{3/2}} {\tau_{\rm s} \sqrt{2 E_0 / m_{\rm fi}} } \nonumber \end{align} \tag{ 66 }$$

which is solved by:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:energydiffsol} f(v) &= \frac{1 }{2 \pi} \frac{S \tau_{\rm s} }{v_0^3+v_{\rm c}^3} \exp{\frac{-(v^2-v^2_0)}{v_{\rm ef\,\!f}^2}} \nonumber \\ & \equiv \frac{1 }{2 \pi} \frac{S \tau_{\rm s} }{E_0^{3/2}+E_{\rm c}^{3/2}} \exp{\frac{-(E-E_0)}{T_{\rm ef\,\!f}}} \nonumber \end{align} \tag{ 67 }$$

$$\begin{align*} \newcommand{\e}{{\rm e}} \displaystyle v_{\rm ef\,\!f}^2 = \frac{D_E/A_E}{ \frac{1}{2}m_\text{fi}}= \frac{D/A\cdot m_{\rm fi}v_0}{\frac{1}{2}m_{\rm fi}} = \frac{2}{m_{\rm fi}} \frac{T_{\rm e} v_0^3+ T_{\rm i}v_{\rm c}^3}{v_0^3+v_{\rm c}^3} \equiv \frac{2 T_{\rm ef\,\!f}}{m_{\rm fi}}. \nonumber \end{align*} \tag{ 68 }$$

This approach is also found e.g. in [44].

In figure 26, we compare both solutions with a NUBEAM calculation, for a ASDEX Upgrade discharge with $T_{\rm e}=3.2$ keV and a particularly high bulk ion temperature of $T_{\rm i}=6.9$ keV, which leads to a relatively strong high-energy-tail due to speed diffusion. It can be seen that equation (67) gives a very good estimate of the slope of the tail, while the agreement between NUBEAM and equation (63) is slightly inferior, which is why we have implemented the former solution.

Zoom In Zoom Out Reset image size

The integral for the contribution of the high-energy-tail to the fast-ion energy density can be solved analytically (owing to the assumption of constant drag and diffusion rates):

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \frac{m_{\rm fi}}{2}\int_{v_0}^{+\infty} v^2 f \cdot 2\pi v^2 {\rm d} v \nonumber \end{align} \tag{ 69 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\Teff}{{v_{\rm ef\,\!f}^2}} \newcommand{\Te}{T_{\rm e}} \displaystyle =& \left(\frac{3}{8} \sqrt{\pi} v_{\rm ef\,\!f}^{5} \exp\frac{v_0^2}{\Teff} \left(1-{\rm erf}\, \frac{v_0}{v_{\rm ef\,\!f}}\right) + \frac{3}{4} v_{\rm ef\,\!f}^4 v_0 + \frac{1}{2} \Teff v_0^3 \right) \nonumber \\ & \times \frac{S \tau_{\rm s}} {v_0^3+v_{\rm c}^3} \frac{m_{\rm fi}}{2} \nonumber \end{align} \tag{ 70 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\Teff}{{v_{\rm ef\,\!f}^2}} \newcommand{\Te}{T_{\rm e}} \displaystyle \approx \left(\frac{3}{8} \frac{v_{\rm ef\,\!f}^{6}}{v_0} + \frac{3}{4} v_{\rm ef\,\!f}^4 v_0 + \frac{1}{2} \Teff v_0^3 \right) \frac{S \tau_{\rm s}} {v_0^3+v_{\rm c}^3} \frac{m_{\rm fi}}{2}. \nonumber \end{align} \tag{ 71 }$$

In the last step, we have used the first term of the Taylor approximation of $(1-{\rm erf}\, x)\approx e^{-x^2} (\frac{1}{\sqrt{\pi} x} - \frac{1}{2\sqrt{\pi} x^3} + \frac{3}{4\sqrt{\pi} x^5} + \dots)$ for large x.

For the correction to the fast-ion density we get:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle & \int_{v_0}^{+\infty} f \cdot 2\pi v^2 {\rm d} v \nonumber \\ & = \frac{v_{\rm ef\,\!f}^2}{4} \left(2 v_0 + \sqrt{\pi} v_{\rm ef\,\!f} \exp\frac{v_0^2}{v_{\rm ef\,\!f}^2}\left(1-{\rm erf}\,\frac{v_0}{v_{\rm ef\,\!f}}\right)\right) \frac{S \tau_{\rm s}} {v_0^3+v_{\rm c}^3} \nonumber \end{align} \tag{ 72 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \approx \frac{v_{\rm ef\,\!f}^2}{4} \left(2 v_0 + \frac{v_{\rm ef\,\!f}^2}{v_0} - \frac{v_{\rm ef\,\!f}^4}{2 v_0^3}+\frac{3 v_{\rm ef\,\!f}^6}{4 v_0^5} \right) \frac{S \tau_{\rm s}} {v_0^3+v_{\rm c}^3} \nonumber \end{align} \tag{ 73 }$$

where it is more beneficial to use the first three terms from the Taylor expansion. In both cases, the considered terms of the expansion are good enough to yield only an error in the order of 10⁻³ for a relatively low $v_0/v_{\rm ef\,\!f}=2$—we consequently use the Taylor approximation for higher $v_0/v_{\rm ef\,\!f}$ and the exact solution for lower $v_0/v_{\rm ef\,\!f}$.

Figure 27 shows in green lines the result of this speed-diffusion correction, which is a significant improvement.

Zoom In Zoom Out Reset image size

The RABBIT code consists of three parts: in the first part, the beam attenuation and fast-ion birth profile is calculated. We use a realistic collisional radiative model resolved for atomic states, which is evaluated in a simplified beam geometry (only along the center-line of the beam). The effect of finite beam-width is included analytically. In the second part, the orbit average of the source term is calculated. This is done with a conventional Runge–Kutta integrator using arbitrary numerical axisymmetric equilibria. To speed up the calculation, only a very low number of orbits are calculated (≈20 per beam energy component) and a suitable interpolation is used in between. The third part calculates the time-dependent solution of the Fokker–Planck equation. Since it is almost purely analytically calculated, this part is clearly the fastest part, while the other two parts take approximately 12 ms each. In total, the RABBIT code needs  ≈25 ms per time-step and beam (serial), which is roughly a factor of 1000 faster than NUBEAM. The calculation time could be further reduced, e.g. by faster hardware, further code optimizations or higher parallelization. The reported 25 ms have been estimated on an Intel Xeon E5-2680 v3 CPU with 2.5 GHz clockrate, using a parallelizing scheme with one thread per beam. A higher parallelization could, e.g. be achieved trivially by using one thread per beam energy component.

A comprehensive benchmark is carried out with the more accurate NUBEAM code, indicating very good agreements for heating profiles as well as profiles of fast-ion density and pressure. The fast-ion current density (necessary for neutral-beam current drive calculations) is more challenging to calculate. In the considered test-cases, we see typical deviations in the order of  ±10% between RABBIT and NUBEAM which should still be acceptable for many applications. The fast-ion current depends directly on $v_{\parallel}$ and thus on the orbital motion during the whole slowing down process. In particular, orbit topology changes (e.g. due to pitch angle scattering) are supposed to play an important role, because initially passing ions carry no longer a significant current when they have become trapped. In a comparison with a Monte-Carlo-based orbit-average (based on a fully realistic 3D beam deposition calculation by FIDASIM) we could check that these deviations mainly stem from the concept of orbit-averaging, and not from the further approximations that RABBIT introduces (namely simplified beam geometry and the orbit interpolation technique). Finding an ad hoc correction for this effect remains a challenging task for the future, but could help to further improve the agreement with NUBEAM.

Several applications of the code are already in planning. An interface with the discharge control system of ASDEX Upgrade and with the RAPTOR code has already been started, and it should allow an improved real-time control of the plasma discharges. Here, the kinetic input profiles will be provided by estimates from RAPTOR, which are constrained by the real-time diagnostics (when available) and complemented by model-based estimates. Currently, real-time T_e and n_e diagnostics are used, and there are plans to make the charge-exchange diagnostics (T_i and rotation) real-time capable as well. The plasma equilibrium is reconstructed in real-time by the JANET code [45]. An extension of these activities and further applications to other tokamaks such as JET, DIII-D and TCV is planned.

It should be noted that RABBIT is also well suited for non-real-time applications, where fast models are desired. For example, it can be used for fast intershot analyzes which are needed to facilitate the decision-making for the next discharge. To this end, RABBIT is foreseen to be used for accurate equilibrium reconstructions at ASDEX Upgrade with the IDE code [46, 47], which allow sophisticated estimates of e.g. the safety factor (q) profile directly after a discharge. This will aid e.g. the development of advanced scenarios, where a fine-tuning of the q-profile is desired. In addition, RABBIT could also be used for scenario optimization or reactor design studies, where the short calculation time will allow large parameter scans.

As outlook, several extensions of the code are possible in the future. The calculation of neutron rates has already been tested in initial benchmarks on ASDEX Upgrade and DIII-D. The inclusion of charge-exchange losses could be relevant for smaller machines such as TCV. In addition, radial transport could be included within our time-dependence numerical scheme, which could allow to model radial transport, e.g. due to MHD activity or plasma turbulence.

We like to thank Mirko Salewski for careful reading of the manuscript and valuable comments. This work has been carried out within the framework of the EUROfusion Consortium and has received funding from the Euratom research and training programme 2014–2018 under grant agreement No 633053. The views and opinions expressed herein do not necessarily reflect those of the European Commission.