Microwave beam broadening due to turbulent plasma density fluctuations within the limit of the Born approximation and beyond

Plasma turbulence, and edge density fluctuations in particular, can under certain conditions broaden the cross-section of injected microwave beams significantly. This can be a severe problem for applications relying on well-localized deposition of the microwave power, like the control of MHD instabilities. Here we investigate this broadening mechanism as a function of fluctuation level, background density and propagation length in a fusion-relevant scenario using two numerical codes, the full-wave code IPF-FDMC and the novel wave kinetic equation solver WKBeam. The latter treats the effects of fluctuations using a statistical approach, based on an iterative solution of the scattering problem (Born approximation). The full-wave simulations are used to benchmark this approach. The Born approximation is shown to be valid over a large parameter range, including ITER-relevant scenarios.


Introduction
Electromagnetic waves in the microwave range of frequencies are widely used in fusionrelevant experiments for heating and diagnostic purposes [1][2][3]. In tokamaks, they are employed among others for control and suppression of MHD instabilities like the sawtooth oscillation and the neoclassical tearing mode (NTM) [4]. These applications require a good localization of the deposited wave power [5]. In particular, NTMs can lead to a degradation of the confinement up to a disruption of the discharge [6]. Since NTMs are driven by small perturbations in the plasma current profile (more precisely in the bootstrap current profile) resulting in the formation of magnetic islands, one way to stabilize them is to drive currents at the islands' positions and restore the original current profile. This can be achieved by injecting microwaves in the electron cyclotron range of frequencies [7,8]. It requires however a precise spatial localization of the place of absorption as the current should be ideally driven within the islands [9].
To reach the magnetic surface on which the NTM develops, the injected microwave beam has to pass the plasma boundary where density fluctuations occur with amplitudes up to 100 % [10]. This can significantly distort the beam and thus spoil the good localization, strongly reducing the efficiency of the NTM stabilization. A sound understanding of this effect is mandatory in order to predict the effectiveness of the microwave beams for the control tasks described above.
The influence of edge plasma density fluctuations on injected microwaves has been studied with geometrical-optics tools in the 1980s in a fusion-relevant context when high-power microwave sources became available [11][12][13]. The topic has been brought back into focus by Tsironis in 2009 [14] which triggered a significant follow-up research looking into this problem using different techniques [15][16][17][18][19][20][21][22]. As a common agreement one can state that (a) substantial broadening of microwave beams due to edge plasma density perturbations is expected, (b) the situation in medium-sized tokamaks differs from large-scale tokamaks like ITER (due to differences in microwave frequency, size of turbulent structures, and propagation length), (c) further and more detailed studies with a minimum of simplifying assumptions are needed for the ITER scenarios which cannot be explored experimentally in today's tokamaks, and (d) the various numerical tools should be cross-benchmarked.
This paper contributes to the understanding of the interaction of microwaves with turbulent plasma density fluctuations with the aid of numerical simulations. Two different codes are used: the full-wave code IPF-FDMC [23] and the WKBeam code [24] which solves the wave kinetic equation in the presence of random fluctuations in the background density. While for the first code an ensemble of different realizations of turbulent density fluctuations is required to reproduce the situation in the experiment, WKBeam allows to directly calculate the average effect by applying a statistical operator. The derivation of this scattering operator is based on the so-called Born approximation [25]. The WKBeam results are thus expected to become invalid for sufficiently high fluctuation levels. These limitations of the latter treatment are explored and quantified.
This paper is the continuation of previous full-wave simulations of scattering from singular blob-like density structures [26] and of first simulations including turbulent electron density fluctuations [27]. It serves as a benchmark for the WKBeam code and the recently published results [28,29]. The paper is organized as follows: both numerical codes and the generation of the electron density fluctuations are described in Sec. 2. Section 3 describes the setup of the simulations, followed by Sec. 4 which explains how the data obtained from both types of simulations is analyzed. In Sec. 5, the influence of the level of the electron density fluctuations on beam broadening is investigated. The role of the absolute value of the background density is then discussed in Sec. 6 and the influence of the thickness of the fluctuation layer is discussed in Sec. 7. Results for changing the injected mode from ordinary (O) to the extra-ordinary (X) are presented in Sec. 8. The summary in Sec. 9 concludes the paper.

Numerics
This section describes the numerical tools which are used throughout the paper. First, the full-wave code is introduced in Sec. 2.1, followed by the WKBeam code in Sec. 2.2. Both codes are only briefly described and the interested reader is referred to the references given in the corresponding sections. Note that a cold plasma model is used here. The most dangerous NTMs in ITER are expected to occur at radial positions corresponding to electron temperatures of approximately 7 keV [28]. The effective refractive index changes only marginally for these temperature, see e.g. Ref. [1], and the corresponding effects on the microwave beam propagation are negligible compared to the effect of density fluctuations.

