Pseudo-differential representation of the metaplectic transform and its application to fast algorithms

TL;DR

This paper introduces a pseudo-differential form of the metaplectic transform, enabling efficient numerical computation of the transform for small-angle rotations and proposing an algorithm with complexity scaling as O(K N^3 N_p).

Contribution

The paper derives a pseudo-differential representation of the metaplectic transform and develops a fast algorithm for its numerical implementation, especially for small-angle rotations.

Findings

Asymptotic differential representations of the MT for small angles.

An efficient algorithm with complexity O(K N^3 N_p) for larger rotations.

Numerical implementation and stability analysis of the algorithm.

Abstract

The metaplectic transform (MT), also known as the linear canonical transform, is a unitary integral mapping which is widely used in signal processing and can be viewed as a generalization of the Fourier transform. For a given function $\psi$ on an $N$-dimensional continuous space $\textbf{q}$, the MT of $\psi$ is parameterized by a rotation (or more generally, a linear symplectic transformation) of the $2N$-dimensional phase space $(\textbf{q},\textbf{p})$, where $\textbf{p}$ is the wavevector space dual to $\textbf{q}$. Here, we derive a pseudo-differential form of the MT. For small-angle rotations, or near-identity transformations of the phase space, it readily yields asymptotic \textit{differential} representations of the MT, which are easy to compute numerically. Rotations by larger angles are implemented as successive applications of $K \gg 1$ small-angle MTs. The algorithm complexity scales as $O(K N^3 N_p)$, where $N_p$ is the number of grid points. We present a numerical implementation of this algorithm and discuss how to mitigate the associated numerical instabilities.

Equations

$$i{\partial }_t|{\psi }_t⟩=\hat{H}|{\psi }_t⟩,\hat{H}\doteq ({\hat{p}}^2+{\hat{q}}^2)/2.$$

$${\hat{M}}_t=\exp \left(-i\hat{H}t\right).$$

$${\partial }_t({\hat{M}}_t^{†}\hat{q}{\hat{M}}_t)$$

$${\partial }_t({\hat{M}}_t^{†}\hat{p}{\hat{M}}_t)$$

$$\hat{Q}=\cos (t)\hat{q}+\sin (t)\hat{p},\hat{P}=-\sin (t)\hat{q}+\cos (t)\hat{p},$$

$$\hat{Q}\doteq {\hat{M}}_t^{†}\hat{q}{\hat{M}}_t,\hat{P}\doteq {\hat{M}}_t^{†}\hat{p}{\hat{M}}_t.$$

$$\int \mathrm{d}x|q(x)⟩⟨q(x)|=\int \mathrm{d}x|Q(x)⟩⟨Q(x)|=\hat{1}$$

$$\Psi (y)$$

$$=\int \mathrm{d}x⟨q(y)|{\hat{M}}_t|q(x)⟩{\psi }_0(x).$$

$$\hat{H}|n⟩=(n+1/2)|n⟩,$$

$${\hat{M}}_t=\exp \left(-it/2\right)\sum _{n=0}^{\infty }\exp \left(-int\right)|n⟩⟨n|.$$

$$\hat{\mathrm{Q}}\doteq {\hat{M}}^{†}\hat{\mathrm{q}}\hat{M},\hat{\mathrm{P}}\doteq {\hat{M}}^{†}\hat{\mathrm{p}}\hat{M},$$

$$\left(\begin{array}{c}\hat{\mathrm{Q}}\\ \hat{\mathrm{P}}\end{array}\right)=\mathsf{S}\left(\begin{array}{c}\hat{\mathrm{q}}\\ \hat{\mathrm{p}}\end{array}\right),\mathsf{S}=\left(\begin{array}{cc}\mathsf{A}& \mathsf{B}\\ \mathsf{C}& \mathsf{D}\end{array}\right),$$

$$\mathsf{S}\left(\begin{array}{cc}{0}_N& {\mathsf{I}}_N\\ -{\mathsf{I}}_N& 0_N\end{array}\right){\mathsf{S}}^{⊺}=\left(\begin{array}{cc}{0}_N& {\mathsf{I}}_N\\ -{\mathsf{I}}_N& 0_N\end{array}\right),$$

$$\mathsf{A}{\mathsf{D}}^{⊺}-\mathsf{B}{\mathsf{C}}^{⊺}$$

$${\mathsf{A}}^{⊺}\mathsf{D}-{\mathsf{C}}^{⊺}\mathsf{B}$$

$$\mathsf{A}{\mathsf{B}}^{⊺}-\mathsf{B}{\mathsf{A}}^{⊺}$$

$${\mathsf{B}}^{⊺}\mathsf{D}-{\mathsf{D}}^{⊺}\mathsf{B}$$

$${\mathsf{C}}^{⊺}\mathsf{A}-{\mathsf{A}}^{⊺}\mathsf{C}$$

$$\mathsf{D}{\mathsf{C}}^{⊺}-\mathsf{C}{\mathsf{D}}^{⊺}$$

$$\Psi (\mathrm{y})$$

$$U(\mathrm{y},\mathrm{x})$$

$$⟨\mathrm{Q}(\mathrm{y})|\hat{\mathrm{p}}|\mathrm{q}(\mathrm{x})⟩={\left[⟨\mathrm{q}(\mathrm{x})|\hat{\mathrm{p}}|\mathrm{Q}(\mathrm{y})⟩\right]}^{\ast }=i{\partial }_{\mathrm{x}}U(\mathrm{y},\mathrm{x})$$

$$\mathrm{y}U(\mathrm{y},\mathrm{x})=\left(\mathsf{A}\mathrm{x}+i\mathsf{B}{\partial }_{\mathrm{x}}\right)U(\mathrm{y},\mathrm{x}),$$

$$U(\mathrm{y},\mathrm{x})=f(\mathrm{y})e^{\frac{i}{2}{\mathrm{x}}^{⊺}{\mathsf{B}}^{-1}\mathsf{A}\mathrm{x}-i{\mathrm{x}}^{⊺}{\mathsf{B}}^{-1}\mathrm{y}}.$$

$${\partial }_{\mathrm{y}}U(\mathrm{y},\mathrm{x})=\left(i\mathsf{C}\mathrm{x}-\mathsf{D}\partial \mathrm{x}\right)U(\mathrm{y},\mathrm{x}).$$

$$f(\mathrm{y})=\alpha e^{\frac{i}{2}{\mathrm{y}}^{⊺}\mathsf{D}{\mathsf{B}}^{-1}\mathrm{y}}.$$

$$\Psi (\mathrm{Q})$$

$$\times \int \mathrm{d}\mathrm{q}{e}^{\frac{i}{2}{\mathrm{q}}^{⊺}{\mathsf{B}}^{-1}\mathsf{A}\mathrm{q}-i{\mathrm{q}}^{⊺}{\mathsf{B}}^{-1}\mathrm{Q}}\psi (\mathrm{q}),$$

$$\Psi (\mathrm{Q})$$

$$\times \int \mathrm{d}\mathrm{u}{e}^{i{\mathrm{u}}^{⊺}{\Lambda }^{-1}\mathrm{u}}\psi ({\mathsf{A}}^{-1}\mathrm{Q}+\mathrm{u}),$$

$$\mathsf{G}\doteq \mathsf{C}{\mathsf{A}}^{-1}/2,\Lambda \doteq 2{\mathsf{A}}^{-1}\mathsf{B}.$$

$$\psi \left(\frac{Q}{A}+u\right)=\sum _{n=0}^{\infty }\frac{u^n}{n!}{\psi }^{(n)}\left(\frac{Q}{A}\right),$$

$${\int }_{-\infty }^{\infty }\mathrm{d}u{e}^{i{\Lambda }^{-1}{u}^2}\psi \left(\frac{Q}{A}+u\right)$$

