2.1. Notation

For any given Hermitian operator $\hat {A}$ on a given Hilbert space $\mathbb{H}$ , one can introduce a complete set of its orthogonal eigenvectors ,

(2.1)

where $\lambda$ are the corresponding eigenvalues. (We use $^*$ to denote complex conjugation, so the second equation states that the eigenvalues of a Hermitian operator are real.) Assuming that the eigenspace of $\hat {A}$ is $\mathbb{R} \equiv ({-}\infty , \infty )$ ,Footnote ⁸ its eigenvectors can always be chosen such that

(2.2)

where $\delta$ is the Dirac delta function. Then, they satisfy the identity

(2.3)

where $\hat {1}$ is a unit operator. (Integrals in this paper are taken over $\mathbb{R}$ except where specified otherwise.) Accordingly, $\hat {A}$ can be expressed as

(2.4)

One can also use the basis of to define the trace of any given operator $\hat {O}$ :

(2.5)

Since the trace is invariant under the change of basis, any $\hat {A}$ generates the same ${\textrm {tr}}\;\hat {O}$ .

Among such operators on $\mathbb{H}$ , let us assign some operator $\hat {x}$ as the position operator. For any given vector $| {\psi }\rangle$ in $\mathbb{H}$ , the eigenvectors of $\hat {x}$ define the $x$ -representation of $| {\psi }\rangle$ , or a field $\psi _{\hat {x}}(x)$ , via

(2.6)

Correspondingly, for any given operator $\hat {O}$ , the field associated with is as follows:

(2.7)

Using (2.3) for $\hat {A} = \hat {x}$ , one can rewrite this as

(2.8)

where the kernel is the ‘matrix element’ of $\hat {O}$ in the $x$ -representation, or simply the $x$ -representation of $\hat {O}$ . In the special case when $\hat {O}_{\hat {x}}$ is a local or differential operator (that is, if $(\hat {O}_{\hat {x}}\psi _{\hat {x}})(x)$ can be expressed through $\psi _{\hat {x}}$ and its derivatives $\psi ^{(m)}_{\hat {x}}$ with finite $m$ ), then $O_{\hat {x}}(x, x')$ can be expressed through the delta function, $\delta (x - x')$ , and its corresponding derivatives, $\delta ^{(m)}(x - x')$ .

Using (2.4), one can further use this function to express $\hat {O}$ as follows:

(2.9)

and its Hermitian adjoint as

(2.10)

By definition,

(2.11)

Then, $(\hat {x}_{\hat {x}} \psi _{\hat {x}})(x) = x\psi _{\hat {x}}(x)$ , so

(2.12) \begin{gather} \hat {x}_{\hat {x}} = x. \end{gather}

Let us also introduce the wavevector operator $\hat {k}$ as a generator of spatial translations:

(2.13)

(here, $a$ is any real constant), whence

(2.14)

At $a \to 0$ , one has $\exp ({-}\mathrm{i} a \hat {k}_{\hat {x}}) \to 1 - \mathrm{i} a \hat {k}_{\hat {x}}$ and $\psi _{\hat {x}}(x - a) \to \psi _{\hat {x}}(x) - a\,\partial _x\psi _{\hat {x}}(x)$ , so (2.14) leads to $-\mathrm{i} \hat {k}_{\hat {x}} \psi _{\hat {x}}(x) \to - \partial _x\psi _{\hat {x}}(x)$ for any $\psi _{\hat {x}}(x)$ . Thus,

(2.15) \begin{gather} \hat {k}_{\hat {x}} = -\mathrm{i}\partial _{x} \end{gather}

and $k_{\hat {x}}(x, x') = -\mathrm{i} \delta '(x - x')$ . (The prime in $\delta '$ denotes the derivative with respect to the argument.) Such $\hat {x}$ and $\hat {k}$ are called Schrödinger operators, and here we will also refer to them as base operators. Note that the base operators satisfy the canonical commutation relation

(2.16) \begin{gather} [\hat {x}, \hat {k}] = \mathrm{i}. \end{gather}

The converse is also true: if given Hermitian operators $\hat {x}$ and $\hat {k}$ satisfy (2.16), then there exists a representation, called the $x$ -representation, in which (2.12) and (2.15) are satisfied; see, e.g. Galindo & Pascual (Reference Galindo and Pascual1990, § 2.12).

For , which is the $x$ -representation of an eigenvector of $\hat {k}$ , (2.15) yields $-\mathrm{i}\partial _{x}\xi _k = k \xi_k$ . Assuming the same normalisation as in (2.2), one obtainsFootnote ⁹

(2.17)

where we chose the integration constant such that $\xi _k(0)$ is real. This convention will also be assumed from now on for the eigenfunctions of all other base operators.

Notably, (2.17) means that the $k$ -representation of any given $| {\psi }\rangle$ is the Fourier transform of the corresponding $\psi _{\hat {x}}$ :

(2.18)

Obviously, the $k$ -representation of $\hat {k}$ is $\hat {k}_{\hat {k}} = k$ . The $k$ -representation of $\hat {x}$ is found from

(2.19)

Then, $\hat {x}_{\hat {k}} = \mathrm{i} \partial _k$ , so $-\hat {x}$ can be understood as a generator of translations in the $k$ -space.

We will also use the notation

(2.20) \begin{gather} \hat {{\mathsf{z}}} \equiv \left ( \begin{array}{c} \hat {x}\\ \hat {k} \end{array} \right ) \equiv (\hat {x}, \hat {k})^{\intercal }, \quad {\mathsf{z}} \equiv \left ( \begin{array}{c} x\\ k \end{array} \right ) \equiv (x, k)^{\intercal }, \end{gather}

where $^{\intercal }$ denotes transposition. In other words, ${\mathsf{z}} \equiv (x, k)^{\intercal }$ will be understood as a column vector and ${\mathsf{z}}^{\intercal } \equiv (x, k)$ will be generally understood as a row vector. (We ignored the difference between the two in the introduction for simplicity.) The upper and lower indices are not distinguished in this paper, but the Einstein summation convention for repeated indices is assumed. For two vectors ${\mathsf{z}}_1 = (a_1, b_1)^{\intercal }$ and ${\mathsf{z}}_2 = (a_2, b_2)^{\intercal }$ , one has ${\mathsf{z}}_1^{\intercal } {\mathsf{ z}}_2 = a_1 a_2 + b_1 b_2 \equiv {\mathsf{z}}_1 \boldsymbol{\cdot} {\mathsf{z}}_2 = {\mathsf{z}}_2 \boldsymbol{\cdot} {\mathsf{ z}}_1$ .Footnote ¹⁰ We will also need a wedge product defined as

(2.21) \begin{gather} {\mathsf{a}} \wedge {\mathsf{b}} \doteq {\mathsf{a}} \boldsymbol{\cdot} {\mathsf{J}} {\mathsf{b}} = - {\mathsf{b}} \wedge {\mathsf{a}}, \end{gather}

and ${\mathsf{J}}$ is the canonical symplectic form, which satisfies ${\mathsf{J}}{\mathsf{ J}} = - {{\mathsf{I}}}$ :

(2.22) \begin{gather} {\mathsf{J}} \doteq \left ( \begin{array}{c@{\quad}c} 0 & 1\\ -1 & 0\\ \end{array} \right )\!, \quad {{\mathsf{I}}} \doteq \left ( \begin{array}{c@{\quad}c} 1 & 0\\ 0 & 1\\ \end{array} \right )\!. \end{gather}

2.2. Weyl symbol calculus

Any given operator $\hat {O}$ can be expressed as a linear superposition of non-commuting phase-space translations $\hat {T}_{{{\mathsf{a}}}} \doteq \mathrm{e}^{-\mathrm{i} {{\mathsf{a}}} \wedge \hat {{\mathsf{z}}}}$ . Specifically, this can be expressed as follows:Footnote ¹¹

(2.23) \begin{gather} \hat {O} = \frac {1}{2\pi } \int {\mathrm{d}}^2{\mathsf{z}}\, O({\mathsf{z}}) \hat {\varDelta }_{{\mathsf{z}}}, \quad O({\mathsf{z}}) = {\textrm {tr}}(\varDelta _{\mathsf{z}} \hat {A}), \end{gather}

where the weight function $O({\mathsf{z}}) \equiv {\textrm {symb}}\; \hat {O}$ is called the Weyl symbol (or just ‘symbol’Footnote ¹² ) of $\hat {O}$ and $\hat {\varDelta }_{{\mathsf{z}}}$ is the Wigner operator defined as

(2.24) \begin{gather} \hat {\varDelta }_{{\mathsf{z}}} \doteq \frac {1}{2\pi } \int {\mathrm{d}}^2{{\mathsf{a}}}\, \mathrm{e}^{\mathrm{i} {{\mathsf{a}}} \wedge {\mathsf{z}}}\,\hat {T}_{{{\mathsf{a}}}}. \end{gather}

Since $\hat {\varDelta }_{{\mathsf{z}}}$ is expressed through $\hat {T}_{{{\mathsf{a}}}}$ , and $\hat {T}_{{{\mathsf{a}}}}$ is expressed through $\hat {{\mathsf{z}}}$ , (2.23) also means that $\hat {O}$ is expressible through $\hat {{\mathsf{z}}}$ and can be symbolically written as

(2.25) \begin{gather} \hat {O} = \mathfrak{O}(\hat {{\mathsf{z}}}) \equiv \mathfrak{O}(\hat {x}, \hat {k}). \end{gather}

In general, (2.25) is just a shorthand for (2.23). If $\mathfrak{O}$ is representable as a sum $\mathfrak{O}_x(\hat {x}) + \mathfrak{O}_k(\hat {k})$ , then $\mathfrak{O}$ is identical to the operator symbol, $O$ , except it is evaluated not on $x$ and $k$ but on the corresponding operators. The $x$ - and $k$ -representations of (2.25) are

(2.26) \begin{gather} \hat {O}_{\hat {x}} = \mathfrak{O}(x, -\mathrm{i} \partial _x), \quad \hat {O}_{\hat {k}} = \mathfrak{O}(\mathrm{i} \partial _k, k). \end{gather}

For differential operators, the corresponding $\mathfrak{O}$ are polynomials in $-\mathrm{i} \partial _x$ (or in $\mathrm{i} \partial _k$ , depending on the context) and are typically known explicitly in a given problem. More generally, (2.26) can be understood as a symbolic pseudo-differential representation of $\hat {Q}$ .

Equation (2.23) can also be represented as

(2.27a)

(2.27b)

and the symbol can be expressed as follows:

(2.28a)

(2.28b)

Note that the symbol of an operator is related to its $x$ -representation via the Fourier transform with respect to the argument difference $x - x'$ at fixed $(x + x')/2$ :

(2.29)

and similarly for the $k$ -representation. This also leads to the following notable properties of Weyl symbols:

(2.30a)

(2.30b)

An operator unambiguously determines its symbol and vice versa. We denote this isomorphism as $\hat {O} \Leftrightarrow O$ . The mapping $\hat {O} \Rightarrow O$ is called the Wigner transform and $O \Rightarrow \hat {O}$ is called the Weyl transform. For uniformity, we will refer to these as the direct and inverse Wigner–Weyl transform. The Wigner–Weyl isomorphism $\Leftrightarrow$ is natural in that it has the following properties:

(2.31) \begin{align} \hat {1} \Leftrightarrow 1, \quad \hat {x} \Leftrightarrow x, \quad \hat {k} \Leftrightarrow k, \quad h(\hat {x}) \Leftrightarrow h(x), \quad h(\hat {k}) \Leftrightarrow h(k), \quad \hat {O}^\dagger \Leftrightarrow O^*, \end{align}

where $h$ is any function and $\hat {O}$ is any operator. For multi-component operators, the Wigner–Weyl transform is applied element-by-element. In particular, if $\hat {O}$ is a matrix operator, one finds that ${\hat {O}}^\dagger \Leftrightarrow O^\dagger$ , where $O$ is the symbol matrix.

The product of two operators maps to the so-called Moyal product, or star product, of their symbols:

(2.32) \begin{eqnarray} \hat {A}\hat {B}\, \Leftrightarrow \, A(x, k) \star B(x, k) \doteq A(x, k)\mathrm{e}^{\mathrm{i}\hat {\mathcal{L}}/2} B(x, k). \end{eqnarray}

Here, $\hat {\mathcal{L}} \doteq \overset {{\scriptscriptstyle \leftarrow }}{\partial }_x \overset {{\scriptscriptstyle \rightarrow }}{\partial }_k - \overset {{\scriptscriptstyle \leftarrow }}{\partial }_k \overset {{\scriptscriptstyle \rightarrow }}{\partial }_x$ and the arrows indicate the directions in which the derivatives act. For example, $A \hat {\mathcal{L}} B$ is just the canonical Poisson bracket:

(2.33) \begin{equation} A \hat {\mathcal{L}} B = \lbrace A, B \rbrace \doteq \frac {\partial A}{\partial x}\frac {\partial B}{\partial k} - \frac {\partial A}{\partial k}\frac {\partial B}{\partial x} \equiv \frac {\partial A}{\partial {\mathsf{z}}^\alpha }\,{\mathsf{J}}^{\alpha \beta }\,\frac {\partial B}{\partial {\mathsf{z}}^\beta }. \end{equation}

These formulae readily yield

(2.34) \begin{eqnarray} h(\hat {x})\hat {k} \Leftrightarrow k h(x) + \frac {\mathrm{i}}{2}\,h'(x), \quad \hat {k} h(\hat {x}) \Leftrightarrow k h(x) - \frac {\mathrm{i}}{2}\,h'(x) \end{eqnarray}

(the prime denotes the derivative with respect to the argument), so

(2.35) \begin{gather} k h(x) \Leftrightarrow \frac{1}{2}(h(\hat {x})\hat {k} + \hat {k}h(\hat {x})) = h(\hat {x})\hat {k} + \frac {1}{2}[\hat {k}, h(\hat {x})] = h(\hat {x})\hat {k} - \frac {\mathrm{i}}{2}\, h'(\hat {x}), \end{gather}

or, in the $x$ -representation,

(2.36) \begin{gather} k h(x) \Leftrightarrow -\mathrm{i} h(x)\partial _x - \mathrm{i} h'(x)/2. \end{gather}

Also, $h(\hat {k}) \mathrm{e}^{\mathrm{i} K \hat {x}} \Leftrightarrow h(k) \star \mathrm{e}^{\mathrm{i} K x} = h(k + K/2) \mathrm{e}^{\mathrm{i} K x}$ , and so on. In particular, the symbol of the commutator $[\hat {A}, \hat {B}]$ can be expressed through the so-called Moyal bracket $\lbrace \! \lbrace A, B \rbrace \! \rbrace$ :

(2.37) \begin{gather} [\hat {A}, \hat {B}] \Leftrightarrow A \star B - B \star A \equiv \mathrm{i} \lbrace \! \lbrace A, B \rbrace \! \rbrace . \end{gather}

The Moyal product and the Moyal bracket are particularly handy when $\partial _{x} \partial _{k} \sim \epsilon \ll 1$ . (Sometimes, $\epsilon$ is equated with the GO parameter $\lambda /L$ , but note that $\star$ is symplectically invariant and the GO parameter is not. One might notice, though, that this $\epsilon$ is similar to the MGO parameter; cf. § 5.4.) Since $\hat {\mathcal{L}} = {\mathcal{O}}(\epsilon )$ , one can express the Moyal product as an asymptotic series in powers of $\epsilon$ :

(2.38) \begin{equation} \star = \hat {1} + \mathrm{i} \hat {\mathcal{L}} /2 - \hat {\mathcal{L}}^2/8 + {\mathcal{O}}(\epsilon ^3). \end{equation}

Accordingly, the Moyal bracket can be approximated with the Poisson bracket:

(2.39) \begin{equation} \lbrace \! \lbrace A, B \rbrace \! \rbrace = \lbrace A, B \rbrace + {\mathcal{O}}(\epsilon ^3) \end{equation}

(assuming that $A$ and $B$ are scalars or commuting matrices). Likewise, the symbol of any operator in the limit $\epsilon \to 0$ can be obtained simply by replacing $\hat {x}$ with $x$ and $\hat {k}$ with $k$ .

2.3. Wigner functions

Of particular interest in the Weyl symbol calculus are ‘density operators’

(2.40)

for various $| {\psi }\rangle$ . (The factor $2\pi$ is sometimes omitted, depending on a convention.) The symbol of such an operator, called the Wigner function,

(2.41) \begin{gather} W^{(\psi )} \doteq {\textrm {symb}} \hat {W}^{(\psi )}, \end{gather}

is a real function and can be expressed as follows:

(2.42a) \begin{align} W^{(\psi )}(x, k) & = \frac {1}{2\pi } \int {\mathrm{d}} s\, \psi _{\hat {x}}(x + s/2)\, \psi _{\hat {x}}^*(x - s/2)\,\mathrm{e}^{-\mathrm{i} k s}\\[-10pt]\nonumber \end{align}

(2.42b) \begin{align} & = \frac {1}{2\pi } \int {\mathrm{d}} r\, \psi _{\hat {k}}(k + r/2)\, \psi _{\hat {k}}^*(k - r/2)\,\mathrm{e}^{\mathrm{i} r x}. \end{align}

(Like in the previous formulae, $2\pi$ here must be replaced with $(2\pi )^n$ if $x$ is $n$ -dimensional.) In particular, note the two special cases:

(2.43)

Any function bilinear in $\psi _{\hat {x}}$ and $\psi ^*_{\hat {x}}$ can be expressed through $W^{(\psi )}$ . For example, for any operators $\hat {L}$ and $\hat {R}$ , one has

(2.44) \begin{align} (\hat {L}\psi _{\hat {x}}(x))(\hat {R}\psi _{\hat {x}}(x))^* = \int {\mathrm{d}} k\,{L}(x, k) \star W^{(\psi )}(x, k) \star {R}^*(x, k), \end{align}

where ${L}$ and ${R}$ are the corresponding symbols. As a corollary, one has

(2.45) \begin{gather} |\psi _{\hat {x}}(x)|^2 = \int {\mathrm{d}} k\,W^{(\psi )}(x, k), \quad |\psi _{\hat {k}}(k)|^2 = \int {\mathrm{d}} x\,W^{(\psi )}(x, k). \end{gather}

If $|\psi _{\hat {x}}(x)|^2$ and $|\psi _{\hat {k}}(k)|^2$ are interpreted as the densities of quanta in the $x$ -space and the $k$ -space, respectively, then one can attribute $W^{(\psi )}$ as a quasiprobability distribution of wave quanta in phase space. The prefix ‘quasi’ is commonly added because $W^{(\psi )}$ can be negative, which is not something that one expects from a probability distribution. That said, a Wigner function averaged over a sufficiently large phase-space volume $\Delta x\,\Delta k \gtrsim 1$ is guaranteed to be non-negative. Such a function can be understood as the spectrum of the two-point correlation function $\langle \psi _{\hat {x}}(x + s/2) \psi _{\hat {x}}^*(x - s/2) \rangle$ over $s$ at given $x$ and represents a local property of the field. For further details, see, for example, Dodin (Reference Dodin2022).

3.1. Definition

Consider a transformation of the base operators to some new base operators:

(3.1) \begin{gather} \hat {{\mathsf{z}}} \equiv \left ( \begin{array}{c} \hat {x}\\[2pt] \hat {k} \end{array} \right ) \mapsto \left ( \begin{array}{c} \hat {q}\\[2pt] \hat {p} \end{array} \right ) \equiv \hat {{\mathsf{y}}}. \end{gather}

Here, $\hat {q}$ can be defined as any Hermitian operator whose eigenvalue space is $\mathbb{R}$ and $\hat {p}$ is defined such that its $q$ -representation is $\hat {p}_{\hat {q}} = -\mathrm{i} \partial _q$ . Accordingly,

(3.2) \begin{gather} [\hat {q}, \hat {p}] = \mathrm{i}, \end{gather}

and thus the corresponding symbols $q(x, k)$ and $p(x, k)$ satisfyFootnote ¹³

(3.3) \begin{gather} \lbrace \! \lbrace {q, p} \rbrace \! \rbrace = 1. \end{gather}

In the GO limit, this corresponds to a variable transformation

(3.4) \begin{gather} {\mathsf{z}} \equiv \left ( \begin{array}{c} x\\[2pt]k \end{array} \right ) \mapsto \left ( \begin{array}{c} q\\[2pt]p \end{array} \right ) \equiv {\mathsf{y}} \end{gather}

with the following Jacobian matrix:

(3.5) \begin{gather} \varXi \doteq \left ( \begin{array}{c@{\quad}c} \partial _x q & \partial _k q\\[2pt] \partial _x p & \partial _k p\\ \end{array} \right )\!. \end{gather}

Since $\varXi {\mathsf{J}} \varXi ^{\intercal } = \lbrace q, p \rbrace {\mathsf{J}}$ and $\lbrace q, p \rbrace \to 1$ by the GO limit of (3.3), this transformation conserves the symplectic form and therefore is canonical. Thus, (3.1) can be interpreted as a canonical transformation for operators.Footnote ¹⁴ Also, as a reminder, canonical transformations conserve the Poisson bracket for any $A$ and $B$ ; i.e.

(3.6) \begin{equation} \lbrace A, B \rbrace \equiv \frac {\partial A}{\partial {\mathsf{z}}^\alpha }\,{\mathsf{J}}^{\alpha \beta }\,\frac {\partial B}{\partial {\mathsf{z}}^\beta } = \frac {\partial A}{\partial {\mathsf{y}}^\alpha }\,{\mathsf{J}}^{\alpha \beta }\,\frac {\partial B}{\partial {\mathsf{y}}^\beta }. \end{equation}

The operator transformation (3.1) induces the following transformation of fields:

(3.7)

or, equivalently,

(3.8) \begin{gather} \psi _{\hat {q}}(q) = \int {\mathrm{d}} x M(q, x)\,\psi _{\hat {x}}(x), \end{gather}

where we have introduced

(3.9)

When $\hat {q} = \hat {k}$ (and thus $\hat {p} = -\hat {x}$ ), the metaplectic transform is simply the Fourier transform, as seen by comparing (3.7) with (2.18). Since the function defines the corresponding operator one can as well express (3.8) asFootnote ¹⁵

(3.10)

We call this operator the MT induced by the operator transformation (3.1).Footnote ¹⁶ Its inverse is the metaplectic operator induced by the inverse transformation $(\hat {q}, \hat {p}) \mapsto (\hat {x}, \hat {k})$ and has the kernel

(3.11)

By (2.10), this means that i.e. $\hat {M}$ is unitary.

A sequence of MTs $\hat {M}_1\hat {M}_0$ corresponding to the operator transformations $\hat {x}_0 \mapsto \hat {x}_1 \mapsto \hat {x}_2$ is an MT $\hat {M}_2$ corresponding to the operator transformation $\hat {x}_0 \mapsto \hat {x}_2$ , because

(3.12) \begin{align} \psi _{\hat {x}_2}(x_2) & = \int {\mathrm{d}} x_1\,M_1(x_2, x_1)\,\psi _{\hat {x}_1}(x_1) \notag \\ & = \int {\mathrm{d}} x_1\,{\mathrm{d}} x_0\,M_1(x_2, x_1)\,M_0(x_1, x_0)\psi _{\hat {x}_0}(x_0) \notag \\ & = \int {\mathrm{d}} x_0\,\mathcal{M}(x_2, x_0)\psi _{\hat {x}_0}(x_0), \end{align}

where

(3.13)

Thus, MTs form a group.

3.2. M-waves

At given $x$ considered as a parameter, the function can be understood as a field on the $q$ -space, specifically, the $q$ -representation of the vector . As such, $M$ can be acted upon by linear operators like any other field. Let us explore the effect of $\hat {x}$ on it. By definition,

(3.14)

This leads to

(3.15) \begin{gather} (\hat {x}_{\hat {q}} M)(q, x) = x M(q, x), \end{gather}

where . Let us symbolically represent $\hat {x}_{\hat {q}}$ through the new base operators $\hat {q}_{\hat {q}} = q$ and $\hat {p}_{\hat {q}} = -\mathrm{i} \partial _q$ as $X(q, -\mathrm{i} \partial _q)$ .Footnote ¹⁷ This leads to the following equation:

(3.16) \begin{gather} X(q, -\mathrm{i} \partial _q) M(q, x) = x M(q, x). \end{gather}

Equation (3.16) is a (pseudo)differential equation for $M$ as a function of $q$ , with $x$ as a parameter. Thus, the solution of (3.16) is defined up only to a function of $x$ and another equation for $M$ is needed to specify this function. To derive such an equation, note that can be understood as the $x$ -representation of . Hence, on the one hand,

(3.17) \begin{gather} (\hat {k}_{\hat {x}} \bar {M})(x, q) = -\mathrm{i} \partial _x \bar {M}(x, q) = - \mathrm{i} \partial _x M^*(q, x). \end{gather}

On the other hand,

(3.18)

Thus, $(\hat {k}_{\hat {q}} M)(q, x) = \mathrm{i} \partial _x M(q, x)$ , or, assuming the notation $\hat {k}_{\hat {q}} = K(q, -\mathrm{i} \partial _q)$ ,

(3.19) \begin{gather} K(q, -\mathrm{i} \partial _q)M(q, x) = \mathrm{i} \partial _x M(q, x). \end{gather}

Equation (3.19) can be interpreted as a Schrödinger equation (SE) for $M$ , with $x$ serving as the time variable. In this sense, $M$ is a wave, which we call an $M$ -wave. Likewise, in retrospect, (3.16) can be interpreted as an SE for a stationary wave with ‘energy’ $x$ .

Together, (3.16) and (3.19) define $M$ up to a constant complex factor. From the normalisation condition, one finds that

(3.20)

This means that, for any given $x$ , $M$ satisfies

(3.21) \begin{gather} \int {\mathrm{d}} q \,{\mathrm{d}} s\, M^*(q, x + s/2) M(q, x - s/2) = 1. \end{gather}

(This is also seen from the definition of $M$ as a representation of a normalised eigenvector.) This defines the normalisation factor up to a constant phase, which can be anything depending how the eigenvectors of $\hat {q}$ are chosen relative to those of $\hat {x}$ . This choice constitutes a gauge freedom of the theory.

Similar equations apply to the inverse-MT kernel $\bar {M}(x, q)$ . That is, assuming the notation $\hat {q}_{\hat {x}} = Q(x, -\mathrm{i} \partial _x)$ and $\hat {p}_{\hat {x}} = P(x, -\mathrm{i} \partial _x)$ , one has

(3.22) \begin{align} Q(x, -\mathrm{i} \partial _x) \bar {M}(x, q) & = q \bar {M}(x, q),\\[-10pt]\nonumber \end{align}

(3.23) \begin{align} P(x, -\mathrm{i} \partial _x) \bar {M}(x, q) & = \mathrm{i} \partial _q \bar {M}(x, q). \end{align}

Equation (3.22) can be interpreted as an SE for a wave $\bar {M}$ , with $q$ serving as the time variable. Because $\bar {M}(x, q) \equiv M^*(x, q)$ , we will also call it an $M$ -wave. Also, (3.22) can be interpreted as an SE for a stationary wave with ‘energy’ $q$ . Like before, one also finds

(3.24)

whence

(3.25) \begin{gather} \int {\mathrm{d}} x\,{\mathrm{d}} s\,\bar {M}^*(x, q + s/2) \bar {M}(x, q - s/2) = 1. \end{gather}

Let us summarise our intermediate results. $M$ -waves are normalised solutions of SEs (3.19) and (3.23), which can also be written in the form (4.1):

(3.26) \begin{gather} (K(q, -\mathrm{i} \partial _q) - \mathrm{i} \partial _x)M = 0, \quad (P(x, -\mathrm{i} \partial _x) - \mathrm{i} \partial _q)\bar {M}=0. \end{gather}

Their solutions can be formally expressed as follows:

(3.27) \begin{gather} M(q, x) = \exp ({-}\mathrm{i} x\hat {k}_{\hat {q}}) M(q, 0), \quad \bar {M}(x, q) = \exp ({-}\mathrm{i} q\hat {p}_{\hat {x}}) \bar {M}(x, 0). \end{gather}

It is generally impossible to find explicit analytic expressions for the propagators $\exp ({-}\mathrm{i} x\hat {k}_{\hat {q}})$ and $\exp ({-}\mathrm{i} q\hat {p}_{\hat {x}})$ for a given variable transformation. However, it is possible to do so for some important special cases and also asymptotically. This will be discussed later, in §§ 3.5 and 4.3.

3.3. Symplectic pseudo-measure

As a function of two variables, an $M$ -wave can be considered as a field on the two-dimensional ‘timespace’ $(x, q)$ . The variable $x$ serves as the time variable in the Schrödinger equation (3.19), so $\mathrm{i} \partial _x$ serves as the frequency operator and $-k$ serves as the frequency variable. Hence, MGO naturally involves Wigner functions of $M$ in the form

(3.28) \begin{gather} W^{(M)}_{\mathsf{z}}(x, q, -k, p) \equiv \mu ({\mathsf{y}}, {\mathsf{z}})/2\pi \end{gather}

(the utility of the added coefficient $2\pi$ will become clear shortly), i.e.

(3.29) \begin{gather} \mu ({\mathsf{y}}, {\mathsf{z}}) = \frac {1}{2\pi } \int {\mathrm{d}} s\,{\mathrm{d}} s'\, M(q + s'\!/2, x + s/2) \, M^*(q - s'\!/2, x - s/2)\,\mathrm{e}^{-\mathrm{i} p s' + \mathrm{i} k s}. \end{gather}

