Quasioptical modeling of wave beams with and without mode conversion: II. Numerical simulations of single-mode beams

TL;DR

This paper introduces a numerical code for quasioptical modeling of electromagnetic beams, validating it against analytical and ray-tracing methods, focusing on single-mode beams in vacuum and plasma.

Contribution

It presents the PARADE code for simulating single-mode electromagnetic beams and compares its results with analytical and ray-tracing solutions.

Findings

PARADE accurately models single-mode beams in vacuum and plasma.

Numerical results agree with Gaussian-beam optics and cold-plasma ray tracing.

The code enables detailed quasioptical simulations for beam analysis.

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. Here, we present a corresponding code PARADE (PAraxial RAy DEscription) and its test applications to single-mode beams in vacuum and also in inhomogeneous magnetized plasma. The numerical results are compared, respectively, with analytic formulas from Gaussian-beam optics and also with cold-plasma ray tracing. Quasioptical simulations of mode-converting beams are reported in the next, third paper of this series.

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 }\hat{\mathrm{D}}{e}^{i\theta },$$

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

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

$${\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,$$

$${\mathrm{D}}_H\psi ={\eta }_s{\Lambda }_s{a}^s,{\mathrm{D}}_H{\eta }_s={\Lambda }_s{\eta }_s,$$

$$\psi =\mathrm{\Xi }a+O({ϵ}_{\perp }),{\mathrm{D}}_H\mathrm{\Xi }=\mathrm{\Xi }\Lambda ,$$

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

$$\frac{\mathrm{d}{x}^{\alpha }}{\mathrm{d}\tau }=\frac{\partial \Lambda (\mathrm{x},\mathrm{p})}{\partial p_{\alpha }},\frac{\mathrm{d}{p}_{\alpha }}{\mathrm{d}\tau }=-\frac{\partial \Lambda (\mathrm{x},\mathrm{p})}{\partial x^{\alpha }},$$

$$V^{\alpha }(\mathrm{x},\mathrm{p})\doteq \frac{\partial \Lambda (\mathrm{x},\mathrm{p})}{\partial p_{\alpha }}.$$

$$\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 {\Lambda }_{\star }}{\partial X^{\alpha }},$$

$${\Lambda }_{\star }=0.$$

$$\mathrm{d}\mathrm{x}={\tilde{\mathrm{e}}}_{\mu }\mathrm{d}{\tilde{x}}^{\mu },$$

$$\tilde{\mathcal{X}}\doteq {\mathcal{X}}^{-1},$$

$${\mathcal{X}}^{\alpha }{}_{\mu }=\frac{\partial x^{\alpha }}{\partial {\tilde{x}}^{\mu }},{\tilde{\mathcal{X}}}^{\mu }{}_{\alpha }=\frac{\partial {\tilde{x}}^{\mu }}{\partial x^{\alpha }}.$$

$${\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 }{}_{\zeta }\approx ({\tilde{\mathrm{e}}}_{\star \zeta }{)}^{\alpha }+[({\tilde{\mathrm{e}}}_{\star \sigma }{)}^{\alpha }{]}^{\prime }{\tilde{\varrho }}^{\sigma },{\mathcal{X}}^{\alpha }{}_{\sigma }\approx ({\tilde{\mathrm{e}}}_{\star \sigma }{)}^{\alpha },$$

$$\mathcal{X}={\mathcal{X}}_{\star }+\Delta \mathcal{X},\Delta \mathcal{X}\approx {\mathcal{Y}}_{\star \sigma }{\tilde{\varrho }}^{\sigma },$$

$$\displaystyle\boldsymbol{\rm\mathcal{Y}}_{\star\sigma}=\left(\begin{array}[]{ccc}[(\tilde{\boldsymbol{\rm e}}_{\star\sigma})^{1}]^{\prime}&0&0\\[3.0pt] [(\tilde{\boldsymbol{\rm e}}_{\star\sigma})^{2}]^{\prime}&0&0\\[3.0pt] [(\tilde{\boldsymbol{\rm e}}_{\star\sigma})^{3}]^{\prime}&0&0\end{array}\right).$$