The full-wave code IPF-FDMC
IPF-FDMC [23] is a finite-difference time-domain (FDTD) [30] code solving Maxwell's equations and the fluid equation of motion of the electrons on a 2D Cartesian grid. It allows to simulate propagation of electromagnetic waves in a cold magnetized plasma. Specifically, the mathematical model considered consists in evolution equations for the magnetic field B, the electric field E and the current density J of the wave, in a plasma equilibrium with background magnetic field B 0 and electron density n e . Specifically, with c the speed of light, ω pe = n e e 2 /(ǫ 0 m e ) the electron plasma frequency, ω ce = |e|B 0 /m e the electron cyclotron frequency, andB 0 the unit vector into the direction of the background magnetic field. An electron collision frequency ν e is included in Eq. (3) as a dissipation mechanism. The code has been successfully benchmarked against cold plasma theory [31] and used to study mode conversion processes [23] and microwave heating in plasmas [32].
For the implementation of Eqs. (1)-(3), the standard FDTD scheme [30] has to be complemented by a discretization scheme for the current equation (3). Here we use a "straight forward" way, that is first advance J x using the old values of J y and J z (in the cross product with the background magnetic field), then advance J y using the updated value of J x and the old value of J z , and finally advance J z using the updated values of J x and J y . As has been shown in previous works, this method is completely sufficient for a rather large set of problems [23,32,33]. More advanced methods exist for situations with extreme density gradients and we would like to refer the interested reader to the detailed and thorough analysis by Heuraux and da Silva [34][35][36].

