Paraxial beams in ﬂuctuating fusion plasmas: diffusive limit and beyond

. A paraxial expansion of the (ensemble-averaged) Wigner function in the relevant wave kinetic equation for electron cyclotron waves in ﬂuctuating plasmas allows the derivation of phase-space equations similar to the equations for the Gaussian beam parameters in the paraxial WKB method [G.V. Pereverzev, Phys. Plasmas 5 , 3529 (1998)]. This is relatively straightforward when the scattering of the wave ﬁeld by density ﬂuctuations can be described by a di ﬀ usion operator in refractive-index space. The general case is rather more complicated, yet we could ﬁnd a heuristic construction of a paraxial Wigner function. Here we use a simple model, which has an analytical solution, to test both the theoretical validity of the di ﬀ usion approximation and the heuristic paraxial approach beyond the di ﬀ usion approximation.


Introduction
The wave kinetic equation (WKE) provides a suitable tool to tackle the problem of the scattering of short-wavelength waves from density fluctuations in fusion plasmas [1,2]. In this approach, the typical correlation length of the fluctuations is not required to be larger than the radiation wavelength. Just the spatial variation of the correlation length and of the root-mean-square density fluctuation is required to satisfy the WKB limit, and this is usually the case, since their magnitude is related to the equilibrium profiles. The solution of the WKE gives the lowest-order approximation of the statistically averaged Wigner function of the wave field in the short-wavelength limit, where the statistical average is taken over an ensemble of density fluctuations. Density fluctuations, typically generated by plasma turbulence, are predicted to lead to a significant beam broadening particularly in large machines, in which the beam path across the plasma can be much longer than in present-day fusion experiments [3,4]. A direct numerical solution of the WKE for electron cyclotron (EC) waves in tokamak geometry has been implemented in the WKBeam code [5]. Comparisons with both fullwave codes [6] and experimental data [7] supported the applicability of the underlying scattering model (based on the Born approximation) to typical fusion-plasma parameters. The numerical approach implemented in WKBeam, although much more efficient than a direct solution of the wave problem (if this can be afforded at all), still leads to a much heavier computational burden as compared to beamtracing codes like TORBEAM [8] or GRAY [9], in which, however, the effect of density fluctuations is not included.
A possibility to reduce the computational cost of the solution of the WKE may still come from a generalization of the paraxial approach used in TORBEAM to the averaged Wigner function. In fact, from the results of the * e-mail: emanuele.poli@ipp.mpg.de WKBeam code, we observe a strong localization of the averaged Wigner function around a trajectory in phase space. In previous work [10], we showed how such a paraxial approximation of the Wigner function can be constructed in a simplified model. The result is a set of ordinary differential equations (ODEs), called phase-space beam-tracing equations. If the diffusion approximation for the wave scattering holds, the derivation of such phase-space beam-tracing equations is relatively straightforward (at least in this simplified model). In the general case, however, a different strategy for the evaluation of the scattering effects must be considered.
In this paper, we discuss the validity of the diffusive limit on the basis of an analytical solution of the WKE. We also use the analytical solution in order to assess the proposed paraxial approach in a particular scattering scenario, namely that of a thin fluctuation layer. In this context, "thin" means that the variation of the beam parameters due to refraction and diffraction can be neglected within the layer, while outside the layer the propagation occurs in the limit of vanishing fluctuations. This setup is similar to the phase-screen model employed in trans-ionospheric communication [11].