$$\Delta \tilde{\mathcal{X}}\approx -{\tilde{\mathcal{X}}}_{\star }(\Delta \mathcal{X}){\tilde{\mathcal{X}}}_{\star }\approx -{\tilde{\mathcal{X}}}_{\star }{\mathcal{Y}}_{\star \sigma }{\tilde{\mathcal{X}}}_{\star }{\tilde{\varrho }}^{\sigma }.$$

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

$$A^{\alpha \beta }={\mathcal{X}}^{\alpha }{}_{\mu }{\mathcal{X}}^{\beta }{}_{\nu }{\tilde{A}}^{\mu \nu },$$

$$A^{\alpha }{}_{\beta }={\mathcal{X}}^{\alpha }{}_{\mu }{\tilde{\mathcal{X}}}^{\nu }{}_{\beta }{\tilde{A}}^{\mu }{}_{\nu },$$

$$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 }},$$

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

$$\begin{gathered} \Lambda a+i\Gamma a-i{\tilde{V}}^{\mu }{\partial }_{\mu }a-\frac{i}{2}({\partial }_{\mu }{\tilde{V}}^{\mu })a-Ua \\ -\frac{1}{2}{\partial }_{\sigma }({\tilde{\Phi }}^{\sigma \bar{\sigma }}{\partial }_{\bar{\sigma }}a)=0. \end{gathered}$$

Full text

I Introduction

Electron-cyclotron heating and current drive in fusion plasmas involve electron-cyclotron waves (ECWs) and require deposition of the wave power with high precision; hence, accurate modeling of ECWs is required. Full-wave simulations of ECWs foot:vdovin06 are computationally-demanding and largely impractical due to the fact that the ECW wavelengths (which are in the mm range) are extremely small compared to the plasma size. Geometrical-optics (GO) ray-tracing simulations ref:tsujimura15 ; ref:marushchenko14 are much faster but cannot adequately describe the structure of the wave beam in the focal region. Many alternatives have been proposed to improve reduced modeling of ECWs, including complex-eikonal methods ref:mazzucato89 ; ref:nowak93 ; ref:peeters96 ; ref:farina07 , beam tracing ref:pereverzev98 ; ref:poli01b ; ref:poli01 ; ref:poli18 , and the quasioptical approach ref:balakin08b ; ref:balakin07a ; ref:balakin07b . However, none of them account for the possible interaction of the two cold-plasma modes (linear mode conversion), which can occur in strongly sheared magnetic field at low density, for example, at the plasma edge my:xo ; ref:kubo15 . Hence, a need remains for a more systematic approach to modeling ECWs.

This work continues a series of papers where we propose a new algorithm for quasioptical modeling of electromagnetic (EM) beams propagating in inhomogeneous plasma with and without mode conversion. Although our project is largely motivated by the need to improve ECW modeling in fusion plasmas, our results are equally applicable also beyond plasma physics, for example, to studying related problems in optics and general relativity ref:bliokh15 ; foot:marius . The general theory was reported in Paper I of our series foot:paper1 . There, a parabolic partial differential equation was derived for a certain projection of the field envelope, whose dimension $N\geq 1$ equals the number of the resonantly-coupled modes. Based on that theory, which is a generalization of extended geometrical optics (XGO) phd:ruiz17 ; my:covar ; my:qdirac ; my:qdiel , we have also developed a code PARADE (PAraxial RAy DEscription) and applied it to several test problems. Paper III of our series foot:paper3 reports PARADE’s applications to quasioptical simulations of mode-converting beams with $N=2$. Here, we present a simplified version of the same code for $N=1$, in which case mode conversion is not considered. This “scalar” version of PARADE serves two goals: (i) it is a stepping stone towards a more general code reported in Paper III; and (ii) it is useful as an independent tool for quasioptical beam tracing without mode conversion in sufficiently dense plasma, especially when the absorption is strongly inhomogeneous and no simple ansatz can be assumed for the beam structure. (In this second capacity, PARADE can also be considered as an alternative to the algorithm proposed in Refs. ref:balakin07a ; ref:balakin07b ; ref:balakin08b . The latter is similar in spirit foot:balakincomp but is limited to single-mode beams by design.) But fusion applications of PARADE foot:fusion will not be discussed in this paper, which is focused on the general method rather than on physics of specific waves. Instead, we shall present examples of PARADE simulations in the basic-physics context, namely, by modeling beam propagation in vacuum and in cold magnetized plasma. Simple geometries will be assumed in order to benchmark the code and illustrate its features most transparently.