The WKBeam code and the Born scattering approximation
The WKBeam code [24] is based on the formalism of the wave kinetic equation [37,38] which describes the average effect of plasma density fluctuations on a traversing microwave beam. The derivation of scattering operator in the wave kinetic equation relies on the Born approximation which is expected to hold for a weakly scattering medium in the sense specified below. No restrictions need to be made on the spatial size of the turbulent density structures, where other methods based on the short-wavelength approximation (e.g. geometrical optics, beam tracing) fail in presence of short-scale fluctuations. The Born approximation imposes, however, a limitation on the amplitude of the turbulent density fluctuations, to be explored in detail in this paper. Since this point is essential for the following discussion, some details about the derivation of the WKBeam model are reported below. The reader is referred to Refs. [24,28] for more details.
In WKBeam, turbulence is described as a time-independent random field of density fluctuations, with the idea that the time average of a physical quantity of interests can be computed as the ensemble average over a sufficiently large number of samples of the random field. The wave beam is modeled again by equations (1)- (3). With the plasma frequency ω pe being a time-independent random field, we can Fourier transform in time.
In the frequency domain Eqs. (1)-(3) can be written as a single equation for the Fourier transformed wave electric fieldÊ =Ê(ω, x), namely, whereε 0 =ε 0 (ω, x) is the cold plasma dielectric tensor [2] computed with the unperturbed electron density n e,0 , I is the identity tensor, and δn e the random fluctuation field. We assume that the expectation value is δn e (x) = 0, and the correlation function δn e (x)δn e (x ′ ) is known. The Born approximation [25] consists in the iterative approximation of a solution of Eq. (4) of the form whereÊ 0 is a solution of the unperturbed problem (Eq. (4) with δn e = 0) and determines the correctors for j ≥ 1. Formally at least, the solution forÊ j is of order (δn e /n e,0 ) j so that we may expect convergence of the series for a small-enough fluctuation level, precisely for [ δn 2 e /n 2 e,0 ] 1 2 ω 2 pe /ω 2 ≪ 1. When the series converges, the wave energy density averaged over random fluctuations is proportional to for we have Ê 1 = 0, which follows from averaging the Eq. (6) with j = 1, while in general Ê j = 0, j ≥ 2. It is expected therefore that the deviation of |Ê| 2 from its unperturbed value |Ê 0 | 2 grows quadratically with the fluctuation strength in the Born scattering regime. The Born expansion (5) has been applied by Karal and Keller [39] in order to obtain an equation for the average wave field Ê (ω, x) and later Mc Donald [38] extended their method to derive an equation for the wave-field correlation function Ê (ω, x)Ê(ω, x ′ ) * , from which the radiative transfer model of WKBeam is obtained. Mc Donald's formal derivation applies to abstract wave equations of the form where D 0 is a linear operator acting on a vector ψ in an abstract Hilbert space, and D 1 is a linear operator with random coefficients. For the specific problem (4), the wave field ψ is the electric fieldÊ(ω, ·), D 0 ψ is the unperturbed operator ∇ × (∇ ×Ê) − k 2 0ε 0Ê , and D 1 ψ amounts to −k 2 0 δne n e,0 (ε 0 − I)Ê and includes random density fluctuations. The wave field is sought in the form ψ = ψ 0 + ψ 1 , where ψ 0 satisfies D 0 ψ 0 = 0. After shifting possible singularities in the complex plane [38], we can construct an operator G such that D 0 G = I. Then Eq. (7) implies ψ = ψ 0 − GD 1 ψ, and iterating, which is the Born expansion (5). We can use this series to evaluate the correlation operator ψψ * , and multiplying on the left by D 0 we have where both D 0 and D 1 are assumed to be Hermitian and the identity D 1 = 0 has been accounted for. At last, we observe that ψ 0 ψ * 0 differs from ψψ * by second-or higher-order terms, hence, The Weyl symbol of the correlation operator ψψ * is by definition the average Wigner matrix W = W (x, N) which is a Hermitian-matrix-valued function of position x and refractive-index vector N. Upon computing the Weyl symbol of the operator equation (9), the relevant equation for W is readily obtained in a form that depends only on the correlation functions of the random density field, and thus lends itself to the asymptotic solution in the short-wavelength limit in spite of the presence of short scale random fluctuations [37]. For the specific case of Eq. (4), the lowest order approximation of the Wigner matrix W in the short wavelength limit is diagonal on the basis of the two cold plasma polarization vectors and the corresponding two real eigenvalues w α are referred to as the Wigner functions of the ordinary (α = O) and extra-ordinary (α = X) modes. The dispersion relation imposes the constraint H α w α = 0 with H α = H α (x, N) being the geometrical optics Hamiltonian for the mode α. Then the equation for W reduces to the wave kinetic equation solved by WKBeam, namely, where the Wigner functions w α are the unknowns and the scattering term S α can be brought to the form S α = β S αβ with where σ αβ is the scattering cross-section [24,28]. Let us remark again that Eq. (10) holds independently on how short the correlation length and thus spatial scale of the fluctuations themselves is. On the other hand, the Born approximation imposes a limit on the fluctuation amplitude. The Born series (8) is controlled by a norm of the operator D 1 which for the case of Eq. (4) can be estimated by [ (δn e /n e,0 ) 2 ] 1 2 ω 2 pe /ω 2 , hence the Born approximation should remain valid even for large values of the relative root-meansquare amplitude of the fluctuations if the wave is propagating sufficiently far away from the cut-off, i.e. at sufficiently low densities or high frequencies. Before concluding this section, it is useful to quote a further result derived in [28], namely that the product Σ α ∆ℓ, where ∆ℓ is the distance travelled in the turbulent region, is found to scale as where n e,cut is the cut-off density, k 0 is the vacuum wave vector and L ⊥ is the perpendicular correlation length of the fluctuations. The quantity Σ α ∆ℓ gives an estimate of the number of scattering events experienced by a given ray. The code WKBeam has been successfully benchmarked against the paraxial WKB code TORBEAM [40] for different, fusion-relevant scenarios [28].

The simulation set-up
The 2D computational domain resembles part of a poloidal cross section in a toroidal magnetic confinement device. The emitting antenna is located in vacuum on the right hand side of the domain. A frequency of f 0 = 50 GHz, corresponding to a vacuum wavelength of λ 0 ≈ 6 mm, is chosen for the microwave beam which is described in detail in Sec. 3.1. After a propagation distance of 5 cm in vacuum, a linearly increasing plasma density profile is encountered, described in detail in Sec. 3.2. The background magnetic field is taken to be homogeneous across the whole domain with a purely toroidal direction and a strength of B tor = 1 T, corresponding to a normalized value of Y = ω ce /ω 0 ≈ 0.56.
Note that the absolute values of the frequency, the plasma density, and of the background magnetic field correspond to the values of the ASDEX-Upgrade tokamak [7,43] reduced by approximately a factor 2.5. The reduced frequency and hence increased vacuum wavelength allows for a coarser numerical grid to be used decreasing the required computational resources. Since the electromagnetic wave equation in a cold plasma depends on plasma density and magnetic field through non-dimensional parameters, X = ω 2 pe /ω 2 0 and Y = ω ce /ω 0 respectively, our simulations models however scenarios at the same X and Y as fusion-relevant scenarios.
The generation of synthetic density fluctuations is outlined in Sec. 3.3.

