Quasioptical modeling of wave beams with and without mode conversion: I. Basic theory

TL;DR

This paper introduces a general quasioptical theory for mode-converting electromagnetic beams in plasma, deriving a simplified parabolic equation for the wave envelope that accounts for mode conversion and polarization effects.

Contribution

It develops a comprehensive theoretical framework and numerical approach for modeling mode-converting electromagnetic beams in plasma, including the derivation of a quasioptical equation.

Findings

Derived a second-order approximate operator for wave envelopes.

Simplified the operator to a parabolic differential equation.

Applicable to scalar and vector mode-converting beams.

Abstract

This work opens a series of papers where we develop a general quasioptical theory for mode-converting electromagnetic beams in plasma and implement it in a numerical algorithm. Here, the basic theory is introduced. We consider a general quasimonochromatic multi-component wave in a weakly inhomogeneous linear medium with no sources. For any given dispersion operator that governs the wave field, we explicitly calculate the approximate operator that governs the wave envelope $\psi$ to the second order in the geometrical-optics parameter. Then, we further simplify this envelope operator by assuming that the gradient of $\psi$ transverse to the local group velocity is much larger than the corresponding parallel gradient. This leads to a parabolic differential equation for $\psi$ ("quasioptical equation") in the basis of the geometrical-optics polarization vectors. Scalar and mode-converting vector beams are described on the same footing. We also explain how to apply this model to electromagnetic waves in general. In the next papers of this series, we report successful quasioptical modeling of radiofrequency wave beams in magnetized plasma based on this theory.

Equations

$$\hat{D}\Psi =0,$$

$$ϵ\doteq \lambda /L\ll 1$$

$$\Psi =\hat{\mathcal{U}}\psi ,\hat{\mathcal{U}}\doteq e^{i\theta (\mathrm{x})}.$$

$$\mathrm{k}(\mathrm{x})\doteq \nabla \theta (\mathrm{x})$$

$$\hat{\mathcal{D}}\psi =0,$$

$$\hat{\mathcal{D}}=e^{-i\theta }\hat{D}{e}^{i\theta }.$$

$$\frac{\mathrm{d}{x}^{\alpha }}{\mathrm{d}\tau }=\frac{\partial H}{\partial k_{\alpha }},\frac{\mathrm{d}{k}_{\alpha }}{\mathrm{d}\tau }=-\frac{\partial H}{\partial x^{\alpha }}$$

$$H(\mathrm{x},\mathrm{k}(\mathrm{x}))=\text{const}.$$

$$(\hat{A}\Psi )(\mathrm{x})\doteq \int {\mathrm{d}}^n{x}^{\prime }\sqrt{g_{\diamond }({\mathrm{x}}^{\prime })}\mathsf{A}(\mathrm{x},{\mathrm{x}}^{\prime })\Psi ({\mathrm{x}}^{\prime }).$$

$$⟨\Psi ⟩\Phi \doteq \int {\mathrm{d}}^{n}x\sqrt{g_{\diamond }(\mathrm{x})}{\Psi }^{\ast }(\mathrm{x})\Phi (\mathrm{x}).$$

$$\hat{\mathrm{x}}=\{{\hat{x}}^0,{\hat{x}}^1,\dots ,{\hat{x}}^{n-1}\},$$

$$\hat{\mathrm{p}}=\{{\hat{p}}_0,{\hat{p}}_1,\dots ,{\hat{p}}_{n-1}\}$$

$${\hat{x}}^{\mu }\Psi =x^{\mu }\Psi ,$$

$${\hat{p}}_{\mu }\Psi =-i{g}_{\diamond }^{-1/4}{\partial }_{\mu }(g_{\diamond }^{1/4}\Psi ).$$

$${\hat{p}}_{\nu }=-i{\partial }_{\nu }-i{q}_{\nu }(\mathrm{x}),q_{\nu }\doteq \frac{1}{4}{\partial }_{\nu }(\ln g_{\diamond }),$$

$$\hat{\mathrm{x}}|\mathrm{x}⟩=\mathrm{x}|\mathrm{x}⟩,\hat{\mathrm{p}}|\mathrm{p}⟩=\mathrm{p}|\mathrm{p}⟩.$$

$$⟨{\mathrm{x}}_1⟩{\mathrm{x}}_2=\mathfrak{G}({\mathrm{x}}_1,{\mathrm{x}}_2)\delta ({\mathrm{x}}_1-{\mathrm{x}}_2),$$

$$⟨{\mathrm{p}}_1⟩{\mathrm{p}}_2=\bar{\mathfrak{G}}({\mathrm{p}}_1,{\mathrm{p}}_2)\delta ({\mathrm{p}}_1-{\mathrm{p}}_2).$$

$$\mathfrak{G}({\mathrm{x}}_1,{\mathrm{x}}_2)\doteq [g_{\diamond }({\mathrm{x}}_1)g_{\diamond }({\mathrm{x}}_2){]}^{-1/4},$$

$$\bar{\mathfrak{G}}({\mathrm{p}}_1,{\mathrm{p}}_2)\doteq [{\bar{g}}_{\diamond }({\mathrm{p}}_1){\bar{g}}_{\diamond }({\mathrm{p}}_2){]}^{-1/4}.$$

$$\Psi (\mathrm{x})=⟨\mathrm{x}⟩\Psi ,\psi (\mathrm{x})=⟨\mathrm{x}⟩\psi .$$

$$⟨\mathrm{x}⟩\mathrm{p}=⟨\mathrm{p}⟩{\mathrm{x}}^{\ast }=\frac{\exp (i\mathrm{p}\cdot \mathrm{x}),}{(2\pi {)}^{n/2}[g_{\diamond }(\mathrm{x}){\bar{g}}_{\diamond }(\mathrm{p}){]}^{1/4}}.$$

$${\hat{\Delta }}_z\doteq \int {\mathrm{d}}^{n}sG(\mathrm{x},\mathrm{s})|\mathrm{x}-\mathrm{s}/2⟩⟨\mathrm{x}+\mathrm{s}/2|e^{-i\mathrm{p}\cdot \mathrm{s}},$$

$$G(\mathrm{x},\mathrm{s})\doteq [g_{\diamond }(\mathrm{x}-\mathrm{s}/2)g_{\diamond }(\mathrm{x}+\mathrm{s}/2){]}^{1/4}.$$

$$A(\mathrm{x},\mathrm{p})\doteq \int {\mathrm{d}}^{n}sG(\mathrm{x},\mathrm{s})⟨\mathrm{x}+\mathrm{s}/2⟩\hat{A}\mathrm{x}-\mathrm{s}/2e^{-i\mathrm{p}\cdot \mathrm{s}}.$$

$${\hat{A}}_H\doteq \frac{1}{2}(\hat{A}+{\hat{A}}^{†}),{\hat{A}}_A\doteq \frac{1}{2i}(\hat{A}-{\hat{A}}^{†}).$$

$$\hat{A}=\frac{1}{(2\pi {)}^n}\int {\mathrm{d}}^{n}x{\mathrm{d}}^{n}pA(\mathrm{x},\mathrm{p}){\hat{\Delta }}_z.$$

$$f(\hat{\mathrm{x}})\iff f(\mathrm{x}),f(\hat{\mathrm{p}})\iff f(\mathrm{p}).$$

$$f(\hat{\mathrm{x}}){\hat{p}}_{\alpha }\iff p_{\alpha }f(\mathrm{x})+\frac{i}{2}{\partial }_{\alpha }f(\mathrm{x}),$$

$${\hat{p}}_{\alpha }f(\hat{\mathrm{x}})\iff p_{\alpha }f(\mathrm{x})-\frac{i}{2}{\partial }_{\alpha }f(\mathrm{x}),$$

$${\hat{p}}_{\alpha }f(\hat{\mathrm{x}}){\hat{p}}_{\beta }\iff$$

$$+\frac{i}{2}{p}_{\alpha }{\partial }_{\beta }f(\mathrm{x})-\frac{i}{2}{p}_{\beta }{\partial }_{\alpha }f(\mathrm{x}).$$

$$\hat{D}|\Psi ⟩=0.$$

$$\hat{\mathcal{D}}|\psi ⟩=0,$$

Full text

I Introduction

I.1 Motivation

Describing the propagation of waves in inhomogeneous media is a classic problem with a long history book:tracy ; book:kravtsov ; book:whitham . It is particularly important in fusion research, where quasistationary beams of electromagnetic (EM) radiation are commonly used for many purposes and need to be modeled with fidelity book:stix . Full-wave modeling, which involves solving the complete Maxwell’s equations, can be impractical at short (cm and mm) wavelengths, especially when multi-dimensional simulations with complex geometries are required, such as those of tokamak and stellarator plasmas. Hence, reduced methods have been widely used in practice. These methods are rooted in geometrical optics (GO) book:tracy and include conventional ray tracing ref:bonoli14 ; tex:genray ; ref:marushchenko14 ; ref:tsujimura15 , complex ray tracing ref:mazzucato89 ; ref:farina07 , beam tracing ref:pereverzev92 ; ref:pereverzev98 ; ref:poli01 ; ref:poli18 , and variations of thereof foot:thereof . There are also other “quasioptical” models, such as in Refs. ref:balakin07a ; ref:balakin07b ; foot:balakin ; foot:phillips86 , that resolve the evolution of the beam transverse structure without adopting any particular ansatz for the intensity profile. Still, they assume that only one branch of the dispersion relation is excited in each given case foot:raycon , i.e., mode conversion does not occur foot:mc . As a result, today’s simulations of mode-converting beams mainly rely on full-wave codes foot:phillips86 and thus have to compromise the fidelity by reducing the number of the dimensions resolved ref:tsujimura15 .

However, simulating mode conversion does not actually require the full-wave approach; in fact, it can be done within the quasioptical approach too, if the latter is generalized properly. Developing the corresponding algorithms would facilitate not only fusion applications. For example, advanced quasioptical models could find use in optics and general relativity, where mode-converting beams are also possible and have been attracting much attention lately ref:bliokh15 ; foot:marius . Hence, it is potentially beneficial to approach the problem of quasioptical modeling as a general-physics problem, without restricting it to fusion applications.

I.2 General idea

Mode conversion can be described quite generally, particularly without restricting it to narrow regions in space, within “extended geometrical optics” (XGO), which was developed recently phd:ruiz17 ; my:covar ; my:qdirac ; my:qdiel ; my:qdirpond . XGO is a theory that calculates the leading-order correction $U$ to the GO dispersion operator of a general vector wave and shows my:qdirac that this correction is analogous to (and a generalization of) the Stern–Gerlach Hamiltonian of a quantum spin-$1/2$ electron. Accordingly, $U$ is responsible for two effects simultaneously: (i) it modifies the ray equations just like spin–orbital interactions affect the electron motion, and (ii) it also governs mode conversion, which appears as a direct analog of spin up–down transitions my:xo . An examination of the quasioptical algorithms such as those in Refs. ref:pereverzev98 ; ref:balakin07a shows that they already involve calculations of terms similar to $U$. Hence, adding the mode-conversion capability to quasioptical codes should not be burdensome and is not expected to slow down the codes considerably. However, formulating the corresponding theory is easier to do using the abstract quantumlike formalism of XGO. One only needs to upgrade the existing XGO by adding diffraction, for which it also helps to introduce more general coordinates with curved metric foot:covar .

In this series of papers (Papers I-III), we propose such an upgrade of XGO and apply it to numerical simulations. Our goal is to develop a general modular framework which later could be applied to a broad variety of problems both in plasma physics and beyond. Within this framework to be presented, one can assume general dispersion, diffraction, and polarization effects, including mode conversion of not just two but arbitrarily many resonant waves. Also importantly, our formulation below is not restricted to EM waves, since we do not specify the dispersion operator in the governing wave equation.

Our series of papers is organized as follows. In Paper I, we introduce the basic theory of waves that diffract and mode-convert simultaneously. In Paper II tex:mypaper2 and Paper III tex:mypaper3 , we apply this theory to perform quasioptical modeling of radiofrequency-wave beams in magnetized plasma as an example. In particular, we consider applications to mode conversion caused by magnetic shear in edge plasma my:xo , which, for example, is a known problem ref:notake05 ; ref:kubo15 in the Large Helical Device ref:iiyoshi99 ; ref:osakabe17 .

I.3 Outline

In this first paper of our series, we consider an arbitrary quasimonochromatic multi-component wave in a weakly inhomogeneous linear medium. Supposing that $\smash{\widehat{D}}$ is some dispersion operator governing the wave dynamics, we simplify $\smash{\widehat{D}}$ and obtain an approximate operator that governs the wave envelope $\psi$. Then, we derive a parabolic differential equation for $\psi$ (“quasioptical equation”) by assuming that the gradient of the wave envelope transverse to the local group velocity is much larger than the corresponding parallel gradient. The resulting theory applies to both scalar and mode-converting vector beams. At the end of the paper, we also discuss how this model can be applied to EM waves in particular. However, readers who are mainly interested in simulations as opposed to the general theory are encouraged to proceed straight to Paper II, where our key equations are overviewed in a simplified form and without derivations.

