Analytical and Statistical Modelling of a Fast Ion Source Formed by Injection of a Neutral Beam into Magnetically Confined Plasma

Abstract

Mathematical modelling of heating and current drive as well as yields and distributions of fusion products in a magnetically confined plasma subject to neutral beam injection requires, in turn, modelling of distributions of fast ions, which is a complex task including calculations of the source of suprathermal particles, i.e., the number of fast ions occurring in unit volume during unit time owing to the injection of fast atoms. The knowledge of the magnetohydrodynamic equilibrium, beam injection geometry and spatial distribution of the magnetic field are the necessary prerequisites. Explicit general analytical formulae for the source of fast ions have been obtained by two different methods. In addition, a method of statistical modelling is presented. Calculations of spatial and angular distributions of the fast ion source for a tokamak and verifications of the obtained results have been performed by a number of methods.

1. Introduction

At an early stage of controlled nuclear fusion research in review paper [1] L.A. Artsimovich mentioned the injection of beams of fast neutral particles into plasma as one of the natural ways of achieving the required level of ion temperatures and predicted that, although little or nothing had then been done to develop this method, probably it would own the future, taking an important place in the overall research programme. In [2], T. Ohkawa suggested the use of beams of neutral atoms as a method of driving the plasma current.

Today neutral beam injection (NBI) is an established method of heating and non-inductive generation of an electric current in many experimental fusion devices with magnetic confinement of plasma [3]. ITER will be equipped with the most powerful neutral beam injectors ever, capable of delivering higher energy beams for longer than any previous systems [4]. Prototypes of neutral beam injectors [5,6] are being developed for the Tokamak with Reactor Technologies (TRT) [7]. NBI is also being considered for demonstrative fusion power plant concepts (DEMO) [8].

In the works on the development of fusion neutron sources and fusion-fission hybrid systems NBI is considered as the main source of high-energy particles generating primary neutrons [9,10]. Calculations of spatial, angular, and energetic distributions of the neutron source in a plasma and the subsequent modelling of neutron transport processes first require distributions of velocities of plasma ions to be obtained, which in turn require calculations of the penetration of fast neutral beams into plasma.

Calculations of the penetration and deposition of neutral beams for large tokamaks such as TFTR [11,12], EAST [13,14], KSTAR [15], and JET [16] were made on the basis of [17,18]. Simulations of the neutral beam injection in ITER were reported in [19]. NBI-induced effects in plasma of a compact fusion neutron source were studied in [20]. Applications of NBI modelling for stellarator/heliotron devices, such as LHD [21] and W7-X [22], were based on the works [23] and [24] correspondingly.

The above-mentioned bibliographic sources are mostly relying on basic descriptions of different sorts of simulation techniques predominantly related to various Monte Carlo methods, rather than explicit general analytical closed-form expressions for the source function of fast ions produced in a magnetically confined plasma due to neutral beam injection. The treatments in [25,26] include a larger mathematical background. A compact overview of commonly adopted approaches to describe the penetration of neutral beams into fusion plasmas is given in [27]. Among newer works dedicated to the modelling of neutral beam injection, articles [28,29,30] can be referenced. As stated in [25,27,31], if the confinement of fast ions is sufficiently good, an ion originating from NBI can be assumed to stay on its magnetic surface. The same assumption is adopted herein. Deviations from such a basic model, i.e., corrections for drift orbits, are beyond the scope of this study and may be introduced on top of the present results.

The purpose of this article is to present rigorous general derivations of the source function by two different methods, providing a cross check, in addition to statistical modelling. The approach presented here is based on explicit analytical formulae, which are reproducible and exact in the frame of the adopted physics model. This approach is preferable from the viewpoint of the development of the physics basis as well as the reliability and the computational speed of the integrated systems codes where the complexity necessitates optimisation.

Two alternative derivations of the analytical formulae are given in Section 2. Statistical modelling based on kernel smoothing is described in Section 3. Examples of calculations of spatial and angular distributions of the source of fast ions originating from neutral beam injection into a tokamak are shown in Section 4, followed by conclusions in Section 5.

2. Analytical Approach