The injected microwave beam
The injected beam is Gaussian as being characteristic for typical fusion experiment [44]. Specifically, for a two-dimensional domain the electric field amplitude of a standard Gaussian beam [45] is given by with z the radial distance to the beam axis, x the axial distance to the beam waist, w 0 the size of the beam waist, R the radius of curvature of the wavefront, and φ 0 the Gouy phase shift [46,47]. Note that w, R and φ 0 are all functions of the axial distance x to the beam waist. (Also note that the radial distance is usually denoted with x, and the axial distance is usually denoted with z.) In the full-wave code, the field distribution is defined in the antenna plane explicitly by Eq. (12) and added to the electric field on the grid resembling a soft source [48]. A focusing beam with the waist located in front of the antenna plane (still inside the computational domain) is considered, and w and R need to be evaluated in the emitting antenna plane using given values of the beam waist w 0 and of the axial distance to the waist x. In WKBeam, in contrast, the parameters w and R in the antenna plane are direct input parameters. Values of w = 1.5 cm and R = 10 cm are used in WKBeam, corresponding to w 0 ≈ 9.7 mm and x ≈ 58.2 mm for the full-wave code as described in detail in the Appendix. Those values were chosen to ensure that the simulation domain contains the beam waist and the subsequently diverging beam and still has a reasonable size with respect to the required computational resources.

The electron plasma density profile
A 1D profile which depends only on the radial coordinate x is chosen for the background electron plasma density. The density starts to increase linearly at x n1 = 2.45 m until x n2 = 2.30 m, where a maximum value of n e,max = 0.2 · 10 20 m −3 ≈ 0.65 · n e,cut is reached, where n e,cut refers to the cut-off density of the injected mode which is, if not explicitly stated otherwise, the O-mode. The density values normalized to n e,cut correspond approximately to those in the scrape-off layer in ASDEX-Upgrade [49]. The linear profile is plotted in Fig. 1 and described by  A layer of turbulent density fluctuations is then added to the background profile, as described in Sec. 3.3. The Gaussian envelope of the fluctuation amplitude is given by the expression where A 0 is the normalized fluctuation strength, R 0 = 1.65 m and a = 0.6 m correspond respectively to the major and minor radius, and w fluct defines the width of the Gaussian envelope. The parameter x shift is used to shift the fluctuation layer radially, where a value of x shift = 0 corresponds to the center of the layer being located at x = 2.40 m. Note that the values used for major and minor radii correspond to the ASDEX Upgrade tokamak [43]. Figure 2 shows as an example one sample for the actual 2D density profile used in the full-wave simulations. The finite spatial extent of the fluctuations in x-direction can be clearly seen. It is remarked that in the region around x shift = 0 the ratio between plasma density and cut-off density closely matches the corresponding value expected in the ITER standard scenario [28].
To ensure statistical significance of the full-wave simulations, the ensemble of density profiles needs to be large enough for each set of turbulence parameters. A size of N = 3000 turbulence realizations has been found to yield relevant results as will be demonstrated in Sec. 5. In WKBeam on the other hand, the statistical parameters of the turbulent density fluctuations are used as an input to directly calculate the average effect on the microwave beam. This leads to a significant reduction of computational time as compared to full-wave simulations.
The perpendicular correlation length of the density structures is set to a value of L ⊥ ≈ 5 mm which is close to the vacuum wavelength of the injected microwave (λ 0 ≈ 6 mm). According to Ref. [27], this can result in pronounced scattering of the microwave. The correlation length is predicted to scale like L ⊥ ≈ 5 − 10ρ s [50], with the drift scale parameter  with deuterium ions, the chosen value of L ⊥ = 5 mm lies within that range.

