The effect of density fluctuations on electron cyclotron beam broadening and implications for ITER

The wave kinetic equation with random fluctuations was originally formulated in [20] adapting a statistical approach proposed by Karal and Keller [21] in combination with Weyl-symbol calculus for the phase-space representation of the wave field [22]. The specific form of the wave kinetic equation relevant to coherent wave beams and implemented in WKBeam [17] is usually written in coordinates normalized to the scale L of the equilibrium plasma profiles so that the semiclassical parameter $\kappa = \omega L /c = k_0 L$ appears explicitly in the theory; here ω is the launched frequency, c the speed of light in free space, and $k_0 = \omega/c$ is the wave number in free space. It is important to note that the scale length L refers to the background plasma equilibrium which does not include turbulent fluctuations. Hence the parameter κ is large even in presence of short-scale fluctuations and the semiclassical limit $\kappa \to +\infty$ can be considered, yet not directly on the solution of Maxwell's equations for the wave field, but rather on the averaged wave energy density, which is represented by the averaged Wigner function. The wave kinetic equation describes the statistically averaged Wigner function of the beam in the limit $\kappa\to+\infty$. We refer to Weber [17] for a sketch of the derivation. Here instead we start from the wave kinetic equation written in physical coordinates as appropriate for the application. It reads

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \{H_{\alpha},w_{\alpha} \} = - 2k_0 \gamma_{\alpha} w_{\alpha} + \sum_{\beta} S_{\alpha\beta} \left(w_\alpha, w_\beta \right) \nonumber \end{align} \tag{ 1 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \; H_{\alpha}w_{\alpha} = 0, \nonumber \end{align} \tag{ 2 }$$

where, for each cold plasma propagation mode (namely the ordinary mode (O-mode, $\alpha = {\rm O}$) or the extra-ordinary mode (X-mode, $\alpha = {\rm X}$)) labeled by Greek indices, $H_\alpha$ is the Hamiltonian of standard ray-tracing theory [23, 24], $\gamma_{\alpha}$ the absorption coefficient, $S_{\alpha\beta}$ the scattering operator to be discussed below in more detail. The unknown is the statistically averaged Wigner function $w_{\alpha} = w_{\alpha}(x, N)$ which depends on the spatial coordinates x and the refractive index N. The left-hand side is written in terms of canonical Poisson brackets on the $(x, N)$ phase space, given by

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \{H_{\alpha},w_{\alpha} \} = \nabla_NH_{\alpha}\cdot \nabla_xw_{\alpha}-\nabla_xH_{\alpha}\cdot \nabla_Nw_{\alpha}. \nonumber \end{align} \tag{ 3 }$$

One should notice that the wave kinetic equation has the form of a steady-state kinetic equation with an additional constraint. The Wigner function can be associated to the electric field energy density. It has an analogous role to distribution function of the kinetic theory as we can derive integral expressions over the Wigner function for any quantity that can be written in terms of $\vert E\vert ^2$, where E is the electrid field. An example of such a quantity is, e.g. the absorption profile. The constraint merely means that there is no energy on phase-space points that do not satisfy the dispersion relation of the mode, $H_\alpha =0$. Further details on the Monte Carlo solution of the wave kinetic equation and on how to construct such integrals of Wigner function are given in [17].

In the WKBeam code, the Hamiltonians $H_\alpha = e^*_\alpha D e_\alpha$ are computed from the Hermitian part of the dispersion tensor $D_{ij} = N^2\delta_{ij} - N_iN_j - \varepsilon_{ij}$ with $\varepsilon = (\varepsilon_{ij})$ being the Hermitian part of the dielectric tensor [25] and $e_\alpha(x, N)$ the polarization unit vector of the mode α.

The absorption coefficient $\gamma_\alpha$, on the other hand, is computed from the hot-plasma dielectric tensor accounting for relativistic electron-cyclotron resonance interaction. More specifically, one can write $\gamma_\alpha = {\rm Im}(N) \cdot \nabla_NH_\alpha$ where ${\rm Im}(N)$ is the imaginary part of the refractive index vector which is computed in the WKBeam code by the same routines used in the complex-geometrical-optics code GRAY [26] and the paraxial WKB code TORBEAM [27].

The general form of the scattering operator includes cross-polarization scattering terms from mode α to mode β and is given by

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:scattop} S_{\alpha\beta} \left(w_\alpha, w_\beta \right) (x, N) =&\, \int \big[ \sigma_{\beta\alpha}(x,N^{\prime},N)w_\beta(x,N^{\prime}) \nonumber \\ & - \sigma_{\alpha\beta}(x,N,N^{\prime})w_\alpha(x,N) \big] {\rm d}N^{\prime}, \nonumber \end{align} \tag{ 4 }$$