Two methods are given in this section to derive an analytical formula for the local source of suprathermal ions supplied by a monoenergetic beam of fast neutral atoms injected into a magnetically confined fusion plasma. The local source $S$ [cm⁻³s⁻¹] is the number of ions, occurring in unit volume during unit time, expressed as a function of coordinates. A phenomenological derivation method and a more formal mathematical method lead to the same result, thus providing a validating cross check.

2.1. Phenomenological Derivation

First, let us consider a mathematical line instead of a finite width beam. Figure 1 shows such a line and a cross-section of nested magnetic surfaces in a device for the magnetic confinement of plasma. It is assumed that once an injected H⁰ (or D⁰, or T⁰) atom loses its electron and becomes an ion, it remains at the same magnetic surface, i.e., that the source $S(\rho )$ is a function of the dimensionless magnetic flux surface label $\rho$, which is a generalisation of radius. This can be e.g., the normalised poloidal radius defined as

$$\rho =\sqrt{\frac{\mathsf{\Psi }-{\mathsf{\Psi }}_a}{{\mathsf{\Psi }}_b-{\mathsf{\Psi }}_a}},$$

(1)

where $\mathsf{\Psi }$ is the local value of the poloidal magnetic flux function, ${\mathsf{\Psi }}_a$ is its value at the magnetic axis, where $\rho =0$, and ${\mathsf{\Psi }}_b$ is its value at the boundary, where $\rho =1$. Let us denote the length along the injection line as $\mathsf{\Lambda }$. The line is crossing the last closed flux surface $\rho =1$ at $\mathsf{\Lambda }={\mathsf{\Lambda }}_{en}$, where it is entering the plasma and at $\mathsf{\Lambda }={\mathsf{\Lambda }}_{ex}$, where it is exiting the plasma. Let us denote the neutral particle current carried by this injection line as $I$ [s⁻¹].

The probability that an atom entering the plasma at ${\mathsf{\Lambda }}_{en}$ penetrates to a position $\mathsf{\Lambda }$ is given by the Poisson exponent

$$\mathrm{P}(\mathsf{\Lambda })=e^{-{\displaystyle \underset{{\mathsf{\Lambda }}_{en}}{\overset{\mathsf{\Lambda }}{\int }}\frac{dl}{{\lambda }_{mfp}}}},$$

(2)

where ${\lambda }_{mfp}$ is the mean free path of an atom with respect to electron loss, so the integral ${\displaystyle \underset{{\mathsf{\Lambda }}_{en}}{\overset{\mathsf{\Lambda }}{\int }}\frac{dl}{{\lambda }_{mfp}}}$ is the average number of electron loss events along the path from ${\mathsf{\Lambda }}_{en}$ to $\mathsf{\Lambda }$. This means that the injected neutral particle density along the line $n_b(\mathsf{\Lambda })$ decreases exponentially as

$$n_b(\mathsf{\Lambda })=n_b({\mathsf{\Lambda }}_{en})e^{-{\displaystyle \underset{{\mathsf{\Lambda }}_{en}}{\overset{\mathsf{\Lambda }}{\int }}\frac{dl}{{\lambda }_{mfp}}}},$$

(3)

where $n_b({\mathsf{\Lambda }}_{en})$ is the initial density at the entrance to the plasma. The electron loss frequency $\nu$ [s⁻¹] is

$$\nu =\frac{{\mathrm{v}}_b}{{\lambda }_{mfp}},$$

(4)

where ${\mathrm{v}}_b$ is the injection velocity. By electron loss throughout this article all the electric charge changing processes of the type $H^0\to H^{+}$ are meant, including charge exchange, impact ionisation, etc.

The above considerations lead to the following qualitatively derived source function of fast ions in the plasma at the position $\mathsf{\Lambda }$ along the fast atom injection line

$$S(\mathsf{\Lambda })=\nu n_b(\mathsf{\Lambda })=n_b({\mathsf{\Lambda }}_{en})\frac{{\mathrm{v}}_b}{{\lambda }_{mfp}}{e}^{-{\displaystyle \underset{{\mathsf{\Lambda }}_{en}}{\overset{\mathsf{\Lambda }}{\int }}\frac{dl}{{\lambda }_{mfp}}}}.$$

