Quasioptical modeling of wave beams with and without mode conversion. III. Numerical simulations of mode-converting beams Free

Geometrical-optics (GO) ray tracing has been widely used to calculate the propagation and absorption of electron-cyclotron waves (ECWs) in inhomogeneous magnetized fusion plasmas in many contexts, including electron-cyclotron heating, electron-cyclotron current drive, and electron-cyclotron emission diagnostics.^1,2 Also, a lot of alternatives,^3–13 which take into account diffraction to describe the beam width of ECWs near the focal region, have been proposed as extensions of the GO approach and have contributed to the improvement of the wave-power deposition-profile simulations. These approaches can treat single-mode waves in sufficiently dense plasmas. However, waves propagating in low-density plasmas often contain two electromagnetic modes, which have close refraction indexes and thus can interact efficiently. Specifically, in fusion plasmas, the mode conversion between the O and X waves that form ECWs is caused by the magnetic-field shear in the peripheral region. One of the ways to properly model this process is to perform one-dimensional full-wave (1DFW) analysis, such as that carried out in Ref. 14. However, one-dimensional analysis becomes inaccurate already when the wave beams experience substantial bending. This can be remedied to some extent by applying “extended geometrical optics” (XGO),^15–19 which is a reduced theory that describes refraction and mode conversion on the same footing. However, XGO, as it is formulated in Refs. 15–19, still does not include diffraction, and so the problem remains open.

The present work continues a series of papers where we propose how to overcome this problem. Our results are quite general; in fact, they are not restricted even to plasma physics, let alone fusion applications. The comprehensive theory generalizing XGO to include diffraction was reported in Paper I of this series.²⁰ In Paper II,²¹ we presented a corresponding quasi-optical code PARADE (PAraxial RAy DEscription) in its reduced version that can simulate single-mode beams without mode conversion. Here, we present a more general version of PARADE and report the first quasi-optical simulations of mode-converting beams. We also demonstrate that PARADE can model splitting of two-mode beams.

Although our work was primarily motivated by the need to improve ECW modeling in fusion plasmas, PARADE can also be useful for studying related problems in optics and general relativity.^22,23 For this reason, we illustrate the PARADE's capabilities below in basic-physics examples. (Fusion-relevant applications of PARADE are discussed in Ref. 24). As seen from these examples, the numerical results produced by PARADE are in good agreement with the corresponding results of ray-tracing and 1DFW simulations to the extent that those are applicable and such comparison is meaningful. More generally, PARADE's simulations surpass ray tracing in that they resolve diffraction and mode conversion. PARADE's simulations can also surpass 1DFW results in that they capture the true three-dimensional structure of wave beams, while the 1DFW approach is restricted to straight rays.

Our paper is organized as follows. In Sec. II, we introduce the key equations derived in Paper I and also adjust them to numerical modeling. In Sec. III, we report simulation results for test problems. In Sec. IV, we summarize our main conclusions.

Here, we outline how the general theory developed in Paper I can be applied, with some adjustments, to describe mode-converting wave beams. The general idea is similar to that presented in Paper II for single-modewaves. We assume a general linear equation for electric field E of a wave,

where $D̂$ is a linear dispersion operator. We also assume that this field can be represented in the eikonal form,

where $ψ$ is a slow complex vector envelope and θ is a fast real “reference phase” to be prescribed. The wave is considered stationary, and so it has a constant frequency ω; then, $ψ$ and θ are the functions of the spatial coordinate x. Correspondingly, the envelope $ψ$ satisfies

where $D̂$ serves as the “envelope dispersion operator” (⁠$≐$ denotes definitions). We introduce $k ≐ ∇θ(x)$ for the local wave vector, $λ≐2π/k$ for the corresponding wavelength, $L∥$ for the inhomogeneity scale of $ψ$ along the beam, and $L⊥$ for the minimum scale of $ψ$ across the beam. The medium-inhomogeneity scale is assumed to be larger than or comparable with $L∥$⁠, and we adopt