where the integration is over the refractive index space. Here the scattering cross-section

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \nonumber \sigma_{\alpha\beta} (x,N,N^{\prime}) =&\, 2\pi k_0 \big\vert e^*_\alpha(x,N) (\mathds{1}-\varepsilon) e_\beta(x,N^{\prime})\big\vert ^2 \nonumber \\ \label{eq:sigma} & \times \Big(\frac{k_0}{2\pi}\Big)^3 \Gamma(k_0,x, N-N^{\prime}) \delta\big(H_\beta(x,N^{\prime})\big), \nonumber \end{align} \tag{ 5 }$$

describes a momentum-preserving three-wave interaction process in which a mode α carrying refractive index N interacts with turbulence exchanging momentum $\Delta N$ and thus scattering into a mode β with refractive index $N^{\prime}$, with momentum conservation $N - N^{\prime} = \Delta N$. The spatial position is not changed in the scattering event. The variation in momentum $\Delta N$ is provided by the turbulence spectrum Γ which is rigorously defined as the Wigner transform of the two-point correlation function of the relative electron density fluctuations $\delta n_e$, namely,

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:gamma} \Gamma(k_0,x,N) = \frac{1}{n^2_e(x)} \int {\rm e}^{-{\rm i}k_0 N \cdot s} \mathbb{E}\big(\delta n_e(x + s/2) \delta n_e (x-s/2)\big){\rm d}s, \nonumber \end{align} \tag{ 6 }$$

where n_e is the average electron density profile of the plasma and s is an integration parameter. Here the term inside $\mathbb{E}$ corresponds to the two-point correlation function of the electron density fluctuations.

In the WKBeam code, the two-point correlation function of the electron density fluctuations is a critical input quantity. Such information is usually not easily recovered, neither from experimental data nor from simulations. Therefore we make use of a model which is based on the current theoretical understanding of turbulence in a magnetized plasma, thus reducing the required inputs to quantities that can be inferred from available data. The model used in this work amounts to

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:basic_model} \frac{1}{n^2_e(x)} \mathbb{E}\big(\delta n_e(x + s/2) \delta n_e (x-s/2)\big) = F(x)^2 {\rm e}^{-\frac{1}{2} s \cdot A(x)s}, \nonumber \end{align} \tag{ 7 }$$

where $F(x)$ is a scalar function accounting for the spatial localization of turbulent density fluctuations and

$$\begin{align*} \newcommand{\e}{{\rm e}} \displaystyle A_{ij}(x) = L_\perp^{-2} \big(\delta_{ij} - b_i(x) b_j(x) \big) + L_{\Vert }^{-2} b_{i}(x)b_{j}(x), \nonumber \end{align*}$$

is a tensor accounting for the anisotropy of turbulence in magnetized plasmas, $L_\perp$ and $L_{\Vert }$ being the correlation lengths in perpendicular and parallel directions with respect to the unit vector b of the equilibrium magnetic field. The Gaussian shape arising from broad-band (white noise) signals is a good approximation for a range of plasma conditions as supported by experimental measurements, e.g. [28–30]. The input parameters for the code are thus reduced to the function F, which is usually specified on the poloidal plane of the tokamak and the two correlation lengths, see section 2.3 below. The function F in particular can be identified with the root-mean-square density fluctuations: In fact, we have, upon setting $s=0$,

$$\begin{align*} \newcommand{\e}{{\rm e}} \displaystyle F(x) = \sqrt{\frac{\mathbb{E}(\delta n^2_e(x))}{n^2_e(x)}} \equiv \Big[\frac{\delta n_e}{n_e}\Big]_{{\rm rms}} (x), \nonumber \end{align*}$$

and the average here is a statistical average (no spatial averaging should be performed). A further advantage of this model given by equation (7) is that the integral in spectrum Γ can be computed analytically.

The effect of the scattering operator on the solution of the wave kinetic equation can be rather complicated and it is in general not just a diffusion process. We can in fact distinguish two radically different regimes (as shown in section 2.5), namely, a superdiffusive regime in which scattering deforms the tails of the wave energy distribution only with little consequences for the quality of the beam, and a diffusive regime in which the whole wave energy distribution spreads with a significant broadening of the beam. An indication of the scattering regime can be obtained from the total scattering cross-section

$$\begin{align*} \newcommand{\e}{{\rm e}} \displaystyle \Sigma_\alpha (x,N) = \sum_\beta \int \sigma_{\alpha\beta}(x,N,N^{\prime}) {\rm d}N^{\prime}, \nonumber \end{align*}$$

which is the average number of scattering events per unit of propagation length for the mode α and near the phase-space point $(x, N)$. If $\newcommand{\e}{{\rm e}} \Delta \ell$ is the distance traveled by a beam in the turbulent plasma, then