$$\sim \sum _{n=0}^{\infty }\frac{1}{n!}{\psi }^{(n)}\left(\frac{Q}{A}\right){\int }_{-\infty }^{\infty }\mathrm{d}u{u}^n{e}^{i{\Lambda }^{-1}{u}^2}.$$

$$\sum _{n=0}^{\infty }\frac{1}{n!}{\psi }^{(n)}\left(\frac{Q}{A}\right){\int }_{-\infty }^{\infty }\mathrm{d}u{u}^n{e}^{i{\Lambda }^{-1}{u}^2}$$

Full text

I Introduction

Suppose a signal described by a square-integrable function $\psi$ of some continuous coordinate ${\boldsymbol{\rm q}}$. Like in quantum mechanics, one can introduce a ‘state vector’ $|\psi\rangle$ such that $\psi$ be the projection of $|\psi\rangle$ onto the coordinate axis. Correspondingly, $\psi$’s Fourier image $\smash{\widetilde{\psi}}$ can be viewed as the projection of $|\psi\rangle$ onto the wavevector axis ${\boldsymbol{\rm p}}$, or equivalently, onto the coordinate axis obtained via rotation of the original phase space $({\boldsymbol{\rm q}},{\boldsymbol{\rm p}})$ by $\pi/2$. But one can also introduce rotations by different angles or, most generally, linear symplectic transformations of the original phase space. Suppose a phase space $({\boldsymbol{\rm Q}},{\boldsymbol{\rm P}})$ obtained via such transformation of $({\boldsymbol{\rm q}},{\boldsymbol{\rm p}})$. One can then obtain $\Psi$, the projection of $|\psi\rangle$ onto the new coordinate space ${\boldsymbol{\rm Q}}$, and relate it to the original projection $\psi$ by a linear unitary mapping. This mapping is called the metaplectic transform (MT) Littlejohn (1986); de Gosson (2006)111It is also sometimes called the linear canonical transform.. It subsumes the Fourier transform as a special case and represents one of the pillars of modern phase space analysis used in many applications Tracy and Kaufman (1993); Tracy et al. (2007); Gopinathan et al. (2008); Camara et al. (2011); Bazarov (2012); Child (2014).

To accommodate these applications, a number of numerical algorithms have been proposed which efficiently compute the MT on both 1-dimensional (1-D) and 2-D configuration spaces Ozaktas et al. (1996); Hennelly and Sheridan (2005); Healy and Sheridan (2010); Koc et al. (2010); Ding et al. (2012); Pei and Huang (2016); Sun and Li (2018). Many of them are reviewed in LABEL:Healy18. Despite this multitude, however, there also exist applications for which suitable MT algorithms have yet to be designed. In particular, consider the modeling of electromagnetic waves in media with slowly-varying parameters. Such waves are usually described by the equations of geometrical optics Tracy et al. (2014), but this approach fails near reflection points, where the local wavenumber goes to zero and its derivative becomes singular. The MT provides a means to reinstate geometrical optics near reflection points, because a simple rotation of the phase space can make the wavenumber nonzero again Littlejohn (1985) (see Sec. V). It is convenient to perform such rotations consecutively along the ray trajectory at small angles; the corresponding MTs will be near-identity. Since the existing algorithms treat the MT as an integral transform, they are not optimal for computing the MT in this limit. A differential representation would be advantageous but remains to be developed.

Here, we propose an algorithm which closes this gap, as it is specifically tailored to computing near-identity MTs. We start by deriving a general pseudo-differential form of the MT. For small-angle phase-space rotations, or more generally, for any near-identity symplectic transformations of the phase space, this readily yields asymptotic differential representations of the MT, which are easy to compute numerically. Rotations by larger angles can be implemented as successive applications of $K\gg 1$ small-angle MTs. We show that the algorithm complexity scales as $O(KN^{3}N_{p})$, where $N$ is the dimension of the configuration space and $N_{p}$ is the number of grid points. This means that our algorithm allows computing the MT in linear time, which is a faster scaling than other published MT algorithms Healy (2018), albeit with a potentially-large prefactor. We then assess the stability of our algorithm, discuss ways to optimize its performance, and present a numerical implementation.

The paper is organized as follows. In Sec. II, we introduce the MT in a familiar setting of elementary quantum mechanics. In Sec. III, we derive the pseudo-differential representation of the MT from its integral representation, and we also discuss its possible truncations. In Sec. IV, we describe how the near-identity MT can be used in an iterative algorithm to perform cumulative MTs which are not near-identity. We also discuss the computational complexity and stability of such an algorithm, and demonstrate how it can be used to simulate quadratic Hamiltonian systems. In Sec. V, we outline briefly how our new algorithm can feature in a ray-tracing code to resolve caustics, using Airy’s equation as an example. In Sec. VI, we present our main conclusions. Auxiliary calculations are presented in appendices.

II Metaplectic transforms and their integral representations

II.1 Special case: a quantum harmonic oscillator and its propagator as an MT

To better understand what the MT is, let us first consider an elementary problem from quantum mechanics, namely, the quantum harmonic oscillator (QHO). The QHO is described by the Schrödinger equation222In the following, we adopt the operator notation that is standard in quantum-mechanical literature and also in optics LABEL:Stoler81. Bold font denotes vectors, sans serif font denotes matrices, and $\doteq$ denotes definitions.

[TABLE]

Equation (1) has the solution $|\psi_{t}\rangle=\smash{\hat{M}}_{t}|\psi_{0}\rangle$, where $|\psi_{0}\rangle$ is an initial wavefunction and the propagator $\smash{\hat{M}}_{t}$ is a unitary operator given by

[TABLE]

An interesting property of $\smash{\hat{M}}_{t}$ is revealed by switching from the Schrödinger representation to the Heisenberg representation, in which the wavefunction is fixed but $\smash{\hat{q}}$ and $\smash{\hat{p}}$ evolve in time as governed by Shankar (1994)

[TABLE]

The coordinate and momentum operators of the QHO are seen to satisfy the same Hamilton’s equations that describe a classical harmonic oscillator Goldstein et al. (2002). The solution to Eqs. (3) is therefore

[TABLE]

where we introduced

[TABLE]

Equations (4) can be considered as a mapping $(\smash{\hat{q}},\smash{\hat{p}})\mapsto(\smash{\hat{Q}},\smash{\hat{P}})$ which is a phase-space rotation by angle $t$. The unitary propagator $\smash{\hat{M}}_{t}$ that effects this rotation is called a metaplectic operator333Here, $\smash{\hat{M}}_{t}$ also acts as the fractional Fourier transform operator, up to a phase..

The metaplectic operator also induces a mapping between the projections of $|\psi_{0}\rangle$ onto the original coordinate axis $q$ and onto the new axis $Q$. The former is defined as $\psi(x)\doteq\langle q(x)|\psi_{0}\rangle$, where $|q(x)\rangle$ is the eigenvector of $\smash{\hat{q}}$ corresponding to the eigenvalue $x$. Likewise, the projection onto $Q$ is $\Psi(y)\doteq\langle Q(y)|\psi_{0}\rangle$, where $|Q(y)\rangle$ is the eigenvector of $\smash{\hat{Q}}$ corresponding to the eigenvalue $y$. We assume the usual normalization, $\langle q(x)|q(y)\rangle=\langle Q(x)|Q(y)\rangle=\delta(x-y)$, so

[TABLE]

and $|Q(x)\rangle=\smash{\hat{M}}_{t}^{\dagger}|q(x)\rangle$. (Here $\smash{\hat{\mathbb{1}}}$ is a unit operator.) Then,

[TABLE]

Note that the right-hand side of (7) is the same as $\psi_{t}(y)\doteq\langle q(y)|\smash{\hat{M}}_{t}|\psi_{0}\rangle$, because in our example $\smash{\hat{M}}_{t}$ is the propagator. Hence, for the QHO considered here, the MT can be equivalently understood as the evolution of the wavefunction in the Schrödinger representation, $|\psi_{0}\rangle\mapsto|\psi_{t}\rangle$, or as the evolution of the projection basis in the Heisenberg representation, $\langle q(y)|\mapsto\langle Q(y)|$.