(5)

However, this qualitative derivation above is not yet properly taking into account the geometry. The desired quantity is the source of fast ions at a given magnetic surface. As Figure 2 illustrates, a given increment of the effective radius $d\rho$ may correspond to different increments $d\mathsf{\Lambda }$ of the length along the injection line, depending on the injection direction.

Let us attribute a cross-section area $dS$ to the injection line, considering it as a narrow infinitesimal beam, as shown in Figure 3. Then $dSd\mathsf{\Lambda }$ is the volume of the cylinder between the magnetic surfaces $\rho$ and $\rho +d\rho$.

The number of atoms undergoing electron loss collisions in the depicted volume $dSd\mathsf{\Lambda }$ at the position $\mathsf{\Lambda }$ per second is the product of this volume $dSd\mathsf{\Lambda }$ and the number of electron loss events per unit volume during unit time $S(\mathsf{\Lambda })$, given by (5). Assuming that all the resulting fast ions are staying at the corresponding magnetic surface $\rho$, we need to divide the mentioned quantity, i.e., the number of atoms undergoing electron loss collisions per second, by the corresponding volume $dV$ enclosed between the magnetic surfaces $\rho$ and $\rho +d\rho$. This will be the desired source function $S(\rho )$ per unit volume during unit time. Using (5) we have

$$S(\rho )=I\frac{e^{-{\displaystyle \underset{{\mathsf{\Lambda }}_{en}}{\overset{\mathsf{\Lambda }}{\int }}\frac{dl}{{\lambda }_{mfp}\left(l\right)}}}}{{\lambda }_{mfp}\left(\mathsf{\Lambda }\right)}\frac{\left(\frac{d\mathsf{\Lambda }}{d\rho }\right)}{\left(\frac{dV}{d\rho }\right)},$$

(6)

because the initial injection current at the entrance equals

$$I=n_b\left({\mathsf{\Lambda }}_{en}\right){\mathrm{v}}_{b}dS.$$

(7)

Such is the phenomenological derivation of the analytic formula for the fast ion source function. Since $\mathsf{\Lambda }(\rho )$ is a multivalued function, the final formula needs to contain the sum over all branches of this function. An example is depicted in Figure 4a showing an injection line crossing the magnetic surfaces, and in Figure 4b showing the variation of the effective radius versus the distance along this injection line. The value ${\rho }^{*}$ designates some particular magnetic surface. There are four branches with alternating monotonicity in this example. These branches are shown by different colours. This illustrates that different values of $\mathsf{\Lambda }$ may contribute to the source function $S(\rho )$ at the same magnetic surface $\rho$.

Therefore, the formula becomes

$$S(\rho )={\displaystyle \sum _k\frac{I}{{\lambda }_{mfp}\left({\mathsf{\Lambda }}_k\right)}}{e}^{-{\displaystyle \underset{{\mathsf{\Lambda }}_{en}}{\overset{{\mathsf{\Lambda }}_k}{\int }}\frac{dl}{{\lambda }_{mfp}\left(l\right)}}}\frac{\left(\frac{d\mathsf{\Lambda }}{d\rho }\right)\left({\mathsf{\Lambda }}_k\right)}{\left(\frac{dV}{d\rho }\right)\left({\mathsf{\Lambda }}_k\right)}.$$

(8)

The argument ${\mathsf{\Lambda }}_k$ is the spatial location. This is the value of $\mathsf{\Lambda }$ corresponding to the given $\rho$ at the k-th branch. The value $\frac{dV}{d\rho }$ does not depend on the branch index k. The mean free path of a fast hydrogen atom in a fusion plasma in accordance with the models [32,33] is

$${\lambda }_{mfp}=\frac{1}{n_e{\sigma }_s},$$

(9)

where $n_e$ is the electron density and ${\sigma }_s$ is the effective stopping cross-section depending on the electron temperature $T_e$, the electron density $n_e$, and impurity densities. The value ${\sigma }_s$ is determined by parameterisations given in [32,33]. The mean free path ${\lambda }_{mfp}$ given by (9) is also a function of the magnetic surface, i.e., does not depend on the branch index k in (8). Taking ${\lambda }_{mfp}$ and $\frac{dV}{d\rho }$ out of the summation over the branches we obtain