$$\begin{align*} \newcommand{\e}{{\rm e}} \displaystyle \lambda_\alpha = \langle \Sigma_\alpha \rangle \Delta \ell, \nonumber \end{align*}$$

where $\langle \Sigma_\alpha \rangle$ is the average of $\Sigma_\alpha(x, N)$ along the beam path (reference ray), gives an estimate of the average number of scattering events for that beam. For $\lambda_\alpha \gg 1$ we expect to have diffusive solutions with significant broadening of the beam. It is therefore useful to have simple estimates for such a parameter.

We begin with the analytical calculation of the spectrum Γ, which in view of equation (7) amounts to

$$\begin{align*} \newcommand{\e}{{\rm e}} \displaystyle \Big(\frac{k_0}{2\pi} \Big)^3 \Gamma(k_0, x, N) = \frac{F^2(x)}{(2\pi)^{3/2} \sigma_\perp^2 \sigma_{\Vert }} \exp\left(-\frac{N_\perp^2}{2\sigma_\perp^2} - \frac{N_{\Vert }^2}{2\sigma_{\Vert }^2}\right), \nonumber \end{align*}$$

where we have introduced the parameters $\sigma_\perp = (k_0L_\perp){\hspace{0pt}}^{-1}$ and $\sigma_{\Vert } = (k_0 L_{\Vert }){\hspace{0pt}}^{-1}$ that play the role of spectral widths and cylindrical coordinates defined by

$$\begin{align*} \newcommand{\e}{{\rm e}} \displaystyle N = N_{\Vert } b(x) + N_\perp (e_1(x) \cos \phi_N + e_2(x) \sin \phi_N), \nonumber \end{align*}$$

with $(e_1, e_2, b)$ being an orthogonal reference frame constructed around the unit vector b of the local magnetic field.

Turbulence in a magnetized plasma is well within the limit $k_0 L_{\Vert } \gg 1$ both because of the high frequency of the wave and the elongated structures of the turbulence along the magnetic field direction ($L_{\Vert }\gg L_{\perp}$). Hence we can consider the limit $\sigma_{\Vert } \to 0^+$, for which one has

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:Gamma_long_Lparallel} \Big(\frac{k_0}{2\pi} \Big)^3 \Gamma(k_0, x, N) \to \frac{F^2(x)}{2\pi \sigma_\perp^2} \exp\Big[-\frac{N_\perp^2}{2\sigma_\perp^2}\Big] \delta(N_{\Vert }). \nonumber \end{align} \tag{ 8 }$$

The presence of the Dirac's delta-function together with two key properties of the cold-plasma dispersion tensor allow us to compute the total cross-section exactly. The key properties are: (i) The dispersion relation $H_\alpha(x, N) = 0$ corresponds to $N_\perp = n_\alpha(x, N_{\Vert })$ with $n_\alpha$ a known function; (ii) The dispersion tensor and thus the polarization vectors $e_\alpha(x, N)$ depend only on $N_\perp, N_{\Vert }$ and not on the polar angle $\phi_N$. Particularly (i) implies

$$\begin{align*} \newcommand{\e}{{\rm e}} \displaystyle \delta(H_\beta(x,N^{\prime})) = \left\vert \frac{\partial H_\beta}{\partial N_\perp} \right\vert ^{-1} \delta(N^{\prime}_\perp - n_\beta(x, N^{\prime}_{\Vert })). \nonumber \end{align*}$$

In view of such properties and in the limit $\sigma_{\Vert } \to 0^+$, we obtain

$$\begin{align*} \newcommand{\e}{{\rm e}} \displaystyle \Sigma_\alpha = \frac{k_0 F^2(x)}{\sigma_\perp^2} \sum_\beta M_{\alpha\beta} \int_0^{2\pi} {\rm e}^{-(n_\alpha n_\beta/\sigma_\perp^2)(1-\cos(\phi_N - \phi^{\prime}_N))} {\rm d}\phi^{\prime}_N, \nonumber \end{align*}$$

where $\Sigma_\alpha$ has been evaluated 'on-shell', i.e. on the dispersion surface $N_\perp = n_\alpha(x, N_{\Vert })$, and we have defined the matrix of coefficients

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:coefficients} M_{\alpha\beta} = n_\beta \left[\frac{\vert e^*_\alpha(x,N) (\mathds{1} - \varepsilon) e_\alpha(x,N^{\prime})\vert ^2} {\vert \partial H_\beta (x,N^{\prime}) / \partial N_\perp\vert } \right] \exp\left(- \frac{(n_\alpha - n_\beta)^2}{2\sigma_\perp^2} \right), \nonumber \end{align} \tag{ 9 }$$