Finally, let us notice the following. As is well-known, the eigenvalues of the QHO Hamiltonian are Shankar (1994)

[TABLE]

with $n$ an integer and $|n\rangle$ the $n$-th eigenstate of $\smash{\hat{H}}$; hence, the specific MT considered in Eq. (2) can also be represented as

[TABLE]

A notable aspect of this formula is that it takes not one but two rotation periods ($t=4\pi$) for $\smash{\hat{M}}_{t}$ to return to its original value $\smash{\hat{M}}_{0}=\smash{\hat{\mathbb{1}}}$. More generally, $\smash{\hat{M}}_{2\pi n}=\smash{\hat{\mathbb{1}}}$ for even $n$ yet $\smash{\hat{M}}_{2\pi n}=-\smash{\hat{\mathbb{1}}}$ for odd $n$. Hence, the same identity transformation on phase space [governed by Eq. (4)] can be effected by two distinct metaplectic operators, $\pm\smash{\hat{\mathbb{1}}}$. This double-valuedness also holds for arbitrary rotation angles, and is in fact a general property of the MT. This is illustrated by analogy with the behavior of the complex function $f(z)\doteq\sqrt{z}$ in Fig. 1.

II.2 General definition of the MT

A more general definition of the MT is as follows. Let $\smash{\hat{{\boldsymbol{\rm q}}}}$ and $\smash{\hat{{\boldsymbol{\rm p}}}}$ be respectively the $N$-dimensional coordinate and momentum operators. Consider

[TABLE]

where $\smash{\hat{M}}$ is a unitary operator such that

[TABLE]

and $\mathsf{S}$ is real and symplectic. The latter means that

[TABLE]

which implies (cf. Appendix A) Luneburg’s relations Luneburg (1964)

[TABLE]

where $\mathsf{0}_{N}$ and $\mathsf{I}_{N}$ denote respectively the $N\times N$ null and identity matrices444Note that at $N=1$, Eqs. (13c)-(13f) are satisfied automatically, and Eqs. (13a) and (13b) are equivalent to $\det\mathsf{S}=1$; hence, a $2\times 2$ matrix is symplectic if and only if it has unit determinant.. Then, $\smash{\hat{M}}$ is called the metaplectic operator corresponding to the chosen $\mathsf{S}$.

Like in the previous section, we now define the MT as the mapping between a given function $\psi$ on the coordinate space associated with $\smash{\hat{{\boldsymbol{\rm q}}}}$ and the projection of the corresponding state vector $|\psi\rangle$ on the coordinate space associated with $\smash{\hat{{\boldsymbol{\rm Q}}}}$. Again, this leads to555Analogous to the Schrödinger and Heisenberg representations of time evolution, there exists in the general case a distinction between whether $\smash{\hat{M}}$ transforms the wavefunction (‘active’ representation) or transforms the projection basis (‘passive’ representation). In our discussion, we assume the passive representation.

[TABLE]

To calculate $U$, let us consider the top row of Eq. (11), $\smash{\hat{{\boldsymbol{\rm Q}}}}=\mathsf{A}\smash{\hat{{\boldsymbol{\rm q}}}}+\mathsf{B}\smash{\hat{{\boldsymbol{\rm p}}}}$, and apply $\langle{\boldsymbol{\rm Q}}({\boldsymbol{\rm y}})|$ from the left and $|{\boldsymbol{\rm q}}({\boldsymbol{\rm x}})\rangle$ from the right. Using the eigenvalue relations along with

[TABLE]

leads to a differential equation Littlejohn (1986); Moshinsky and Quesne (1971)

[TABLE]

which can be solved to yield

[TABLE]

Doing the same with the bottom row of Eq. (11) leads to

[TABLE]

Using Eqs. (13), (17), and (18) determines $f({\boldsymbol{\rm y}})$ up to a multiplicative constant:

[TABLE]

Normalization determines the constant $\alpha$ up to a phase. The phase requires more involved analysis to determine, and the result is not unique: there exist two possible phases which differ by $\pi$. This ambiguity is required to ensure that the metaplectic operators form a group, but results in a one-to-two correspondence between the symplectic and the metaplectic groups Littlejohn (1986). In other words, changing the overall sign of a metaplectic operator does not change the resulting phase-space transformation, which Eqs. (4) and (9) demonstrate for the QHO example. As discussed in Sec. II.1 (and also related to the general Bohr-Sommerfeld rule Shankar (1994)), the sign ambiguity becomes important when one considers a family of transformations parameterized by some path variable $t$. A closed trajectory in the space of symplectic matrices, $\mathsf{S}_{t}$, results in a closed trajectory in the space of metaplectic operators only for even winding numbers. In contrast, for odd winding numbers $\smash{\hat{M}}_{t}$ changes sign, just like the function $f(z)\doteq\sqrt{z}$ changes sign each time $z$ encircles the origin in the complex plane Littlejohn (1986) (see Fig. 1).

Including the phase and sign ambiguity, the final result for the transformation is Collins (1970); Moshinsky and Quesne (1971); Littlejohn (1986)

[TABLE]

where $\mathsf{B}^{-1}\mathsf{A}$ and $\mathsf{D}\mathsf{B}^{-1}$ are symmetric due to Eqs. (13c) and (13d). Equation (20) defines $\Psi({\boldsymbol{\rm Q}})$ as the MT image of $\psi({\boldsymbol{\rm q}})$. In writing Eq. (20), we have dropped the ${\boldsymbol{\rm x}}$ and ${\boldsymbol{\rm y}}$ notation in favor of ${\boldsymbol{\rm q}}$ and ${\boldsymbol{\rm Q}}$, as there is no longer any risk of ambiguity, and our branch cut convention restricts all complex phases to the interval $(-\pi,\pi]$.

III Pseudo-differential representation of the Metaplectic Transform

Here, we develop a pseudo-differential representation of Eq. (20). This representation is particularly useful when $\mathsf{A}^{-1}\mathsf{B}$ is small, because then the MT can be approximated by a finite-order differential transform, which is easier to evaluate than the integral transform of Eq. (20). Specifically, we proceed as follows. Using the substitution ${\boldsymbol{\rm u}}\doteq{\boldsymbol{\rm q}}-\mathsf{A}^{-1}{\boldsymbol{\rm Q}}$, Eq. (20) can be re-written as

[TABLE]

where we have defined the matrices

[TABLE]

Notably, both $\mathsf{G}$ and $\Lambda$ are symmetric per Eqs. (13c) and (13e). In the following, we shall assume that $\Lambda$ is small. This assumption is not strictly necessary, since the final result is convergent for all values of $\|\Lambda\|$ and thereby possesses a natural analytic continuation; however, it aids intuition in the forthcoming derivation.

III.1 1-D case

Let us first consider the $1$-D case ($N=1$) for simplicity. Since $\Lambda^{-1}$ is assumed large, only small values of $u$ will contribute to the integral of Eq. (21). Therefore, we can expand the function $\psi\left(Q/A+u\right)$ around the point $u=0$ as

[TABLE]

where $\psi^{(n)}(Q/A)$ is the $n$-th derivative of $\psi(q)$ evaluated at $q=Q/A$. (Here, we assume that $\psi$ is smooth, but we shall revisit this assumption below.) Hence,

[TABLE]

By parity, all integrals with odd powers of $u$ are identically zero, so the sum can be written solely in terms of even powers as

[TABLE]

Let us introduce a dummy multiplicative variable $s$, which will eventually be taken to unity. Then, since

[TABLE]

we obtain

[TABLE]

where the first line invokes Leibniz’s rule, the final equality follows from the binomial theorem Olver et al. (2010), and $\Gamma(z)$ is the gamma function Olver et al. (2010). By combining Eqs. (24), (25), and (27), we obtain the asymptotic representation

[TABLE]

