Cross-polarization scattering of diffracting electron-cyclotron beams in a turbulent plasma with the WKBeam code

In turbulent plasmas, density fluctuations are expected to scatter radiofrequency wave beams, causing a degradation of the beam quality and thus reducing their performance. The WKBeam code is a Monte-Carlo solver for the wave kinetic equation, which describes such an effect, so far limited to a single wave mode so that it cannot account for cross-polarization scattering. In this work a new feature of the WKBeam code is presented, which allows the analysis of cross-polarization scattering of electron-cyclotron (EC) wave beams in realistic tokamak scenarios. We prove the convergence of the numerical scheme and give a preliminary overview of such effects in ITER scenarios.


Introduction
It has been suggested that in next-generation tokamak machines such as ITER, scattering by turbulent fluctuations can deteriorate the quality of electron-cyclotron (EC) beams and thus impair applications relying on high localization of the power deposition [1,2,3]. The WKBeam code [4] provides numerical simulations for the propagation of EC waves in a magnetically confined plasma, taking into account scattering phenomena due to density fluctuations. The code is based on a Monte-Carlo solution of the wave kinetic equation in the geometrical optics phase-space [5] for a single wave mode. Relevant quantities -e.g., the power deposition profile -can be analyzed in realistic tokamak scenarios for either the ordinary (O) or extraordinary (X) mode. The method can also treat diffracting beams and short-scale fluctuations, essentially being limited only by the applicability of the Born scattering approximation, which imposes an upper bound on the density spatial correlation function.
Under certain circumstances, however, plasma turbulence can induce mode-to-mode transitions, i.e., part of the injected power, initially assigned to the desired wave mode, can be transferred to another mode with different polarization. Cross-polarization scattering could be significant, in particular, in presence of large amplitude turbulent fluctuations in the tokamak edge, where the cold-plasma wave modes are nearly degenerate and the turbulence spectrum might induce a "jump" in wavevector space between the two branches. The aim of this work is to design and implement a numerical scheme for WKBeam in order to account for such cross-polarization scattering events. Specifically, we design a Monte-Carlo algorithm for the wave kinetic equation (WKE) in the constrained steady-state form which is relevant to time-harmonic coherent wave fields, namely, where the unknown is the statistically averaged Wigner function w α = w α (x, N ) representing the wave energy density carried by the mode α (α = X and α = O for the extra-ordinary and ordinary mode, respectively). More specifically, w α is a (generalized) function defined on the standard geometrical optics phase-space with coordinates (x, N ), where x = (x 1 , . . . , x d ) is the position in the physical space normalized to the scale length L of the plasma equilibrium density profile, and N = (N 1 , . . . , N d ) is the refractive index vector (N = ck/ω, where c is the speed of light in the phase-space and ω the angular frequency of the wave beam). Eq. (1) describes the transport of the Wigner function in the phase-space. Particularly, are the Poisson brackets of Hamiltonian ray theory, γ α is the absorption coefficient and, S α ({w β }) accounts for scattering due to density fluctuations. Specifically, where Σ α (x, N ) = β R d σ αβ (x, N, N )dN . The scattering cross-section σ αβ (x, N, N ) is related to the probability of the transition (α; x, N ) → (β; x, N ), i.e., a ray in the mode α at (x, N ) being scattered to a ray of mode β and refractive index N (α = β means scattering to the same mode); the position x never changes. Eq. (2) plays the role of a dispersion relation, while Eq. (3) gives a boundary condition on the surface Σ, which represents the launching mirror in the EC system. A statistical description of turbulence combined with tools from the calculus of pseudodifferential operators provides the theoretical framework for the derivation of the WKE from the relevant wave equation for the electric field of the beam, accounting for density fluctuations and weak absorption, in the short-wavelength limit [6,5].

Description of the algorithm
The code WKBeam [4] provides a numerical solution of the WKE in realistic tokamak scenarios, together with a set of post-processing tools which allow an analysis of the with initial condition u 0 α chosen so that w α = u α dt satisfies (1)-(3). The path followed is very similar to standard ray tracing, coupled with a random generation of scattering events. The beam is modeled as an ensemble of markers in the phase space which follow the rays and randomly modify their trajectory in presence of the fluctuations. The code as presented in [7] was able to solve the problem for a single wave mode (O-mode or X-mode), neglecting the terms σ αβ with α = β in (4). In this work a modification of the algorithm is presented, such that those terms are actually taken into account. With the new algorithm, the possibility of energy transfer to undesired polarization modes can be investigated. Following the theory of Monte-Carlo methods for transport processes [8], we design a numerical scheme which approximates a stochastic process {(X(t), N (t), α(t))}, constructed so that its probability density solves the Cauchy problem (5) with γ α = 0. The parameter t ∈ [0, T ] along each stochastic ray trajectory is discretized with constant step ∆t = T /n s , thus defining the intervals the left and right limits of the stochastic ray path at each node t = t j . Given a number N m of markers labeled by the index i ∈ {0, . . . , N m − 1}, the algorithm consists of the following operations: The number of scattering events in the interval is extracted from a Poisson distribution of parameterλ∆t, with (v) For each of such events, scattering from mode α to mode β occurs with probabilitỹ The target mode α i,j+1 is selected accordingly and the new value N + i,j+1 is determined through the Metropolis-Hastings algorithm from the probability densitỹ (vi) Iterate to the next interval I j+1 , with initial condition (x + i,j+1 , N + i,j+1 , α i,j+1 ). By standard techniques [8], one can prove that the probability density of such process approximates the density of {X(t), N (t), α(t)} which solves the auxiliary Cauchy problem (5), with γ α = 0. Absorption is then introduced by weighting each marker with the weight exp −2 t 0 γ α dt .