Then, Eq. (3) becomes

where $D$ serves as the “effective dispersion tensor” found from $D̂$ and the operator $L̂ϵ=O(ϵ⊥)$ is specified in Paper I (also see below). We suppose the following ordering:

where the indices H and A denote the Hermitian and anti-Hermitian parts of $D̂$⁠, respectively. Assuming that the spatial dispersion is weak, $D$ can be replaced with the homogeneous-plasma dispersion tensor,

where $1$ is a unit matrix and $ε$ is the homogeneous-plasma dielectric tensor;²⁵ its dependence on ω is assumed but not emphasized since ω is constant. (Here, p denotes any given wave vector, as opposed to k, which is the specific wave vector determined by θ; see above.) Note that Eq. (7) assumes the Euclidean metric. Although we shall also use curvilinear coordinates below, expressing D in those coordinates will not be necessary as explained in Sec. VI D of Paper I.

With $D≈D$ and D given by Eq. (7), Eq. (6) implies $εA=O(ϵ∥)$⁠. Hence, $DH$ is the dominant part of D, and Eq. (5) yields

Since the Hermitian matrix $DH=O(1)$ is the dominant part of D and has enough eigenvectors $ηs$ to form a complete orthonormal basis, it is convenient to represent the envelope $ψ$ on this basis,

where a^s are the complex amplitudes. [Summation over repeating indices is assumed. For all functions derived from D, such as $ηs$⁠, the notation convention $f≡f(x)≡f(x,k(x))$ will also be assumed by default.] Then,

where Λ_s are the corresponding eigenvalues. Due to Eq. (8) and the mutual orthogonality of all $ηs$⁠, $Λsas$ is small for every given s individually; hence, either a^s is small or Λ_s is small. In the latter case, $ηs$ approximately satisfies the eigenmode equation, $DHηs=Λsηs≈0$⁠, and so it can be viewed as the local polarization vector of a GO mode. Then, the corresponding a^s can be understood as the local scalar amplitude of an actual GO mode, and $as=O(1)$ is allowed.

In this paper, we consider the situation when there are two modes, henceforth called O and X modes, whose eigenvalues $Λo≈Λx$ are both small within the region of interest. (This assumption is specified in Sec. II D. Also note that the single-mode case is discussed in Paper II.) In other words, the two modes are approximately in resonance with each other. Then,

where $ηo$ and $ηx$ are the O- and X-mode polarization vectors, $η¯$ is some third eigenvector of D that is orthogonal to both of them, and

The small amplitude $a¯$ can be easily calculated as a perturbation and is included in our theory²⁰ but does not need to be considered below explicitly. Instead, we introduce a two-dimensional amplitude vector

and the 3 × 2 “polarization matrix” $Ξ$ that contains the vectors $ηo$ and $ηx$ as its columns,

Then, $ψ$ can be expressed as follows:

We also introduce the dual basis vectors $ηo$ and $ηx$ and an auxiliary polarization matrix,

As seen easily, this matrix satisfies $Ξ+Ξ=1$⁠, and

Also, one can express the amplitude vector as

[here, there is no $O(ϵ⊥)$ correction, unlike in Eq. (15)], whose squared length $|a|2≡|ao|2+|ax|2$ satisfies

up to $O(ϵ∥)$⁠. Since we consider the beam dynamics in coordinates that are close to Euclidean, it will be sufficient, within the accuracy of our model, to adopt²⁰ 

Then, $Ξ+$ is simply the Hermitian conjugate of $Ξ$⁠.

We consider the wave evolution in curvilinear coordinates that are linked to a “reference ray” (RR), which is governed by

Here, ζ is the path along the ray, X and K are the RR coordinate and wave vector,

is the group velocity, $V⋆ ≐ |V⋆|$⁠, and the index $⋆$ denotes that the corresponding quantity is evaluated on $(X,K)$⁠. In particular, $H⋆ ≐ H(X,K)$⁠, and H is defined as follows:

We require $H⋆$ to be exactly zero initially, which is ensured by choosing an appropriate K; then, $H⋆$ remains zero at all ζ, as seen from Eq. (21).

The RR-based coordinates are introduced as $x̃≡{ζ,ϱ̃1,ϱ̃2}$⁠, where $ϱ̃σ$ are orthogonal coordinates transverse to the RR as specified in Paper II. (Here and further, the indices σ and $σ¯$ span from 1 to 2; other Greek indices span from 1 to 3.) The basis vectors $ẽμ$ of the new coordinates (⁠$dx=ẽμdx̃μ$⁠) are defined such that

Then,

For the transformation, matrices are defined as

which leads to

The specific choice of $ẽ⋆μ$ is described in Paper II.

We shall use tilde to denote the components of vectors and tensors measured in the RR-based coordinates $x̃$⁠. These components can be mapped to those in the laboratory coordinates using the standard formulas.²¹ For example, for any vector or covector A, one has

We also introduce first-order partial derivatives

and the second-order derivatives are denoted as follows:

Accordingly, $f̃|μ ≐ ∂f(x,k)/∂x̃μ$⁠, and so on.

Let us split the matrix (17) into its scalar part H and its traceless part M; this gives $Λ=H1+M$⁠. We assume $M(x,k(x)) ≲ O(ϵ⊥)$ and $M⋆ ≲ O(ϵ∥)$⁠; then,

where $π̃μ ≐ k̃μ(x)−K̃μ(ζ)$ is given by $π̃μ=−ϱ̃σH̃⋆|σṼ⋆μ/V⋆2$⁠.²¹ We can also rewrite this as

This leads to

Let us also introduce the matrices $uα ≐ M|α≲O(ϵ⊥)$⁠,

and a rescaled amplitude vector

By combining the results of Papers I and II, one readily finds that $ϕ$ is governed by the following parabolic partial differential equation, which is the main “quasi-optical” equation used in PARADE

Here, we introduced the notation

and the following coefficients:

(As a reminder, the index A in the latter formula denotes the anti-Hermitian part.) Also,

the matrix $ϑ̃⋆σσ$ is obtained from here by setting $σ¯=σ$ and summing over σ accordingly, and

Finally,

Note that the structure of the two-mode quasi-optical equation (39) is the same as that of the single-mode equation in Paper II except for the following: (i) there are additional terms $ũ⋆$⁠, $M̃⋆$⁠, and $M⋆$⁠, and (ii) the coefficients are matrices rather than scalars. In particular, U and $Γ$ are generally non-diagonal and thus can cause mode conversion. Also note that, like in Paper II, one finds the following corollary:

where $P ≐ ∫|ϕ|2 d2ϱ̃$ equals, up to a constant coefficient, the energy flux carried by the beam. This shows that the total energy flux of the beam is conserved when $Γ=0$⁠.

To facilitate benchmarking of PARADE, we also introduce the electromagnetic-field parametrization in terms of the polarization angles²⁶ used in the 1DFW code presented in Ref. 14. The two components of the field envelope $ψ$ projected on the transverse space ${ϱ̃1,ϱ̃2}$ can be expressed as follows:

The geometrical meaning of the polarization angles α and β can be understood from Fig. 1. Explicitly, these angles can be expressed through $ψ$ as²⁶ 

(note the difference in the signs in the denominators), where $δ ≐ arg(ψ̃⋆2)−arg(ψ̃⋆1)$⁠. In the following PARADE simulations, the initial a is prescribed by prescribing α and β. Specifically, α and β determine $ψ̃μ$ via Eq. (45), and Eq. (18) yields $as=ηs*αXαμψ̃μ$⁠.