For the inverse transformation, the corresponding function is

(3.30) \begin{align} \bar {\mu }({\mathsf{z}}, {\mathsf{y}}) & = \frac {1}{2\pi } \int {\mathrm{d}} s\,{\mathrm{d}} s'\, \bar {M}(x + s/2, q + s'\!/2) \, \bar {M}^*(x - s/2, q - s'\!/2)\,\mathrm{e}^{-\mathrm{i} k s + \mathrm{i} p s'} \notag \\[2pt] & = \frac {1}{2\pi } \int {\mathrm{d}} s\,{\mathrm{d}} s'\, M^*(q + s'\!/2, x + s/2) \, M(q - s'\!/2, x - s/2)\,\mathrm{e}^{\mathrm{i} p s' - \mathrm{i} k s} \notag \\ & = \mu ^*({\mathsf{y}}, {\mathsf{z}}).\vphantom {\int }\\[-20pt]\nonumber \end{align}

Like any Wigner function, $\mu$ is real, so one can also express (3.30) simply as

(3.31) \begin{gather} \bar {\mu }({\mathsf{z}}, {\mathsf{y}}) = \mu ({\mathsf{y}}, {\mathsf{z}}). \end{gather}

Also note that

(3.32) \begin{align} \int \mu ({\mathsf{y}}, {\mathsf{z}})\,{\mathrm{d}} {\mathsf{z}} & = \int {\mathrm{d}} x\,{\mathrm{d}} s\,{\mathrm{d}} s'\, M(q + s'\!/2, x + s/2) \, M^*(q - s'\!/2, x - s/2)\,\delta (s)\,\mathrm{e}^{-\mathrm{i} p s'} \notag \\ & = \int {\mathrm{d}} s'\,\mathrm{e}^{-\mathrm{i} p s'}\int {\mathrm{d}} x\, M(q + s'\!/2, x) \, M^*(q - s'\!/2, x) \notag \\ & = \int {\mathrm{d}} s'\,\mathrm{e}^{-\mathrm{i} p s'}\,\delta (s') \notag \\ & = 1 \end{align}

(here we used (3.24)) and, similarly, $\int \mu ({\mathsf{y}}, {\mathsf{z}})\,{\mathrm{d}}{\mathsf{y}} = 1$ too. In summary then,

(3.33) \begin{gather} \int \mu ({\mathsf{y}}, {\mathsf{z}})\,{\mathrm{d}} {\mathsf{z}} = \int \mu ({\mathsf{y}}, {\mathsf{ z}})\,{\mathrm{d}} {\mathsf{y}} = 1. \end{gather}

As to be seen later (§ 3.4), these serve as pseudo-measures on the ${\mathsf{z}}$ -space and the ${\mathsf{y}}$ -space, respectively. (The prefix ‘pseudo’ is added because $\mu$ can be negative.) From here, the arguments of $\mu$ will be occasionally omitted where they are obvious from the context.

3.4. Transformation of operators

A metaplectic transform induces a transformation of the operator representation $\hat {O}_{\hat {x}} \mapsto \hat {O}_{\hat {q}}$ , which is found as follows. Note that

(3.34)

Using (2.9), one can also represent this as follows:

(3.35)

whence

(3.36) \begin{gather} \hat {O}_{\hat {q}} = \hat {M}\hat {O}_{\hat {x}} \hat {M}^\dagger . \end{gather}

Let us use the following formula for the Weyl symbol of a given operator (cf. (2.28a )):

(3.37)

The subindex in the symbol’s notation refers to the fact that this symbol is defined on the phase space ${\mathsf{y}} \doteq (q, p)^{\intercal }$ . We will call it the ${\mathsf{y}}$ -symbol of $\hat {O}$ , and the symbol defined on the phase space ${\mathsf{z}} \doteq (x, k)^{\intercal }$ will be called the ${\mathsf{z}}$ -symbol of $\hat {O}$ . Recall that

(3.38)

Then,

(3.39)

Expressing the brackets through $M$ , one obtains the following expression for the ${\mathsf{ y}}$ -symbol in terms of the ${\mathsf{z}}$ -symbol:

(3.40) \begin{align} O_{\mathsf{y}}(q, p) & = \frac {1}{2\pi } \int {\mathrm{d}} x\,{\mathrm{d}} k\,{\mathrm{d}} s\,{\mathrm{d}} s'\, M(q + s'\!/2, x + s/2) \, M^*(q - s'\!/2, x - s/2) \notag \\ &\quad \times O_{\mathsf{z}}(x, k)\,\mathrm{e}^{-\mathrm{i} p s' + \mathrm{i} k s}, \end{align}

and the inverse transformation can be written similarly. Equivalently, this result can be expressed in the following compact form:

(3.41) \begin{gather} O_{\mathsf{y}}({\mathsf{y}}) = \int O_{\mathsf{z}}({\mathsf{z}})\,\mu ({\mathsf{y}}, {\mathsf{ z}})\,{\mathrm{d}} {\mathsf{z}}, \quad O_{\mathsf{z}}({\mathsf{z}}) = \int O_{\mathsf{y}}({\mathsf{y}})\,\mu ({\mathsf{y}}, {\mathsf{ z}})\,{\mathrm{d}} {\mathsf{y}}. \end{gather}

(If $O_{\mathsf{z}} \equiv 1$ , then (3.40) readily yields $O_{\mathsf{y}} \equiv 1$ , and vice versa, by (3.33).) In particular, these allow one to recalculate the Wigner functions from one representation to another.

3.5. Examples

In this section, we present several examples of MTs that are used as building blocks of the theory presented in later sections.

3.5.1. Shift

For example, let us consider a variable transformation that is a phase-space shift by a constant vector $\varDelta \doteq (\varDelta _q, \varDelta _p)^{\intercal }$ :

(3.42) \begin{gather} \hat {q} = \hat {x} + \varDelta _q \hat {1}, \quad \hat {p} = \hat {k} + \varDelta _p \hat {1}. \end{gather}

This corresponds to

(3.43) \begin{gather} X(q, -\mathrm{i} \partial _q) = q - \varDelta _q, \quad K(q, -\mathrm{i} \partial _q) = -\mathrm{i} \partial _q - \varDelta _p. \end{gather}

Then, (3.16) requires that $M(q, x) = C(q)\delta (q - x - \varDelta _q)$ , where $C$ is some function. From (3.19), one finds that $({-}\mathrm{i} \partial _q - \mathrm{i} \partial _x - \varDelta _p)M(q, x) = 0$ , whence $\partial _q C = \mathrm{i} \varDelta _p C$ . With (3.21) for the normalisation, this gives

(3.44) \begin{gather} M(q, p) = \mathrm{e}^{\mathrm{i} \varDelta _p q}\,\delta (q - x - \varDelta _q) \end{gather}

(up to an arbitrary constant phase factor, which we choose equal to one) and, therefore,

(3.45) \begin{gather} \psi _{\hat {q}}(q) = \mathrm{e}^{\mathrm{i} \varDelta _p q}\psi _{\hat {x}}(q - \varDelta _q). \end{gather}

Notice that even at $\varDelta _q = 0$ , when $q$ is the same as $x$ , the function $\psi _{\hat {q}}$ is not the same as $\psi _{\hat {x}}$ due to non-zero $\varDelta _p$ . Also,

