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

Dodin, I. Y.; Ruiz, D. E.; Yanagihara, K.; Zhou, Y.; Kubo, S.

doi:10.1063/1.5095076

This work opens a series of papers where we develop a general quasi-optical 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 multicomponent 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 ψ to the second order in the geometrical-optics parameter. Then, we further simplify this envelope operator by assuming that the gradient of ψ transverse to the local group velocity is much larger than the corresponding parallel gradient. This leads to a parabolic differential equation for ψ (“quasioptical equation”) on 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 radio frequency wave beams in magnetized plasma based on this theory.

I. INTRODUCTION

A. Motivation

Describing the propagation of waves in inhomogeneous media is a classic problem with a long history.^1–3 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.⁴ Full-wave modeling, which involves solving complete Maxwell's equations, can be impractical at short (cm and mm) wavelengths, especially when multidimensional 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)¹ and include conventional ray tracing,^5–8 complex ray tracing,^9,10 beam tracing,^11–14 and variations thereof.¹⁵ There are also other “quasioptical” models, such as in Refs. 16–19, which 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,²⁰ i.e., mode conversion does not occur.²¹ As a result, today's simulations of mode-converting beams mainly rely on full-wave codes¹⁹ and thus have to compromise the fidelity by reducing the number of dimensions resolved.⁸

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.^22,23 Hence, it is potentially beneficial to approach the problem of quasioptical modeling as a general-physics problem, without restricting it to fusion applications.

B. 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.^24–28 XGO is a theory which calculates the leading-order correction U to the GO dispersion operator of a general vector wave and shows²⁶ 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 affecting the electron motion and (ii) it also governs mode conversion, which appears as a direct analog of spin up–down transitions.²⁹ An examination of the quasioptical algorithms such as those in Refs. 12 and 16 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.³⁰

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 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 many arbitrary 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³¹ and Paper III,³² 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,²⁹ which, for example, is a known problem^33,34 in the Large Helical Device.^35,36

C. Outline

In this first paper of our series, we consider an arbitrary quasimonochromatic multicomponent wave in a weakly inhomogeneous linear medium. Supposing that $\hat{D}$ is some dispersion operator governing the wave dynamics, we simplify $\hat{D}$ and obtain an approximate operator that governs the wave envelope ψ. Then, we derive a parabolic differential equation for ψ (“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 the supplementary material. 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ⁿ with coordinates $x \equiv {x^{0}, x^{1}, \dots, x^{n - 1}}$ and some general metric tensor $g (x)$ ⁠. For simplicity, assume that Mⁿ is diffeomorphic to $ℝ^{n}$ ⁠, i.e., the n-dimensional Euclidean space or pseudo-Euclidean space with the same metric signature as Mⁿ. (For more information on why the diffeomorphism with $ℝ^{n}$ is needed, see Sec. IV A 2.) Suppose that the wave field $Ψ \equiv Ψ (x)$ ⁠, which may have multiple components, is governed by a linear equation with no source terms,

\hat{D} Ψ = 0,

(1)

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

ϵ ≐ λ / L ≪ 1

(2)

(the symbol $≐$ denotes definitions), where λ 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.³⁷ Below, we propose a systematic reduction of Eq. (1) using the smallness of ϵ and eventually obtain a quasioptical model based on this equation. The idea of quasioptical modeling will be formalized later (Secs. IV E and V F).

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

Ψ = \hat{U} ψ, \hat{U} ≐ e^{i θ (x)} .

(3)

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

k (x) ≐ \nabla θ (x)

(4)

is the wave vector identical or close to that predicted by the GO approximation. This implies that the envelope ψ and also k are slow functions [and that the wavelength substituted in Eq. (2) is $λ \approx 2 π / k$ ]. Then, Eq. (1) becomes

\hat{D} ψ = 0,

(5)

where the “envelope dispersion operator” is $\hat{D} ≐ {\hat{U}}^{†} \hat{D} \hat{U}$ (the dagger denotes the adjoint, as usual) or, more explicitly,

\hat{D} = e^{- i θ} \hat{D} e^{i θ} .

(6)

The reference phase θ is treated as a prescribed function. As it will become clear later, knowing θ per se is not actually needed for our purposes; instead, it is k that matters. The latter can be calculated on some “reference rays” using the conventional ray equations,¹

\frac{d x^{α}}{d τ} = \frac{\partial H}{\partial k_{α}}, \frac{d k_{α}}{d τ} = - \frac{\partial H}{\partial x^{α}}

(7)

(τ is any parameter along the ray), which also leads to

H (x, k (x)) = const .

(8)

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

Also note that our approach equally applies to stationary and nonstationary waves. In the case of a stationary wave, we assume that x is a coordinate in physical space (“spatial problem”) and the wave frequency ω serves as a constant parameter. In the case of a nonstationary wave, we assume that x is a coordinate in spacetime (“spacetime problem”), and then, ω is a part of k. In spacetime problems, we assume coordinates such that $x^{0} = ct$ ⁠, where c is the speed of light, t is the time, and the metric signature is $(-, +, +, \dots)$ ⁠; then, $k_{0} = - ω / c$ (in the case of the Minkowski metric, this implies that $k^{0} = ω / c$ ⁠), and so $ω = - \partial_{t} θ$ ⁠, 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 $\hat{D}$ can be asymptotically expanded in ϵ for any $\hat{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. (For recent overviews, see Refs. 1 and 24. Specific key papers include, but are not limited to, Refs. 39–44. In particular, Ref. 43 is one of the first papers which explicitly discuss applications of the Weyl calculus to plasma problems; although also see Refs. 45–47.) In order to present our framework and notation in a self-contained manner, we restate some basics in Sec. IV A and the supplementary material. Also note that our calculations in Sec. IV B are similar in spirit to those in Ref. 43, except that we perform a higher-order expansion and assume a more general metric. We shall also comment throughout this paper on how our equations reproduce other known results in the corresponding limits.

A. Weyl calculus

1. Basic definitions

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

(\hat{A} Ψ) (x) ≐ \int d^{n} x' \sqrt{g_{⋄} (x')} A (x, x') Ψ (x') .

(9)

Here, the integral is taken over $ℝ^{n}$ (and so are all integrals below, up to dimensions), $g_{⋄} ≐ | \det g |$ ⁠, and $A$ is some kernel function that determines $\hat{A}$ ⁠. Also, consider a family of all unitary operators ${\hat{A}}_{u}$ ⁠, which is a subset of all possible $\hat{A}$ ⁠. For a given Ψ, all image functions ${\hat{A}}_{u} Ψ$ are mutually equivalent up to an isomorphism, and so their family ${{\hat{A}}_{u} Ψ}$ can be viewed as a single object, a “state vector” $| Ψ ⟩$ ⁠, which belongs to a Hilbert space $H_{1}^{n}$ with the inner product

⟨ Ψ | Φ ⟩ ≐ \int d^{n} x \sqrt{g_{⋄} (x)} Ψ^{*} (x) Φ (x) .

(10)

Then, Eq. (9) can be viewed as the “x representation” of $\hat{A}$ ⁠, while the operator itself can be understood more generally as a transformation of $| Ψ ⟩$ ⁠, i.e., of the whole family ${{\hat{A}}_{u} Ψ}$ ⁠. Using this invariant notation, one can formulate a machinery, called the Weyl calculus,⁴⁸ which allows efficient asymptotic approximation of operators using the smallness of ϵ. Below, we overview the key theorems of the Weyl calculus which are used in our paper. Readers who are interested in details and proof of these theorems can find them in the supplementary material.

2. Coordinate and momentum operators

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

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

(11)

\hat{p} = {{\hat{p}}_{0}, {\hat{p}}_{1}, \dots, {\hat{p}}_{n - 1}}

(12)

such that the x representations of ${\hat{x}}^{μ}$ and ${\hat{p}}_{μ}$ are as follows:

{\hat{x}}^{μ} Ψ = x^{μ} Ψ,

(13)

{\hat{p}}_{μ} Ψ = - i g_{⋄}^{- 1 / 4} \partial_{μ} (g_{⋄}^{1 / 4} Ψ) .

(14)

Here, $\partial_{μ} ≐ \partial / \partial x^{μ}, x^{μ}$ and p_μ are the corresponding eigenvalues, and the factors $g_{⋄}^{\pm 1 / 4}$ are introduced to keep $\hat{p}$ self-adjoint under the inner product (10).⁴⁹ (This would not be the case if Mⁿ were not diffeomorphic to $ℝ^{n}$ ⁠.⁵⁰) Since

{\hat{p}}_{ν} = - i \partial_{ν} - i q_{ν} (x), q_{ν} ≐ \frac{1}{4} \partial_{ν} (ln g_{⋄}),

(15)

and q_ν commutes with ${\hat{x}}^{μ}$ ⁠, the usual commutation relation $[{\hat{x}}^{μ}, {\hat{p}}_{ν}] = i δ_{ν}^{μ}$ is satisfied.

Let us consider the eigenvectors $| x ⟩$ and $| p ⟩$ of the coordinate and momentum operators, which are defined as

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

(16)

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

⟨ x_{1} | x_{2} ⟩ = G (x_{1}, x_{2}) δ (x_{1} - x_{2}),

(17)

⟨ p_{1} | p_{2} ⟩ = \bar{G} (p_{1}, p_{2}) δ (p_{1} - p_{2}) .

(18)

Here, we introduced

G (x_{1}, x_{2}) ≐ {[g_{⋄} (x_{1}) g_{⋄} (x_{2})]}^{- 1 / 4},

(19)

\bar{G} (p_{1}, p_{2}) ≐ {[{\bar{g}}_{⋄} (p_{1}) {\bar{g}}_{⋄} (p_{2})]}^{- 1 / 4} .

(20)

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

Ψ (x) = ⟨ x | Ψ ⟩, ψ (x) = ⟨ x | ψ ⟩ .

(21)

The p representations of these state vectors are introduced similarly as $⟨ p | Ψ ⟩$ and $⟨ p | ψ ⟩$ ⁠. One can also show that $⟨ x_{1} | \hat{A} | x_{2} ⟩ = A (x_{1}, x_{2})$ and

⟨ x | p ⟩ = {⟨ p | x ⟩}^{*} = \frac{exp (i p \cdot x),}{{(2 π)}^{n / 2} {[g_{⋄} (x) {\bar{g}}_{⋄} (p)]}^{1 / 4}} .

(22)

Here, $p \cdot x ≐ p_{μ} x^{μ}$ ⁠, and summation over repeated indices is assumed from now on.

3. Wigner–Weyl transform

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

{\hat{Δ}}_{z} ≐ \int d^{n} s G (x, s) | x - s / 2 ⟩ ⟨ x + s / 2 | e^{- i p \cdot s},

(23)

G (x, s) ≐ {[g_{⋄} (x - s / 2) g_{⋄} (x + s / 2)]}^{1 / 4} .

(24)

(Note that G = 1 if the x space is Euclidean or pseudo-Euclidean.) Using ${\hat{Δ}}_{z}$ ⁠, we define the Wigner–Weyl transform $W_{z} : \hat{A} \mapsto A$ ⁠, which maps a given operator $\hat{A}$ on $H_{1}^{n}$ to a function A (“Weyl symbol”) on the z space. Specifically, the Weyl symbol of a given operator $\hat{A}$ is $A (x, p) ≐ tr ({\hat{Δ}}_{z} \hat{A})$ (“tr” stands for trace); i.e.,

A (x, p) ≐ \int d^{n} s G (x, s) ⟨ x + s / 2 | \hat{A} | x - s / 2 ⟩ e^{- i p \cdot s} .

(25)

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 $\hat{A}$ as $\hat{A} = {\hat{A}}_{H} + i {\hat{A}}_{A}$ ⁠, where the subscripts denote the Hermitian and anti-Hermitian parts,

{\hat{A}}_{H} ≐ \frac{1}{2} (\hat{A} + {\hat{A}}^{†}), {\hat{A}}_{A} ≐ \frac{1}{2 i} (\hat{A} - {\hat{A}}^{†}) .

(26)