This 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 outline the numerical algorithm used in PARADE. In Sec. IV, we report simulation results for test problems. In Sec. V, 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 describing single-mode wave beams. We start by assuming 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 discussed below. We also assume that the 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 (Sec. II.1). 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$ for the local wave vector, $\lambda\doteq 2\pi/k$ for the wavelength, $L_{\parallel}$ for the characteristic scale of the beam field $\boldsymbol{\rm\psi}$ along the group velocity at the beam center, and $L_{\perp}$ for the minimum scale of the field in the plane transverse to this group velocity. The medium-inhomogeneity scale is assumed to be larger than or comparable with $L_{\parallel}$.

We adopt the quasioptical ordering, specifically,

[TABLE]

Using this ordering, the “envelope dispersion operator” $\smash{\widehat{\boldsymbol{\rm\mathcal{D}}}}$ can be expressed as

[TABLE]

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

[TABLE]

where the indices $H$ and $A$ denote the Hermitian and anti-Hermitian parts, correspondingly. Assuming that spatial dispersion is weak enough, $\boldsymbol{\rm\mathsf{D}}$ can be replaced foot:paper1 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]

As a Hermitian matrix, $\boldsymbol{\rm D}_{H}$ has enough (three) eigenvectors $\boldsymbol{\rm\eta}_{s}$ to form a complete orthonormal basis. Then, it is convenient to represent the envelope in this basis,

[TABLE]

where $a^{s}$ are the corresponding complex amplitudes. [Summation over repeated indices is assumed. For all functions derived from $\boldsymbol{\rm D}$, such as $\boldsymbol{\rm\eta}_{s}$ here, 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}=O(\epsilon_{\perp})$ 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, so that $a^{s}=O(1)$ is allowed. In this paper, we assume conditions such that only one mode can be excited at given $\boldsymbol{\rm x}$ near the fixed frequency of interest. Then,

[TABLE]

where the polarization vector $\boldsymbol{\rm\Xi}$ is the eigenvector $\boldsymbol{\rm\eta}_{s}$ corresponding to the excited mode and $\Lambda$ is the corresponding eigenvalue. (The notation is chosen here such that our formulas extend to mode-converting waves with generalized definitions of $\boldsymbol{\rm\Xi}$, $\Lambda$, and $a$ as in Refs. foot:paper1 ; foot:paper3 .) The $O(\epsilon_{\perp})$ term in Eq. (11) can be easily calculated from $a$ and is included in our theory as a perturbation foot:paper1 but will not be considered below explicitly. Hence, we shall describe the wave in terms of the scalar $a$, which can also be expressed as follows:

[TABLE]

Here, the row vector $\boldsymbol{\rm\Xi}^{+}$ is the dual polarization vector, so $\boldsymbol{\rm\Xi}^{+}\boldsymbol{\rm\Xi}=1$. Since we consider the beam dynamics in coordinates that are close to Euclidean, it will be sufficient, within the accuracy of our model, to consider the row vector $\boldsymbol{\rm\Xi}^{+}$ simply as the conjugate transpose of the column vector $\boldsymbol{\rm\Xi}$.

II.2 Reference ray

We shall consider the beam evolution in the coordinate system that is determined by a family of rays. These rays are governed by Hamilton’s equations foot:paper1

[TABLE]

where $\tau$ serves as an effective time variable, $x^{\alpha}$ denote the ray location in the laboratory coordinates, and $p_{\alpha}$ denote the ray wave-vector components in the associated momentum space. The initial conditions are assumed such that $\Lambda(\boldsymbol{\rm x}(0),\boldsymbol{\rm p}(0))$ is of order $\epsilon_{\parallel}$ or smaller foot:paper1 . (A more specific definition will not be needed.) Then, it also remains small indefinitely, because $\Lambda$ is conserved by Eqs. (13).