Here, we present PARADE simulations of mode-converting wave beams with the focus on two effects: (i) the mode-amplitude transformation during the O–X conversion caused by the magnetic-field shear and (ii) splitting of beams that consist of multiple modes. The simulation algorithm is the same as the one used in PARADE for single-mode waves in Paper II. Also like in Paper II, all simulations were performed on a laptop with Intel Core^TM i7-7660U processor and took only a few seconds to run, as further specified in the figure captions. For D, we assume the dispersion tensor of collisionless cold electron plasma for simplicity, and so there is no dissipation.

As mentioned in Sec. II D, mode conversion in the vector equation (39) is governed by non-diagonal matrices $U⋆$ and $Γ$⁠. In the cold-plasma model, $Γ$ is zero, but $U⋆$ is generally non-negligible. In order to illustrate the effect of $U⋆$ on the polarization state of wave beams in low-density plasmas, we compared PARADE simulations using Eq. (39) with those using the single-mode scheme. This scheme is described in Paper II, and it is also similar to other existing quasi-optical models^3–13 in that it ignores the coupling between the electromagnetic modes.

As an example, we chose the initial wave to be a pure O mode. Also, the assumed geometry is as follows. We introduce the standard notation {x, y, z} for the laboratory coordinates. The origin is chosen to be the RR starting point, and $ez$⁠, which is the unit vector along the z axis, is chosen to be the orientation of the RR initial wave vector. Then, we assume a slab geometry with electron density

and magnetic field

The parameters are $n0=1.0×1018$ m⁻³, $z0=0.9$ m, $Ln=0.9$ m, $B0=0.4$ T, $θo=80.0°$⁠, $θs=80.0°$⁠, and $Lb=0.9$ m. The polarization angles are chosen to be $α=80.0°$ and $β=−27.0°$⁠. We also adopt the initial beam profile as Gaussian^21,27 with the focal lengths of $Z1=Z2=4.0$ m and waist sizes of $w0,1=w0,2=5.0$ cm. The wave frequency is f = 77.0 GHz, which corresponds to the vacuum wavelength of $λ0≈4$ mm.

The simulation results are presented in Fig. 2. It is seen that, when the vector model (39) is used (solid lines), the variation of the polarization angles is slow and overall small [Fig. 2(d)] because the plasma density is low ( $fpc<0.2$ in units f ). Correspondingly, the relative intensities of the O and X waves, which are defined as

vary rapidly [Fig. 2(e)], because of the strong shear of the magnetic field with the scale length of $Lb=0.9$ m. This variation illustrates the shear-driven mode conversion. In contrast, in the single-mode scheme (dashed lines), where the mode conversion is not taken into account, the initially excited O mode remains pure. Then, h^s values are fixed and, accordingly, the polarization angles vary rapidly although the ambient plasma has low density.

As the second example, we compared predictions of PARADE with predictions of the 1DFW code presented in Ref. 14, so as to verify that the shear-driven mode conversion and smooth variation of the polarization state are modeled accurately. Here, the initial polarization angles are chosen to be $α=35.0°$ and $β=−10.0°$⁠. The simulations are performed in a slab geometry [Eqs. (48) and (49)] with $n0=3.0× 1018$ m⁻³ and $Lb=5.4$ m; the other parameters are the same as in Sec. III A 1. The comparison shows good agreement of the two codes (Fig. 3).

As another example, we consider how the mode conversion is influenced by the magnetic-field shear and by the plasma density. For simplicity, we adopt that the magnetic field is given by Eq. (49), as earlier, whereas the density is constant, $n=n0$⁠. First, the O–X mode conversion is simulated with fixed $n0=1.0×1018$ m⁻³ and different scales of the shear, ranging from $Lb=0.3$ m to $Lb=3.6$ m (Fig. 4). The parameters other than n₀ and L_b are kept the same as in Sec. III A 2. Since the low density is assumed (⁠$fpc<0.2$ in units f), the polarization state does not change notably with L_b [Figs. 4(c) and 4(d)]. However, the relative intensities of the O and X waves [Eq. (50)] do because a^s are defined as the projections of $ψ$ on the polarization vectors [Eq. (18)] that are linked to B. As seen in Figs. 4(e) and 4(f), smaller L_b corresponds to shorter periods of the energy exchange between the O and X waves. Specifically, the mode conversion occurs twice per rotation of the magnetic field in the (x, y) plane.