where both $n_\alpha$ and $n_\beta$ should be evaluated at $(x, N_{\Vert })$, while the term in square brackets should be evaluated at $N^{\prime}_{\Vert } = N_{\Vert }$ and $N^{\prime}_\perp = n_\beta(x, N_{\Vert })$.

The remaining integral can be computed by means of the Jacobi–Anger expansion,

$$\begin{align*} \newcommand{\e}{{\rm e}} \displaystyle {\rm e}^{{\rm i}z\cos\theta} = \sum_{n \in \mathbb{Z}} i^n J_n(z){\rm e}^{{\rm i}n\theta}, \nonumber \end{align*}$$

where $J_n(z)$ are the Bessel functions and the sum is over the set $\mathbb{Z}$ of integer numbers. We can use this identity with $z=-{\rm i}n_\alpha n_\beta / \sigma_\perp^2$ and $\theta = \phi_N - \phi^{\prime}_N$ and only the term $n=0$ gives a contribution. At last we obtain

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:exactSigma} \Sigma_\alpha = \frac{2\pi k_0}{\sigma_\perp^2} \Big[\frac{\delta n_e}{n_e}\Big]^2_{{\rm rms}} \sum_\beta M_{\alpha\beta} {\rm e}^{-n_\alpha n_\beta /\sigma_\perp^2} I_0(n_\alpha n_\beta/\sigma_\perp^2), \nonumber \end{align} \tag{ 10 }$$

and $I_n(z) = i^{-n} J_n({\rm i}z)$ are the modified Bessel functions of first kind and, in particular, I₀ is even so that $I_0(-n_\alpha n_\beta / \sigma_\perp) = I_0(n_\alpha n_\beta / \sigma_\perp)$.

Expression (10) is strictly only valid under the assumption of long parallel correlation length $k_0L_{\Vert } \gg 1$ which is usually satisfied in cases of interest. If in addition, we consider the limit $\sigma_\perp = (k_0 L_\perp){\hspace{0pt}}^{-1} \to 0^+$, which is often the case for beams in the EC frequency range, the Hankel's large-argument expansion of the modified Bessel function,

$$\begin{align*} \newcommand{\e}{{\rm e}} \displaystyle I_0(z) \sim e^z /(2\pi z)^{1/2}, \nonumber \end{align*}$$

yields the simpler expression

$$\begin{align*} \newcommand{\e}{{\rm e}} \displaystyle \Sigma_\alpha \sim \sqrt{2\pi} k_0 \Big[\frac{\delta n_e}{n_e}\Big]^2_{{\rm rms}} \left[ \sum_\beta M_{\alpha\beta} (n_\alpha n_\beta)^{-1/2} \right] k_0L_\perp, \nonumber \end{align*}$$

which exhibits the linear dependence with $k_0L_\perp$.

With the aim of practical estimate of scattering, however, it is convenient to further simplify the expression for $\Sigma_\alpha$. In most cases the toroidal injection angles are modest (<20 degrees), hence we consider the simple case of perpendicular propagation. In the reference frame in which the local magnetic field defines the third axis and the refractive index belongs to the plane spanned by the first and third axis (Stix reference frame [31]), perpendicular propagation implies $N = (N_\perp, 0, 0)$. For the ordinary mode $\alpha=O$ we find

$$\begin{align*} \newcommand{\e}{{\rm e}} \displaystyle e_O(x,N) = \left(\begin{array}{@{}c@{}} 0 \nonumber \\ 0 \nonumber \\ 1 \end{array}\right), \qquad H_{O} = N_\perp^2 - 1 + \frac{\omega_{\rm pe}^2}{\omega^2}, \qquad e_O^* (\mathds{1} - \varepsilon) e_O = \frac{\omega_{\rm pe}^2}{\omega^2}, \nonumber \end{align*}$$

where $\omega_{\rm pe}$ is the electron plasma frequency. As a result, for perpendicular propagation,

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:Moo} M_{OO} = \frac{1}{2} \frac{\omega_{\rm pe}^4}{\omega^4}. \nonumber \end{align} \tag{ 11 }$$

Even under such simple conditions, the calculations for the extra-ordinary mode $\alpha=X$, although doable, yield rather cumbersome results. In fact one has

$$\begin{align*} \newcommand{\e}{{\rm e}} \displaystyle e_X(x,N) \propto \left(\begin{array}{@{}c@{}} \frac{D}{S+H_X} \nonumber \\ i \nonumber \\ 0 \end{array}\right), \quad H_{X} = \frac{N_\perp^2}{2} - S + \sqrt{\frac{N_\perp^4}{4} + D^2}, \nonumber \end{align*}$$

with the proportionality constant for e_X being determined by $\vert e_X\vert ^2 = 1$ and