Both ${\hat{A}}_{H}$ and ${\hat{A}}_{A}$ (not to be confused with $i {\hat{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 $W^{- 1} : A \mapsto \hat{A}$ ⁠, which maps a given function A to the corresponding operator $\hat{A}$ via

\hat{A} = \frac{1}{{(2 π)}^{n}} \int d^{n} x d^{n} p A (x, p) {\hat{Δ}}_{z} .

(27)

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

f (\hat{x}) \Leftrightarrow f (x), f (\hat{p}) \Leftrightarrow f (p) .

(28)

However, the Weyl symbols of operators that cannot be represented as $f_{1} (\hat{x}) + f_{2} (\hat{p})$ are generally more complicated. In particular, one can show that (see the supplementary material)

f (\hat{x}) {\hat{p}}_{α} \Leftrightarrow p_{α} f (x) + \frac{i}{2} \partial_{α} f (x),

(29)

{\hat{p}}_{α} f (\hat{x}) \Leftrightarrow p_{α} f (x) - \frac{i}{2} \partial_{α} f (x),

(30)

and also

\begin{matrix} {\hat{p}}_{α} f (\hat{x}) {\hat{p}}_{β} \Leftrightarrow p_{α} p_{β} f (x) + \frac{1}{4} \partial_{α β}^{2} f (x) \\ + \frac{i}{2} p_{α} \partial_{β} f (x) - \frac{i}{2} p_{β} \partial_{α} f (x) . \end{matrix}

(31)

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

B. Approximate $\hat{D}$

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

\hat{D} | Ψ ⟩ = 0.

(32)

Accordingly, Eq. (5) becomes

\hat{D} | ψ ⟩ = 0,

(33)

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

\hat{D} = e^{- i θ (\hat{x})} \hat{D} e^{i θ (\hat{x})} .

(34)

We shall now approximate $\hat{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 ϵ [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 $\hat{D}$ as $D (x, p) = e^{- i θ (x)} ⋆ D (x, p) ⋆ e^{i θ (x)}$ ⁠, where D is the Weyl symbol of $\hat{D}$ and $⋆$ is the Moyal product (see the supplementary material). By expanding $⋆$ in the GO parameter, one can then obtain an approximate symbol of $\hat{D}$ formally. But here we adopt a different approach, which is more transparent. We start by expressing $\hat{D}$ through D explicitly,

\hat{D} = \frac{1}{{(2 π)}^{n}} \int d^{n} x d^{n} p d^{n} s G (x, s) | x - s / 2 ⟩ e^{- i p \cdot s} D (x, p) ⟨ x + s / 2 |,

(35)

where G is a metric factor given by Eq. (24). Using the fact that $θ (\hat{x}) | x \pm s / 2 ⟩ = θ (x \pm s / 2) | x \pm s / 2 ⟩$ ⁠, one can rewrite Eq. (34) as follows:

\hat{D} = \frac{1}{{(2 π)}^{n}} \int d^{n} x d^{n} p d^{n} s G (x, s) e^{i [θ (x + s / 2) - θ (x - s / 2) - p \cdot s]} | x - s / 2 ⟩ D (x, p) ⟨ x + s / 2 | .

(36)

Consider a formal Taylor expansion of the reference phase θ in s,

θ (x \pm s / 2) = θ (x) \pm \frac{s^{α}}{2} \frac{\partial θ (x)}{\partial x^{α}} + \frac{1}{2} \frac{s^{α} s^{β}}{4} \frac{\partial^{2} θ (x)}{\partial x^{α} \partial x^{β}} \pm \frac{1}{6} \frac{s^{α} s^{β} s^{γ}}{8} \frac{\partial^{3} θ (x)}{\partial x^{α} \partial x^{β} \partial x^{γ}} + \dots

(37)

This gives

θ (x + s / 2) - θ (x - s / 2) = s^{α} \frac{\partial θ (x)}{\partial x^{α}} + \frac{s^{α} s^{β} s^{γ}}{24} \frac{\partial^{3} θ (x)}{\partial x^{α} \partial x^{β} \partial x^{γ}} + \dots = k (x) \cdot s + \frac{s^{α} s^{β} s^{γ}}{24} \frac{\partial^{2} k_{γ} (x)}{\partial x^{α} \partial x^{β}} + \dots,

(38)

and so

\begin{matrix} e^{i [θ (x + s / 2) - θ (x - s / 2) - p \cdot s]} = [1 + \frac{i}{24} s^{α} s^{β} s^{γ} \frac{\partial^{2} k_{γ} (x)}{\partial x^{α} \partial x^{β}} + \dots] e^{i [k (x) - p] \cdot s} \\ = [1 + \frac{1}{24} \frac{\partial^{2} k_{γ} (x)}{\partial x^{α} \partial x^{β}} \frac{\partial^{3}}{\partial p_{α} \partial p_{β} \partial p_{γ}} + \dots] e^{i [k (x) - p] \cdot s} . \end{matrix}

(39)

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

\begin{matrix} \hat{D} = \frac{1}{{(2 π)}^{n}} \int d^{n} x d^{n} p d^{n} s G (x, s) {[1 + \frac{1}{24} \frac{\partial^{2} k_{γ} (x)}{\partial x^{α} \partial x^{β}} \frac{\partial^{3}}{\partial p_{α} \partial p_{β} \partial p_{γ}} + \dots] e^{i [k (x) - p] \cdot s}} | x - s / 2 ⟩ D (x, p) ⟨ x + s / 2 | \\ = \frac{1}{{(2 π)}^{n}} \int d^{n} x d^{n} p d^{n} s G (x, s) e^{i [k (x) - p] \cdot s} | x - s / 2 ⟩ [D (x, p) - \frac{1}{24} \frac{\partial^{2} k_{γ} (x)}{\partial x^{α} \partial x^{β}} \frac{\partial^{3} D (x, p)}{\partial p_{α} \partial p_{β} \partial p_{γ}} + \dots] ⟨ x + s / 2 | \\ = \frac{1}{{(2 π)}^{n}} \int d^{n} x d^{n} p d^{n} s G (x, s) e^{- i p \cdot s} | x - s / 2 ⟩ D' (x, k (x) + p) ⟨ x + s / 2 |, \end{matrix}

(40)

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

D' (x, p) ≐ D (x, p) - \frac{1}{24} \frac{\partial^{2} k_{γ} (x)}{\partial x^{α} \partial x^{β}} \frac{\partial^{3} D (x, p)}{\partial p_{α} \partial p_{β} \partial p_{γ}} + \dots

(41)

The first term on the right-hand side of Eq. (41) is O(1) because at $ϵ \to 0$ ⁠, it becomes the GO dispersion function, which is assumed order-one.⁵¹ (At nonzero ϵ, the Weyl symbol D generally does not coincide with the dispersion function that governs the wave in a homogeneous medium; see Sec. IV A 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' (x, p)$ at p close to k. Then,

\frac{\partial^{2} k_{γ}}{\partial x^{α} \partial x^{β}} \frac{\partial^{3} D}{\partial p_{α} \partial p_{β} \partial p_{γ}} \sim \frac{k}{L^{2}} \frac{D}{k^{3}} \sim ϵ^{2} D,

(42)

meaning that the second term in Eq. (41) is roughly ϵ² times smaller than the first one. This fact will be used below.

The inverse Wigner–Weyl transform will turn the coordinate p into the operator $\hat{p}$ ⁠. The latter will act on the wave envelope [Eq. (33)], which is considered slow in the coordinate representation, and so $\hat{p} | ψ ⟩ = O (ϵ)$ ⁠. In this sense, $p = O (ϵ)$ ⁠, and so by Taylor-expanding $D'$ in p, we are effectively expanding $\hat{D}$ in ϵ. It is sufficient for our purposes to adopt the second-order expansion,

\hat{D} \approx \frac{1}{{(2 π)}^{n}} \int d^{n} x d^{n} p d^{n} s G (x, s) e^{- i p \cdot s} | x - s / 2 ⟩ [D' (x) + p_{μ} V^{μ} (x) + \frac{1}{2} p_{μ} p_{ν} Θ^{μ ν} (x)] ⟨ x + s / 2 |,

(43)

where we introduced $D' (x) ≐ D' (x, k (x))$ and

V^{μ} (x) ≐ \frac{\partial D' (x, k (x))}{\partial k_{μ}},

(44)

Θ^{μ ν} (x) ≐ \frac{\partial^{2} D' (x, k (x))}{\partial k_{μ} \partial k_{ν}} .

(45)

By properties of the inverse Wigner–Weyl transform $W^{- 1}$ (Sec. IV A 3), $\hat{D}$ can be written as follows:

\begin{matrix} \hat{D} \approx W^{- 1} [D' (x)] + W^{- 1} [p_{μ} V^{μ} (x)] + \frac{1}{2} W^{- 1} [p_{μ} p_{ν} Θ^{μ ν} (x)] . \end{matrix}

(46)

The inverse Wigner–Weyl transforms here can be calculated using Eqs. (28)–(31) and $Θ^{μ ν} = Θ^{ν μ}$ ⁠. This leads to

\begin{matrix} \hat{D} \approx D' (\hat{x}) + \frac{1}{2} [{\hat{p}}_{μ} V^{μ} (\hat{x}) + V^{μ} (\hat{x}) {\hat{p}}_{μ}] \\ + \frac{1}{2} {\hat{p}}_{μ} Θ^{μ ν} (\hat{x}) {\hat{p}}_{ν} - \frac{1}{8} (\partial_{μ ν}^{2} Θ^{μ ν}) (\hat{x}) . \end{matrix}

(47)

One can further simplify Eq. (47) as follows. Suppose that the medium, including the metric, is characterized by some parameters $P$ ⁠. Within the accuracy of the theory developed in this paper, we shall neglect terms involving $\partial^{2} P$ and $(\partial P) (\partial ψ)$ while generally retaining terms such as $\partial^{2} ψ$ (as to be explained Sec. IV E 2), where ∂ denotes a generic spatial derivative. Hence, we can ignore the difference between D and $D'$ ⁠; we can also omit $\partial_{μ ν}^{2} Θ^{μ ν}$ in Eq. (47), and so $\hat{D}$ can be expressed as follows:

\hat{D} \approx D (\hat{x}) + \frac{1}{2} [{\hat{p}}_{μ} V^{μ} (\hat{x}) + V^{μ} (\hat{x}) {\hat{p}}_{μ}] + \frac{1}{2} {\hat{p}}_{μ} Θ^{μ ν} (\hat{x}) {\hat{p}}_{ν} .

(48)

Also, $V^{μ}$ and $Θ^{μ ν}$ can be calculated through the zeroth-ϵ limit of D, denoted D₀, which is just the dispersion function of the homogeneous medium,

V^{μ} (x) \approx \frac{\partial D_{0} (x, k (x))}{\partial k_{μ}}, Θ^{μ ν} (x) \approx \frac{\partial^{2} D_{0} (x, k (x))}{\partial k_{μ} \partial k_{ν}} .

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

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

X_{μ} = g_{μ ν} X^{ν}, X^{μ} = g^{μ ν} X_{ν},

(49)

where the matrix $g^{μ ν}$ is the inverse of $g_{μ ν}$ ⁠. Similar rules apply to manipulating tensors, as usual.

C. Approximate envelope equation

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

[{\hat{p}}_{μ} f (\hat{x})] ψ = - i g_{⋄}^{- 1 / 4} \partial_{μ} [g_{⋄}^{1 / 4} f (x) ψ],

(50)

where g and ψ are also functions of x. Then,

\begin{matrix} \hat{D} ψ \approx D ψ - i V^{μ} \partial_{μ} ψ - \frac{i}{2} {V^{μ}}_{; μ} ψ - \frac{1}{2 g_{⋄}^{1 / 4}} \partial_{μ} [Θ^{μ ν} \partial_{ν} (g_{⋄}^{1 / 4} ψ)], \end{matrix}

(51)

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

{X^{μ}}_{; μ} ≐ \frac{1}{\sqrt{g_{⋄}}} \frac{\partial}{\partial x^{μ}} (\sqrt{g_{⋄}} X^{μ}) .

(52)

Except for Dψ, 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 only the zeroth-order dispersion function D₀ but also the leading-order correction $D - D_{0} = O (\partial P)$ that may appear for inhomogeneous media], even though this may not be obvious directly from the definition of the Weyl transform (Sec. IV A). To emphasize this fact, and also to ensure consistency of notation with the vector-wave theory to be discussed in Secs. V and VI, we denote this covariant function as $D$ ⁠, and we also adopt the notation

D_{H} ≐ Re D, D_{A} ≐ Im D .

(53)

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

({\hat{D}}_{H} + i {\hat{D}}_{A}) ψ = 0.

(54)

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

\begin{matrix} {\hat{D}}_{H} ψ = D_{H} ψ - i V^{μ} ψ_{, μ} - \frac{i}{2} {V^{μ}}_{; μ} ψ - \frac{1}{2 g_{⋄}^{1 / 4}} {[{D_{H}}^{| μ ν} {(g_{⋄}^{1 / 4} ψ)}_{, ν}]}_{, μ}, \end{matrix}

(55)

\begin{matrix} {\hat{D}}_{A} ψ = D_{A} ψ - i V_{A}^{μ} ψ_{, μ} - \frac{i}{2} {V_{A}^{μ}}_{; μ} ψ - \frac{1}{2 g_{⋄}^{1 / 4}} {[D_{A}^{| μ ν} {(g_{⋄}^{1 / 4} ψ)}_{, ν}]}_{, μ} \end{matrix}

(56)

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

V^{μ} ≐ V_{H}^{μ} = D_{H}^{| μ},

(57)

where we assume a standard notation

f_{, μ} ≐ \frac{\partial f}{\partial x^{μ}}, f_{, μ ν} ≐ \frac{\partial^{2} f}{\partial x^{μ} \partial x^{ν}} .

(58)

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

f_{, μ} (x, k (x)) = f_{| μ} + f^{| ν} k_{ν | μ} .

(59)

The true partial derivatives are introduced as follows:

f_{| μ} ≐ \frac{\partial f (x, k)}{\partial x^{μ}}, f_{| μ ν} ≐ \frac{\partial^{2} f (x, k)}{\partial x^{μ} \partial x^{ν}},

(60)

f^{| μ} ≐ \frac{\partial f (x, k)}{\partial k_{μ}}, f^{| μ ν} ≐ \frac{\partial^{2} f (x, k)}{\partial k_{μ} \partial k_{ν}} .

(61)

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

k_{ν | μ} = k_{ν, μ} = k_{μ, ν} = θ_{, μ ν} = θ_{, ν μ} .

(62)

D. Leading-order approximation: Geometrical optics

For the envelope approximation to hold, i.e., for $ψ (x)$ to remain a slow function, $k (x)$ must be chosen such that $D (x, k (x)) ψ ≲ O (ϵ)$ ⁠. [The other terms in Eq. (5) are automatically small on the account of Eq. (2).] Let us adopt the usual GO ordering⁵²

D_{H} = O (1), D_{A} = O (ϵ),

(63)

or more rigorously, $D_{A} ≲ O (ϵ)$ ⁠; i.e., $D_{A}$ much smaller than O(ϵ) is also allowed.⁵³ Then, one can define $k (x)$ such that

D_{H} (x, k (x)) = 0.

(64)

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

D_{H} (x^{(0)}, k^{(0)}) = 0.

(65)

To the first order in ϵ, Eq. (5) is

V^{μ} ψ_{, μ} + \frac{1}{2} {V^{μ}}_{; μ} ψ = D_{A} ψ,

(66)

and as a corollary,

{J^{μ}}_{; μ} = - 2 D_{A} | ψ |^{2}, J^{μ} ≐ - V^{μ} | ψ |^{2} .

(67)

At zero $D_{A}$ ⁠, the dynamics is conservative, and Eq. (67) reflects conservation of the wave action, or quanta.⁵⁴ Accordingly, $D_{A}$ determines the dissipation rate, $J^{μ}$ can be identified (at least up to a constant factor) as the action flux density, and V^μ 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.^1,54 Below, we extend them by retaining the second-order terms neglected in Eq. (66).

E. Quasioptical model

1. Ray-based coordinates

Let us start by introducing ray-based coordinates as follows.⁵⁵ Suppose multiple rays launched with different $x^{(0)}$ within the beam. Their trajectories can be found using the ray equations as discussed in Sec. IV D; this determines the path ζ along each ray as a function of the initial coordinate. We treat ζ as the longitudinal coordinate along the wave beam, and each of its isosurfaces is considered as the transverse space $M_{⊥}^{n - 1} (ζ)$ ⁠. The coordinates on this space, $ϱ \equiv {ϱ^{1}, ϱ^{2}, \dots, ϱ^{n - 1}}$ ⁠, can be introduced arbitrarily (yet they will be specified below). Then, we define n – 1 independent vector fields $e_{σ}$ via

e_{σ} ≐ \partial x / \partial ϱ^{σ}

(68)

and another vector $e_{0}$ such that it is linearly independent of all $e_{σ}$ ⁠. Hence, any $d x$ can be decomposed as follows:

d x = e_{0} d ζ + e_{σ} d ϱ^{σ},

(69)

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

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

e^{μ} \cdot e_{ν} = δ_{ν}^{μ},

(70)

and so we adopt $e^{0} ≐ \nabla ζ$ ⁠. Then, the general form of $e_{0}$ is

e_{0} = \pm α^{2} e^{0} + β,

(71)

where the first sign is determined by the metric signature and $α ≐ | e^{0} \cdot e^{0} |^{- 1 / 2}, β ≐ e_{σ} β^{σ}$ ⁠, and $β^{σ}$ are arbitrary coefficients. This leads to the following metric representation in the coordinates ${ζ, ϱ}$ ⁠:

g = (\begin{matrix} \pm α^{2} + β^{T} h β & {(h β)}^{T} \\ h β & h \end{matrix}) .

(72)

Here, h is the matrix with elements $h_{σ \tilde{σ}} ≐ e_{σ} \cdot e_{\tilde{σ}}$ and $^{T}$ denotes the transposition. By using a known theorem for the determinant of a block matrix,⁵⁶ one obtains

g_{⋄} = α^{2} h_{⋄}, h_{⋄} ≐ | \det h | .

(73)

Let us also introduce the transverse projection of $d ϱ$ ⁠,

d ϱ_{⊥} ≐ (1 \mp α^{2} e^{0} e^{0}) d ϱ, d ϱ_{⊥}^{σ} = d ϱ^{σ} + β^{σ} d ζ,

(74)

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

d x \cdot d x = \pm α^{2} {(d ζ)}^{2} + h_{σ \tilde{σ}} d ϱ_{⊥}^{σ} d ϱ_{⊥}^{\tilde{σ}},

(75)

and so $h_{σ \tilde{σ}}$ serves as the transverse metric. Below, we assume $β = 0$ ⁠, and so $ϱ_{⊥}^{σ}$ and $ϱ^{σ}$ do not need to be distinguished, and we can define the inner product on $M_{⊥}^{n - 1}$ as follows:

{⟨ ψ | ϕ ⟩}^{⊥} ≐ \int d^{n - 1} ϱ \sqrt{h_{⋄} (ζ, ϱ)} ψ^{*} (ζ, ϱ) ϕ (ζ, ϱ) .

(76)

We also choose transverse coordinates specifically such that $ϱ = const$ ⁠. In other words, the coordinate $ϱ$ 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} (0)$ ⁠. Then, the transverse group velocity is zero,

V^{σ} \equiv e^{σ} \cdot \frac{d x}{d τ} = \frac{d ϱ^{σ}}{d τ} = 0,

(77)

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

2. Quasioptical equation

Suppose a wave beam such that its longitudinal scale $L_{| |}$ ⁠, which is defined via $ψ_{, ζ} \sim ψ / L_{| |}$ ⁠, is much larger than its perpendicular scale $L_{⊥}$ ⁠, which is defined via $ψ_{, σ} \sim ψ / L_{⊥}$ and may or may not be the beam global width. (We assume the notation $ψ_{, ζ} ≐ \partial ψ / \partial ζ$ and $ψ_{, σ} ≐ \partial ψ / \partial ϱ^{σ}$ ⁠.) We introduce two GO parameters,

ϵ_{| |} ≐ λ / L_{| |}, ϵ_{⊥} ≐ λ / L_{⊥}, ϵ_{| |} ≪ ϵ_{⊥}

(78)

[so the original parameter (2) is $ϵ = ϵ_{⊥}$ ], and we shall neglect terms smaller than $O (ϵ_{| |} ϵ_{⊥})$ from now on. Also, we require that either α is chosen as independent of $ϱ$ ⁠, or, more generally, the relative variation of α across the wave beam is small enough such that $α_{, σ}$ is negligible. We also assume

{(D_{H})}_{| μ} ≲ O (ϵ_{| |}), g_{α β, μ} ≲ O (ϵ_{| |}),

(79)

where the derivatives can be taken over the longitudinal or transverse coordinates. (Having $g_{α β, μ} \sim 1 / L_{| |}$ is indeed common in typical applications.³⁰)

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

{\hat{D}}_{H} ψ = - iV ψ_{, ζ} - \frac{i}{2 \sqrt{g_{⋄}}} {(\sqrt{g_{⋄}} V)}_{, ζ} ψ + \hat{G} ψ .

(80)

Here, $V ≐ V^{0}$ ⁠, and we used $V^{μ} ψ_{, μ} = V ψ_{, ζ}$ due to Eq. (77). The operator $\hat{G}$ is self-adjoint under the inner product (76) and given by

\hat{G} ψ ≐ - \frac{1}{2 h_{⋄}^{1 / 4}} {[{D_{H}}^{| σ \tilde{σ}} {(h_{⋄}^{1 / 4} ψ)}_{, \tilde{σ}}]}_{, σ} .

(81)

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

{\hat{D}}_{A} ψ = D_{A} - i V_{A}^{σ} ψ_{, σ} - \frac{i}{2 \sqrt{h_{⋄}}} {(\sqrt{h_{⋄}} V_{A}^{σ})}_{, σ} ψ,

(82)

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

iV ψ_{, ζ} + \frac{i}{2 \sqrt{g_{⋄}}} {(\sqrt{g_{⋄}} V)}_{, ζ} ψ - \hat{G} ψ = i {\hat{D}}_{A} ψ,

(83)

and, as a corollary,

\frac{1}{α} \frac{d}{d ζ} {⟨ ψ | α V | ψ ⟩}^{⊥} = 2 {⟨ ψ | {\hat{D}}_{A} | ψ ⟩}^{⊥},

(84)

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

{⟨ ψ | α V | ψ ⟩}^{⊥} = const .

(85)

Equation (84) indicates that ${\hat{D}}_{A}$ cannot be much larger than $\partial / \partial ζ$ ⁠, and so the correct scaling to assume for the dissipation term is $D_{A} ≲ O (ϵ_{| |})$ ⁠. Then, Eq. (83) leads to the “quasioptical scaling”

ϵ_{| |} \sim ϵ_{⊥}^{2} .

(86)

Hence, the difference between ${\hat{D}}_{A}$ and $D_{A}$ in Eq. (82) is $≲ O (ϵ^{3})$ ⁠, and so it is beyond the accuracy of our theory, which is $O (ϵ^{2})$ ⁠. For this reason, we henceforth adopt⁵⁷

{\hat{D}}_{A} = D_{A},

(87)

which will also simplify our notation. [Also, as a reminder, $\partial_{μ} g_{α β} ≲ O (ϵ_{| |}) \sim ϵ_{⊥}^{2}$ ⁠, and the parameters of the medium vary similarly.] This leads to

iV ψ_{, ζ} + \frac{i}{2 \sqrt{g_{⋄}}} {(\sqrt{g_{⋄}} V)}_{, ζ} ψ - \hat{G} ψ = i D_{A} ψ .

(88)

Equation (88) is the 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 ζ), Eq. (88) is much easier to solve than Eq. (5) with the original expressions for ${\hat{D}}_{H}$ and ${\hat{D}}_{A}$ given by Eqs. (55) and (56).

3. Simplified equations

A typical metric of interest for quasioptical modeling has the form $g_{μ ν} = δ_{μ ν} + R_{μ ν σ} ϱ^{σ}$ with $R_{μ ν σ} = O (L_{| |}^{- 1})$ ⁠.³⁰ This satisfies the assumed scaling (79) and corresponds to $ln g_{⋄} = O (L_{⊥} / L_{| |})$ ⁠. Then,

{(ln g_{⋄})}_{, ζ} \sim L_{⊥} / L_{| |}^{2} = O (ϵ_{| |}^{3 / 2}),

(89)

which is negligible. [In contrast, ${(ln g_{⋄})}_{, σ} = O (ϵ_{| |})$ ⁠.] The derivatives of $g_{⋄}$ in Eq. (81) are negligible too. Hence, the metric factors in the quasioptical equation (88) 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. (88) becomes

iV ψ_{, ζ} + \frac{i}{2} V_{, ζ} ψ + \frac{1}{2} {({D_{H}}^{| σ \tilde{σ}} ψ_{, \tilde{σ}})}_{, σ} = i D_{A} ψ .

(90)

Using the variable transformation $ϕ ≐ \sqrt{V} ψ$ ⁠, this equation can be further simplified as follows:

i ϕ_{, ζ} + \frac{1}{2} {(Σ^{σ \tilde{σ}} ϕ_{, \tilde{σ}})}_{, σ} = i ϒ ϕ,

(91)

which is just a dissipative Schrödinger equation with

Σ^{σ \tilde{σ}} ≐ {D_{H}}^{| σ \tilde{σ}} / V, ϒ ≐ D_{A} / V .

(92)

[In Eq. (91), we used that $V_{, σ} ϕ_{, \tilde{σ}}$ is negligible compared to $V ϕ_{, \tilde{σ} σ} = O (ϵ_{⊥})$ ⁠.] Also, Eq. (84) becomes

\frac{d}{d ζ} {⟨ ϕ | ϕ ⟩}^{⊥} = 2 {⟨ ϕ | ʋ | ϕ ⟩}^{⊥} .

(93)

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

D_{H} (x, k_{σ}, k_{0} (x, k_{σ})) = 0.

(94)

By differentiating this with respect to $k_{σ}$ ⁠, we obtain ${k_{0}}^{| σ} = - D_{H}^{| σ} / D_{H}^{| ζ}$ ⁠. By differentiating the latter formula once again, now with respect to $k_{\tilde{σ}}$ ⁠, we also obtain

\begin{matrix} {k_{0}}^{| σ \tilde{σ}} = - \frac{1}{D_{H}^{| ζ}} (D_{H}^{| ζ ζ} {k_{0}}^{| σ} {k_{0}}^{| \tilde{σ}} + D_{H}^{| ζ \tilde{σ}} {k_{0}}^{| σ} + D_{H}^{| ζ σ} {k_{0}}^{| \tilde{σ}} + D_{H}^{| σ \tilde{σ}}) . \end{matrix}

(95)

Recall that $D_{H}^{| σ} = V^{σ}$ and $V^{σ} = 0$ [Eq. (77)]. Then, ${k_{0}}^{| σ} = 0$ ⁠, and by comparing Eq. (95) with Eq. (92) and using Eq. (57), one finds that ${k_{0}}^{| σ \tilde{σ}} = - Σ^{σ \tilde{σ}}$ ⁠. Thus, Eq. (91) can be rewritten as

i ϕ_{, ζ} - \frac{1}{2} {({k_{0}}^{| σ \tilde{σ}} ϕ_{, \tilde{σ}})}_{, σ} = i ϒ ϕ .

(96)

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

V. VECTOR WAVES

A. Hilbert space for vector waves

Now, let us generalize the above results to vector waves. Suppose an m-component wave field $Ψ \equiv Ψ (x)$ ⁠, so

Ψ = (\begin{matrix} Ψ^{1} \\ Ψ^{2} \\ ⋮ \\ Ψ^{m} \end{matrix}), | Ψ ⟩ = (\begin{matrix} | Ψ^{1} ⟩ \\ | Ψ^{2} ⟩ \\ ⋮ \\ | Ψ^{m} ⟩ \end{matrix}) .

(97)

Here, $| Ψ ⟩$ is a vector on some Hilbert space $H_{m}^{n}$ for which we assume the following general inner product:

{⟨ Ψ | Φ ⟩}_{m} ≐ \int d^{n} x \sqrt{g_{⋄} (x)} γ_{ab} (x) Ψ^{a *} (x) Φ^{b} (x) .

(98)

Here, the Latin indices that characterize the field components span from 1 to m. (This is in contrast to the Greek indices introduced earlier, which characterize the components of x and span from 0 to n – 1.) The matrix $γ (x)$ can be any m × 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, Ref. 27 describes a six-dimensional field on a three-dimensional spacetime; also see Sec. VI B.) This means that the space to which $Ψ (x)$ belongs is not necessarily the space tangent to Mⁿ. With that said, having $γ = g$ is possible as a special case; see Secs. VI A and VI C.

Using γ as a metric, we introduce the standard rules for manipulating the Latin indices,

Ψ_{a} ≐ γ_{ab} Ψ^{b}, Ψ^{a} = γ^{ab} Ψ_{b},

(99)

where $γ^{ab}$ are the elements of $γ^{- 1}$ ⁠. In particular, in the x representation, one has

Ψ^{†} (x) ≐ ⟨ Ψ | x ⟩ = (Ψ_{1}^{*}, Ψ_{2}^{*}, \dots, Ψ_{m}^{*}),

(100)

and we shall also use the following local dot product:

Ψ \cdot Φ ≐ Ψ^{†} (x) Φ (x) = Ψ_{a}^{*} (x) Φ^{a} (x) .

(101)

Under this dot product, a matrix A with mixed-index elements ${A^{a}}_{b}$ is self-adjoint [ $(A Ψ) \cdot Φ = Ψ \cdot (A Φ)$ ] 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 γ is Euclidean or pseudo-Euclidean.

B. Operators on $H_{m}^{n}$

Any operator $\hat{A}$ on $H_{m}^{n}$ can be understood as an m × m matrix with elements ${\hat{A}}^{a}_{b}$ which are operators on $H_{1}^{n}$ ⁠. For any given ${\hat{A}}^{a}_{b}$ ⁠, we introduce ${\underline{\hat{A}}}_{ab} ≐ γ_{ac} {\hat{A}}^{c}_{b}$ ⁠, which is also an operator on $H_{1}^{n}$ ⁠. Correspondingly, two types of Weyl images can be defined, namely, $W_{z} [{\hat{A}}^{a}_{b}]$ and $W_{z} [{\underline{\hat{A}}}_{ab}]$ ⁠. Below, we assume the notation

{\underline{A}}_{ab} ≐ W_{z} [{\underline{\hat{A}}}_{ab}], {A^{a}}_{b} ≐ γ^{ac} {\underline{A}}_{cb} .

(102)

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 $W_{z} [{\hat{A}}^{a}_{b}]$ because the Weyl symbols do not transform as tensors,

W_{z} [{\hat{A}}^{a}_{b}] \neq γ^{ac} W_{z} [{\underline{\hat{A}}}_{cb}],

(103)

except for those which are independent of p. Symbols $W_{z} [{\hat{A}}^{a}_{b}]$ will not be used below, and so no special notation is introduced for them.

Note that by Eq. (98), the following equality holds for any $Φ, Ψ$ ⁠, and $\hat{A}$ on $H_{m}^{n}$ ⁠:

\begin{matrix} ⟨ Ψ | \hat{A} Φ ⟩ m = \int d^{n} x \sqrt{g_{⋄}} γ_{ab} Ψ^{a *} ({\hat{A}}^{b}_{c} Φ^{c}) \\ = \int d^{n} x \sqrt{g_{⋄}} Ψ^{a *} ({\underline{\hat{A}}}_{ab} Φ^{b}) \\ = ⟨ Ψ^{a} | {\underline{\hat{A}}}_{ab} Φ^{b} ⟩, \end{matrix}

(104)

where we invoked the definition of the inner product on $H_{1}^{n}$ [Eq. (10)]. By the definition of the adjoint operator ${\hat{A}}^{†}$ ⁠,

\begin{matrix} ⟨ Ψ | \hat{A} Φ ⟩ m = {⟨ {\hat{A}}^{†} Ψ | Φ ⟩}_{m} \\ = \int d^{n} x \sqrt{g_{⋄}} γ_{cb} {[{({\hat{A}}^{†})}^{c}_{a} Ψ^{a}]}^{*} Φ^{b} \\ = \int d^{n} x \sqrt{g_{⋄}} {[{(\underline{{\hat{A}}^{†}})}_{ba} Ψ^{a}]}^{*} Φ^{b} . \end{matrix}

(105)

In order to express ${\hat{A}}^{†}$ through $\hat{A}$ ⁠, let us represent ${\underline{\hat{A}}}_{ab}$ in terms of the corresponding Weyl image ${\underline{A}}_{ab}$ and the Wigner operator ${\hat{Δ}}_{z}$ on $H_{1}^{n}$ ⁠, with $z \equiv (x, p)$ (Sec. IV A 3). Then, Eq. (105) leads to

\begin{matrix} ⟨ Ψ | \hat{A} Φ ⟩ m = \int \frac{d^{2 n} z'}{{(2 π)}^{n}} {\underline{A}}_{ab} (z') ⟨ Ψ^{a} | {\hat{Δ}}_{z'} Φ^{b} ⟩ \\ = \int \frac{d^{2 n} z'}{{(2 π)}^{n}} ⟨ {\underline{A}}_{ab}^{*} (z') {\hat{Δ}}_{z'} Ψ^{a} | Φ^{b} ⟩, \end{matrix}

(106)

where we used that ${\hat{Δ}}_{z'}$ is self-adjoint on $H_{1}^{n}$ ⁠. By comparing this with Eq. (105), one finds that, on one hand,

{(\underline{{\hat{A}}^{†}})}_{ba} = \int \frac{d^{2 n} z'}{{(2 π)}^{n}} {\underline{A}}_{ab}^{*} (z') {\hat{Δ}}_{z'} .

(107)

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

{(\underline{{\hat{A}}^{†}})}_{ba} = \int \frac{d^{2 n} z'}{{(2 π)}^{n}} W_{z'} [{(\underline{{\hat{A}}^{†}})}_{ba}] {\hat{Δ}}_{z'} .

(108)

Hence, the Weyl image of ${\hat{A}}^{†}$ is simply the conjugate transpose of the Weyl image of $\hat{A}$ in the sense that

W_{z} [{(\underline{{\hat{A}}^{†}})}_{ab}] = {\underline{A}}_{ba}^{*} (z) .

(109)

In other words, for operators with both indices lowered, the operations $W_{z}$ and $^{†}$ commute.

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

{(A_{H})}^{a}_{b} = \frac{1}{2} γ^{ac} ({\underline{A}}_{cb} + {\underline{A}}_{bc}^{*}),

(110)

{(A_{A})}^{a}_{b} = \frac{1}{2 i} γ^{ac} ({\underline{A}}_{cb} - {\underline{A}}_{bc}^{*}) .

(111)

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

{({\underline{A}}_{H})}_{ab} = \frac{1}{2} ({\underline{A}}_{ab} + {\underline{A}}_{ba}^{*}),

(112)

{({\underline{A}}_{A})}_{ab} = \frac{1}{2 i} ({\underline{A}}_{ab} - {\underline{A}}_{ba}^{*}) .

(113)

C. Weyl expansion of the dispersion operator

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

{\underline{\hat{D}}}_{ab} Ψ^{b} = 0, a = 1, 2, \dots m,

(114)

where ${\underline{\hat{D}}}_{ab} ≐ γ_{ac} {\hat{D}}^{c}_{b}$ ⁠. Equivalently, this can be written as the envelope equation,

{\underline{\hat{D}}}_{ab} ψ^{b} = 0, {\underline{\hat{D}}}_{ab} ≐ e^{- i θ} {\underline{\hat{D}}}_{ab} e^{i θ} .

(115)

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

{\hat{\underline{D}}}_{ab} \approx {\underline{D}}_{ab} + \frac{1}{2} [{\hat{p}}_{μ} {({\underline{V}}^{μ})}_{ab} + {({\underline{V}}^{μ})}_{ab} {\hat{p}}_{μ}] + \frac{1}{2} {\hat{p}}_{μ} {({\underline{Θ}}^{μ ν})}_{ab} {\hat{p}}_{ν} .

(116)

Here, the coefficients are x-dependent matrices; specifically, ${\underline{D}}_{ab} (x) ≐ {\underline{D}}_{ab} (x, k (x))$ ⁠, and

{[{\underline{V}}^{μ} (x)]}_{ab} ≐ \frac{\partial {\underline{D}}_{ab} (x, k (x))}{\partial k_{μ}},

(117)

{[{\underline{Θ}}^{μ ν} (x)]}_{ab} ≐ \frac{\partial^{2} {\underline{D}}_{ab} (x, k (x))}{\partial k_{μ} \partial k_{ν}} .

(118)

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

{\hat{D}}^{a}_{b} ψ^{b} = 0, {\hat{D}}^{a}_{b} ≐ γ^{ac} {\hat{\underline{D}}}_{cb} .

(119)

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

Unlike in scalar waves (Sec. IV C), ${\underline{D}}_{ab} (x, 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

{\underline{D}}_{ab} (x, k) = {\underline{D}}_{ab} (x, k) + {\underline{C}}_{ab} (x, k),

(120)

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

{[{\underline{V}}^{μ} (x)]}_{ab} ≐ \frac{\partial {\underline{D}}_{ab} (x, k (x))}{\partial k_{μ}},

(121)

{[Θ^{μ ν} (x)]}_{ab} ≐ \frac{\partial^{2} {\underline{D}}_{ab} (x, k (x))}{\partial k_{μ} \partial k_{ν}},

(122)

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:

\hat{D} \approx D_{H} + \hat{D},

(123)

\hat{D} = C_{H} + i D_{A} + {\hat{D}}_{H 1} + {\hat{D}}_{H 2} .

(124)

Here, $D \equiv D (x)$ and $D = γ^{- 1} \underline{D}$ ⁠, and so ${D^{a}}_{b} (x) = γ^{ac} {\underline{D}}_{cb} (x)$ ⁠; similar conventions are assumed for C and D. Accordingly, the matrices $D_{H}$ and D_A are self-adjoint and so are the operators ${\hat{D}}_{H 1}$ and ${\hat{D}}_{H 2}$ ⁠, which are given by

{\hat{D}}_{H 1} ≐ \frac{1}{2} γ^{- 1} ({\hat{p}}_{μ} {\underline{V}}_{H}^{μ} + {\underline{V}}_{H}^{μ} {\hat{p}}_{μ}),

(125)

{\hat{D}}_{H 2} ≐ \frac{1}{2} γ^{- 1} {\hat{p}}_{μ} {\underline{Θ}}_{H}^{μ ν} {\hat{p}}_{ν} .

(126)

Using $[γ^{- 1}, {\hat{p}}_{μ}] = i {(γ^{- 1})}_{, μ}$ and ${(γ^{- 1} γ)}_{, μ} = 0$ to obtain

{(γ^{- 1})}_{, μ} γ = - γ^{- 1} γ_{, μ} \equiv - ℓ_{μ},

(127)

one can also express these as

{\hat{D}}_{H 1} = \frac{1}{2} [{\hat{p}}_{μ} V_{H}^{μ} (x) + V_{H}^{μ} (x) {\hat{p}}_{μ} - i ℓ_{μ} V_{H}^{μ}],

(128)

{\hat{D}}_{H 2} \approx \frac{1}{2} {\hat{p}}_{μ} Θ_{H}^{μ ν} (x) {\hat{p}}_{ν} .

(129)

(The term $\propto ℓ_{μ}$ must be kept in ${\hat{D}}_{H 1}$ ⁠, but a similar term in ${\hat{D}}_{H 2}$ can be neglected.) Also,

V_{H}^{μ} (x) ≐ γ^{- 1} {\underline{V}}_{H}^{μ} (x), Θ_{H}^{μ ν} (x) ≐ γ^{- 1} {\underline{Θ}}_{H}^{μ ν} (x),

(130)

and since γ is k-independent, this leads to

{[V_{H}^{μ} (x)]}^{a}_{b} ≐ \frac{\partial {(D_{H})}^{a}_{b} (x, k (x))}{\partial k_{μ}},

(131)

{[Θ_{H}^{μ ν} (x)]}^{a}_{b} ≐ \frac{\partial^{2} {(D_{H})}^{a}_{b} (x, k (x))}{\partial k_{μ} \partial k_{ν}} .

(132)

If ${\hat{D}}_{H 2}$ is neglected, the envelope equation (119) becomes a Dirac-type equation similar to those considered in the context of XGO^24–28 and also, for instance, in Refs. 59 and 60. We shall also discuss this limit in Sec. V E.

D. Active and passive modes

1. Basis from the eigenvectors of $D_{H}$

Let us assume the GO ordering as in Eq. (63), except that $D_{H}$ and $D_{A}$ are now matrices. Then, $\hat{D} = O (ϵ)$ ⁠, and so the envelope equation [Eq. (5)] acquires the form

D_{H} (x) ψ = O (ϵ) .

(133)

Hence, it is convenient to express ψ through the eigenvectors of $D_{H} (x)$ ⁠. Let us denote these eigenvectors as $η_{s} (x)$ and the corresponding eigenvalues as $Λ_{s} (x)$ ⁠, and so

D_{H} η_{s} = Λ_{s} η_{s} .

(134)

[Note that $η_{s} (x) = η_{s} (x, k (x))$ ⁠, where $η_{s} (x, p)$ are the eigenvectors of $D_{H} (x, p)$ ⁠. Likewise, $Λ_{s} (x) = Λ_{s} (x, k (x))$ ⁠, where $Λ_{s} (x, p)$ are the eigenvalues of $D_{H} (x, p)$ ⁠.] Since $D_{H}$ is self-adjoint, it has m eigenvectors η_s which form a complete basis ${η_{s}}$ ⁠. Let us also introduce the corresponding dual basis ${η^{s}}$ as usual, i.e., such that

η^{s} \cdot η_{s'} = δ_{s^{'}}^{s} .

(135)

Then, we can represent ψ as

ψ = η_{s} a^{s}, a^{s} = η^{s} \cdot ψ \equiv {(η^{s})}_{a}^{*} ψ^{a},

(136)

where a^s serves as the components of ψ in the basis ${η_{s}}$ ⁠. [Remember that the dot product (101) includes conjugation.] Since $D_{H}$ is self-adjoint, it is always possible to make the basis ${η_{s}}$ orthonormal, and this choice is assumed below. Then, $η_{s} \cdot η_{s'} = δ_{ss'}$ ⁠; thus, $η_{s} = η^{s}$ ⁠, i.e.,

{(η^{s})}_{a} = {(η_{s})}_{a} = γ_{ab} {(η_{s})}^{b} .

(137)

Also, for any two fields $ψ = η_{s} a^{s}$ and $ϕ = η_{s} b^{s}$ ⁠, the inner product (98) can be written as follows:

{⟨ ψ | ϕ ⟩}_{m} = \int d^{n} x \sqrt{g_{⋄} (x)} a_{s}^{*} (x) b^{s} (x),

(138)

where $a_{s} ≐ δ_{ss'} a^{s^{'}}$ ⁠. The matrix δ with elements $δ_{ss'}$ serves as a Euclidean metric for manipulating the mode indices s and $s'$ ⁠. 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 γ [Eq. (99)].

2. Amplitude vectors

The physical meaning of the expansion coefficients a^s is understood as follows. Consider multiplying Eq. (133) by η^s from the left. This gives $Λ_{s} a^{s} = O (ϵ)$ ⁠, where no summation over s is assumed. This shows that, for a given s, there are two possibilities: either a^s is small or Λ_s is small. In the first case, the polarization η_s does not correspond to a propagating wave mode per se; the small nonzero projection of ψ on η_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 $k (x)$ approximately satisfies the local dispersion relation $Λ_{s} (x, k (x)) \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-mode 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, \dots, N$ and the remaining $\bar{N} = m - N$ modes with $s = (N + 1), (N + 2), \dots, m$ are passive. Let us also adopt the notation

{\bar{a}}^{s} = a^{s + N}, {\bar{η}}_{s} = η_{s + N}, {\bar{Λ}}_{s} = Λ_{s + N}

(139)

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

a ≐ (\begin{matrix} a^{1} \\ ⋮ \\ a^{N} \end{matrix}), \bar{a} ≐ (\begin{matrix} {\bar{a}}^{1} \\ ⋮ \\ {\bar{a}}^{\bar{N}} \end{matrix})

(140)

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

a^{†} = (a_{1}^{*}, a_{2}^{*}, \dots, a_{N}^{*}), {\bar{a}}^{†} = ({\bar{a}}_{1}^{*}, {\bar{a}}_{2}^{*}, \dots, {\bar{a}}_{\bar{N}}^{*}),

(141)

a_{s} ≐ δ_{ss'} a^{s^{'}}, {\bar{a}}_{s} ≐ δ_{ss'} {\bar{a}}^{s^{'}} .

(142)

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

3. Polarization matrices

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

Ξ = (η_{1}, η_{2}, \dots, η_{N}), \bar{Ξ} = ({\bar{η}}_{1}, {\bar{η}}_{2}, \dots, {\bar{η}}_{\bar{N}}) .

(143)

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

{Ξ^{a}}_{s} = {(η_{s})}^{a}, {\bar{Ξ}}^{a}_{s} = {({\bar{η}}_{s})}^{a} .

(144)

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

ψ = Ξ a + \bar{Ξ} \bar{a} .

(145)

Also, consider the auxiliary polarization matrices

Ξ^{+} ≐ (\begin{matrix} η^{1 *} \\ ⋮ \\ η^{N *} \end{matrix}), {\bar{Ξ}}^{+} ≐ (\begin{matrix} {\bar{η}}^{1 *} \\ ⋮ \\ {\bar{η}}^{\bar{N} *} \end{matrix}),

(146)

which have $η^{s *}$ as their rows,

{(Ξ^{+})}^{s}_{a} = {(η^{s})}_{a}^{*}, {({\bar{Ξ}}^{+})}^{s}_{a} = {({\bar{η}}^{s})}_{a}^{*} .

(147)

Since

{(ψ^{†})}_{a} = {(η_{s})}_{a}^{*} a^{s *} + {({\bar{η}}_{s})}_{a}^{*} {\bar{a}}^{s *} = {(η^{s})}_{a}^{*} a_{s}^{*} + {({\bar{η}}^{s})}_{a}^{*} {\bar{a}}_{s}^{*},

(148)

one obtains

ψ^{†} = a^{†} Ξ^{+} + {\bar{a}}^{†} {\bar{Ξ}}^{+} .

(149)

However, note that in general, $^{+}$ does not mean to represent the adjoint in the common sense [cf. Eq. (154)].

Next, notice that

{(Ξ^{+} Ξ)}^{s}_{s'} = {(η^{s})}_{a}^{*} {(η_{s^{'}})}^{a} = η^{s} \cdot η_{s'} = δ_{s^{'}}^{s},

(150)

where the indices s and $s'$ span from 1 to N. Analogous formulas apply to $\bar{Ξ}$ ⁠. Then,

{({\bar{Ξ}}^{+} Ξ)}^{s}_{s'} = {({\bar{η}}^{s})}_{a}^{*} {(η_{s'})}^{a} = {\bar{η}}^{s} \cdot η_{s'} = 0,

(151)

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

Ξ^{+} Ξ = 1, {\bar{Ξ}}^{+} \bar{Ξ} = 1, {\bar{Ξ}}^{+} Ξ = 0, Ξ^{+} \bar{Ξ} = 0,

(152)

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

a = Ξ^{+} ψ, \bar{a} = {\bar{Ξ}}^{+} ψ .

(153)

Another property that we shall use later on is

{(Ξ^{+})}^{s}_{a} = δ^{{ss}^{'}} {(η_{s^{'}})}_{a}^{*} = δ^{{ss}^{'}} {(η_{s^{'}})}^{b *} γ_{ab} = {(δ^{- 1} Ξ_{γ}^{H})}^{s}_{a},

