Parametric dependence of microwave beam broadening by plasma density turbulence

High-power microwave beams used for heating and current drive in magnetically confined fusion plasmas can be broadened significantly by plasma turbulence, negatively impacting the efficiency of the machine. The dependence of this beam broadening on plasma and beam parameters is not yet fully understood, particularly where the dependence on one parameter is not separable from the dependence on the other parameters, meaning the dependence must be expressed via functions of linear combinations of parameters, rather than functions of single parameters. The aim of this work is to develop an empirical model for how the broadening depends on plasma and beam parameters, allowing for the easy estimation of beam broadening by turbulence without the need for computationally expensive full-wave simulations. In this paper, a microwave beam is simulated propagating through a turbulent layer of plasma using the 2D full-wave cold plasma code EMIT-2D. The dependence of beam broadening on background plasma density, fluctuation amplitude, turbulence correlation lengths in the radial and poloidal direction, thickness of the turbulence layer, and microwave beam waist are considered. The parameter scans are conducted in pairwise combinations of the parameters in order to determine the separability of the dependencies. We find that the dependence on the radial and poloidal correlation lengths are not separable from each other, and neither are the dependences on the fluctuation level and the background density, but all other dependencies are separable. Ignoring this inseparability in the correlation lengths will usually result in an over-prediction of the broadening in tokamak plasmas. An empirical formula for the beam broadening based on the turbulence and beam parameters is found for fusion-relevant scenarios, making prediction of the effect possible in microseconds, instead of the hours required for full-wave simulation. This could then be of use for integrated modelling of heating and current drive systems.


Introduction
In magnetically confined fusion devices, microwaves are often used in diagnostic tools and to inject power into the plasmas for the purposes of heating [1,2], non-inductive start-up [3,4], and current drive [5]. In conventional tokamaks, global heating and current drive can be achieved via a range of mechanisms, including ion cyclotron resonance heating, neutral beam injection, and electron cyclotron resonance heating. Localized electron cyclotron current drive (ECCD) can also be used to stabilize MHD instabilities such as neoclassical tearing modes [6]. Microwave techniques have the advantages of requiring less space on the vessel wall due to a small launching antenna which can be placed far away from the plasma without loss of efficiency, and the gyrotrons used to generate the microwaves can be housed outside the bio-shield.
In order to inject power, microwave beams must cross the plasma boundary, a region where density fluctuation levels can reach 100% on length scales comparable to that of the microwave wavelength [7]. These fluctuations cause scattering of incident microwaves, leading to the broadening of microwave beams travelling through the plasma. For ITERlike scenarios, it has been predicted that this broadening could result in the doubling of the beam width compared to if no turbulence was present [8]. This broadening has the potential to significantly impact the efficiency of both global and local power injection via microwave beams. It's therefore important to be able to predict this effect when designing and running microwave power injection systems. Whilst full-wave simulations can achieve this, their computational expense makes them impractical when a wide range of plasma and launch scenarios are considered. The determination of the parametric dependence of this broadening could allow the effect to be accounted for in integrated modelling and optimisation studies for future tokamaks in a way that is computationally cheap.
There has been previous investigation of this effect, including development of analytical models for both O-mode [9] and X-mode [10]. This approach has the benefit of allowing for very fast prediction of beam broadening based on plasma and beam parameters, making it ideal for optimisation studies or inter-shot analysis. However, these models rely on methods which will not be applicable in the case of strong turbulence locally creating plasma densities close to or even above the cut-off density in the wave's path, meaning that they will not always be applicable for some fusion relevant scenarios.
Beam broadening has also been investigated using ray tracing methods where the scattering effect is described using a Fokker-Planck solver [8,11,12]. There have also been studies using a beam tracing code based on the wave kinetic equation, incorporating the effect of turbulence via a scattering operator derived under the Born approximation [13]. However, these methods are also valid only within certain limits, such as when the turbulence amplitude is small compared to the cut-off density, or when the refractive index does not vary over length scales that are small compared to the wavelength.
In order to simulate particular fusion relevant scenarios where the turbulent scale length is comparable to the microwave wavelength and the fluctuation level is high, directly simulating the beam path through turbulent plasma using full-wave codes is often required. Using this method, 1D parameter scans have previously been carried out [14]. However, the range of parameters is yet to be extended to encompass fusion relevant scenarios. Furthermore, it has not been previously investigated if the impact of changing one parameter is independent of changes in other parameters i.e. if the dependencies are 'separable' as defined in section 2.6.
The work presented here investigates how the broadening of microwave beams by a layer of turbulent plasma depends on plasma and beam parameters. The parameter ranges were set to cover fusion relevant scenarios such as beams used for ECCD, including ranges where turbulence scale length is comparable to microwave wavelength, and fluctuation level is as high as 50%. We considered pairwise combinations of parameters to determine whether the dependence on each parameter is independent of the others and conducted a point-wise fit to the data-set in order to determine an empirical formula for the beam broadening. By determining the dependence of broadening on turbulence and beam parameters in fusion relevant scenarios, a predictive model can be developed which does not require full-wave simulations but is still applicable in parameter regimes that are not analytically tractable. This would be of great use in development of future tokamaks by allowing the quick prediction of beam broadening, and could be incorporated into an integrated model used to optimise heating and current drive efficiency. It could also be used to make predictions on timescales useful for inter-shot analysis as well as real-time predictions based on plasma measurements.