This paper is organized as follows. In Sec. II, we introduce the general problem. In Sec. III, we formalize the concept of the envelope dispersion operator. In Sec. IV, we derive an approximation for the envelope dispersion operator for scalar waves, and we also derive its quasioptical approximation. In Sec. V, we extend this model to vector waves. In Sec. VI, we explain how to apply the resulting theory to EM waves in particular. In Sec. VII, we present our main conclusions. In Appendix A, we summarize some of our notations. In Appendix B, some auxiliary calculations are presented. Our paper also contains Supplementary Material foot:supp . There, we overview the Weyl calculus on a curved configuration space, which is used in this work.

II General problem

Consider a wave propagating on an $n$-dimensional configuration space $M^{n}$ with coordinates ${\boldsymbol{\rm x}}\equiv\{x^{0},x^{1},\ldots,x^{n-1}\}$ and some general metric tensor $g({\boldsymbol{\rm x}})$. For simplicity, assume that $M^{n}$ is diffeomorphic to $\mathbb{R}^{n}$, i.e., the $n$-dimensional Euclidean space or pseudo-Euclidean space with the same metric signature as $M^{n}$. (For more information on why the diffeomorphism with $\mathbb{R}^{n}$ is needed, see Sec. IV.1.2.) Suppose that the wave field $\Psi\equiv\Psi({\boldsymbol{\rm x}})$, which may have multiple components, is governed by a linear equation with no source terms,

[TABLE]

where $\smash{\widehat{D}}$ is a differential or, most generally, integral dispersion operator. (For vector waves, $\smash{\widehat{D}}$ is a matrix whose elements are operators; see Sec. V.) We shall assume that the GO parameter $\epsilon$ is small; namely,

[TABLE]

(the symbol $\doteq$ denotes definitions), where $\lambda$ is the characteristic wave period, or wavelength, and $L$ is the least characteristic scale among those of the wave envelope and of the medium, including the metric foot:GOappr . Below, we propose a systematic reduction of Eq. (1) using the smallness of $\epsilon$ and eventually obtain a quasioptical model based on this equation. The idea of quasioptical modeling will be formalized later (Secs. IV.5 and V.6).

As a side note, we emphasize that the theory to be developed does not describe wave transformations near cutoffs, where the GO parameter (2) is not small. It is possible to waive this limitation, but formulating such generalized theory is left to future publications.

III Envelope dispersion operator

As the first step, let us introduce a unitary variable transformation

[TABLE]

(Other $\smash{\widehat{\mathcal{U}}}$ may also be justified in some cases, e.g., for dealing with caustics or quasiperiodic media my:qponder , but we shall not consider this possibility in the present work.) The phase $\theta$, which we call the “reference phase”, serves as a gauge potential. It is a real function such that

[TABLE]

is the wave vector identical or close to that predicted by the GO approximation. This implies that the envelope $\psi$ and also ${\boldsymbol{\rm k}}$ are slow functions [and that the wavelength entering Eq. (2) is $\lambda\approx 2\pi/k$]. Then, Eq. (1) becomes

[TABLE]

where the “envelope dispersion operator” is $\smash{\widehat{\mathcal{D}}}\doteq\smash{\widehat{\mathcal{U}}}^{\dagger}\smash{\widehat{D}}\smash{\widehat{\mathcal{U}}}$ (the dagger denotes the adjoint, as usual), or more explicitly,

[TABLE]

The reference phase $\theta$ is treated as a prescribed function. As will become clear later, knowing $\theta$ per se is not actually needed for our purposes; instead, it is ${\boldsymbol{\rm k}}$ that matters. The latter can be calculated on some “reference rays” using the conventional ray equations book:tracy

[TABLE]

($\tau$ is any parameter along the ray), which also lead to

[TABLE]

Note that, in general, $\theta$ and ${\boldsymbol{\rm k}}$ are defined uniquely only to the leading order, so there exists some freedom in choosing reference rays and their Hamiltonian $H$. This means that more than one $\smash{\widehat{\mathcal{D}}}$ is possible. Still, envelope equations that (slightly) differ in the choice of ${\boldsymbol{\rm k}}$ are equivalent in the sense that the total field $\Psi$ that they describe is the same by construction. We shall discuss this in more detail in Secs. V.5 and V.7.

Also note that our approach equally applies to stationary and nonstationary waves. In the case of a stationary wave, we assume that ${\boldsymbol{\rm x}}$ is a coordinate in physical space (“spatial problem”), and the wave frequency $\omega$ serves as a constant parameter. In the case of a nonstationary wave, we assume that ${\boldsymbol{\rm x}}$ is a coordinate in spacetime (“spacetime problem”), and then $\omega$ is a part of ${\boldsymbol{\rm k}}$. In spacetime problems, we assume coordinates such that $x^{0}=ct$, where $c$ is the speed of light, $t$ is time, and the metric signature is $(-,+,+,\ldots)$; then, $k_{0}=-\omega/c$ (in case of the Minkowski metric, this implies that $k^{0}=\omega/c$), so $\omega=-\partial_{t}\theta$, as usual. In other respects, spatial and spacetime problems are described on the same footing and will be distinguished only in Sec. VI.

Having defined this terminology, we shall now discuss how $\smash{\widehat{\mathcal{D}}}$ can be expanded in $\epsilon$ asymptotically for any $\smash{\widehat{D}}$.

IV Scalar waves

Our asymptotic theory of mode-converting beams is based on the well-known phase-space methods in quantum and classical wave theory. (The key papers include, but are not limited to, Refs. foot:bookweyl ; ref:wigner32 ; ref:moyal49 ; ref:littlejohn86 ; ref:friedland87b ; ref:mcdonald88 ; ref:littlejohn91 ; for recent overviews, see Refs. book:tracy ; phd:ruiz17 .) In order to present our framework and notation in a self-contained manner, we restate some basics in Sec. IV.1 and Supplementary Material foot:supp . Also note that our calculations in Sec. IV.2 are similar in spirit to those in LABEL:ref:mcdonald88, except we perform a higher-order expansion and assume a more general metric. We shall also comment throughout the paper on how our equations reproduce other known results in the corresponding limits.

IV.1 Weyl calculus

IV.1.1 Basic definitions

Until Sec. V, we shall assume that $\Psi$ is a scalar function. Any given operator $\smash{\widehat{A}}$ acting on it maps $\Psi$ to a new scalar function $\smash{\widehat{A}}\Psi$ that can be expressed as follows:

[TABLE]

Here the integral is taken over $\mathbb{R}^{n}$ (and so are all integrals below, up to dimension), $g_{\diamond}\doteq|\det g\,|$, and $\mathsf{A}$ is some kernel function that determines $\smash{\widehat{A}}$. Consider also a family of all unitary operators $\smash{\widehat{A}}_{u}$, which is a subset of all possible $\smash{\widehat{A}}$. For a given $\Psi$, all image functions $\smash{\widehat{A}}_{u}\Psi$ are mutually equivalent up to an isomorphism, so their family $\smash{\{\smash{\widehat{A}}_{u}\Psi\}}$ can be viewed as a single object, a “state vector” $\ket{\Psi}$, which belongs to a Hilbert space $\mathcal{H}_{1}^{n}$ with inner product

[TABLE]

Then, Eq. (9) can be viewed as the “${\boldsymbol{\rm x}}$ representation” of $\smash{\widehat{A}}$, while the operator itself can be understood more generally as a transformation of $\ket{\Psi}$, i.e., of the whole family $\smash{\{\smash{\widehat{A}}_{u}\Psi\}}$. Using this invariant notation, one can formulate a machinery, called the Weyl calculus foot:wc , that allows efficient asymptotic approximation of operators using the smallness of $\epsilon$. Below, we overview the key theorems of the Weyl calculus that are used in our paper. Readers who are interested in details and proofs of these theorems can find them in the Supplementary Material foot:supp .

IV.1.2 Coordinate and momentum operators

We start by defining the coordinate and momentum (wave-vector) operators

[TABLE]

such that the ${\boldsymbol{\rm x}}$ representations of $\smash{\widehat{x}}^{\mu}$ and $\smash{\widehat{p}}_{\mu}$ be as follows:

[TABLE]

Here $\partial_{\mu}\doteq\partial/\partial x^{\mu},$ $x^{\mu}$ and $p_{\mu}$ are the corresponding eigenvalues, and the factors $\smash{g_{\diamond}^{\pm 1/4}}$ are introduced to keep $\smash{\widehat{{\boldsymbol{\rm p}}}}$ self-adjoint under the inner product (10) ref:dewitt52 . (This would not be the case if $M^{n}$ were not diffeomorphic to $\mathbb{R}^{n}$ ref:domingos84 .) Since

[TABLE]

and $q_{\nu}$ commutes with $\smash{\widehat{x}}^{\mu}$, one arrives at the usual commutation relation $[\smash{\widehat{x}}^{\mu},\smash{\widehat{p}}_{\nu}]=i\delta^{\mu}_{\nu}$.

Let us consider the eigenvectors $\ket{{\boldsymbol{\rm x}}}$ and $\ket{{\boldsymbol{\rm p}}}$ of the coordinate and momentum operators, which are defined as

[TABLE]

Since the operators are self-adjoint, these eigenvectors can be chosen as mutually orthogonal, and we shall assume the following normalization:

[TABLE]

Here we introduced

[TABLE]

The function $\bar{g}_{\diamond}$ can be chosen arbitrarily as long as it is kept positive. It plays a role similar to that of $g_{\diamond}$ in the Weyl calculus, but note that this is just a normalization factor, and we introduced it only to maintain the symmetry between $\smash{\widehat{{\boldsymbol{\rm x}}}}$ and $\smash{\widehat{{\boldsymbol{\rm p}}}}$. One can show then foot:supp that our original function $\Psi$ and the envelope $\psi$ can be expressed through the corresponding state vectors as

[TABLE]

The ${\boldsymbol{\rm p}}$ representations of these state vectors are introduced similarly as $\braket{{\boldsymbol{\rm p}}}{\Psi}$ and $\braket{{\boldsymbol{\rm p}}}{\psi}$. One can also show foot:supp that $\braket{{\boldsymbol{\rm x}}_{1}}{\smash{\widehat{A}}}{{\boldsymbol{\rm x}}_{2}}=\mathsf{A}({\boldsymbol{\rm x}}_{1},{\boldsymbol{\rm x}}_{2})$, and

[TABLE]

Here, ${\boldsymbol{\rm p}}\cdot{\boldsymbol{\rm x}}\doteq p_{\mu}x^{\mu}$, and summation over repeated indices is assumed here and further.

IV.1.3 Wigner–Weyl transform

The set of all eigenvalues of the coordinate and momentum operators form a $2n$-dimensional “phase space” $z\equiv({\boldsymbol{\rm x}},{\boldsymbol{\rm p}})$. For every given $z$, we introduce the so-called Wigner operator $\smash{\widehat{\Delta}}_{z}$, which is self-adjoint and defined as

[TABLE]

(Note that $G=1$ if the ${\boldsymbol{\rm x}}$ space is Euclidean or pseudo-Euclidean.) Using $\smash{\widehat{\Delta}}_{z}$, we define the Wigner–Weyl transform $\mathscr{W}_{z}:\smash{\widehat{A}}\mapsto A$, which maps a given operator $\smash{\widehat{A}}$ on $\mathcal{H}_{1}^{n}$ to a function $A$ (“Weyl symbol”) on the $z$ space. Specifically, the Weyl symbol of a given operator $\smash{\widehat{A}}$ is $A({\boldsymbol{\rm x}},{\boldsymbol{\rm p}})\doteq\text{tr}(\smash{\widehat{\Delta}}_{z}\smash{\widehat{A}})$ (“tr” stands for trace); i.e.,

[TABLE]

As can be seen easily from this definition, if a given operator is self-adjoint, then its Weyl symbol is real. This also leads to the following corollary. Consider splitting a given $\smash{\widehat{A}}$ as $\smash{\widehat{A}}=\smash{\widehat{A}}_{H}+i\smash{\widehat{A}}_{A}$, where the subscripts denote the Hermitian and anti-Hermitian parts,

[TABLE]

Both $\smash{\widehat{A}}_{H}$ and $\smash{\widehat{A}}_{A}$ (not to be confused with $i\smash{\widehat{A}}_{A}$) are self-adjoint by definition. Thus, the corresponding Weyl images $A_{H}$ and $A_{A}$ are real.

We also define the inverse Wigner–Weyl transform $\mathscr{W}^{-1}:A\mapsto\smash{\widehat{A}}$, which maps a given function $A$ to the corresponding operator $\smash{\widehat{A}}$ via

[TABLE]

The direct and inverse transforms set the “Weyl correspondence” between operators and functions on the $({\boldsymbol{\rm x}},{\boldsymbol{\rm p}})$ space, $\smash{\widehat{A}}\Leftrightarrow A({\boldsymbol{\rm x}},{\boldsymbol{\rm p}})$. As can be checked by a direct calculation, for any function $f$, one has