and similarly for ${\bar{Ξ}}^{+}$ ⁠. Here, $^{H}$ denotes the conjugate transpose (⁠ $Ξ^{H} ≐ Ξ^{T *}$ ⁠) and $δ^{- 1}$ is the Kronecker matrix with upper-index elements $δ^{{ss}^{'}}$ ⁠. Hence,

Ξ^{+} = δ^{- 1} Ξ^{H} γ, {\bar{Ξ}}^{+} = δ^{- 1} {\bar{Ξ}}^{H} γ .

(154)

The factor $δ^{- 1}$ can be omitted if one ignores the difference between upper and lower indices s and $s'$ that refer to the mode number. If γ is Euclidean, one can also ignore the difference between upper and lower coordinate indices; then, $Ξ^{+} = Ξ^{H}$ ⁠.

As a side remark, note that Ξ and $\bar{Ξ}$ are generally nonsquare, and so $Π ≐ Ξ Ξ^{+}$ and $\bar{Π} ≐ \bar{Ξ} {\bar{Ξ}}^{+}$ are not unit matrices but rather projectors. Indeed, due to Eqs. (152), one has $Π^{2} = Π$ and ${\bar{Π}}^{2} = \bar{Π}$ ⁠. Also, by applying Π and $\bar{Π}$ to Eq. (145), one obtains

Ξ a = Π ψ, \bar{Ξ} a = \bar{Π} ψ .

(155)

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

4. Equation for the active modes

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

D_{H} = η_{s} Λ_{s} η^{s} = Ξ Λ Ξ^{+} + \bar{Ξ} \bar{Λ} {\bar{Ξ}}^{+},

(156)

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

Λ ≐ diag {Λ_{1}, Λ_{2}, \dots, Λ_{N}} = O (ϵ),

(157)

\bar{Λ} ≐ diag {{\bar{Λ}}_{1}, {\bar{Λ}}_{2}, \dots, {\bar{Λ}}_{\bar{N}}} = O (1) .

(158)

These and Eq. (152) also lead to

D_{H} Ξ = Ξ Λ = O (ϵ), D_{H} \bar{Ξ} = \bar{Ξ} \bar{Λ} = O (1) .

(159)

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

\begin{matrix} 0 = D_{H} Ξ a + D_{H} \bar{Ξ} \bar{a} + \hat{D} Ξ a + \hat{D} \bar{Ξ} \bar{a} \\ = Ξ Λ a + \bar{Ξ} \bar{Λ} \bar{a} + \hat{D} Ξ a + \hat{D} \bar{Ξ} \bar{a}, \end{matrix}

(160)

where we used Eq. (159). Let us multiply this by $Ξ^{+}$ and, separately, by ${\bar{Ξ}}^{+}$ ⁠. Then, due to Eq. (152), one obtains

Λ a + Ξ^{+} \hat{D} Ξ a + Ξ^{+} \hat{D} \bar{Ξ} \bar{a} = 0,

(161)

\bar{Λ} \bar{a} + {\bar{Ξ}}^{+} \hat{D} Ξ a + {\bar{Ξ}}^{+} \hat{D} \bar{Ξ} \bar{a} = 0.

(162)

Let us also use Eq. (162) to express $\bar{a}$ through a and substitute the result into Eq. (161). Since both $\bar{a}$ and $\hat{D}$ are of order ϵ, it is sufficient to solve Eq. (162) for $\bar{a}$ approximately to the leading order,

\bar{a} \approx - {\bar{Λ}}^{- 1} {\bar{Ξ}}^{+} \hat{D} Ξ a .

(163)

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

(Λ + Ξ^{+} \hat{D} Ξ - Ξ^{+} \hat{D} \bar{Ξ} {\bar{Λ}}^{- 1} {\bar{Ξ}}^{+} \hat{D} Ξ) a = 0.

(164)

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

E. Leading-order approximation: Extended geometrical optics

Upon substituting Eqs. (124)–(126) into Eq. (164), we obtain the following equation to lowest (first) order in ϵ:

(Λ + i Γ + \hat{K}) a = 0,

(165)

where we introduced

Γ ≐ Ξ^{+} D_{A} Ξ,

(166)

\hat{K} ≐ Ξ^{+} ({\hat{D}}_{H 1} + C_{H}) Ξ .

(167)

As shown in Appendix B 1,

\hat{K} a = - i V^{μ} a_{, μ} - \frac{i}{2} {V^{μ}}_{; μ} a - Ua,

(168)

V^{μ} ≐ Ξ^{+} V_{H}^{μ} Ξ \approx Λ^{| μ},

(169)

U = δ^{- 1} {(Ξ_{, μ}^{H} {\underline{V}}_{H}^{μ} Ξ)}_{A} - Ξ^{+} C_{H} Ξ,

(170)

where Ξ is considered as a function of x. [As a reminder, $Ξ^{H} ≐ Ξ^{T *}$ is the conjugate transpose of Ξ; ${\underline{V}}_{H}^{μ}$ is given by Eq. (121); and $δ^{- 1}$ is a unit matrix which only raises the mode index.] Alternatively, Ξ can be considered as a function of $(x, k)$ ⁠. Then, as shown in Appendix B 2,

\begin{matrix} U \approx δ^{- 1} {(Λ_{| μ} Ξ^{H} γ Ξ^{| μ} - Λ^{| μ} Ξ^{H} γ Ξ_{| μ} + Ξ^{H | μ} {\underline{D}}_{H} Ξ_{| μ})}_{A} \\ - Ξ^{+} C_{H} Ξ, \end{matrix}

(171)

where the partial derivatives $_{| μ}$ and $^{| μ}$ are defined in Sec. IV C. 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.^24–28 These effects are discussed in further detail in Sec. V G.

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

Λ a - i V^{μ} a_{, μ} - \frac{i}{2} {V^{μ}}_{; μ} a = (U - i Γ) a .

(172)

This generalizes the XGO equation derived in Refs. 24–29 to dissipative waves and curved coordinates. Since the vector a belongs to the space where the metric is Euclidean by definition (Sec. V D 2), the elements of this vector are covariant, i.e., invariant with respect to coordinate transformations in the physical space. The same applies to the eigenvalues comprising Λ; thus, the whole left-hand side of Eq. (172) 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 $Γ$ are covariant.

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

{J^{μ}}_{; μ} = - 2 a^{†} Γ a, J^{μ} ≐ - a^{†} V^{μ} a .

(173)

Using Eqs. (100), (130), (145), (149), (155), and (169), one obtains

\begin{matrix} - J^{μ} = a^{†} Ξ^{+} V_{H}^{μ} Ξ a = {(Π ψ)}^{†} V_{H}^{μ} (Π ψ) \approx ψ^{†} V_{H}^{μ} ψ = ψ^{*} {\underline{V}}_{H}^{μ} ψ, \end{matrix}

(174)

and so Eq. (173) can be viewed as a generalization of Eq. (67). At zero $Γ$ ⁠, the dynamics is conservative, and Eq. (173) reflects conservation of the wave action or quanta.⁵⁴ Accordingly, $Γ$ determines the dissipation rate, $J^{μ}$ can be identified (up to a constant factor) as the action flux density summed over all active modes (for example, see Sec. VI B), and the elements of the diagonal matrix $Λ^{| μ}$ 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 F.

For scalar waves studied in Sec. IV, we had $Λ = D_{H}$ and we defined $k (x)$ such that this term be zero [Eq. (64)]. Now, Λ is a diagonal matrix with $N \geq 1$ nonzero elements, and so it cannot be zeroed entirely by imposing just one scalar constraint on $k (x)$ ⁠. Hence, there can be more than one natural way to define $k (x)$ ⁠. One esthetically pleasing and convenient²⁹ option is to require that Λ or $Λ - U$ be traceless. This amounts to choosing the reference-ray Hamiltonian in Sec. III as $H = tr Λ / N$ or $H = tr (Λ - U) / N$ ⁠, respectively. [Arbitrary constant factors can be introduced instead of N, for that only redefines τ in Eq. (7).] Another option is to adopt $Λ_{s} = 0$ for some single $s \leq N$ ⁠, since all active modes have close wave vectors anyway; then, $H = Λ_{s}$ ⁠. As mentioned in Sec. III, all such choices of $k (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 = Λ - U$ is preferable, as discussed in Sec. V G 2.

F. Quasioptical model

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 $V_{avr}$ ⁠, which can be defined via $N V_{avr}^{μ} ≐ tr Λ^{| μ}$ ⁠. Then,

V^{μ} = V_{avr}^{μ} 1 + Δ V^{μ}, Δ V^{μ} ≐ V^{μ} - V_{avr}^{μ} 1,

(175)

where $Δ V^{μ}$ are matrices with eigenvalues much smaller than $V_{avr}$ ⁠. Let us assume a ray-based coordinate system aligned with $V_{avr}$ with the inner product on the transverse space $M_{⊥}^{n - 1}$ defined similar to Eq. (76) [cf. also Eq. (138)], namely,

{⟨ a | b ⟩}_{N}^{⊥} ≐ \int d^{n - 1} ϱ \sqrt{h_{⋄} (ζ, ϱ)} a_{s}^{*} (ζ, ϱ) b^{s} (ζ, ϱ) .

(176)

Let us also adopt the quasioptical ordering as we did for scalar waves in Sec. IV E. In particular, we allow $Δ V^{μ} = O (ϵ_{⊥})$ ⁠. Then, Eq. (164) becomes

(Λ + i Γ + \hat{K} + \hat{G}) a = 0.

(177)

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

\hat{K} a = - iV a_{, ζ} - \frac{i}{2 \sqrt{g_{⋄}}} {(\sqrt{g_{⋄}} V)}_{, ζ} a + Δ \hat{K} a,

(178)

Δ \hat{K} a ≐ - Ua - i Δ V^{σ} a_{, σ} - \frac{i}{2 \sqrt{h_{⋄}}} {(\sqrt{h_{⋄}} Δ V^{σ})}_{, σ} a,

(179)

Δ V^{σ} = Ξ^{+} {D_{H}}^{| σ} Ξ,

(180)

where $V ≐ V_{avr}^{0}$ (we used the fact that $V_{avr}^{σ} = 0$ by definition of the ray-based coordinates), $U = O (ϵ_{| |})$ is given by Eq. (171), and $Γ$ is given by Eq. (166). Finally,

\hat{G} ≐ Ξ^{+} {\hat{D}}_{H 2} Ξ - Ξ^{+} {\hat{D}}_{H 1} \bar{Ξ} {\bar{Λ}}^{- 1} {\bar{Ξ}}^{+} {\hat{D}}_{H 1} Ξ .

(181)

As shown in Appendix B 3, $\hat{G}$ can also be simplified as

\hat{G} a \approx - \frac{1}{2 h_{⋄}^{1 / 4}} {[Λ^{| σ \tilde{σ}} {(h_{⋄}^{1 / 4} a)}_{, \tilde{σ}}]}_{, σ},

(182)

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

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

Λ a + i Γ a - iV a_{, ζ} - \frac{i}{2 \sqrt{g_{⋄}}} {(\sqrt{g_{⋄}} V)}_{, ζ} a + \hat{G} a + Δ \hat{K} a = 0.

(183)

Equation (183) is the main result of our paper. Like Eq. (172), this equation is covariant within the accuracy of our theory. As a reminder, a is generally a vector (⁠ $dim a = N \geq 1$ ⁠), whose elements are the scalar amplitudes of active modes (Sec. V D 2); Λ is the diagonal matrix formed by the eigenvalues of $D_{H}$ ⁠; $Γ$ is given by Eq. (166); $V = N^{- 1} tr Λ^{| ζ}$ is a scalar; $^{| ζ} ≐ \partial / \partial k_{0}$ ⁠; $\hat{G}$ is given by Eq. (182); and $Δ \hat{K}$ is given by Eq. (179). There, U is given by Eq. (171) and $Δ V$ is given by Eq. (180). The derivatives of $h_{⋄}$ could be neglected within the assumed accuracy, but a purist might want to retain them in order to keep $Δ \hat{K}$ and $\hat{G}$ precisely self-adjoint under the inner product (176).

Like in Sec. IV E 2, we require the α parameter of the ray-based coordinates to be defined as $ϱ$ -independent. Then, Eq. (183) has the following corollary:

\frac{1}{α} \frac{d}{d ζ} {⟨ a | α V | a ⟩}_{N}^{⊥} = 2 {⟨ a | Γ | a ⟩}_{N}^{⊥},

(184)

which is similar to Eq. (84). For the case of zero $Γ$ ⁠, Eq. (184) predicts conservation of the wave-action flux through the beam cross section,

{⟨ ψ | α V | ψ ⟩}_{N}^{⊥} = const .

(185)

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

2. Simplified equations

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

\begin{matrix} Λ a + i Γ a - iV a_{, ζ} - \frac{i}{2} V_{, ζ} a - \frac{1}{2} {(Λ^{| σ \tilde{σ}} a_{, \tilde{σ}})}_{, σ} \\ - Ua - i Δ V^{σ} a_{, σ} - \frac{i}{2} Δ {V^{σ}}_{, σ} a = 0. \end{matrix}

(186)

Using the variable transformation $ϕ ≐ \sqrt{V} a$ ⁠, this equation can be further simplified as

i ϕ_{, ζ} = \hat{χ} ϕ + i ϒ ϕ .

(187)

Here, $\hat{χ}$ is an operator self-adjoint under the inner product (176) and given by

\hat{χ} ϕ \approx Q ϕ - \frac{1}{2} {(Σ^{σ \tilde{σ}} ϕ_{, \tilde{σ}})}_{, σ} - \frac{Δ V^{σ}}{V} i ϕ_{, σ} - \frac{i}{2} {(\frac{Δ V^{σ}}{V})}_{, σ} ϕ .

Also, Q, Σ, and $ϒ$ are the self-adjoint matrices given by

Q ≐ (Λ - U) / V, Σ^{σ \tilde{σ}} ≐ Λ^{| σ \tilde{σ}} / V, ϒ ≐ Γ / V .

(188)

Accordingly, Eq. (184) becomes

\frac{d}{d ζ} {⟨ ϕ | ϕ ⟩}_{N}^{⊥} = 2 {⟨ ϕ | ϒ | ϕ ⟩}_{N}^{⊥},

(189)

and so ${⟨ ϕ | ϕ ⟩}_{N}^{⊥}$ is conserved if $ϒ$ is zero.

G. Summary and discussion

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

1. Scalar beams

In the case of a scalar beam (Sec. IV), one has $Ξ = 1$ ⁠, and so

Λ = D_{H}, Γ = D_{A}, U = 0, Δ V^{σ} = 0.

(190)

Then, the equations from Sec. IV E are reproduced. In particular, Λ is made zero by the choice of k [Eq. (64)]. This implies that the ray Hamiltonian is $H = Λ$ ⁠, which leads to the following ray equations:

\frac{d x^{α}}{d τ} = \frac{\partial Λ}{\partial k_{α}}, \frac{d k_{α}}{d τ} = - \frac{\partial Λ}{\partial x^{α}} .

(191)

2. Single-mode vector beams

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

Λ = η^{†} D_{H} η, Γ = η^{†} D_{A} η,

(192)

and so they are scalars. Also, $Δ V^{σ} = 0$ ⁠, but U is generally a nonzero scalar function. There are two natural ways to define k in this case. One is to require that $Λ = 0$ ⁠, i.e., to adopt $H = Λ$ ⁠. 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 multimode 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 k by requiring $Λ - U = 0$ ⁠, i.e., by adopting $H = Λ - 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,

\frac{d x^{α}}{d τ} = \frac{\partial (Λ - U)}{\partial k_{α}}, \frac{d k_{α}}{d τ} = - \frac{\partial (Λ - U)}{\partial x^{α}} .

(193)

The effect of U on the ray trajectories is known as the (spin) Hall effect of light in optics.^{22,23,61–64} 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 Λ is a scalar, one can rewrite Eq. (171) for $U (x, k)$ as follows:

\begin{matrix} U = U_{0} - Λ^{| μ} Im (η^{†} η_{| μ}) + Λ_{| μ} Im (η^{†} η^{| μ}) \\ \approx U_{0} - {\dot{x}}^{μ} A_{μ}^{(x)} - {\dot{k}}_{μ} A_{(k)}^{μ} . \end{matrix}

(194)

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

U_{0} ≐ Im (η^{H | μ} {\underline{D}}_{H} η_{| μ}) - η^{H} {\underline{C}}_{H} η,

(195)

A_{μ}^{(x)} ≐ Im (η^{†} η_{| μ}), A_{(k)}^{μ} ≐ Im (η^{†} η^{| μ}) .

(196)

Then, the phase-space Lagrangian of a ray, $L [x, k] = k_{μ} {\dot{x}}^{μ} - H$ ⁠, has a noncanonical structure, namely,

L = [k_{μ} + A_{μ}^{(x)}] {\dot{x}}^{μ} + A_{(k)}^{μ} {\dot{k}}_{μ} - (Λ - U_{0}) .

(197)

The equations originating from a special case of this noncanonical phase-space Lagrangian were used, for example, in Ref. 62 to study the Hall effect for light propagating in a nonbirefringent material. A comparison of the resulting noncanonical ray equations with the canonical ray equations (193) can be found in Ref. 27. Also notably, the terms $A_{μ}^{(x)}$ and $A_{(k)}^{μ}$ are known as the Berry connection and are linked to the Berry phase or the geometrical phase of light.^27,44,62

3. Multimode 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 $Γ, Δ V^{σ}$ ⁠, and U are generally nondiagonal, and 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. (187), which then becomes

i ϕ_{, ζ} = Q ϕ,

(198)

with Q being a self-adjoint matrix. The modes decouple when Q is close to diagonal. More generally, the dynamics governed by Eq. (198) 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 ζ⁶⁵). Alternatively, it can be mapped to a precession equation for a real $(N^{2} - 1)$ -dimensional “spin” vector.²⁶ 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 Ref. 29.

As a reminder, the model presented here (Sec. V F) 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_{, ζ} ≪ a_{, σ}$ does not hold. The quasioptical description is inapplicable to such beams; however, the first-order XGO model described in Sec. V E can be used.

VI. ELECTROMAGNETIC WAVES

A. Covariant formulation

Here, we shall explain how the above theory applies to EM waves. We start with Maxwell's equation⁶⁸

\frac{1}{\sqrt{g_{⋄}}} \partial_{α} (\sqrt{g_{⋄}} F^{α β}) = - \frac{4 π}{c} J^{β},

(199)

where $F^{α β}$ is the EM tensor, namely,

F^{α β} = g^{α μ} g^{β ν} F_{μ ν}, F_{α β} ≐ \partial_{α} A_{β} - \partial_{β} A_{α},

(200)

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

{\hat{D}}_{vac}^{(A)} A = - \frac{4 π}{c} J .

(201)

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

\begin{matrix} {[{\hat{D}}_{vac}^{(A)}]}^{α}_{β} A^{β} = \frac{1}{\sqrt{g_{⋄}}} \partial_{λ} {\sqrt{g_{⋄}} g^{α μ} g^{λ ν} \\ \times [\partial_{ν} (g_{μ β} A^{β}) - \partial_{μ} (g_{ν β} A^{β})]}, \end{matrix}

(202)

which is self-adjoint under the inner product (98). By assuming $J = c \hat{ϑ} A$ ⁠, where $\hat{ϑ}$ is some linear operator, one can cast the equation for A in the form (1), namely,

[{\hat{D}}_{vac}^{(A)} + 4 π \hat{ϑ}] A = 0.

(203)

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 $γ = g$ ⁠. The corresponding Weyl symbols are calculated as follows. It can be shown straightforwardly that ${[{\hat{\underline{D}}}_{vac}^{(A)}]}_{α β} ≐ g_{α γ} {[{\hat{D}}_{vac}^{(A)}]}^{γ}_{β}$ is approximately given by

\begin{matrix} {[{\hat{\underline{D}}}_{vac}^{(A)}]}_{α β} \approx {\hat{p}}_{α} {\hat{p}}_{β} - {\hat{p}}_{μ} g^{μ ν} g_{α β} {\hat{p}}_{ν} \\ + i (q_{α} {\hat{p}}_{β} - q_{β} {\hat{p}}_{α}) + i [{(ℓ_{β})}^{μ}_{α} - {(ℓ_{α})}^{μ}_{β}] {\hat{p}}_{μ} . \end{matrix}

(204)

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 $ℓ_{μ} ≐ g^{- 1} g_{, μ}$ ⁠, which is in agreement with Eq. (127) because $g = γ$ ⁠. Finally, in the form introduced in Sec. V C, the above result is

\begin{matrix} {[{\underline{D}}^{(A)}]}_{α β} \approx p_{α} p_{β} - g_{α β} g^{μ ν} p_{μ} p_{ν} + 4 π {({\underline{ϑ}}_{0})}_{α β}, \\ C_{α β}^{(A)} \approx i (q_{α} k_{β} - q_{β} k_{γ}) + i [{(ℓ_{β})}^{μ}_{α} - {(ℓ_{α})}^{μ}_{β}] k_{μ} + 4 π {(Δ \underline{ϑ})}_{α β} \end{matrix}

where ${({\underline{ϑ}}_{0})}_{α β}$ is the GO limit of $\underline{ϑ}, {(Δ \underline{ϑ})}_{α β}$ is the remaining part of $\underline{ϑ}$ ⁠, and $\underline{ϑ}$ itself is the symbol of $\underline{\hat{ϑ}}$ ⁠.

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. 23, 63, and 64.)

B. Noncovariant formulation

As a special case, suppose a metric of the form

g = (\begin{matrix} - 1 & 0 \\ 0 & h \end{matrix}),

(205)

with time-independent spatial metric 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. (203) and replace the vector potential with the electric field,

E^{a} ≐ F^{0 a} = - i {\hat{p}}_{0} A^{a} .

(206)

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

{\hat{D}}^{(E)} E = 0,

(207)

{\hat{D}}^{(E)} ≐ {({\hat{p}}_{0})}^{- 1} [{\hat{D}}_{vac}^{(A)} + 4 π \hat{ϑ}] {({\hat{p}}_{0})}^{- 1},

(208)

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

Let us introduce the conductivity operator $\hat{σ}$ via $J = \hat{σ} E$ and notice that $4 π \hat{σ} = ic {\hat{p}}_{0} \hat{χ}$ by the definition of the susceptibility operator $\hat{χ}$ ⁠. Then, $4 π \hat{ϑ} = {\hat{p}}_{0} \hat{χ} {\hat{p}}_{0}$ ⁠, and so ${\hat{D}}^{(E)}$ can be expressed as follows:

{\hat{D}}^{(E)} ≐ {({\hat{p}}_{0})}^{- 1} {\hat{D}}_{vac}^{(A)} {({\hat{p}}_{0})}^{- 1} + \hat{χ},

(209)

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

{\underline{D}}_{ab}^{(E)} = {(p_{0})}^{- 2} [p_{a} p_{b} - p_{c} p_{d} h^{cd} h_{ab} + {\underline{C}}_{ab}^{(A)}] + {\underline{ε}}_{ab},

(210)

where ${\underline{ε}}_{ab} ≐ h_{ab} + {({\underline{χ}}_{0})}_{ab}$ serves as the dielectric tensor and ${\underline{χ}}_{0}$ is the GO limit of $\underline{χ}$ ⁠.

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

{\underline{D}}_{ab}^{(E)} = {\underline{D}}_{ab}^{(E)} + {\underline{C}}_{ab}^{(E)},

(211)

{\underline{D}}_{ab}^{(E)} = \frac{1}{16 π p_{0}^{2}} (p_{a} p_{b} - p_{c} p_{d} h^{cd} h_{ab}) + \frac{{\underline{ε}}_{ab}}{16 π},

(212)

{\underline{C}}_{ab}^{(E)} ≐ \frac{1}{16 π p_{0}^{2}} {\underline{C}}_{ab}^{(A)} .

(213)

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

q_{a} = \frac{1}{4} \partial_{a} ln | \det h |, ℓ_{a} = h^{- 1} \partial_{a} h,

(214)

and so they both can be of order $ϵ_{| |}$ ⁠; hence, $C = O (ϵ_{| |})$ ⁠. Note that such C is generally non-negligible even for near-Euclidean metrics such as those discussed in Sec. IV E 3.

From Eq. (174), we obtain

J^{α} = - E^{a *} E^{b} \frac{\partial}{\partial k_{α}} {[{\underline{D}}_{H}^{(E)} (t, x, ω, k)]}_{ab},

(215)

which acts as the true action-flux density.⁵⁴ As a reminder, the action density $I ≐ J^{0}$ can be written as

I \approx \frac{E^{a *} E^{b}}{16 π ω^{2}} \frac{\partial}{\partial ω} {ω^{2} {[{\underline{ε}}_{H} (t, x, ω, k)]}_{ab}},

(216)

where we used ${[{\underline{D}}_{H}^{(E)}]}_{ab} E^{b} \approx 0$ ⁠. Using the latter and Faraday's law $B \approx c k \times E / ω$ for the magnetic field B, one can also rewrite $I$ in another common form,⁴

I \approx \frac{E^{a *} E^{b}}{16 π ω} \frac{\partial}{\partial ω} {ω {[{\underline{ε}}_{H} (t, x, ω, k)]}_{ab}} + \frac{B_{a}^{*} B^{a}}{16 π ω} .

(217)

One can also show⁵⁴ that the spatial component of the action flux density (215) can be expressed as

J \approx \frac{S}{ω} - \frac{E^{a *} E^{b}}{16 π} \frac{\partial}{\partial k} {[{\underline{ε}}_{H} (t, x, ω, k)]}_{ab},

(218)

where S is the (time- or space-averaged) Poynting vector,

S = \frac{c}{8 π} Re (E \times B^{*}) .

(219)

The GO equation (173) can be written as

\frac{\partial I}{\partial t} + \frac{1}{\sqrt{h_{⋄}}} \frac{\partial}{\partial x^{a}} (\sqrt{h_{⋄}} J^{a}) = - \frac{E^{a *} E^{b}}{8 π} {[{\underline{ε}}_{A} (t, x, ω, k)]}_{ab} .

The corresponding density of the wave energy is $ω I$ ⁠, and the density of the energy flux is $ω J$ ⁠, as flows from the variational principle.⁵⁴

C. Stationary waves

For stationary waves with fixed frequency, one has ${\hat{p}}_{0} = - ω / c$ and ω can be treated as a constant parameter. Such waves can be studied on the three-dimensional configuration space, namely, the physical space x. In this case, $n = m = 3$ and $γ = h$ ⁠. If there is no spatial dispersion, then $\hat{χ} = χ (\hat{x})$ ⁠, where the dependence on ω is assumed but not emphasized. Assuming that the underlying space is flat (which is always the case in laboratory applications), the symbol ${χ^{a}}_{b}$ is then 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²⁷ and for waves in cold plasmas.^24,25,67)

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