Finally, using well-known properties of the gamma function yields the pseudo-differential representation of the MT in $1$-D:

[TABLE]

or symbolically,

[TABLE]

We can also express Eq. (29b) in an equivalent vector form $|\Psi\rangle=\smash{\hat{M}}|\psi\rangle$, where $\smash{\hat{M}}$ is the manifestly-unitary MT operator given as

[TABLE]

and $\smash{\hat{D}}_{A}$ is the inverse dilation operator defined via its effect in the spatial representation, $\langle q|\smash{\hat{D}}_{A}|\psi\rangle\doteq\psi(q/A)$.

We call Eqs. (29) the 1-D pseudo-differential metaplectic transform (PMT). Although the above derivation assumes smooth $\psi$, the final result can be understood more generally, which is why the asymptotic relation has been replaced with an exact equality. As shown in Appendix B, the operator (30) has exactly the same kernel as the original integral MT (20) and exists on the space of all functions which have a well-defined Fourier transform; i.e., smoothness of $\psi$ is not required. In this sense, Eq. (29b) should be understood not as a symbolic representation of the series (29a) (whose convergence depends on details of $\psi$) but rather as a symbolic representation of the integral MT (20). This new representation is advantageous in that it is compact, and facilitates asymptotic expansions of the MT to any order of $\Lambda$.

Let us also discuss the case when $\Lambda$ is small and $\psi$ is smooth enough so Eq. (29a) can be approximated with a truncated series. We define the $m$-th order near-identity metaplectic transform (NIMT) as the truncation of Eq. (29a) that neglects all terms with $n>m$. This nomenclature is chosen because up to a phase, the limit $B\to 0$ reduces Eqs. (29) to a scaled-identity operation. Also, to be connected with the identity, we explicitly choose the overall $+$ sign when performing NIMT truncations. Decreasing $m$ will increase the locality of the truncated transformation, because the necessary stencil width to compute the $m$-th order NIMT will decrease. This enables the $m$-th order NIMT to be performed pointwise, as the transformed function evaluated at some point $Q=Q_{0}$ depends only on the original function and its first $2m$ derivatives evaluated at the corresponding point $q=q_{0}(Q_{0})$.

When the order is not specified, the ‘NIMT’ refers solely to the first-order NIMT,

[TABLE]

as it is the lowest-order truncation that remains practical. (The truncation at $m=0$ is too simplified to yield an accurate representation of the MT, regardless the smoothness of $\psi$.) We shall make use of Eq. (31) in Sec. IV.

III.2 N-D case

The generalization from 1-D to the arbitrary $N$-D case is straightforward. We consider again the integral of Eq. (21). Since $\Lambda$ is a symmetric matrix, by the spectral theorem it can be diagonalized. Let us enumerate with subscripts $j\in\{1,\ldots,N\}$ vector components with respect to the diagonalizing basis of $\Lambda$. Then,

[TABLE]

where $\lambda_{j}$ is the $j$-th eigenvalue of $\Lambda$. As before, $\psi(\mathsf{A}^{-1}{\boldsymbol{\rm Q}}+{\boldsymbol{\rm u}})$ is expanded around ${\boldsymbol{\rm u}}=0$. In multiple dimensions, this expansion is written as

[TABLE]

with the shorthand notation

[TABLE]

denoting the derivatives of $\psi$ along the eigenvectors of $\mathsf{A}^{-1}\mathsf{B}$. The integral therefore becomes

[TABLE]

where once again, the summation has been restricted to even integers by parity considerations.

The remaining integrals are of the same form as those in Eq. (27). Hence, the $N$-D PMT is ultimately obtained as

[TABLE]

or symbolically,

[TABLE]

where the notation $\Lambda\text{\large:}\nabla\nabla$ denotes the double dot product between $\Lambda$ and the Hessian operator $\nabla\nabla$; i.e., $\Lambda:\nabla\nabla=\Lambda_{ab}(\partial_{x_{a}})(\partial_{x_{b}})$ summed over common indices. In this case, the equivalent operator representation for Eq. (36b) uses

[TABLE]

and $\langle{\boldsymbol{\rm q}}|\smash{\hat{D}}_{\mathsf{A}}|\psi\rangle\doteq\psi(\mathsf{A}^{-1}{\boldsymbol{\rm q}})$. Retaining only the terms corresponding to $\sum_{j=1}^{N}n_{j}=0$ and $\sum_{j=1}^{N}n_{j}=1$ in (36), the $N$-D NIMT is

[TABLE]

where the term in brackets is evaluated at ${\boldsymbol{\rm q}}=\mathsf{A}^{-1}{\boldsymbol{\rm Q}}$, and the overall $+$ sign is assumed, as in Sec. III.1. Since matrix operations can be computationally expensive when $N$ is large, Appendix C provides some low-order approximations for $\det{\mathsf{A}}$, $\mathsf{A}^{-1}$, $\mathsf{A}^{-1}\mathsf{B}$, and $\mathsf{C}\mathsf{A}^{-1}$ for use when $\mathsf{S}$ is near-identity. We also provide auxiliary calculations when $\psi({\boldsymbol{\rm q}})$ is eikonal in Appendix D.

IV Finite transformations via iterated Near-Identity transformations

The pseudo-differential representation of the MT naturally gives rise to an iterative algorithm: successive applications of the NIMT can compute a finite transformation from a sequence of near-identity transformations. To see this, consider the MT of a function $\psi$ that results from the symplectic matrix $\mathsf{S}$, which may be the result of a single optical operation or a cascade of operations. As the symplectic group is topologically connected, it is always possible to find a smooth trajectory of symplectic matrices $\mathsf{S}_{t}$ with parameterization $t$ such that $\mathsf{S}_{0}=\mathsf{I}_{2N}$ and $\mathsf{S}_{T}=\mathsf{S}$ at some final $T$.

Let us discretize $\mathsf{S}_{t}$ with a uniform step size $\Delta t\doteq T/K$ such that $\forall~{}k\in\{1,2,\ldots,K\}$, the matrix $\mathsf{S}_{(k-1)\Delta t}^{-1}\mathsf{S}_{k\Delta t}$ is near-identity. Then, since

[TABLE]

one can compute the MT associated with $\mathsf{S}$ by iteratively applying the NIMT: first the NIMT associated with $\mathsf{S}_{\Delta t}$ (which is near-identity by definition), next the NIMT associated with $\mathsf{S}_{\Delta t}^{-1}\,\mathsf{S}_{2\Delta t}$, and so forth until finally, the NIMT associated with $\mathsf{S}_{T-\Delta t}^{-1}\,\mathsf{S}_{T}$. Hence,

[TABLE]

where $\mathcal{N}_{\mathsf{S}}$ is the NIMT associated with symplectic matrix $\mathsf{S}$.

Note that the discretization of $\mathsf{S}_{t}$ by itself does not incur any errors, so the accuracy of Eq. (40) depends solely on the truncation order of the NIMT. Another advantage of this approach is that the algorithm is independent of the dimensionality. One only needs to adjust the size of $\mathsf{S}$ when changing from, say, a $1$-D application to a $3$-D application. This is not true for other numerical MT algorithms in the literature, which can only handle up to $2$-D and are explicitly different depending on whether $\mathsf{S}$ is ‘separable’ or ‘nonseparable’ Koc et al. (2010); Ding et al. (2012); Pei and Huang (2016). Such restrictions do not arise with the iterated NIMT.

IV.1 Computational efficiency