$$S(\rho )=\frac{I}{{\lambda }_{mfp}\left(\rho \right)\frac{dV}{d\rho }\left(\rho \right)}{\displaystyle \sum _k{\left(\frac{d\mathsf{\Lambda }}{d\rho }\right)}_k{e}^{-{\displaystyle \underset{{\mathsf{\Lambda }}_{en}}{\overset{{\mathsf{\Lambda }}_k}{\int }}\frac{dl}{{\lambda }_{mfp}\left(l\right)}}}}.$$

(10)

Such is the source function of fast ions in the plasma due to the injection of fast atoms along a mathematical line representing a narrow beam.

2.2. Mathematical Derivation

A formula identical to (10) can be obtained formally using the language of mathematical statistics as follows. The probability for an atom to undergo no electron loss collisions along the path inside the plasma from ${\mathsf{\Lambda }}_{en}$ to the point $\mathsf{\Lambda }$ on the injection line is given by (2). The value $\mathsf{\Lambda }$, where the electron loss takes place, is a random number. The cumulative distribution function for $\mathsf{\Lambda }$ is then

$$\mathrm{F}(\mathsf{\Lambda })=1-\mathrm{P}(\mathsf{\Lambda }),$$

(11)

which is the probability that the electron loss does occur somewhere along the path inside the plasma from ${\mathsf{\Lambda }}_{en}$ to the point $\mathsf{\Lambda }$, i.e., that the atom does not penetrate as deep as to the point $\mathsf{\Lambda }$. For the probability density function (p.d.f.) for $\mathsf{\Lambda }$, using (2), we therefore have

$$f(\mathsf{\Lambda })=\frac{d\mathrm{F}}{d\mathsf{\Lambda }}=\frac{1}{{\lambda }_{mfp}}{e}^{-{\displaystyle \underset{{\mathsf{\Lambda }}_{en}}{\overset{\mathsf{\Lambda }}{\int }}\frac{dl}{{\lambda }_{mfp}\left(l\right)}}}.$$

(12)

The probability that the electron loss occurs between $\mathsf{\Lambda }$ and $\mathsf{\Lambda }+d\mathsf{\Lambda }$ is

$$d\mathrm{P}=f(\mathsf{\Lambda })d\mathsf{\Lambda }.$$

(13)

Let $d𝒫$ be the probability that the electron loss occurs between $\rho$ and $\rho +d\rho$. To calculate $d𝒫$ we need the summation over all ${\mathsf{\Lambda }}_k$ values corresponding to the given $\rho$ value, i.e., over all branches of the function $\mathsf{\Lambda }(\rho )$, thus

$$d𝒫={\displaystyle \sum _{k}f\left({\mathsf{\Lambda }}_k\right)\frac{d{\mathsf{\Lambda }}_k}{d\rho }d\rho }=\frac{1}{{\lambda }_{mfp}\left(\rho \right)}\left({\displaystyle \sum _k\frac{d{\mathsf{\Lambda }}_k}{d\rho }{e}^{-{\displaystyle \underset{{\mathsf{\Lambda }}_{en}}{\overset{{\mathsf{\Lambda }}_k}{\int }}\frac{dl}{{\lambda }_{mfp}\left(l\right)}}}}\right)d\rho .$$

(14)

This means that the probability density $f\left(\rho \right)$ for the $\rho$ coordinate of the electron loss location is

$$f\left(\rho \right)=\frac{1}{{\lambda }_{mfp}\left(\rho \right)}{\displaystyle \sum _k{\left(\frac{d\mathsf{\Lambda }}{d\rho }\right)}_k{e}^{-{\displaystyle \underset{{\mathsf{\Lambda }}_{en}}{\overset{{\mathsf{\Lambda }}_k}{\int }}\frac{dl}{{\lambda }_{mfp}\left(l\right)}}}}.$$

(15)

The source function of fast ions at the magnetic surface $\rho$ is

$$S\left(\rho \right)=I\frac{d𝒫\left(\rho \right)}{dV\left(\rho \right)},$$