EMIT-2D
EMIT-2D is a 2D full-wave, cold plasma code [15,16]. It uses a finite difference time domain (FDTD) scheme [17] to solve Maxwell's equations (1) and (2) with a cold plasma response equation (3) on a 2D Cartesian grid with absorbing boundary conditions: where c is the speed of light, ϵ 0 is vacuum permittivity, e is the elementary charge, m e is the mass of an electron, and the electron plasma frequency, ω p,e = (ne 2 )/(ϵ 0 m e ), where n is the electron number density. In equation (3) B 0 is the background magnetic field, and it is assumed that there is no background electric field or current density. The equation is then linearised in terms of perturbed quantities. It is also assumed that there are no finite temperature effects, and neglects the contribution of the ions due to their slow response time. The code is therefore valid only in cold plasmas where wave frequencies are sufficiently high, and background electric field and current density are sufficiently small. This is a remarkably good approximation for fusion reactor conditions unless we are close to a resonance.

Launched beam
The code launches a 2D Gaussian beam. The z-direction is the direction of propagation and the background magnetic field points in the y-direction. In order to excite the beam, a softsource antenna is used where the field added is given by: where w 0 is the beam waist radius at the focal point, z is the distance along the beam path, and r is the distance from the beam centre. The beam waist radius along the beam line is given by: and is defined as the distance from the beam centre at which the electric field amplitude falls to 1/e of its peak value. R(z) is the radius of curvature of the wave fronts and ψ(z) is the Guoy phase. By setting B 0 to point in the y-direction, we launched an Xmode beam in vacuum. We kept the normalized field strength constant throughout the simulation domain at Y = ω ce /ω = 0.5 where ω ce = eB/m e is the electron cyclotron frequency.

