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

TL;DR

This paper presents numerical simulations of mode-converting electromagnetic beams using a quasioptical modeling algorithm, demonstrating the PARADE code's ability to accurately simulate beam splitting and mode conversion phenomena.

Contribution

It introduces the first quasioptical simulations of mode-converting beams and validates the PARADE code against full-wave and ray-tracing simulations.

Findings

PARADE accurately models mode conversion and beam splitting.

Simulation results agree with full-wave and ray-tracing methods.

The approach advances quasioptical modeling capabilities.

Abstract

This work continues a series of papers where we propose an algorithm for quasioptical modeling of electromagnetic beams with and without mode conversion. The general theory was reported in the first paper of this series, where a parabolic partial differential equation was derived for the field envelope that may contain one or multiple modes with close group velocities. In the second paper, we presented a corresponding code PARADE (PAraxial RAy DEscription) and its test applications to single-mode beams. Here, we report quasioptical simulations of mode-converting beams for the first time. 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 that one-dimensional and ray-tracing simulations are applicable).

Figures

Equations

$$\hat{\mathrm{D}}\mathrm{E}=0,$$

$$\mathrm{E}=\psi e^{i\theta },$$

$$\hat{\mathcal{D}}\psi =0,\hat{\mathcal{D}}\doteq e^{-i\theta (\mathrm{x})}\hat{\mathrm{D}}{e}^{i\theta (\mathrm{x})},$$

$${ϵ}_{\parallel }\doteq \lambda /L_{\parallel },{ϵ}_{\perp }\doteq \lambda /L_{\perp },{ϵ}_{\parallel }\sim {ϵ}_{\perp }^2\ll 1.$$

$$\mathsf{D}\psi +{\hat{\mathcal{L}}}_{ϵ}\psi =0,$$

$${\mathsf{D}}_H=O(1),{\mathsf{D}}_A=O({ϵ}_{\parallel }),$$

$$\mathrm{D}(\mathrm{x},\mathrm{p})=\frac{c^2}{16\pi {\omega }^2}[\mathrm{p}\mathrm{p}-(\mathrm{p}\cdot \mathrm{p})1]+\frac{1}{16\pi }\epsilon (\mathrm{x},\mathrm{p}),$$

$${\mathrm{D}}_H\psi =O({ϵ}_{\perp }).$$

$$\psi ={\eta }_s{a}^s,$$

$$\displaystyle\boldsymbol{\rm D}_{H}\boldsymbol{\rm\psi}=\boldsymbol{\rm\eta}_{s}\Lambda_{s}a^{s},\quad\boldsymbol{\rm D}_{H}\boldsymbol{\rm\eta}_{s}=\Lambda_{s}\boldsymbol{\rm\eta}_{s}{\color[rgb]{0,0,0},}$$

$$\psi ={\eta }_{\mathrm{o}}{a}^{\mathrm{o}}+{\eta }_{\mathrm{x}}{a}^{\mathrm{x}}+\bar{\eta }\bar{a},$$

$$a^{\mathrm{o}}=O(1),a^{\mathrm{x}}=O(1),\bar{a}=O({ϵ}_{\perp }).$$

$$\displaystyle\boldsymbol{\rm a}=\left(\begin{array}[]{c}a^{\rm o}\\ a^{\rm x}\end{array}\right)$$

$$\displaystyle\boldsymbol{\rm\Xi}=\left(\begin{array}[]{cc}\boldsymbol{\rm\eta}_{\rm o}&\boldsymbol{\rm\eta}_{\rm x}\end{array}\right).$$

$$\psi =\mathrm{\Xi }\mathrm{a}+O({ϵ}_{\perp }).$$

$$\displaystyle\boldsymbol{\rm\Xi}^{+}=\left(\begin{array}[]{c}\boldsymbol{\rm\eta}^{{\rm o}*}\\ \boldsymbol{\rm\eta}^{{\rm x}*}\end{array}\right).$$