(3.46) \begin{align} \mu & = \frac {1}{2\pi } \int {\mathrm{d}} s\,{\mathrm{d}} s'\, \mathrm{e}^{\mathrm{i} \varDelta _p s' - \mathrm{i} p s' + \mathrm{i} k s}\,\delta (q + s'\!/2 - x - s/2 - \varDelta _q)\,\delta (q - s'\!/2 - x + s/2 - \varDelta _q) \notag \\[3pt] & = \frac {1}{2\pi } \int {\mathrm{d}} s\,{\mathrm{d}} s'\, \mathrm{e}^{\mathrm{i} (\varDelta _p - p) s' + \mathrm{i} k s}\,\delta (q - x - \varDelta _q)\,\delta (s - s') \notag \\[3pt] & = \frac {1}{2\pi } \int {\mathrm{d}} s\, \mathrm{e}^{\mathrm{i} (k + \varDelta _p - p) s}\,\delta (q - x - \varDelta _q) \notag \\[3pt] & = \delta (x - q + \varDelta _q)\,\delta (k - p + \varDelta _p). \end{align}

Then, by (3.41), one obtains

(3.47) \begin{gather} O_{\mathsf{y}}(q, p) = O_{\mathsf{z}}(q - \varDelta _q, p - \varDelta _p) = O_{\mathsf{z}}(x(q, p), k(q, p)). \end{gather}

Simply put, this means that $O_{\mathsf{y}}({\mathsf{y}}) = O_{\mathsf{z}}({\mathsf{z}})$ ; i.e. the symbol of an operator is invariant with respect to phase-space shifts.

3.5.2. Rescaling

Let us consider a rescaling canonical transformation:

(3.48) \begin{gather} \hat {q} = \hat {x}/\alpha , \quad \hat {p} = \alpha \hat {k}, \end{gather}

where $\alpha$ is a non-zero real constant. (Both equations contain the same factor $\alpha$ to keep the transformation canonical.) Clearly, $M(q, x) = C\delta (q - x/\alpha )$ , where $C$ is a constant. From (3.21), one finds that $C = |\alpha |^{-1/2}$ up to an arbitrary phase, so

(3.49) \begin{gather} M(q, x) = |\alpha |^{-1/2}\,\delta (q - x/\alpha ), \quad \psi _{\hat {q}}(q) = |\alpha |^{1/2} \psi _{\hat {x}}(\alpha q). \end{gather}

3.5.3. Linear symplectic transformation

Let us also consider a linear symplectic transformation (LST), that is,

(3.50) \begin{gather} \left ( \begin{array}{c} \hat {q} \\[3pt] \hat {p} \end{array} \right ) = \underbrace { \left ( \begin{array}{c@{\quad}c} A & B\\[3pt] C & D \end{array} \right ) }_{{\mathsf{S}}} \left ( \begin{array}{c} \hat {x} \\[3pt] \hat {k} \end{array} \right )\!, \quad \end{gather}

where the (constant) coefficients satisfy

(3.51) \begin{gather} AD - BC = 1, \end{gather}

i.e. $\det {{\mathsf{S}}} = 1$ .Footnote ¹⁸ This subsumes relabelling coordinates and momenta; specifically, $A = D = 0$ and $B = - C = \pm 1$ corresponds to the transformation $(\hat {q}, \hat {p}) = (\pm \hat {k}, \mp \hat {x})$ . (Note that $\hat {k}$ and $\hat {x}$ necessarily enter the latter formula with opposite signs to ensure symplecticity.) The rescaling transformation considered in § 3.5.2 is also subsumed as a limit ( $B \to 0$ , $C \to 0$ and $A \to 1/D$ ).

Equivalently, (3.50) can be represented in the following inverted form:

(3.52) \begin{gather} \left ( \begin{array}{c} \hat {x} \\[3pt] \hat {k} \end{array} \right ) = \underbrace {\left ( \begin{array}{c@{\quad}c} \bar {A} & \bar {B}\\[3pt] \bar {C} & \bar {D} \end{array} \right )}_{\bar {{{\mathsf{S}}}}} \left ( \begin{array}{c} \hat {q} \\[3pt] \hat {p} \end{array} \right )\!, \quad \end{gather}

where $\bar {{{\mathsf{S}}}} = {{\mathsf{S}}}^{-1}$ is also symplectic; specifically,

(3.53) \begin{gather} \left ( \begin{array}{c@{\quad}c} \bar {A} & \bar {B}\\[3pt] \bar {C} & \bar {D} \end{array} \right ) = \left ( \begin{array}{c@{\quad}c} D & -B\\[3pt] -C & A \end{array} \right )\!, \end{gather}

so $\det \bar {{{\mathsf{S}}}} = \det {{\mathsf{S}}} = 1$ . Then, (3.16) and (3.19) yields

(3.54) \begin{gather} (\bar {A}q - \mathrm{i} \bar {B} \partial _q - x) M(q, x) = 0, \quad (\bar {C}q - \mathrm{i} \bar {D} \partial _q - \mathrm{i} \partial _x) M(q, x) = 0. \end{gather}

With the normalisation (3.21) taken into account, a straightforward calculation yields the well-known result (Littlejohn Reference Littlejohn1986):

(3.55) \begin{gather} M(q, x) = (2\pi \mathrm{i} \bar {B})^{-1/2}\, \exp \left ({-} \frac {\mathrm{i}(\bar {A}q^2 - 2xq + \bar {D}x^2)}{2\bar {B}}\right )\!. \end{gather}

As usual, an arbitrary constant phase can be added to the overall expression, depending on a convention. The phase of the square root as a function of $\bar {B}$ must be chosen such that MTs remain a group. Loosely speaking, one can require that a phase-space rotationFootnote ¹⁹ by $\pm 2\pi$ causes the phase of $\bar {B}$ to change by $\pm 2\pi$ too; then, the phase of $\bar {B}^{1/2}$ changes by $\pm \pi$ . (For a detailed discussion, see Littlejohn (Reference Littlejohn1986).) Then, the MT corresponding to a complete phase-space rotation in either direction is not an identity operator $\hat {1}$ but $-\hat {1}$ .

The corresponding transformation of symbols is calculated as follows. First, note that

(3.56) \begin{align} &M(q + s'\!/2, x + s/2)\, M^*(q - s'\!/2, x - s/2) \,\mathrm{e}^{- \mathrm{i} p s' + \mathrm{i} k s} \notag \\ &\quad = \frac {1}{2\pi |\bar {B}|}\, \exp \left (\frac {\mathrm{i} s(q - \bar {D}x + \bar {B}k) + \mathrm{i} s'(x - \bar {A} q - \bar {B}\!p)}{\bar {B}}\right )\!, \end{align}

and

(3.57) \begin{align} \mu & = \int \frac {{\mathrm{d}} s'}{2\pi |\bar {B}|}\,\exp \left (\frac {\mathrm{i} s'(x - \bar {A} q - \bar {B}\!p)}{\bar {B}}\right ) \int \frac {{\mathrm{d}} s}{2\pi }\, \exp \left (\mathrm{i} s\left (\frac {q - \bar {D}x}{\bar {B}} + k\right )\right ) \notag \\[2pt] & = \delta (x - \bar {A} q - \bar {B}\!p)\, \delta \!\left (k - \frac {\bar {D} x - q}{\bar {B}}\right ) \notag \\[2pt] & = \delta (x - \bar {A} q - \bar {B}\!p)\, \delta \!\left (k - \frac {\bar {D}(\bar {A} q + \bar {B}\!p) - q}{\bar {B}}\right ) \notag \\[2pt] & = \delta (x - \bar {A} q - \bar {B}\!p)\, \delta (k - \bar {C}q - \bar {D}\!p), \end{align}

where we used (3.51). Then, by (3.41), one obtains

(3.58) \begin{gather} O_{\mathsf{y}}(q, p) = O_{\mathsf{z}}(\bar {A}q + \bar {B}\!p, \bar {C}q + \bar {D}x) = O_{\mathsf{z}}(x(q, p), k(q, p)). \end{gather}

This shows that the symbols of operators are invariant with respect to LSTs. In particular, Wigner functions are LST-invariant. Also, for $\hat {O} = \hat {{\mathsf{y}}}$ , (3.58) gives ${\mathsf{z}}_{\mathsf{y}} = \bar {{{\mathsf{ S}}}} {\mathsf{y}}$ and, similarly, ${\mathsf{y}}_{\mathsf{z}} = {{\mathsf{S}}} {\mathsf{z}}$ ; i.e. the symbols of the base operators transform at LSTs exactly as the base operators themselves.

3.5.4. Eikonal transform

Let us also consider the ‘eikonal transform’ (this term will become clear shortly):

(3.59) \begin{gather} \hat {q} = \hat {x}, \quad \hat {p} = \hat {k} - \mathcal{K}(\hat {x}) = \hat {k} - \mathcal{K}(\hat {q}), \end{gather}

which satisfies the commutation relation (3.2). Like in § 3.5.1, one finds that

(3.60) \begin{gather} M(q, p) = \mathrm{e}^{-\mathrm{i} \theta (q)}\,\delta (q - x), \end{gather}

where $\theta \doteq \int \mathcal{K}(q)\,{\mathrm{d}} q = \int \mathcal{K}(x)\,{\mathrm{d}} x$ . Then, $\psi _{\hat {q}}$ is given by

(3.61) \begin{gather} \psi _{\hat {q}}(q) = \psi _{\hat {q}}(x) = \mathrm{e}^{-\mathrm{i} \theta (x)}\psi _{\hat {x}}(x), \end{gather}

and

(3.62) \begin{align} \mu & = \frac {1}{2\pi } \int {\mathrm{d}} s\,{\mathrm{d}} s'\, \mathrm{e}^{-\mathrm{i} \theta (q + s'\!/2) + \mathrm{i} \theta (q - s'\!/2) - \mathrm{i} p s' + \mathrm{i} k s}\, \delta (q + s'\!/2 - x - s/2)\delta (q - s'\!/2 - x + s/2) \notag \\[3pt] & = \frac {1}{2\pi } \int {\mathrm{d}} s\,{\mathrm{d}} s'\, \mathrm{e}^{-\mathrm{i} \theta (q + s'\!/2) + \mathrm{i} \theta (q - s'\!/2) - \mathrm{i} p s' + \mathrm{i} k s}\, \delta (s' - s)\delta (q - x) \notag \\[3pt] & = \frac {1}{2\pi } \int {\mathrm{d}} s\, \mathrm{e}^{-\mathrm{i} \theta (q + s/2) + \mathrm{i} \theta (q - s/2) + \mathrm{i} (k - p) s}\, \delta (q - x). \end{align}

Then, by (3.41), the operator symbols are transformed as follows:

(3.63) \begin{gather} O_{\mathsf{y}}(q, p) = \frac {1}{2\pi } \int {\mathrm{d}} k\,{\mathrm{d}} s\, O_{\mathsf{z}}(q, k)\,\mathrm{e}^{\mathrm{i} (k - p)s -\mathrm{i} \theta (q + s/2) + \mathrm{i} \theta (q - s/2)}. \end{gather}

Using Taylor expansion, one has

(3.64) \begin{gather} \theta \left (q + \frac {s}{2}\right ) - \theta \left (q - \frac {s}{2}\right ) = \mathcal{K}(q)\, s + \frac {1}{24}\frac {\partial ^2 \mathcal{K}(q)}{\partial q^2}\, s^3 + \ldots , \end{gather}

where we used $\partial _q \theta (q) = \mathcal{K}(q)$ . Thus,

(3.65) \begin{align} O_{\mathsf{y}}(q, p) & = \frac {1}{2\pi } \int {\mathrm{d}} k\,{\mathrm{d}} s\, O_{\mathsf{z}}(q, k)\, \left (1 - \frac {\mathrm{i}}{24}\frac {\partial ^2 \mathcal{K}(q)}{\partial q^2}\,s^3 + \ldots \right )\mathrm{e}^{\mathrm{i} (k - p - \mathcal{K}(q))s} \notag \\[2pt] & = \frac {1}{2\pi } \left (1 - \frac {1}{24}\frac {\partial ^2 \mathcal{K}(q)}{\partial q^2}\,\frac {\partial ^3}{\partial p^3} + \ldots \right )\int {\mathrm{d}} k\,{\mathrm{d}} s\, O_{\mathsf{z}}(q, k)\, \mathrm{e}^{\mathrm{i} (k - p - \mathcal{K}(q))s} \notag \\[2pt] & = \frac {1}{2\pi } \int {\mathrm{d}} k\,{\mathrm{d}} s\, O_{\mathsf{z}}(q, k)\, \mathrm{e}^{\mathrm{i} (k - p - \mathcal{K}(q))s} + {\mathcal{O}}(\epsilon ^2) \notag \\ & = O_{\mathsf{z}}(q, p + \mathcal{K}(q)) + {\mathcal{O}}(\epsilon ^2). \vphantom {\int } \end{align}

In particular, if $\theta (q)$ is the eikonal and $\psi _{\hat {q}}$ is the wave envelope, then $\mathcal{K}$ serves as the reference wavevector and the symbol transformation reproduces the one derived, for example, by Dodin et al. (Reference Dodin, Ruiz, Yanagihara, Zhou and Kubo2019). The importance of (3.65) in this case is in that, within GO, where ${\mathcal{O}}(\epsilon ^2)$ is negligible, one may use

(3.66) \begin{gather} O_{\mathsf{y}}(q, p) \approx O_{\mathsf{z}}(q, p + \mathcal{K}(q)) = O_{\mathsf{z}}(x(q, p), k(q, p)). \end{gather}

4.1. Basic equations

The Weyl symbol calculus and MTs provide a natural framework for formulating GO equations (Tracy et al. Reference Tracy, Brizard, Richardson and Kaufman2014). Here, we restate the conventional GO (McDonald Reference McDonald1988; Dodin et al. Reference Dodin, Ruiz, Yanagihara, Zhou and Kubo2019), which will be an important reference point for our further discussion. Let us start with a generic linear wave equationFootnote ²⁰

(4.1) \begin{gather} \hat {D}_{\hat {x}}\psi _{\hat {x}} = 0. \end{gather}

Here, $\psi _{\hat {x}}$ is a field on some $x$ -space (generally, timespace), which we assume one-dimensional for simplicity, and $\hat {D}_{\hat {x}}$ is some dispersion operator, which can be an integral operator or a differential operator as a special case. Suppose that $\psi _{\hat {x}}$ is quasimonochromatic, i.e. representable as

(4.2) \begin{gather} \psi _{\hat {x}}(x) = \mathrm{e}^{\mathrm{i}\theta (x)}\varPsi (x), \end{gather}

where the (real) phase $\theta$ is fast compared with its gradientFootnote ²¹ $k(x) \doteq \partial _x \theta (x)$ and to the envelope $\varPsi$ . The evolution of the latter can be described by the following equation:

(4.3) \begin{gather} \hat {\mathcal{D}}_{\hat {x}}\varPsi = 0, \quad \hat {\mathcal{D}} \doteq \mathrm{e}^{-\mathrm{i}\theta (\hat {x})}\hat {D}\mathrm{e}^{\mathrm{i}\theta (\hat {x})}. \end{gather}

Our goal is to approximate the operator $\hat {\mathcal{D}}$ with a low-order differential operator.

As we have just discussed (§ 3.5.4), the mapping from $\psi _{\hat {x}}$ to $\varPsi$ can be considered as an MT that corresponds to the variable transformation (3.59), so

(4.4) \begin{gather} \hat {\mathcal{D}}_{\hat {x}} = \hat {D}_{\hat {q}}, \quad \hat {x} = \hat {q}, \quad \hat {p} = \hat {k} - k(\hat {x}), \end{gather}

i.e. the envelope-evolution operator is just the $q$ -representation of the original ‘full-wave’ dispersion operator $\hat {D}$ . Then, its symbol can be calculated using (3.66). Since $\varPsi$ that $\hat {\mathcal{D}}_{\hat {x}}$ acts upon is smooth, only small values of $p$ matter in this case, so one can further approximate (3.66) with its first-order Taylor expansion:Footnote ²²

(4.5) \begin{gather} D_{\mathsf{y}}(q, p) \approx D_{\mathsf{z}}(q, k(q)) + p V(q), \quad V(q) \doteq (\partial _p D_{\mathsf{z}}(q, p))_{p = k(q)}. \end{gather}

Loosely, the expansion parameter here is the ratio of the characteristic wavenumber $p$ of the envelope and $k(q)$ . Denoting the former as the inverse inhomogeneity scale, $1/L$ , and introducing the wavelength $\lambda \sim 1/k(q)$ , this small parameter can be roughly estimated as

(4.6) \begin{gather} \epsilon \sim \lambda /L. \end{gather}

Using the analogue of (2.36) for base operators $\hat {q}$ and $\hat {p}$ , one immediately obtains

(4.7) \begin{gather} \hat {D}_{\hat {q}} \approx D_{\mathsf{z}}(q, k(q)) - \mathrm{i} V(q) \partial _q - \frac {\mathrm{i}}{2}\,V'(q), \end{gather}

where the prime denotes the derivative with respect to the argument, so $V'(q) = (\partial ^2_{qp} D_{\mathsf{z}}(q, p) + p'(q)\,\partial ^{2}_{pp} D_{\mathsf{z}}(q, p))_{p = p(q)}$ . Using (4.4), one can then rewrite (4.3) as follows:

(4.8) \begin{gather} (H_{\mathsf{z}}(x) + \mathrm{i} \varGamma _{\mathsf{z}}(x))\varPsi - \mathrm{i} V(x) \partial _x \varPsi - \frac {\mathrm{i}}{2}\,V'(x)\varPsi \approx 0, \end{gather}

where we introduced $H_{\mathsf{z}}(x) \doteq {\textrm {re}} D_{\mathsf{z}}(x, k(x))$ and $\varGamma _{\mathsf{z}}(x) \doteq {\textrm {im}} D_{\mathsf{z}}(x, k(x))$ . (For vector waves, these would be the Hermitian and anti-Hermitian parts of $D_{\mathsf{z}}(q, k(q))$ , correspondingly.) It is common to assume that ${\textrm {re}} D_{\mathsf{z}}$ is an order-one function, while ${\textrm {im}} D_{\mathsf{z}}$ is order- $\epsilon$ . Also, remember that $\partial _x \varPsi$ and $V'$ are order- $\epsilon$ as well, so one can redefine $V$ as

(4.9) \begin{gather} V(x) \doteq (\partial _k H_{\mathsf{z}}(x, k))_{k = k(x)} \end{gather}

without loss of accuracy. Then, all coefficients in (4.8) are real, except $\mathrm{i}$ per se.

Note that the function $\theta$ has been unspecified so far. One can choose it such thatFootnote ²³

(4.10) \begin{gather} 0 = H_{\mathsf{z}}({\mathsf{z}}) \equiv H_{\mathsf{z}}(x, \partial _x \theta ) \end{gather}

(here, ${\mathsf{z}} \equiv (x, k)^{\intercal }$ , as usual), which can be considered as a local dispersion relation or a Hamilton–Jacobi equation for $\theta$ . For given initial conditions, this determines the ‘reference ray’ trajectory $k(x)$ , which can also be described by Hamilton’s ray equations (Appendix A)

(4.11) \begin{gather} \dot {{\mathsf{z}}} = {\mathsf{J}} \partial _{\mathsf{z}} H_{\mathsf{z}} \equiv \lbrace {\mathsf{z}}, H_{\mathsf{z}} \rbrace . \end{gather}