[TABLE]

However, the Weyl symbols of operators that cannot be represented as $f_{1}(\smash{\widehat{{\boldsymbol{\rm x}}}})+f_{2}(\smash{\widehat{{\boldsymbol{\rm p}}}})$ are generally more complicated. In particular, one can show that foot:supp

[TABLE]

and also foot:supp

[TABLE]

Overall, the Weyl symbol of an operator that is any given combination $f(\smash{\widehat{{\boldsymbol{\rm x}}}},\smash{\widehat{{\boldsymbol{\rm p}}}})$ of $\smash{\widehat{{\boldsymbol{\rm x}}}}$ and $\smash{\widehat{{\boldsymbol{\rm p}}}}$ approaches $f({\boldsymbol{\rm x}},{\boldsymbol{\rm p}})$ in the GO limit, when $[\smash{\widehat{x}}^{\mu},\smash{\widehat{p}}_{\nu}]$ is negligible. However, in general, $f(\smash{\widehat{{\boldsymbol{\rm x}}}},\smash{\widehat{{\boldsymbol{\rm p}}}})$ does not map simply to $f({\boldsymbol{\rm x}},{\boldsymbol{\rm p}})$.

IV.2 Approximate $\boldsymbol{\smash{\widehat{\mathcal{D}}}}$

Using the notation introduced in Sec. IV.1, one can express Eq. (1) in the following invariant form:

[TABLE]

Accordingly, Eq. (5) becomes

[TABLE]

where the invariant form of the envelope dispersion operator [Eq. (6)] is as follows:

[TABLE]

We shall now approximate $\smash{\widehat{\mathcal{D}}}$ in three steps: (i) we map the right-hand side of Eq. (34) onto a function, or symbol, using the Wigner–Weyl transform; (ii) we approximate this symbol using the smallness of $\epsilon$ [Eq. (2)]; and (iii) we produce an operator out of the approximated symbol using the inverse Wigner–Weyl transform.

The first two steps can be done by expressing the Weyl symbol of $\smash{\widehat{\mathcal{D}}}$ as $\mathcal{D}({\boldsymbol{\rm x}},{\boldsymbol{\rm p}})=e^{-i\theta({\boldsymbol{\rm x}})}\star D({\boldsymbol{\rm x}},{\boldsymbol{\rm p}})\star e^{i\theta({\boldsymbol{\rm x}})}$, where $D$ is the Weyl symbol of $\smash{\widehat{D}}$ and $\star$ is the Moyal product foot:supp . By expanding $\star$ in the GO parameter, one can then obtain an approximate symbol of $\smash{\widehat{\mathcal{D}}}$ formally. But here we adopt a different approach, which is more transparent. We start by expressing $\smash{\widehat{D}}$ through $D$ explicitly,

[TABLE]

where $G$ is a metric factor given by Eq. (24). Using the fact that $\theta(\smash{\widehat{{\boldsymbol{\rm x}}}})\ket{{\boldsymbol{\rm x}}\pm{\boldsymbol{\rm s}}/2}=\theta({\boldsymbol{\rm x}}\pm{\boldsymbol{\rm s}}/2)\ket{{\boldsymbol{\rm x}}\pm{\boldsymbol{\rm s}}/2}$, one can rewrite Eq. (34) as follows:

[TABLE]

Consider a formal Taylor expansion of the reference phase $\theta$ in ${\boldsymbol{\rm s}}$:

[TABLE]

This gives

[TABLE]

so

[TABLE]

This leads to the following expression for the envelope dispersion operator:

[TABLE]

where we integrated by parts and introduced an “effective dispersion function”

[TABLE]

The first term on the right-hand side of Eq. (41) is $O(1)$, because at $\epsilon\to 0$, it becomes the GO dispersion function, which is assumed order-one foot:orderone . (At nonzero $\epsilon$, the Weyl symbol $D$ generally does not coincide with the dispersion function that governs the wave in a homogeneous medium; see Sec. IV.1.3.) The second term can be estimated by assuming $\partial/\partial x\sim 1/L$ and $\partial/\partial p\sim 1/k$, since we are interested in $D^{\prime}({\boldsymbol{\rm x}},{\boldsymbol{\rm p}})$ at ${\boldsymbol{\rm p}}$ close to ${\boldsymbol{\rm k}}$. Then,

[TABLE]

meaning that the second term in Eq. (41) is roughly $\epsilon^{2}$ times smaller than the first one. This fact will be used below.

The inverse Wigner–Weyl transform will turn the coordinate ${\boldsymbol{\rm p}}$ into the operator $\smash{\widehat{{\boldsymbol{\rm p}}}}$. The latter will act on the wave envelope [Eq. (33)], which is considered slow in the coordinate representation, so $\smash{\widehat{{\boldsymbol{\rm p}}}}\ket{\psi}=O(\epsilon)$. In this sense, $p=O(\epsilon)$, so by Taylor-expanding $D^{\prime}$ in ${\boldsymbol{\rm p}}$, we are effectively expanding $\smash{\widehat{\mathcal{D}}}$ in $\epsilon$. It is sufficient for our purposes to adopt the second-order expansion,

[TABLE]

where we introduced $D^{\prime}({\boldsymbol{\rm x}})\doteq D^{\prime}\textbf{(}{\boldsymbol{\rm x}},{\boldsymbol{\rm k}}({\boldsymbol{\rm x}})\textbf{)}$ and

[TABLE]

By properties of the inverse Wigner–Weyl transform $\mathscr{W}^{-1}$ (Sec. IV.1.3), $\smash{\widehat{\mathcal{D}}}$ can be written as follows:

[TABLE]

The inverse Wigner–Weyl transforms here can be calculated using Eqs. (28)-(IV.1.3) and $\Theta^{\mu\nu}=\Theta^{\nu\mu}$. This leads to

[TABLE]

One can further simplify Eq. (47) as follows. Suppose that the medium, including the metric, is characterized by some parameters $\mathsf{P}$. Within the accuracy of the theory developed in this paper, we shall neglect terms involving $\partial^{2}\mathsf{P}$ and $(\partial\mathsf{P})(\partial\psi)$ while generally retaining terms such as $\partial^{2}\psi$ (as to be explained Sec. IV.5.2), where $\partial$ denotes a generic spatial derivative. Hence, we can ignore the difference between $D$ and $D^{\prime}$, and we can also omit $\smash{\partial_{\mu\nu}^{2}\Theta^{\mu\nu}}$ in Eq. (47), so $\smash{\widehat{\mathcal{D}}}$ can be expressed as follows:

[TABLE]

Also, $\mathcal{V}^{\mu}$ and $\Theta^{\mu\nu}$ can be calculated through the zeroth-$\epsilon$ limit of $D$, denoted $D_{0}$, which is just the dispersion function of the homogeneous medium:

[TABLE]

Importantly, $D_{0}({\boldsymbol{\rm x}},{\boldsymbol{\rm k}})$ can depend only on ${\boldsymbol{\rm k}}$ and on local parameters of the medium (but not on their derivatives). As a gradient of a scalar [Eq. (4)], ${\boldsymbol{\rm k}}$ is a true covector; hence, $D_{0}({\boldsymbol{\rm x}},{\boldsymbol{\rm k}})$ is a true scalar. This makes $\mathcal{V}^{\mu}$ a vector and $\Theta^{\mu\nu}$ a tensor within the assumed accuracy.

As a side remark, note that vectors like $\mathcal{V}^{\mu}$, which are denoted with upper indices, belong to the space tangent to $M^{n}$. (Such vectors must not be confused with multi-component fields denoted with Latin indices in Sec. V.) Similarly, covectors like $k_{\mu}$, which are denoted with lower indices, belong to the space cotangent to $M^{n}$. For any given $X$, one does not have to distinguish $X^{\mu}$ and $X_{\mu}$ only if the metric is Euclidean. In the case of a general metric $g_{\mu\nu}$, one has

[TABLE]

where the matrix $g^{\mu\nu}$ is the inverse of $g_{\mu\nu}$. Similar rules apply to manipulating tensors, as usual.

IV.3 Approximate envelope equation

Let us derive the ${\boldsymbol{\rm x}}$ representation of Eq. (48), which we need to obtain an explicit form of Eq. (5). From the ${\boldsymbol{\rm x}}$ representation of $\smash{\widehat{{\boldsymbol{\rm x}}}}$ and $\smash{\widehat{{\boldsymbol{\rm p}}}}$ given by Eqs. (13) and (14), terms like $\smash{\widehat{{\boldsymbol{\rm p}}}}f(\smash{\widehat{{\boldsymbol{\rm x}}}})$ (for any $f$) must be interpreted as

[TABLE]

where $g$ and $\psi$ are also functions of ${\boldsymbol{\rm x}}$. Then,

[TABLE]

where all terms are functions of ${\boldsymbol{\rm x}}$, and ;μ denotes the covariant derivative with respect to $x^{\mu}$, so $X^{\mu}{}_{;\mu}$ (for any given vector $X^{\mu}$) is the divergence; namely,

[TABLE]

Except for $D\psi$, all terms in Eq. (51) are covariant (i.e., invariant with respect to coordinate transformations) within the accuracy of our theory. Hence, the symbol $D$ is also covariant within the accuracy of our theory [this includes not just the zeroth-order dispersion function $D_{0}$ but also the leading-order correction $D-D_{0}=O(\partial\mathsf{P})$ that may appear in inhomogeneous medium], even though this may not be obvious directly from the definition of the Weyl transform (Sec. IV.1). To emphasize this fact, and also to ensure consistency of notation with the vector-wave theory to be discussed in the later sections, we denote this covariant function as $\mathsf{D}$, and we also adopt the notation

[TABLE]

Then, the approximate form of the envelope equation (5) can be summarized as follows:

[TABLE]

Here, $\smash{\widehat{\mathcal{D}}}_{H}$ and $\smash{\widehat{\mathcal{D}}}_{A}$ are, respectively, the Hermitian and anti-Hermitian parts of the approximated $\smash{\widehat{\mathcal{D}}}$, namely,

[TABLE]

[TABLE]

[from now on, we do not emphasize the approximate nature of Eq. (51), so $\approx$ will be replaced with $=$], and

[TABLE]

where we assume a standard notation

[TABLE]

Note that ,μ is the “full” derivative in the sense that, for any given $f\textbf{(}{\boldsymbol{\rm x}},{\boldsymbol{\rm k}}({\boldsymbol{\rm x}})\textbf{)}$, it applies both to the first and the second argument of $f$; namely, the chain rule leads to

[TABLE]

The true partial derivatives are introduced as follows:

[TABLE]

The derivatives ,μ and |μ are equivalent for functions that depend only on ${\boldsymbol{\rm x}}$. In particular, ${k_{\nu}\equiv\partial_{\nu}\theta({\boldsymbol{\rm x}})}$ satisfies

[TABLE]

IV.4 Leading-order approximation:

geometrical optics

For the envelope approximation to hold, i.e., for $\psi({\boldsymbol{\rm x}})$ to remain a slow function, ${\boldsymbol{\rm k}}({\boldsymbol{\rm x}})$ must be chosen such that $\mathsf{D}\textbf{(}{\boldsymbol{\rm x}},{\boldsymbol{\rm k}}({\boldsymbol{\rm x}})\textbf{)}\psi\lesssim O(\epsilon)$. [The other terms in Eq. (5) are automatically small on the account of Eq. (2).] Let us adopt the usual GO ordering foot:dissip

[TABLE]

or more rigorously, $\mathsf{D}_{A}\lesssim O(\epsilon)$; i.e., $\mathsf{D}_{A}$ much smaller than $O(\epsilon)$ is also allowed foot:herm . Then, one can define ${\boldsymbol{\rm k}}({\boldsymbol{\rm x}})$ such that

[TABLE]

This can be achieved by calculating ${\boldsymbol{\rm k}}$ as discussed in Sec. III with the ray Hamiltonian ${H=\mathsf{D}_{H}}$ and the initial condition [at any chosen $\smash{{\boldsymbol{\rm x}}^{(0)}}$ and $\smash{{\boldsymbol{\rm k}}^{(0)}}$] such that

[TABLE]

To the first order in $\epsilon$, Eq. (5) is

[TABLE]

and as a corollary,

[TABLE]

At zero $\mathsf{D}_{A}$, the dynamics is conservative, and Eq. (67) reflects conservation of the wave action, or quanta my:amc . Accordingly, $\mathsf{D}_{A}$ determines the dissipation rate, $\mathcal{J}^{\mu}$ can be identified (at least, up to a constant factor) as the action flux density, and $V^{\mu}$ is proportional to the group velocity (in space or in spacetime, depending on the problem). Equations (64)-(67) coincide with the equations of the traditional GO theory book:tracy ; my:amc . Below, we extend them by retaining the second-order terms neglected in Eq. (66).

IV.5 Quasioptical model

IV.5.1 Ray-based coordinates