$$\begin{align*} \newcommand{\e}{{\rm e}} \displaystyle e_X^* (\mathds{1} - \varepsilon) e_X = 1 - S \frac{S^2 - D^2}{S^2 + D^2}, \nonumber \end{align*}$$

where

$$\begin{align*} \newcommand{\e}{{\rm e}} \displaystyle S = 1 - \frac{\omega_{\rm pe}^2}{\omega^2 - \omega_{\rm ce}^2}, \qquad D = \frac{\omega_{\rm ce}}{\omega} \frac{\omega_{\rm pe}^2}{\omega^2 - \omega_{\rm ce}^2}, \nonumber \end{align*}$$

are the standard parameters of the cold-plasma dielectric tensor and $\omega_{\rm ce} >0$ is the electron cyclotron frequency (defined with the elementary charge so that it is a positive number). The dispersion relation $H_X=0$ in this case corresponds to

$$\begin{align*} \newcommand{\e}{{\rm e}} \displaystyle N_\perp^2 = (S^2 - D^2) / S, \nonumber \end{align*}$$

as expected, but the Hamiltonian H_X has been computed as the eigenvalue of the dispersion tensor. Then we have

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:Mxx} M_{XX} = \frac{1}{2} \frac{[S^2 + D^2 - S(S^2 - D^2)]^2}{S^2(S^2 + D^2)}. \nonumber \end{align} \tag{ 12 }$$

In this specific case of perpendicular propagation, the cross-polarization terms $M_{OX}$ and $M_{XO}$ are both exactly zero, no matter how close the mode dispersion surfaces are. This is due to the block-diagonal structure of the matrix $\mathds{1} - \varepsilon$ with the polarization e_O and e_X being in the subspace of the two different diagonal blocks of the matrix. The fact that modes are exactly decoupled for perpendicular propagation independently of the fluctuation level can also be inferred directly by inspection of Maxwell's equations.

The expression for $M_{XX}$ is still rather complicated for practical use. Instead of an exact calculation one could consider an easier lower bound. This can be computed on noting that the maximum and minimum values of any quadratic form $F(z) = z^* Az$, where z is a complex vector satisfying $z^*z=1$ and A a complex positive semi-definite Hermitian matrix, amount to the maximum and minimum eigenvalues, respectively. We can therefore obtain a bound for $e_\alpha^*(\mathds{1} - \varepsilon)e_\beta$ as well as for $\partial H_\alpha / \partial N_\perp = e_\alpha^* \partial D/\partial N_\perp e_\alpha$. Since we are interested in the X-mode for perpendicular propagation, we can limit our calculation to the blocks of the matrices orthogonal to the O-mode polarization and computing the eigenvalues we obtain

$$\begin{align*} \newcommand{\e}{{\rm e}} \displaystyle \frac{\omega_{\rm pe}^2}{\omega (\omega + \omega_{\rm ce})} \leqslant e_X^*(\mathds{1} - \varepsilon)e_X \leqslant \frac{\omega_{\rm pe}^2}{\omega (\omega - \omega_{\rm ce})}, \qquad 0 \leqslant \frac{\partial H_X }{\partial N_\perp} \leqslant 2N_\perp. \nonumber \end{align*}$$

The lower bound on the derivative of H_X cannot be used as we need $\vert \partial H_X / \partial N_\perp\vert ^{-1}$. That however allows us to obtain a lower bound for $\Sigma_X$. In summary we have

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:OmodeExact} \Sigma_O = \sqrt{\frac{\pi}{2}} k_0 \Big[\frac{\delta n_e}{n_e}\Big]^2_{{\rm rms}} \frac{\omega_{\rm pe}^4}{\omega^4} \frac{k_0L_\perp}{n_O}, \nonumber \end{align} \tag{ 13 }$$

for the ordinary mode while, with the lower bound for $M_{XX}$,

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:XmodeLowerBound} \Sigma_X \geqslant \sqrt{\frac{\pi}{2}} k_0 \Big[\frac{\delta n_e}{n_e}\Big]^2_{{\rm rms}} \frac{\omega_{\rm pe}^4}{\omega^2 (\omega + \omega_{\rm ce})^2} \frac{k_0L_\perp}{n_X}, \nonumber \end{align} \tag{ 14 }$$

for the extra-ordinary mode, both results being valid for perpendicular propagation under the limits $k_0 L_{\Vert }, k_0 L_\perp \gg 1$.