Density fluctuations
In order to study the effect of plasma density fluctuations on a traversing microwave beam by means of full-wave simulations it is necessary to let the microwave beam interact with an ensemble of density profiles, each being a sample of the same random field, and then average over the resulting wave electric fields. To generate a large number of individual samples, we use synthetic turbulence, as it allows us to generate the required large ensembles in a reasonable time as opposed to using large-scale plasma turbulence codes. It also allows us to ensure that the statistics of the random field is the same as that assumed in the WKBeam code.
The computational domain is defined on a 2D grid which is a reasonable simplification as turbulence in magnetized plasmas is highly anisotropic with typically very small wave numbers parallel to the background magnetic field [10]. The 2D domain corresponds approximately to a poloidal cross section in a toroidal magnetic confinement device. The full electron plasma density profile in the 2D simulation domain, described in detail in Sec. 3.2, can be written as where x and z are the radial and vertical coordinates, respectively, n e,0 (x) is the unperturbed background profile, F (x) an envelope of the fluctuations' amplitude basically defining their spatial location, and δn(x, z) is a random field such that δn e (x, z)/n e,0 (x) = F (x)δn(x, z) is the density fluctuation. Note that there is no dependence on time here as the density fluctuations appear to be frozen in the time frame of the microwave (also referred to as frozen plasma assumption): typical frequencies of the density fluctuations lie in the kHz range [10], whereas the microwave oscillates in the GHz range. In addition, for the densities considered here, the group velocity of the microwave is several orders of magnitude above the propagation speed of the density structures [10]. The fluctuations themselves are generated by a truncated sum of Fourier-like modes: with A i,j the amplitudes of the modes and ϕ i,j independent random phases uniformly distributed in the interval [0, 2π). The correlation length of the two-point autocorrelation function of the density fluctuations correspond to the average perpendicular structure size L ⊥ . Although the corresponding spectra in the experiments exhibit usually some kind of power law (see e.g. Refs. [41,42]) in contrast to the Gaussian shape used here, this is not expected to lead to significantly different scattering of the microwave beam for the parameters used here as the power laws differ most strongly from the Gaussian at large k-values for which, according to previous investigations [27], strongly reduced scattering is expected. In WKBeam, the effect of plasma density fluctuations is included via a scattering operator. The input parameters in the current model can be reduced to the spatial localization of the fluctuation layer, F (x), and the two correlation lengths, L ⊥ and L || (for details, see Ref. [28]). This ensures that both codes use the same plasma density profiles (including fluctuations) as input.

Data analysis
The full-wave simulations, based on a time-dependent scheme, start with the excitation of the microwave beam. They are stopped when a steady state solution is achieved. Including a safety margin in computational time, this corresponds to a value of T = 100 T wave , where T wave denotes the oscillation period. At various radial positions, the time-averaged squared wave electric field is recorded across the whole z-range of the computational domain: where t is the time coordinate and the tilde indicates that a scenario with turbulent density fluctuations is used (as opposed to a scenario without fluctuations where just the linear profile as described by Eq. (13) is used). Such detector antenna signals are acquired for each sample at a given set of radial positions x. The ensemble-averaged signals of the full-wave simulations can then be compared with the output of WKBeam (which yields directly the squared wave electric field).
As will be demonstrated in the following section, the ensemble-averaged beam cross-section can be approximately described by a Gaussian. In order to quantify the broadening of the injected microwave beam, it is thus convenient to fit a Gaussian of the shape to the (ensemble-averaged) scattered signal in the detector antenna planes with a 0 , a 1 , and a 2 being the fit parameters. The value obtained for the (averaged) beam width w, corresponding to the fit parameter a 2 , is then compared with the width w of the same beam propagating in the unperturbed scenario (i.e. without fluctuations but with the background profile). A normalized beam broadening b =w/w is obtained in this way for each set of turbulence parameters allowing to compare and benchmark the WKBeam results with the full-wave simulations. Note that another example for FDTD full-wave simulations of electromagnetic waves passing through random media consists in calculating the scattering coefficient, see e.g. Ref. [54]. The relevant physical quantity in our case is however the beam broadening as outlined in the introduction. Although the transverse ensemble-averaged beam profile can be well approximated by a Gaussian in most cases, there will be a few scenarios where a Cauchy distribution is more suitable (as will be discussed in the following Sections). Therefore, a general Cauchy distribution of the form will also be fitted to the detector antenna signals (via a non-linear least square fit), where 2a 1 corresponds to the full width at half maximum (FWHM) and a 2 to the median (a 0 fits the amplitude). The normalised beam broadening for those cases will be analysed in terms of the FWHM, i.e. the value of a 1 obtained from the ensemble-averaged signal is normalised with the corresponding value for the case without fluctuations.
The plasma density in the fluctuation layer can locally reach values close to the cutoff density of the microwave if large fluctuation levels are considered. This can result in microwave power being reflected. The corresponding quantity is routinely measured in the full-wave simulations and its value, averaged over the whole ensemble of fluctuations, is found to be at maximum 1 % for the worst case scenario, and at least one order of magnitude lower for most cases.

Influence of the fluctuation level on beam broadening
In this section, the results from full-wave simulations and WKBeam calculations are compared first for the case without turbulent density fluctuations (the zeroth case, comparison in vacuum, yielded excellent agreement and is not included in this paper). For this case we find that the two codes are indeed in very good agreement. Having established a reference solution, density fluctuations are included with their envelope localized at a radial position, as described in Sec. 3.2. The fluctuation amplitude is varied in a series of scans and the resulting values for the beam broadening are compared for the two codes.