Let us start by introducing ray-based coordinates as follows foot:mymetric . Suppose multiple rays launched with different ${\boldsymbol{\rm x}}^{(0)}$ within the beam. Their trajectories can be found using the ray equations as discussed in Sec. IV.4; this determines the path $\zeta$ along each ray as a function of the initial coordinate. We treat $\zeta$ as the longitudinal coordinate along the wave beam, and each of its isosurfaces is considered as the transverse space $M^{n-1}_{\perp}(\zeta)$. The coordinates on this space, ${\boldsymbol{\rm\varrho}}\equiv\{\varrho^{1},\varrho^{2},\ldots,\varrho^{n-1}\}$, can be introduced arbitrarily (yet they will be specified below). Then, we define $n-1$ independent vector fields ${\boldsymbol{\rm e}}_{\sigma}$ via

[TABLE]

and another vector ${\boldsymbol{\rm e}}_{0}$ such that it is linearly-independent from all ${\boldsymbol{\rm e}}_{\sigma}$. Hence, any $\mathrm{d}{\boldsymbol{\rm x}}$ can be decomposed as follows:

[TABLE]

and summation over repeated indices $\sigma$ (and $\tilde{\sigma}$) is henceforth assumed from 1 to $n-1$.

Let us also construct the dual basis $\{{\boldsymbol{\rm e}}^{\mu}\}$. (We treat vectors and covectors on the same footing; namely, ${\boldsymbol{\rm X}}\cdot{\boldsymbol{\rm Y}}=X_{\mu}Y^{\mu}=X^{\mu}Y_{\mu}$, and ${\boldsymbol{\rm X}}={\boldsymbol{\rm Y}}$ equally means $X^{\mu}=Y^{\mu}$ and $X_{\mu}=Y_{\mu}$.) By definition,

[TABLE]

so we adopt ${\boldsymbol{\rm e}}^{0}\doteq\nabla\zeta$. Then, the general form of ${\boldsymbol{\rm e}}_{0}$ is

[TABLE]

where the first sign is determined by the metric signature, $\alpha\doteq|{\boldsymbol{\rm e}}^{0}\cdot{\boldsymbol{\rm e}}^{0}|^{-1/2}$, ${\boldsymbol{\rm\beta}}\doteq{\boldsymbol{\rm e}}_{\sigma}\beta^{\sigma}$, and $\beta^{\sigma}$ are arbitrary coefficients. This leads to the following metric representation in the coordinates $\{\zeta,{\boldsymbol{\rm\varrho}}\}$:

[TABLE]

Here, $\boldsymbol{h}$ is the matrix with elements $h_{\sigma\tilde{\sigma}}\doteq{\boldsymbol{\rm e}}_{\sigma}\cdot{\boldsymbol{\rm e}}_{\tilde{\sigma}}$, and ⊺ denotes transposition. By using a known theorem for the determinant of a block matrix foot:det , one obtains

[TABLE]

Let us also introduce the transverse projection of $\mathrm{d}{\boldsymbol{\rm\varrho}}$,

[TABLE]

where $\mathbb{1}$ is a unit matrix and ${\boldsymbol{\rm e}}^{0}{\boldsymbol{\rm e}}^{0}$ is a dyad formed out of ${\boldsymbol{\rm e}}^{0}$. Then, one obtains from Eq. (69) that

[TABLE]

so $h_{\sigma\tilde{\sigma}}$ serves as the transverse metric. Below, we assume ${\boldsymbol{\rm\beta}}=0$, so $\varrho_{\perp}^{\sigma}$ and $\varrho^{\sigma}$ do not need to be distinguished and we can define the inner product on $M^{n-1}_{\perp}$ as follows:

[TABLE]

We also choose transverse coordinates specifically such that ${\boldsymbol{\rm\varrho}}=\text{const}$. In other words, the coordinate ${\boldsymbol{\rm\varrho}}$ of a given point in space is the initial location of the ray that arrives at this point from the initial transverse surface $M^{n-1}_{\perp}(0)$. Then, the transverse group velocity is zero,

[TABLE]

where we used Eqs. (69) and (70). This simplifies our equations below. A different definition of ${\boldsymbol{\rm\varrho}}$ is assumed in Papers II and III, leading to straightforward modifications of the equations, which are not discussed below.

IV.5.2 Quasioptical equation

Suppose a wave beam such that its longitudinal scale $L_{\lVert}$, which is defined via $\psi_{,\zeta}\sim\psi/L_{\lVert}$, is much larger than its perpendicular scale $L_{\perp}$, which is defined via $\psi_{,\sigma}\sim\psi/L_{\perp}$ and may or may not be the beam global width. (We assume the notation $\psi_{,\zeta}\doteq\partial\psi/\partial\zeta$ and $\psi_{,\sigma}\doteq\partial\psi/\partial\varrho^{\sigma}$.) We introduce two GO parameters,

[TABLE]

[so the original parameter (2) is $\epsilon=\epsilon_{\perp}$], and we shall neglect terms smaller than $O(\epsilon_{\lVert}\epsilon_{\perp})$ from now on. Also, we require that either $\alpha$ is chosen as independent of ${\boldsymbol{\rm\varrho}}$ or, more generally, the relative variation of $\alpha$ across the wave beam is small enough such that $\alpha_{,\sigma}$ is negligible. We also assume

[TABLE]

where the derivatives can be taken over the longitudinal or transverse coordinates. (Having $g_{\alpha\beta,\mu}\sim 1/L_{\lVert}$ is indeed typical in typical applications tex:mypaper2 .)

Under these assumptions, Eq. (55) can be approximated as follows [assuming also Eq. (64)]:

[TABLE]

Here, $V\doteq V^{0}$, and we used that $V^{\mu}\psi_{,\mu}=V\psi_{,\zeta}$ due to Eq. (79). The operator $\smash{\widehat{\mathcal{G}}}$ is self-adjoint under the inner product (78) and given by

[TABLE]

(We kept only the transverse derivatives in this term because the longitudinal ones are beyond the assumed accuracy.) Similarly, Eq. (56) becomes

[TABLE]

where the higher-order terms are neglected because $\mathsf{D}_{A}$ is assumed small [Eq. (63)]. Although the last term in Eq. (84) is negligible within the assumed accuracy, it is retained anyway to make $\smash{\widehat{\mathcal{D}}}_{A}$ (not to be confused with $i\smash{\widehat{\mathcal{D}}}_{A}$) exactly self-adjoint under the inner product (78). Then, Eq. (5) becomes

[TABLE]

and, as a corollary,

[TABLE]

where we used the fact that $\smash{\widehat{\mathcal{G}}}$ and $\smash{\widehat{\mathcal{D}}}_{A}$ are self-adjoint. In particular, for zero $\mathsf{D}_{A}$, Eq. (86) predicts conservation of the wave-action flux through the beam cross section,

[TABLE]

Equation (86) indicates that $\smash{\widehat{\mathcal{D}}}_{A}$ cannot be much larger than $\partial/\partial\zeta$, so the correct scaling to assume for the dissipation term is $\mathsf{D}_{A}\lesssim O(\epsilon_{\lVert})$. Then, Eq. (85) leads to the “quasioptical scaling”

[TABLE]

Hence, the difference between $\smash{\widehat{\mathcal{D}}}_{A}$ and $\mathsf{D}_{A}$ in Eq. (84) is $\lesssim O(\epsilon^{3})$, so it is beyond the accuracy of our theory, which is $O(\epsilon^{2})$. For this reason, we henceforth adopt foot:deff

[TABLE]

which will also simplify our notation. [Also, as a reminder, $\partial_{\mu}g_{\alpha\beta}\lesssim O(\epsilon_{\lVert})\sim\epsilon_{\perp}^{2}$, and the parameters of the medium vary similarly.] This leads o

[TABLE]

Equation (90) is a general quasioptical equation for a scalar-wave beam in an inhomogeneous medium. Since it is a parabolic equation (it contains only the first-order derivative with respect to $\zeta$), Eq. (90) is much easier to solve than Eq. (5) with the original expressions for $\smash{\widehat{\mathcal{D}}}_{H}$ and $\smash{\widehat{\mathcal{D}}}_{A}$ given by Eqs. (55) and (56).

IV.5.3 Simplified equations

A typical metric of interest for quasioptical modeling has the form $\smash{g_{\mu\nu}=\delta_{\mu\nu}+R_{\mu\nu\sigma}\varrho^{\sigma}}$ with $\smash{R_{\mu\nu\sigma}=O(L_{\lVert}^{-1})}$ tex:mypaper2 . This satisfies the assumed scaling (81) and corresponds to $\ln g_{\diamond}=O(L_{\perp}/L_{\lVert})$. Then,

[TABLE]

which is negligible. [In contrast, $(\ln g_{\diamond})_{,\sigma}=\smash{O(\epsilon_{\lVert})}$.] The derivatives of $g_{\diamond}$ in Eq. (83) are negligible too. Hence, the metric factors in the quasioptical equation (90) can be replaced with unity. (However, the effect of the metric on the dispersion matrix may remain important in this case; see Sec. VI.) Then, Eq. (90) becomes

[TABLE]

Using the variable transformation $\phi\doteq\sqrt{V}\psi$, this equation can be further simplified as follows:

[TABLE]

which is just a dissipative Schrödinger equation with

[TABLE]

[In Eq. (93), we used that $V_{,\sigma}\phi_{,\tilde{\sigma}}$ is negligible compared to $V\phi_{,\tilde{\sigma}\sigma}=O(\epsilon_{\perp})$.] Also, Eq. (86) becomes

[TABLE]

Note that $\Sigma$ can also be expressed alternatively as follows. Suppose one finds the group velocity and constructs a ray-based metric (as described in Sec. IV.5.1) in the vicinity of a given ${\boldsymbol{\rm x}}$. Then, Eq. (64) can be considered as an equation for the local $k_{0}$ as a function of the transverse wave-vector components $k_{\sigma}$ (and ${\boldsymbol{\rm x}}$) in this prescribed metric; i.e.,

[TABLE]

By differentiating this with respect to $k_{\sigma}$, we obtain $\smash{k_{0}{}^{|\sigma}=-\smash{\mathsf{D}_{H}^{|\sigma}/\mathsf{D}_{H}^{|\zeta}}}$. By differentiating the latter formula once again, now with respect to $k_{\tilde{\sigma}}$, we also obtain

[TABLE]

Recall that $\smash{\mathsf{D}_{H}^{|\sigma}=V^{\sigma}}$ and $V^{\sigma}=0$ [Eq. (79)]. Then, $\smash{k_{0}{}^{|\sigma}=0}$, and by comparing Eq. (97) with Eq. (94) and using Eq. (57), one finds that $\smash{k_{0}{}^{|\sigma\tilde{\sigma}}=-\Sigma^{\sigma\tilde{\sigma}}}$. Thus, Eq. (93) can be rewritten as

[TABLE]

In the case of a spacetime problem, $\smash{k_{0}{}^{|\sigma\tilde{\sigma}}}$ can also be expressed through the wave frequency in the laboratory frame. Then, the familiar Schrödinger equation for the wave envelope foot:karpman can be reproduced. We do not discuss it further for this subject is not essential to our paper.

V Vector waves

V.1 Hilbert space for vector waves

Now, let us generalize the above results to vector waves. Suppose an $m$-component wave field $\Psi\equiv\Psi({\boldsymbol{\rm x}})$, so

[TABLE]

Here, $\ket{\Psi}$ is a vector on some Hilbert space $\mathcal{H}_{m}^{n}$ for which we assume the following general inner product:

[TABLE]

Here, the Latin indices that characterize the field components span from 1 to $m$. (This is in contrast with the Greek indices introduced earlier, which characterize the components of ${\boldsymbol{\rm x}}$ and span from 0 to $n-1$.) The matrix $\gamma({\boldsymbol{\rm x}})$ can be any $m\times m$ symmetric matrix. It can be considered as an additional metric. This metric does not have to be the same as $g$, as seen already from the fact that $m$ and $n$ do not have to be the same. (For example, LABEL:my:qdiel describes a six-dimensional field on a three-dimensional spacetime; also see Sec. VI.2.) This means that the space to which $\Psi({\boldsymbol{\rm x}})$ belongs is not necessarily the space tangent to $M^{n}$. That said, having $\gamma=g$ is possible as a special case; see Secs. VI.1 and VI.3.

Using $\gamma$ as a metric, we introduce the standard rules for manipulating the Latin indices,

[TABLE]

where $\gamma^{ab}$ are elements of $\gamma^{-1}$. In particular, in the ${\boldsymbol{\rm x}}$ representation, one has

[TABLE]

and we shall also use the following local dot product:

[TABLE]

Under this dot product, a matrix $A$ with mixed-index elements $A^{a}{}_{b}$ is self-adjoint [$(A\Psi)\cdot\Phi=\Psi\cdot(A\Phi)$] if the matrix with the corresponding lower-index elements, $A_{ab}$, is Hermitian. The difference between Hermitian and self-adjoint matrices can be ignored if the metric $\gamma$ is Euclidean or pseudo-Euclidean.

V.2 Operators on $\boldsymbol{\mathcal{H}_{m}^{n}}$