We also simulated the O–X mode conversion for n₀ scanned in the range of $1.0×1018$ m⁻³ to $5.0×1018$ m⁻³ at fixed $Lb=0.9$ m. The other parameters were kept the same as in Sec. III A 2. As seen in Fig. 5, larger n₀ corresponds to stronger variations of the polarization angles [Figs. 5(c) and 5(d)] and slower variations of the relative intensities [Figs. 5(e) and 5(f)]. This tendency predicted by PARADE is anticipated because the two electromagnetic eigenmodes are nonresonant at sufficiently dense plasma but become resonant at low densities, where they both have the refraction indexes close to unity. The influence of the shear [Figs. 4(e) and 4(f)] and plasma density [Figs. 5(e) and 5(f)] on the relative intensities is particularly relevant in practice due to the fact that magnetically confined fusion plasmas have inhomogeneous shear and nonzero inhomogeneous density even outside the edge. Specific applications of PARADE and comparison with 1DFW simulations (where refraction is neglected) are reported in Ref. 24.

PARADE is also advantageous in that it can efficiently model splitting of multi-mode beams. To demonstrate this capability, we performed numerical simulations in a slab geometry with

Here, $n0=1.0×1019$ m⁻³, $x0=4.0$ m, $Ln=4.0$ m, $B0=1.0$ T, $Lb=4.0$ m, and $ez$ is the unit vector along the z axis. (The orientation of the magnetic field $B/|B|$ is chosen to be constant here to suppress the shear-driven mode conversion.) The origin is chosen to be the RR starting point, and ${x,y,z}={1.0,0.0,0.2}$ m is chosen to be the target point, which is used to fix the orientation of the RR initial wave vector. The initial polarization angles are $α=10.0°$ and $β=−30.0°$⁠. Also, for the initial beam profile, we adopt a Gaussian profile^21,27 with the focal lengths $Z1=Z2=3.0$ m and waist sizes $w0,1=w0,2=5.0$ cm. The wave frequency is chosen to be f = 77.0 GHz. Figure 6 shows the evolution of the transverse profile of $|a|$ at different locations along the RR. One can see the gradual splitting of the original beam into O-mode and X-mode beams propagating along separate ray trajectories. (Note that a single RR is used for this simulation, in contrast to single-mode simulations, where each mode would have its own RR.) Figure 7 shows the trajectories of the locations of the amplitude maxima. For comparison, this figure also shows the trajectories obtained from ray-tracing simulations, where O and X waves are modeled as independent. It is seen that PARADE's quasi-optical simulations are in good agreement with conventional ray tracing. Importantly, such quasi-optical modeling of a two-mode beam is adequate only as long as the group velocities of the O and X waves remain close enough to each other; otherwise, the ordering (4) cannot be maintained. (That is why we consider only weak splitting here.) However, by the time the two group velocities become very different, the O and X waves also become nonresonant and thus independent. Such waves can also be modeled with PARADE, except that the single-mode algorithm²¹ must be used instead.

This work continues a series of papers where we propose a new code PARADE for quasi-optical modeling of electromagnetic beams with and without mode conversion. The general theoretical model underlying PARADE and its application to single-mode beams were presented earlier.^20,21 Here, we apply PARADE to produce the first quasi-optical simulations of mode-converting beams. We also demonstrate that PARADE can model splitting of two-mode beams. The numerical results produced by PARADE show good agreement with those of one-dimensional full-wave simulations and also with conventional ray tracing to the extent those are applicable.

This work was supported by the U.S. DOE through Contract No. DE-AC02-09CH11466. This work was also supported by JSPS KAKENHI Grant No. JP17H03514.