Model problem and analytical solution
The Wigner function description of EC wave beams is formulated in spatial coordinates x, normalized to the scale L of variation of the background medium. Here, the scale length L is defined with the equilibrium quantities, that are assumed independent of time, and therefore does not capture the scale of fluctuations, which can be arbitrary. Further, we have the refractive-index vector N = ck/ω, where c is the speed of light in free space, and ω the angular frequency of the wave field. The statistically averaged Wigner matrix of the wave electric field [1] is defined on the position-refractive-index phase space with coordinates (x, N) and depends on the parameter κ ωL/c. Since L is the scale length of the equilibrium, for EC waves in plasmas we have κ 1. In the limit κ → +∞, each eigenvalue W of the leading-order term of Wigner matrix satisfies the constrained transport equation [1,5] where {F, G} = ∇ N F · ∇ x G − ∇ x F · ∇ N G denotes the canonical Poisson bracket of two smooth functions F, G in the (x, N) phase space, H is the geometrical optics Hamiltonian, hence H = 0 gives the dispersion relation of the wave mode associated to W, γ is the absorption coefficient, and with η a given function depending on the mode polarization (the precise definition [5] is not needed here), δ is the Dirac's distribution, and where n e (x) = n e,0 (x)+δn e (x)/ √ κ is the total electron density split into the equilibrium part n e,0 , and small random fluctuations δn e / √ κ. The reasons for the scaling 1/ √ κ are discussed in earlier work [1,5]. The expectationvalue operator E is the average over an ensemble of fluctuations δn e . Equation (1) complemented by appropriate boundary condition at the launching mirror is the WKE solved by the WKBeam code.
The simplified test problem addressed in this paper is obtained by choosing a two-dimensional domain (d = 2) with coordinates x = (x, y) = (τ, y), N = (N x , N y ), and where F = F(τ) = E(δn 2 e ) 1/2 /n e,0 is the root-mean-square of density fluctuations relative to the equilibrium density, and L C is the correlation length, here assumed to be constant. The dependence on κ is understood in the rest of the paper.
Upon using δ(H) = (2N x,0 ) −1 δ(N x − N x,0 ), and the ap- is a distribution defined on the reduced (y, N y ) phase space and satisfies the linear Boltzmann equation with One should note that simplification (4) is somewhat unphysical: On the one hand, one assumes that the frequency is so large that propagation can be described with the vacuum Hamiltonian; on the other hand, fluctuations are retained. Therefore model (4) should be regarded as a computational test. Nonetheless, such simplified model captures the essential parts of the WKE: Indeed, equation (5) is essentially the same as the more realistic model studied by Chellaï et al. [7], the only difference being simpler coefficients. Particularly, we can obtain an analytical solution of (5) on the line of the derivation sketched in the appendix of the paper by Chellaï et al. [7]. This analytical solution is then written in Fourier space, that is, where andv is the Fourier transform of the initial condition, with A 0 being the initial electric-field amplitude, w 0 the initial Gaussian beam width, and R 0 is the initial phase-front radius of curvature. The energy densities in both real and refractiveindex space are given by respectively, and their calculation requires the function v in the (y, N y )-space. However, their moments can be evaluated directly fromv [7]. In particular, from the secondand fourth-order moments we can compute the width and the kurtosis [12] of both distributions. The skewness is identically zero for both distributions. Upon defining the quantities for any integer n ≥ 0, the results can be written as where w x and w N are the widths in real and refractive-index space, respectively, whereas kurt(ρ) denotes the kurtosis of the density ρ. We recall that any value in excess of 3 indicates that the distribution decays to zero more slowly than a Gaussian, and thus have "fat tails". In addition, and w 2 0 /L 2 R are the squared beam width and spectral width of a standard Gaussian beam, with 1/L 2 R = 1/R 2 0 + 4/(κ 2 w 4 0 ).

Diffusive limit and diffusive scattering
For propagation distances longer than the transport meanfree path [13, section 5.1 and references therein] the energy density E x in physical space can be described by a diffusion equation. This regime is referred to as diffusive limit of the linear Boltzmann equation (5). A validity criterion for the diffusive regime has been discussed in [4] and can be re-derived here on the basis of kurt(E x ). For instance, with a piecewise-constant profile F 2 (τ) = F 2 max > 0 for 0 ≤ τ ≤ τ F and F 2 (τ) = 0 otherwise, from (10) one finds that as τ → +∞, is the parameter introduced by Snicker et al. [4] and It is remarkable that the kurtosis of both E x and E N have the same asymptotic behavior. We deduce that both distribution observed at large τ, i.e., well after propagation through the layer of fluctuations, have negligible excess kurtosis if either λ 1, in agreement with the criterion established earlier [4], or ξ 1. We shall now show that the parameter ξ controls the regime in which the scattering operator (6) can be approximated by a diffusion operator acting in refractive-index space only. This should not be confused with the diffusive limit [13] in which one finds a diffusion operator in position space. It is also somewhat stronger than the condition of negligible kurtosis, since a diffusion operator preserves the Gaussian form of the initial condition. Formally this "diffusive scattering approximation" is obtained by a Taylor expansion of the solution, where S r (v) involves the fourth-order remainder term in Taylor's formula for v (odd-order terms vanish due to symmetry of the integrand). Neglecting S r in equation (5) leads to with refractive-index diffusion coefficient The validity of (14) can be established on noting that the Taylor's expansion of ϕ 0 for q 2 2κ 2 L 2 C gives, cf. (8), and the exact solution (7) of equation (5) reduces to the exact solution of the diffusive-scattering approximation (14). Therefore, we expect that (14) yields a good approximation of the solution of (5) if the initial conditionv 0 is concentrated in a region of its domain where q 2 2κ 2 L 2 C . For a Gaussian beam, the width in q of the initial condition (9) is given by 8L 2 R /w 2 0 and thus the condition amounts to ξ 1 as anticipated. In the special case of flat initial phase front 1/R 0 → 0, we have L R = κw 2 0 /2 and this reduces to ξ = w 2 0 /L 2 C 1. We note that in realistic application, typically ξ > 1 and the full scattering operator (6) should be used.