Any operator $\smash{\widehat{A}}$ on $\mathcal{H}_{m}^{n}$ can be understood as an $m\times m$ matrix with elements $\smash{\widehat{A}}^{a}{}_{b}$ which are operators on $\mathcal{H}_{1}^{n}$. For any given $\smash{\widehat{A}}^{a}{}_{b}$, we introduce $\underline{\smash{\widehat{A}}}_{ab}\doteq\gamma_{ac}\smash{\widehat{A}}^{c}{}_{b}$, which is also an operator on $\mathcal{H}_{1}^{n}$. Correspondingly, two types of Weyl images can be defined, namely, $\mathscr{W}_{z}[\smash{\widehat{A}}^{a}{}_{b}]$ and $\mathscr{W}_{z}[\underline{\smash{\widehat{A}}}_{ab}]$. Below, we assume the notation

[TABLE]

Also, $A$ will be an index-free notation of the matrix with elements $A^{a}{}_{b}$, and $\underline{A}$ will be an index-free notation of the matrix with elements $\underline{A}_{ab}$. Importantly, $A^{a}{}_{b}$ should not be confused with $\mathscr{W}_{z}[\smash{\widehat{A}}^{a}{}_{b}]$, because the Weyl symbols do not transform as tensors,

[TABLE]

except for those which are independent of ${\boldsymbol{\rm p}}$. Symbols $\mathscr{W}_{z}[\smash{\widehat{A}}^{a}{}_{b}]$ will not be used below, so no special notation is introduced for them.

Note that by Eq. (108), the following equality holds for any $\Phi$, $\Psi$, and $\smash{\widehat{A}}$ on $\mathcal{H}_{m}^{n}$:

[TABLE]

where we invoked the definition of the inner product on $\mathcal{H}_{1}^{n}$ [Eq. (10)]. By definition of the adjoint operator $\smash{\widehat{A}}^{\dagger}$,

[TABLE]

In order to express $\smash{\widehat{A}}^{\dagger}$ through $\smash{\widehat{A}}$, let us represent $\underline{\smash{\widehat{A}}}_{ab}$ in terms of the corresponding Weyl image $\underline{A}_{ab}$ and the Wigner operator $\smash{\widehat{\Delta}}_{z}$ on $\mathcal{H}_{1}^{n}$, with $z\equiv({\boldsymbol{\rm x}},{\boldsymbol{\rm p}})$ (Sec. IV.1.3). Then, Eq. (115) leads to

[TABLE]

where we used that $\smash{\widehat{\Delta}}_{z^{\prime}}$ is self-adjoint on $\mathcal{H}_{1}^{n}$. By comparing this with Eq. (115), one finds that, on one hand,

[TABLE]

On the other hand, by Eq. (27) for the inverse Wigner–Weyl transform, one has

[TABLE]

Hence, the Weyl image of $\smash{\widehat{A}}^{\dagger}$ is simply the conjugate transpose of the Weyl image of $\smash{\widehat{A}}$ in the sense that

[TABLE]

In other words, for operators with both indices lowered, the operations $\mathscr{W}_{z}$ and † commute.

Equation (119) implies that, if $\smash{\widehat{A}}$ is self-adjoint on $\mathcal{H}_{m}^{n}$, then $\underline{A}$ is a Hermitian matrix, and thus the corresponding matrix $A$ is self-adjoint (Sec. V.1). In particular, this means that the Hermitian and anti-Hermitian parts [Eq. (26)] of an operator are determined by, respectively, the Hermitian and anti-Hermitian parts of its lower-index Weyl symbol. We summarize this as follows:

[TABLE]

For lower-index matrices, the Hermitian and anti-Hermitian parts are defined as usual:

[TABLE]

V.3 Weyl expansion of the dispersion operator

Now, let us consider the dispersion operator $\smash{\widehat{D}}$ in particular. According to the above definition, it is viewed as an $m\times m$ matrix with elements $\smash{\widehat{D}}^{a}{}_{b}$. By lowering the index in the envelope equation (1), one can write this equation as follows:

[TABLE]

where $\underline{\smash{\widehat{D}}}_{ab}\doteq\gamma_{ac}\smash{\widehat{D}}^{c}{}_{b}$. Equivalently, this can be written as the envelope equation

[TABLE]

Then, the approximate operators $\underline{\smash{\widehat{\mathcal{D}}}}_{ab}$ can be expressed just like $\smash{\widehat{\mathcal{D}}}$ for scalar waves (Sec. IV.2),

[TABLE]

Here, the coefficients are ${\boldsymbol{\rm x}}$-dependent matrices; specifically, $\underline{D}_{ab}({\boldsymbol{\rm x}})\doteq\underline{D}_{ab}\textbf{(}{\boldsymbol{\rm x}},{\boldsymbol{\rm k}}({\boldsymbol{\rm x}})\textbf{)}$, and

[TABLE]

One can also multiply Eq. (125) by $\gamma^{-1}$ to raise the first index. Then, one obtains

[TABLE]

Note that if $\smash{\widehat{D}}$ is self-adjoint, then $\underline{D}_{ab}$ is Hermitian; then, $\smash{\widehat{\mathcal{D}}}^{a}{}_{b}$ is self-adjoint too.

Unlike in scalar waves (Sec. IV.3), $\underline{D}_{ab}({\boldsymbol{\rm x}},{\boldsymbol{\rm k}})$ is not a covariant object, i.e., in this case, not a true tensor. (For example, see Sec. VI.) It is convenient to split it as

[TABLE]

where $\underline{\mathsf{D}}$ is a true tensor, and $\underline{C}=O(\partial\mathsf{P})$ is a small remainder. [Here $\partial\mathsf{P}$ means the same as in Sec. IV.3, so $C=O(\epsilon)$ under the GO ordering, and $C=O(\epsilon_{\parallel})$ under the quasioptical ordering.] The splitting Eq. (130) is not unique and is a matter of convenience; for example, one can define $\underline{\mathsf{D}}({\boldsymbol{\rm x}},{\boldsymbol{\rm k}})$ simply as the zeroth-order dispersion tensor. Hence, we can redefine

[TABLE]

which is equivalent to the above definitions within the accuracy of our theory. Then, assuming the same orderings as in the scalar-wave case (also see below), the following approximation is enough both for the first-order theory and for the quasioptical theory:

[TABLE]

Here, $\mathsf{D}\equiv\mathsf{D}({\boldsymbol{\rm x}})$, and $\mathsf{D}=\gamma^{-1}\underline{\mathsf{D}}$, so $\mathsf{D}^{a}{}_{b}({\boldsymbol{\rm x}})=\gamma^{ac}\underline{\mathsf{D}}_{cb}({\boldsymbol{\rm x}})$; similar conventions are assumed for $C$ and $D$. Accordingly, the matrices $\mathsf{D}_{H}$ and $D_{A}$ are self-adjoint, and so are the operators $\smash{\widehat{\mathscr{D}}}_{H1}$ and $\smash{\widehat{\mathscr{D}}}_{H2}$, which are given by

[TABLE]

Using $[\gamma^{-1},\smash{\widehat{p}}_{\mu}]=i(\gamma^{-1})_{,\mu}$ and $(\gamma^{-1}\gamma)_{,\mu}=0$ to obtain

[TABLE]

one can also express these as

[TABLE]

(The term $\propto\ell_{\mu}$ must be kept in $\smash{\widehat{\mathscr{D}}}_{H1}$ but a similar term in $\smash{\widehat{\mathscr{D}}}_{H2}$ can be neglected.) Also,

[TABLE]

and since $\gamma$ is ${\boldsymbol{\rm k}}$-independent, this leads to

[TABLE]

If $\smash{\widehat{\mathscr{D}}}_{H2}$ is neglected, the envelope equation (129) becomes a Dirac-type equation similar to those considered in the context of XGO phd:ruiz17 ; my:covar ; my:qdiel ; my:qdirac ; my:qdirpond and also, for instance, in Refs. ref:bridges01 ; ref:friedland87 . We shall also discuss this limit in Sec. V.5.

V.4 Active and passive modes

V.4.1 Basis from the eigenvectors of $\mathsf{D}_{H}$

Let us assume the GO ordering as in Eq. (63), except that $\mathsf{D}_{H}$ and $\mathsf{D}_{A}$ are now matrices. Then, $\smash{\widehat{\mathscr{D}}}=O(\epsilon)$, so the envelope equation [Eq. (5)] acquires the form

[TABLE]

Hence, it is convenient to express $\psi$ through the eigenvectors of $\mathsf{D}_{H}({\boldsymbol{\rm x}})$. Let us denote these eigenvectors as $\eta_{s}({\boldsymbol{\rm x}})$ and the corresponding eigenvalues as $\Lambda_{s}({\boldsymbol{\rm x}})$, so

[TABLE]

[Note that $\eta_{s}({\boldsymbol{\rm x}})=\eta_{s}\textbf{(}{\boldsymbol{\rm x}},{\boldsymbol{\rm k}}({\boldsymbol{\rm x}})\textbf{)}$, where $\eta_{s}({\boldsymbol{\rm x}},{\boldsymbol{\rm p}})$ are the eigenvectors of $\mathsf{D}_{H}({\boldsymbol{\rm x}},{\boldsymbol{\rm p}})$. Likewise, $\Lambda_{s}({\boldsymbol{\rm x}})=\Lambda_{s}\textbf{(}{\boldsymbol{\rm x}},{\boldsymbol{\rm k}}({\boldsymbol{\rm x}})\textbf{)}$, where $\Lambda_{s}({\boldsymbol{\rm x}},{\boldsymbol{\rm p}})$ are the eigenvalues of $\mathsf{D}_{H}({\boldsymbol{\rm x}},{\boldsymbol{\rm p}})$.] Since $\mathsf{D}_{H}$ is self-adjoint, it has $m$ eigenvectors $\eta_{s}$ that form a complete basis $\{\eta_{s}\}$. Let us also introduce the corresponding dual basis $\{\eta^{s}\}$ as usual, i.e., such that

[TABLE]

Then, we can represent $\psi$ as

[TABLE]

where $a^{s}$ serve as the components of $\psi$ in the basis $\{\eta_{s}\}$. [Remember that the dot product (111) includes conjugation.] Since $\mathsf{D}_{H}$ is self-adjoint, it is always possible to make the basis $\{\eta_{s}\}$ orthonormal, and this choice is assumed below. Then, $\eta_{s}\cdot\eta_{s^{\prime}}=\delta_{ss^{\prime}}$; thus $\eta_{s}=\eta^{s}$, i.e.,

[TABLE]

Also, for any two fields $\psi=\eta_{s}a^{s}$ and $\phi=\eta_{s}b^{s}$, the inner product (108) can be written as follows:

[TABLE]

where $a_{s}\doteq\delta_{ss^{\prime}}a^{s^{\prime}}$. The matrix $\delta$ with elements $\delta_{ss^{\prime}}$ serves as a Euclidean metric for manipulating the mode indices $s$ and $s^{\prime}$. Those should not to be confused with the coordinate indices denoted with Greek letters and also with other Latin indices that are manipulated by the metric $\gamma$ [Eq. (109)].

V.4.2 Amplitude vectors

The physical meaning of the expansion coefficients $a^{s}$ is understood as follows. Consider multiplying Eq. (143) by $\eta^{s}$ from the left. That gives $\Lambda_{s}a^{s}=O(\epsilon)$, where no summation over $s$ is assumed. This shows that, for a given $s$, there are two possibilities: either $a^{s}$ is small or $\Lambda_{s}$ is small. In the first case, the polarization $\eta_{s}$ does not correspond to a propagating wave mode per se; the small nonzero projection of $\psi$ on $\eta_{s}$ is only due to the fact that the wave field is not strictly sinusoidal in inhomogeneous medium. We call such “modes” passive. In the second case, the local ${\boldsymbol{\rm k}}({\boldsymbol{\rm x}})$ approximately satisfies the local dispersion relation $\Lambda_{s}\textbf{(}{\boldsymbol{\rm x}},{\boldsymbol{\rm k}}({\boldsymbol{\rm x}})\textbf{)}\approx 0$. Then, $a^{s}$ can be understood as the local scalar amplitude of an actual GO mode, so that $a^{s}=O(1)$ is allowed. We call such modes active, and mode conversion occurs when more than one active mode exists. In other words, for a given boundary or initial conditions, active modes are those that are excited resonantly, while passive modes are those that are nonresonant and, thus, adiabatically isolated. (However, this does not mean that the passive-modes amplitudes are simply negligible; see below.)

Assuming that there are $N\geq 1$ active modes, we shall order them such that they correspond to $s=1,2,\ldots,N$ and the remaining $\bar{N}=m-N$ modes with $s={(N+1)},{(N+2)},\ldots,m$ are passive. Let us also adopt the notation

[TABLE]

($s=1,2,\ldots,\bar{N}$) for the passive-mode amplitudes, polarizations, and eigenvalues. Then, it is convenient to introduce the following “amplitude vectors”

[TABLE]