$$\displaystyle\boldsymbol{\rm\Lambda}\doteq\boldsymbol{\rm\Xi}^{+}\boldsymbol{\rm D}_{H}\boldsymbol{\rm\Xi}=\left(\begin{array}[]{cc}\Lambda_{\rm o}&0\\ 0&\Lambda_{\rm x}\end{array}\right).$$

$$\mathrm{a}={\mathrm{\Xi }}^{+}\psi$$

$$|\mathrm{a}{|}^2\approx |\psi {|}^2\equiv |{\psi }^x{|}^2+|{\psi }^y{|}^2+|{\psi }^z{|}^2$$

$${\eta }^{\mathrm{o}}\approx {\eta }_{\mathrm{o}}{\eta }^{\mathrm{x}}\approx {\eta }_{\mathrm{x}}.$$

$$\frac{\mathrm{d}{X}^{\alpha }}{\mathrm{d}\zeta }=\frac{V_{\star }^{\alpha }}{V_{\star }},\frac{\mathrm{d}{K}_{\alpha }}{\mathrm{d}\zeta }=-\frac{1}{V_{\star }}\frac{\partial H_{\star }}{\partial X^{\alpha }}.$$

$$V_{\star }^{\alpha }(\zeta )\doteq \frac{\partial H(\mathrm{X},\mathrm{K})}{\partial K_{\alpha }}$$

$$H(\mathrm{x},\mathrm{p})\doteq \frac{1}{2}[{\Lambda }_{\mathrm{o}}(\mathrm{x},\mathrm{p})+{\Lambda }_{\mathrm{x}}(\mathrm{x},\mathrm{p})].$$

$${\tilde{\mathrm{e}}}_{\star \mu }\cdot {\tilde{\mathrm{e}}}_{\star \nu }={\delta }_{\mu \nu },{\left[\partial {\tilde{\mathrm{e}}}_{\sigma }(\tilde{x})/\partial {\tilde{\varrho }}^{\bar{\sigma }}\right]}_{\star }=0.$$

$$\displaystyle\boldsymbol{\rm x}\approx\boldsymbol{\rm X}(\zeta)+\left(\begin{array}[]{cc}\tilde{\boldsymbol{\rm e}}_{\star 1}&\tilde{\boldsymbol{\rm e}}_{\star 2}\end{array}\right)\left(\begin{array}[]{c}\tilde{\varrho}^{1}\\ \tilde{\varrho}^{2}\end{array}\right).$$

$${\mathcal{X}}^{\alpha }{}_{\mu }\doteq \frac{\partial x^{\alpha }}{\partial {\tilde{x}}^{\mu }}=({\tilde{\mathrm{e}}}_{\mu }{)}^{\alpha },$$

$${\tilde{\mathcal{X}}}^{\mu }{}_{\alpha }\doteq \frac{\partial {\tilde{x}}^{\mu }}{\partial x^{\alpha }}=({\mathcal{X}}^{-1}{)}^{\mu }{}_{\alpha },$$

$${\mathcal{X}}^{\alpha }{}_{\zeta }\approx {\mathcal{X}}_{\star }^{\alpha }{}_{\zeta }+[({\tilde{\mathrm{e}}}_{\star \sigma }{)}^{\alpha }{]}^{\prime }{\tilde{\varrho }}^{\sigma },$$

$${\mathcal{X}}^{\alpha }{}_{\sigma }\approx {\mathcal{X}}_{\star }^{\alpha }{}_{\sigma },$$

$$\displaystyle\mathcal{X}^{\alpha}_{\star}{}_{\mu}=(\tilde{\boldsymbol{\rm e}}_{\star\mu})^{\alpha}{\color[rgb]{0,0,0}.}$$

$$A^{\alpha }={\mathcal{X}}^{\alpha }{}_{\mu }{\tilde{A}}^{\mu },A_{\alpha }={\tilde{A}}_{\mu }{\tilde{\mathcal{X}}}^{\mu }{}_{\alpha }.$$