{\underline{χ}}_{ab} = {\underline{χ}}_{ab}^{(0)} (x) + {\underline{χ}}_{ab}^{(k)} (x, p), {\underline{χ}}_{ab}^{(k)} ≪ {\underline{χ}}_{ab}^{(0)} .

(220)

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

In principle, our general method is also applicable to waves in media with strong spatial dispersion. However, $\hat{χ}$ 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 $\hat{χ}$ ⁠. A related discussion can also be found in Ref. 43.

D. Waves in flat space

Let us also discuss the special case when the configuration space is flat and Eq. (220) applies, as in typical laboratory applications. Although the ray-based metric 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 ${\tilde{h}}_{ab} = δ_{ab}$ and $\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, tildes are quantities in the ray-based coordinates, and nontilded 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 C). Then, since U and $Γ$ are frame-invariant (Sec. V E), one can calculate them in the laboratory coordinates using

\begin{matrix} U = \tilde{U} \approx {({\tilde{Λ}}_{| μ} {\tilde{Ξ}}^{H} {\tilde{Ξ}}^{| μ} - {\tilde{Λ}}^{| μ} {\tilde{Ξ}}^{H} {\tilde{Ξ}}_{| μ} + {\tilde{Ξ}}^{H | μ} {\tilde{D}}_{H} {\tilde{Ξ}}_{| μ})}_{A}, \end{matrix}