Let us estimate the computational efficiency of the iterated NIMT. We should first emphasize that although the NIMT appears to require interpolation, this is not strictly necessary. Suppose that $\psi({\boldsymbol{\rm q}})$ is only known on a discrete set of points $\{{\boldsymbol{\rm q}}_{k}\}$. The discretization of $\psi({\boldsymbol{\rm q}})$ can be used to inform the discretization of $\Psi({\boldsymbol{\rm Q}})$ by evaluating the NIMT only at the corresponding points $\{{\boldsymbol{\rm Q}}_{k}\doteq\mathsf{A}{\boldsymbol{\rm q}}_{k}\}$. No interpolation is required, unless, one needs to evaluate $\Psi({\boldsymbol{\rm Q}})$ off-grid. In that case, either the discrete set $\{\psi({\boldsymbol{\rm q}}_{k})\}$ must be interpolated and transformed, or the discrete set $\{\Psi({\boldsymbol{\rm Q}}_{k})\}$ must be interpolated. For this reason, and because interpolation efficiency is highly implementation-specific, we do not account for interpolation in our runtime estimate.

From Eq. (31), evaluating $\Psi(Q)$ at $N_{p}$ discrete points using the $1$-D NIMT requires only $O(N_{p})$ floating-point operations (FLOPs). For the $N$-D case, this estimate becomes $O(N^{3}N_{p})$, since each evaluation includes a matrix multiplication Trefethen and Bau, III (1997). Thus, the NIMT always scales linearly with the number of sample points, independent of dimensionality. The iterated NIMT remains ‘fast’ with respect to the number of sample points, since the FLOP count scales as $O(KN^{3}N_{p})$, with $K$ the number of iterations. The linear scaling is faster than many of the other MT algorithms found in the literature Healy (2018), which scale as $O(N_{p}\log N_{p})$.

IV.2 Computational stability

Although the iterated NIMT scales faster than other published MT algorithms, it may not be as stable. Intuitively, one would expect that refining the discretization of $\mathsf{S}_{t}$ would increase the accuracy of the iterated NIMT, since the magnitude of $\|\mathsf{A}^{-1}_{j}\mathsf{B}_{j}\|$ for each successive $j$-th application of the NIMT would decrease. As the magnitude of $\|\mathsf{A}^{-1}_{j}\mathsf{B}_{j}\|$ decreases, however, the number of iterations required to generate a fixed final transformation increases. Careful analysis is needed to determine if the truncation errors of the iterated NIMT accumulate coherently, which we accomplish by estimating the parameter regimes in which the iterated NIMT is non-unitary. For simplicity, the forthcoming analysis is restricted to $1$-D.

Let us consider how the PMT and the iterated NIMT transform the single-parameter family of exponential functions $\psi_{\kappa}(q)\doteq e^{\kappa q}$, with $\kappa$ being complex. Generally speaking, we define an MT algorithm as stable, or non-magnifying, if the norms of the transformed function $\Psi_{\kappa}(Q)$ and the original function $\psi_{\kappa}(q)$ satisfy $\|\Psi_{\kappa}(Q)\|\leq\|\psi_{\kappa}(q)\|$; conversely, we define an MT algorithm as unstable, or magnifying, if $\|\Psi_{\kappa}(Q)\|>\|\psi_{\kappa}(q)\|$. Unitarity corresponds to a strict equality. The ratio $\|\Psi_{\kappa}(Q)\|/\|\psi_{\kappa}(q)\|$ is referred to as the magnification factor. Additionally, we define an MT algorithm as either $L^{2}$-stable or, respectively, $L^{2}$-unitary if the algorithm is stable or unitary along the entire imaginary $\kappa$ axis. This is because any $L^{2}$ function can be expanded into Fourier modes; thus, an $L^{2}$-unitary MT algorithm will be exactly unitary for any $L^{2}$ function. In our analysis, we shall only consider the class of function norms where $\|e^{ig(Q)}f(Q/A)\|=\sqrt{A}\,\|f(q)\|$ for $g(Q)$ real, an example of which being the $L^{2}$ norm.

Since $\psi_{\kappa}^{\prime}(q)=\kappa\psi_{\kappa}(q)$, the PMT of $\psi_{\kappa}(q)$ is

[TABLE]

where $\mathcal{P}_{\mathsf{S}}$ is the PMT for symplectic matrix $\mathsf{S}$. Let us define the rescaled variable $w\doteq\kappa\sqrt{B/A}$. Then, the PMT is stable when

[TABLE]

This region of the complex $w$ plane is shown in Fig. 2. The PMT is stable within the first and third quadrants of the complex plane, and is unitary along the real and imaginary $w$ axes. Hence, the PMT is $L^{2}$-unitary. Interestingly, the PMT is not unitary on its entire domain. This is because the domain of the PMT includes both square-integrable functions and functions where the integral of Eq. (20) does not converge, such as $e^{q}$. The cost of this expanded domain is the loss of global unitarity, albeit for functions whose $L^{2}$ norms are undefined.

We proceed to analyze the NIMT. Applied once, the NIMT of $\psi_{\kappa}(q)$ is

[TABLE]

Reintroducing $w$, the NIMT is stable where

[TABLE]

which is shown alongside the stability region of the PMT in Fig. 2. This region makes up only a small subset of the stability region of the PMT.

Notably, the NIMT is no longer $L^{2}$-stable; as such, square-integrable functions will be magnified. There are three ways to minimize the magnification: (i) reduce the step size to $B/A\lesssim 1/2\kappa_{i}$, with $\kappa_{i}$ the largest Fourier mode number; (ii) apply a low-pass filter to remove fast-growing Fourier modes; or (iii) increase the truncation order. However, we show shortly that increasing the truncation order of the NIMT increases its vulnerability to numerical noise, so only (i) and (ii) are recommended.

Let us now assess how subsequent iterations of the NIMT affect its stability. We first observe that each iteration of the NIMT adds the overall phase $Q^{2}C/2A$ that contributes to the derivatives of subsequent NIMT iterations. This sequence will very quickly become unwieldy as the iteration number increases. To achieve an analytical estimate of the iterated NIMT stability, we shall therefore neglect the contributions of the phase to all derivatives. This is consistent with the near-identity limit, where $C/A$ is vanishingly small.

In this approximation, the norm of the $K$-iterated NIMT is

[TABLE]

where for $n=1$, we define $\prod_{j=1}^{n-1}A^{2}_{j}=1$. When the iteration is uniform, i.e., $A_{n}=A$ and $B_{n}=B$, Eq. (45) simplifies to

[TABLE]

where we have reintroduced $w\doteq\kappa\sqrt{B/A}$. Hence, the $K$-iterated NIMT is stable where

[TABLE]

where $(a;q)_{K}\doteq\prod_{n=0}^{K-1}\left(1-aq^{n}\right)$ is the $q$-Pochhammer symbol Olver et al. (2010), i.e., the $q$-analog of the rising factorial.

Figure 3 shows the stability region at four different iteration numbers: $K=2$, $K=5$, $K=10$, and $K=20$. As Eq. (47) indicates, the stability of the iterated NIMT now explicitly depends on $A$, so each subplot of Fig. 3 includes stability diagrams for $A=0.9$, $A=1$, $A=1.1$, and $A=2$. These values were chosen to emphasize the near-identity behavior of the iterated NIMT, when $A\approx 1$. There are two notable observations. First, the stability region for $A=1$ is independent of $K$. For other values of $A$, the stability region changes significantly with $K$, decreasing for $A<1$ and increasing for $A>1$. Second, the sensitivity of the iterated NIMT increases with $K$, as seen by considering the rate at which the $A=1.1$ and $A=0.9$ contours separate.

Consequently, a step size $B/A$ that is initially stable, but with $A<1$, will become quickly and increasingly unstable as the NIMT is iterated. This introduces an interesting tradeoff consideration when computing a finite transformation: is it better to use a coarse discretization with a large step size but few iterations, or a fine discretization with a small step size but many iterations? The answer depends largely on implementation specifics; we find in the following subsection that a fine discretization is preferable for our chosen example, but this is not necessarily indicative of a general principle.

Although we shall not dwell much on implementation details, we must make one cautionary remark regarding the finite-difference scheme used to discretize the NIMT. Because discrete differentiation is poorly conditioned, any noise in the original function $\psi(q)$ will be magnified when its derivatives are computed. Since derivatives are computed with each iteration of the NIMT, this noise will grow geometrically. We call this instability the d-instability (with ‘d’ standing for discretization). As shown in Fig. 4, it is particularly disastrous for iterated NIMTs with large truncation order.