$$f_{|\alpha }\doteq \frac{\partial f(\mathrm{x},\mathrm{k})}{\partial x^{\alpha }},f^{|\alpha }\doteq \frac{\partial f(\mathrm{x},\mathrm{k})}{\partial k_{\alpha }},$$

Full text

I Introduction

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 ref:tsujimura15 ; ref:marushchenko14 . Also, a lot of alternatives ref:mazzucato89 ; ref:nowak93 ; ref:peeters96 ; ref:farina07 ; ref:pereverzev98 ; ref:poli01b ; ref:poli01 ; ref:poli18 ; ref:balakin08b ; ref:balakin07a ; ref:balakin07b , 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 at 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 LABEL:ref:kubo15. 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) my:xo ; phd:ruiz17 ; my:covar ; my:qdirac ; my:qdiel , which is a reduced theory that describes refraction and mode conversion on the same footing. However, XGO, as it is formulated in Refs. my:xo ; phd:ruiz17 ; my:covar ; my:qdirac ; my:qdiel , still does not include diffraction, 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 foot:paper1 . In Paper II foot:paper2 , we presented a corresponding quasioptical 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 quasioptical 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 ref:bliokh15 ; foot:marius . For this reason, we illustrate the PARADE’s capabilities below on basic-physics examples. (Fusion-relevant applications of PARADE are discussed in LABEL:foot:itc27). 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.

II Theoretical model

II.1 Basic equations

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-mode waves. We assume a general linear equation for the electric field $\boldsymbol{\rm E}$ of a wave,

[TABLE]

where $\smash{\widehat{\boldsymbol{\rm D}}}$ is a linear dispersion operator. We also assume that this field can be represented in the eikonal form,

[TABLE]

where $\boldsymbol{\rm\psi}$ is a slow complex vector envelope and $\theta$ is a fast real “reference phase” to be prescribed. The wave is considered stationary, so it has a constant frequency $\omega$; then $\boldsymbol{\rm\psi}$ and $\theta$ are functions of the spatial coordinate $\boldsymbol{\rm x}$. Correspondingly, the envelope $\boldsymbol{\rm\psi}$ satisfies

[TABLE]

where $\smash{\widehat{\boldsymbol{\rm\mathcal{D}}}}$ serves as the “envelope dispersion operator” ($\doteq$ denotes definitions). We introduce $\boldsymbol{\rm k}\doteq\nabla\theta(\boldsymbol{\rm x})$ for the local wave vector, $\lambda\doteq 2\pi/k$ for the corresponding wavelength, $L_{\parallel}$ for the inhomogeneity scale of $\boldsymbol{\rm\psi}$ along the beam, and $L_{\perp}$ for the minimum scale of $\boldsymbol{\rm\psi}$ across the beam. The medium-inhomogeneity scale is assumed to be larger than or comparable with $L_{\parallel}$, and we adopt

[TABLE]

Then, Eq. (3) becomes

[TABLE]

where $\boldsymbol{\rm\mathsf{D}}$ serves as the “effective dispersion tensor” found from $\smash{\widehat{\boldsymbol{\rm D}}}$ and the operator $\smash{\widehat{\boldsymbol{\rm\mathcal{L}}}}_{\epsilon}=O(\epsilon_{\perp})$ is specified in Paper I (also see below). We suppose the following ordering:

[TABLE]

where the indices $H$ and $A$ denote the Hermitian and anti-Hermitian parts of $\smash{\widehat{\boldsymbol{\rm D}}}$, respectively. Assuming that the spatial dispersion is weak, $\boldsymbol{\rm\mathsf{D}}$ can be replaced with the homogeneous-plasma dispersion tensor,

[TABLE]

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

With $\boldsymbol{\rm\mathsf{D}}\approx\boldsymbol{\rm D}$ and $\boldsymbol{\rm D}$ given by Eq. (7), Eq. (6) implies $\boldsymbol{\rm\varepsilon}_{A}=O(\epsilon_{\parallel})$. Hence, $\boldsymbol{\rm D}_{H}$ is the dominant part of $\boldsymbol{\rm D}$, and Eq. (5) yields

