Wigner-function-based solution schemes for electromagnetic wave beams in fluctuating media

Electromagnetic waves are described by Maxwell’s equations together with the constitutive equation of the considered medium. The latter equation in general may introduce complicated operators. As an example, for electron cyclotron (EC) waves in a hot plasma, an integral operator is present. Moreover, the wavelength and computational domain may differ by orders of magnitude making a direct numerical solution unfeasible, with the available numerical techniques. On the other hand, given the scale separation between the free-space wavelength λ0\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\lambda _0$$\end{document} and the scale L of the medium inhomogeneity, an asymptotic solution for a wave beam can be constructed in the limit κ=2πL/λ0→∞\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\kappa = 2\pi L / \lambda _0 \rightarrow \infty$$\end{document}, which is referred to as the semiclassical limit. One example is the paraxial Wentzel-Kramer-Brillouin (pWKB) approximation. However, the semiclassical limit of the wave field may be inaccurate when random short-scale fluctuations of the medium are present. A phase-space description based on the statistically averaged Wigner function may solve this problem. The Wigner function in the semiclassical limit is determined by the wave kinetic equation (WKE), derived from Maxwell’s equations. We present a paraxial expansion of the Wigner function around the central ray and derive a set of ordinary differential equations (phase-space beam-tracing equations) for the Gaussian beam width along the central ray trajectory.