A simple estimate of the parameter $\lambda_\alpha$ can now be readily written by using either equation (13) or (14) as relevant to the considered mode. In order to compute $\lambda_\alpha$ rigorously we should average $\Sigma_\alpha$ over the propagation path of the beam in the turbulent plasma. For simplicity however we can evaluate $\Sigma_\alpha$, which for perpendicular propagation depends on x only, at a given point in space chosen to represent the typical value of plasma parameter in the turbulent region. Then we can write an estimate for $\lambda_\alpha$ which can be evaluated for a given tokamak scenario without need of ray- or beam-tracing calculations. That is,

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:lambda} \lambda_\alpha \approx \sqrt{\frac{\pi}{2}} \frac{1}{n_\alpha} \Big[\frac{\delta n_e}{n_c}\Big]^2_{{\rm rms}} (k_0^2 L_{\perp} \Delta \ell) \left\{\begin{array}{@{}ccc@{}} 1, & \; & \alpha = O, \nonumber \\ \omega^2 /(\omega + \omega_{\rm ce})^2, & \; & \alpha = X, \end{array} \right. \nonumber \end{align} \tag{ 15 }$$

and we have introduced the cut-off density $n_c=n_e(\omega/\omega_{\rm pe}){\hspace{0pt}}^2$. The authors would like to remind that this chain of approximations and simplifications were made to obtain a simplified expression for the categorization of the transport regimes. In the implementation of the code, the general expression (4) is evaluated during a scattering event as described in [32].

The turbulence model employed in WKBeam simulations consists in a model profile for the function F defined in equation (7), which is the root-mean-square relative electron-density fluctuations. We model F as a function of the radial coordinate $\rho_p=\sqrt{\psi_p}$, where $\psi_p$ is the normalized poloidal flux ($\psi_p=0$ and $\psi_p = 1$ corresponds to magnetic axis and the separatrix, respectively), and the geometrical poloidal angle θ. At first we shall assume that F is a purely radial profile, i.e.

$$\begin{align*} \newcommand{\e}{{\rm e}} \displaystyle F(\rho_p, \theta) = F_r(\rho_p), \nonumber \end{align*}$$

with F_r depending only on the radial coordinate $\rho_p$. The specific profile used for F_r was selected based on the experimental work in [33, 34] and it is shown in figure 1(a). It is a piecewise defined function combining core region until $\rho_p=0.97$, a (linear) ramp-up region at $\rho_p=0.97-1.0$ and a constant edge region at $\rho_p>1.0$. The core value in H-modes is a few percent (and has a negligible effect on the beam broadening), here it is chosen to be 2%. While the edge level dominates the beam broadening and here it was selected to be 20%. A Gaussian centred at $\rho_p=1$ is shown for comparison, as it has been used in the earlier works [11, 12]. Note that even though the function has rather high values outside separatrix, the fluctuation density itself is rather small due to small electron density. This is further illustrated in figure 1(b). As this extrapolation from current experiments to ITER is definitely uncertain, a scan for the edge value is presented in section 3.3. The correlation length was assumed to scale like $L_c\approx 5$–$10 \rho_{s}$ [35], where $\rho_{s}$ is the sound gyroradius. The authors are aware that this scaling might break down when going from open to closed field lines. As the temperature of the plasma is a radial function, so is the expected correlation length. In this paper, however, a constant value of the correlation length was used as the main effect of the fluctuations can be isolated to a radial domain close to the edge (we note, however, that the theory allows a radial dependency on the correlation length). Using ITER parameters, the correlation length in the studied scenario at the location of interest is 1-2 cm. Consequently, a value of 2 cm was used and a scan around this value is presented in section 3.3.

Zoom In Zoom Out Reset image size

In the plasma core, the turbulence exhibits strong ballooning structure leading to poloidal dependency in the density fluctuation amplitude. According to recent 2D measurements combined with gyrokinetic modeling [36] suggest this ballooning behavior can be modeled (apart from fine-structure assumed to be not particularly important for the current application) using

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:ballooning} F(\rho_p, \theta) = F_r(\rho_p)\frac{1}{2}\left(1+\cos\left(\theta\right)\right), \nonumber \end{align} \tag{ 16 }$$

where θ is the geometrical poloidal angle. For ITER, this ballooning leads to the shape illustrated in figure 2. In this figure the radial function shown in figure 1 was used. At the location where the UL EC beam enters the plasma, the value of the fluctuation amplitude is roughly half of the midplane value. However, only very marginal amount of information exists about the ballooning in open field lines either experimentally or theoretically [37–39], and even these observations are often of limited interest for us (limiter plasmas, L-mode). To compensate for these uncertainties, a scan over different fluctuation amplitude levels is performed, as detailed in the following sections.

Zoom In Zoom Out Reset image size