(16)

where $I$ [s⁻¹] is the ‘input’ neutral current attributed to the injection line under consideration, $dV\left(\rho \right)$ is the volume where the ions reside, occurring due to the electron loss by the injected atoms, $d𝒫\left(\rho \right)$ is the probability that the electron loss occurs particularly at the given $\rho$ value. Using (14), (15) and (16) we arrive at the source function

$$S\left(\rho \right)=\frac{I}{dV/d\rho }f\left(\rho \right)=\frac{I}{{\lambda }_{mfp}\left(\rho \right)\left(\frac{dV}{d\rho }\right)\left(\rho \right)}{\displaystyle \sum _k{\left(\frac{d\mathsf{\Lambda }}{d\rho }\right)}_k{e}^{-{\displaystyle \underset{{\mathsf{\Lambda }}_{en}}{\overset{{\mathsf{\Lambda }}_k}{\int }}\frac{dl}{{\lambda }_{mfp}\left(l\right)}}}},$$

(17)

which coincides with (10), i.e., the formal result is the same as the phenomenologically derived expression.

3. Statistical Approach

An alternative way to calculate $f\left(\rho \right)$ instead of Formula (15) is to use the so-called kernel probability density function (p.d.f.) estimator

$$f\left(\rho \right)\approx \frac{1}{Nh}{\displaystyle \sum _{i=1}^{N}K\left(\frac{\rho -{\rho }_i}{h}\right)},$$

(18)

where $K$ is called a kernel function, $h$ is an arbitrary positive number called the kernel bandwidth, and $N$ is the number of random values of the electron loss locations ${\rho }_i$. The selection of an optimum bandwidth and an optimum kernel function is a separate subject in mathematical statistics. With $I_{(-1{,}_{}1)}\left(z\right)$ being equal to unity inside the $\left(-1{,}_{}1\right)$ interval and equal to nought outside, kernel functions typically in use include the rectangular window function

$$K(z)=\frac{1}{2}{I}_{(-1{,}_{}1)}\left(z\right),$$

(19)

Gaussian kernel

$$K(z)=\frac{1}{\sqrt{2\pi }}{e}^{-z^2/2},$$

(20)

Epanechnikov kernel

$$K(z)=\frac{3}{4}\left(1-z^2\right)I_{(-1{,}_{}1)}\left(z\right),$$

(21)

and others.

The algorithm to calculate $f\left(\rho \right)$ via the kernel p.d.f. estimator (18) is as follows. First of all, $N$ uniformly distributed random numbers $u_i$ in the $\left(0{,}_{}1\right)$ interval need to be generated, where $i\in \overline{1,N}$. Each $u_i$ will represent the value of the probability (2) for an atom to reach the location ${\mathsf{\Lambda }}_i$ along the injection line in plasma in accordance with

$$u_i=e^{-{\displaystyle \underset{{\mathsf{\Lambda }}_{en}}{\overset{{\mathsf{\Lambda }}_i}{\int }}\frac{dl}{{\lambda }_{mfp}\left(l\right)}}}.$$

(22)

Next, the corresponding ${\mathsf{\Lambda }}_i$ values need to be calculated from (22) for each of the $N$ numbers $u_i$. The knowledge of the local mean free path ${\lambda }_{mfp}$ of an atom in the plasma is required, as a function of the kinetic energy of the atom and the plasma parameters depending on the spatial coordinates. After that, the obtained values ${\mathsf{\Lambda }}_i$ of the length along the injection line need to be recalculated into the corresponding values of the effective radius ${\rho }_i$ where the electron loss occurs. For this, the equation of the injection line is needed and the transformation of the coordinates in configuration space (Cartesian of cylindrical) to flux coordinates using the magnetohydrodynamic equilibrium data determining the geometry of magnetic surfaces. Finally, the obtained ${\rho }_i$ values need to be substituted to (18).

Calculations show that both the statistical approach (18) and the analytical Formula (15) for the probability density $f\left(\rho \right)$ for the $\rho$ coordinate of the electron loss location lead to the same source function $S(\rho )$.

4. Calculations of Spatial and Angular Distributions of the Fast Ion Source