Introduction: a model wave equation
In this work, we analyze electromagnetic waves of high frequency (order of 100GHz , with the associated vacuum wavelength of the order of few mm) propagating through a fluctuating plasma. Such waves are applied in nuclear fusion devices like tokamaks for heating and current drive purposes, and they are referred to as electron cyclotron (EC) waves [1]. For this kind of waves in tokamaks usually the wavelength is much smaller than the scale length L of the background medium (typically tens of centimeters). In such limit, semiclassical methods apply. We start from a precise definition of the medium scale L in terms of the gradient of the background medium, that is, for every spatial coordinate r and for every function f describing the background medium (for a plasma, f can represent the particle density, temperature, magnetic field). Neglecting for a moment the fluctuations we consider a stationary medium, that is, f depends on position, but not on time. For simplicity, we absorb the scale length L in a normalized and dimensionless spatial coordinate with the consequence that gradients of any quantity f with respect to x are of the same order as f itself. The semiclassical limit then is defined by where 0 is the free-space wavelength (in physical units).
Generally, the properties of a non-uniform dispersive medium enter the wave equation via the dielectric operator (1) for the electric displacement ⃗ D . The operator ̂ is composed by a Hermitian contribution ̂h (describing the dispersive properties of the medium) and an anti-Hermitian contribution îa (describing dissipation), In the case of EC waves in fusion plasmas, it usually makes sense to apply the cold plasma approximation to the Hermitian part ̂h [1]. In this approximation, the thermal motion of plasma particles is neglected, and the Hermitian part ̂h of the operator ̂ reduces to a Hermitian matrix. The focus of this paper is on the effect of random fluctuations of the medium on the propagation of the beam. Therefore, we consider a simplified model, and dissipative effects are not addressed, i.e., ̂a = 0.
It should be noted in this respect that the influence of the anti-Hermitian (dissipative) part of ̂ is usually limited to narrow regions in the plasma where the resonance condition between the wave and the (cyclotron) motion of the electrons is satisfied, and significant damping can occur. Outside these regions, it is justified (and a standard approach) to treat the propagation in the dissipationless limit.
We also consider an isotropic medium so that ̂h = ( , x)1 , where 1 is the identity matrix, and = ( , x) is the frequency-dependent real valued dielectric function. The dependence on is implied in the remaining part of this paper.
Electromagnetic waves are vector valued. For simplicity, we consider the case of a linearly polarized wave with constant polarization unit vector , i.e., the electric field is written as This is possible e.g., where is a symmetry direction ( ⋅ ∇ = 0 ), and we look for symmetric solutions ⋅ ∇E = 0 . The choices made simplify the mathematical model, retaining the essential physics effects that we are interested in (refraction, diffraction and scattering) without the technical difficulties. As an example, the assumptions hold for the O-mode of EC waves in slab geometry. Then the equation for electromagnetic waves in frequency domain reduces to the Helmholtz equation, for the scalar wave field E ≡ E(x) . The dielectric function can be expressed as and the function n = n(x) is referred to as refractive index and describes the local properties of the medium at a given frequency and spatial position x. In the case of EC waves in a plasma, the function n = n(x) can be computed from the plasma parameters [2].
In a turbulent fusion plasma, the correlation time C of fluctuations is typically in the order of tens of microseconds [3], whereas the beam is turned on for a time of hundreds of milliseconds or larger [4], and both time scales are orders of magnitude larger than the time of propagation of the beam (estimated by the typical size of the machine divided by speed of light). Heuristically, such time-scale separation can be invoked in order to justify a statistical approach to the description of the effects of random fluctuations on the wave field: Instead of physical time-dependent fluctuations of plasma parameters f(t, x), we consider time-independent random fields f ( , x) where varies in an abstract probability space with probability density ℙ . Then, the refractive index becomes a random field n(x) = n( , x) (the variable is implied when not explicitly needed). For each point , one can solve the frequency-domain problem (7) with n = n( , x) . The solution E = E( , x) thus obtained is a random electric field, with the same probability density ℙ . The relation to physical quantities is provided by the ergodic hypothesis. E.g., for the wave field, that is, the time average of the physical wave field E(t, x) can be computed as an ensemble-average of the random wave field E( , x) . The same hypothesis is made for any derived physical quantity, such as the electric-field energy density |E| 2 ∕8 .
We always assume that the fluctuations are weak in the sense that with a random field n 2 of order O(1) and n 2 (x) = 0 . Further, we have the two-point correlation function This function enters the subsequent analysis and is proportional to the two-point correlation function of the refractive index perturbation n 2 − n 2 0 which is of order O −1 . Finally, we will be looking for the ensemble averaged Wigner function associated to the electric field which is a function of position x and refractive-index vector N, together forming coordinates in the phase-space z = (x, N) . Unlike the refractive index defined as a function of the spatial position x in (8), here N is an independent coordinate. A connection between both is drawn in Sect. 2.
Special attention must be paid to the fact that for wave beams (without dissipation) the wave field is not in L 2 (squared integrable). This is in contrast to the case of wave packets usually analyzed in quantum mechanics, e.g., in the work of Graefe et al. [5]. However, in literature a Wigner function description of beams can be found for various applications when time is a parameter in the wave equation. For example, Hermite-Laguerre-Gaussian Wigner beams are described in [6], and a decomposition into multiple Gaussian beams is employed in [7]. In the field of plasma physics, a robust description of drift waves and zonal flows as well as geodesic acoustic modes is possible upon a Wigner function description [8,9]. In the present paper, as an additional feature "mixing of the Wigner function" is induced by plasma fluctuations, as suggested in the work of McDonald [10].
In Sect. 2, we exploit the well-known fact that under simplifying assumptions of a straight beam trajectory the wave equation may be reduced to a Schrödinger-type equation. From this starting point, an initial value problem for the Wigner function can be formulated [10]. A paraxial solution is proposed in Sect. 3. In Sect. 4, mixing is quantified by means of an entropy functional [11]. A benchmark and examples are presented in Sect. 5 for standard test cases. Finally, a numerical scheme which allows for super-diffusive scattering is introduced in Sect. 6.