A basic description of the d-instability is afforded by the transformation of a constant function. Suppose one attempts to transform a function that is identically zero everywhere except at a single grid point, where the function is erroneously non-zero by some unspecified noise source. When the grid spacing $h$ is uniform, the growth rate of the d-instability, $\gamma$, can be estimated analytically. Let $\Delta_{k}$ be a $k$-th order finite-difference matrix such that $h^{-k}\Delta_{k}{\boldsymbol{\rm f}}$ equals the $k$-th discrete derivative of ${\boldsymbol{\rm f}}$. In this specific test problem, any non-zero norm is due to noise; hence, the error of the $K$-th iterated, $m$-th order NIMT is bounded with the triangle inequality as

[TABLE]

where $\boldsymbol{\psi}$ and $\boldsymbol{\Psi}$ are the discretized versions of $\psi(q)$ and $\Psi(Q)$ respectively, and $\|\Delta_{2n}\|$ is the subordinate matrix norm of $\Delta_{2n}$. There is a freedom to choose the norm with which Eq. (48) is evaluated; we choose the $\infty$-norm, denoted $\|\ldots\|_{\infty}$, as it yields the readily-evaluated matrix row norm as its subordinate Trefethen and Bau, III (1997). Considering only the leading order in $1/h$, $\gamma$ is estimated as

[TABLE]

Equation (49) has been purposefully separated: for increasing $m$, the factor $\|\Delta_{2m}\|/h^{2m}$ increases while the factor $|B/2A|^{m}/m!$ decreases. In fact, as defined, $\gamma\to 0$ as $m\to\infty$ for any reasonable class of $\Delta_{2m}$; this does not mean the d-instability disappears for high truncation orders, but rather that the d-instability is not dominated by the leading order in $1/h$ when $m$ is large. Instead, a subset of intermediate-order terms dominate, which are not included in Eq. (49). For central finite-difference schemes with homogeneous boundary conditions, $\|\Delta_{2n}\|_{\infty}=2^{2n}$ Abramowitz and Stegun (1970), and the growth rate is uniformly estimated to be

[TABLE]

where $\Gamma\left(s,x\right)$ is the incomplete gamma function Olver et al. (2010). Notably, $\Gamma\left(s,x\right)\to\Gamma\left(s\right)$ as $s\to\infty$. Equation (49) is sufficient for the error estimation of low truncation order schemes; for large $m$, however, Eq. (50) should be used instead.

Thus far, our discussion of the d-instability has been contingent on a maliciously designed initial condition. Such a specific state will not likely arise in practical applications; nevertheless, local d-instabilities can certainly arise. For example, consider the NIMT of a function $\psi(q)$ that asymptotes to [math] at the domain edge. Near the domain edge, $\psi(q)$ is nearly constant, but a source of error, interpolation or otherwise, will inevitably cause at least one data point to deviate. The local d-instability will then grow rapidly, and will propagate inward from the domain edge until the transformed function is entirely dominated by noise. Since the d-instability growth rate scales with truncation order, using low $m$ schemes will minimize its deleterious effects. Marginally smoothing the input data before taking derivatives will also delay its onset.

IV.3 Numerical example: Time evolution of the QHO

To demonstrate the iterated NIMT, let us consider once again the time evolution of the $1$-D QHO, introduced in Sec. II. There, the symplectic matrix $\mathsf{S}$ can be expressed as

[TABLE]

where we have added the index $t$ to emphasize the dependence on time. This matrix can be represented as $\mathsf{S}_{t}=(\mathsf{S}_{\Delta t})^{K}$, where $\mathsf{S}_{\Delta t}$ evolves the system by $\Delta t\ll 1$ and $K=t/\Delta t$. From Eq. (11), the scalar functions $A_{\Delta t}$, $B_{\Delta t}$, $C_{\Delta t}$, and $D_{\Delta t}$ to be used in the iterated NIMT are

[TABLE]

For visualization, it is useful to introduce the Wigner function that corresponds to $\psi(q)$, defined as Wigner (1932):

[TABLE]

As shown in Refs. de Gosson (2006); Littlejohn (1986); Lohmann (1993), the Wigner function of the metaplectic image $\Psi$ is simply the Wigner function of the original $\psi$ correspondingly rotated. This is also readily understood from the physical meaning of $W_{\psi}$. Specifically, if $\psi$ is a wave field, then $W_{\psi}$ can be interpreted as the phase-space quasiprobability distribution function of the wave quanta. The prefix ‘quasi’ marks the fact that $W_{\psi}$ is not positive-definite unless averaged over a phase space volume of size $\Delta q\,\Delta p\gtrsim 2\pi$ Cartwright (1976); O’Connell and Wigner (1981); nonetheless, $W_{\psi}$ is always real by definition, even for complex $\psi$.

For our example, we consider the time evolution of three initial states: (i) a chirped Gaussian profile, $\psi(q)=e^{(i-1)q^{2}}$, which is relevant for bit-flip operations in chirp-modulated communication Kaminsky and Simanjuntak (2005), (ii) a squeezed coherent state, $\psi(q)=e^{-(q+3)^{2}+2iq}$, which is relevant for high-sensitivity detectors Xiao et al. (1987), and (iii) the QHO eigenstate corresponding to $n=3$, namely, $\psi(q)=(48\sqrt{\pi})^{-1/2}e^{-q^{2}/2}H_{3}(q)$, where $H_{n}(q)$ is the $n$-th degree Hermite polynomial. For these choices, the exact metaplectic image of $\psi(q)$ can be found explicitly from Eq. (20), which facilitates benchmarking of our algorithm. They are given respectively by

[TABLE]

The overall sign is chosen based on the winding number of $\mathsf{S}_{t}$ as discussed in Fig. 1: an odd winding number corresponds to the $-$ sign, while an even winding number corresponds to the $+$ sign. Each of these three example functions is evolved in time using the iterated NIMT with a uniform step size of $\Delta t=\pi/2000$, and the resulting Wigner functions are shown in Fig. 5. Here, $\Psi_{t}(Q)$ is discretized on an equally-spaced grid ranging from $[-10,10]$, and in the final example, the second-order NIMT was used in place of the first-order NIMT.

As the NIMT is sequentially applied, Fig. 5 shows that the resultant Wigner functions indeed rotate in phase space as expected. (In the third example, $W_{\psi}$ is rotationally-symmetric, so it is preserved by the MT.) This shows that the iterated NIMT can indeed perform finite transformations with high accuracy. For computing the Fourier transform, which corresponds to $t=\pi/2$, the iterated NIMT is robust to changes in the step size; discretizing the trajectory into $10^{2}$, $10^{3}$, and $10^{4}$ steps all yielded well-behaved solutions. The same is not true for changes in grid resolution, nor in changes of truncation order. Indeed, $\Psi_{t}(Q)$ quickly succumbed to amplified noise when (i) a Chebyshev-spaced grid was used in place of the equally-spaced grid, and (ii) the truncation order was increased beyond third-order.

Recall from Fig. 3 that the iterated NIMT is typically a magnifying transformation whose magnification factor depends in a complicated manner on both the path discretization and the input function. For our chosen examples, the magnification is reduced by refining the discretization of the path $\mathsf{S}_{t}$ (Fig. 6). When a step size of $\pi/500$ is used, $\Psi_{t}(Q)$ quickly disrupts and becomes dominated by noise. However, refining the discretization by a factor of $10$ avoids the numerical instability and leads to a well-behaved solution.

We reiterate that the magnification of the NIMT is not reduced for every input function by refining the discretization; a rigorous profiling should be performed to determine how the magnification scales with path discretization when using the iterated NIMT in a new application. Alternatively, since the magnification scales with Fourier mode number, occasionally smoothing the signal between NIMT iterations will suppress high-frequency growth. This approach is shown in the final column of Fig. 6, where a third-degree Savitzky–Golay filter Savitzky and Golay (1964) with a window size of $5$ is applied every $50$ iterations.