Case without turbulence
The linear density profile given by Eq. (13)    of x = 2.35 m from both, the full-wave simulations and the WKBeam calculations, agree within the resolution of the plot shown in Fig. 3. The Gaussians fitted to the signals (see Eq. (18)) are also included. The value of the beam size w obtained from the Gaussian for WKBeam is larger by 0.2 % than the corresponding value for the full-wave simulations. It is in principle possible to further reduce this difference by adjusting the spatial resolution in both codes since the detector antenna signals need to be interpolated to the respective other code. This would however increase the total demand for computational resources and the introduction of the turbulent fluctuations leads anyway to an increase of at least an order of magnitude in the difference (see Sec. 5.2). We thus decided to use the values yielding the (already very good) agreement shown in Fig. 3.

Including a layer of turbulent density fluctuations
As a next step, a layer of turbulent plasma density fluctuations is added to the background profile. The position and the width of the layer are kept constant, using values of x shift = 0 and w fluct = 2 cm, respectively. To illustrate the effect of the fluctuations, a snapshot of the absolute wave electric field obtained from full-wave simulations for a single sample is shown in Fig. 4. The fluctuation layer clearly perturbs the injected beam, leading to a splitting into multiple beams which destroys the intended spatial localization of the absorption.
As an example, full-wave and WKBeam beam profiles are compared for a fluctuation amplitude A 0 = 0.3 in Fig. 5. On average, i.e. averaging over the full ensemble, a small broadening of the beam as compared to the case without turbulence is found. The ensemble-averaged signal resembles a smooth Gaussian-like beam, illustrating the sufficient size of the ensemble. Note that the signal is symmetric around z = 0 and thus does on average not change its original direction of propagation. Comparing with WKBeam, no differences are noticeable in this representation, proving the validity of the WKBeam calculations for this set of parameters. A tiny deviation to a Gaussian fit, also included in the figure, is only found in the tail of signals, which is slightly elevated. The elevated tails can become more pronounced, as will be discussed in Sec. 6. From the full-wave simulations, a broadening normalized to the case without fluctuations of b fw = 1.070 (±0.002) is obtained and for the WKBeam calculations it is b WB = 1.086 which is larger by approximately 1 %. Normalizing the broadening to the case without plasma (i.e. the vacuum case), b fw and b WB need to be multiplied with an additional factor of approximately 1.119, resulting in b fw,vac ≈ 1.197 and b WB,vac ≈ 1.215. If not mentioned explicitly otherwise, the beam broadening is normalized to the case without fluctuations (but with plasma) in the rest of this paper.

Scanning the fluctuation level
One of the main goals of this paper is to investigate the deviation of the beam broadening as predicted by the WKBeam code with respect to the reference solution provided by the full-wave solver. In particular, as discussed in Sec. 2.2, the Born approximation is expected to become inaccurate with increasing fluctuation level and background density. Figure 6 (left) shows the resulting scaling which yields an increased broadening with increasing fluctuation level. Error bars are not shown since the standard deviation of the ensemble-averaged beam broadening is smaller than the symbol size used in the plots.
The broadening for WKBeam is consistently larger than the full-wave solution and the absolute difference increases with increasing fluctuation strength. In WKBeam, the beam broadening b as a function of the fluctuation amplitude A 0 follows a power law. The solid line included in the plot in Fig. 6 (left) corresponds to a fit to the WKBeam values, obtaining a functional dependence of b = 0.96 · A 1.99 0 . For small fluctuation levels, a quadratic dependence is expected as the scenario resembles a phase screen [55]. The full-wave simulations, in contrast, exhibit a reduced increase for large values of A 0 . WKBeam is thus overestimating the broadening for large fluctuation levels which can be illustrated by the relative difference plotted in Fig. 6 (right). The maximum overestimation is with 6 % still considered to be small. It can become more significant if the background density is larger, as presented in the following section.