Paraxial wave equation
In order to illustrate the proposed phase-space paraxial expansion in a simple context free of technical difficulties, we consider equation (7) under the following conditions on the refractive index of the medium: sian coordinates, and z-axis is a symmetry direction; then n 2 (x, y) = n 2 0 (x, y) + 1 √ n 2 (x, y).
(ii) n 0 (x, 0) > 0 , so that no reflection from cut-off can occur; (iii) y n 0 (x, y) Under these conditions, we first simplify the wave equation (7) via the standard paraxial approximation. This preliminary step is however not essential. We consider a beam which is localized around the x -axis (to be referred to as the central ray) and propagating straight in x-direction. Existence of such solution requires the assumption made above in point (iii) that n 0 (x, 0)∕ y = 0.
For the wave field E, we introduce a change of variable with a phase factor (x) which will be chosen appropriately below [12][13][14] In addition, we assume that the new field a(x, y) is slowly varying in x-direction ( n a∕ x n = O(1) ) and rapidly decreasing in y-direction, i.e., y n a(x, y) = O −n∕2 and n a∕ y n = O +n∕2 . A physically relevant example of such a function is the Gaussian a(x, y) = a 0 (x)e − 2 y 2 w 2 (x) . As a consequence, from the wave equation (7) we find where the term proportional to 2 a∕ x 2 is of O −2 . Here, the refractive index function has been Taylor expanded, with the integral remainder q(x, y) . By assumption, the firstorder derivative vanishes.
As an asymptotic series in −1 , the leading order term in equation (14) is zero if we choose In equation (14), we have the first derivative of a with respect to x only and, thus, can consider the x-coordinate as a parameter of the evolution of a. However, we introduce the change of variable x = x( ) , defined by the transformation where we have assumed that n 2 0 (x, 0) > 0 . This absorbs the prefactor d d x in (14). Further, we define a rescaled field function At last, to the lowest order we are left with where, with some abuse of notation, we write q( , y) and n 2 ( , y) for q(x( ), y) and n 2 (x( ), y) , respectively. Equation (19) is a one-dimensional Schrödinger equation, where −1 plays the role of the Planck's constant, the first term on the r.h.s. represents the kinetic energy operator and the second term the potential.
The derivation of the WKE for this type of equations is one of the examples in the paper of McDonald [10], which can be applied to equation (19). In this specific case, the Hamiltonian is which is the Weyl symbol of the operator in equation (19). The dispersion relation H = 0 with H given in (20) has been introduced. The Wigner function W A associated to the amplitude field A is related to the averaged Wigner function W E introduced in equation (12) by which follows from the respective definitions.
One should notice that the substitution of (20) yields the dispersion relation which is the standard paraxial approximation of the Helm- The WKE (23) is the starting point for the derivation of a set of phase-space beam-tracing equations in the next section.

Paraxial solution scheme
We search for approximate solutions of (23) such that the cross section in phase-space, that is w at constant , is a Gaussian around y = N y = 0, with the amplitude c and the quadratic form in the exponential being the unknowns to be determined. The matrix G associated to the quadratic form is assumed to be strictly positive definite.
Initial conditions are given at = 0 . For EC beams usually an initially Gaussian shape for the wave field is assumed, with A 0 a constant measuring the initial field amplitude, with R 0 the initial wave-front curvature radius and Explicit computation shows that the Wigner transformation (12) in y-direction is given by (28) at = 0 and with This set of equations yields the initial parameters c(0), G rr (0) , G rN (0) and G NN (0) in terms of the initial beam parameters.
From ansatz (28) and the fact that G is strictly positive definite, it can be proven that y m N n y w = O − (m+n)∕2 , justifying a paraxial expansion of the WKE around y = N y = 0 in the semiclassical limit.
In substituting (28) into (23), we observe that The derivatives of the Hamiltonian are where H is given in equation (20). The fluctuation spectrum to lowest order reads where Γ 0 ( , N y ) = Γ( , 0, 0, N y ) , and the reminder for O N has been estimated by means of the relation (21) between N , y and N y . We consider the diffusive limit of the scattering operator [15,16]. Heuristically, this holds when the scale ΔN y of w in direction N y satisfies L C ≫ ΔN y [17,18] where L C is the correlation length defined as the width of the two-point correlation function C(x + s 2 , x − s 2 ) with respect to s. Then, Taylor expansion of the Wigner function w( , y, N � y ) in the (30b) scattering operator around N y up to second order is appropriate. Under the integral, zero-order terms cancel out, the firstorder term does not contribute due to symmetry of the Gaussian, and the resulting second order term yields Here, the diffusion coefficient has been defined, and the remainder arises from the fourth order term of the Taylor expansion. The second derivative of w in (28) amounts to Using (31-34) and (36) in equation (23) and collecting powers of (y, N y ) yields This system of ordinary differential equations describes the evolution of the beam parameters c, G rr , G rN and G NN and is referred to as phase-space beam tracing equations. One should notice that when D = 0 , the solution of (37) is independent of , but in general the diffusion term introduces a residual dependence on .