Let us introduce the path along each ray, $\zeta\doteq\int_{0}^{\tau}V\,d\tau$, where $V\doteq|\boldsymbol{\rm V}|$ and $\boldsymbol{\rm V}$ is given by

[TABLE]

At any given $\zeta$, the coordinates of rays with different initial $\boldsymbol{\rm x}$ form a “transverse space”. Let us define a two-dimensional transverse coordinate $\tilde{\boldsymbol{\rm\varrho}}\equiv\{\tilde{\varrho}^{1},\tilde{\varrho}^{2}\}$ on this space such that its origin corresponds to the location of some “reference ray” (RR). This ray can be chosen arbitrarily as long as it is representative of the rays within the beam. The coordinate and the wave vector of the RR are denoted as $\boldsymbol{\rm X}$ and $\boldsymbol{\rm K}$, respectively. Like the coordinate and momentum of other rays [Eqs. (13)], they can be found by solving the ray equations

[TABLE]

where the index $\star$ denotes evaluation at $(\boldsymbol{\rm X},\boldsymbol{\rm K})$. Note that we require $\Lambda_{\star}$ to be exactly zero initially; then, it also remains zero at all $\zeta$,

[TABLE]

II.3 Ray-based coordinates

Since there is a one-to-one correspondence between the laboratory coordinates $x^{\alpha}\equiv\{x^{1},x^{2},x^{3}\}$ and $\tilde{x}^{\mu}\equiv\{\zeta,\tilde{\varrho}^{1},\tilde{\varrho}^{2}\}$, one can choose the latter as the new coordinates. The two sets of coordinates are connected via

[TABLE]

where $\tilde{\boldsymbol{\rm e}}_{\mu}$ are the basis vectors of the “tilded” coordinate system. Accordingly, the transformation matrices $\boldsymbol{\rm\mathcal{X}}$ and

[TABLE]

where $\mathcal{X}^{\alpha}{}_{\mu}\doteq(\tilde{\boldsymbol{\rm e}}_{\mu})^{\alpha}$ are given by

[TABLE]

We choose the transverse coordinates such that

[TABLE]

(Here and further, the indices $\sigma$ and ${\bar{\sigma}}$ span from 1 to 2; other Greek indices span from 1 to 3. The specific choice of $\tilde{\boldsymbol{\rm e}}_{\mu}$ used in PARADE is discussed in Sec. III.) Then,

[TABLE]

and

[TABLE]

where the prime denotes the derivative with respect to $\zeta$. This leads to

[TABLE]

where $\boldsymbol{\rm\mathcal{Y}}_{\star\sigma}$ is a $3\times 3$ matrix given by

[TABLE]

Using that $\boldsymbol{\rm\mathcal{X}}\tilde{\boldsymbol{\rm\mathcal{X}}}=\mathbb{1}$, we also obtain $\tilde{\boldsymbol{\rm\mathcal{X}}}=\tilde{\boldsymbol{\rm\mathcal{X}}}_{\star}+\Delta\tilde{\boldsymbol{\rm\mathcal{X}}}$,

[TABLE]

For future references, note that we shall use tilde also to mark the components of vectors and tensors measured in these new coordinates. As usual, these components are connected to those in the laboratory coordinates via book:landau2