Influence of the turbulence radial location on beam broadening
In this section, the background density in the fluctuation layer is varied by shifting the layer radially, i.e. along the x-direction. To this end, the parameter x shift (see Eq. (14)) is varied in the range x shift = 0, . . . , 4 cm in steps of 1 cm, corresponding to background density values at the center of the fluctuation layer of n e /n e,cut =  Figure 7 shows as an example a full-wave simulation for a value of x shift = 4 cm and A 0 = 0.30 with the fluctuation layer centered around a position of x = 2.36 m. The increased background density at the fluctuation layer is expected to result in stronger scattering as refraction effects become more pronounced: power scattered by the turbulent density structures off the original direction of propagation experiences stronger refraction resulting on average in an increased beam broadening.
As a further example, the detector antenna signals for x shift = 4 cm and A 0 = 0.50 are plotted in Fig. 8. Since the antenna position used in the previous section (x = 2.35 m) would be situated inside the fluctuation layer, a position of x = 2.31 m is chosen here. As a first observation, the ensemble-averaged full-wave signal and the WKBeam signal are both broader (and correspondingly with a reduced peak amplitude) than the example shown in Fig. 5, where (a) the background density was lower and (b) the normalized fluctuation amplitude was with a value of A 0 = 0.3 also lower. Although both these differences make the scenario considered here a harder test for the Born approximation, the agreement between WKBeam and full-wave solution is still remarkably good.
Although a disagreement between full-wave and WKBeam signals can be seen, it is not considered to be significant.
Another observation is related to the shape of the signals: a Gaussian seems no longer be the adequate function to describe them, a pronounced elevation at the tails can be clearly seen. Therefore, a Cauchy distribution as described by Eq. (19) has also been fitted to the signals. Comparing the fitted Cauchy distribution with the fitted Gaussian, see Fig. 8, the first one seems to be more suitable to describe the broadened microwave beam for these fluctuation parameters. This finding corresponds to the results presented in Ref. [28]: a Cauchy distribution thus corresponds to a scattering process of superdiffusive nature, whereas a Gaussian shape corresponds to a diffusive process.
Note that the beam broadening deduced from either the Gaussian or the Cauchy  does not differ much: for the Gaussian fit it is b fw = 1.652 (±0.001) and b WB = 1.838 for full-wave and WKBeam, respectively, and for the Cauchy fit the values are b fw ≈ 1.60 and b WB ≈ 1.74. To get the broadening normalized to the vacuum case a slightly different factor as in Sec. 5 is required here due to the different detector antenna position: the factor is approximately 1.211, resulting in b fw,vac ≈ 2.00 and b WB,vac ≈ 2.23 for the Gaussian fit.
The broadening is no longer the only measure of interest. Instead, an increasing amount of energy is located in the tails. This could create a problem for an actual experiment, as it basically means that more scattering events far off the original direction of beam propagation occur, threatening diagnostics or other wall components.
The resulting beam broadening deduced from the Gaussians and the Cauchy distributions fitted to the detector antenna signals is shown, respectively, in Fig. 9  values of x shift (and increasing values of A 0 ) for both codes. Significant broadening by more than a factor of two is found. One can also see that the WKBeam values can still be represented by power laws and that the full-wave values are consistently smaller. Not much differences can be seen between using a Gaussian or the Cauchy distribution in this representation, the broadening seems to be very similar. An asymptotic behavior is observed towards larger values of x shift , i.e. higher background densities: the slope of the fitted power law for x shift = 4 cm is only slightly larger than for x shift = 3 cm whereas there is a substantial difference going from x shift = 0 cm to x shift = 1 cm. The overestimation of the beam broadening of WKBeam can be described in a more quantitative way by calculating the percentage deviation d of the WKBeam values to the full-wave values. Polynomials of 2 nd order can be fitted to the deviation as a function of fluctuation level A 0 and x shift (representing the background density). As shown in Fig. 10, maximum deviations of 18 % are found for the parameters used in this paper, where for fluctuation levels below 50 % the overestimation of WKBeam stays below 10 %.

Influence of the width of the fluctuation layer on beam broadening
The width of the fluctuation layer is varied in this section in order to investigate the influence of the propagation length inside of the fluctuating density area on beam broadening. The parameter w fluct , see Eq. (14), is varied, where values of w fluct = 1, 2, 3 cm are used. The parameter x shift is kept constant at a value of x shift = 0 and the fluctuation amplitude A 0 is varied over the same range as in the previous cases.
With increasing width of the fluctuation layer the fluctuating density extends to regions of higher density and thus lower x values. Therefore, the position of the detector antennas is, as in the last section, set to a position of x = 2.31 m. Figure 11  shows the beam broadening as obtained from fitting a Gaussian to the detector antenna signals as a function of A 0 with w fluct as parameter. Similar to the results presented so far, the WKBeam calculations yield larger values than the full-wave simulations, with increasing absolute deviation for increasing width of the fluctuation layer. The dependence of the broadening on A 0 is found to be stronger for larger value of w fluct , i.e. with increasing propagation length in the fluctuation layer, the broadening increases further with maximum beam broadening values of approximately two.
The percentage deviation of the WKBeam results with respect to the full-wave results does not show a strong dependence on w fluct , see Fig. 11 (right). The deviations for w fluct = 1 cm and w fluct = 2 cm as a function of A 0 are very similar. Only for w fluct = 3 cm a slightly stronger increase with increasing value of A 0 , i.e. a slightly steeper slope, is observed. This is, however, thought to be caused by the increased spatial extension of the fluctuation layer into regions with higher background densities (see Sec. 6) instead of an increased propagation distance in the fluctuation layer alone.
For most of the cases in this scan, the Gaussian provides the better fit to the detector antenna signals, only for large fluctuation amplitude cases (A 0 ≥ 0.5) at w fluct = 3 cm, the Cauchy distribution is the better approximation of the transverse beam profile.