and the corresponding row vectors that are dual to the amplitude vectors under the Euclidean complex dot product; namely,

[TABLE]

Below, we seek to derive an approximate envelope equation in terms of these amplitude vectors.

V.4.3 Polarization matrices

Before we proceed, let us introduce the following notation. First, consider the “polarization matrices”

[TABLE]

These are non-square matrices that have active- and passive-mode polarizations as their columns; namely,

[TABLE]

Using these, Eq. (146) can be rewritten compactly as

[TABLE]

Also consider the auxiliary polarization matrices

[TABLE]

which have the $\eta^{s*}$ as their rows,

[TABLE]

Since

[TABLE]

one obtains

[TABLE]

However, note that in general, + does not mean to represent the adjoint in the common sense [cf. Eq. (176)].

Next, notice that

[TABLE]

where the indices $s$ and $s^{\prime}$ span from 1 to $N$. Analogous formulas apply to $\bar{\Xi}$. Then,

[TABLE]

and similarly for $\Xi^{+}\bar{\Xi}$. In other words, one has

[TABLE]

where $\mathbb{1}$ is the unit square matrix and $\mathbb{0}$ is the zero square matrix, correspondingly. [The dimensions of $\mathbb{1}$ and $\mathbb{0}$ are different in different equations.] Hence, if one multiplies Eq. (161) by $\Xi^{+}$ and, separately, by $\bar{\Xi}^{+}$, one obtains concise formulas for the amplitude vectors $a$ and $\bar{a}$,

[TABLE]

Another property that we shall use later on is

[TABLE]

and similarly for $\bar{\Xi}^{+}$. Here, H denotes the conjugate transpose ($\Xi^{\sf H}\doteq\Xi^{\intercal*}$), and $\delta^{-1}$ is the Kronecker matrix with upper-index elements $\smash{\delta^{ss^{\prime}}}$. Hence,

[TABLE]

The factor $\delta^{-1}$ can be omitted if one ignores the difference between upper and lower indices $s$ and $s^{\prime}$ that refer to the mode number. If $\gamma$ is Euclidean, one can also ignore the difference between upper and lower coordinate indices; then, $\Xi^{+}=\Xi^{\sf H}$.

As a side remark, note that $\Xi$ and $\bar{\Xi}$ are generally non-square, so $\Pi\doteq\Xi\Xi^{+}$ and $\bar{\Pi}\doteq\bar{\Xi}\bar{\Xi}^{+}$ are not unit matrices but rather projectors. Indeed, due to Eqs. (174), one has $\Pi^{2}=\Pi$ and $\bar{\Pi}^{2}=\bar{\Pi}$. Also, by applying $\Pi$ and $\bar{\Pi}$ to Eq. (161), one obtains

[TABLE]

Thus, $\Pi\psi$ is the projection of $\psi$ on the active-mode space, and $\bar{\Pi}\psi$ is the projection of $\psi$ on the passive-mode space.

V.4.4 Equation for the active modes

Using the eigendecomposition theorem and Eq. (174), one obtains

[TABLE]

where $\Lambda$ and $\bar{\Lambda}$ are the diagonal eigenvalue matrices,

[TABLE]

These and Eq. (174) also lead to

[TABLE]

Hence, the envelope equation (5) can be written as

[TABLE]

where we used Eq. (181). Let us multiply this by $\Xi^{+}$ and, separately, by $\bar{\Xi}^{+}$. Then, due to Eq. (174), one obtains

[TABLE]

Let us also use Eq. (184) to express $\bar{a}$ through $a$ and substitute the result into Eq. (183). Since $\bar{a}$ and $\smash{\widehat{\mathscr{D}}}$ are both of order $\epsilon$, it is sufficient to solve Eq. (184) for $\bar{a}$ approximately to the leading order,

[TABLE]

Then, Eq. (183) becomes an equation for just the $N$-dimensional active-mode amplitude vector,

[TABLE]

As a reminder, this equation is valid up to $O(\epsilon^{2})$.

V.5 Leading-order approximation:

extended geometrical optics

Upon substituting Eqs. (134) into Eq. (186), we obtain the following equation to lowest (first) order in $\epsilon$:

[TABLE]

where we introduced

[TABLE]

As shown in Appendix B.1,

[TABLE]

where $\Xi$ is considered as a function of ${\boldsymbol{\rm x}}$. [As a reminder, $\Xi^{\sf H}\doteq\Xi^{\intercal*}$ is the conjugate transpose of $\Xi$; $\underline{\mathcal{V}}_{H}^{\mu}$ is given by Eq. (131); and $\delta^{-1}$ is a unit matrix that only raises the mode index.] Alternatively, $\Xi$ can be considered as a function of $({\boldsymbol{\rm x}},{\boldsymbol{\rm k}})$. Then, as shown in Appendix B.2,

[TABLE]

where the partial derivatives |μ and |μ are defined in Sec. IV.3. The term $U$ is the Stern–Gerlach Hamiltonian mentioned in the introduction (Sec. I), except here it is generalized to an arbitrary metric. This term causes polarization-driven bending of the ray trajectories, which is missed in traditional GO; also, it causes mode conversion, if more than one active mode is present phd:ruiz17 ; my:covar ; my:qdiel ; my:qdirac ; my:qdirpond . These effects are is discussed in further detail in Sec. V.7.

More explicitly, Eq. (187) can be written as

[TABLE]

This generalizes the XGO equation derived in Refs. phd:ruiz17 ; my:covar ; my:qdiel ; my:qdirpond ; my:qdirac ; my:xo to dissipative waves and curved coordinates. Since the vector $a$ belongs to the space where the metric is Euclidean by definition (Sec. V.4.2), the elements of this vector are covariant, i.e., invariant with respect to coordinate transformations in the physical space. Same applies to the eigenvalues comprising $\Lambda$; thus, the whole left-hand-side of Eq. (194) is covariant too. This means that the operator on the right-hand side is also covariant, which includes its Hermitian and anti-Hermitian parts separately. Hence, $U$ and $\Gamma$ are covariant.

Also note that Eq. (194) is similar to Eq. (66) for scalar waves and has a similar corollary,

[TABLE]

Using Eqs. (110), (140), (161), (171), (177), and (191) one obtains

[TABLE]

so Eq. (195) can be viewed as a generalization of Eq. (67). At zero $\Gamma$, the dynamics is conservative, and Eq. (195) reflects conservation of the wave action, or quanta my:amc . Accordingly, $\Gamma$ determines the dissipation rate, $\mathcal{J}^{\mu}$ can be identified (up to a constant factor) as the action flux density summed over all active modes (for example, see Sec. VI.2), and the elements of the diagonal matrix $\Lambda^{|\mu}$ are proportional to the group velocities of the corresponding active modes. We shall call them the group velocities (without “proportional to”) for brevity. We also emphasize that this model applies even when the group velocities of the active modes are very different, unlike the quasioptical model discussed in Sec. V.6.

For scalar waves studied in Sec. IV, we had $\Lambda=\mathsf{D}_{H}$ and we defined ${\boldsymbol{\rm k}}({\boldsymbol{\rm x}})$ such that this term be zero [Eq. (64)]. Now, $\Lambda$ is a diagonal matrix with $N\geq 1$ nonzero elements, so it cannot be zeroed entirely by imposing just one scalar constraint on ${\boldsymbol{\rm k}}({\boldsymbol{\rm x}})$. Hence, there can be more than one natural way to define ${\boldsymbol{\rm k}}({\boldsymbol{\rm x}})$. One aesthetically pleasing and convenient my:xo option is to require that $\Lambda$ or $\Lambda-U$ be traceless. This amounts to choosing the reference-ray Hamiltonian in Sec. III as ${H=\text{tr}\,\Lambda/N}$ or $H={\text{tr}\,(\Lambda-U)/N}$, correspondingly. [Arbitrary constant factors can be introduced instead of $N$, for that only redefines $\tau$ in Eq. (7).] Another option is to adopt $\Lambda_{s}=0$ for some single ${s\leq N}$, since all active modes have close wave vectors anyway; then $H=\Lambda_{s}$. As mentioned in Sec. III, all such choices of ${\boldsymbol{\rm k}}({\boldsymbol{\rm x}})$ are equally justified. Although they lead to slightly different envelope equations, those equations are equivalent in the sense that they all describe the same total field by construction. (One might call this a gauge freedom.) However, for the case $N=1$, choosing $H=\Lambda-U$ is preferable, as discussed in Sec. V.7.2.

V.6 Quasioptical model

V.6.1 Quasioptical equation

Suppose that a wave propagates largely as a single beam. This implies that all active modes have group velocities close to their average group velocity ${\boldsymbol{\rm V}}_{\rm avr}$, which can be defined via $NV_{\rm avr}^{\mu}\doteq\text{tr}\,\Lambda^{|\mu}$. Then,

[TABLE]

where $\Delta V^{\mu}$ are matrices with eigenvalues much smaller than $V_{\rm avr}$. Let us assume a ray-based coordinate system aligned with ${\boldsymbol{\rm V}}_{\rm avr}$ with the inner product on the transverse space $M^{n-1}_{\perp}$ defined similarly to Eq. (78) [cf. also Eq. (148)], namely,

[TABLE]

Let us also adopt the quasioptical ordering as we did for scalar waves in Sec. IV.5. In particular, we allow $\Delta V^{\mu}=O(\epsilon_{\perp})$. Then, Eq. (186) becomes

[TABLE]

Here, $\smash{\widehat{K}}$ is defined as in Eq. (189) and can be approximated as follows:

[TABLE]

where $V\doteq V_{\rm avr}^{0}$ (we used the fact that $V_{\rm avr}^{\sigma}=0$ by definition of the ray-based coordinates), $U=O(\epsilon_{\lVert})$ is given by Eq. (193), and $\Gamma$ is given by Eq. (188). Finally,

[TABLE]

As shown in Appendix B.3, $\smash{\widehat{\mathcal{G}}}$ can also be simplified as

[TABLE]

so it is also self-adjoint under the inner product (198).

In summary then, the quasioptical equation for vector waves can be written as follows:

[TABLE]

Equation (205) is the main result of our paper. Like Eq. (194), this equation is covariant within the accuracy of our theory. As a reminder, $a$ is generally a vector (${\text{dim}\,a=N\geq 1}$), whose elements are the scalar amplitudes of active modes (Sec. V.4.2); $\Lambda$ is the diagonal matrix formed by the eigenvalues of $\mathsf{D}_{H}$; $\Gamma$ is given by Eq. (188); $V=N^{-1}\text{tr}\,\Lambda^{|\zeta}$ is a scalar; ${}^{|\zeta}\doteq\partial/\partial k_{0}$; $\smash{\widehat{\mathcal{G}}}$ is given by Eq. (204); and $\Delta\smash{\widehat{K}}$ is given by Eq. (201). There, $U$ is given by Eq. (193), and $\Delta V$ is given by Eq. (202). The derivatives of $\smash{h_{\diamond}}$ could be neglected within the assumed accuracy, but a purist might want to retain them in order to keep $\Delta\smash{\widehat{K}}$ and $\smash{\widehat{\mathcal{G}}}$ precisely self-adjoint under the inner product (198).

Like in Sec. IV.5.2, we require the $\alpha$ parameter of the ray-based coordinates to be defined as ${\boldsymbol{\rm\varrho}}$-independent. Then, Eq. (205) has the following corollary:

[TABLE]

which is similar to Eq. (86). For the case of zero $\Gamma$, Eq. (206) predicts conservation of the wave-action flux through the beam cross section,

[TABLE]

Unlike in Eq. (87), this is the flux of all active modes combined, and the fluxes of individual active modes may not be conserved separately.

V.6.2 Simplified equations

If the metric is nearly Euclidean (or pseudo-Euclidean) as in Sec. IV.5.3, we can drop the metric factors and rewrite Eq. (205) as follows:

[TABLE]

Using the variable transformation $\phi\doteq\sqrt{V}a$, this equation can be further simplified as

[TABLE]

Here, $\smash{\widehat{\chi}}$ is an operator self-adjoint under the inner product (198) and given by

[TABLE]

Also, $Q$, $\Sigma$, and $\Upsilon$ are self-adjoint matrices given by

[TABLE]

Accordingly, Eq. (206) becomes

[TABLE]

so $\smash{\braket{\phi}{\phi}^{\perp}_{N}}$ is conserved if $\Upsilon$ is zero.

V.7 Summary and discussion

The quasioptical model proposed above is equally applicable to scalar beams, single-mode vector beams (${N=1}$), and multi-mode vector beams (${N>1}$).

V.7.1 Scalar beams

In the case of a scalar beam (Sec. IV), one has $\Xi=1$, so

[TABLE]

Then, the equations from Sec. IV.5 are reproduced. In particular, $\Lambda$ is made zero by the choice of ${\boldsymbol{\rm k}}$ [Eq. (64)]. This implies that the ray Hamiltonian is $H=\Lambda$, which leads to the following ray equations:

[TABLE]

V.7.2 Single-mode vector beams