[TABLE]

Since the Hermitian matrix $\boldsymbol{\rm D}_{H}=O(1)$ is the dominant part of $\boldsymbol{\rm D}$ and has enough eigenvectors $\boldsymbol{\rm\eta}_{s}$ to form a complete orthonormal basis, it is convenient to represent the envelope $\boldsymbol{\rm\psi}$ in this basis,

[TABLE]

where $a^{s}$ are the complex amplitudes. [Summation over repeating indices is assumed. For all functions derived from $\boldsymbol{\rm D}$, such as $\boldsymbol{\rm\eta}_{s}$, the notation convention $f\equiv f(\boldsymbol{\rm x})\equiv f(\boldsymbol{\rm x},\boldsymbol{\rm k}(\boldsymbol{\rm x}))$ will also be assumed by default.] Then,

[TABLE]

where $\Lambda_{s}$ are the corresponding eigenvalues. Due to Eq. (8) and the mutual orthogonality of all $\boldsymbol{\rm\eta}_{s}$, this means that $\Lambda_{s}a^{s}$ is small for every given $s$ individually; hence, either $a^{s}$ is small or $\Lambda_{s}$ is small. In the latter case, $\boldsymbol{\rm\eta}_{s}$ approximately satisfies the eigenmode equation, $\boldsymbol{\rm D}_{H}\boldsymbol{\rm\eta}_{s}=\Lambda_{s}\boldsymbol{\rm\eta}_{s}\approx 0$, 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 $a^{s}=O(1)$ is allowed.

II.2 Polarization matrices