O-and X-mode comparison
In the cases presented so far, the injected microwave beam was in O-mode polarization.
In this section, we present simulations results of WKBeam's capability of injecting a beam in X-mode polarization (by comparing and benchmarking with the corresponding full-wave simulations). Instead of repeating all scans presented so far, from which no further knowledge would be gained, we restrict ourselves to a scan of the fluctuation amplitude A 0 for fixed values of w fluct = 2 cm and x shift = 0.
To avoid the right-hand cut-off of the X-mode [56], we have reduced the background magnetic field to B tor = 0.25 T which results in a very similar background density profile when units normalized to the respective cut-offs are considered. Figure 12 shows the beam broadening deduced from the Gaussian fits to the detector antenna signals at a position of x = 2.35 m. The equivalent case for O-mode injection was shown in Fig. 6 and one can see that they are very similar. For the X-mode case, the power law fitted to the full-wave simulations yields for the normalized beam broadening as a function of the fluctuation amplitude b = 1.05 · A 2.05 0 which is, again, very similar to the O-mode case. The small difference is due to the slightly different normalized background density resulting in slightly different refraction.

Summary
We have investigated the broadening of a microwave beam passing through a layer of turbulent plasma density fluctuations, resembling the situation of a fusion edge plasma. The results from the WKBeam code were compared with full-wave simulations, performed with IPF-FDMC, over a large parameter range in order to benchmark WKBeam and explore the ranges of validity of the underlying approximations, and specifically the Born approximation. This approximation allows to directly calculate the effect of fluctuations in WKBeam by applying a scattering operator whereas the full-wave simulations require an ensemble-average. For the scenarios presented here, this leads to a speed-up of WKBeam of approximately a factor of 4 as compared to the full-wave simulations (scaling the actual wall-clock times of the computations down to a single process). This value will become significantly larger when increasing the wave frequencies which requires a higher spatial resolution and thus larger numerical grids in the full-wave simulations.
Substantial broadening of the injected microwave beam up to a factor of 2 was found in the scenarios considered. For all cases, WKBeam yielded larger broadening than the full-wave simulations. Up to fluctuation levels of approximately 50 %, however, the overestimation remains below 10 %. If the background density in the fluctuation layer exceeds values of 30 − 40 % of the cut-off density of the corresponding mode, the overestimation of WKBeam becomes more pronounced, reaching values of 20 % at about 70 % fluctuation level.
The relative deviation from the full-wave solution was found to depend only weakly on the propagation length through the fluctuation layer.
An important observation is the change of the transverse profile of the scattered beam from a Gaussian to a Cauchy distribution for strong scattering. One consequence are the elevated tails of the profile which means more power is scattered into directions far off the original propagation direction.
The parameter range investigated in this paper also includes the ITER scenario recently analyzed with WKBeam: values of δn e /n e = 20 % at a normalized background density of X = 0.2 were assumed [28]. According to the results presented in this paper, only a small overestimation on the percentage level is expected from WKBeam for these parameters, strengthening the main result from Ref. [28] that significant beam broadening should be expected for ITER (NTM stabilization should nevertheless still be achievable within the capabilities of the EC upper launcher system). One can furthermore conclude that the interaction of microwaves in the EC range of frequencies with edge density fluctuations can be well described within the limit of the Born approximation in large-scale fusion-relevant tokamak experiments.
No significant difference was found when changing the beam polarization from Oto X-mode. Since both codes are able to investigate X-mode polarized beams, a project to study cross-polarization scattering due to density fluctuations was started [57]. A thorough benchmark and analysis of this problem is however beyond the scope of this paper and will be published in a following paper.
in the antenna plane are input parameters. According to e.g. Ref. [1], w and R in the antenna plane are given by the following equations: with x R = πw 2 0 /λ 0 the Rayleigh range. Doing some algebra yields the required expressions for the full-wave code antenna input: x = w 4 π 2 R w 4 π 2 + λ 2 R 2 , . (A.4) Using the WKBeam input parameters w = 1.5 cm and R = 10 cm, values of w 0 ≈ 9.7 mm and x ≈ 58.2 mm are then obtained for the full-wave input parameters.