The turbulent density profiles
We chose to study the influence of the following key parameters: background density, fluctuation level, turbulence structure size, width of the turbulence layer, and beam waist. The fluctuations were only in the density profiles, making the turbulence purely electrostatic.
To launch the beam in vacuum, we transitioned the density from zero to its constant background value via a hyperbolic tangent function. In the constant background density region, we transitioned the fluctuation level from zero to a constant value via the same hyperbolic tangent function, before transitioning it back to zero. Though this sharp increase in density is not necessarily realistic for fusion plasmas, we decided that introducing a pedestal like background density would increase the dimensionality of the problem too much. As such, this simplified approach was taken. By considering the simulations run through background density profiles, we were able to determine that this sharp increase in density did not perturb the beam. An example of the background density profile and fluctuation envelope, with an image of a turbulent snapshot are shown in figure 1.
The simulation domain shown in figure 1 is in 2D. This is an acceptable simplification as the turbulent structures are elongated along field lines [7], making the scattering in the toroidal direction small compared to the radial and poloidal direction. The simulation domain is therefore approximately equivalent to a poloidal cross-section in a tokamak.
The turbulent density profile shown in figure 1 represents a snapshot of the plasma. In reality, the plasma would be moving at a speed significantly slower than the wave speed. This allows us to make the assumption that the turbulence is 'frozen', and the overall effect on the beam can be calculated by averaging over an ensemble of turbulence profiles.
In order to generate an ensemble of turbulent density profiles, we used synthetic turbulence. This allowed full control over the turbulence parameters and, for the large number of profiles needed, was less computationally demanding than fluid turbulence generation from a turbulence code.
The density fluctuations were generated using a truncated sum of Fourier-like modes given by [18]: where A ij are the amplitudes of the modes, and ϕ i,j are independent random phases uniformly distributed between 0 and 2π. The amplitudes are related to the structure size by [18]: where σ x = π L bin /a x and σ z = π L perp /a z , a x and a z are the box size, and L bin and L perp are the turbulence correlation lengths in the x and z directions respectively. Within the geometry of the simulations, the correlation length perpendicular to the beam path is L perp , and the correlation length bi-normal to the beam path and the background magnetic field is L bin . In a tokamak scenario, these correspond approximately to the radial and poloidal correlation lengths respectively for perpendicular injection along the radial coordinate.
We generated an ensemble of twenty turbulence profiles for each combination of turbulence correlation lengths. We then scaled the turbulence to the required fluctuation amplitude, applied an envelope function to achieve the required turbulence layer thickness, and added it to the background density. Any negative density values were truncated to zero, resulting in the negative fluctuation amplitude being slightly lower than the positive fluctuation amplitude.

Calculating beam broadening
In order to determine the broadening, Γ, for each combination of parameters in a scan, we found the ensemble average RMS electric field (E RMS ) at the back-plane of the simulation for each turbulent snapshot. We then calculated the ensemble average of these E RMS profiles for the ensemble of turbulence snapshots, and compared it to a beam which propagated through an equivalent background profile with no turbulence present. For all simulations, we used an ensemble   Parameters that are varied, how they are varied, and the 'base' value each parameter takes in scans where it is not varied. Note that the background density is only varied through nine values instead of ten, as at the timescales needed for higher densities, the stability condition proved prohibitively computationally expensive.

Parameter
Varied as Min Max Step Base

Parameter scans
The parameters we investigated are listed in table 1, with the value they take in the 'base case' included in every parameter scan, and the ranges they are varied over. All parameters are normalised. Length scales are normalised to the vacuum wavelength, λ 0 , background density is normalised to the O-mode cut-off density n crit , and fluctuation amplitude is normalised to background density. We varied the turbulence correlation lengths logarithmically, as we anticipated that broadening will be much smaller when the ratio between correlation length and wavelength becomes either small or large, as seen in previous work [14]. We therefore wanted to treat the ratio symmetrically, which is done by varying its logarithm. From previous work in the field, we expected scattering to be maximal when the correlation lengths are close to the vacuum wavelength, and to decrease as correlation length moves from this value in either direction [14].
For low density gradients we expected scattering to increase quadratically with fluctuation amplitude [10,12,14]. We also expected it to increase with background density and width of the turbulence layer [10,14]. For the microwave beam waist, we expected an initial increase followed by a gradual decrease, as seen in previous studies [14]. This is due to the dependence of the microwave beam width on its initial waist varying similarly non-monotonically, as can be seen from equation (5).
In order to determine the separability of the dependence on each of these parameters, we carried out a series of 2D scans considering each of the pairwise combinations. Each of these 2D scans required 100 ensembles of 20 simulations to be run, each at an approximate computational cost of 80 CPU hours.