V Phase-space rotation for cutoff removal

In addition to the time integration of quadratic Hamiltonian systems, the NIMT is also naturally suited for modeling caustics that arise in geometrical optics near cutoffs. As motivated in the introduction, such caustics can be resolved by rotating the phase space using metaplectic operators. Although this idea of using phase-space rotations to avoid caustics is not entirely new Tracy et al. (2014); Littlejohn (1985), it is not yet a common practice, and therefore merits brief discussion.

Consider a wave field incident on an isolated cutoff in a $1$-D inhomogeneous medium. As is well-known, the corresponding wave field $\psi$ near the cutoff is approximately described by Airy’s equation Tracy et al. (2014)

[TABLE]

with the cutoff located at $a$. Applying the Wentzel–Kramers–Brillouin (WKB) approximation to Eq. (55) yields the dispersion surface $p(q)=\sqrt{a-q}$ on which the wave ‘quanta’ is asymptotically confined, as well as the divergent wave envelope $\phi(q)\sim(a-q)^{-1/4}$. Thus, the caustic at $q=a$ manifests as a singularity in the WKB envelope, as illustrated in the first column of Fig. 7.

Let us now rotate the phase space using the MT corresponding to Eq. (51) as

[TABLE]

with $t$ now specifying the (negative) rotation angle rather than time. Then, Eq. (55) becomes

[TABLE]

where $\Psi(Q)$ is the metaplectic image of $\psi(q)$. Applying the WKB approximation to Eq. (57) yields

[TABLE]

with $\alpha_{\pm}$ constants determined by boundary conditions, which should be matched on either side of the caustic separately due to Stokes phenomenon Heading (1962).

In Fig. 7, the WKB result is compared to the exact result, which can be computed via Eq. (20) as

[TABLE]

where $\mathrm{Ai}$ is the Airy function Olver et al. (2010). As the phase space is rotated, the caustic moves steadily towards increasing $Q$. At a rotation angle of $\pi/2$, Eq. (57) becomes

[TABLE]

In this case, the caustic disappears entirely and the WKB approximation, obtained from the ($+$) solution to Eqs. (58) as

[TABLE]

becomes exact. Importantly, the WKB approximation to Eq. (60) holds at any $a$, even though for $a\geq 0$ there are values of $Q$ at which the wavenumber $P$ approaches zero. (For example, see Fig. 8 for the case $a=0$.) This is to be expected because: (i) $a$ can be removed from Eq. (55) by a simple variable transformation, and hence has no fundamental meaning, and (ii) it is $\mathrm{d}P/\mathrm{d}Q$ that determines the validity of geometrical optics, not the value of $P$ per se.

For other equations, a single MT is not sufficient to reinstate geometrical optics for the entire field. However, multiple MTs applied sequentially can. Specifically, a phase space can be continually rotated using the NIMT such that $\mathrm{d}P/\mathrm{d}Q$ always remains finite (Fig. 9). In that frame, the WKB approximation will hold indefinitely and there will be no caustics.

For example, consider the wave equation

[TABLE]

with $a$ constant. The WKB approximation is applicable only far enough from the cutoffs located at $q=\pm a$, and there is no single MT for which the image of Eq. (62) will be free of cutoffs. However, note that Eq. (62) is the time-independent limit of the QHO (1), whose solutions are eigenfunctions of the phase-space rotation operator (Sec. II). Therefore, in the appropriately-rotating frame which maintains the wavenumber constant (say, $P=a$), the WKB approximation holds perfectly. More general wave equations can be handled in a similar manner, but require introducing additional machinery. Hence, the corresponding discussion is postponed until future publications.

VI Conclusion

In this work, we derive a pseudo-differential representation of the MT in arbitrary dimensions. This is a general result which can be useful for both analytical and numerical applications. An important example is the simulation of a wavepacket evolving in a quadratic potential, whose propagator is a metaplectic operator. Evolving the system by $\Delta t$ would invoke an MT that is near-identity, which is not a common consideration in MT-algorithm design. In contrast, the pseudo-differential representation that we propose here readily displays the simplicity of the MT in the near-identity limit, suggestive of a new algorithm.

Specifically, in the near-identity limit the pseudo-differential series can be accurately truncated; the correspondingly finite stencil width then enables local, pointwise transformations. This is useful when transforming ‘incomplete’ functions, e.g., signals measured over finite intervals; it also leads naturally to a linear time algorithm called the NIMT. When applied once, the NIMT performs a fast, near-identity transformation; when iterated, the NIMT can perform an arbitrary MT by synthesizing a series of near-identity transformations. With a computational efficiency of $O(KN^{3}N_{p})$, the NIMT is potentially faster than existing MT algorithms, which often scale as $O(N_{p}\log N_{p})$ from their similarity with the fast Fourier transform. Moreover, unlike these other algorithms, the NIMT is the same algorithm regardless the number of dimensions and the structure of $\mathsf{S}$. Hence, the NIMT is flexible in its application, and should thereby complement the existing collection of MT algorithms.

We assess the stability of the iterated NIMT and identify two dominant instabilities: the loss of unitarity via truncation error (magnification), and the poor conditioning of discrete derivatives (d-instability). One might expect the NIMT magnification to be suppressed by reducing the transformation ‘step size’, i.e., its deviation from identity, or by increasing the number of terms retained; however, this is not true. Reducing the step size increases the number of iterations needed to perform a finite transformation, and it is not clear whether this tradeoff is beneficial in the general case. Increasing the truncation order indeed decreases the NIMT magnification, but also increases its susceptibility to the d-instability. The most robust avenue to NIMT stability therefore appears to be the combined use of a low-order truncation with occasional smoothing, which we demonstrate in a numerical example.

Acknowledgements

This work was supported by the U.S. DOE through Contract No. DE-AC02-09CH11466.

Appendix A Derivation of Eqs. (13)

Here, we present the derivation of Eqs. (13) from Eq. (12), which is known Tracy and Kaufman (1993) but included here for completeness. Consider

[TABLE]

Since $\mathsf{J}\mathsf{J}=-\mathsf{I}_{2N}$, Eq. (12) implies that

[TABLE]

and also that $\mathsf{S}^{-1}$ is symplectic, i.e.,

[TABLE]

Using

[TABLE]

together with Eq. (12) leads to Eqs. (13a), (13c), and (13f). Likewise, since

[TABLE]

Eq. (65) readily yields Eqs. (13b), (13d), and (13e).

Appendix B Deriving the metaplectic transform from its pseudo-differential representation

Here, we show the pseudo-differential representation (36b) leads to the original integral representation (20) regardless the size of $\|\Lambda\|$. This proves that the PMT is in fact exact, even though it was originally derived in Sec. III using an expansion in $\|\Lambda\|$.

As a starting point, let us rewrite Eq. (36b) as

[TABLE]

where we have replaced $\nabla$ with $\partial_{{\boldsymbol{\rm q}}^{\prime}}$ to avoid ambiguities. We introduce the Fourier representation of $\psi({\boldsymbol{\rm q}})$ as

[TABLE]

which, when substituted into Eq. (68), yields

[TABLE]

The Gaussian integral can be performed explicitly,

[TABLE]

with the branch cut chosen to restrict all complex phases to the interval $(-\pi,\pi]$, and with $\Delta{\boldsymbol{\rm q}}\doteq{\boldsymbol{\rm q}}^{\prime}-{\boldsymbol{\rm q}}$. Then, performing the trivial integration over $\mathrm{d}{\boldsymbol{\rm q}}^{\prime}$ yields Eq. (20).

Note that neither the smoothness nor even the differentiability of $\psi$ is invoked in the above argument; the existence of the Fourier image of $\psi$ is enough.

Appendix C Asymptotic parameterization of near-identity symplectic matrices