Paraxial solution of the WKE
The Fokker-Planck equation (14), which holds for ξ 1, admits exact solutions of the form of a Gaussian distribution where the amplitude c and the Gaussian parameters G rr , G rN and G NN at the launch position can be obtained in terms of electric field amplitude, beam width and focusing, cf. (30) of [14]. Substitution into (14) gives ordinary differential equations for the parameters of the Gaussian, that is, [14] dc dτ = −2κDG NN c, The general form (5) of the WKE is more complicated since the operator S 0 does not preserve the Gaussian form of the function v. However, equation (5) can be written as an abstract evolution equation of the form [14] dv dτ = Av + Bv, where Av = −N y ∂ y v and Bv = S 0 (v). A time-discretized solution in terms of the Lie-Trotter operator splitting method [15] is proposed [14,16]. If 0 = τ 0 < τ 1 < · · · < τ k < · · · is a sequence of points in τ, the exact solution at τ = τ k formally can be written as with the exact evolution operator Φ A+B (τ k , τ k−1 ) accounting for the effect of both operators A and B during the time interval ∆τ k = τ k − τ k−1 and v(τ k−1 ) the solution before this time step. After Fourier transform, Φ A+B can be constructed in the same way as the exact solution (7). Analogously, let Φ A and Φ B be the evolution operators for the two sub-problems respectively. A first-order approximation of Φ A+B (τ k , τ k−1 ) is then obtained by the Lie-Trotter formula where h τ = max k ∆τ k is the maximum step in τ.
For solutions of the form (16) the evolution operator Φ A is exactly constructed by solving equations (17) with D = 0 for the parameters of the Gaussian. For the evolution operator Φ B associated to the full scattering operator, we apply the explicit Euler scheme which gives When applied to a Gaussian, the result of the action of this operator is an amplitude decay of the original Gaussian function (first term) plus the convolution of the (Gaussian) fluctuation spectrum with the Gaussian initial Wigner function. The result of such computation, again, has the Gaussian shape as in (16) with the parameters (c, G rr , G rN , G NN ) to be determined by evaluation of the convolution integral. As a result of the procedure described above, starting from a Gaussian Wigner function after one time step we have the sum of two Gaussian functions. Owing to the linearity of equation (21) the strategy can be applied on a superposition of Gaussians as well, doubling the number of such Gaussian components of the beam after each time step. Thus, at time τ k the obtained paraxial approximation v Euler of the solution can be written as superposition with an exponentially growing number of scattered contributions m = 2 k (starting with one mode at τ 0 ). Here, each Gaussian mode v i is described by ansatz (16), with a set of parameters (c i , G i ). The development of these modes is illustrated in figure 1. However, the fast growth of the number of Gaussian modes might require a large numerical effort when the solution is constructed in practice and makes it impossible to provide an analytical interpretation of the results. On the other hand, many effects may be retained when a simplifying model is considered. Here, we introduce the model of a thin fluctuation layer, described by the piecewise-constant profile introduced above, with F 2 max = D/τ F . The overall effect of fluctuations for τ > τ F is determined by the integral of the fluctuation amplitude function τ F 0 F 2 (τ)dτ = D and, thus, in this model measured by the parameter D. We consider layers so thin that the variation of the parameters in each Gaussian component v i of (22) due to the operator Φ A can be neglected, hence Φ A (τ k , τ k−1 ) ≈ 1. Within the layer, this simplification allows us to replace the paraxial solution v Euler with an adapted simpler approximation given by a linear combinations of functions that are constructed by a repeated application of the convolution operator in (6). This construction of the functions is reminiscent of the Krylov subspace method. The first element is chosen to be the initial Gaussian launched at τ = 0 with parameters where Φ = 2/κw 2 0 computed from the initial beam width w 0 and the phase front is flat (1/R 0 → 0) for simplicity. Then the functions ψ α = ψ α (y, N y ) for α = 0, 1, . . . are defined recursively by with constant γ α determined by the normalization condition ψ α (0, 0) = 1 and with K being the convolution in N y by the Gaussian exp[−κ 2 L 2 C N 2 y /2] so that Bψ = (κ 2 L 2 C F 2 /4)K(ψ) − Σψ. Explicitly, (24) yields ψ α (y, N y ) = e −κG rN with parameters defined recursively by 4 EPJ Web of Conferences 277, 01003 (2023) https://doi.org/10.1051/epjconf/202327701003 EC21 The recurrence relation for G NN α , in particular, leads to the closed equation which can be proven by induction. The solution within the layer is then approximated by a linear combination of such Gaussian functions, where C α = C α (τ) are coefficients to be determined. Outside the layer, each Gaussian component is propagated with the exact evolution operator Φ A (τ, τ F ), i.e., by solving (17) with D = 0.
The subspace spanned by a finite number of ψ α is not closed under the action of the operator K, or equivalently B. Therefore a truncation is applied, that is, where the term containing ψ m+1 has been dropped. The error of this approximation depends on how fast the coefficients C α decay to zero for α → +∞. The substitution of (26) into the equation for Φ B in (20) together with the truncated form of Bv (m) leads to with the coefficients and σ −1 = 0 by definition (no contribution to the nonscattered Gaussian function). Equation (27) can be readily integrated with the result that When (25), (28), and the definition of D are inserted, again by induction one can show that The projection onto the configuration space of the Wigner function (26) can be computed from (30), with the result After inserting the amplitudes C α (τ F ) and letting m → +∞ the sum becomes a series which can be computed analytically, with the result which is equivalent to the initial beam cross section in configuration space E x (0, ·). The conservation of E x is a property of the exact evolution operator Φ B , since S 0 (v)dN y = 0, and we find that it is preserved by the approximated solution (26) in the limit m → +∞.
As for the spectrum E N , with the approximated solution v (m) and in the limit m → +∞, we obtain  where the parameters λ and ξ are defined in section 3, cf. equations (12) and (13) with L R = κ 2 w 2 0 /2 for an initially flat phase-front, and G NN α = Φ −1 /(1 + αξ). This result reflects the fact that the scattering operator Φ B has a direct impact on the spectrum of the beam only, leaving its spatial structure unaffected. In turn, application of Φ A leads to a spatial evolution of the beam which preserves its spectrum. The spectrum E N computed with v (m) , equation (33), becomes a superposition of Gaussian functions of different widths. With every application of the convolution part of the scattering operator, G NN α decreases leading to eventually broad tails in the spectrum arising from higher α modes. We can compute the kurtosis of the distribution (32). The even-order moments are N 2n y E N (τ F , N y )dN y = π κ Φ n (2n − 1)!! (2κ) n C 0 (0) × e −λ/2 ∞ α=0 (λ/2) α α! (1 + αξ) n .
The remaining sums can be traced back to the Taylor series of the function f (λ) = e λ/2 and evaluated analytically, with the result kurt E N (τ F ) = 3 + 3 1 This result is identical to the exact asymptotic limit obtained in equation (11). The evolution of the spatial beam width, spatial kurtosis, spectral beam width and spectral kurtosis is shown versus the spatial position x = τ along the beam axis for two scenarios in the figures 2 to 5. For scenario 1, the parameters represent a relatively thin fluctuation layer (x F = τ F = 1 cm), whereas scenario 2 includes a broader fluctuation layer with x F = τ F = 30 cm. As expected, in all figures multi-Gaussian solution and the analytical solution coincide very well when a thin fluctuation layer is considered. In figures 2 and 3 one can see that within the fluctuation layer the quantities describing the spatial distribution are exactly conserved by the multi-Gaussian solution. As discussed above, the effect of a disturbed spectrum (due to the action of Φ B ) on the spatial distribution arises from the evolution of the beam via the transport operator Φ A . This evolution is deferred in the multi-Gaussian solution for which the operator Φ A is applied only after the fluctuation layer. In turn, the transport operator leaves the spectral distribution unchanged as can be seen in figures 4 and 5. Broadening of the spectrum inside the fluctuation layer can be observed in figure 4. The kurtosis in figure 5 starts from the Gaussian value kurt = 3 and, after reaching a maximum, matches the exact value analytically if m → +∞, as shown above.