(221)

\begin{matrix} Γ = \tilde{Γ} = {\tilde{Ξ}}^{H} {\tilde{D}}_{A} \tilde{Ξ}, \end{matrix}

(222)

where both $δ^{- 1}$ and $\tilde{γ}$ 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.⁶⁸ In contrast, the O–X conversion caused by the magnetic shear in dilute plasmas²⁹ 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 G 2. (In special cases like the X–B transformation,⁴ when the group velocity becomes zero on a ray, while k remains large, one may want to replace the coordinate ζ with τ, where $d τ = d ζ / V$ ⁠, to remove the singularity in the amplitude equation.) In all these cases, our theory is expected to yield results which agree with full-wave modeling up to local errors of order $O (ϵ_{∥})$ ⁠, because this is the accuracy within which the passive-mode contribution is calculated (Sec. V D 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 $\hat{D}$ that governs the original wave field Ψ, we explicitly calculate the approximate operator $\hat{D}$ that governs the wave envelope ψ to the second order in the GO parameter ϵ. Then, we further simplify this envelope operator by assuming that the gradient of ψ transverse to the local group velocity is much larger than the corresponding parallel gradient. This leads to a parabolic differential equation for ψ (“quasioptical equation”) on the basis of the GO polarization vectors [Eq. (136)]. Our main results can be found in Sec. V F, which includes the general quasioptical equation (183) and its special case (187). Scalar and mode-converting vector beams are described on the same footing (Sec. V G). 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 radio frequency wave beams in magnetized plasma based on this theory.

SUPPLEMENTARY MATERIAL

See the supplementary material for an overview of the Weyl calculus on a curved configuration space. There, we present proofs of the statements summarized in Sec. IV A.

ACKNOWLEDGMENTS

The first author (I.Y.D.) 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. This work was supported by JSPS KAKENHI Grant No. JP17H03514. This 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 No. DE-NA-0003525. This paper describes objective technical results and analysis. Any subjective views or opinions that might be expressed in this paper do not necessarily represent the views of the U.S. Department of Energy or the United States Government.

APPENDIX A: SUMMARY OF THE SELECTED NOTATIONS

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

1. Basic symbols

• $≐$ denotes a definition.

• ^̂ denotes an operator.

• $^{†}$ denotes either a dual vector [Eqs. (100) and (141)] or an adjoint operator.

• $^{T}$ denotes the matrix transpose.

• $^{H} =^{T *}$ denotes the conjugate matrix transpose.

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

• _H and _A denote the Hermitian and anti-Hermitian parts of an operator [Eqs. (26)] or those of a matrix (Sec. V B). When applied to a scalar, _H and _A denote the real and imaginary parts, respectively.

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

• $W_{z}$ and $W^{- 1}$ denote the Wigner–Weyl transform and its inverse (Sec. IV A 3), respectively.

2. Index manipulation

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

X_{μ} = g_{μ ν} X^{ν}, μ, ν = 0, 1, \dots, (n - 1) .

(A1)

(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ⁿ (Sec. V A). On that space, an additional metric γ is introduced, and the indices are manipulated as follows:

Ψ_{a} = γ_{ab} Ψ^{b}, a, b = 1, 2, \dots, m .

(A2)

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

We also introduce the Euclidean metric δ on the space of amplitude vectors a (Sec. V D 2),

a_{s} = δ_{ss'} a^{s^{'}} .

(A3)

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

3. Derivatives

For spatial derivatives, we use the notation that is standard, for example, in general relativity.⁶⁹ In particular, ${X^{μ}}_{; μ}$ is the divergence, namely,

{X^{μ}}_{; μ} = \frac{1}{\sqrt{g_{⋄}}} {(\sqrt{g_{⋄}} X^{μ})}_{, μ} .

(A4)

Here, $g_{⋄} ≐ | \det g |$ ⁠,

f_{, μ} \equiv \partial_{μ} f ≐ \frac{\partial f}{\partial x^{μ}}, f_{, μ ν} ≐ \frac{\partial^{2} f}{\partial x^{μ} \partial x^{ν}} .

(A5)

For functions of the form $f (x, k)$ ⁠, we also introduce the following partial derivatives:

f_{| μ} ≐ \frac{\partial f (x, k)}{\partial x^{μ}}, f_{| μ ν} ≐ \frac{\partial^{2} f (x, k)}{\partial x^{μ} \partial x^{ν}},

(A6)

f^{| μ} ≐ \frac{\partial f (x, k)}{\partial k_{μ}}, f^{| μ ν} ≐ \frac{\partial^{2} f (x, k)}{\partial k_{μ} \partial k_{ν}} .

(A7)

Since $k = \nabla θ (x)$ ⁠, one has

f_{, μ} (x, k (x)) = f_{| μ} + f^{| ν} k_{ν | μ},

(A8)

k_{ν | μ} = k_{ν, μ} = k_{μ, ν} = θ_{, μ ν} = θ_{, ν μ} .

(A9)

4. Inner products

The inner product for scalar waves on Mⁿ is

⟨ Ψ | Φ ⟩ = \int d^{n} x \sqrt{g_{⋄} (x)} Ψ^{*} (x) Φ (x),

(A10)

and for m-dimensional vector waves on Mⁿ,

{⟨ Ψ | Φ ⟩}_{m} ≐ \int d^{n} x \sqrt{g_{⋄} (x)} γ_{ab} (x) Ψ^{a *} (x) Φ^{b} (x) .

(A11)

In particular, in the x representation,

Ψ^{†} (x) ≐ ⟨ Ψ | x ⟩ = (Ψ_{1}^{*}, Ψ_{2}^{*}, \dots, Ψ_{m}^{*}) .

(A12)

We also use the dot product

Ψ \cdot Φ ≐ Ψ^{†} (x) Φ (x) = Ψ_{a}^{*} (x) Φ^{a} (x) .

(A13)

The inner product on the transverse space $M_{⊥}^{n - 1}$ is defined for scalar fields as

{⟨ ψ | ϕ ⟩}^{⊥} ≐ \int d^{n - 1} ϱ \sqrt{h_{⋄} (ζ, ϱ)} ψ^{*} (ζ, ϱ) ϕ (ζ, ϱ)

(A14)

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

{⟨ a | b ⟩}_{N}^{⊥} ≐ \int d^{n - 1} ϱ \sqrt{h_{⋄} (ζ, ϱ)} a_{s}^{*} (ζ, ϱ) b^{s} (ζ, ϱ) .

(A15)

APPENDIX B: AUXILIARY CALCULATIONS

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

1. Calculation of $\hat{K}$ and $U (x)$

We start with Eq. (167) for $\hat{K}$ ⁠. This equation can be rewritten as follows:

\begin{matrix} \hat{K} a = \frac{1}{2} Ξ^{+} ({\hat{p}}_{μ} V_{H}^{μ} + V_{H}^{μ} {\hat{p}}_{μ} - i ℓ_{μ} V_{H}^{μ}) Ξ a + Ξ^{+} C_{H} Ξ \\ = - \frac{i}{2 g_{⋄}^{1 / 4}} [Ξ^{+} {(g_{⋄}^{1 / 4} V_{H}^{μ} Ξ a)}_{, μ} + Ξ^{+} V_{H}^{μ} {(g_{⋄}^{1 / 4} Ξ a)}_{, μ}] - \frac{i}{2} Ξ^{+} ℓ_{μ} V_{H}^{μ} Ξ a + Ξ^{+} C_{H} Ξ \\ = - \frac{i}{g_{⋄}^{1 / 4}} {(g_{⋄}^{1 / 4})}_{, μ} Ξ^{+} V_{H}^{μ} Ξ a - i Ξ^{+} V_{H}^{μ} Ξ a_{, μ} - \frac{i}{2} [Ξ^{+} {(V_{H}^{μ} Ξ)}_{, μ} + Ξ^{+} V_{H}^{μ} Ξ_{, μ}] a - \frac{i}{2} Ξ^{+} ℓ_{μ} V_{H}^{μ} Ξ a + Ξ^{+} C_{H} Ξ \\ = - \frac{i}{2} {(ln \sqrt{g_{⋄}})}_{, μ} V^{μ} a - i V^{μ} a_{, μ} - \frac{i}{2} [{V^{μ}}_{, μ} - {(Ξ^{+})}_{, μ} V_{H}^{μ} Ξ + Ξ^{+} V_{H}^{μ} Ξ_{, μ} + Ξ^{+} ℓ_{μ} V_{H}^{μ} Ξ] a + Ξ^{+} C_{H} Ξ \\ = - i V^{μ} a_{, μ} - \frac{i}{2} {V^{μ}}_{; μ} a - \frac{i}{2} [Ξ^{+} V_{H}^{μ} Ξ_{, μ} - {(Ξ^{+})}_{, μ} V_{H}^{μ} Ξ + Ξ^{+} ℓ_{μ} V_{H}^{μ} Ξ] a + Ξ^{+} C_{H} Ξ, \end{matrix}

(B1)

where we used Eq. (52) and introduced

V^{μ} ≐ Ξ^{+} V_{H}^{μ} Ξ = {(Ξ^{+} D_{H} Ξ)}^{| μ} - Ξ^{+ | μ} D_{H} Ξ - Ξ^{+} D_{H} Ξ^{| μ} = Λ^{| μ} - Ξ^{+ | μ} Ξ Λ - Λ Ξ^{+} Ξ^{| μ} \approx Λ^{| μ},

(B2)

where, at the end, we used $Λ = O (ϵ)$ (Sec. V D 4), so that the last two terms in Eq. (B2) can be neglected.

From Eq. (154), we have ${(Ξ^{+})}_{, μ} = δ^{- 1} Ξ_{, μ}^{H} γ + Ξ^{+} ℓ_{μ}$ ⁠. Using this and also Eq. (130), one obtains

\hat{K} a = - i V^{μ} a_{, μ} - \frac{i}{2} {V^{μ}}_{; μ} a - Ua, U = U_{0} - Ξ^{+} C_{H} Ξ,

(B3)

U_{0} = \frac{i}{2} δ^{- 1} (Ξ^{H} {\underline{V}}_{H}^{μ} Ξ_{, μ} - Ξ_{, μ}^{H} {\underline{V}}_{H}^{μ} Ξ) = δ^{- 1} {(Ξ_{, μ}^{H} {\underline{V}}_{H}^{μ} Ξ)}_{A},

(B4)

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

2. Calculation of $U_{0} (x, k)$

Let us also derive an alternative expression for U₀ [Eq. (B4)] in terms of the partial derivatives $_{| μ}$ instead of the full derivatives $_{, μ}$ ⁠. Using Eq. (59), one obtains $Ξ_{, μ} = Ξ^{| ν} k_{ν, μ} + Ξ_{| μ}$ ⁠. Then, Eq. (B4) can be rewritten as follows:

U_{0} = \frac{i}{2} δ^{- 1} (Ξ^{H} {\underline{V}}_{H}^{μ} Ξ^{| ν} - Ξ^{H | ν} {\underline{V}}_{H}^{μ} Ξ) k_{ν, μ} + \frac{i}{2} δ^{- 1} (Ξ^{H} {\underline{V}}_{H}^{μ} Ξ_{| μ} - Ξ_{| μ}^{H} {\underline{V}}_{H}^{μ} Ξ) .

(B5)

By differentiating Eq. (159) with respect to k_μ, we obtain

V_{H}^{μ} Ξ = Ξ^{| μ} Λ + Ξ Λ^{| μ} - D_{H} Ξ^{| μ} \approx Ξ Λ^{| μ} - D_{H} Ξ^{| μ},

(B6)

where we used that Λ (but not its partial derivatives) is small. Let us multiply this by γ to obtain

{\underline{V}}_{H}^{μ} Ξ \approx γ Ξ Λ^{| μ} - {\underline{D}}_{H} Ξ^{| μ}, Ξ^{H} {\underline{V}}_{H}^{μ} \approx Λ^{| μ} Ξ^{H} γ - Ξ^{H | μ} {\underline{D}}_{H},

(B7)

where the latter equality is the complex conjugate of the former. Using these and $k_{ν, μ} = θ_{, ν μ} = θ_{, μ ν} = O (ϵ)$ ⁠, we further obtain, to the leading order,

\begin{matrix} (Ξ^{H} {\underline{V}}_{H}^{μ} Ξ^{| ν} - Ξ^{H | ν} {\underline{V}}_{H}^{μ} Ξ) k_{ν, μ} \approx (Λ^{| μ} Ξ^{H} γ Ξ^{| ν} - Ξ^{H | μ} {\underline{D}}_{H} Ξ^{| ν} - Ξ^{H | ν} γ Ξ Λ^{| μ} + Ξ^{H | ν} {\underline{D}}_{H} Ξ^{| μ}) k_{ν, μ} \\ = (Λ^{| μ} Ξ^{H} γ Ξ^{| ν} - Ξ^{H | ν} γ Ξ Λ^{| μ}) k_{ν, μ} \\ = (Λ_{, μ} Ξ^{H} γ Ξ^{| μ} - Ξ^{H | μ} γ Ξ Λ_{, μ}) - (Λ_{| μ} Ξ^{H} γ Ξ^{| μ} - Ξ^{H | μ} γ Ξ Λ_{| μ}) \\ \approx - (Λ_{| μ} Ξ^{H} γ Ξ^{| μ} - Ξ^{H | μ} γ Ξ Λ_{| μ}), \end{matrix}

(B8)

where we neglected terms of the second order in ϵ. We also used the fact that $Λ_{s} (x, k (x)) ≲ O (ϵ)$ by the definition of active modes, in addition, Λ_s are changing slowly, and so $Λ_{, μ} ≲ O (ϵ^{2})$ ⁠. (Note that this is only true for $Λ_{, μ}$ but not for $Λ_{| μ}$ ⁠.) Using Eqs. (B7), we also obtain similarly that

\begin{matrix} Ξ^{H} {\underline{V}}_{H}^{μ} Ξ_{| μ} - Ξ_{| μ}^{H} {\underline{V}}_{H}^{μ} Ξ \approx Λ^{| μ} Ξ^{H} γ Ξ_{| μ} - Ξ^{H | μ} {\underline{D}}_{H} Ξ_{| μ} \\ - Ξ_{| μ}^{H} γ Ξ Λ^{| μ} + Ξ_{| μ}^{H} {\underline{D}}_{H} Ξ^{| μ} . \end{matrix}

(B9)

Then,

\begin{matrix} U_{0} = δ^{- 1} (Λ_{| μ} Ξ^{H} γ Ξ^{| μ} - Ξ^{H | μ} γ Ξ Λ_{| μ} - Λ^{| μ} Ξ^{H} γ Ξ_{| μ} + Ξ^{H | μ} {\underline{D}}_{H} Ξ_{| μ} + Ξ_{| μ}^{H} γ Ξ Λ^{| μ} - Ξ_{| μ}^{H} {\underline{D}}_{H} Ξ^{| μ}) / (2 i) \\ = δ^{- 1} {(Λ_{| μ} Ξ^{H} γ Ξ^{| μ} - Λ^{| μ} Ξ^{H} γ Ξ_{| μ} + Ξ^{H | μ} {\underline{D}}_{H} Ξ_{| μ})}_{A} . \end{matrix}

(B10)

As a reminder, the subscript A denotes the anti-Hermitian part, and $δ^{- 1}$ is a unit matrix which only raises the mode index. Equation (B10) is similar to the corresponding result derived in Refs. 24–26. Extra terms are included in those papers in the expressions for U₀, but they are of the second order in ϵ and could have been neglected.

3. Calculation of $\hat{G}$

Consider Eq. (181) for $\hat{G}$ derived under the assumption of the quasioptical ordering. Here, we simplify that equation and derive the explicit formula (182). First, note that

Ξ^{+} {\hat{D}}_{H 2} Ξ a \approx \frac{1}{2} Ξ^{+} {\hat{p}}_{μ} Θ_{H}^{μ ν} {\hat{p}}_{ν} Ξ a \approx - \frac{1}{2 h_{⋄}^{1 / 4}} {[Ξ^{+} {D_{H}}^{| σ \tilde{σ}} Ξ {(h_{⋄}^{1 / 4} a)}_{, \tilde{σ}}]}_{, σ},

(B11)

where we used Eq. (129). Also, using Eq. (128) for ${\hat{D}}_{H 1}$ ⁠, one obtains

\begin{matrix} Ξ^{+} {\hat{D}}_{H 1} \bar{Ξ} {\bar{Λ}}^{- 1} {\bar{Ξ}}^{+} {\hat{D}}_{H 1} Ξ a \approx - (Ξ^{+} V_{H}^{σ} \bar{Ξ} {\bar{Λ}}^{- 1} {\bar{Ξ}}^{+} V_{H}^{\tilde{σ}} Ξ) a_{, σ \tilde{σ}} \\ = - [Ξ^{+} V_{H}^{(σ} \bar{Ξ} {\bar{Λ}}^{- 1} {\bar{Ξ}}^{+} V_{H}^{\tilde{σ})} Ξ] a_{, σ \tilde{σ}} \\ \approx - h_{⋄}^{- 1 / 4} {[Ξ^{+} V_{H}^{(σ} \bar{Ξ} {\bar{Λ}}^{- 1} {\bar{Ξ}}^{+} V_{H}^{\tilde{σ})} Ξ {(h_{⋄}^{1 / 4} a)}_{, \tilde{σ}}]}_{, σ} . \end{matrix}

(B12)

Here, in the first line, we substituted Eq. (128) for ${\hat{D}}_{H 1}$ and also Eq. (14) for ${\hat{p}}_{μ}$ ⁠. Then, we only kept the terms involving derivatives of a along the perpendicular direction of the wave beam. The parentheses in the indices denote symmetrization, namely, for any A and B,

A^{(σ} B^{\tilde{σ})} = \frac{1}{2} (A^{σ} B^{\tilde{σ}} + A^{\tilde{σ}} B^{σ}) .

(B13)

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

\hat{G} \approx - \frac{1}{2 h_{⋄}^{1 / 4}} {[Φ^{σ \tilde{σ}} {(h_{⋄}^{1 / 4} a)}_{, \tilde{σ}}]}_{, σ}, Φ^{σ \tilde{σ}} = Ξ^{+} {D_{H}}^{| σ \tilde{σ}} Ξ - 2 Ξ^{+} {D_{H}}^{(| σ} \bar{Ξ} {\bar{Λ}}^{- 1} {\bar{Ξ}}^{+} {D_{H}}^{| \tilde{σ})} Ξ .

(B14)

In order to simplify the expression for $Φ^{σ \tilde{σ}}$ ⁠, note that

{\bar{Ξ}}^{+} {D_{H}}^{| \tilde{σ}} Ξ = {({\bar{Ξ}}^{+} D_{H} Ξ)}^{| \tilde{σ}} - {\bar{Ξ}}^{+ | \tilde{σ}} D_{H} Ξ - {\bar{Ξ}}^{+} D_{H} Ξ^{| \tilde{σ}} = {({\bar{Ξ}}^{+} Ξ Λ)}^{| \tilde{σ}} - {\bar{Ξ}}^{+ | \tilde{σ}} Ξ Λ - \bar{Λ} {\bar{Ξ}}^{+} Ξ^{| \tilde{σ}} \approx \bar{Λ} {\bar{Ξ}}^{+ | \tilde{σ}} Ξ,

(B15)

where we used ${\bar{Ξ}}^{+} Ξ = 0$ (and thus ${\bar{Ξ}}^{+} Ξ^{| \tilde{σ}} = - {\bar{Ξ}}^{+ | \tilde{σ}} Ξ$ ⁠) and neglected terms proportional to $Λ = O (ϵ_{| |})$ ⁠. Hence,

Ξ^{+} {D_{H}}^{| σ} \bar{Ξ} {\bar{Λ}}^{- 1} {\bar{Ξ}}^{+} {D_{H}}^{| \tilde{σ}} Ξ = Ξ^{+} {\bar{Ξ}}^{| σ} \bar{Λ} {\bar{Λ}}^{- 1} \bar{Λ} {\bar{Ξ}}^{+ | \tilde{σ}} Ξ = Ξ^{+} {\bar{Ξ}}^{| σ} \bar{Λ} {\bar{Ξ}}^{+ | \tilde{σ}} Ξ .

(B16)

Using Eqs. (152) for Ξ and $\bar{Ξ}$ like we did above, Eq. (156) for $D_{H}$ ⁠, and the fact that $Λ = O (ϵ)$ ⁠, we finally obtain

\begin{matrix} Φ^{σ \tilde{σ}} = Ξ^{+} {(Ξ Λ Ξ^{+} + \bar{Ξ} \bar{Λ} {\bar{Ξ}}^{+})}^{| σ \tilde{σ}} Ξ - 2 Ξ^{+} {\bar{Ξ}}^{(| σ} \bar{Λ} {\bar{Ξ}}^{+ | \tilde{σ})} Ξ \\ = Ξ^{+} {(Ξ^{| σ} Λ Ξ^{+} + Ξ Λ^{| σ} Ξ^{+} + Ξ Λ Ξ^{+ | σ} + {\bar{Ξ}}^{| σ} \bar{Λ} {\bar{Ξ}}^{+} + \bar{Ξ} {\bar{Λ}}^{| σ} {\bar{Ξ}}^{+} + \bar{Ξ} \bar{Λ} {\bar{Ξ}}^{+ | σ})}^{| \tilde{σ}} Ξ - 2 Ξ^{+} {\bar{Ξ}}^{(| σ} \bar{Λ} {\bar{Ξ}}^{+ | \tilde{σ})} Ξ \\ \approx Ξ^{+} (Ξ^{| σ} Λ^{| \tilde{σ}} Ξ^{+} + Ξ^{| \tilde{σ}} Λ^{| σ} Ξ^{+} + Ξ Λ^{| σ \tilde{σ}} Ξ^{+} + Ξ Λ^{| σ} Ξ^{+ | \tilde{σ}} + Ξ Λ^{| \tilde{σ}} Ξ^{+ | σ} + {\bar{Ξ}}^{| σ} \bar{Λ} {\bar{Ξ}}^{+ | \tilde{σ}} + {\bar{Ξ}}^{| \tilde{σ}} \bar{Λ} {\bar{Ξ}}^{+ | σ}) Ξ - 2 Ξ^{+} {\bar{Ξ}}^{(| σ} \bar{Λ} {\bar{Ξ}}^{+ | \tilde{σ})} Ξ \\ = Λ^{| σ \tilde{σ}} + Ξ^{+} Ξ^{| σ} Λ^{| \tilde{σ}} + Ξ^{+} Ξ^{| \tilde{σ}} Λ^{| σ} - Λ^{| σ} Ξ^{+} Ξ^{| \tilde{σ}} - Λ^{| \tilde{σ}} Ξ^{+} Ξ^{| σ} \\ \approx Λ^{| σ \tilde{σ}} + [Ξ^{+} Ξ^{(| σ}, V^{\tilde{σ})}], \end{matrix}

(B17)

where we substituted $Λ^{| σ} \approx V^{σ}$ [Eq. (B2)]. The quasioptical approximation implies that $V^{\tilde{σ}}$ is close to a scalar matrix, and so the commutator $[Ξ^{+} Ξ^{(| σ}, V^{\tilde{σ})}]$ can be neglected. Hence, $Φ^{σ \tilde{σ}} \approx Λ^{| σ \tilde{σ}}$ ⁠.

References

1.

E. R.

Tracy

,

A. J.

Brizard

,

A. S.

Richardson

, and

A. N.

Kaufman

,

Ray Tracing and Beyond: Phase Space Methods in Plasma Wave Theory

(

Cambridge University Press

,

New York

,

2014

).

Google Scholar

Crossref

2.

Yu. A.

Kravtsov

and

Yu. I.

Orlov

,

Geometrical Optics of Inhomogeneous Media

(

Springer-Verlag

,

New York

,

1990

).

Google Scholar

Crossref

3.

G. B.

Whitham

,

Linear and Nonlinear Waves

(

Wiley

,

New York

,

1974

).

Google Scholar

4.

T. H.

Stix

,

Waves in Plasmas

(

AIP

,

New York

,

1992

).

Google Scholar

5.

P. T.

Bonoli

, “

Review of recent experimental and modeling progress in the lower hybrid range of frequencies at ITER relevant parameters

,”

Phys. Plasmas

21

,

061508

(

2014

).

https://doi.org/10.1063/1.4884360

Google Scholar

Crossref

6.

See http://www.compxco.com/genray.html for the GENRAY code by A. P. Smirnov and R. W. Harvey.

7.

N. B.

Marushchenko

,

Y.

Turkin

, and

H.

Maassberg

, “

Ray-tracing code TRAVIS for ECR heating, EC current drive and ECE diagnostic

,”

Comput. Phys. Commun.

185

,

165

(

2014

).

https://doi.org/10.1016/j.cpc.2013.09.002

Google Scholar

Crossref

8.

T. I.

Tsujimura

,

S.

Kubo

,

H.

Takahashi

,

R.

Makino

,

R.

Seki

,

Y.

Yoshimura

,

H.

Igami

,

T.

Shimozuma

,

K.

Ida

,

C.

Suzuki

,

M.

Emoto

,

M.

Yokoyama

,

T.

Kobayashi

,

C.

Moon

,

K.

Nagaoka

,

M.

Osakabe

,

S.

Kobayashi

,

S.

Ito

,

Y.

Mizuno

,

K.

Okada

,

A.

Ejiri

,

T.

Mutoh

, and

LHD Experiment Group,

“

Development and application of a ray-tracing code integrating with 3D equilibrium mapping in LHD ECH experiments

,”

Nucl. Fusion

55

,

123019

(

2015

).

https://doi.org/10.1088/0029-5515/55/12/123019

Google Scholar

Crossref

9.

E.

Mazzucato

, “

Propagation of a Gaussian beam in a nonhomogeneous plasma

,”

Phys. Fluids B

1

,

1855

(

1989

).

https://doi.org/10.1063/1.858917

Google Scholar

Crossref

10.

D.

Farina

, “

A quasi-optical beam-tracing code for electron cyclotron absorption and current drive: GRAY

,”

Fusion Sci. Technol.

52

,

154

(

2007

).

https://doi.org/10.13182/FST07-A1494

Google Scholar

Crossref

11.

G. V.

Pereverzev

, “

Use of the multidimensional WKB method to describe propagation of lower hybrid waves in tokamak plasma

,”

Nucl. Fusion

32

,

1091

(

1992

).

https://doi.org/10.1088/0029-5515/32/7/I01

Google Scholar

Crossref

12.

G. V.

Pereverzev

, “

Beam tracing in inhomogeneous anisotropic plasmas

,”

Phys. Plasmas

5

,

3529

(

1998

).

https://doi.org/10.1063/1.873070

Google Scholar

Crossref

13.

E.

Poli

,

A. G.

Peeters

, and

G. V.

Pereverzev

, “

TORBEAM, a beam tracing code for electron-cyclotron waves in tokamak plasmas

,”

Comput. Phys. Commun.

136

,

90

(

2001

).

https://doi.org/10.1016/S0010-4655(01)00146-1

Google Scholar

Crossref

14.

E.

Poli

,

A.

Bock

,

M.

Lochbrunner

,

O.

Maj

,

M.

Reich

,

A.

Snicker

,

A.

Stegmeir

,

F.

Volpe

,

N.

Bertelli

,

R.

Bilato

,

G. D.

Conway

,

D.

Farina

,

F.

Felici

,

L.

Figini

,

R.

Fischer

,

C.

Galperti

,

T.

Happel

,

Y. R.

Lin-Liu

,

N. B.

Marushchenko

,

U.

Mszanowski

,

F. M.

Poli

,

J.

Stober

,

E.

Westerhof

,

R.

Zille

,

A. G.

Peeters

, and

G. V.

Pereverzev

, “

TORBEAM 2.0, a paraxial beam tracing code for electron-cyclotron beams in fusion plasmas for extended physics applications

,”

Comput. Phys. Commun.

225

,

36

(

2018

).

https://doi.org/10.1016/j.cpc.2017.12.018

Google Scholar

Crossref

15.

It is not a purpose of this paper to provide a comprehensive list of existing approaches to reduced modeling of EM waves in plasmas, so the above reference list is not intended as exhaustive.

16.

A. A.

Balakin

,

M. A.

Balakina

,

G. V.

Permitin

, and

A. I.

Smirnov

, “

Quasi-optical description of wave beams in smoothly inhomogeneous anisotropic media

,”

J. Phys. D: Appl. Phys.

40

,

4285

(

2007

).

https://doi.org/10.1088/0022-3727/40/14/028

Google Scholar

Crossref

17.

A. A.

Balakin

,

M. A.

Balakina

,

G. V.

Permitin

, and

A. I.

Smirnov

, “

Scalar equation for wave beams in a magnetized plasma

,”

Plasma Phys. Rep.

33

,

302

(

2007

).

https://doi.org/10.1134/S1063780X07040058

Google Scholar

Crossref

18.

For applications and comparison with other codes, see also

A. A.

Balakin

,

M. A.

Balakina

, and

E.

Westerhof

, “

ECRH power deposition from a quasioptical point of view

,”

Nucl. Fusion

48

,

065003

(

2008

);

https://doi.org/10.1088/0029-5515/48/6/065003

Crossref

Google Scholar

O.

Maj

,

A. A.

Balakin

, and

E.

Poli

, “

Effects of aberration on paraxial wave beams beam tracing versus quasioptical solutions

,”

Plasma Phys. Controlled Fusion

52

,

085006

(

2010

).

https://doi.org/10.1088/0741-3335/52/8/085006

Crossref

Google Scholar

19.

There have been works on quasioptical modeling for systems where wave beams are naturally aligned with an axis of symmetry. For example, see

C. K.

Phillips

,

F. W.

Perkins

, and

D. Q.

Hwang

, “

Parabolic approximation method for fast magnetosonic wave propagation in tokamaks

,”

Phys. Fluids

29

,

1608

(

1986

). However, the methods used there do not easily extend to beam modeling in more general settings.

https://doi.org/10.1063/1.865677

Google Scholar

Crossref

20.

This does not extend to ray-tracing codes. For example, a ray-tracing code raycon¹ captures mode conversion by using an approximate yet somewhat generic analytic solution that applies in the limit of a narrow region of resonant coupling. A related ray-tracing model with mode conversion was also proposed in Refs. 24 and 25.

21.

Mode conversion is understood here as the exchange of action between near-resonant but nonintersecting branches of the dispersion relation.¹ Mathematically, this process is similar to quantum tunneling. Adiabatic transformations of a wave mode along continuous dispersion curves, such as the X–B transformation,⁴ are not mode conversion in the true sense of the term and can also be simulated using already existing reduced models.

22.

K. Y.

Bliokh

,

F. J.

Rodríguez-Fortuño

,

F.

Nori

, and

A. V.

Zayats

, “

Spin-orbit interactions of light

,”

Nat. Photonics

9

,

796

(

2015

).

https://doi.org/10.1038/nphoton.2015.201

Google Scholar

Crossref

23.

M. A.

Oancea

,

C. F.

Paganini

,

J.

Joudioux

, and

L.

Andersson

, “

An overview of the gravitational spin Hall effect

,” e-print arXiv:1904.09963.

24.

D. E.

Ruiz

, “

Geometric theory of waves and its applications to plasma physics

,” Ph.D. thesis (Princeton University,

2017

);

Google Scholar

D. E.

Ruiz

, e-print arXiv:1708.05423.

25.

D. E.

Ruiz

and

I. Y.

Dodin

, “

Extending geometrical optics: A Lagrangian theory for vector waves

,”

Phys. Plasmas

24

,

055704

(

2017

).

https://doi.org/10.1063/1.4977537

Google Scholar

Crossref

26.

D. E.

Ruiz

and

I. Y.

Dodin

, “

Lagrangian geometrical optics of nonadiabatic vector waves and spin particles

,”

Phys. Lett. A

379

,

2337

(

2015

).

https://doi.org/10.1016/j.physleta.2015.07.038

Google Scholar

Crossref

27.

D. E.

Ruiz

and

I. Y.

Dodin

, “

First-principles variational formulation of polarization effects in geometrical optics

,”

Phys. Rev. A

92

,

043805

(

2015

).

https://doi.org/10.1103/PhysRevA.92.043805

Google Scholar

Crossref

28.

D. E.

Ruiz

,

C. L.

Ellison

, and

I. Y.

Dodin

, “

Relativistic ponderomotive Hamiltonian of a Dirac particle in a vacuum laser field

,”

Phys. Rev. A

92

,

062124

(

2015

).

https://doi.org/10.1103/PhysRevA.92.062124

Google Scholar

Crossref

29.

I. Y.

Dodin

,

D. E.

Ruiz

, and

S.

Kubo

, “

Mode conversion in cold low-density plasma with a sheared magnetic field

,”

Phys. Plasmas

24

,

122116

(

2017

).

https://doi.org/10.1063/1.5003931

Google Scholar

Crossref

30.

This also facilitates the aforementioned applications to related problems in general relativity,²³ as will be discussed in future publications.

31.

K.

Yanagihara

,

I. Y.

Dodin

, and

S.

Kubo

, “

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

,”

Phys. Plasmas

26

,

072111

(

2019

).

https://doi.org/10.1063/1.5095173

Google Scholar

Crossref

32.

K.

Yanagihara

,

I. Y.

Dodin

, and

S.

Kubo

, “

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

,”

Phys. Plasmas

26

,

072112

(

2019

).

https://doi.org/10.1063/1.5095174

Google Scholar

Crossref

33.

T.

Notake

,

S.

Kubo

,

T.

Shimozuma

,

H.

Idei

,

Y.

Yoshimura

,

S.

Inagaki

,

K.

Ohkubo

,

S.

Kobayashi

,

Y.

Mizuno

,

S.

Ito

,

Y.

Takita

,

T.

Watari

,

K.

Narihara

,

T.

Morisaki

,

I.

Yamada

,

Y.

Nagayama

,

K.

Tanaka

,

S.

Sakakibara

,

R.

Kumazawa

,

T.

Seki

,

K.

Saito

,

T.

Mutoh

,

A.

Shimizu

,

A.

Komori

, and

LHD Experiment Group,

“

Optimization of incident wave polarization for ECRH in LHD

,”

Plasma Phys. Controlled Fusion

47

,

531

(

2005

).

https://doi.org/10.1088/0741-3335/47/3/009

Google Scholar

Crossref

34.

S.

Kubo

,

H.

Igami

,

T. I.

Tsujimura

,

T.

Shimozuma

,

H.

Takahashi

,

Y.

Yoshimura

,

M.

Nishiura

,

R.

Makino

, and

T.

Mutoh

, “

Plasma interface of the EC waves to the LHD peripheral region

,”

AIP Conf. Proc.

1689

,

090006

(

2015

).

https://doi.org/10.1063/1.4936543

Google Scholar

Crossref

35.

A.

Iiyoshi

,

A.

Komori

,

A.

Ejiri

,

M.

Emoto

,

H.

Funaba

,

M.

Goto

,

K.

Ida

,

H.

Idei

,

S.

Inagaki

,

S.

Kado

,

O.

Kaneko

,

K.

Kawahata

,

T.

Kobuchi

,

S.

Kubo

,

R.

Kumazawa

,

S.

Masuzaki

,

T.

Minami

,

J.

Miyazawa

,

T.

Morisaki

,

S.

Morita

,

S.

Murakami

,

S.

Muto

,

T.

Mutoh

,

Y.

Nagayama

,

Y.

Nakamura

,

H.

Nakanishi

,

K.

Narihara

,

K.

Nishimura

,

N.

Noda

,

S.

Ohdachi

,

N.

Ohyabu

,

Y.

Oka

,

M.

Osakabe

,

T.

Ozaki

,

B. J.

Peterson

,

A.

Sagara

,

S.

Sakakibara

,

R.

Sakamoto

,

H.

Sasao

,

M.

Sasao

,

K.

Sato

,

M.

Sato

,

T.

Seki

,

T.

Shimozuma

,

M.

Shoji

,

H.

Suzuki

,

Y.

Takeiri

,

K.

Tanaka

,

K.

Toi

,

T.

Tokuzawa

,

K.

Tsumori

,

K.

Tsuzuki

,

K. Y.

Watanabe

,

T.

Watari

,

H.

Yamada

,

I.

Yamada

,

S.

Yamaguchi

,

M.

Yokoyama

,

R.

Akiyama

,

H.

Chikaraishi

,

K.

Haba

,

S.

Hamaguchi

,

M.

Iima

,

S.

Imagawa

,

N.

Inoue

,

K.

Iwamoto

,

S.

Kitagawa

,

J.

Kodaira

,

Y.

Kubota

,

R.

Maekawa

,

T.

Mito

,

T.

Nagasaka

,

A.

Nishimura

,

C.

Takahashi

,

K.

Takahata

,

Y.

Takita

,

H.

Tamura

,

T.

Tsuzuki

,

S.

Yamada

,

K.

Yamauchi

,

N.

Yanagi

,

H.

Yonezu

,

Y.

Hamada

,

K.

Matsuoka

,

K.

Murai

,

K.

Ohkubo

,

I.

Ohtake

,

M.

Okamoto

,

S.

Satoh

,

T.

Satow

,

S.

Sudo

,

S.

Tanahashi

,

K.

Yamazaki

,

M.

Fujiwara

, and

O.

Motojima

, “

Overview of the Large Helical Device project

,”

Nucl. Fusion

39

,

1245

(

1999

).

https://doi.org/10.1088/0029-5515/39/9Y/313

Google Scholar

Crossref

36.

M.

Osakabe

,

Y.

Takeiri

,

T.

Morisaki

,

G.

Motojima

,

K.

Ogawa

,

M.

Isobe

,

M.

Tanaka

,

S.

Murakami

,

A.

Shimizu

,

K.

Nagaoka

,

H.

Takahashi

,

K.

Nagasaki

,

H.

Takahashi

,

T.

Fujita

,

Y.

Oya

,

M.

Sakamoto

,

Y.

Ueda

,

T.

Akiyama

,

H.

Kasahara

,

S.

Sakakibara

,

R.

Sakamoto

,

M.

Tokitani

,

H.

Yamada

,

M.

Yokoyama

,

Y.

Yoshimura

, and

LHD Experiment Group

, “

Current status of large helical device and its prospect for deuterium experiment

,”

Fusion Sci. Technol.

72

,

199

(

2017

).

https://doi.org/10.1080/15361055.2017.1335145

Google Scholar

Crossref

37.

For a more rigorous definition of the GO approximation, see, for instance, Ref. 2 or

I. Y.

Dodin

,

A. I.

Zhmoginov

, and

D. E.

Ruiz

, “

Variational principles for dissipative (sub)systems, with applications to the theory of linear dispersion and geometrical optics

,”

Phys. Lett. A

381

,

1411

(

2017

).

https://doi.org/10.1016/j.physleta.2017.02.023

Google Scholar

Crossref

38.

D. E.

Ruiz

and

I. Y.

Dodin

, “

Ponderomotive dynamics of waves in quasiperiodically modulated media

,”

Phys. Rev. A

95

,

032114

(

2017

).

https://doi.org/10.1103/PhysRevA.95.032114

Google Scholar

Crossref

39.

H.

Weyl

,

The Theory of Groups and Quantum Mechanics

(

Dover

,

New York

,

1931

).

Google Scholar

40.

E.

Wigner

, “

On the quantum correction for the thermodynamic equilibrium

,”

Phys. Rev.

40

,

749

(

1932

).

https://doi.org/10.1103/PhysRev.40.749

Google Scholar

Crossref

41.

J. E.

Moyal

, “

Quantum mechanics as a statistical theory

,”

Proc. Cambridge Philos. Soc.

45

,

99

(

1949

).

https://doi.org/10.1017/S0305004100000487

Google Scholar

Crossref

42.

R. G.

Littlejohn

, “

The semiclassical evolution of wave packets

,”

Phys. Rep.

138

,

193

(

1986

).

https://doi.org/10.1016/0370-1573(86)90103-1

Google Scholar

Crossref

43.

S. W.

McDonald

, “

Phase-space representations of wave equations with applications to the eikonal approximation for short-wavelength waves

,”

Phys. Rep.

158

,

337

(

1988

).

https://doi.org/10.1016/0370-1573(88)90012-9

Google Scholar

Crossref

44.

R. G.

Littlejohn

and

W. G.

Flynn

, “

Geometric phases in the asymptotic theory of coupled wave equations

,”

Phys. Rev. A

44

,

5239

(

1991

).

https://doi.org/10.1103/PhysRevA.44.5239

Google Scholar

Crossref

PubMed

45.

I. B.

Bernstein

, “

Geometric optics in space and time varying plasmas. I

,”

Phys. Fluids

18

,

320

(

1975

).

https://doi.org/10.1063/1.861140

Google Scholar

Crossref

46.

H. L.

Berk

and

D.

Pfirsch

, “

WKB method for systems of integral equations

,”

J. Math. Phys.

21

,

2054

(

1980

).

https://doi.org/10.1063/1.524716

Google Scholar

Crossref

47.

L.

Friedland

,

G.

Goldner

, and

A. N.

Kaufman

, “

Four-dimensional eikonal theory of linear mode conversion

,”

Phys. Rev. Lett.

58

,

1392

(

1987

).

https://doi.org/10.1103/PhysRevLett.58.1392

Google Scholar

Crossref

PubMed

48.

For Euclidean coordinates, a detailed discussion of the Weyl calculus can be found in Refs. 1 and 24. A generalization to curved metric was considered also in

C.

Gneiting

,

T.

Fischer

, and

K.

Hornberger

, “

Quantum phase-space representation for curved configuration spaces

,”

Phys. Rev. A

88

,

062117

(

2013

) and, to some extent, Ref. 66.

https://doi.org/10.1103/PhysRevA.88.062117

Google Scholar

Crossref

49.

B. S.

DeWitt

, “

Point transformations in quantum mechanics

,”

Phys. Rev.

85

,

653

(

1952

).

https://doi.org/10.1103/PhysRev.85.653

Google Scholar

Crossref

50.

J. M.

Domingos

and

M. H.

Caldeira

, “

Self-adjointness of momentum operators in generalized coordinates

,”

Found. Phys.

14

,

147

(

1984

).

https://doi.org/10.1007/BF00729971

Google Scholar

Crossref

51.

Eventually (Sec. V E), we shall choose $k (x)$ such that $D (x, k (x))$ be small. But that will characterize $k (x)$ rather than the function D itself, which is considered order-one when evaluated on general $(x, p)$ ⁠.

52.

Reduced codes typically assume that the anti-Hermitian part of the dielectric tensor is small. That is not always the case; for example, see

M. D.

Tokman

,

E.

Westerhof

, and

M. A.

Gavrilova

, “

Wave power flux and ray-tracing in regions of resonant absorption

,”

Plasma Phys. Controlled Fusion

42

,

91

(

2000

) and references therein. However, this subject will not be discussed below.

https://doi.org/10.1088/0741-3335/42/2/302

Google Scholar

Crossref

53.

Strictly speaking, Eq. (1) can always be multiplied by a constant matrix to ensure that the new dispersion operator satisfies the ordering (63) at a given location. Hence, the only nontrivial requirement that we actually impose here is that such $D_{A}$ does not grow too much as the wave propagates.

54.

I. Y.

Dodin

and

N. J.

Fisch

, “

Axiomatic geometrical optics, Abraham-Minkowski controversy, and photon properties derived classically

,”

Phys. Rev. A

86

,

053834

(

2012

).

https://doi.org/10.1103/PhysRevA.86.053834

Google Scholar

Crossref

55.

We use an approach similar to the 3 + 1 approach in general relativity. For example, see

I. Y.

Dodin

and

N. J.

Fisch

, “

Vlasov equation and collisionless hydrodynamics adapted to curved spacetime

,”

Phys. Plasmas

17

,

112118

(

2010

) including the supplementary material there.

https://doi.org/10.1063/1.3497005

Google Scholar

Crossref

56.

For example, see

P. D.

Powell

, “

Calculating determinants of block matrices

,” e-print arXiv:1112.4379.

57.

In principle, ${\hat{D}}_{A}$ produces effects qualitatively different from those produced by ${\hat{D}}_{H}$ ⁠, so one might argue that the accuracy of its approximation does not necessarily need to equal that of ${\hat{D}}_{H}$ ⁠. Besides, $D_{A}$ is often more abrupt than $D_{H}$ ⁠, so retaining $V_{A}^{σ}$ in Eq. (82) may in fact lead to a better model in principle. However, in such cases, the Weyl symbol D of the original dispersion operator is typically known only to the zeroth order (Sec. VI C), and thus so is $D_{A}$ ⁠. Hence, retaining first-order corrections the damping coefficient does not seem warranted.

58.

For example, see the derivation of the linear part of the nonlinear Schrödinger equation for a wave in a homogeneous medium in

V. I.

Karpman

,

Non-Linear Waves in Dispersive Media

(

Pergamon Press

,

New York

,

1974

), Sec. 27.

Google Scholar

59.

T. J.

Bridges

and

S.

Reich

, “

Multi-symplectic integrators: Numerical schemes for Hamiltonian PDEs that conserve symplecticity

,”

Phys. Lett. A

284

,

184

(

2001

).

https://doi.org/10.1016/S0375-9601(01)00294-8

Google Scholar

Crossref

60.

L.

Friedland

and

A. N.

Kaufman

, “

Congruent reduction in geometric optics and mode conversion

,”

Phys. Fluids

30

,

3050

(

1987

).

https://doi.org/10.1063/1.866480

Google Scholar

Crossref

61.

M.

Onoda

,

S.

Murakami

, and

N.

Nagaosa

, “

Hall effect of light

,”

Phys. Rev. Lett.

93

,

083901

(

2004

).

https://doi.org/10.1103/PhysRevLett.93.083901

Google Scholar

Crossref

PubMed

62.

K. Y.

Bliokh

,

A.

Niv

,

V.

Kleiner

, and

E.

Hasman

, “

Geometrodynamics of spinning light

,”

Nat. Photonics

2

,

748

(

2008

).

https://doi.org/10.1038/nphoton.2008.229

Google Scholar

Crossref

63.

P.

Gosselin

,

A.

Bérard

, and

H.

Mohrbach

, “

Spin Hall effect of photons in a static gravitational field

,”

Phys. Rev. D

75

,

084035

(

2007

).

https://doi.org/10.1103/PhysRevD.75.084035

Google Scholar

Crossref

64.

C.

Duval

and

T.

Schücker

, “

Gravitational birefringence of light in Robertson–Walker cosmologies

,”

Phys. Rev. D

96

,

043517

(

2017

).

https://doi.org/10.1103/PhysRevD.96.043517

Google Scholar

Crossref

65.

I. Y.

Dodin

, “

Geometric view on noneikonal waves

,”

Phys. Lett. A

378

,

1598

(

2014

).

https://doi.org/10.1016/j.physleta.2014.04.004

Google Scholar

Crossref

66.

L. D.

Landau

and

E. M.

Lifshitz

,

The Classical Theory of Fields

(

Pergamon Press

,

New York

,

1971

).

Google Scholar

67.

L.

Friedland

, “

Gyroresonant absorption from congruent reduction of an anisotropic pressure fluid model

,”

Phys. Fluids

31

,

2615

(

1988

).

https://doi.org/10.1063/1.866539

Google Scholar

Crossref

68.

J.

Preinhaelter

and

V.

Kopecky

, “

Penetration of high-frequency waves into a weakly inhomogeneous magnetized plasma at oblique incidence and their transformation to Bernstein modes

,”

J. Plasma Phys.

10

,

1

(

1973

).

https://doi.org/10.1017/S0022377800007649

Google Scholar

Crossref

69.

C. W.

Misner

,

K. S.

Thorne

, and

J. A.

Wheeler

,

Gravitation

(

Freeman

,

San Francisco

,

1973

).

Google Scholar

2019

Author(s)

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

I. INTRODUCTION

A. Motivation

B. General idea

C. Outline

II. GENERAL PROBLEM

III. ENVELOPE DISPERSION OPERATOR

IV. SCALAR WAVES

A. Weyl calculus

1. Basic definitions

2. Coordinate and momentum operators

3. Wigner–Weyl transform

B. Approximate D̂

C. Approximate envelope equation

D. Leading-order approximation: Geometrical optics

E. Quasioptical model

1. Ray-based coordinates

2. Quasioptical equation

3. Simplified equations

V. VECTOR WAVES

A. Hilbert space for vector waves

B. Operators on Hmn

C. Weyl expansion of the dispersion operator

D. Active and passive modes

1. Basis from the eigenvectors of DH

2. Amplitude vectors

3. Polarization matrices

4. Equation for the active modes

E. Leading-order approximation: Extended geometrical optics

F. Quasioptical model

1. Quasioptical equation

2. Simplified equations

G. Summary and discussion

1. Scalar beams

2. Single-mode vector beams

3. Multimode vector beams

VI. ELECTROMAGNETIC WAVES

A. Covariant formulation

B. Noncovariant formulation

C. Stationary waves

D. Waves in flat space

VII. CONCLUSIONS

SUPPLEMENTARY MATERIAL

ACKNOWLEDGMENTS

APPENDIX A: SUMMARY OF THE SELECTED NOTATIONS

1. Basic symbols

2. Index manipulation

3. Derivatives

4. Inner products

APPENDIX B: AUXILIARY CALCULATIONS

1. Calculation of K̂ and U(x)

2. Calculation of U0(x,k)

3. Calculation of Ĝ

References

Supplementary Material

Citing articles via

Resources

Explore

pubs.aip.org

Connect with AIP Publishing

This Feature Is Available To Subscribers Only

B. Approximate $\hat{D}$

B. Operators on $H_{m}^{n}$

1. Basis from the eigenvectors of $D_{H}$

1. Calculation of $\hat{K}$ and $U (x)$

2. Calculation of $U_{0} (x, k)$

3. Calculation of $\hat{G}$