In the case of a single-mode vector beam ($N=1$), one has $\Xi=\eta$, where $\eta$ is the polarization vector. Then,

[TABLE]

so they are scalars. Also, $\Delta V^{\sigma}=0$, but $U$ is generally a nonzero scalar function. There are two natural ways to define ${\boldsymbol{\rm k}}$ in this case. One is to require that $\Lambda=0$, i.e., to adopt $H=\Lambda$. In this case, $U$ is left in the envelope equation and can intensify the envelope inhomogeneity in the direction perpendicular to the group velocity. This may eventually result in the violation of the envelope approximation. (For multi-mode beams, this is less of a concern because they typically split into single-mode beams before the issue becomes significant.) Hence, it is potentially advantageous to define ${\boldsymbol{\rm k}}$ by requiring $\Lambda-U=0$, i.e., by adopting $H=\Lambda-U$. In this case, the amplitude equation is identical to that of a scalar wave while the effect of $U$ is absorbed in the reference phase, leading to modified ray equations,

[TABLE]

The effect of $U$ on the ray trajectories is known as the (spin) Hall effect of light in optics ref:onoda04 ; ref:bliokh08 ; ref:bliokh15 ; ref:gosselin07 ; ref:duval17 ; foot:marius . In plasma physics, this effect is typically neglected. (To our knowledge, all existing ray-tracing codes ignore $U$ entirely.)

Finally, note that the ray dynamics can also be represented in an alternative form. Since $\Lambda$ is a scalar, one can rewrite Eq. (193) for $U({\boldsymbol{\rm x}},{\boldsymbol{\rm k}})$ as follows:

[TABLE]

Since there is only one mode, there is no mode index to manipulate, so the “metric” factor $\delta$ has been dropped. Likewise, the anti-Hermitian parts are simply the imaginary parts in this case. We also substituted the GO ray equations (213), where we replaced $\mathrm{d}/\mathrm{d}\tau$ with dots, and adopted

[TABLE]

Then, the phase-space Lagrangian of a ray, $\mathsf{L}[{\boldsymbol{\rm x}},{\boldsymbol{\rm k}}]=k_{\mu}\dot{x}^{\mu}-H$, has a non-canonical structure, namely,

[TABLE]

The equations originating from a special case of this non-canonical phase-space Lagrangian were used, for example, in LABEL:ref:bliokh08 to study the Hall effect for light propagating in non-birefringent material. A comparison of the resulting non-canonical ray equations with the canonical ray equations (215) can be found in LABEL:my:qdiel. Also notably, the terms $\smash{\mathsf{A}_{\mu}^{(x)}}$ and $\smash{\mathsf{A}^{\mu}_{(k)}}$ are known as the Berry connection and are linked to the Berry phase, or the geometrical phase of light ref:bliokh08 ; my:qdiel ; ref:littlejohn91 .

V.7.3 Multi-mode vector beams

In the case of multiple active modes ($N>1$), $a$ is an $N$-dimensional vector, and the coefficients in the amplitude equation are matrices (except for $V$, which is a scalar). In particular, the matrices $\Gamma$, $\Delta V^{\sigma}$, and $U$ are generally nondiagonal, so they can cause mode conversion. In the simplest case, when both the transverse gradients and dissipation are negligible, this process is most transparently described by Eq. (209), which then becomes

[TABLE]

with $Q$ being a self-adjoint matrix. The modes decouple when $Q$ is close to diagonal. More generally, the dynamics governed by Eq. (220) is similar to that of an $N$-level quantum system with a Hamiltonian $Q$ (and to the dynamics of $N$ coupled classical oscillators whose parameters do not change too rapidly with $\zeta$ my:wkin ). Alternatively, it can be mapped to a precession equation for a real $(N^{2}-1)$-dimensional “spin” vector my:qdirac . This approach is particularly intuitive at $N=2$, when the spin is three-dimensional. A detailed analysis of mode conversion for this case can be found in LABEL:my:xo.

As a reminder, the model presented here (Sec. V.6) relies on the assumption that the group velocities of the active modes are close to each other. Otherwise, beam splitting occurs rapidly, and the assumption that $a_{,\zeta}\ll a_{,\sigma}$ does not hold. The quasioptical description is inapplicable to such beams; however, the first-order XGO model described in Sec. V.5 can be used.

VI Electromagnetic waves

VI.1 Covariant formulation

Here, we shall explain how the above theory applies to EM waves. We start with Maxwell’s equation book:landau2

[TABLE]

where $F^{\alpha\beta}$ is the EM tensor, namely,

[TABLE]

$A$ is the vector potential, and $J$ is the current density. This equation can also be represented as follows:

[TABLE]

Here $\smash{\widehat{D}}_{\rm vac}^{(A)}$ is the vacuum dispersion operator,

[TABLE]

which is self-adjoint under the inner product (108). By assuming $J=c\smash{\widehat{\vartheta}}A$, were $\smash{\widehat{\vartheta}}$ is some linear operator, one can cast the equation for $A$ in the form (1), namely,

[TABLE]

Hence, the theory developed in the previous sections readily applies, specifically, with the dimension of the configuration space being $n=4$, the dimension of the vector field $A$ being $m=4$ (so $m=n$), and $\gamma=g$. The corresponding Weyl symbols are calculated as follows. It can be shown straightforwardly that $\smash{[\smash{\widehat{\underline{D}}}_{\rm vac}^{(A)}]{}_{\alpha\beta}\doteq g_{\alpha\gamma}[\smash{\widehat{D}}_{\rm vac}^{(A)}]^{\gamma}{}_{\beta}}$ is approximately given by

[TABLE]

Here, we omitted second-order derivatives of $g$, which are negligible within the accuracy of our theory. Also, $q$ is defined as in Eq. (15), and $\ell_{\mu}\doteq g^{-1}g_{,\mu}$, which is in agreement with Eq. (137) because $g=\gamma$. Finally, in the form introduced in Sec. V.3, the above result is

[TABLE]

where $(\underline{\vartheta}_{0})_{\alpha\beta}$ is the GO limit of $\underline{\vartheta}$, $(\Delta\underline{\vartheta})_{\alpha\beta}$ is the remaining part of $\underline{\vartheta}$, and $\underline{\vartheta}$ itself is the symbol of $\underline{\smash{\widehat{\vartheta}}}$.

As a side remark, these equations can be used to calculate the Hall effect of light in gravitational fields, i.e., the deviation of vacuum light rays from geodesics in curved spacetime. (Special cases of this effect are discussed in Refs. ref:gosselin07 ; ref:duval17 ; foot:marius .)

VI.2 Non-covariant formulation

As a special case, suppose a metric of the form

[TABLE]

with time-independent spatial metric $\boldsymbol{h}$. Then, the above equations can be simplified as follows. Let us assume the Weyl gauge ($A_{0}=0$). Then, it is sufficient to consider just the spatial part of the four-dimensional Eq. (225) and replace the vector potential with the electric field

[TABLE]

(As usual, $\smash{\widehat{p}}_{0}\doteq-ic^{-1}\partial_{t}$, which is a self-adjoint operator. Assuming we work with wave fields that have zero time average, $\smash{\widehat{p}}_{0}$ can also be considered reversible.) Unlike in the previous sections, we now use the Latin indices to denote spatial components, and we shall adopt the bold font for spatial vectors and matrices when using index-free notation. Specifically, one obtains

[TABLE]

where we also multiplied the equation by $(\smash{\widehat{p}}_{0})^{-1}$. This has the form (1) with $n=4$, $m=3$, and $\gamma=\boldsymbol{h}$.

Let us introduce the conductivity operator $\smash{\widehat{{\boldsymbol{\rm\sigma}}}}$ via ${\boldsymbol{\rm J}}=\smash{\widehat{{\boldsymbol{\rm\sigma}}}}{\boldsymbol{\rm E}}$ and notice that $4\pi\smash{\widehat{{\boldsymbol{\rm\sigma}}}}=ic\smash{\widehat{p}}_{0}\smash{\widehat{{\boldsymbol{\rm\chi}}}}$ by definition of the susceptibility operator $\smash{\widehat{{\boldsymbol{\rm\chi}}}}$. Then, $4\pi\smash{\widehat{{\boldsymbol{\rm\vartheta}}}}=\smash{\widehat{p}}_{0}\smash{\widehat{{\boldsymbol{\rm\chi}}}}\smash{\widehat{p}}_{0}$, so $\smash{\widehat{{\boldsymbol{\rm D}}}}^{(E)}$ can be expressed as follows:

[TABLE]

or equivalently, $\smash{\widehat{{\boldsymbol{\rm D}}}}^{(E)}=(\smash{\widehat{p}}_{0})^{-2}\smash{\widehat{{\boldsymbol{\rm D}}}}_{\rm vac}^{(A)}+\smash{\widehat{{\boldsymbol{\rm\chi}}}}$. The corresponding Weyl image is $\smash{\underline{D}^{(E)}_{ab}=(p_{0}){}^{-2}[\underline{D}_{\rm vac}^{(A)}]_{ab}+\underline{\chi}_{ab}}$, or

[TABLE]

where $\smash{\underline{\varepsilon}_{ab}\doteq h_{ab}+(\underline{\chi}_{0})_{ab}}$ serves as the dielectric tensor and $\underline{\chi}_{0}$ is the GO limit of $\underline{\chi}$.

Using the fact that the dispersion operator is defined only up to a constant factor, let us introduce an additional factor $1/(16\pi)$. Then, Eq. (234) becomes

[TABLE]

The quantity $\underline{C}^{(A)}_{ab}$ was derived in Sec. VI.1, and the vector $q$ [Eq. (15)] and the matrices $\ell_{a}$ [Eq. (137)] can be written as

[TABLE]

so they both can be of order $\epsilon_{\lVert}$; hence, $C=O(\epsilon_{\lVert})$. Note that such $C$ is generally non-negligible even for near-Euclidean metrics such as those discussed in Sec. IV.5.3.

From Eq. (196), we obtain

[TABLE]

becomes as the true action-flux density my:amc . As a reminder, the action density $\mathcal{I}\doteq\mathcal{J}^{0}$ can be written as

[TABLE]

where we used $[\underline{D}_{H}^{(E)}]_{ab}E^{b}\approx 0$. Using the latter and Faraday’s law ${\boldsymbol{\rm B}}\approx c{\boldsymbol{\rm k}}\times{\boldsymbol{\rm E}}/\omega$ for the magnetic field ${\boldsymbol{\rm B}}$, one can also rewrite $\mathcal{I}$ in another common form book:stix ,

[TABLE]

One can also show my:amc that the spatial component of the action flux density (239) can be expressed as

[TABLE]

where ${\boldsymbol{\rm S}}$ is the (time- or space-averaged) Poynting vector,

[TABLE]

The GO equation (195) can be written as

[TABLE]

The corresponding density of the wave energy is $\omega\mathcal{I}$, and the density of the energy flux is $\omega{\boldsymbol{\rm\mathcal{J}}}$, as flows from the variational principle my:amc .

VI.3 Stationary waves

For stationary waves with fixed frequency, one has $\smash{\smash{\widehat{p}}_{0}}=-\omega/c$ and $\omega$ can be treated as a constant parameter. Such waves can be studied on the three-dimensional configuration space, namely, the physical space ${\boldsymbol{\rm x}}$. In this case, $n=m=3$ and $\gamma=\boldsymbol{h}$. If there is no spatial dispersion, then $\smash{\smash{\widehat{{\boldsymbol{\rm\chi}}}}={\boldsymbol{\rm\chi}}(\smash{\widehat{{\boldsymbol{\rm x}}}})}$, where the dependence on $\omega$ is assumed but not emphasized. Assuming that the underlying space is flat (which is always the case in laboratory applications), the symbol $\chi^{a}{}_{b}$ in this case is identical to the susceptibility of the homogeneous medium. (Other formulations of the susceptibility have also been proposed in this case, e.g., for light in nondispersive dielectric media my:qdiel and for waves in cold plasmas my:covar ; phd:ruiz17 ; ref:friedland88 .)

These properties can also be extended, approximately, to media with weak spatial dispersion, i.e., such that

[TABLE]

For example, in plasma, $\smash{\underline{\chi}_{ab}^{(k)}}$ is proportional to the temperature, and for EM waves, thermal effects often can be considered as small perturbations book:stix . Although the Weyl symbol $\smash{\underline{\chi}_{ab}^{(k)}}$ is, strictly speaking, different from the corresponding part of the susceptibility in homogeneous medium, the difference is of order $\smash{\epsilon_{\parallel}\underline{\chi}_{ab}^{(k)}}$, so it is much smaller than $\epsilon_{\parallel}$ and thus can be neglected. In other words, one can simply replace the small $\smash{\underline{\chi}_{ab}^{(k)}}$ with its GO limit, which is usually well known book:stix .

In principle, our general method is also applicable to waves in media with strong spatial dispersion. However, $\smash{\smash{\widehat{{\boldsymbol{\rm\chi}}}}}$ may be hard to calculate in that case unless a medium is homogeneous, and no extrapolation of a known homogeneous-medium model is guaranteed to yield the true general $\smash{\smash{\widehat{{\boldsymbol{\rm\chi}}}}}$.