The dot in $\dot {{\mathsf{z}}}$ denotes a derivative with respect to the ray time (which can be different from the physical time $t$ ; see Appendix A) and $H_{\mathsf{z}}$ serves as the ray Hamiltonian. Accordingly, we henceforth address $\hat {H}$ , i.e. the Hermitian part of $\hat {D}$ , as the wave Hamiltonian and $\hat {\varGamma }$ , i.e. the anti-Hermitian part of $\hat {D}$ , as the dissipation operator.

The reference ray determines the coefficients in (4.8), which becomes

(4.12) \begin{gather} V(x) \partial _x \varPsi + \frac {1}{2}\,V'(x)\varPsi = \varGamma _{\mathsf{z}}(x) \varPsi , \quad \partial _x (V(x) |\varPsi |^2) = 2 \varGamma _{\mathsf{z}}(x) |\varPsi |^2. \end{gather}

(As a reminder, $V'(x)$ in these equations is generally not just $\partial ^2_{xk} H_{\mathsf{ z}}$ , but also includes the term $k'(x)\,\partial ^2_{kk} H_{\mathsf{z}}$ .) Since the coefficients in (4.12) are real, the phase of $\varPsi$ is conserved and one also readily finds $|\varPsi |$ :

(4.13) \begin{gather} |\varPsi |^2 = \frac {C}{V(x)} \exp \left (\int \frac {2\varGamma _{\mathsf{z}}(x)}{V(x)}\,{\mathrm{d}} x\right )\!, \end{gather}

where the constant $C$ is determined by the initial conditions. A straightforward generalisation to the multidimensional space $\boldsymbol{x}$ gives

(4.14) \begin{gather} {\nabla } \boldsymbol{\cdot} (\boldsymbol{V}\, |\varPsi |^2) = 2 \varGamma _{\mathsf{z}} |\varPsi |^2, \quad \boldsymbol{V} \doteq \partial _{\boldsymbol{k}} H_{\mathsf{z}}(\boldsymbol{x}, \boldsymbol{k})_{\boldsymbol{k} = \boldsymbol{k}(\boldsymbol{x})}, \end{gather}

with $\boldsymbol{k}(\boldsymbol{x}) \doteq \boldsymbol{\nabla }\theta$ . In particular, if $\hat {D}$ is Hermitian, (4.12) becomes a conservation law:

(4.15) \begin{gather} {\nabla } \boldsymbol{\cdot} (\boldsymbol{V}\,|\varPsi |^2) = 0. \end{gather}

For stationary waves, this can be understood as energy conservation. More generally, when $\boldsymbol{x}$ includes time, (4.15) represents conservation of the wave action, which is conserved even when the energy is not. For details, see, for example, Dodin (Reference Dodin2022, § 7). A generalisation to vector waves is discussed by Dodin et al. (Reference Dodin, Ruiz, Yanagihara, Zhou and Kubo2019) and, in the broader context of MGO, in § 7.

Although often convenient, this model significantly relies on the smallness of the GO parameter (4.6). Near cutoffs, where the local $\lambda$ is large (infinite) and the group velocity turns to zero, (4.13) predicts $|\varPsi |^2 \sim 1/V(x) \to \infty$ , which indicates limitations of said approximation. In this case, it is useful to promote (4.1) to an abstract vector equation,

(4.16)

and seek an alternative representation of (4.16) in which the corresponding wavelength changes with the new coordinate slowly or not at all. We will discuss this in § 5, after we have introduced some more analytical tools in the next sections.

4.2. Geometrical optics of M-waves: special case

It is instructive to apply the previous approach to $M$ -waves in particular, which can be studied using GO methods like any other waves. In this section, we consider an important special case where the new variables $(q, p)$ are linked to $H_{\mathsf{z}}$ as follows. (Related transformations will be considered in §§ 5.3 and 5.4.)

Like in § 4.1, suppose that the reference-ray trajectory is given by some $k = k(x)$ . Then, near the reference ray, one hasFootnote ²⁴

(4.17) \begin{gather} H_{\mathsf{z}}({\mathsf{z}}) \approx V(x)(k - k(x)), \end{gather}

where $V(x) \doteq (\partial _k H_{\mathsf{z}})_{k = k(x)}$ can be interpreted as the group velocity on the reference ray. Let us choose the new variables $(q, p)$ such that $p = H_{\mathsf{ z}}(x, k)$ , whence $k$ can be readily expressed as a function of $(x, p)$ :

(4.18) \begin{gather} k(x, p) = k(x) + p/V(x). \end{gather}

To find the corresponding canonical coordinate $q$ , we look for a generating function in the type-2 form, $F = F(x, p)$ , such that (Goldstein, Poole & Safko Reference Goldstein, Poole and Safko2011)

(4.19) \begin{gather} k(x, p) = \partial _x F(x, p), \quad q(x, p) = \partial _p F(x, p). \end{gather}

The former gives

(4.20) \begin{gather} F(x, p) = \int ^{x}_0 k(\tilde {x}, p)\,{\mathrm{d}} \tilde {x}, \end{gather}

where we set the integration constant to zero by choice. (The tilde is used to distinguish the dummy integration variable $\tilde {x}$ from the actual canonical coordinate $x$ .) Then,

(4.21) \begin{gather} q(x, p) = \int ^{x}_0 \frac {\partial k(\tilde {x}, p)}{\partial p}\,{\mathrm{d}} \tilde {x} = \int ^{x}_0 \frac {{\mathrm{d}} \tilde {x}}{V(\tilde {x})} \equiv q(x). \end{gather}

Now, we can derive the corresponding approximation for $M$ -waves. One way to do this is to perform an eikonal MT on $M$ and directly apply the results of § 4.1. However, since we have already introduced the expansion (4.17), one might as well notice that the Weyl transform of (4.17) readily yields the following approximation of $\hat {H}_{\hat {x}}$ :

(4.22) \begin{gather} \hat {H}_{\hat {x}} \approx -\mathrm{i} V(x) \partial _x - \frac {\mathrm{i}}{2}\,V'(x) - V(x)k(x). \end{gather}

Then, (3.23) leads to the following equations for $\bar {M} \equiv \bar {M}(x, q)$ :

(4.23) \begin{gather} \mathrm{i} \partial _{q} \bar {M} = \left ({-}\mathrm{i} V(x) \partial _x - \frac {\mathrm{i}}{2}\,V'(x) - V(x)k(x)\right )\bar {M}. \end{gather}

Also, (3.22) gives $(q - q(x))\bar {M}(x, q) = 0$ , whence $\bar {M}(x, q)\, \propto \, \delta (q - q(x))$ . Thus,

(4.24) \begin{gather} \bar {M}(x, q) = \frac {\mathrm{e}^{\mathrm{i} \vartheta (x)}}{\sqrt {V(x)}}\,\delta (q - q(x)), \quad \vartheta (x) \doteq \int ^{x}_0 k(\tilde {x})\,{\mathrm{d}} \tilde {x}, \end{gather}

where the normalisation is chosen such that (3.24) is satisfied:

(4.25) \begin{align} \int {\mathrm{d}} x \,\bar {M}^*(x, q') \bar {M}(x, q'') & = \int \frac {{\mathrm{d}} x}{|V(x)|} \,\delta (q' - q(x))\,\delta (q'' - q(x)) \notag \\[2pt] & = \delta (q' - q'') \int \frac {{\mathrm{d}} x}{|V(x)|} \,\delta (q' - q(x)) \notag \\[2pt] & = \delta (q' - q'') \int {\mathrm{d}} q\,\delta (q' - q) \notag \\[2pt] & = \delta (q' - q''). \end{align}

Note that (4.24) readily provides a general GO solution to (4.1). Indeed, according to (4.16), a wave $| {\psi }\rangle$ is an eigenstate of $\hat {H} = \hat {p}$ corresponding to the eigenvalue $p = 0$ . Then, the $q$ -representation of this field is

(4.26)

where $a$ is a constant amplitude and we have used (2.17). Accordingly,

(4.27) \begin{gather} \psi _{\hat {x}}(x) = \int {\mathrm{d}} q\,\bar {M}(x, q)\,\psi _{\hat {q}}(q) = \psi _{\hat {q}} \int {\mathrm{d}} q\,\bar {M}(x, q) = \psi _{\hat {q}}\,\frac {\mathrm{e}^{\mathrm{i} \vartheta (x)}}{\sqrt {V(x)}}. \end{gather}

This coincides with the well-known WKB solution for eikonal waves, with $\psi _{\hat {q}}$ serving as a constant factor determined by the initial conditions.

Of course, the same result can as well be obtained by directly applying GO to $\psi _{\hat {x}}$ . However, the advantage of using $M$ -waves is that doing so separates the problem of calculating the evolution of $| {\psi }\rangle$ from calculating its representation. Errors that may emerge from the inaccuracy of the mapping $\psi _{\hat {q}} \mapsto \psi _{\hat {x}}$ do not accumulate (except, possibly, in the phase). For example, even though (4.27) fails when $V$ turns to zero, its validity gets reinstated as soon as $V$ becomes large enough. In other words, one does not need to integrate the amplitude equation through singularity. One needs to solve continuously only for $\psi _{\hat {q}}$ (which is constant here but can be more complicated in general) while mapping it to $\psi _{\hat {x}}$ only occasionally, when that is easy to do with a desired accuracy.

4.3. Geometrical optics of M-waves: generic case

Let us now consider a generic case when GO $M$ -waves are not delta-shaped, but quasimonochromatic in both $x$ and $q$ .Footnote ²⁵ Then, one can search for a solution of (3.26) in the eikonal form, as usual:

(4.28) \begin{gather} \bar {M}(x, q) \equiv M^*(q, x) = \mathrm{e}^{\mathrm{i} \varTheta (x, q)}\mathfrak{M}(x, q), \end{gather}

where $\varTheta$ is a real phase (eikonal) and $\mathfrak{M}$ is a real amplitude. Such $M$ -waves are well studied in quantum mechanics; see, for example, Miller (Reference Miller1974, Section II). Here, we offer a concise derivation of those known results from an alternative perspective.

Assuming that $\varTheta$ is fast compared with $\mathfrak{M}$ , it satisfies the two ‘local dispersion relations’, or Hamilton–Jacobi equations, that flow from (3.26):

(4.29a) \begin{align} K(q, -\partial _{q}\varTheta ) & = \partial _x\varTheta ,\\[-10pt]\nonumber \end{align}

(4.29b) \begin{align} P(x, \partial _x\varTheta ) & = -\partial _{q}\varTheta . \end{align}

Since $K(q, -\partial _{q}\varTheta ) = K(q, p) = k$ and $P(x, \partial _x\varTheta ) = P(x, k) = p$ , (4.29) can be expressed as

(4.30) \begin{gather} k = \partial _x\varTheta (x, q), \quad p = -\partial _q\varTheta (x, q). \end{gather}

As a function that satisfies (4.30), $\varTheta (x, q)$ can be recognised as the type-1 generating function of the canonical transformation (3.4). Generating functions of other types emerge if one re-attributes the old or (and) new coordinates as the momenta, and vice versa, with the appropriate sign changes to preserve symplecticity (§ 3.5.3). For example, in canonical variables $(q', p') \doteq (p, -q)$ , (4.30) becomes

(4.31) \begin{gather} k = \partial _x\varTheta (x, p'), \quad q' = \partial _{p'}\varTheta (x, p'), \end{gather}

which makes $\varTheta$ a type-2 generating function. Likewise, using $(x', k') \doteq (k, -x)$ , one has

(4.32) \begin{gather} x' = -\partial _{k'}\varTheta (k', q), \quad p = -\partial _q\varTheta (k', q), \end{gather}

so $\varTheta$ becomes a type-3 generating function. Finally, using both $(q', p')$ and $(x', k')$ , one can write

(4.33) \begin{gather} x' = -\partial _{k'}\varTheta (k', p'), \quad q' = \partial _{p'}\varTheta (k', p'), \end{gather}

in which case, $\varTheta$ serves as a type-4 generating function. A reader interested in brushing up on this topic is referred to Goldstein et al. (Reference Goldstein, Poole and Safko2011).

Since $K(q, -\mathrm{i} \partial _q)$ serves as the Hamiltonian for the dynamics in ‘time’ $x$ , the symbol $K(q, p)$ serves as the ray Hamiltonian. Since this Hamiltonian does not depend on $x$ explicitly, it is conserved on rays, ${\mathrm{d}}_x K = 0$ , where ${\mathrm{d}}_x$ is the convective derivative associated with the ‘group velocity’ $K_p(x, q) \doteq (\partial _p K)(q, p(x, q))$ as a field on the ‘timespace’ $(x, q)$ . Similarly, $P(x, -\mathrm{i} \partial _x)$ serves as the Hamiltonian for the dynamics in ‘time’ $q$ . Since this Hamiltonian does not depend on $q$ explicitly, it is conserved on rays, ${\mathrm{d}}_q P = 0$ , where ${\mathrm{d}}_q$ is now the convective derivative associated with the ‘group velocity’ $P_k(x, q) \doteq (\partial _k P)(x, k(x, q))$ as a field on the ‘timespace’ $(q, x)$ . From

(4.34a) \begin{align} & 0 = {\mathrm{d}}_x K = \partial _x k(x, q) + K_p \partial _{q} k(x, q),\\[-10pt]\nonumber \end{align}

(4.34b) \begin{align} & 0 = {\mathrm{d}}_q P = \partial _{q} p(x, q) + P_k \partial _x p(x, q), \end{align}

one obtains the following equations that will be useful in the next section:

(4.35a) \begin{align} & K_p(x, q) = - \frac {\partial _x k(x, q)}{\partial _q k(x, q)} = - \frac {\varTheta _{xx}}{\varTheta _{xq}},\\[-10pt]\nonumber \end{align}

(4.35b) \begin{align} & P_k(x, q) = - \frac {\partial _q p(x, q)}{\partial _x p(x, q)} = - \frac {\varTheta _{qq}}{\varTheta _{xq}}. \end{align}

The lower indices in $\varTheta \equiv \varTheta (x, q)$ denote the corresponding partial derivatives and the argument $(x, q)$ will now be occasionally omitted for brevity.

To calculate the envelope $\mathfrak{M}$ , let us apply the action conservation theorem (4.15) to $M$ -waves. For the dynamics on timespace $\boldsymbol{x} = (x, q)$ , one has $\boldsymbol{V} = (1, K_p)$ and $|\varPsi |^2 = \mathfrak{M}^2$ , so (4.15) leads to $\partial _x \mathfrak{M}^2 + \partial _q (K_p \mathfrak{M}^2) = 0$ . Similarly, for the dynamics on timespace $\boldsymbol{x} = (q, x)$ , one has $\boldsymbol{V} = (1, P_k)$ and $|\varPsi |^2 = \mathfrak{M}^2$ , so (4.15) leads to $\partial _q \mathfrak{M}^2 + \partial _x (P_k \mathfrak{M}^2) = 0$ . Let us rewrite these equations as follows:

(4.36a) \begin{align} \partial _x \mathfrak{M}^2 + (\partial _q K_p) \mathfrak{M}^2 + K_p \partial _q \mathfrak{M}^2 & = 0,\\[-10pt]\nonumber \end{align}

(4.36b) \begin{align} \partial _q \mathfrak{M}^2 + (\partial _x P_k) \mathfrak{M}^2 + P_k \partial _x \mathfrak{M}^2 & = 0. \end{align}

Substituting $\partial _q \mathfrak{M}^2$ from the second equation into the first one yields

(4.37) \begin{gather} \ln \mathfrak{M}^2 = \int {\mathrm{d}} x\,\frac {K_p(\partial _x P_k) - \partial _q K_p}{1 - P_k K_p} = \ln |\varTheta _{xq}| + \text{const.}, \end{gather}

where the latter equality is obtained by direct calculation using (4.35). The integration constant is independent of $x$ . Since one can equally arrive at (4.37) by excluding $\partial _x \mathfrak{M}^2$ instead of $\partial _q \mathfrak{M}^2$ and integrating over $q$ instead of $x$ , this constant cannot depend on $q$ either. Thus, the GO solution for $M$ -waves is

(4.38) \begin{gather} \bar {M}(x, q) = C \sqrt {\varTheta _{xq}}\,\mathrm{e}^{\mathrm{i}\varTheta }, \end{gather}

where $C$ is a constant. This factor is determined by the normalisation condition (3.25):

(4.39) \begin{align} 1 \approx |C|^2 \int {\mathrm{d}} x\,{\mathrm{d}} s |\varTheta _{xq}(x, q)|\,\mathrm{e}^{-\mathrm{i} \varTheta _{q}(x, q) s} = |C|^2 \int {\mathrm{d}} \xi \,{\mathrm{d}} s\,\mathrm{e}^{-\mathrm{i} \xi s} = 2\pi |C|^2, \end{align}

so $|C| = (2\pi )^{-1/2}$ . (Here, we introduced a variable transformation $x \mapsto \xi \doteq \varTheta _q(x, q)$ at fixed $q$ and used that ${\mathrm{d}}\xi = \varTheta _{qx}{\mathrm{d}} x$ .) Then, up to a piecewise-constant phase, one has

(4.40) \begin{gather} \bar {M}(x, q) = \sqrt {\frac {\varTheta _{xq}}{2\pi }}\,\mathrm{e}^{\mathrm{i}\varTheta }, \end{gather}

and, as a reminder, $\varTheta _{xq} = \partial _q k = - \partial _x p$ .Footnote ²⁶ The initial phase is determined by the preferred gauge (§ 3.2). Once it is set, the remaining phase is determined by the branch of the square root and changes by $\pm \pi /2$ when $\varTheta _{xq}$ goes through zero. This is the standard Maslov-index problem, which we will not revisit here, but see Miller (Reference Miller1974, pp. 85–86).

An example illustrating (4.40) will be discussed in Appendix B.2, and this formula is also useful for deriving the quantisation condition for closed orbits (Appendix D). However, keep in mind that (4.40) is constructed only as the leading-order approximation and we will also need alternative approximations of $M$ -waves for our purposes (§ 5.4).

Finally, one might wonder why (4.40) is so different from (4.24), which seems to be based on similar approximations. The explanation is as follows. In § 4.2, we used a type-2 generating function $F(x, p)$ to obtain $q = q(x, p)$ and $k = k(x, p)$ . Normally, one can invert the former to obtain $p = p(x, q)$ and $k = k(x, p(x, q))$ ; then, one can reformulate this as a type-1 transformation and the results of the present section apply. However, the specific transformation in § 4.2 is singular in that, by (4.21), our specific $q(x, p)$ happens to be independent of $p$ . This makes the inversion impossible and leads to a delta-shaped $M$ -wave (4.24).

5.1. Preliminaries

5.1.1. Tangent space