To ensure that the implementation of the theory is handled properly, several benchmarks have been carried out. These included a benchmark study without the fluctuations against well benchmarked and validated TORBEAM code [27]. Several different scenarios were studied, including various ASDEX Upgrade discharges and ITER. These cases involve variation in kinetic profiles, magnetic fields and also injection angles. Both the beam shape and absorption profiles were compared. In most of the cases, for which an ASDEX Upgrade example is shown in figure 3 and ITER examples in figures 4 and 5, both the beam and the profiles are extremely close to each other. However, some discrepancy in the absorption profiles were observed in particular for large toroidal injection angles like shown in figure 5. The reason for these discrepancies is thought to be a combination of physical differences between the codes (TORBEAM adopts the paraxial approximation, so that all the relevant information is computed on the central ray only) and the projection to flux surfaces necessary on the TORBEAM side when constructing the absorption profile.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Furthermore, the code was benchmarked against full-wave code IPF-FDMC [40]. In these studies, the effect of fluctuations on the beam shape was studied under ITER-relevant conditions, i.e. $n_e/n_c \approx 0.2$. The main goal of this benchmark was to find the threshold above which the Born approximation, as introduced in the seminal paper by Karal and Keller [21] and applied to the wave-kinetic equation by McDonald [20], breaks down. See [17] for details in the implementation in WKBeam. This threshold is expected to be a function of the fluctuation amplitude and, therefore, a scan against this variable was carried out until the codes disagreed. The main result is that the WKBeam method agrees well with the full-wave calculations in the levels of fluctuations amplitudes expected in ITER, i.e. $\delta n_e/n_e< 40$% as shown in figure 6. The beam broadening is here defined as the ratio of the fitted Gaussian widths of the beams with and without scattering. The WKBeam solution fits nicely on a second order polynomial fit represented by a solid curve in this figure. Therefore the basic physics assumption upon which the WKBeam code relies appears valid under ITER-like conditions.

Zoom In Zoom Out Reset image size

As an example of present-day experiment, we take ASDEX Upgrade discharge #25 485. There the EC beam was in X-mode polarization with a frequency of 140 GHz. Using the information of the plasma quantities, we may construct a 2D map of the scattering parameter in equation (15) as a function of $\newcommand{\e}{{\rm e}} k_0L_{\perp}\Delta \ell$ and $\delta n_e/n_c$ as motivated by this equation. For this analysis, we have set a Gaussian layer around the separatrix, as illustrated in figure 1 and calculated the distance s the beam travels in this layer. The nominal values for the correlation length was 0.4 cm for AUG and 2 cm for ITER and the nominal amplitude was $\delta n_e/n_e=10\%$. In figure 7 we show the $\lambda=1$ curve together with the nominal values of λ for the parameters discussed above. Clearly ITER is well in the diffusive regime, while AUG nominal values are close to boundary but on the super-diffusive side.

Zoom In Zoom Out Reset image size

For illustrative purposes, we considered also modified cases with $L_{\perp}=2$ cm and $\delta n_e/n_e=20\%$ for AUG and $L_{\perp}=0.2$ cm for ITER. With these parameters, the scattering parameter ends up clearly on the diffusive side in the AUG case and close to the boundary in ITER case. As the derivation of the scattering parameter includes several approximations, the boundary cannot be considered to be sharp. Moreover, the real transport, within the limits of our model, can be obtained by running the code and it was done for the two cases discussed above. The resulting statistically averaged electric energy density $\propto \mathbb{E}(\vert E\vert ^2)$ after the fluctuation layer as a function of a distance along the perpendicular distance to the EC beam is illustrated in figures 8 and 9. The diffusive beam shape in both cases can accurately be fitted with a Gaussian while the superdiffusive case has characteristic high tails, that suggest a Kappa or Lévy distribution [41, 42]. The ASDEX case is shown in figure 8. Red curve is a solution with out fluctuations, blue curve the solution using nominal parameters that produce superdiffusive behavior and green curve a modified turbulence case that produce a diffusive solution. For the latter two, points represent a fit (as described below). The ITER case is shown in figure 9. Here blue line is the case with out fluctuations, red line is the nominal diffusive case and green one is the superdiffusive example obtained using modified turbulence parameters. Again, the dots corresponds to fits.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

To further illustrate the transport regimes, we fitted for the superdiffusive cases a generalized Cauchy distribution with

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle E^2_{{\rm fit}}(x)=\frac{a^2}{\left(x^2+b^2\right)^c}, \nonumber \end{align} \tag{ 17 }$$

where $a, b$ and c are free parameters of the fit. Using the beam shape coming from WKBeam to generate the least-square fit, c has a value of $2.18\pm 0.20$ (95% confidence level), which is very close to Cauchy distribution ($c=2$). We have also counted the number of scattering events along each ray. In the diffusive regime, almost all rays ($>$99%) have scattered at least once and most of them several times, whilst in the superdiffusive regime only a fraction of the rays undergo a scattering event (depends on how deep in the superdiffusive regime one lies, but here for example roughly 20%) and having more than one event for a ray is very unlikely. Last remark is the difference between the superdiffusive beam shapes of ITER and AUG examples, in ITER the high tails are more visible than in the AUG case. This can be explained by the fact that in AUG case the fluctuation layer is so short, and the scattering probability correspondingly low, that not enough scattering events occur to raise the tails, and the original Gaussian beam shape is retained rather well.