VI.4 Waves in flat space

Let us also discuss the special case when the configuration space is flat and Eq. (244) applies, as in typical laboratory applications. Although the ray-based metric $\boldsymbol{h}$ is still non-Euclidean in this case (unless the reference rays are straight), one can use the flatness of the configuration space in the following sense. In the Euclidean laboratory coordinates $\tilde{x}$, the spatial metric has the form $\smash{\tilde{h}_{ab}=\delta_{ab}}$ and $\smash{\tilde{C}}$ vanishes. (Here and further, the tildes denote that the corresponding quantities are evaluated in the laboratory coordinates. Note that Papers II and III assume a different notation; there, tilded are quantities in the ray-based coordinates, and non-tilded are those in the laboratory coordinates.) This simplifies the expression for the dispersion matrix and makes $D_{H}$ a true tensor (modulo the insignificant corrections caused by weak spatial dispersion discussed in Sec. VI.3). Then, since $U$ and $\Gamma$ are frame-invariant (Sec. V.5), one can calculate them in the laboratory coordinates using

[TABLE]

where both $\delta^{-1}$ and $\tilde{\gamma}$ are omitted for simplicity, since they are unit matrices. This also extends to the quasioptical model. Since calculations in the laboratory coordinates are generally easier than the corresponding calculations in the ray-based coordinates, this simplifies applications of our theory to numerical simulations, as discussed further in Paper II.

Specific applications of this theory, including those to waves in fusion plasmas, will be discussed in future publications. Here, we only note that our formulation relies on the GO ordering and thus is inapplicable to wave transformations near caustics, for example, as in the O–X conversion caused by the density gradient in dense magnetized plasmas ref:preinhaelter73 . In contrast, the O–X conversion caused by the magnetic shear in dilute plasmas my:xo can be modeled naturally, as also demonstrated explicitly in Paper III. Any adiabatic transformations of waves along continuous dispersion curves can also be modeled naturally, specifically, within the single-mode approximation summarized in Sec. V.7.2. (In special cases like the X–B transformation book:stix , when the group velocity becomes zero on a ray while $k$ remains large, one may want to replace the coordinate $\zeta$ with $\tau$, where $\mathrm{d}\tau=\mathrm{d}\zeta/V$, to remove the singularity in the amplitude equation.) In all these cases, our theory is expected to yield results that agree with full-wave modeling up to local errors of order $O(\epsilon_{\parallel})$, because this is the accuracy within which the passive-mode contribution is calculated (Sec. V.4.4).

VII Conclusions

In summary, we propose a quasioptical theory of mode-converting wave beams in inhomogeneous media such as plasma. This includes the following. For any given dispersion operator $\smash{\widehat{D}}$ that governs the original wave field $\Psi$, we explicitly calculate the approximate operator $\smash{\widehat{\mathcal{D}}}$ that governs the wave envelope $\psi$ to the second order in the GO parameter $\epsilon$. Then, we further simplify this envelope operator by assuming that the gradient of $\psi$ transverse to the local group velocity is much larger than the corresponding parallel gradient. This leads to a parabolic differential equation for $\psi$ (“quasioptical equation”) in the basis of the GO polarization vectors [Eq. (146)]. Our main results can be found in Sec. V.6, which includes the general quasioptical equation (205) and its special case (209). Scalar and mode-converting vector beams are described on the same footing (Sec. V.7). We also explain how to apply this model to EM waves considered as a special case (Sec. VI). In the follow-up papers, we report successful quasioptical modeling of radiofrequency wave beams in magnetized plasma based on this theory.

VIII Acknowledgments

The first author (IYD) acknowledges the support and hospitality of the National Institute of Fusion Science, Japan. The authors are also thankful to M. A. Oancea for stimulating discussions regarding the XGO generalization to non-Euclidean metrics. The work was supported by JSPS KAKENHI Grant Number JP17H03514. The work was also supported by the U.S. DOE through Contract No. DE-AC02-09CH11466 and by the Laboratory Directed Research and Development program at Sandia National Laboratories. Sandia National Laboratories is a multimission laboratory managed and operated by National Technology and Engineering Solutions of Sandia, LLC., a wholly owned subsidiary of Honeywell International, Inc., for the U.S. DOE National Nuclear Security Administration under contract DE-NA-0003525. This paper describes objective technical results and analysis. Any subjective views or opinions that might be expressed in the paper do not necessarily represent the views of the U.S. Department of Energy or the United States Government.

Appendix A Summary of selected notations

Here, we present a summary of selected notations used in the main text.

A.1 Basic symbols

• •

$\doteq$ denotes a definition.

• •

$\smash{\widehat{\phantom{A}}}$ denotes an operator.

• •

† denotes either a dual vector [Eqs. (110) and (157)] or an adjoint operator.

• •

⊺ denotes the matrix transpose.

• •

${}^{\sf H}={}^{\intercal*}$ denotes the conjugate matrix transpose.

• •

• is used only in $\Xi^{+}$ and $\bar{\Xi}^{+}$ defined in Eqs. (168).

• •

H and A denote the Hermitian and anti-Hermitian parts of an operator [Eqs. (26)] or those of a matrix (Sec. V.2). When applied to a scalar, H and A denote the real and imaginary parts, correspondingly.

• •

The underline notation $\underline{\phantom{A}}$ is explained in Sec. V.2 and is used for matrices with lower indices only.

• •

$\mathscr{W}_{z}$ and $\mathscr{W}^{-1}$ denote the Wigner–Weyl transform and its inverse (Sec. IV.1.3), respectively.

A.2 Index manipulation

On the $n$-dimensional configuration space $M^{n}$, we assume a general metric tensor $g$ with components $g_{\mu\nu}$. For (co)vector fields on the space (co)tangent to $M^{n}$, the indices are manipulated as usual, namely, via

[TABLE]

(Summation over repeated indices is always assumed unless specified otherwise.) For vector waves, we also introduce an additional vector space that may or may not be the same as the space tangent to $M^{n}$ (Sec. V.1). On that space, an additional metric $\gamma$ is introduced, and the indices are manipulated as follows:

[TABLE]

In special cases, Eqs. (245) and (246) can be equivalent, but that is not a generic situation.

We also introduce the Euclidean metric $\delta$ on the space of amplitude vectors $a$ (Sec. V.4.2),

[TABLE]

For active modes, which are of our primary interest, the mode indices $s$ and $s^{\prime}$ range from 1 to $N$. The distinction between $a^{s}$ and $a_{s}$ is made only for aesthetic reasons (consistency of notation) and can be neglected otherwise.

A.3 Derivatives

For spatial derivatives, we use the notation that is standard, for example, in general relativity book:misner77 . In particular, $X^{\mu}{}_{;\mu}$ is the divergence; namely,

[TABLE]

Here, $g_{\diamond}\doteq|\det g\,|$, and

[TABLE]

For functions of the form $f({\boldsymbol{\rm x}},{\boldsymbol{\rm k}})$, we also introduce the following partial derivatives:

[TABLE]

Since ${\boldsymbol{\rm k}}=\nabla\theta({\boldsymbol{\rm x}})$, one has

[TABLE]

A.4 Inner products

The inner product for scalar waves on $M^{n}$ is

[TABLE]

and for $m$-dimensional vector waves on $M^{n}$,

[TABLE]

In particular, in the ${\boldsymbol{\rm x}}$ representation,

[TABLE]

We also use the dot product

[TABLE]

The inner product on the transverse space $M^{n-1}_{\perp}$ is defined for scalar fields as

[TABLE]

and for the $N$-dimensional “active” parts of the vector fields (Sec. V.4.2) as

[TABLE]

Appendix B Auxiliary calculations

Here we present some auxiliary calculations whose results are used in Sec. V.

B.1 Calculation of $\boldsymbol{\smash{\widehat{K}}}$ and $\boldsymbol{U}({\boldsymbol{\rm x}})$

We start with Eq. (189) for $\smash{\widehat{K}}$. This equation can be rewritten as follows:

[TABLE]

where we used Eq. (52) and introduced

[TABLE]

where, at the end, we used $\Lambda\sim O(\epsilon)$ (Sec. V.4.4), so that the last two terms in Eq. (261) can be neglected.

From Eq. (176), we have $(\Xi^{+})_{,\mu}=\delta^{-1}\Xi^{\sf H}_{,\mu}\gamma+\Xi^{\sf+}\ell_{\mu}$. Using this and also Eq. (140), one obtains

[TABLE]

where the subscript $A$ denotes the anti-Hermitian part, as usual.

B.2 Calculation of $\boldsymbol{U}_{0}({\boldsymbol{\rm x}},{\boldsymbol{\rm k}})$

Let us also derive an alternative expression for $U_{0}$ [Eq. (263)] in terms of the partial derivatives |μ instead of the “full” derivatives ,μ. Using Eq. (59), one obtains $\Xi_{,\mu}=\Xi^{|\nu}k_{\nu,\mu}+\Xi_{|\mu}$. Then, Eq. (263) can be rewritten as follows:

[TABLE]

By differentiating Eq. (181) with respect to $k_{\mu}$, we obtain

[TABLE]

where we used that $\Lambda$ (but not its partial derivatives) is small. Let us multiply this by $\gamma$ to obtain

[TABLE]

where the latter equality is the complex conjugate of the former. Using these and $k_{\nu,\mu}=\theta_{,\nu\mu}=\theta_{,\mu\nu}=O(\epsilon)$, we further obtain, to the leading order,

[TABLE]

where we neglected terms of the second order in $\epsilon$. We also used the fact that ${\Lambda_{s}\textbf{(}{\boldsymbol{\rm x}},{\boldsymbol{\rm k}}({\boldsymbol{\rm x}})\textbf{)}\lesssim O(\epsilon)}$ by definition of active modes and, in addition, $\Lambda_{s}$ are changing slowly, so $\Lambda_{,\mu}\lesssim O(\epsilon^{2})$. (Note that this is only true for $\Lambda_{,\mu}$ but not for $\Lambda_{|\mu}$.) Using Eqs. (266), we also obtain similarly that

[TABLE]

Then,

[TABLE]

As a reminder, the subscript $A$ denotes the anti-Hermitian part, and $\delta^{-1}$ is a unit matrix that only raises the mode index. Equation (269) is similar to the corresponding result derived in Refs. my:qdirac ; my:covar ; phd:ruiz17 . Extra terms are included in those papers in the expressions for $U_{0}$, but they are of the second order in $\epsilon$ and could be neglected.

B.3 Calculation of $\boldsymbol{\smash{\widehat{\mathcal{G}}}}$

Consider Eq. (203) for $\smash{\widehat{\mathcal{G}}}$ derived under the assumption of the quasioptical ordering. Here, we simplify that equation and derive the explicit formula (204). First, note that

[TABLE]

where we used Eq. (139). Also, using Eq. (138) for $\smash{\widehat{\mathscr{D}}}_{H1}$, one obtains

[TABLE]

Here, in the first line, we substituted Eq. (138) for $\smash{\widehat{\mathscr{D}}}_{H1}$ and also Eq. (14) for $\smash{\widehat{p}}_{\mu}$. Then, we only kept the terms involving derivatives of $a$ along the perpendicular direction of the wave beam. The parenthesis in the indices denote symmetrization; namely, for any $A$ and $B$,

[TABLE]

The metric factor $h_{\diamond}$ has been added to keep the operator self-adjoint under the inner product (198). (Accounting for the inhomogeneity of $\boldsymbol{h}$ in $\smash{\widehat{\mathcal{G}}}$ is beyond the accuracy of our theory. However, keeping $\smash{\widehat{\mathcal{G}}}$ self-adjoint is convenient and physically meaningful, for it is known that the exact operator is not responsible for dissipation.) Hence,

[TABLE]

In order to simplify the expression for $\Phi^{\sigma\tilde{\sigma}}$, note that

[TABLE]

where we used $\bar{\Xi}^{+}\Xi=\mathbb{0}$ (and thus $\bar{\Xi}^{+}\Xi^{|\tilde{\sigma}}=-\bar{\Xi}^{+|\tilde{\sigma}}\Xi$) and neglected terms proportional to $\Lambda=O(\epsilon_{\lVert})$. Hence,

[TABLE]

Using Eqs. (174) for $\Xi$ and $\bar{\Xi}$ like we did above, Eq. (178) for $\mathsf{D}_{H}$, and the fact that $\Lambda=O(\epsilon)$, we finally obtain

[TABLE]

where we substituted $\Lambda^{|\sigma}\approx V^{\sigma}$ [Eq. (261)]. The quasioptical approximation implies that $V^{\tilde{\sigma}}$ is close to a scalar matrix, so the commutator $[\Xi^{+}\Xi^{(|\sigma},V^{\tilde{\sigma})}]$ can be neglected. Hence, $\Phi^{\sigma\tilde{\sigma}}\approx\Lambda^{|\sigma\tilde{\sigma}}$.