The $N$-D NIMT involves computing the quantities $\det{\mathsf{A}}$, $\mathsf{A}^{-1}$, $\mathsf{A}^{-1}\mathsf{B}$, and $\mathsf{C}\mathsf{A}^{-1}$. However, when $\mathsf{S}$ is near-identity, one can derive approximate asymptotic formulas for these quantities which help calculate them more efficiently. In particular, calculating the lowest-order terms does not require any explicit matrix multiplications.

Generally speaking, the near-identity behavior of a group is governed by its Lie algebra. For the group of $2N\times 2N$ real symplectic matrices, denoted $\text{Sp}(2N,\mathbb{R})$, the Lie algebra is the space of all $2N\times 2N$ real Hamiltonian matrices Dragt (2005). Note that a matrix $\mathsf{H}$ is Hamiltonian if and only if $\mathsf{J}\mathsf{H}$ is symmetric, with $\mathsf{J}$ defined in Eq. (63).

By the connectivity of $\text{Sp}(2N,\mathbb{R})$ and the polar decomposition, any symplectic matrix $\mathsf{S}$ can be parameterized as Littlejohn (1986); Dragt (1982); Hall (2015)

[TABLE]

where $\mathsf{H}_{s}$ and $\mathsf{H}_{a}$ are symmetric and antisymmetric Hamiltonian matrices, respectively. The formal parameter $\epsilon$ has been introduced to aid with ordering the forthcoming expansions when $\mathsf{H}_{s}$ and $\mathsf{H}_{a}$ are small. Note that if $\mathsf{H}$ is Hamiltonian, then $\mathsf{H}^{\intercal}$ also is; hence, $\mathsf{H}_{s}$ and $\mathsf{H}_{a}$ can be uniquely represented as

[TABLE]

for some Hamiltonian matrix $\mathsf{H}$ Dragt (1982). In this sense, $\mathsf{S}$ is parameterized by a single Hamiltonian matrix $\mathsf{H}$.

Let us consider the case when $\mathsf{S}$ is near-identity, meaning $\mathsf{H}$ is close to $\mathsf{0}_{2N}$. Expanding Eq. (72) in $\epsilon$ yields

[TABLE]

Since any Hamiltonian matrix can be decomposed as

[TABLE]

with $\mathsf{U}$ and $\mathsf{W}$ being symmetric matrices, we obtain the following expansions from Eq. (74):

[TABLE]

where

[TABLE]

One can show that

[TABLE]

satisfies both $\mathsf{A}^{-1}\mathsf{A}=\mathsf{I}_{N}$ and $\mathsf{A}\mathsf{A}^{-1}=\mathsf{I}_{N}$ to $O(\epsilon^{3})$. By direct multiplication one also obtains

[TABLE]

where the subscript s denotes the symmetric part. Notably, the expansions of both $\mathsf{A}^{-1}\mathsf{B}$ and $\mathsf{C}\mathsf{A}^{-1}$ are symmetric at each order of $\epsilon$, as required by Eqs. (13c) and (13e). Finally, let us approximate $\det{\mathsf{A}}$ as

[TABLE]

Up to the factor $\epsilon^{N}$, the right-hand side of Eq. (80a) is simply the characteristic polynomial of $-\mathsf{M}$. Using, for example, Faddeev–LeVerrier’s method leads to

[TABLE]

Appendix D Reducing the PMT to an envelope equation for eikonal functions

Often, the function $\psi({\boldsymbol{\rm q}})$ can be characterized by a rapidly-varying phase $\theta({\boldsymbol{\rm q}})$, and a complex envelope $\phi({\boldsymbol{\rm q}})$ which varies much slower than $\theta({\boldsymbol{\rm q}})$. If such a partition is defined, then we call $\psi({\boldsymbol{\rm q}})$ an eikonal function. Eikonal solutions to physical systems are frequently sought as a means to develop approximate, reduced models; an example is the WKB approximation for quantum particles Heading (1962). In reduced models, phase and envelope dynamics are typically governed by separate equations, which often makes it convenient to consider the phase and envelope as separate entities Tracy et al. (2014). Let us therefore explore how the PMT partitions eikonal functions.

Let $\psi({\boldsymbol{\rm q}})=\phi({\boldsymbol{\rm q}})e^{i\theta({\boldsymbol{\rm q}})}$, and let ${\boldsymbol{\rm k}}({\boldsymbol{\rm q}})\doteq\nabla\theta({\boldsymbol{\rm q}})$ with component functions $\{k_{j}({\boldsymbol{\rm q}})\}$. Then, by induction

[TABLE]

An analogous result is obtained in the case of mixed partial derivatives, which implies that $\nabla$ and $\widetilde{\nabla}\doteq i{\boldsymbol{\rm k}}({\boldsymbol{\rm q}})+\nabla$ have the same commutation relations among their vector components; hence, the phase function effects a formal mapping from a differential operator acting on the full function $\psi({\boldsymbol{\rm q}})$ to the differential operator acting solely on the envelope $\phi({\boldsymbol{\rm q}})$. For example, see the definition of the envelope dispersion operator in LABEL:Dodin19a.

For an eikonal function, the PMT is

[TABLE]

At least for near-identity transformations, $\Psi$ can also be cast in the eikonal form. Let $\Psi({\boldsymbol{\rm Q}})=\Phi({\boldsymbol{\rm Q}})e^{i\Theta({\boldsymbol{\rm Q}})}$, then

[TABLE]

Since $\Phi({\boldsymbol{\rm Q}})$ is generally complex, the definition of $\Theta({\boldsymbol{\rm Q}})$ is not unique, so choosing it is a matter of convenience (as long as $\Theta$ remains fast compared to $\Phi$). Here, we choose to define $\Theta({\boldsymbol{\rm Q}})$ such that it is (i) real, (ii) independent of $\phi({\boldsymbol{\rm q}})$, and (iii) simplifies the resultant expression for $\Phi({\boldsymbol{\rm Q}})$ as much as possible. Then, the first-order truncation of Eq. (84) yields the eikonal partition

[TABLE]

where $\otimes$ is the tensor product. If one prefers, additional approximations can be placed on Eqs. (85) that are consistent with the eikonal ordering ansatz, such as neglecting $\nabla\nabla\phi$ in favor of the terms involving ${\boldsymbol{\rm k}}$.

Let us also calculate the local wavevector in the new coordinates, ${\boldsymbol{\rm K}}\doteq\partial_{{\boldsymbol{\rm Q}}}\Theta$. From Eq. (85) one obtains

[TABLE]

where

[TABLE]

When ${\boldsymbol{\rm Q}}$ is obtained as ${\boldsymbol{\rm Q}}=\mathsf{A}{\boldsymbol{\rm q}}+\mathsf{B}{\boldsymbol{\rm k}}({\boldsymbol{\rm q}})$, Eq. (86) becomes

[TABLE]

Assuming that $\epsilon\doteq\|\mathsf{A}^{-1}\mathsf{B}\|$ is small, then ${\boldsymbol{\rm R}}\left[{\boldsymbol{\rm q}}+\mathsf{A}^{-1}\mathsf{B}{\boldsymbol{\rm k}}({\boldsymbol{\rm q}})\right]\approx{\boldsymbol{\rm k}}({\boldsymbol{\rm q}})+O(\epsilon^{2})$. Substituting this into Eq. (88) yields

[TABLE]

where ${\boldsymbol{\rm P}}$ is defined in Eq. (11). This shows that the transform (85) maps (${\boldsymbol{\rm q}}$,${\boldsymbol{\rm k}}({\boldsymbol{\rm q}})$) to (${\boldsymbol{\rm Q}}$,${\boldsymbol{\rm K}}({\boldsymbol{\rm Q}})$) with $O(\epsilon^{2})$ accuracy, which is consistent with the accuracy of Eqs. (85). In this sense, this transform is natural and can be useful for modeling the propagation of eikonal waves, as we shall discuss in a separate paper.