In this paper, we consider the situation when there are two modes, henceforth called O and X modes, whose eigenvalues $\Lambda_{\rm o}\approx\Lambda_{\rm x}$ are both small within the region of interest. (This assumption is specified in Sec. II.4. 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,

[TABLE]

where $\boldsymbol{\rm\eta}_{\rm o}$ and $\boldsymbol{\rm\eta}_{\rm x}$ are the O- and X-mode polarization vectors, $\bar{\boldsymbol{\rm\eta}}$ is some third eigenvector of $\boldsymbol{\rm D}$ that is orthogonal to both of them, and

[TABLE]

The small amplitude $\bar{a}$ can be easily calculated as a perturbation and is included in our theory foot:paper1 but does not need to be considered below explicitly. Instead, we introduce a two-dimensional amplitude vector

[TABLE]

and the $3\times 2$ “polarization matrix” $\boldsymbol{\rm\Xi}$ that contains the vectors $\boldsymbol{\rm\eta}_{\rm o}$ and $\boldsymbol{\rm\eta}_{\rm x}$ as its columns,

[TABLE]

Then, $\boldsymbol{\rm\psi}$ can be expressed as follows:

[TABLE]

We also introduce the dual basis vectors $\boldsymbol{\rm\eta}^{\rm o}$ and $\boldsymbol{\rm\eta}^{\rm x}$ and an auxiliary polarization matrix

[TABLE]

As seen easily, this matrix satisfies $\boldsymbol{\rm\Xi}^{+}\boldsymbol{\rm\Xi}=\mathbb{1}$, and

[TABLE]

Also, one can express the amplitude vector as

[TABLE]

[here, there is no $O(\epsilon_{\perp})$ correction, unlike in Eq. (18)], whose squared length $|\boldsymbol{\rm a}|^{2}\equiv|a^{\rm o}|^{2}+|a^{\rm x}|^{2}$ satisfies

[TABLE]

up to $O(\epsilon_{\parallel})$. 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 foot:paper1

[TABLE]

Then, $\boldsymbol{\rm\Xi}^{+}$ is simply the Hermitian conjugate of $\boldsymbol{\rm\Xi}$.

II.3 Reference ray and new coordinates

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

[TABLE]

Here, $\zeta$ is the path along the ray, $\boldsymbol{\rm X}$ and $\boldsymbol{\rm K}$ are the RR coordinate and wave vector,

[TABLE]

is the group velocity, $V_{\star}\doteq|\boldsymbol{\rm V}_{\star}|$, and the index $\star$ denotes that the corresponding quantity is evaluated on $(\boldsymbol{\rm X},\boldsymbol{\rm K})$. In particular, $H_{\star}\doteq H(\boldsymbol{\rm X},\boldsymbol{\rm K})$, and $H$ is defined as follows:

[TABLE]

We require $H_{\star}$ to be exactly zero initially, which is ensured by choosing an appropriate $\boldsymbol{\rm K}$; then, $H_{\star}$ remains zero at all $\zeta$, as seen from Eqs. (28).

The RR-based coordinates are introduced as $\tilde{x}\equiv\{\zeta,\tilde{\varrho}^{1},\tilde{\varrho}^{2}\}$, where $\tilde{\varrho}^{\sigma}$ are orthogonal coordinates transverse to the RR as specified in Paper II. (Here and further, the indices $\sigma$ and ${\bar{\sigma}}$ span from 1 to 2; other Greek indices span from 1 to 3.) The basis vectors $\tilde{\boldsymbol{\rm e}}_{\mu}$ of the new coordinates ($\mathrm{d}\boldsymbol{\rm x}=\tilde{\boldsymbol{\rm e}}_{\mu}\mathrm{d}\tilde{x}^{\mu}$) are defined such that

[TABLE]

Then,

[TABLE]

For the transformation matrices defined as

[TABLE]

this leads to

[TABLE]

The specific choice of $\tilde{\boldsymbol{\rm e}}_{\star\mu}$ is described in Paper II.

We shall use tilde to denote the components of vectors and tensors measured in the RR-based coordinates $\tilde{x}$. These components can be mapped to those in the laboratory coordinates using the standard formulas foot:paper2 . For example, for any vector or covector $\boldsymbol{\rm A}$, one has

[TABLE]

We also introduce first-order partial derivatives

[TABLE]

and the second-order derivatives are denoted as follows:

[TABLE]

Accordingly, $\tilde{f}_{|\mu}\doteq\partial f(\boldsymbol{\rm x},\boldsymbol{\rm k})/\partial\tilde{x}^{\mu}$, and so on.

II.4 Quasioptical equation

Let us split the matrix (24) into its scalar part $H$ and its traceless part $\boldsymbol{\rm M}$; this gives $\boldsymbol{\rm\Lambda}=H\mathbb{1}+\boldsymbol{\rm M}$. We assume $\boldsymbol{\rm M}(\boldsymbol{\rm x},\boldsymbol{\rm k}(\boldsymbol{\rm x}))\lesssim O(\epsilon_{\perp})$ and $\boldsymbol{\rm M}_{\star}\lesssim O(\epsilon_{\parallel})$; then,

[TABLE]

where $\tilde{\pi}_{\mu}\doteq\tilde{k}_{\mu}(\boldsymbol{\rm x})-\tilde{K}_{\mu}(\zeta)$ is given by $\tilde{\pi}_{\mu}=-\tilde{\varrho}^{\sigma}\tilde{H}_{\star|\sigma}\tilde{V}_{\star\mu}/V_{\star}^{2}$ foot:paper2 . We can also rewrite this as

[TABLE]

This leads to

[TABLE]

Let us also introduce the matrices $\boldsymbol{\rm u}^{\alpha}\doteq\boldsymbol{\rm M}^{|\alpha}\lesssim O(\epsilon_{\perp})$,

[TABLE]

and a rescaled amplitude vector

[TABLE]

By combining the results of Papers I and II, one readily finds that $\boldsymbol{\rm\phi}$ is governed by the following parabolic partial differential equation, which is the main “quasioptical” equation used in PARADE:

[TABLE]

Here, we introduced the notation

[TABLE]

and the following coefficients:

[TABLE]

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

[TABLE]

the matrix $\tilde{\boldsymbol{\rm\vartheta}}_{\star}^{\sigma}{}_{\sigma}$ is obtained from here by setting ${\bar{\sigma}}=\sigma$ and summing over $\sigma$ accordingly, and

[TABLE]

Finally,

[TABLE]

Note that the structure of the two-mode quasioptical equation (49) is the same as that of the single-mode equation in Paper II except for the following: (i) there are additional terms $\tilde{\boldsymbol{\rm u}}_{\star}$, $\smash{\tilde{\boldsymbol{\rm\mathfrak{M}}}_{\star}}$, and $\boldsymbol{\rm M}_{\star}$; and (ii) the coefficients are matrices rather than scalars. In particular, $\boldsymbol{\rm U}$ and $\boldsymbol{\rm\Gamma}$ are generally non-diagonal and thus can cause mode conversion. Also note that, like in Paper II, one finds the following corollary:

[TABLE]

where $P\doteq\int|\boldsymbol{\rm\phi}|^{2}\,d^{2}\tilde{\varrho}$ 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 $\boldsymbol{\rm\Gamma}=0$.

II.5 Polarization angles

To facilitate benchmarking of PARADE, we also introduce the electromagnetic-field parametrization in terms of the polarization angles book:born65 used in the 1DFW code presented in LABEL:ref:kubo15. The two components of the field envelope $\boldsymbol{\rm\psi}$ projected on the transverse space $\{\tilde{\varrho}^{1},\tilde{\varrho}^{2}\}$ can be expressed as follows:

[TABLE]

The geometrical meaning of the polarization angles $\alpha$ and $\beta$ can be understood from Fig. 1. Explicitly, these angles can be expressed through $\boldsymbol{\rm\psi}$ as book:born65

[TABLE]

(note the difference in the signs in the denominators), where $\delta\doteq\arg(\tilde{\psi}^{2}_{\star})-\arg(\tilde{\psi}^{1}_{\star})$. In the following PARADE simulations, the initial $\boldsymbol{\rm a}$ is prescribed by prescribing $\alpha$ and $\beta$. Specifically, $\alpha$ and $\beta$ determine $\tilde{\psi}^{\mu}$ via Eq. (62), and Eq. (25) yields $a^{s}=\eta^{s*}{}_{\alpha}\mathcal{X}^{\alpha}{}_{\mu}\tilde{\psi}^{\mu}$.

III Simulation results

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 done on a laptop with Intel CoreTM i7-7660U processor and took only a few seconds to run, as further specified in the figure captions. For $\boldsymbol{\rm D}$, we assume the dispersion tensor of collisionless cold electron plasma for simplicity, so there is no dissipation.

III.1 Mode conversion in a sheared magnetic field

III.1.1 Comparison with uncoupled-mode simulations

As mentioned in Sec. II.4, mode conversion in the vector equation (49) is governed by non-diagonal matrices $\boldsymbol{\rm U}_{\star}$ and $\boldsymbol{\rm\Gamma}$. In the cold-plasma model, $\boldsymbol{\rm\Gamma}$ is zero but $\boldsymbol{\rm U}_{\star}$ is generally non-negligible. In order to illustrate the effect of $\boldsymbol{\rm U}_{\star}$ on the polarization state of wave beams in low-density plasmas, we compared PARADE simulations using Eq. (49) with those using the single-mode scheme. This scheme is described in Paper II, and it is also similar to other existing quasioptical models ref:mazzucato89 ; ref:nowak93 ; ref:peeters96 ; ref:farina07 ; ref:pereverzev98 ; ref:poli01b ; ref:poli01 ; ref:poli18 ; ref:balakin08b ; ref:balakin07a ; ref:balakin07b 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 $\boldsymbol{\rm e}_{z}$, 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

[TABLE]

and magnetic field

[TABLE]

The parameters are $n_{0}=1.0\times 10^{18}$ m*-3*, $z_{0}=0.9$ m, $L_{n}=0.9$ m, $B_{0}=0.4$ T, $\theta_{o}=80.0^{\circ}$, $\theta_{s}=80.0^{\circ}$, and $L_{b}=0.9$ m. The polarization angles are chosen to be $\alpha=80.0^{\circ}$ and $\beta=-27.0^{\circ}$. We also adopt the initial beam profile as Gaussian foot:paper2 ; book:yariv with the focal lengths $\mathcal{Z}_{1}=\mathcal{Z}_{2}=4.0$ m and waist sizes $w_{0,1}=w_{0,2}=5.0$ cm. The wave frequency is $f=77.0$ GHz, which corresponds to the vacuum wavelength $\lambda_{0}\approx 4$ mm.

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

[TABLE]

vary rapidly [Fig. 2(e)], because of the strong shear of the magnetic field with scale length $L_{b}=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}$ are fixed and, accordingly, the polarization angles vary rapidly, although the ambient plasma has low density.