Our goal here is to construct a phase-space coordinate system aligned with the reference ray. We will assume that the ray is governed by the Hamiltonian $H_{\mathsf{z}}$ and thus satisfies $H_{\mathsf{ z}}({\mathsf{z}}) = 0$ , like in § 4.1. (However, other Hamiltonians can be advantageous or even necessary in some cases, as will be discussed in § 7.) Consider a narrow region of a fixed phase-space point ${\mathsf{z}}_0$ that satisfies $H_{\mathsf{z}}({\mathsf{z}}_0) = 0$ . Assuming that $H_{\mathsf{z}}$ is sufficiently smooth (the quantitative condition for this will be discussed in § 5.1.2), let us consider its Taylor expansion in ${\mathsf{z}} - {\mathsf{z}}_0$ :

(5.1) \begin{gather} H_{\mathsf{z}}({\mathsf{z}}) = ({\mathsf{z}} - {\mathsf{z}}_0) \boldsymbol{\cdot} \underbrace {\partial _{\mathsf{z}} H_{\mathsf{z}}({\mathsf{z}}_0)}_{-{\mathsf{J}} {\mathsf{v}}} +\, {\varUpsilon }, \quad {\mathsf{v}} \doteq {\mathsf{J}} \partial _{\mathsf{z}} H_{\mathsf{z}}({\mathsf{z}}_0), \end{gather}

where we used ${\mathsf{J}}{\mathsf{J}} = -{{\mathsf{I}}}$ . Here, ${\mathsf{v}}$ is understood as the phase-space velocity at ${\mathsf{z}}_0$ , so we will also write it as ${\mathsf{v}} = \dot {{\mathsf{z}}}_0$ . (The index 0 denotes that the corresponding quantity is evaluated at ${\mathsf{z}}_0$ .) The term ${\varUpsilon }$ subsumes the second- and higher-order terms of the Taylor series. Albeit non-negligible, it is small at small ${\mathsf{z}} - {\mathsf{ z}}_0$ , so GO applicability is determined entirely by the first term, $-({\mathsf{z}} - {\mathsf{ z}}_0) \boldsymbol{\cdot} {\mathsf{J}} {\mathsf{v}}$ . Because this term is linear in ${\mathsf{z}}$ , a linear variable transformation is generally enough to reinstate GO at least locally. This is done as follows.

Assuming the same notation as earlier, let us adopt the operator transformation (3.1) in the form

(5.2) \begin{gather} \hat {{\mathsf{y}}} = {{\mathsf{S}}} (\hat {{\mathsf{z}}} - {\mathsf{z}}_0), \end{gather}

where ${{\mathsf{S}}}$ is a symplectic matrix. As discussed in § 3.5.3, this corresponds to the phase-space transformation ${\mathsf{y}}({\mathsf{z}}) \equiv {\mathsf{y}}_{\mathsf{z}} = {{\mathsf{S}}} ({\mathsf{z}} - {\mathsf{z}}_0)$ or, equivalently, ${\mathsf{z}}({\mathsf{y}}) \equiv {\mathsf{z}}_{\mathsf{y}} = {\mathsf{z}}_0 + \bar {{{\mathsf{S}}}}{\mathsf{y}}$ , where $\bar {{{\mathsf{S}}}}$ is a symplectic matrix that is the inverse of ${{\mathsf{S}}}$ . Let us represent $\bar {{{\mathsf{S}}}}$ through its columns ${\mathsf{s}}_q$ and ${\mathsf{s}}_p$ ,

(5.3) \begin{gather} \bar {{{\mathsf{S}}}} = ( \begin{array}{c@{\quad}c} {\mathsf{s}}_q & {\mathsf{s}}_p \end{array} ), \end{gather}

and assume the notation ${\mathsf{y}} = (q, p)^{\intercal }$ , as before. Then, $q$ and $p$ can be understood as the components of ${\mathsf{z}} - {\mathsf{z}}_0$ in the basis formed by the vectors ${\mathsf{s}}_q$ and ${\mathsf{s}}_p$ :

(5.4) \begin{gather} {\mathsf{z}} - {\mathsf{z}}_0 = {\mathsf{s}}_q q + {\mathsf{s}}_p p. \end{gather}

Also, the transformation (5.2) can be expressed through (5.3) using symplecticity of ${{\mathsf{S}}}$ ; that is, inverting ${{\mathsf{S}}} {\mathsf{J}} {{\mathsf{S}}}^{\intercal } = {\mathsf{J}}$ yields ${{\mathsf{S}}} = - {\mathsf{J}}\bar {{{\mathsf{S}}}}^{\intercal }{\mathsf{J}}$ (see also (3.53)).

Let us choose ${\mathsf{s}}_q$ such that it is parallel to ${\mathsf{v}}$ , i.e. the $q$ -axis is tangent to the ray,

(5.5) \begin{gather} {\mathsf{s}}_q = {\mathsf{v}}/v. \end{gather}

Because of this, we will call ${\mathsf{y}}$ a tangent phase space. The non-zero scalar $v$ can be anything for now (but this freedom will be removed in § 5.1.2 for nonlinear coordinate transformations, where we also comment on the physical meaning of $v$ ). It may be tempting to choose ${\mathsf{s}}_p$ such that it is perpendicular to ${\mathsf{s}}_q$ ; however, angles are not well defined on a symplectic space for the lack of a metric.Footnote ²⁷ Instead, the only natural constraint on ${\mathsf{s}}_p$ is that $\bar {{{\mathsf{S}}}}$ must be symplectic, i.e. ${\mathsf{J}} = \bar {{{\mathsf{ S}}}}^{\intercal } {\mathsf{J}} \bar {{{\mathsf{S}}}}$ :

(5.6) \begin{gather} {\mathsf{J}} = \left ( \begin{array}{c} {\mathsf{s}}_q^{\intercal } \\[6pt] {\mathsf{s}}_p^{\intercal } \end{array} \right ) {\mathsf{J}} \left ( \begin{array}{c@{\quad}c} {\mathsf{s}}_q & {\mathsf{s}}_p \end{array} \right ) = \left ( \begin{array}{c@{\quad}c} {\mathsf{s}}_q \boldsymbol{\cdot} {\mathsf{J}} {\mathsf{s}}_q & {\mathsf{s}}_q \boldsymbol{\cdot} {\mathsf{J}} {\mathsf{s}}_p\\[2pt] {\mathsf{s}}_p \boldsymbol{\cdot} {\mathsf{J}} {\mathsf{s}}_q & {\mathsf{s}}_p \boldsymbol{\cdot} {\mathsf{J}} {\mathsf{s}}_p\\[2pt] \end{array} \right ) = ({\mathsf{s}}_q \wedge {\mathsf{s}}_p) {\mathsf{J}}. \end{gather}

Hence, the requirement that $\bar {{{\mathsf{S}}}}$ must be symplectic is equivalent to the requirement that ${\mathsf{s}}_q \wedge {\mathsf{s}}_p = 1$ , whence

(5.7) \begin{gather} {\mathsf{s}}_p = {\mathsf{u}}/u, \quad u \doteq {\mathsf{v}} \wedge {\mathsf{u}}/v. \end{gather}

We choose ${\mathsf{u}}$ such that $u$ is non-zero. Assuming the notation ${\mathsf{ v}} = (a_1, a_2)^{\intercal }$ and ${\mathsf{u}} = (b_1, b_2)^{\intercal }$ , one has ${\mathsf{v}} \wedge {\mathsf{u}} = a_1 b_2 - a_2 b_1$ , so having non-zero $u$ is equivalent to having ${\mathsf{u}}$ linearly independent from ${\mathsf{v}}$ , i.e. having (5.3) non-degenerate. Other than that, ${\mathsf{u}}$ can be chosen arbitrarily (but see also § 5.1.2).

Because the transformation (5.3) is linear, it preserves the symbol of $\hat {H}$ , so $H_{\mathsf{y}}({\mathsf{y}}) = H_{\mathsf{ z}}({\mathsf{z}}) = {\varUpsilon } - (\bar {{{\mathsf{S}}}} {\mathsf{y}})^{\intercal } {\mathsf{J}} {\mathsf{ v}}$ . Notice that

(5.8) \begin{align} - {\mathsf{y}}^{\intercal } \left ( \begin{array}{c} {\mathsf{s}}_q^{\intercal }\\[7pt] {\mathsf{s}}_p^{\intercal } \end{array} \right ) {\mathsf{J}} {\mathsf{v}} = - {\mathsf{y}}^{\intercal } \left ( \begin{array}{c} {\mathsf{v}} \wedge {\mathsf{v}}/v\\ {\mathsf{u}} \wedge {\mathsf{v}}/u \end{array} \right ) = \left ( \begin{array}{c@{\quad}c} q & p \end{array} \right ) \left ( \begin{array}{c} 0\\ v \end{array} \right ) = vp, \end{align}

so, in summary, $H_{\mathsf{y}}({\mathsf{y}}) = {\varUpsilon } + vp$ . The corresponding ray equations are

(5.9) \begin{gather} \dot {{\mathsf{y}}} = {\mathsf{J}} \partial _{\mathsf{y}} H_{\mathsf{y}} \equiv \lbrace {\mathsf{y}}, H_{\mathsf{y}} \rbrace, \end{gather}

or, explicitly,

(5.10) \begin{gather} \dot {q} = v + \partial _p {\varUpsilon }, \quad \dot {p} = - \partial _q {\varUpsilon }. \end{gather}

Since ${\varUpsilon } = {\mathcal{O}}(({\mathsf{z}} - {\mathsf{z}}_0)^2)$ , these lead to $\dot {q}_0 = v$ and $\dot {p}_0 = 0$ . In this sense, $v$ can be understood as the ‘symplectic speed’ at ${\mathsf{z}}_0$ . However, $v$ must not be confused with the length of the phase-space-velocity vector ${\mathsf{v}}$ , because the vector length is generally undefined on a symplectic space (see the footnote Footnote 10).

Also note that, since ${\varUpsilon }$ is small, the wavelength $2\pi /p$ evolves slowly. This justifies the GO approximation and (5.9) a posteriori. As one can also see easily, (5.9) are equivalent to the ray equations in the $x$ -representation, (4.11). In other words, the ray equations (4.11) remain applicable even near cutoffs, where GO is inapplicable in the $x$ -representation. In practice, this fact is often taken for granted, but it is non-trivial and requires justification. As seen from the previous argument, such ‘unreasonable effectiveness’ of the ray equations is due to the existence of a symplectic transformation ${\mathsf{z}} \mapsto {\mathsf{y}}$ that reinstates GO locally for any smooth $H_{\mathsf{z}}$ while conserving the Poisson bracket and the dispersion-operator symbol.

With GO reinstated in the $q$ -representation, one can also search for the wavefield in the eikonal form $\psi _{\hat {q}} = \mathrm{e}^{\mathrm{i}\theta (q)}\varPsi (q)$ and derive the envelope operator $\hat {\mathcal{H}}_{\mathsf{y}}$ for the envelope $\varPsi$ in the $q$ -space the same way the envelope equation was derived in § 4.1 in the $x$ -space. If ${\varUpsilon }$ is small globally, i.e. if $H_{\mathsf{z}}$ is mostly linear in ${\mathsf{z}}$ (locally) at all ${\mathsf{z}}$ , this envelope equation will remain sufficient indefinitely. However, such dispersion operators are rarely found in practice. A typical $H_{\mathsf{z}}$ is a nonlinear function of $\mathsf{z}$ , so the choice of the $q$ -representation that is convenient at some ${\mathsf{z}}$ is likely to become inconvenient at other ${\mathsf{z}}$ . Then, a single linear transformation is not enough.

There are two ways to deal with this problem. One is to define different linear transformations at different points ${\mathsf{z}}_0^{(i)}$ along the ray and use MTs to map the wavefield from the $q$ -space at ${\mathsf{z}}_0^{(i)}$ to the $q$ -space at ${\mathsf{z}}_0^{(i + 1)}$ (figure 1 a). For example, one can choose ${\mathsf{z}}_0^{(i + 1)}$ to be infinitesimally close to ${\mathsf{z}}_0^{(i)}$ ; then, the MT for remapping the field will be a near-identity transformation, which can be convenient. This approach is summarised in (Lopez & Dodin Reference Lopez and Dodin2022). The alternative is to use just one but nonlinear transformation such that the new coordinate space is aligned with the ray indefinitely (figure 1 b). No field remapping is needed in this case other than to initialise the field in the new coordinate space and to output the final results back to the $x$ -space. This is the approach that we discuss below.

Figure 1.

Two ways of constructing the ${\mathsf{y}}$ -coordinate grids around a reference ray (blue; coincides with the dispersion surface) in the ${\mathsf{z}}$ -space. (a) Grids are constructed in patches near predefined locations (red dots) on the ray via linear variable transformations. The $q$ -axes (red) are tangent to the ray. (b) A single grid is constructed via a nonlinear variable transformation. The coordinate axis (red) coincides with the ray indefinitely.

5.1.2. Ray path as a coordinate

Let us consider a generic ray that is curved, i.e. has $\ddot {{\mathsf{z}}}$ linearly independent from $\dot {{\mathsf{z}}}$ in the region of interest. (Straight rays can be handled as a limit; see later.) Then, at any given ${\mathsf{z}}_0$ on the ray, the vectors

(5.11) \begin{gather} {\mathsf{v}} = \dot {{\mathsf{z}}}_0, \quad {\mathsf{u}} = \ddot {{\mathsf{z}}}_0 \end{gather}

form a symplectic basisFootnote ²⁸

(5.12) \begin{gather} {\mathsf{s}}_q = {\mathsf{v}}/v, \quad {\mathsf{s}}_p = {\mathsf{u}}/u, \quad {\mathsf{s}}_q \wedge {\mathsf{s}}_p = 1, \end{gather}

where $u \doteq {\mathsf{v}} \wedge {\mathsf{u}}/v$ and the scalar $v$ remains to be defined. This basis can be used to construct new ray-aligned coordinates ${\mathsf{r}} = (\tau , h)$ (here, ‘ ${\mathsf{r}}$ ’ stands for ‘ray’) such that the ray itself would be an isoline of the new momentum $h$ . Then, one can derive the envelope equation as a partial differential equation (PDE) in the new-coordinate representation. The envelope remains slow in this representation because the reference wavevector satisfies the dispersion relation indefinitely and there are no cutoffs in this representation by construction of ${\mathsf{ r}}$ .

Since we are interested only in the local dynamics on the ray, it will be sufficient for us to find ${\mathsf{r}}({\mathsf{z}})$ and said PDE near a predefined generic ${\mathsf{z}}_0$ on the ray. Then, it is convenient to search for the mapping ${\mathsf{z}} \mapsto {\mathsf{r}}$ as a superposition of two transformations:

(5.13) \begin{gather} \underbrace {(x, k)^{\intercal }}_{{\mathsf{z}}} \mapsto \underbrace {(q, p)^{\intercal }}_{{\mathsf{y}}} \mapsto \underbrace {(\tau , h)^{\intercal }}_{\mathsf{r}}. \end{gather}

The first one, ${\mathsf{y}} = \bar {{{\mathsf{S}}}} ({\mathsf{z}} - {\mathsf{z}}_0)$ , is the same as in § 5.1.1, except now ${\mathsf{ u}}$ is specified by (5.11); i.e. the $q$ -axis is tangent to the ray at ${\mathsf{z}} = {\mathsf{z}}_0$ . The second transformation, ${\mathsf{y}} \mapsto {\mathsf{r}}$ , is nonlinear and modifies the coordinate grid such that the new coordinate axis coincides with the ray. (The momentum axis is also modified accordingly, such that the canonical commutation relation (3.2) is preserved.) The advantage of the splitting (5.13) is that the mapping ${\mathsf{z}} \mapsto {\mathsf{y}}$ is easy to construct exactly and, for smooth rays, ${\mathsf{y}} \mapsto {\mathsf{ r}}$ can also be found analytically, at least asymptotically. In the following, we show how to do this step by step.

5.2. ‘Osculating’ harmonic oscillator

Let us start with a simplified model, where $H_{\mathsf{z}}$ can be adequately approximated (locally) with its second-order Taylor expansion in ${\mathsf{z}} - {\mathsf{z}}_0$ :Footnote ²⁹

(5.14) \begin{gather} H_{\mathsf{z}}({\mathsf{z}}) = - ({\mathsf{z}} - {\mathsf{z}}_0) \boldsymbol{\cdot} {\mathsf{J}} {\mathsf{v}} + \frac {1}{2}\,({\mathsf{z}} - {\mathsf{z}}_0) \boldsymbol{\cdot} {\mathsf{g}}({\mathsf{z}} - {\mathsf{z}}_0). \end{gather}

Here, ${\mathsf{v}} \doteq {\mathsf{J}} \partial _{\mathsf{z}} H_{\mathsf{z}}({\mathsf{z}}_0) \equiv \dot {{\mathsf{z}}}_0$ , as earlier, and ${\mathsf{g}}$ is a symmetric matrix given by ${\mathsf{g}} \doteq (\partial ^2_{\mathsf{z}} H)_0$ . A related description of this model can also be found in Tracy et al. (Reference Tracy, Brizard, Richardson and Kaufman2014, § 5.3.3) for a specific form of $H_{\mathsf{z}}$ . Here, we allow for a general $H_{\mathsf{z}}$ to keep the theory symplectically invariant, at least in its final form discussed in § 5.4.

Equation (5.14) leads to the following Hamilton’s equations:

(5.15) \begin{gather} \dot {{\mathsf{z}}} = {\mathsf{v}} + {\mathsf{J}} {\mathsf{g}}({\mathsf{z}} - {\mathsf{z}}_0), \quad \ddot {{\mathsf{z}}} = {\mathsf{J}} {\mathsf{g}}\dot {{\mathsf{z}}}, \end{gather}

so $\ddot {{\mathsf{z}}}_0 = ({\mathsf{J}} {\mathsf{g}}\dot {{\mathsf{z}}})_0 = {\mathsf{J}} {\mathsf{ g}}{\mathsf{v}}$ . Then, ${\mathsf{u}} = {\mathsf{J}} {\mathsf{g}}{\mathsf{v}}$ and $u = - \varepsilon /v$ . Here, $v$ is a normalisation factor that was introduced in (5.5) and remains to be specified. Also,

(5.16) \begin{gather} \varepsilon \doteq (\ddot {{\mathsf{z}}} \wedge \dot {{\mathsf{z}}})_0 = ({\mathsf{J}} {\mathsf{g}} {\mathsf{v}})^{\intercal } {\mathsf{J}} {\mathsf{v}} = {\mathsf{v}} \boldsymbol{\cdot} {\mathsf{g}}{\mathsf{v}}, \end{gather}

where we used ${\mathsf{J}}^{\intercal }{\mathsf{J}} = {{\mathsf{I}}}$ and ${\mathsf{ g}}^{\intercal } = {\mathsf{g}}$ . Also note that