Fitting
To determine the dependence of the broadening, Γ = FWHM ensemble /FWHM background , on each parameter and confirm the separability of the dependencies, we performed pointwise fits to the data-set from each 2D parameter scan.
By separable dependence we mean that, for two parameters a and b, the overall broadening can be expressed as a product of two independent functions, f and g, such that: Within the context of a point-wise fit, the broadening for a particular combination of parameters can then be calculated as: where C is a normalisation factor, fixing the value of the 'base case' broadening to be Γ basecase = 1 + C. To allow for the possibility that the dependence is not separable, we included a rotation of co-ordinate axis in the fit, adding an additional fit parameter, θ. The fit parameters f i and g j then depend on (a cos θ − b sin θ) and (a sin θ + b cos θ) as: As each parameter was varied across ten values, for each 2D parameter scan this resulted in 20 fitting parameters for 100 data points. The exception is scans where background density is varied where, due to computational constraints, only nine values were included. This is owing to the fact that the size of the simulation domain required increases for increased broadening, and that the wave speed in higher density plasma decreases, making the simulations take significantly longer to run. It would be possible to expand upon the range included here, but for this work they were deemed computationally too expensive.
Once a θ value was found for each 2D scan, we knew which dependencies were separable and which were not. We then performed a point-wise fit to the whole data-set, only allowing for inseparability where it had been found. Dependence of broadening on microwave beam waist (w 0 ) and width of the turbulence layer(W turb ). Colour represents broadening factor compared to a beam propagating through the same background density profile, with no turbulence present.

Parametric dependence
Performing the fits to the 2D scans individually, the only pairwise combinations where we found the dependence to be inseparable were those for L perp /λ 0 vs L bin /λ 0 , and for δn e /n e,0 vs n e,0 /n crit .
An example plot of a 2D dependence of Γ on a pair of separable parameters (w 0 /λ 0 vs W turb /λ 0 ) is shown in figure 3. Plots of the 2D dependence of Γ on both pairs of inseparable parameters (δn e /n e,0 vs n e,0 /n crit and L perp /λ 0 vs L bin /λ 0 ) are shown in figures 4 and 5 respectively.
The results of a global pointwise fit, as described in section 2.6, are shown in figure 6. These show how the broadening depends on each parameter (or linear combination of parameters in the cases was the dependence is not separable). It should be noted that for the linear combinations of parameters, the extremal points only correspond to one ensemble of simulations. Every other point corresponds to multiple ensembles. This is why the errors on the extremal points are noticeably larger.
The uncertainty in the broadening measurements from the simulations comes from the standard error on the ensemble average E RMS . There is also an uncertainty associated with the standard error on the ensemble average correlation lengths, fluctuation level, and background density. This contributes to the error on the fit parameters, given by their covariance.
From figure 6, we can see that there is only a weak dependence on the sum of the logarithms of the correlation lengths and a stronger dependence on their difference. This corresponds to having a stronger dependence on the ratio (L perp /λ 0 ) 0.4 /(L bin /λ 0 ) 0.9 , meaning that having a larger L perp compared to L bin increases broadening. In fusion relevant plasmas, the poloidal correlation length (which corresponds to L bin ) is usually longer than the radial correlation length (which corresponds to L perp ). This would result in less broadening. Dependence of broadening on background density (n e,0 /n crit ) and fluctuation amplitude (δne/n e,0 ). Colour represents broadening factor compared to a beam propagating through the same background density profile, with no turbulence present. From this we see that broadening increases with both background density and fluctuation amplitude, and the tilted nature indicates the dependence is not separable. Figure 5. Dependence of broadening on turbulence correlation length in the direction perpendicular to the magnetic field (Lperp) and in the direction binormal to both the magnetic field and direction of beam propagation (L bin ). Colour represents broadening factor compared to a beam propagating through the same background density profile, with no turbulence present. From this we see that both correlation lengths affect broadening differently, and the tilted nature indicates the dependence is not separable. This rough trend matches the analytical result that (in the low fluctuation level limit) broadening is proportional to L perp /L 2 bin [10].
We can also see that the linear combination of n e,0 /n crit and δn e /n e,0 based on their difference has a weaker effect on the broadening than the linear combination based on their sum. This suggests tha perhaps a key consideration is the peak Figure 6. The normalised dependence of broadening on each parameter considered. For the case where the dependence of two parameters is not separable, the dependence on a linear combination of those parameters is shown instead. The orange points correspond to the base value of each parameter, which was used in scans where other parameters were varied. The fit parameters are normalised to this point. The solid lines serve as a guide to the eye. density of n e,0 /n crit + δn e /n e,0 and how close it is to a cut-off. Due to the parameter ranges chosen, most points of the scan are outside of the limit where analytical theory can be applied, where we find the dependence of n e,0 /n crit and δn e /n e,0 is no longer separable and the broadening no longer scales with (δn e /n e,0 ) 2 as predicted by theory [10].
The form of the dependence on microwave beam waist matches that found in [14], where, as w 0 /λ 0 is increased, there is first a sharp increase in the broadening, followed by a gradual decrease. We believe this is due to the dependence of beam width on w 0 /λ 0 even in the case where turbulence is not present, which can be seen in equation (5).
Finally, the dependence of beam broadening on width of turbulence layer shows that as W turb /λ 0 increases, so does broadening, until it reaches a maximal value and appears to level off. We believe this is due to a kind of saturation effect, where the beam can no longer become any broader without simply appearing as background noise.
The functional form of the dependence of broadening on plasma and beam parameters, based on separability of parameters, is then: The numerical values of these fit parameters can be found in appendix. The model presented in this paper is also available as a Python script [19].