Verification and convergence
In order to test the algorithm and its convergence, we consider a stand-alone test problem, in which the particular form of the scattering operator allows the use of a pseudospectral solver to compute a reference solution. The simplified model consists in considering an equivalent Cauchy problem (5) in a 2D phase-space (one coordinate x in space, one coordinate y in momentum), with the model cross-section and Hamiltonians with associated Hamiltonian fields X α = (∂ y H α , −∂ x H α ) given by The numbers s αβ ∈ R + quantify the scattering from mode α to mode β, and L c plays the role of the turbulence correlation length normalized to the scale length L; the spectral width in (7) is ∝ 1/L c . The absorption term, which does not take parte to the definition of the relevant stochastic process, is omitted in this test-model. Under these assumptions, the integral operator S α amounts to a convolution, making the pseudospectral solver a suitable tool to build a reference deterministic solution. Fig. 1 shows a test run with an initial Gaussian distribution in the mode α = 0. Panels 1a and 1b show the pseudospectral solution, while panels 1c and 1d the Monte-Carlo solution reconstructed by means of a kernel density estimation (KDE). Each panel shows the reconstructed solution u α (t, x, y) at various values of the parameter t = 0, 1, . . . , 5. The contours of the pseudospectral solution show excellent qualitative agreement with those of the reconstructed Monte-Carlo distribution. The solution evolves from the initial position, set at (x, y) = (5.5, 3.4), towards the left-hand side and scatters in momentum y, into both the launched mode α = 0 (panels 1a/1b) and the secondary mode α = 1 (panels 1c/1d), which is not present at t = 0. The streamlines of the two Hamiltonian fields X α for the corresponding mode are represented (black arrows). One can notice a deflection of the center of mass of the solution from the streamline trajectory. This effect is entirely due to cross-polarization scattering: energy is transferred to the "passive" mode α = 1, flows along the streamlines of X 1 and has a finite probability of being scattered back into α = 0, effectively producing a mechanism to transfer energy perpendicularly to the Hamiltonian field. goes to zero as 1/ √ N m . Here the max is computed for K running over the cells of the dual grid of the pseudospectral solution, µ Nm is the measure that counts the number of markers in a certain set, and ν is the measure associated to the pseudospectral solution.

Preliminary physics results
We present here some examples of simulations performed with WKBeam in realistic tokamak scenarios, in order to provide a taste of what kind of situations could be investigated thanks to the cross-polarization scattering algorithm. A crucial role in this framework is played by the model used for the density fluctuations. In all cases the two-point density correlation function -which largely determines the cross section σ αβis taken in the form where n e (ρ) is the local equilibrium density, F (ρ, θ) is an envelope, ρ = √ ψ, with ψ being the normalized poloidal flux of the tokamak equilibrium, and θ the geometrical poloidal angle. The functions ρ and θ are evaluated at the middle point (x + x )/2. At last the correlation is Gaussian in x−x , with the matrix A accounting for the turbulence anisotropy along the equilibrium magnetic field direction. One has maximum flexibility in the choice of the envelope F (ρ, θ). Models used in this paper are sketched in Fig. 3. A possible choice for F (ρ, θ) is a Gaussian with center at ρ = ρ 0 , amplitude (δn e /n e ) 0 , standard deviation δρ, and uniform in θ. We observe different behaviours for different values of such parameters. In Table 1  In ASDEX the effect of fluctuations is actually negligible, as already pointed out in [7]. Fig. 4 (blue squares, red circles and green triangles) shows the power deposition profiles for ITER, together with the curve in absence of fluctuations for a comparison (black stars). All ITER scenarios refer to the upper launcher, low stearing mirror, with power deposition on the q = 2 resonant surface. The beam frequency is 170 GHz, and the correlation length of the fluctuations is 2cm in the direction perpendicular to the local magnetic field, while it is set to a large number in the parallel direction in order to take into account the flute-mode character of turbulence. For the scenario "ITER 1", cross-polarization scattering is negligible, as in the region where the fluctuations are located (ρ = 1.0) the density is high enough, and the dipersion relations of the modes are well separated [7]. Nonetheless, simple scattering to the same launched mode produces a significant broadening of the power deposition profile. In the scenario "ITER 2", the Gaussian layer of turbulence is moved far out in the scrape-off layer (ρ = 1.04), and its relative strength (δn e /n e ) 0 is increased by a factor 10. This is an extremely severe fluctuation profile, which is articially manufactured in order to obtain a scenario in which a considerable fraction of the beam is scattered into the secondary mode (X-mode), and is afterwards reflected at the cutoff, located near the plasma separatrix. This produces a loss in terms of power deposited in the target region. On the other hand, the width of the power deposition region is not degraded much, in contrast to "ITER 1" where simple scattering is dominant. Together with the two standard Gaussian models, we have constructed a third scenario, ITER 3, in which (δn e /n e ) 0 = 0.02 and the envelope is represented in Fig. 3a (green triangles curve). This should reflect better the expected shape of the correlation function at the edge of ITER plasmas. With such a model, about 6% of the injected power is scattered to the X-mode, which is reflected at the cut-off and does not reach the target resonance. Moreover, a significant broadening of the beam can also be observed (Fig. 4, ITER 3). A poloidal projection of the beam in this scenario is shown in Fig. 5. Figure 5: Scenario ITER 3 (poloidal section). The beam is injected from the upper launcher, lower stearing mirror, on the q = 2 resonance surface. About 6% of the power is scattered to the X-mode and reflected. A considerable broadening is observed for the O-mode near the deposition region. Again, the inaccuracy in the reconstruction of the X-point region is an equilibrium post-processing artifact which does not affect our simulations.