(5.17) \begin{align}&\qquad\qquad\quad {\mathsf{v}} \boldsymbol{\cdot} {\mathsf{g}}{\mathsf{u}} = {\mathsf{v}} \boldsymbol{\cdot} {\mathsf{g}} {\mathsf{J}} {\mathsf{g}}{\mathsf{v}} = ({\mathsf{g}}{\mathsf{v}}) \wedge ({\mathsf{g}}{\mathsf{v}}) = 0,\\[-10pt]\nonumber \end{align}

(5.18) \begin{align}& {\mathsf{u}} \boldsymbol{\cdot} {\mathsf{g}} {\mathsf{u}} = ({\mathsf{J}} {\mathsf{g}}{\mathsf{v}})^{\intercal } {\mathsf{g}} ({\mathsf{J}} {\mathsf{g}}{\mathsf{v}}) = {\mathsf{v}}^{\intercal } {\mathsf{g}}^{\intercal } {\mathsf{J}}^{\intercal } {\mathsf{g}} {\mathsf{J}} {\mathsf{g}} {\mathsf{v}} = {\mathsf{v}} \boldsymbol{\cdot} {\mathsf{g}} ({-}{\mathsf{J}} {\mathsf{g}} {\mathsf{J}} {\mathsf{g}}) {\mathsf{v}} = g \varepsilon , \end{align}

where we used ${\mathsf{J}}^{\intercal } = -{\mathsf{J}}$ , and, for $2 \times 2$ matrices, $-{\mathsf{J}} {\mathsf{g}} {\mathsf{J}} {\mathsf{g}} = g{{\mathsf{I}}}$ , where $g \doteq \det {\mathsf{g}}$ . Similarly, the time derivative of $\ddot {{\mathsf{z}}} \wedge \dot {{\mathsf{z}}}$ is

(5.19) \begin{gather} \dddot {{\mathsf{z}}} \wedge \dot {{\mathsf{z}}} = ({\mathsf{J}} {\mathsf{g}} \ddot {{\mathsf{z}}}) \wedge \dot {{\mathsf{z}}} = \ddot {{\mathsf{z}}}^{\intercal } {\mathsf{g}}^{\intercal } {\mathsf{J}}^{\intercal } {\mathsf{J}} \dot {{\mathsf{z}}} = \ddot {{\mathsf{z}}}^{\intercal } {\mathsf{g}} \dot {{\mathsf{z}}} = \ddot {{\mathsf{z}}}^{\intercal } {\mathsf{J}}^{\intercal } ({\mathsf{J}} {\mathsf{g}} \dot {{\mathsf{z}}}) = -\ddot {{\mathsf{z}}} \wedge \ddot {{\mathsf{z}}} = 0, \end{gather}

so $\ddot {{\mathsf{z}}} \wedge \dot {{\mathsf{z}}} = (\ddot {{\mathsf{z}}} \wedge \dot {{\mathsf{ z}}})_0$ . This means that, to the extent that the model (5.14) is applicable, one can replace the definition (5.16) with $\varepsilon = \ddot {{\mathsf{z}}} \wedge \dot {{\mathsf{z}}}$ .

Let us use these results to see how (5.14) is transformed by the mapping ${\mathsf{z}} \mapsto {\mathsf{y}}$ . The first term in (5.14) becomes $p v$ , like before. The second term becomes

(5.20) \begin{align} \frac {1}{2}\,({\mathsf{z}} - {\mathsf{z}}_0) \boldsymbol{\cdot} {\mathsf{g}}({\mathsf{z}} - {\mathsf{z}}_0) & = \frac {1}{2}\,{\mathsf{y}}^{\intercal } \left ( \begin{array}{c} {\mathsf{v}}^{\intercal }/v\\[4pt] {\mathsf{u}}^{\intercal }/u \end{array} \right ) {\mathsf{g}} \left ( \begin{array}{c@{\quad}c} {\mathsf{v}}/v & {\mathsf{u}}/u \end{array} \right ) {\mathsf{y}} \notag \\[4pt] & = \frac {1}{2}\,{\mathsf{y}} \boldsymbol{\cdot} \left ( \begin{array}{c@{\quad}c} {\mathsf{v}} \boldsymbol{\cdot} {\mathsf{g}} {\mathsf{v}}/v^2 & {\mathsf{v}} \boldsymbol{\cdot} {\mathsf{g}} {\mathsf{u}}/vu\\[4pt] {\mathsf{u}} \boldsymbol{\cdot} {\mathsf{g}} {\mathsf{v}}/vu & {\mathsf{u}} \boldsymbol{\cdot} {\mathsf{g}} {\mathsf{u}}/u^2 \end{array} \right ) {\mathsf{y}} \notag \\[4pt] & = \frac {1}{2}\,{\mathsf{y}} \boldsymbol{\cdot} \left ( \begin{array}{c@{\quad}c} \varepsilon /v^2 & 0\\[4pt] 0 & gv^2/\varepsilon \end{array} \right ) {\mathsf{y}}. \end{align}

Assuming the notation $\sigma _f \doteq {\textrm {sgn}} f$ for any $f$ , let us introduce $\varOmega \doteq \sigma _\varepsilon \sqrt {|g|}$ and choose $v = (\varepsilon ^2/|g|)^{1/4}$ , so that the absolute values of the diagonal elements in the matrix in (5.20) are the same. Then, $\varepsilon /v^2 = |g|v^2/\varepsilon = \varOmega$ , so $v = |\varepsilon /\varOmega |^{1/2}$ , $u = - \sigma _\varepsilon |\varepsilon \varOmega |^{1/2}$ and

(5.21) \begin{gather} \left ( \begin{array}{c@{\quad}c} \varepsilon /v^2 & 0\\[4pt] 0 & gv^2/\varepsilon \end{array} \right ) = \varOmega \left ( \begin{array}{c@{\quad}c} 1 & 0\\[2pt] 0 & \sigma _g \end{array} \right )\!. \end{gather}

Then,

(5.22) \begin{gather} H_{{\mathsf{y}}}({\mathsf{y}}) = H_{\mathsf{z}}({\mathsf{z}}) = \frac {\varOmega }{2}(q^2 + \sigma _g (p + R)^2 - \sigma _g R^2), \end{gather}

where $R \doteq \sigma _g v/\varOmega$ , so $|R| = |v/\varOmega | = |\varepsilon |^{1/2}|g|^{-3/4}$ .Footnote ³⁰ One can recognise this procedure as the Williamson diagonalisation of ${\mathsf{g}}$ (at least, for $g \gt 0$ ), with $\varOmega$ being the symplectic eigenvalue of ${\mathsf{g}}$ (Nicacio Reference Nicacio2021). Note that $\varOmega$ and $R$ are symplectic invariants by definition, and so are $\varepsilon$ , $g$ and $v$ . Also, it is easily seen that $R$ is dimensionless and the other quantities have the following units (the square brackets here stand for ‘units of ’):

(5.23) \begin{gather} [\varepsilon ] = [H_{\mathsf{z}}]^3, \quad [g] = [H_{\mathsf{z}}]^2, \quad [\varOmega ] = [H_{\mathsf{z}}], \quad [v] = [H_{\mathsf{z}}]. \end{gather}

Also useful for us will be the following formula:

(5.24) \begin{gather} |\varepsilon |^{1/3} = |\varOmega \| R|^{2/3} = v/|R|^{1/3}. \end{gather}

The phase-space variables $(q, p)^{\intercal } \equiv {\mathsf{y}}$ can be understood as a local approximation to the natural canonical variables of $\hat {H}$ . The matrix ${\mathsf{g}}$ serves as a natural (pseudo)metric induced by $\hat {H}$ on the symplectic space; then, $\varepsilon = \|{\mathsf{v}}\|_{\mathsf{g}}^2$ can be understood as the square of ${\mathsf{ v}}$ in this metric. The ray trajectories in the ${\mathsf{y}}$ -space satisfy $H_{{\mathsf{ y}}}({\mathsf{y}}) = 0$ , so, within the approximation (5.22), they are conic sections: circles at $g \gt 0$ , parabolas at $g = 0$ and hyperbolas at $g \lt 0$ . However, we are interested only in small $q$ and $p$ , in which case $\sigma _g$ has little effect on the ray trajectories (figure 2).Footnote ³¹ Hence, locally, the ray dynamics can be thought of as the dynamics of a harmonic oscillator with the Hamiltonian

(5.25) \begin{gather} H_{{\mathsf{y}}}({\mathsf{y}}) = \frac {\varOmega }{2}\,(q^2 + (p + R)^2 - R^2). \end{gather}

We call it the ‘osculating’ harmonic oscillator (OHO). (To a dynamical system, an OHO basically is what an osculating circle is to a curve.) The phase-space orbit of OHO is a circle with radius $|R|$ . Since our approximation (5.14) is generally adequate only on scales small compared with this scale, our model is applicable only when both $q$ and $p$ are much less than $R$ .

Figure 2.

Ray trajectories $H_{{\mathsf{y}}}({\mathsf{y}}) = 0$ for various signs of $g$ and $\varepsilon$ . Darker colours mark, loosely, the areas of where the assumed model is adequate. Lighter colours mark the areas that are beyond the validity domain of the model. The arrows mark the direction of the ray propagation, which is always towards positive $q$ . In the two right columns, dashed are the osculating circles. The radius of each circle is $|R| = |v/\varOmega |$ . Dotted are the asymptotes of the hyperbolas.

The parameter $\varOmega$ can be understood as OHO’s angular frequency. At $\varOmega \gt 0$ , or $\varepsilon \gt 0$ , the phase-space trajectory is oriented clockwise, as usual. (Note that the canonical phase of the oscillator normally grows when the polar angle on the $(q, p)$ plane decreases.) At $\varOmega \lt 0$ , or $\varepsilon \lt 0$ , the trajectory is counter-clockwise. At $\varOmega = 0$ , which corresponds to $\varepsilon = 0$ , the trajectory is flat, i.e. $R = \infty$ . Hence, it is convenient to classify the ray trajectories in terms of their (dimensionless) ‘symplectic curvature’

(5.26) \begin{gather} \kappa \doteq \frac {1}{R} = \frac {\varOmega }{v} = \sigma _\varepsilon \,\frac {|g|^{3/4}}{|\varepsilon |^{1/2}}, \end{gather}

and $R$ will be called the symplectic radius, accordingly. (Here and further, ${\textrm {sgn}} R \equiv {\textrm {sgn}}\, \kappa = \sigma _\varepsilon$ .) Note that the symplectic curvature can be very different from the regular curvature of a trajectory in a metric space. In particular, one can easily show that

(5.27) \begin{gather} \varepsilon = - \frac {{\mathrm{d}}^2 x}{{\mathrm{d}} k^2}\left (\frac {\partial H_{\mathsf{z}}}{\partial x}\right )^3 = - \frac {{\mathrm{d}}^2 k}{{\mathrm{d}} x^2}\left (\frac {\partial H_{\mathsf{z}}}{\partial k}\right )^3, \end{gather}

where the full derivatives are taken along the ray trajectory, $H_{{\mathsf{z}}}({\mathsf{z}}) = 0$ . Thus, at inflection points (where $\varepsilon \to 0$ ), $\kappa$ is infinite, whereas the regular curvature would be zero.

Finally, notice the following. The OHO approximation (5.25) for $H_{\mathsf{z}}$ corresponds to the following approximation of $\hat {H}$ itself: $\hat {H} \approx (\varOmega /2)(\hat {q}^2 + (\hat {p} + R)^2 - R^2)$ . This means that, in the absence of dissipation, the wave (4.16) becomes the equation of a quantum harmonic oscillator (QHO) with energy $\varOmega R^2/2$ , in the natural units that correspond to unit Planck constant, $\hbar = 1$ . Since the QHO energy can be expressed through its quantum mode number $n$ as $\varOmega (n + 1/2)$ , one obtains

(5.28) \begin{gather} n = \frac {R^2}{2} - \frac {1}{2}. \end{gather}

This readily determines applicability conditions for GO in the $q$ -representation in the region $q \ll R$ in which we are interested. Indeed, as commonly known, one can adequately approximate the dynamics of such QHO using GO at $n \gg 1$ , i.e. $R \gg 1$ . This makes $\kappa$ a natural small parameter for MGO. One might be surprised that this parameter becomes infinite at inflection points, where rays become straight and thus, supposedly, MGO should work only better. However, notice that, near inflection points, $R$ evolves faster than the ray moves along the OHO trajectory per se and, thus, the very concept of the orbit radius becomes irrelevant. In other words, the quadratic approximation (5.14) is inapplicable in such regimes and one should either include cubic terms or simply ignore the effects produced by the quadratic terms, in which case MGO is reinstated. What matters then is not the local $R$ per se, but the characteristic $R$ . We will return to this subject in § 5.4.

5.3. Parabolic approximation

There is more than one way to define the transformation ${\mathsf{y}} \mapsto {\mathsf{r}}$ near a given ${\mathsf{z}}_0$ . One set of obvious variables to try are the angle–action variables of the OHO. However, as discussed in Appendix B, this approach is not optimal, so let us try a different one.

5.3.1. Operator transformation

As pointed out in § 5.2, the OHO approximation is valid only at $q \ll R$ , where $p \sim q^2/R \ll q$ . Because of this, one can neglect $p^2$ compared with $q^2$ in (5.22) and obtain the following ‘parabolic’ approximation of the ray Hamiltonian:Footnote ³²

(5.29) \begin{gather} H_{{\mathsf{y}}}({\mathsf{y}}) = H_{\mathsf{z}}({\mathsf{z}}) \approx pv + \frac {\varOmega }{2}\,q^2. \end{gather}

(Note that this holds for any signs of $g$ and $\varepsilon$ , so the discontinuous function ${\textrm {sgn}}$ never even appears in this formulation.) Let us now introduce new canonical variables $(\tau , h)^{\intercal } \equiv {\mathsf{r}}$ such that $h = pv + \varOmega q^2/2$ . To find a suitable $\tau$ , let us search for the corresponding generating function $F$ in the type-2 form, $F = F(q, h)$ (Goldstein et al. Reference Goldstein, Poole and Safko2011):

(5.30) \begin{gather} p(q, h) = \partial _q F(q, h), \quad \tau (q, h) = \partial _h F(q, h). \end{gather}

Since $p(q, h) = (h - \varOmega q^2/2)/v$ , one has $F(q, h) = (h q - \varOmega q^3/6)/v$ (up to an arbitrary function of $h$ , which we choose to be zero), so

(5.31) \begin{gather} \tau = q/v, \quad h = pv + \varOmega q^2/2, \end{gather}

where we have repeated the expression for $h$ for completeness. Since $v$ and $w$ are defined as values at a fixed ${\mathsf{z}}_0$ , the corresponding operators are given simply by

(5.32) \begin{gather} \hat {\tau } = \hat {q}/v, \quad \hat {h} = \hat {p}v + \varOmega \hat {q}^2/2. \end{gather}

(Basically, (5.32) is just a rescaled and truncated small- $p$ approximation of (B4).) One can readily see that they exactly satisfy the canonical commutation relation $[\hat {\tau }, \hat {h}] = \mathrm{i}$ .

5.3.2. M-waves

In addition to avoiding infinite series, the previous model is advantageous in that it allows one to explicitly calculate $M$ -waves in a closed form. Indeed, (3.23) for becomes an equation solvable analytically:

(5.33) \begin{gather} \mathrm{i} \partial _\tau \bar {M} = -\mathrm{i} v \partial _q \bar {M} + \frac {1}{2}\,w q^2 \bar {M} \end{gather}

(and (3.19) is just the complex conjugate thereof). Since $\hat {\tau } = \hat {q}/v$ , one can adopt ; then, the initial condition for (5.33) is

(5.34)

The corresponding solution of (5.33) is

(5.35) \begin{gather} \bar {M}(q, \tau ) = \frac {1}{\sqrt {v}}\,\delta \left (\tau - \frac {q}{v}\right ) \exp \left ({-}\frac {\mathrm{i} \varOmega q^3}{6v}\right )\!. \end{gather}

Accordingly, the pseudo-measure (3.29), with $M \equiv \bar {M}^*$ , is given by