The initial angle between the fast ion velocity and the magnetic field, i.e., the angle between the injection line and the magnetic field at the spatial point where the fast neutral beam atom undergoes the electron loss, is obviously not constant along the injection line. The source of fast ions has a certain angular distribution. This section illustrates the approach to calculate this distribution.

In addition, the source function $S(\rho )$ considered in Section 2 was obtained for the injection along a mathematical line, which is an abstraction, necessary as a building block for the modelling of realistic beams. This section also describes the approach to calculate the source of fast ions in the plasma for the case of a finite width of a neutral beam.

Figure 5a shows the poloidal cross-section of magnetic surfaces, i.e., the nested isolines of the normalised poloidal radius calculated on the basis of magnetohydrodynamic equilibrium data for the DEMO-FNS project [34] for a candidate operating regime used herein as a sample case. Either analytical equilibria such as the Soloviev solution [35] of the Grad–Shafranov equation or numerical magnetohydrodynamic codes such as [36] can be employed. The values $R_{ax}$ and $Z_{ax}$ in Figure 5a refer to the magnetic axis position, $a_h$ is the horizontal minor radius, the values $\epsilon (1)$, $\tau (1)$, and $\Delta (1)$ correspond to the elongation, triangularity and the Shfranov shift calculated for the last closed flux surface (LCFS), $V(1)$ denotes the volume enclosed by the LCFS, and $S(1)$ is the LCFS area. The functions of the volume enclosed by magnetic surfaces and its derivative, with respect to the normalised poloidal radius, are shown in Figure 5b versus the normalised poloidal radius defined by (1).

Details of the implementation of the analytical geometry to calculate coordinates along injection lines and the transformation between cylindrical coordinates $(R,Z,\phi )$ and flux coordinates $(\rho ,\vartheta ,\phi )$ are not discussed here since these both are rather common techniques. Suffice it to mention that the volume enclosed by the magnetic surface labeled $\rho$ is given by

$$V(\rho )={\displaystyle \underset{0}{\overset{\rho }{\int }}{\displaystyle \underset{0}{\overset{2\pi }{\int }}{\displaystyle \underset{0}{\overset{2\pi }{\int }}R(\tilde{\rho },\vartheta ,\phi )\left|J\right|d\tilde{\rho }d\vartheta d\phi }}},$$

(23)

the azimuthal (toroidal) angle $\phi \in \left[0,2\pi \right)$ preserves when azimuthally symmetric magnetic surfaces are assumed for tokamaks, and the Jacobian determinant is

$$J=\left|\begin{array}{ccc}\frac{\partial R}{\partial \rho }& \frac{\partial R}{\partial \vartheta }& 0\\ \frac{\partial Z}{\partial \rho }& \frac{\partial Z}{\partial \vartheta }& 0\\ 0& 0& 1\end{array}\right|.$$

(24)

Jacobian determinant (24) calculated using the DEMO-FNS magnetohydrodynamic equilibrium data is depicted in Figure 6 as a surface plot versus the normalised poloidal radius $\rho$ and the poloidal angle $\vartheta$ to illustrate a typical example of the qualitative behaviour for a tokamak equilibrium where the volume derivative is

$$\frac{dV}{d\rho }=2\pi {\displaystyle \underset{0}{\overset{2\pi }{\int }}R(\rho ,\vartheta ,\phi )\left|J\right|d\vartheta }.$$

(25)

The coordinate transformation procedure is involved in such calculations.

A 3D mesh in Figure 7a illustrates the last closed flux surface corresponding to the normalised poloidal radius $\rho =1$. The three lines in Figure 7a show three different injection directions corresponding to three different values of tangency radius. The values of the effective radius (dashed lines) and the angle (solid lines) between the injection line and the magnetic field vector are shown in Figure 7b for each of the injection lines by the corresponding colours. The radial profiles of the plasma parameters used in calculations are shown in Figure 8. The profiles were fixed, and no self-consistent iterative procedures were involved herein to simulate the effect of neutral beam injection on the plasma parameters. The specified impurity densities represent a sample model example used for illustrative purposes. In practice either numerical modelling, or experimental data are required.