Entropy
In presence of fluctuations, the statistically averaged Wigner function does not necessarily belong to the range of the Wigner transform, since in general we cannot guarantee that there exists a deterministic wave field Ẽ such that (E(x)E(y)) =Ẽ(x)Ẽ(y) . When this is the case, e.g., when E =Ẽ with probability one, then we say that the wave field is in a pure state. In general, we speak of mixed states. Following [19], we use an entropy functional which measures how distant w is from a pure state: Indeed, it can be shown that we have S = 0 for the initial cross section with parameters (30). In this paper, we focus on the paraxial approximation, and with ansatz (28), the integrals can be computed analytically, with the result that From (39), the time derivative can be evaluatd upon making use of the phase-space beam tracing equations (37), with the result that This shows that with fluctuations suppressed (i.e., D = 0 ) the derivative of the entropy vanishes identically and, in consequence, the beam remains pure. By contrast, for D > 0 we have a monotonously increasing entropy with S → 1 for → +∞ , which shows the increasing impact of fluctuations.

Examples
Convergence of the paraxial solution is demonstrated for the cases of free space and a lens-like medium, addressed in Sects. 5.1 and 5.2, respectively. In both cases, we choose with F the (in general -dependent) fluctuation strength which corresponds to the root mean square of the refractive index fluctuations.
The test scenarios are chosen with a focus on simplicity of the resulting expressions. However, background media with more complicated refractive index functions n would also be possible as far as the points (i)-(iii) of Sect. 2 are fulfilled. For the application of EC beams in tokamaks, one could e.g., imagine a beam propagating through a layer of increasing refractive index n = n(x) , following the example of [20].
For the example results, the Wigner function itself is considered as well as the projections on the configuration space and refractive-index space, given by (39) , respectively. With the paraxial ansatz, one has that both E y and E N y are Gaussians, of widths respectively. As a reference solution, we have analytical expressions (no fluctuations only) and the WKBeam solution [17] (including the effect of fluctuations).

Free space results
Here, we make the idealized assumption that the beam propagates in free space ( n 2 0 = 1 and q = 0 in equation 15) but that fluctuations of the refractive index around n 2 0 = 1 are possible. This reduces the phase-space beam-tracing equations (37) to We first discuss a scenario without fluctuations, i.e., D = 0 , and a focused beam propagating in positive x-direction. Figure 1 shows the phase-space structure for several cross sections. Due to the Gaussian ansatz (28), elliptical contours are observed.
In general, during the propagation of the beam the ellipse rotates and gets squeezed. Via (30c) it is seen that the offdiagonal term in the exponential of the Gaussian ansatz vanishes if and only if the beam curvature is infinite, and this occurs at the focal point of a Gaussian beam. Accordingly the principal axes are aligned to the coordinate axes in that case (top right of Fig. 1). From equation (45a), it follows that the maximum value of the Wigner function at y = N y = 0 remains constant. Along with equations (45b) to (45d), it can be proven that the phase-space integral ∫ w d y d N y is a constant of motion, as it should due to energy conservation. Further, it can be checked that S( ) ≡ 0 . This shows that, indeed, the Wigner function under study corresponds to a pure state as expected if fluctuations are not considered.
The projection on configuration space (42a) is proportional to the energy density of the beam resolved in y and, thus, establishes a link to a wave field description. The result of equation (44a) along with (43a) is in agreement with the standard results on the width of a Gaussian beam in free space [21]. For the beam parameters as used in Fig. 1, the width is plotted in Fig. 2 (no fluctuations case). It shows a minimum at the position of the focal point.
In order to understand the effect of fluctuations on the beam, we consider a fluctuation layer with constant fluctuation strength of F(x) = 0.02 if 0 < x < 10cm and zero elsewhere and a two-point correlation length L C = 3cm.
The entropy as computed from formula (39) increases inside the fluctuation layer and then remains constant, cf. Figure 2. This confirms that mixing of the wave field is induced by fluctuations.
Due to the special structure of the diffusion operator (34), which modifies the shape of the cross section in N y -direction only, fluctuations do not directly induce beam broadening in configuration space (top panel of Fig. 2). By contrast, diffusion affects the refractive index projection directly (middle panel of the figure). Fluctuation-induced beam broadening in configuration space is observed after further propagation and must hence be understood as an indirect effect, provoked by the broadened spectrum.
In order to demonstrate the validity of the paraxial solution in Figs. 3 and 4, we show a comparison of the computed beam width to the analytical solution disregarding fluctuations and to the WKBeam solution.
In both cases, a clear effect of fluctuations is found. We see good agreement of the paraxial solution with the  = 120cm and x = 160cm ). Beam frequency f = 140GHz , initial beam width Δ 0 = 3cm , initial curvature radius R 0 = 100cm , no fluctuations WKBeam solution (not based on the diffusive approximation) in Fig. 3, whereas the paraxial approximation cannot reproduce the WKBeam results in Fig. 4. This can be explained by the small two-point correlation length in the latter case which invalidates the diffusive approximation, see also discussion in Sect. 3.
For the application of EC waves in tokamak plasmas, the fluctuation parameters are within the diffusive regime for an ITER scenario and beyond for an ASDEX Upgrade scenario [18].