[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

II.4.1 Basic formulas

Using the fact that the metric in the assumed ray-based coordinates is Euclidean on the RR and overall smooth, the envelope equation can be represented as follows foot:paper1 :

[TABLE]

Here, we introduced

[TABLE]

and also the following notation for the derivatives:

[TABLE]

Note that all coefficients in Eq. (33) are considered as functions of $\boldsymbol{\rm x}$; namely, $f\equiv f(\boldsymbol{\rm x},\boldsymbol{\rm k}(\boldsymbol{\rm x}))$ for any $f$. Accordingly, the derivatives (37) act on both arguments. Also note that $\Gamma$ and $U$ are calculated in the laboratory coordinates; otherwise, additional terms would emerge foot:paper1 .

Within the assumed accuracy, $\tilde{\Phi}^{\sigma{\bar{\sigma}}}$ and $U$ can be replaced with those evaluated on the RR,

[TABLE]

In contrast, $\Lambda$, $\Gamma$, and $\tilde{V}^{\mu}$ generally must be calculated with a higher accuracy. (An example of a numerical simulation that ignores this fact is presented in Sec. IV.3 for a reference.) In order to do that, we need to find $\boldsymbol{\rm k}(\boldsymbol{\rm x})$, so we start by specifying $\theta$, through which $\boldsymbol{\rm k}$ is defined (Sec. II.1). Let us adopt

[TABLE]

(This choice of $\theta$ is different from the one assumed in Paper I but similar to that in LABEL:ref:balakin07a.) Hence,

[TABLE]

To simplify this, note that

[TABLE]

by definition. Hence, the transformation $(X^{\alpha},K^{\alpha})\mapsto(\tilde{X}^{\mu},\tilde{K}^{\mu})$ is symplectic (canonical), so the “tilded” coordinates satisfy Hamilton’s equations similar to Eq. (15),

[TABLE]

By combining Eqs. (41), (42), and (44), one arrives at

[TABLE]

where we used $\tilde{V}_{\star\mu}=V_{\star}\delta_{\mu}^{\zeta}$, and $\delta_{\mu}^{\zeta}$ is a Kronecker symbol.

II.4.2 Calculating $\Lambda$ and $\Gamma$

Expanding $\Lambda(\boldsymbol{\rm x},\boldsymbol{\rm k}(\boldsymbol{\rm x}))$ around the RR, one obtains

[TABLE]

Using Eqs. (16) and (45), one can reduce this down to

[TABLE]

which is of the order of $O(\epsilon_{\perp}^{2})\sim O(\epsilon_{\parallel})$, as needed. Since the metric on the RR is Euclidean, the second-order derivatives entering this equation are transformed as true tensors. Hence, they can be mapped to the laboratory coordinates with $\boldsymbol{\rm\mathcal{X}}_{\star}$ and $\tilde{\boldsymbol{\rm\mathcal{X}}}_{\star}$. This leads to $\Lambda(\boldsymbol{\rm x},\boldsymbol{\rm k}(\boldsymbol{\rm x}))\approx\tilde{\mathfrak{L}}_{\star\sigma{\bar{\sigma}}}\tilde{\varrho}^{\sigma}\tilde{\varrho}^{{\bar{\sigma}}}$ with

[TABLE]

where we substituted [cf. Eq. (45)]

[TABLE]

(Note that, to the leading order, $\boldsymbol{\rm\varrho}\equiv\{\varrho^{1},\varrho^{2}\}$ can be considered as a true vector and transforms accordingly.)

Using Eqs. (48), one can also adopt the following approximation of Eq. (34) for $\Gamma$:

[TABLE]

Here, the polarization vectors $\boldsymbol{\rm\Xi}$ and $\boldsymbol{\rm\Xi}^{+}$ are evaluated at the RR (hence the stars), because they are determined by $\boldsymbol{\rm D}_{H}$, which is a slow function. In contrast, $\boldsymbol{\rm D}_{A}$ can vary substantially within the beam cross-section, so it can be important to calculate $\boldsymbol{\rm D}_{A}$ separately at each point.

II.4.3 Calculating $\tilde{V}$

In order to calculate the terms involving $\tilde{V}^{\mu}=\tilde{\mathcal{X}}^{\mu}{}_{\alpha}V^{\alpha}$, let us represent this function as $\tilde{V}^{\mu}=\tilde{V}_{\star}^{\mu}+\tilde{v}^{\mu}$,

[TABLE]

where $\Delta$ denotes deviations from the corresponding values of the RR. Using Eq. (31), we obtain

[TABLE]

Also note that

[TABLE]

Hence, $\tilde{v}^{\mu}(\boldsymbol{\rm x})\approx\tilde{\vartheta}_{\star}^{\mu}{}_{{\bar{\sigma}}}\tilde{\varrho}^{{\bar{\sigma}}}$, where

[TABLE]

Then, since $\tilde{V}^{\sigma}_{\star}=0$ and $\tilde{v}^{\zeta}\partial_{\zeta}=O(\epsilon_{\parallel})$, one obtains

[TABLE]

where $\tilde{\vartheta}_{\star}^{\sigma}{}_{\sigma}$ can be calculated from Eq. (53) by taking $\mu=\bar{\sigma}=\sigma$ and summing over $\sigma$ accordingly. Also,

[TABLE]

II.4.4 Final result

Using the above results, one can express Eq. (33) as

[TABLE]

which also has the following corollary:

[TABLE]

Equation (57) shows that $P$ is conserved when $\Gamma$ is zero. As explained in Paper I, this represents the conservation of the wave action, and $P$ equals the energy flux up to a constant coefficient. Assuming the notation

[TABLE]

so $P=\int|\phi|^{2}\,\mathrm{d}^{2}\tilde{\varrho}$, one can also rewrite Eq. (56) as

[TABLE]

Equation (59) is the main quasioptical equation used in PARADE. As a reminder, we assume Eq. (55) for $V_{\star}$, Eq. (49) for $\Gamma$, Eq. (47) for $\tilde{\mathfrak{L}}_{\star\sigma{\bar{\sigma}}}$, Eq. (39) for $U_{\star}$, and Eq. (38) for $\tilde{\Phi}^{\sigma{\bar{\sigma}}}_{\star}$. We also assume Eq. (53) for $\tilde{\vartheta}_{\star}^{\sigma}{}_{{\bar{\sigma}}}$, with Eq. (30) for $\boldsymbol{\rm\mathcal{Y}}_{\star}$. The quantity $\tilde{\vartheta}_{\star}^{\sigma}{}_{\sigma}$ can be calculated from Eq. (53) by taking $\mu=\bar{\sigma}=\sigma$ and summing over $\sigma$ accordingly. The matrix $\boldsymbol{\rm\mathcal{X}}$ is given by Eq. (25), and the matrix $\tilde{\boldsymbol{\rm\mathcal{X}}}$ is given by Eq. (18). The index $\star$ denotes that the corresponding quantities are evaluated at the RR (i.e., at $\tilde{\boldsymbol{\rm\varrho}}=0$) and thus are independent of the transverse coordinates.

III Numerical algorithm

The code PARADE is based on the model described in Sec. II. The RR starting point $\boldsymbol{\rm X}_{0}$, the orientation of the RR initial wave vector, the frequency $f\doteq\omega/2\pi$, the initial profile of the complex amplitude $a_{0}$, the magnetic field profile $\boldsymbol{\rm B}(\boldsymbol{\rm x})$, and the plasma density profile $n(\boldsymbol{\rm x})$ are considered pre-determined. Currently, they are entered as analytic functions, but it is also possible to use experimental data or profiles generated numerically. The value of $|\boldsymbol{\rm K}_{0}|$ is determined by the local dispersion relation at $\boldsymbol{\rm X}_{0}$. The collisionless cold-electron-plasma model is assumed for the dielectric tensor book:stix here, but the code can also accept a general dispersion tensor $\boldsymbol{\rm D}$. In particular, warm-plasma simulations will be reported separately.

Figure 1 outlines the code operation. At each point on the RR, starting from the initial point $\boldsymbol{\rm X}_{0}$, the dispersion tensor $\boldsymbol{\rm D}=\boldsymbol{\rm D}_{H}+i\boldsymbol{\rm D}_{A}$ is calculated. Then, the eigenvalues of $\boldsymbol{\rm D}_{H}$ are calculated, and the eigenvalue $\Lambda_{\star}$ is found that has the least absolute value. This $\Lambda_{\star}$ is then numerically differentiated and used as the Hamiltonian in Eqs. (15) to propagate the RR using the fourth-order Runge–Kutta method. Based on the knowledge of the RR trajectory $\boldsymbol{\rm X}(\zeta)$, the ray-based coordinates $\tilde{x}^{\mu}$ are constructed by choosing the basis vectors $\tilde{\boldsymbol{\rm e}}_{\star\mu}$, which is done as follows.

We adopt $\tilde{\boldsymbol{\rm e}}_{\star\zeta}$ to be the unit vector along the RR group velocity. The other two vectors $\tilde{\boldsymbol{\rm e}}_{\star\sigma}$ (where $\sigma=1,2$) can, in principle, be chosen arbitrarily as long as they satisfy Eqs. (20) and the resulting metric is smooth. The traditional Frenet–Serret frame can be used but sometimes fails the latter requirement. Instead, we choose

[TABLE]

at each given time step, where the vectors marked with ∙ are those from the previous time step. (The prime denotes the derivative with respect to $\zeta$, so $|\boldsymbol{\rm X}^{\prime}|=1$.) For the very first time step, $\tilde{\boldsymbol{\rm e}}_{\star 1}^{\bullet}$ and $\tilde{\boldsymbol{\rm e}}_{\star 2}^{\bullet}$ can be chosen arbitrarily as long as they are orthogonal to each other and to the initial group velocity on the RR. This ensures that the basis vectors $\tilde{\boldsymbol{\rm e}}_{\star\mu}$ transform continuously along the RR and result in a smooth ray-based metric.

As the next step, the polarization eigenvector $\boldsymbol{\rm\Xi}_{\star}$ is calculated as the eigenvector of $\boldsymbol{\rm D}_{H}$ corresponding to the eigenvalue $\Lambda_{\star}$; also, its complex conjugate $\boldsymbol{\rm\Xi}^{+}_{\star}$ is determined. After these vectors are numerically differentiated, the coefficients in Eq. (59) are calculated as outlined in Sec. II.4.4. The field $a$ is mapped to $\phi$ using Eq. (58), and $\phi$ is differentiated with respect to $\tilde{\varrho}^{\sigma}$ at fixed $\zeta$ using the finite-difference approach. Then, $\phi$ is propagated by solving Eq. (59) using the fourth-order Runge–Kutta method. Finally, $\phi$ is mapped back to $a$, and the three-dimensional profile $a(\boldsymbol{\rm x})$ is yielded as a data array.

IV Simulation results

In this section, test simulation results are reported. These results were obtained by applying PARADE to sample problems that fit the assumptions of the model described above. (Simulations that involve mode conversion are presented in Paper III.) 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.

IV.1 Beams in vacuum

As the first example, we simulated the propagation of a wave beam in vacuum. For the initial beam profile, we assumed a Gaussian profile

[TABLE]

where

[TABLE]

and $k=2\pi f/c$, with the focal lengths $\mathcal{Z}_{1}=3.0$ m and $\mathcal{Z}_{2}=4.0$ m, the waist sizes $w_{0,1}=2.0$ cm and $w_{0,2}=3.0$ cm, and the wave frequency $f=154.0$ GHz, which corresponds to the vacuum wavelength $\lambda_{0}\approx 2$ mm. The results are shown in Fig. 2. At $\zeta>0$, the beam width $w_{\sigma}$ is determined numerically as the distance on which $|a|$ drops one $e$-fold from its maximum along the $\tilde{\varrho}^{\sigma}$-axis. As seen from Fig. 2(d), PARADE demonstrates good agreement with the analytical formula from Gaussian-beam optics book:yariv .

IV.2 Straight beams in cold plasma

As the second example, we simulated the propagation of waves in collisionless magnetized cold electron plasma. From now on, we introduce the standard notation $\{x,y,z\}$ for the laboratory coordinates and $\boldsymbol{\rm e}_{\alpha}$ for the unit vectors along the corresponding axes. In particular, the $x$ axis is chosen to be along the initial direction of the wave vector. We assume a slab geometry with

[TABLE]

Here, $n_{0}=2.0\times 10^{19}$ m*-3*, $L_{n}=1.0$ m, and $B_{0}=2.0$ T. The initial field has the Gaussian profile (63) with the focal lengths $\mathcal{Z}_{1}=\mathcal{Z}_{2}=3.0$ m, the waist sizes $w_{0,1}=w_{0,2}=2.0$ cm, and the wave frequency $f=154.0$ GHz. The simulation results are presented in Fig. 3. Since the waves propagate along the gradient of $n$, the RR trajectory is straight and does not depend on the polarization, just like in the vacuum case. However, the intensity profiles of the waves with different polarizations, which are the O and X waves here, are different from each other and from those in vacuum. The X wave is affected by the plasma more significantly. This is explained by the fact that the wave frequency $f$ is closer to the X-wave cutoff frequency $f_{\rm rc}$ than to the O-wave cutoff frequency $f_{\rm pc}$. In addition, although $\mathcal{Z}_{1}$ and $\mathcal{Z}_{2}$ are assumed to be equal initially [dashed line in Fig. 3(b) and (c)], they start to differ within the plasma [solid lines in Fig. 3(b) and (c)] due to the magnetic field. Also, while the focal lengths are different from those in vacuum, the waist sizes are, in fact, the same. This is because $\boldsymbol{\rm k}$ is uniform in the transverse plane for the lack of the density gradient in that plane. Simulations with transverse gradients are presented in Sec. IV.3. Then, we have also tested PARADE for the energy conservation [Eq. (57)], which is expected here because waves in collisionless cold plasma have $\Gamma=0$. While the wave intensity and group velocity change considerably along the RR path, the wave-energy flux is conserved by the code with high accuracy, as seen in Figs. 3(d) and (e).

IV.3 Bended beams in cold plasma

As another example, we applied PARADE to simulate the propagation of bended wave beams launched at generic angles relative to the density gradient and the magnetic field. The assumed density and magnetic-field profiles are as follows:

[TABLE]

where $n_{0}=5.0\times 10^{19}$ m*-3*, $x_{c}=2.0$ m, $L_{n}=2.0$ m, $B_{0}=2.0$ T, and $L_{b}=2.0$ m. The RR starting point is $\{x,y,z\}=\{0.1,0.0,3.0\}$ m, and the target point used to fix the orientation of the RR initial wave vector is $\{x,y,z\}=\{3.1,0.0,0.0\}$ m. For the initial profile, we adopted the Gaussian profile (63) with the focal lengths $\mathcal{Z}_{1}=\mathcal{Z}_{2}=1.5$ m, the waist sizes $w_{0,1}=w_{0,2}=2.7$ cm, and the wave frequency $f=77.0$ GHz (or $\lambda_{0}\approx 4$ mm). A simulation of the X-mode propagation is shown in Fig. 4(a). Since PARADE accounts for diffraction, the beam retains a nonzero width on the whole trajectory. This is different from conventional ray tracing, where diffraction is neglected and each ray propagates independently. A comparison is presented in Fig. 4(b). The difference is essential for resonant-plasma heating and current drive, where the wave energy deposition must be calculated accurately, as also discussed in Refs. ref:mazzucato89 ; ref:nowak93 ; ref:peeters96 ; ref:farina07 ; ref:pereverzev98 ; ref:poli01b ; ref:poli01 ; ref:poli18 ; ref:balakin08b ; ref:balakin07a ; ref:balakin07b .

The corresponding evolution of the key frequencies along the RR is presented in Fig. 5(a), and the corresponding beam widths are presented in Fig. 5(b). For a reference, Fig. 5(c) shows the same widths calculated with neglected $\smash{\tilde{\boldsymbol{\rm\varrho}}}$-dependence of $\boldsymbol{\rm k}$. (In the latter case, the waist sizes are almost the same as for vacuum waves.) These results show that capturing the transverse inhomogeneity of $\boldsymbol{\rm k}$ in practical simulations can be important.

We also note that in our cold-plasma PARADE simulations, the beams remain Gaussian with high accuracy (Fig. 6). However, it is generally not the case when resonant dissipation is included, for such dissipation can be strongly inhomogeneous in the beam cross-section. PARADE simulations of dissipative beams in hot plasma will be reported in a separate publication.

V Conclusions

We report the first quasioptical simulations of wave beams in cold magnetized plasma using a new code PARADE. This code does not assume any particular beam structure but instead solves a parabolic equation for the wave envelope; hence, it is particularly well-suited to modeling strongly inhomogeneous resonant dissipation in the future. The general theoretical model underlying PARADE was largely derived in the previous paper foot:paper1 and describes multi-mode beams that can experience mode conversion. Here, we consider a simplified version of this theory in which mode conversion is not considered and a beam can be described as a scalar field. We adjust the theory to numerical modeling and simulate the propagation of wave beams in vacuum and inhomogeneous magnetized plasma. Our numerical results show good agreement with Gaussian-beam optics and surpass conventional ray tracing. Quasioptical simulations of mode-converting beams are reported in the next, third paper of this series.

VI 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.