(5.36) \begin{align} \mu & = \frac {1}{2\pi } \int {\mathrm{d}} s'\,{\mathrm{d}} s\, M(\tau + s/2, q + s'\!/2) \, M^*(\tau - s/2, q - s'\!/2)\,\mathrm{e}^{-\mathrm{i} h s + \mathrm{i} p s'} \notag \\ & = \frac {1}{2\pi v} \int {\mathrm{d}} s\,{\mathrm{d}} s'\, \delta \left (\tau - \frac {q}{v} + \frac {s - s'\!/v}{2}\right ) \delta \left (\tau - \frac {q}{v} - \frac {s - s'\!/v}{2}\right ) \notag \\ & \quad \times \exp \left ( \frac {\mathrm{i} \varOmega }{6v}\left (q + \frac {s'}{2}\right )^3 - \frac {\mathrm{i} \varOmega }{6v}\left (q - \frac {s'}{2}\right )^3 -\mathrm{i} h s + \mathrm{i} p s' \right ) \notag \\ & = \frac {1}{2\pi v} \int {\mathrm{d}} s\,{\mathrm{d}} s'\, \delta \left (\tau - \frac {q}{v}\right ) \delta \left (s - \frac {s'}{v}\right ) \notag \\ & \quad \times \exp \left ( \frac {\mathrm{i} \varOmega }{6v}\left (q + \frac {s'}{2}\right )^3 - \frac {\mathrm{i} \varOmega }{6v}\left (q - \frac {s'}{2}\right )^3 -\mathrm{i} h s + \mathrm{i} p s' \right ) \notag \\ & = \frac {1}{2\pi } \int {\mathrm{d}} s\, \delta \left (\tau - \frac {q}{v}\right ) \notag \\ & \quad \times \exp \left ( \frac {\mathrm{i} \varOmega }{6v}\left (q + \frac {s v}{2}\right )^3 - \frac {\mathrm{i} \varOmega }{6v}\left (q - \frac {s v}{2}\right )^3 + \mathrm{i} s (p v - h) \right ) \notag \\ & = \frac {1}{2\pi } \int {\mathrm{d}} s\, \delta \left (\tau - \frac {q}{v}\right ) \exp \left ( \frac {\mathrm{i} \varepsilon s^3}{24} + \mathrm{i} s (H_{\mathsf{y}}({\mathsf{y}}) - h) \right )\!, \end{align}

where we used $\varOmega v^2 = \varepsilon$ . This result can also be expressed as

(5.37) \begin{gather} \mu ({\mathsf{r}}, {\mathsf{y}}) = \textrm {Ai}_\varepsilon (H_{\mathsf{y}}({\mathsf{y}}) - h)\,\delta (\tau - q/v), \end{gather}

where we have introduced the following (real) function, with $\gamma$ as a real parameter:

(5.38) \begin{gather} \textrm {Ai}_\gamma (z) \doteq \frac {1}{2\pi }\int _{-\infty }^{\infty } {\mathrm{d}} t\,\mathrm{e}^{\mathrm{i} z t + \mathrm{i} \gamma t^3/24} = \frac {2}{|\gamma |^{1/3}}\,\textrm {Ai}\left (\frac {2}{|\gamma |^{1/3}}\,\sigma _\gamma z\right )\!. \end{gather}

Here, $\sigma _\gamma \doteq {\textrm {sgn}} \gamma$ , and $\textrm {Ai}$ is the Airy function of the first kind. Importantly (Appendix C),

(5.39) \begin{gather} \lim _{\gamma \to 0} \textrm {Ai}_\gamma (z) = \delta (z). \end{gather}

As to be seen, the function $\textrm {Ai}_\gamma$ , which we will call a rescaled Airy function (RAF), plays a fundamental role in MGO and emerges in multiple contexts.

5.4. General case

Let us now present a formulation that is not restricted to the quadratic approximation (5.14) and also subsumes the standard WKB approximation (4.27) as a special case.

5.4.1. Operator transformation

Let us start with assuming a generic $H_{\mathsf{z}}$ in the form

(5.40) \begin{gather} H_{\mathsf{z}}({\mathsf{z}}) = \mathcal{H}(x/\Delta x, k/\Delta k), \end{gather}

where $\Delta x$ and $\Delta k$ are the characteristic scales of the system in $x$ and $k$ , respectively, and $\mathcal{H}$ is some function with order-one scales. Consider new dimensionless variables

(5.41) \begin{gather} X \doteq R x/\Delta x, \quad K \doteq R k/\Delta k, \end{gather}

where $R$ is a dimensionless constant. We require that the transformation ${\mathsf{z}} \mapsto {\mathsf{Z}} \equiv (X, K)^{\intercal }$ be symplectic; then, $H_{\mathsf{Z}}({\mathsf{Z}}) = H_{\mathsf{z}}({\mathsf{z}})$ and

(5.42) \begin{gather} R = \sqrt {\Delta x \Delta k}. \end{gather}

We will call this quantity a symplectic scale. To the extent that the quadratic approximation (5.14) is adequate at least locally, $R$ coincides with the absolute value of the symplectic radius that was introduced in § 5.2. For this reason, in the following, we occasionally mix the terms ‘symplectic scale’ and ‘symplectic radius’, and also refer to the inverse symplectic scale as symplectic curvature.

Much like in § 4.1, our goal now is to derive field equations that are accurate to the first (not zeroth) order in the ‘MGO parameter’

(5.43) \begin{gather} \epsilon \doteq R^{-2} \equiv (\Delta x \Delta k)^{-1} \ll 1. \end{gather}

It is similar to the GO parameter introduced in § 4.1, but the absolute value of $k$ is now replaced with the spectral-inhomogeneity scale $\Delta k$ ; i.e. the smallness of the local wavelength per se is irrelevant. To the extent that the OHO model is applicable, $\epsilon$ coincides with $\kappa ^2$ , where $\kappa$ is the symplectic curvature introduced in § 5.2.

Let us perform a linear canonical transformation ${\mathsf{Z}} \mapsto {\mathsf{y}} \equiv (q, p)^{\intercal }$ , as in § 5.1.1, so that the $q$ -axis is along ${\mathsf{v}}$ that is tangent to the ray. Then, $H_{\mathsf{ y}}({\mathsf{y}}) = \mathcal{H}(q/R, p/R)$ , so the ray trajectory satisfies $\mathcal{H}(q/R, p/R) = 0$ . Its solution (for a given branch) can be written in the form $p = R\mathcal{P}(q/R) \equiv p(q)$ , where $\mathcal{P}$ is an order-one function that has order-one scales. (If $H_{\mathsf{z}}$ is linear in ${\mathsf{z}}$ , then $\mathcal{P} = 0$ by construction. If $H_{\mathsf{z}}$ is quadratic, then $p(q) = -R + \sqrt {R^2 - q^2}$ , so $p(q) \approx q^2/2R$ at $q \ll R$ , as in § 5.3.) By Taylor-expanding $H_{\mathsf{ y}}({\mathsf{ y}})$ in $p - p(q)$ , one obtains

(5.44) \begin{gather} H_{\mathsf{y}}({\mathsf{y}}) \approx V(q)(p - p(q)), \end{gather}

where $V(q) \doteq (\partial _p H_{\mathsf{y}})_{p = p(q)}$ , so $V(0) = v$ . Like (4.17), (5.44) can be understood as Hayes’s representation of $H_{\mathsf{y}}$ . For GO, this representation is precise in that it leads to the correct ray equations and correctly captures the leading-order amplitude dynamics, because those are entirely determined by the first-order derivatives of the symbol. The MT that follows from (5.44) will not be exact, because ${\mathcal{O}}([p - p(q)]^2)$ terms have been neglected. Still, assuming $p - p(q) = {\mathcal{O}}(q^2)$ , (5.44) differs from the exact $H_{\mathsf{y}}$ by at most ${\mathcal{O}}(q^4)$ , which is assumed negligible (at not-too-large $q$ that we are interested in locally).

Let us perform a canonical transformation ${\mathsf{y}} \,{\mapsto}\, {\mathsf{r}} \equiv (\tau , h)^{\intercal }$ such that $h = H_{\mathsf{y}}$ , i.e.Footnote ³³

(5.45) \begin{gather} p(q, h) = p(q) + h/V(q). \end{gather}

Like in § 5.3, to find what the new coordinate $\tau$ should be to keep the transformation canonical, we look for a generating function $F$ in the type-2 form:

(5.46) \begin{gather} p(q, h) = \partial _q F(q, h), \quad \tau (q, h) = \partial _h F(q, h). \end{gather}

The former gives

(5.47) \begin{gather} F(q, h) = \int ^{q}_0 p(\tilde {q}, h)\,{\mathrm{d}} \tilde {q}, \end{gather}

where we set the integration constant to zero by choice. (The tilde is used to distinguish the dummy integration variable $\tilde {q}$ from the actual canonical coordinate $q$ .) Then,

(5.48) \begin{gather} \tau (q, h) = \int ^{q}_0 \frac {\partial p(\tilde {q}, h)}{\partial h}\,{\mathrm{d}} \tilde {q} = \int ^{q}_0 \frac {{\mathrm{d}} \tilde {q}}{V(\tilde {q})} \equiv \tau (q). \end{gather}

In other words, $\tau$ and $h$ are the ray time and the ray energy. (By the ray energy, we mean the value of the ray Hamiltonian, not the wave energy that depends on the wave amplitude.) The corresponding operators are

(5.49) \begin{gather} \hat {\tau } = \tau (\hat {q}), \quad \hat {h} = \frac {1}{2}(\hat {V} \hat {p} + \hat {p} \hat {V}) - \hat {V} p(\hat {q}), \end{gather}

where $\hat {V} \doteq V(\hat {q})$ (note, though, that $\hat {p}$ is not the same as $p(\hat {q})$ ), so $[\hat {\tau }, \hat {V}] = 0$ and

(5.50) \begin{align} [\hat {\tau }, \hat {h}] & = [\hat {\tau }, \hat {V} \hat {p} + \hat {p} \hat {V}]/2 \notag \\ & = (\hat {\tau }\hat {V} \hat {p} - \hat {V} \hat {p}\hat {\tau } + \hat {\tau } \hat {p} \hat {V} - \hat {p} \hat {V} \hat {\tau })/2 \notag \\ & = (\hat {\tau }\hat {V} \hat {p} - \hat {V} [\hat {p}, \hat {\tau }] - \hat {V}\hat {\tau } \hat {p} + \hat {\tau } [\hat {p}, \hat {V}] + \hat {\tau } \hat {V} \hat {p} - \hat {p} \hat {V} \hat {\tau })/2 \notag \\ & = ([\hat {\tau }\hat {V}, \hat {p}] - \hat {V} [\hat {p}, \hat {\tau }] + \hat {\tau } [\hat {p}, \hat {V}])/2 \notag \\ & = \mathrm{i}((\tau (q) V(q))' + V(q) \tau '(q) - \tau (q) V'(q))/2 \notag \\ & = \mathrm{i} V(q)\tau '(q) \notag \\ & = \mathrm{i}. \end{align}

5.4.2. M-waves

The approximation (5.44) corresponds to the following approximation of the operator $\hat {H}_{\hat {q}}$ (cf. § 4.1):

(5.51) \begin{gather} \hat {H}_{\hat {q}} \approx -\mathrm{i} V(q) \partial _q - \frac {\mathrm{i}}{2}\,V'(q) - V(q)p(q). \end{gather}

Then, (3.22) and (3.23) lead to the following equations for $\bar {M} \equiv \bar {M}(q, \tau )$ :

(5.52) \begin{align}&\qquad\qquad\qquad (\tau - \tau (q))\bar {M} = 0,\\[-10pt]\nonumber \end{align}

(5.53) \begin{align}& \mathrm{i} \partial _\tau \bar {M} = \left ({-}\mathrm{i} V(q) \partial _q - \frac {\mathrm{i}}{2}\,V'(q) - V(q)p(q)\right )\bar {M}. \end{align}

These can be readily solved:

(5.54) \begin{gather} \bar {M}(q, \tau ) = \frac {\mathrm{e}^{\mathrm{i} \vartheta (q)}}{\sqrt {V(q)}}\,\delta (\tau - \tau (q)), \quad \vartheta (q) \doteq \int ^{q}_0 p(\tilde {q})\,{\mathrm{d}} \tilde {q}, \end{gather}

which generalises (5.35) to non-constant $V(q)$ and arbitrary $p(q)$ . To the extent that (5.44) is valid, this result is exact. One can also notice parallels with (4.24), which is the same formula applied to different coordinates.

To calculate the generalisation of the pseudo-measure (5.37), notice that

(5.55) \begin{gather} \tau \left (q \pm \frac {s'}{2}\right ) \approx \tau (q) \pm \frac {s'}{2V} - \frac {V'}{8V^2}\,s'^2, \end{gather}

where $V \equiv V(q)$ and $V'(q) \equiv {\mathrm{d}} V(q)/{\mathrm{d}} q$ , while the prime in $s'$ is only a notation that distinguishes said variable from $s$ . Then,

(5.56) \begin{align} & \delta \left ( \tau + \frac {s}{2} - \tau \left (q + \frac {s'}{2}\right ) \right ) \delta \left ( \tau - \frac {s}{2} - \tau \left (q - \frac {s'}{2}\right ) \right ) \notag \\[2pt] & \quad \approx \delta \left ( \tau - \tau (q) + \frac {s - s'\!/V}{2} + \frac {V'{s}^{\prime 2}}{8V^2} \right ) \delta \left ( \tau - \tau (q) - \frac {s - s'\!/V}{2} + \frac {V'{s}^{\prime 2}}{8V^2} \right ) \notag \\[2pt] & \quad = \delta \left (\tau - \tau (q) + \frac {V'{s}^{\prime 2}}{8V^2}\right ) \delta \left (s - \frac {s'}{V}\right ) \notag \\[2pt] & \quad = \delta \left (\tau - \tau (q) + \frac {V's^2}{8}\right ) \delta \left (s - \frac {s'}{V}\right ) \notag \\[2pt] & \quad \approx \left (1 + \frac {V's^2}{8}\,\frac {\partial }{\partial \tau }\right ) \delta (\tau - \tau (q)) \delta \left (s - \frac {s'}{V}\right )\!, \end{align}

where the last equality is based on the first-order Taylor expansion. (Admittedly, Taylor-expanding a delta function is a questionable procedure, but it makes sense as a shorthand for Taylor-expanding what actually matters, namely, integrals of delta functions.) Then,

(5.57) \begin{gather} \mu = \big [1 - (V'/8)\,\partial ^3_{\tau h h}\big ]\mu _0, \quad \mu _0 = \mathfrak{A}\, \delta (\tau (q) - \tau ), \end{gather}

where $\mathfrak{A}$ is a generalisation of (5.38):

(5.58) \begin{gather} \mathfrak{A} \doteq \frac {1}{2\pi } \int {\mathrm{d}} s\,\mathrm{e}^{\mathrm{i}\varPhi }, \quad \varPhi \doteq (p V - h) s - \int _{q - s V/2}^{q + s V/2} p(\tilde {q})\,{\mathrm{d}}\tilde {q}. \end{gather}

Remember that $p$ is an independent variable, $V \equiv V(q)$ and $p(q)$ is a function introduced in § 5.4.1.

5.4.3. Symplectically invariant model for $\mu$

Since we are already neglecting terms of order $V''$ , the factor $V'$ in (5.57) can be considered commuting with the derivatives. Then, using (3.41) and integrating by parts, one obtains the following formula for remapping symbols between the ${\mathsf{y}}$ - and ${\mathsf{r}}$ - representations:

(5.59) \begin{gather} O_{\mathsf{y}}({\mathsf{y}}) = \int \left (O_{\mathsf{r}}({\mathsf{r}}) + \frac {V'}{8} \,\partial ^3_{\tau h h} O_{\mathsf{r}}({\mathsf{ r}}) \right )\mu _0 \,{\mathrm{d}}\tau \,{\mathrm{d}} h. \end{gather}

In the following, we will be interested in $O_{\mathsf{r}}$ that are either $\tau$ -independent or smooth, namely, have symplectic scales of order $R$ or larger. In either case, the second term in (5.59) is ${\mathcal{O}}(R^{-4})$ , i.e. ${\mathcal{O}}(\epsilon ^2)$ , which is considered negligible. Thus, in the following, we adopt $\mu \approx \mu _0$ , i.e.

(5.60) \begin{gather} \mu \approx \mathfrak{A}\, \delta (\tau (q) - \tau ). \end{gather}

One can also show (Appendix C) that, at $R \gg 1$ , one can approximate $\mathfrak{A}$ as

(5.61) \begin{gather} \mathfrak{A} \approx \textrm {Ai}_\varepsilon (H_{\mathsf{y}}({\mathsf{y}}) - h) \end{gather}

with $\varepsilon = - V^3 p''$ . In fact, even when (5.43) is not satisfied, yet all higher-order derivatives of $p(q)$ are small enough, $\mathfrak{A}$ can still be approximated with a delta function, $\mathfrak{A} \approx \delta (H_{\mathsf{y}}({\mathsf{y}}) - h)$ , which may be sufficiently accurate for practical purposes. In this case, $\mathfrak{A}$ can be safely replaced with $\textrm {Ai}_\varepsilon$ because the latter can be approximated with a delta function as well. Hence, we will assume (5.61) either literally, when (5.43) is satisfied, or in the sense that both $\mathfrak{A}$ and $\textrm {Ai}_\varepsilon$ are close enough to $\delta (H_{\mathsf{y}}({\mathsf{y}}) - h)$ . Then, (5.60) becomes

(5.62) \begin{gather} \mu \approx \textrm {Ai}_\varepsilon (H_{\mathsf{y}}({\mathsf{y}}) - h)\, \delta (\tau (q) - \tau ). \end{gather}

Because ${\mathsf{y}}$ and ${\mathsf{z}}$ are connected by a linear transformation, which preserves $\mu$ , one can also consider (5.62) as the symplectic measure for the whole transformation ${\mathsf{r}} \leftrightarrow {\mathsf{z}}$ :

(5.63) \begin{gather} \mu = \textrm {Ai}_\varepsilon (H_{\mathsf{z}}({\mathsf{z}}) - h)\,\delta (\tau ({\mathsf{z}}) - \tau ). \end{gather}

At $\varepsilon \to 0$ , when the transformation ${\mathsf{z}} \mapsto {\mathsf{r}}$ is linear,Footnote ³⁴ (5.63) yields $\mu ({\mathsf{r}}, {\mathsf{z}}) = \delta (H_{\mathsf{z}}({\mathsf{z}}) - h)\,\delta (\tau ({\mathsf{z}}) - \tau )$ , in agreement with (3.57). Application of (5.63) in the general case requires that one expresses $\varepsilon$ as a function of ${\mathsf{z}}$ (as opposed to $q$ ). Since our calculation is valid only for small $q$ , multiple coordinate charts would have to be introduced and properly merged for this, which is inconvenient. However, as long as $\varepsilon$ is smooth (or indistinguishable from zero), one can use (5.16) to extrapolate it to all ${\mathsf{z}}$ via

(5.64) \begin{gather} \varepsilon = ({\mathsf{J}} \partial _{\mathsf{z}} H_{\mathsf{z}}) \boldsymbol{\cdot} (\partial _{\mathsf{z}}^2 H_{\mathsf{ z}}) ({\mathsf{J}} \partial _{\mathsf{z}} H_{\mathsf{z}}). \end{gather}

This ensures that $\varepsilon$ is close to the desired value $\varepsilon = - V^3 p''$ in the ray vicinity, while far from the ray, our result is not expected to be applicable in any case. Then, one can readily transform the symbols of operators with

(5.65) \begin{gather} O_{\mathsf{z}}({\mathsf{z}}) \approx \int \textrm {Ai}_\varepsilon (H_{\mathsf{z}}({\mathsf{z}}) - \tilde {h})\,O_{\mathsf{r}}(\tau ({\mathsf{z}}), \tilde {h})\, {\mathrm{d}} \tilde {h}, \end{gather}

i.e. $O_{\mathsf{z}}({\mathsf{z}})$ is the Airy transform (Widder Reference Widder1979) of $O_{\mathsf{r}}$ . Equation (5.65) is one of the key new results of this paper. As will be seen in the following, it is an important theorem that makes MGO a practical theory.

Assuming the symplectic scale of $O_{\mathsf{r}}$ is of order $R$ , one can take this integral using (C2). This leads to $O_{\mathsf{z}}({\mathsf{z}}) \approx O_{\mathsf{r}}(\tau ({\mathsf{z}}), H_{\mathsf{z}}({\mathsf{ z}})) + {\mathcal{O}}(\varepsilon \partial _h^3 O_{\mathsf{r}})$ . Based on (5.64), the second term on the right-hand side can be estimated as $\varepsilon \partial _h^3 O_{\mathsf{r}} \sim O_{\mathsf{r}}/R^4$ . Thus,

(5.66) \begin{gather} O_{\mathsf{z}}({\mathsf{z}}) = O_{\mathsf{r}}({\mathsf{r}}({\mathsf{z}})) + {\mathcal{O}}(\epsilon ^2). \end{gather}

The importance of this result is in that the approximate equality $O_{\mathsf{z}}({\mathsf{z}}) \approx O_{\mathsf{r}}({\mathsf{r}}({\mathsf{z}}))$ is satisfied with accuracy not $\mathcal{O}(\epsilon )$ but ${\mathcal{O}}(\epsilon ^2)$ , which will be negligible for our purposes (§ 6). This is another new key result of this paper.Footnote ³⁵ Keep in mind, though, that (5.66) is invalid for small-scale symbols, such as the delta-shaped symbols discussed in § 6.2. In that case, the more general (5.65) should be used instead.

6.1. Dynamics in the $\tau$ -space

Armed with (5.66), we are now ready to formulate the MGO envelope equation. Let us return to our original wave equation (4.16), which accounts for dissipation, and rewrite it in the $\tau$ -representation:

(6.1) \begin{gather} \big(\hat {H}_{\hat {\tau }} + \mathrm{i}\hat {\varGamma }_{\hat {\tau }}\big)\psi _{\hat {\tau }} = 0. \end{gather}

(As earlier, we assume that $\psi _{\hat {\tau }}$ is a scalar function, but see § 7.2 for a more general case.) Since $\hat {H}_{\hat {\tau }} = \hat {h}_{\hat {\tau }} = -\mathrm{i} \partial _\tau$ , this leads to

(6.2) \begin{gather} \big({-} \partial _\tau + \hat {\varGamma }_{\hat {\tau }}\big)\psi _{\hat {\tau }} = 0. \end{gather}

Assuming that $\varGamma _{\mathsf{r}} = {\mathcal{O}}(\epsilon )$ , the operator $\hat {\varGamma }_{\hat {\tau }}$ can be approximated simply by its symbol on the reference ray, $\varGamma _{\mathsf{r}}(\tau )$ . Then, (6.2) becomes

(6.3) \begin{gather} \partial _\tau \psi _{\hat {\tau }} = \varGamma _{\mathsf{r}}(\tau ) \psi _{\hat {\tau }}. \end{gather}

Assuming that $\varGamma _{\mathsf{r}}$ is sufficiently smooth, (5.66) predicts that $\varGamma _{\mathsf{r}}({\mathsf{r}}) \approx \varGamma _{\mathsf{z}}({\mathsf{z}}({\mathsf{r}}))$ , which is a known function when the dispersion operator is given. Then, (6.3) can be readily integrated,

(6.4) \begin{gather} \psi _{\hat {\tau }} = a \exp \left (\int {\mathrm{d}}\tau \varGamma _{\mathsf{r}}(\tau )\right )\!, \quad a = \text{const.}, \end{gather}

and, unlike in the GO (4.12), this solution exhibits no singularities. Also notably, the ray equations in the ${\mathsf{ r}}$ -representation are particularly simple:

(6.5) \begin{gather} \dot {\tau } = 1, \quad \dot {h} = 0. \end{gather}

A generalisation to vector waves is discussed in § 7. One can also generalise the MGO envelope equation to quasioptical beams that experience transverse diffraction. This can be done similarly to how quasioptical equations were derived by Dodin et al. (Reference Dodin, Ruiz, Yanagihara, Zhou and Kubo2019) and Yanagihara, Dodin & Kubo (Reference Yanagihara, Dodin and Kubo2019a , Reference Yanagihara, Dodin and Kubob ) for conventional GO. However, presenting these calculations in detail would require introducing additional cumbersome machinery that is not specific to MGO, so we postpone discussing this subject until further publications.

6.2. Mapping solutions to the $x$ -space

Mapping solutions to the $x$ -space is generally the most challenging part in MGO applications (Donnelly et al. Reference Donnelly, Lopez and Dodin2021; Højlund Marholt et al. Reference Højlund Marholt, Senstius and Nielsen2024; Lopez, Højlund & Senstius Reference Lopez, Højlund and Senstius2024). However, knowing the field at all $x$ is usually unnecessary. One might not even need to leave the $\tau$ -space except for initialising the field and also calculating the final state, and those mappings may be easy to do.Footnote ³⁶ In other cases, one can proceed as follows.

6.2.1. General case

To calculate how the wave evolves in the original $x$ -space, one needs to remap the $\tau$ -representation to the $x$ -representation:

(6.6)

Since this leads to

(6.7)

where $\vartheta (q) = \int ^{q}_0 p(\tilde {q})\,{\mathrm{d}} \tilde {q}$ . Since consists of a phase-space shift (§ 3.5.1) and an LST (§ 3.5.3), the integrand in (6.7) is easy to write explicitly in any given case.

How to approximate such an integral depends on the problem and can be done using various techniques, with accuracy requirements determined by specific applications. For example, stationary-phase or steepest-descent methods might suffice, or one might be able to reduce (6.7) to a known representation of a special function. In particular, to the extent that (6.7) can be expressed through the Airy function (or, equivalently, Bessel functions of order $1/3$ ), it is similar to previous results (Langer Reference Langer1937; Chester, Friedman & Ursell Reference Chester, Friedman and Ursell1957). Similar integrals have been studied also in the context of semiclassical calculations in atomic and molecular physics; for example, see Kay (Reference Kay1994b ).

Remember also that the above-mentioned result is derived under the assumption that only small values of $q$ contribute to the integral (6.7). This allows one to derive local approximations $\psi _{\hat {x}}$ using different definitions of $(q, p)$ for different regions of interest. A global representation is not necessarily desirable for practical applications. When it is, though, it can be constructed in one of two ways.

One way to do this is to numerically find the type-1 generating function $\varTheta (x, \tau )$ of the transformation $(x, k) \mapsto (\tau , h)$ using

(6.8) \begin{gather} k = \partial _x\varTheta (x, \tau ), \quad h = -\partial _\tau \varTheta (x, \tau ), \quad {\mathrm{d}}\varTheta = k\,{\mathrm{d}} x - h\,{\mathrm{d}}\tau , \end{gather}

and then invoke (6.6) in combination with the approximate formula (4.40). For the specific coordinates that we assume here, the latter takes the form

(6.9)

This formula is also useful for determining the quantisation condition for periodic orbits, as discussed in Appendix D.

Alternatively, one may notice that $\psi _{\hat {x}}$ per se is rarely needed in applications. Instead, of interest are bilinear functionals of $\psi _{\hat {x}}$ and $\psi _{\hat {x}}^*$ , such as the energy density or the dissipated power. Any such functional can be expressed through the field’s Wigner function.Footnote ³⁷ By (2.43) and (6.4), the Wigner function of $\psi _{\hat {\tau }}$ in the ${\mathsf{r}}$ -representation is

(6.10) \begin{gather} W_{\mathsf{r}}(\tau , h) \approx W_0(\tau ) \delta (h), \quad W_0(\tau ) = |a|^2\exp \left (2 \int {\mathrm{d}}\tau \varGamma _{\mathsf{r}}(\tau )\right )\!. \end{gather}

By (3.41), its ${\mathsf{ z}}$ -representation is $W_{\mathsf{z}}({\mathsf{z}}) = \int {\mathrm{d}} \tau \,{\mathrm{d}} h\,W_{\mathsf{r}}(\tau , h)\,\mu ({\mathsf{r}}, {\mathsf{z}})$ . This leads to a remarkably concise expression for $W_{\mathsf{z}}$ :

(6.11) \begin{gather} W_{\mathsf{z}}({\mathsf{z}}) = W_0(\tau ({\mathsf{z}}))\, \textrm {Ai}_\varepsilon (H_{\mathsf{z}}({\mathsf{z}})). \end{gather}

Figure 3.

Approximation (6.16), with (5.64) for $\varepsilon$ , for the Wigner functions $W_{\mathsf{z}}$ of fields satisfying (4.16): (a) $H_{\mathsf{z}}(x, k) = (x^2 + k^2 - 11)/2$ , which corresponds to a QHO with $n = 5$ ; (b) $H_{\mathsf{z}}(x, k) = - 2 - \cos (x/10) + k^2/2$ . The intensity denotes the magnitude of $W_{\mathsf{z}}$ (arbitrary units). The white areas in panel (b) correspond to $\varepsilon \approx 0$ . The approximation is accurate near the ray trajectories given by $H_{\mathsf{z}}(x, k) = 0$ (dashed). In panel (b), there are two such trajectories, one for each sign of $k$ . Strictly speaking, the Wigner function shown corresponds to a standing wave, to which both these branches contribute equally. Otherwise, although still applicable locally, (6.16) cannot be continued across the region where the contributions of the two disconnected branches are comparable (i.e. at $|k| \lesssim 1$ in this case).

6.2.2. Eigenwaves

Let us now consider an important special case when dissipation is negligible, at least locally ( $\hat {\varGamma } = 0$ ), and a wave satisfies $\hat {H}_{\hat {\tau }}\psi _{\hat {\tau }} \approx 0$ . Such a wave can be understood as an eigenstate of $\hat {H}$ corresponding to the zero eigenvalue, $h = 0$ ; i.e. , where $a_h \equiv \sqrt {2\pi }a$ is some constant amplitude. Then, using (2.17), one obtains

(6.12)

which, of course, can also be obtained directly from (6.3). For $\psi _{\hat {x}}$ , (6.7) readily yields

(6.13)

In the special case when $q \equiv x$ , one has so the standard WKB solution is obtained:

(6.14) \begin{gather} \psi _{\hat {x}}(x) = \psi _{\hat {\tau }}\,\frac {\mathrm{e}^{\mathrm{i} \vartheta (x)}}{\sqrt {V(x)}}. \end{gather}

More generally, (6.13) represents, expectedly, the same WKB solution in the tangent space, which is then mapped to the $x$ -space using a linear MT. One can approach the integral (6.13) in the same way as discussed in § 6.2.1.

Alternatively, one can use (6.10) and (6.11), which, in this case, become

(6.15) \begin{align}&\quad W_{\mathsf{r}}(\tau , h) = W_0 \delta (h),\\[-10pt]\nonumber \end{align}

(6.16) \begin{align}& W_{\mathsf{z}}({\mathsf{z}}) \approx W_0 \textrm {Ai}_\varepsilon (H_{\mathsf{z}}({\mathsf{z}})), \end{align}

with $W_0 = |a|^2 = \text{const}$ . Examples of such $W_{\mathsf{z}}({\mathsf{z}})$ are presented in figure 3 for rays in a quadratic potential, which corresponds to a QHO, and a sinusoidal potential. Similar patterns are also seen in direct wave simulations; for example, see Zhu & Dodin (Reference Zhu and Dodin2021, figure 3a2) or Weinbub & Ferry (Reference Weinbub and Ferry2018, figure 9). Note that these complicated patterns are merely a result of mapping the Wigner function to the original coordinates ${\mathsf{z}}$ , as opposed to the natural coordinates ${\mathsf{r}}$ , where the Wigner functions have a simple form (6.15). Also note that (6.15) is more precise than (6.16) in that it does not rely on approximations of the MT and is valid at all ${\mathsf{z}}$ . In contrast, (6.16) is accurate only in the vicinity of the ray trajectories, where $H_{\mathsf{z}}$ is small. This is illustrated in figure 4, which shows a comparison of (6.16) with the exact analytical result for a QHO, $H_{\mathsf{z}}(x, k) = (x^2 + k^2)/2 - (n + 1/2)$ (Mostowski & Pietraszewicz Reference Mostowski and Pietraszewicz2021):

(6.17) \begin{gather} W_{\mathsf{z}}(x, k) = \frac {({-}1)^n}{\pi }\, e^{-\xi ^2} L_n(2\xi ), \quad \xi \doteq x^2 + k^2. \end{gather}

Here, $n$ is a non-negative integer that corresponds to the state number, $L_n$ is $n$ th Laguerre polynomial, and the standard normalisation $a_h = 1$ is assumed, which corresponds to $W_0 = (2\pi )^{-1}$ .

Figure 4.

Wigner functions $W(x) \equiv W_{\mathsf{z}}(x, k = 0)$ of a QHO, $H_{\mathsf{z}} = (x^2 + k^2)/2 - (n + 1/2)$ : (a) $n = 5$ ; (b) $n = 10$ . Red – exact analytical result (6.17); blue – approximation (6.16), with $W_0 = (2\pi )^{-1}$ . The latter is valid only close to the ray (dashed curve in figure 3 a), which, at $k = 0$ considered here, corresponds to the reflection points $x = \sqrt {2n + 1}$ .

To the extent that $\textrm {Ai}_\varepsilon$ can be replaced with the delta function, (6.16) yields

(6.18) \begin{gather} |\psi _{\hat {x}} (x)|^2 = \int {\mathrm{d}} k\, W_{\mathsf{z}}(x, k) \approx W_0 \int {\mathrm{d}} k\,\delta (H_{\mathsf{z}}(x, k)) = \frac {W_0}{|v_{\text{g}}|}, \end{gather}

where $v_{\text{g}} \doteq \partial _k H_{\mathsf{z}}$ and the derivative is evaluated on the dispersion curve at a given $x$ and the corresponding $k$ on the dispersion curve. This coincides with the result that flows from the action conservation for stationary waves (see, e.g. Dodin Reference Dodin2022, § 7.3.2),

(6.19) \begin{gather} \partial _x(v_{\text{g}} |\psi _{\hat {x}}|^2) = 0, \end{gather}

and $|\psi _{\hat {x}}|^2$ depends on $x$ only due to inhomogeneity of $v_{\text{g}}$ , i.e. slowly. If a wave experiences reflection, one might want to account for the multiple branches, i.e. two or more values of $k$ at which the delta-function argument in (6.18) turns to zero at given $x$ .Footnote ³⁸ Then,

(6.20) \begin{gather} |\psi _{\hat {x}} (x)|^2 = \sum _i \frac {W_0}{|v_{\text{g}}|_i}, \end{gather}

where the summation is taken over said branches. Like in the previous case, the dependence on $x$ (due to inhomogeneity of $|v_{\text{g}}|_i$ ) is slow in this case. This is because the assumed delta approximation for $W_{\mathsf{z}}$ ignores the interference between multiple branches; i.e. (6.20) can be understood as $|\psi _{\hat {x}} (x)|^2$ averaged over their relative phases.

One can also integrate the actual Wigner function (6.16) without resorting to the delta approximation, at least as long as the parabolic approximation (5.29) holds. For example, in the simplest case when $(q, p) \equiv (x, k)$ , one readily obtains that

(6.21) \begin{gather} |\psi _{\hat {x}} (x)|^2 = W_0 \int {\mathrm{d}} k\,\textrm {Ai}_\varepsilon \left (v k + \frac {v x^2}{2R}\right ) = \frac {W_0}{v}, \end{gather}

which is in agreement with (6.18). More generally, when $(q, p)$ are connected with $(x, k)$ via (3.50) (we assume ${\mathsf{z}}_0 = 0$ for brevity), one can write

(6.22) \begin{gather} H_{\mathsf{z}}({\mathsf{z}}) = \underbrace {- \frac {v}{B}\left (x + \frac {R D^2}{2B}\right )}_{U(x)} + \frac {v B^2}{2R}\bigg (\underbrace {\vphantom {\bigg (}k + \frac {A x}{B} + \frac {R D}{B^2}}_K\bigg )^2, \end{gather}

where we used (3.51). Then, using (Vallée et al. Reference Vallée, Soares and de Izarra1997)

(6.23) \begin{gather} \int _{-\infty }^\infty \textrm {Ai}(X + P^2)\,{\mathrm{d}} P = 2^{2/3} \pi \, \textrm {Ai}^2(X/2^{2/3}) \end{gather}

and (5.24), one finds that

(6.24) \begin{align} |\psi _{\hat {x}} (x)|^2 & = W_0 \int {\mathrm{d}} K\, \textrm {Ai}_\varepsilon \left (U(x) + \frac {v B^2}{2R}\,K^2\right ) \notag \\ & = W_0\,\frac {\pi v}{2^{1/3}|B|}\, \textrm {Ai}_\varepsilon ^2\left (\frac {U(x)}{2^{2/3}}\right )\!. \end{align}

In particular, the case when $B = 1$ and $D = 0$ corresponds to having a reflection point at $x = 0$ , when the ${\mathsf{y}}$ -space is rotated by $\pi /2$ relative to the ${\mathsf{z}}$ -space. Then, one obtains

(6.25) \begin{gather} |\psi _{\hat {x}} (x)|^2 = 2\pi W_0\,\frac {2^{2/3}}{|R|^{1/3}}\, \textrm {Ai}^2({-}\sigma _\varepsilon |2R|^{1/3}x), \end{gather}

which corresponds to approximating (4.16) with the (appropriately rescaled) Airy equation near the reflection point. Although this is a commonly known result, a comparison of (6.25) with the exact solutions for a QHO is shown in figure 5 for completeness.

Figure 5.

Normalised eigenstates $|\psi |^2 \equiv |\psi _{\hat {x}}(x)|^2$ of a QHO for two sample states: (a) $n = 15$ and (b) $n = 50$ . Blue – exact analytical result; orange – approximation (6.25), with $W_0 = (2\pi )^{-1}$ and $R = \sqrt {2n + 1}$ . Unlike in (6.21)–(6.26), which assume ${\mathsf{z}}_0 = 0$ for simplicity, the $x$ -axis origin here is at the centre of the potential well, as usual; i.e. to reproduce these plots, one should replace $x$ with $x - R$ in (6.25).

Let us also consider the opposite limit, that is, $B \to 0$ and $D \to 1$ . One might expect to recover (6.21) in this case, yet (6.24) leads to a different result. Assuming the notation $\xi \doteq |R|^{2/3}/(2^{1/3}|B|) \gg 1$ , it can be written as follows:

(6.26) \begin{gather} |\psi _{\hat {x}} (x)|^2 \to \frac {2W_0}{v}\,2\pi \xi \, \textrm {Ai}^2({-}\xi ^2) \approx \frac {2W_0}{v}\left (1 + \sin \frac {4\xi ^{3}}{3} \right )\!. \end{gather}

This is understood as follows. Notice that now the problem is posed such that $|\psi _{\hat {x}} (x)|^2$ automatically includes the contributions from two branches, that is, the two roots $k(x)$ of the parabolic Hamiltonian (6.22). At $B \to 0$ , $D \to 1$ and $x \to 0$ , one of these roots approaches zero, while the other one tends to infinity, $k \approx$ $-2R/B^2 \to -\infty$ . To the extent that the interference with this second branch can be ignored, one can average out the second term in (6.26) to zero and obtain a result consistent with (6.20), $|\psi _{\hat {x}} (x)|^2 \to 2W_0/v$ . Likewise, if the second branch is ignored entirely,Footnote ³⁹ then

(6.27) \begin{gather} |\psi _{\hat {x}} (x)|^2 \to W_0/v, \end{gather}

as expected. This can equivalently be done by replacing the function $\textrm {Ai}^2$ in (6.26) with $|\textrm {Ai} + \mathrm{i}\,\textrm {Gi}|^2/4$ , where $\textrm {Gi}$ is the Scorer function, to isolate the contribution from a single branch (Lopez Reference Lopez2024). The asymptotic form of $|\textrm {Ai} + \mathrm{i}\,\textrm {Gi}|$ for large $\xi$ then gives (6.27).