Conclusions and further work
We used a 2D full-wave code to simulate microwave beams propagating through turbulent plasma density profiles to determine the dependence of the beam broadening on plasma and beam parameters. We carried out a series of 2D scans in order to identify which dependencies were separable. We found that the dependencies on L perp /λ 0 and L bin /λ 0 were not separable from each other, and neither were the dependencies on n e,0 /n crit and δn e /n e,0 . We found that all other dependencies were separable. We then found the dependence of each parameter, or linear combination of inseparable parameters, by means of a point-wise fit to the whole data-set.
Where applicable, agreement with previous studies was found. However, it would be beneficial to compare this to other methods such as the analytical models [9,10], ray tracing simulations [8,11,12], and beam tracing simulations [13,18] in order to determine the regions of agreement and disagreement.
The inseparability of the dependence of broadening on the two turbulence correlation lengths is an important result. Currently when simulating beam broadening in fusion relevant scenarios it is common only to have good data for the correlation length in one direction, resulting in both correlation lengths being set equal to each other. This has the potential to significantly under or over predict the broadening of the microwave beam. The sensitivity of broadening to correlation length emphasises the need for diagnostic tools which can measure turbulence parameters to high degrees of precision. For example, from figure 5 the effect of uncertainty in correlation length measurement can be seen. If we take an example where only L bin = 2λ 0 is measured, and for simulation purposes the correlation length are set equal to each other despite L perp = 1λ 0 , the broadening would be over-predicted by a factor of 1.2 resulting in a large error on the predicted deposition profile.
A detailed understanding of the parametric dependence of beam broadening on plasma and beam parameters also introduces the possibility of being able to use beam broadening as a turbulence diagnostic. For example, if the other parameters are known well enough, the ratio of the turbulence correlation lengths could be deduced from the measured beam broadening.
Using the parameter dependencies provided here, it should now be possible to predict how much a beam will be broadened by without further simulation, given the scenario falls within the parameter ranges investigated here. This can be done almost instantaneously, rather than the hours that would be required for an ensemble of full-wave simulations to be carried out. This ability to rapidly predict beam broadening could allow for its inclusion in integrated modelling and calculation during experiment. Table A1. Look up table for beam broadening by plasma turbulence based on plasma and beam parameters. p is the parameter (or combination of parameters where the dependence is not separable), f is the value of the fit function, and e is the error on that fit. p 1 = ( Lperp/λ 0 ) 0.4 / ( L bin /λ 0 ) 0.9 , p 2 = ( Lperp/λ 0 ) 0.9 ( L bin /λ 0 ) 0.4 , p 3 = 0.2 ( n e,0 /n crit ) − 1.0 ( δne/n e,0 ) , p 4 = 1.0 ( n e,0 /n crit ) + 0.2 ( δne/n e,0 ) , p 5 = w 0 /λ 0 , and p 6 = W turb /λ 0 .