III.1.2 Comparison with the 1DFW code

As the second example, we compared predictions of PARADE with predictions of the 1DFW code presented in LABEL:ref:kubo15, 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 $\alpha=35.0^{\circ}$ and $\beta=-10.0^{\circ}$. The simulations are performed in a slab geometry [Eqs. (65) and (72)] with $n_{0}=3.0\times 10^{18}$ m*-3* and $L_{b}=5.4$ m; the other parameters are the same as in Sec. III.1.1. The comparison shows good agreement of the two codes (Fig. 3).

III.1.3 Parameter scan

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. (72), as earlier, whereas the density is constant, $n=n_{0}$. First, the O–X mode conversion is simulated with fixed $n_{0}=1.0\times 10^{18}$ m*-3* and different scales of the shear, ranging from $L_{b}=0.3$ m to $L_{b}=3.6$ m (Fig. 4). The parameters other than $n_{0}$ and $L_{b}$ are kept the same as in Sec. III.1.2. Since the low density is assumed ($f_{\rm pc}<0.2$ in units $f$), the polarization state does not change notably with $L_{b}$ [Figs. 4(c) and (d)]. However, the relative intensities of the O and X waves [Eq. (73)] do, because $a^{s}$ are defined as the projections of $\boldsymbol{\rm\psi}$ on the polarization vectors [Eq. (25)] that are linked to $\boldsymbol{\rm B}$. As seen in Figs. 4(e) and (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_{0}$ scanned in the range between $1.0\times 10^{18}$ m*-3* to $5.0\times 10^{18}$ m*-3* at fixed $L_{b}=0.9$ m. The other parameters were kept the same as in Sec. III.1.2. As seen in Fig. 5, larger $n_{0}$ corresponds to stronger variations of the polarization angles [Figs. 5(c) and (d)] and slower variations of the relative intensities [Figs. 5(e) and (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 (f)] and plasma density [Figs. 5(e) and (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 LABEL:foot:itc27.

III.2 Weak splitting of mode-converting beams

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

[TABLE]

Here, $n_{0}=1.0\times 10^{19}$ m*-3*, $x_{0}=4.0$ m, $L_{n}=4.0$ m, $B_{0}=1.0$ T, $L_{b}=4.0$ m, and $\boldsymbol{\rm e}_{z}$ is the unit vector along the $z$ axis. (The orientation of the magnetic field $\boldsymbol{\rm B}/|\boldsymbol{\rm 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 $\alpha=10.0^{\circ}$ and $\beta=-30.0^{\circ}$. Also, for the initial beam profile, we adopt a Gaussian profile foot:paper2 ; book:yariv with the focal lengths $\mathcal{Z}_{1}=\mathcal{Z}_{2}=3.0$ m and waist sizes $w_{0,1}=w_{0,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 $|\boldsymbol{\rm 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 quasioptical simulations are in good agreement with conventional ray tracing. Importantly, such quasioptical 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 the single-mode algorithm foot:paper2 must be used instead.

IV Conclusions

This work continues a series of papers where we propose a new code PARADE for quasioptical 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 foot:paper1 ; foot:paper2 . Here, we apply PARADE to produce the first quasioptical 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.

V Acknowledgments

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