We have concentrated in this paper on the 15 MA $Q=10$ standard H-mode ITER scenario. It is the most important mission of ITER and thus it deserves special attention. Here we have modelled the flat-top phase according to the simulations presented in [43]⁴. The relevant 1D profiles of electron density and temperature are shown in figure 10 and flux surfaces of the equilibrium in figure 11 (the distorted shape around the X-point is due to an equilibrium re-processing and does not affect the results in the top region, where the beam propagates).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Relevant parameters for the EC beams are given in table 1. Four different configurations are listed, namely the upper steering mirror (USM) and lower steering mirror (LSM) and both are separately optimized for current drive at $q=2$ and $q=3/2$ surfaces, where the most dangerous NTMs are expected to occur in the given scenario.

Table 1. Most important beam parameters for USM and LSM including both (2, 1) and (3, 2) resonance surfaces.

Quantity	USM ($q=2$)	LSM ($q=2$)	USM ($q=3/2$)	LSM ($q=3/2$)
Frequency (GHz)	170	170	170	170
Mode	O	O	O	O
Input power (MW)	1	1	1	1
Beam width (cm)	5.047	4.813	5.047	4.813
Focal length (cm)	318.6	200.1	318.6	200.1
Antenna R (cm)	699.9	705.4	699.9	705.4
Antenna z (cm)	441.4	417.8	441.4	417.8
Toroidal angle (deg)	20	20	20	20
Poloidal angle (deg)	46.8	40.5	54.7	50.1
Resonance surface ($\rho_p$)	0.87	0.87	0.77	0.77

In figure 12 we show the absorption profiles of electron cyclotron beams. The set of profiles on the left-hand side of the figure (around $\rho_p=0.77$) corresponds to deposition on the $q=3/2$ surface while the set on the right-hand side (around $\rho_p=0.87$) corresponds to the $q=2$ surface. Blue/black/red colour represents the case with 0/10/20% density fluctuation amplitudes, respectively, while the solid line shows the USM configuration and the dashed line the LSM configuration. An immediate observation is that the beams are broadened by a factor of  ≈3 in all cases. Another peculiar feature of the scattered profiles is the non-Gaussian shape. This can be explained by a combination of beam-broadening-induced geometrical factors and aberration effects (these are naturally included in the model as each ray carries information of the refractive index N and also of the temperature that affects the absorption). The effect of density fluctuations is similar for both mirrors and both resonance surfaces, while the geometric distortion of the Gaussian shape is slightly stronger for the inner $q=3/2$ surface compared to the outer $q=2$. Each of these simulations was performed using 4$\times 10^5$–2$\times10^6$ rays, which allows a statistically well-converged determination of the power deposition profiles.

Zoom In Zoom Out Reset image size

One major uncertainty of the density fluctuations on open field lines is related to ballooning. In the core region with nested flux surfaces, the fluctuations are strongly damped at the high-field side but only very scarce information is available for the open field line region of H-mode plasmas that is interesting for the present study. To evaluate the effect of ballooning on the profile broadening, we have assumed a core-like ballooning given by equation (16). A simulation was carried out for the USM $q=(2, 1)$ configuration, assuming 20% fluctuation amplitude at the outer midplane. For comparison, the results are shown together with 20% and 10% fluctuation amplitudes without ballooning in figure 13. As expected, the ballooning case is very close to the 10% fluctuation amplitude case as the level of fluctuations at the location where the EC beam enters the plasma is roughly half of the midplane value, in this case 10%.

Zoom In Zoom Out Reset image size

Since some of the key input parameters have large uncertainties, a scan has been performed to assess the sensitivity of the results on these parameters. As a basis for this scan, we have selected the USM configuration for the $q=2$ resonance surface. The scanned variables are the fluctuation amplitude $\delta n_e/n_e$ and the correlation length $L_c=L_{\perp}$ shown in figure 14. In this figure we show the full width at half-maximum (FWHM) of the deposition profile compared to the quiescent case. It is readily noticed that the effect of fluctuations is to broaden the beam by a factor of 1.7 to 5.1, in the range of most likely fluctuation amplitudes by a factor of 2.5 to 3.5. In this scan no ballooning was assumed. In the case of correlation $L_{\perp}$ scan, an inverse proportionality is observed. Different properties of the turbulence in the open-field-line region are mimicked through different values of the correlation length. An investigation of the different statistical properties of the edge turbulence requires a separate study and is left for future work.

Zoom In Zoom Out Reset image size