Lens-like medium results
A background medium with and a central refractive index n 0 (x, 0) = 1 exhibits cut-off layers at y = ±1 . The results for such a profile of the refractive index are presented in Fig. 5.
The reference solution [22] (solution of the Helmholtz equation (7)

Multi-Gaussian solution scheme
In this section, we propose a heuristic approach to the approximation of the solution for w in terms of Gaussian functions beyond the diffusion limit (34). The difficulty is due to the We apply to equation (47) the Lie-Trotter splitting scheme [23]. Such solution method requires a time discretization k = k ⋅ Δ with a small time step Δ . Formally, the solution w( k ) can be computed from the solution at the previous time w( k−1 ) as  where Φ A+B is the time evolution operator. We consider the sub-problems solved by the evolution operators Φ A and Φ B , respectively. In the Lie-Trotter scheme, the time evolution operator of the complete problem is approximated by a composition of the separate operators, For the operator A defined in (48), the effect of Φ A on a beam mode has been derived in Sect. 3: The matrix G solves equations (37) with D = 0 , which then can be integrated along the time step.
Integration of equation (53) in terms of the explicit Euler method yields The factor in round brackets multiplying the function w leads to an amplitude decay of the propagating mode. The second summand on the r.h.s. amounts to the convolution of the Gaussian fluctuation spectrum Γ 0 with the Gaussian w. The result of such convolution is again a Gaussian function with parameters (c, G rr , G rN , G NN ) to be determined explicitly from the convolution term. Thus, operator Φ B applied to a Gaussian generates a new Gaussian which is summed to w( k−1 ) . The total amount of energy carried by the beam (measured by the y-N y -integral of the Wigner function) is conserved when both the amplitude decay and generation of a new mode are accounted for.
The number of Gaussian contributions is exponentially growing during time propagation, as illustrated in Fig. 6.
The scheme presented here can easily hit the wall due to insufficient computer resources, given the fast increase of the number of equations to be integrated. Even though a direct solution of the equations can be provided with a limited number of time steps rather than an out-of-the-box technique, it should be understood as a basis for more efficient numerical approaches which will be investigated in the future.

Summary
We have addressed paraxial solutions of the WKE thus reducing the problem of beam propagation to integration of ordinary differential equations. The results produced by such method have been found accurate within the limits of validity of the paraxial approximation for a beam propagating in free space and in a lens-like medium. The evolution of the shape of the beam cross section in phase-space has been discussed, and the entropy has been introduced. As expected, the paraxial approach is inaccurate whenever the diffusion approximation of the scattering operator breaks. For such cases, however, we have shown how the paraxial approximation can be combined with a Lie-Trotter splitting scheme which in principle can also deal with such scenarios. While a naive implementation of the splitting method is computationally too expensive due to an exponential increase of Gaussian terms in the solution, more efficient generalizations of such an idea are currently under investigation. Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.