Actual neutral beams of a finite width are approximated by a large number of mathematical lines, i.e., infinitesimal beams. Injection geometry and beam divergence are taken into account. A certain neutral current is attributed to each of the injection lines in accordance with the beam current density distribution in the beam cross-section plane. Figure 9a shows such approximations of beams with cross-section areas comparable with the poloidal cross-section area of the plasma column. Three beams are shown, each of which is approximated by a large number of mathematical lines. The central lines of the beams are the lines shown in Figure 7a by the corresponding colours.

The resultant source function of fast ions in the plasma is calculated for a neutral beam of a finite width as the sum of the source functions obtained as described in Section 2 for each of the mathematical lines all together approximating the beam. An example of such a calculation is shown in Figure 9b for three neutral beams shown in Figure 9a assuming the magnetohydrodynamic equilibrium depicted in Figure 5 and Figure 6.

The distribution of the source of fast ions in the plasma noticeably depends on the beam injection geometry. Different colours in Figure 9b correspond to different tangency radii of the beam central lines. Mathematical modelling of the interaction of fast ions with the background plasma is described in [37] and references therein.

The initial beam current density distribution in the beam cross-section plane is shown in Figure 10a, and the residual beam current density distribution where the “remainder” of the beam “exits” the plasma is shown in Figure 10b. This is the so-called shine-through which is rather small in this example being three orders of magnitude weaker than the initial beam current density. For a more peripheral injection, i.e., for the beam with tangency radius R = 2906 mm shown by the yellow colour in Figure 9a, the shine-through current density is greater reaching several per cent of the initial current density.

The angular distribution of the source function at each magnetic surface is calculated as a histogram of the values of the angle between the injection line and the magnetic field at the point where the beam atom undergoes the electron loss. Figure 11a,b show the obtained angular distributions at $\rho =0.4$ and at $\rho =0.8$ respectively. Angular distributions of the source of fast ions in the plasma noticeably depend on the beam injection geometry as well.

Figure 12 shows verifications of the spatial and angular distributions of the source of fast ions in the plasma due to neutral beam injection. The lowermost curves are those obtained for finite width beams with nominal cross-section areas. The other coloured curves were calculated for narrower beams with smaller and smaller cross-section areas. Figure 12a shows that as the beam cross-section is decreasing, the source function is approaching the radial distribution, shown by the grey colour, obtained for a mathematical line corresponding to the central line of the beam. Similarly, Figure 12b shows that as the beam cross-section is decreasing, the angular distribution is approaching the curve, shown by the grey colour, obtained for a mathematical line corresponding to the central line of the beam, as it ought to be.

It should be noted that although an azimuthally symmetric sample case is illustrated for simplicity, the presented general formalism for the fast ion source due to neutral beam injection is certainly applicable equally well for stellarator/heliotron devices, however the treatment of the equilibrium and coordinates will be different, e.g., based on [38].

5. Conclusions

Explicit general analytical formulae have been obtained for the source of fast ions originating from the injection of fast neutral beams into magnetically confined fusion plasma. Two derivation methods for the analytical formulae have been demonstrated. Both the phenomenological and the formal mathematical approach lead to identical results, which fact is one of successfully performed verifications of the obtained results. In addition, an alternative method has been described for calculations of the source of fast ions in a beam heated plasma. This method is based on mathematical statistics and the use of kernel estimation of the probability density, also called kernel smoothing. For this approach related to statistical modelling, an explicit mathematical formulation is given as well. This method of statistical modelling of the penetration of neutral beams into plasma produces the same results as calculations by the analytical formulae, which is an additional verification.

The obtained results have been applied to a sample case input data corresponding to a candidate operating regime of a classical tokamak DEMO-FNS [34]. Spatial and angular distributions of the source of fast ions in plasma subject to neutral beam injection have been obtained. As an extra verification, a number of test calculations have been performed with a gradual decrease in the cross-section area of the injected beam. As it ought to be, in this case the results approach the distributions calculated independently for the simple limiting case of a mathematical line.

The use of a rigorous analytical approach is beneficial for the development of the physics basis of controlled fusion and the reliability, reproducibility, and operation speed of integrated numerical codes.