Restoring geometrical optics near caustics using sequenced metaplectic transforms

Geometrical optics (GO) is often used to model wave propagation in weakly inhomogeneous media and quantum-particle motion in the semiclassical limit. However, GO predicts spurious singularities of the wavefield near reflection points and, more generally, at caustics. We present a new formulation of GO, called metaplectic geometrical optics (MGO), that is free from these singularities and can be applied to any linear wave equation. MGO uses sequenced metaplectic transforms of the wavefield, corresponding to symplectic transformations of the ray phase space, such that caustics disappear in the new variables, and GO is reinstated. The Airy problem and the quantum harmonic oscillator are studied analytically using MGO for illustration. In both cases, the MGO solutions are remarkably close to the exact solutions and remain finite at cutoffs, unlike the usual GO solutions.


I. INTRODUCTION
The geometrical-optics (GO) approximation, sometimes called the Wentzel-Kramers-Brillouin (WKB) approximation, is commonly used to model wave propagation in weakly inhomogeneous media, a special case being the semiclassical motion of quantum particles [1][2][3][4]. It is applicable when, loosely speaking, kL ≫ 1, where k is the local wavenumber and L is the smallest scale among those that characterize the local properties of the medium, the wave envelope, and k itself. However, even for initially smooth fields, GO often predicts the appearance of caustics, where kL → 0. Examples of simple caustics include cutoffs (where k → 0) and focal points (where L → 0). The GO wavefield has spurious singularities at caustics, which is an issue for many applications of GO, such as calculating scattering cross sections [5,6] and modeling thermonuclear fusion experiments [7][8][9][10]. Hence, developing reduced models for handling caustics is an important practical problem.
In some cases, this problem is solved by locally reducing a given wave equation to a simpler one with a known solution, such as Airy's equation [11][12][13]. In other cases, cutoffs are modeled as discrete interfaces, such as in specular-reflection or perfect-conductor approximations [14][15][16]. However, such approaches assume that the spatial structure of the caustic is known a priori, which is not often the case. A more fundamental description of caustics was developed by Maslov [17] based on geometrical properties of GO solutions in the ray phase space. By occasionally rotating the phase space by π/2 using the Fourier transform in one or more spatial variables, one can remove a caustic and locally reinstate GO. Still, this approach is inconvenient for simulations because the rotation points have to be introduced ad hoc, requiring the simulations to be supervised. (Maslov's approach is discussed in more detail in Sec. II C.) Improvements have been made in the specific context of the Helmholtz equation [18], but not in the general case.
Here, we propose an alternate formulation of GO that can handle caustics without encountering these issues. We use the same general idea as in Maslov's method; however, instead of rotating the phase space by π/2 at select locations, we rotate the phase space continually along a GO ray. This is accomplished with the metaplectic transform (MT) [19,20], or more precisely, the near-identity metaplectic transform (NIMT) [21]. Importantly, we assume neither the caustic structure nor a specific wave equation. Using the Weyl symbol calculus [4], we formulate this procedure for any linear wave, including those governed by integro-differential equations. Hence, our approach can be useful for a wide variety of applications, such as in optics and in plasma physics.
This paper is organized as follows. In Sec. II, we review the basic equations of GO and introduce caustics. We also discuss Maslov's method in more detail. In Sec. III, our new approach, called metaplectic GO or MGO, is derived. Subsidiary results include an algorithm to explicitly construct the rotation matrices that our method employs and also the GO equations in an arbitrarily rotated phase space. In Sec. IV, we discuss two examples of one-dimensional (1-D) waves governed by Airy's equation and by Weber's equation. We study these systems analytically and show that MGO yields an accurate approximation at all locations, including near the cutoffs. In Sec. V, we summarize our main conclusions. Auxiliary calculations are presented in appendices.

A. The geometrical-optics equations
Consider an undriven linear scalar wave equation on an N -D configuration space with coordinates q, or q-space, which is assumed to be Euclidean 1 . In general, such an equation can be written as where ψ(q) is the wavefield on q, and D(q, q ′ ) is some dispersion kernel. Through use of Dirac δ-functions, Eq. (1) can describe local differential wave equations (but it can also describe integro-differential wave equations such as those common in plasma physics [4]). For example, choosing D(q, q ′ ) = ∇ ′2 δ(q ′ − q) + n 2 (q ′ )δ(q ′ − q) leads to the Helmholtz equation with a spatially varying index of refraction n(q). Here, ∇ and ∇ ′ are the gradient with respect to q and q ′ respectively. Let us introduce 2 a Hilbert space of state vectors |ψ such that ψ(q) is the projection of a given |ψ onto the eigenbasis {|q } of the coordinate operatorq. We adopt the usual normalization such that where1 is the identity operator. Then, where the symbol . = denotes definitions. We defineq through its action on the Hilbert space asq|q ′ = q ′ |q ′ , or equivalently, through its matrix elements The canonically conjugate momentum operatorp is similarly defined through its matrix elements as Let us further define the dispersion operatorD through its matrix elements q|D|q ′ = D(q, q ′ ). Then, Eq. (1) can be represented aŝ with Eq. (1) being simply the projection of Eq. (6) onto the coordinate eigenbasis. (Here, |0 is the null vector.) Note thatD is expressed as a function ofq andp. When the dispersion kernel D(q, q ′ ) describes a local differential wave equation, the construction ofD(q,p) is trivial. For example, the aforementioned Helmholtz equation haŝ D(q,p) = −p 2 + n 2 (q). However, constructingD(q,p) for integro-differential wave equations requires a pseudodifferential representation of D(q, q ′ ). Such a representation can be obtained using the Weyl symbol calculus, which we shall introduce momentarily. GO is the asymptotic model of Eq. (6) for the shortwavelength limit (Sec. I). In this limit, the wavefield can be partitioned into a rapidly varying phase and a slowly varying envelope. Following Ref. [22], we define the envelope state vector |φ via the unitary transformation where θ(q) is a hermitian operator representing the phase of ψ(q). Under this transformation, Eq. (6) becomes We shall now approximate the envelope dispersion operator in the GO limit. This is readily accomplished using the Weyl symbol calculus, which provides a mapping between functions and operators [24]. Let A(z) be a function on classical phase space with coordinates z . = (q, p) ⊺ . The Weyl transform W maps A(z) to an operatorÂ(ẑ) on the Hilbert space aŝ whereẑ . = (q,p) ⊺ and we have introduced the matrix with 0 N and I N being respectively the N -D null and identity matrices. (Here, ⊺ denotes the matrix transpose, which also denotes the scalar dot product for vectors, i.e., a ⊺ b ≡ a·b.) The function A(z) is called the Weyl symbol ofÂ(ẑ). The inverse Weyl transform, sometimes called the Wigner transform, maps operators back to phasespace functions as The Weyl symbol calculus is reviewed in Appendix A.
With the Weyl symbol calculus, approximating operators becomes as easy as approximating functions: one performs a Wigner transform to obtain the operator's Weyl symbol, approximates the symbol in the desired limit using, say, familiar Taylor expansions, then performs a Weyl transform to obtain the correspondingly approximated operator. As shown in Appendix B, applying this procedure to Eq. (8) ultimately yields where are interpreted respectively as the local dispersion relation, the local wavevector, and the local group velocity. Importantly, note that Projecting Eq. (12) onto q-space then yields the GO equations, where φ(q) . = q|φ . Specifically, Eq. (15a) represents a local dispersion relation, and Eq. (15b) represents an envelope equation. Note that for vector fields, additional polarization terms can emerge [22,25].

B. Ray trajectories and caustics
Equations (15) are commonly solved along the 1-D characteristic ray trajectories of Eq. (15a). These rays are described by Hamilton's equations where τ 1 is some variable used for ray parameterization.
Since the value of D(z) [Eq. (15a)] and ∇×p(q) [Eq. (14)] are conserved along each ray, the solutions to Eq. (16) are confined to an N -D surface in the 2N -D phase space called the dispersion manifold. A parametric representation of the dispersion manifold is readily obtained by integrating Eq. (16). Indeed, the formal solution to Eq. (16) can be written as z(τ ), where τ . = (τ 1 , τ ⊥ ) and τ ⊥ are coordinates on the (N − 1)-D sub-surface of initial conditions z 0 (τ ⊥ ) . = z(0, τ ⊥ ). Importantly, since Eq. (16) is a first-order autonomous system, rays can never cross in phase space; however, their projections onto q-space, q(τ ), have no such restriction. This can be problematic for the GO model [Eqs. (15)], which is constructed in q-space using q(τ ). To see why, let us integrate Eq. (15b) along a ray. The first term in Eq. (15b) is clearly the directional derivative of φ along the ray trajectory. Following Ref. [2], the second term is simplified upon noting that where tr denotes the matrix trace and v(τ ) by the chain rule, Jacobi's formula for the derivative of the matrix determinant implies that where J(τ ) . = det ∂ τ q(τ ) is the Jacobian determinant of the ray evolution in q-space, and log(x) is the natural logarithm. , det ∂τ q = 0 typically occurs at boundaries between regions where q(τ ) = q 0 has differing numbers of roots. This is illustrated above in 1-D, where q ′ (τ ) = 0 at q1, q2, and q3. An exception is when det ∂τ q = 0 corresponds to a degenerate saddlepoint (an inflection point in 1-D); however, these structures are not stable to small perturbations, so they are not often seen in physical systems.
To better understand where and why caustics occur, let us consider the extended 'ray parameter' space (q, τ ). In this space, the ray trajectories are represented by the graph τ = τ (q), which is obtained by a formal inversion of q(τ ). The condition J(τ ) = 0 has a geometric interpretation in this space: J(τ ) = 0 where the projection of τ (q) onto q-space becomes singular. Importantly, caustics do not occur every time rays cross in q-space, but rather, when the number of rays crossing in q-space changes abruptly. In this sense, caustics appear as topological boundaries (see Fig. 1).
Since τ are coordinates on the dispersion manifold, the same geometric interpretation of caustics must hold in phase space as well. Indeed, since it follows that Hence, J(τ ) = 0 where the dispersion manifold has a singular projection onto q-space as well. Formulating caustics as projection singularities in phase space is advantageous because it highlights the arbitrariness in the initial choice to project Eq. (6) onto q-space. More generally, a caustic occurs wherever the dispersion manifold has a singular projection onto the chosen projection plane. Consequently, caustics can be removed by rotating the projection plane.
C. Maslov's method for caustic removal A popular paradigm for performing such phase-space rotations is Maslov's method [17,26]. This method takes advantage of the fact that a well-behaved dispersion manifold cannot have a singular projection onto both q-space and p-space simultaneously. A caustic that appears in q-space is thus absent in p-space, and vice versa. By repeatedly switching between q-space and p-space as caustics are approached, one can construct a GO framework that does not produce singularities along a given ray.
To illustrate this method qualitatively, let us consider the dispersion manifold shown in Fig. 2. A q-space caustic occurs at q(τ 1 ), while p-space caustics occur at p(τ 2 ) and p(τ 3 ). Consider a wave initially located at q(τ 0 ). In region A, that is, for τ between τ 0 and some τ a to be specified momentarily, q-space can be used as the projection plane. Consequently, the wave envelope is evolved using Eq. (15b). Near the q-space caustic at q(τ 1 ), however, Eq. (15b) breaks down and cannot be used. Instead, the projection plane should be switched from q-space to p-space prior to encountering the caustic at τ 1 . The switching location, τ a , must be far enough from q-space and p-space caustics such that GO is accurate in both representations near τ a , but is otherwise arbitrary.
Next, the wavefield ψ(q) is transformed to its p-space representation Ψ(p) at τ a . This is achieved using the Fourier transform (FT) subsequently evaluated via the stationary phase approximation (SPA) [27]. Indeed, since the phase of the FT integrand is stationary where ∂ q θ(q) = p, which is satisfied along the dispersion manifold by definition. When τ a is chosen sufficiently far from both q-space caustics (such that φ is not singular) and p-space caustics (such that ∂ 2 q θ ≡ ∂ q p is nonzero), the SPA of Eq. (24) is manifold with coordinate τ that exhibits a q-space caustic at τ = τ1 and two p-space caustics at τ = τ2 and τ = τ3. Region A (τ ≤ τA) is far from the caustics, so both the q-space and p-space GO solutions are well-behaved. However, region B (τA < τ < τB) is close to the q-space caustic, so the q-space GO solution is singular while the p-space GO solution is well-behaved. Similarly, region C (τ ≥ τB) is close to the p-space caustics, so the q-space GO solution is well-behaved while the p-space GO solution is singular.
Thus, the SPA has the important role in Maslov's method of localizing the FT to become a pointwise mapping from Being absent from p-space caustics, Ψ(p) is evolved in region B from τ a to τ b using GO formulated in p-space. We shall derive the GO equations in various projection planes, including p-space, in the following section; for the moment, let us simply note that the p-space GO equations are not obtained by projecting Eq. (12) onto {|p }, but instead, more sophisticated machinery must be introduced. After propagating Ψ(p) through region B, the projection plane must be switched back to q-space to avoid the p-space caustic at p(τ 2 ). This is accomplished by using the inverse FT evaluated via SPA. Since there are no remaining q-space caustics, ψ(q) can be evolved using q-space GO for all subsequent τ > τ b .
Maslov's method has been very successful for the theoretical analysis of caustics. However, the lack of rigorous criteria for choosing when to switch between q-space and p-space is unsatisfying and can be awkward for practical calculations. For a code, this selection must be performed using an external module that supervises the envelope evaluation, detects when a caustic is becoming 'close' using some ad hoc cost function, then triggers a switch in representation. A framework that could proceed unsupervised would be more desirable. In the remainder of the paper, we shall present such a framework.

III. RESTORING GEOMETRICAL OPTICS USING PHASE-SPACE ROTATIONS
Maslov's method for caustic removal uses only a small subset of all the possible projective planes; indeed,

FIG. 3:
New framework for GO that is free from caustic singularities. This new approach consists of three steps: for a given point on the dispersion manifold, (i) the GO solution is calculated on the local tangent plane using metaplectic geometrical optics (MGO), (ii) the MGO solution is projected onto the tangent plane of the next point on the dispersion manifold using a near-identity metaplectic transform (NIMT) to initialize the next MGO calculation, (iii) the MGO solution is projected onto configuration space using a metaplectic transform (MT) subsequently evaluated using stationary phase or steepest-descent method. The process is repeated for all points on the dispersion manifold.
only q-space, p-space, and simple combinations such as (q x , p y ) are considered. This restriction ultimately results in an algorithm which must be supervised. However, allowing for a wider variety of possible projective planes can eliminate this need for supervision. This is the basis for our approach, which is outlined in Fig. 3.
In short, rather than sometimes switching between qspace and p-space as in Maslov's method, we propose to always switch between q-space and the local tangent plane of the dispersion manifold at a desired query point z(τ ). Each point on the dispersion manifold is thus treated equally, and there is no need to arbitrarily designate specific points as 'switching' points. Also, there will never be a caustic near z(τ ) by definition. For these reasons, our approach should be easy to implement in a code. Let us now develop this idea in more detail.

A. Phase-space rotations via metaplectic operators
Let us first discuss how to transform between q-space and the local tangent plane of the dispersion manifold. Generally speaking, for a given plane in phase space to be a valid choice for a GO projective plane, it must be related to q-space by a linear symplectic transformation. This is because linear symplectic transformations pre-serve the Poisson bracket, and consequently define an equivalency class on phase space [28].
A linear symplectic transformation is defined by a 2N × 2N matrix S which satisfies This matrix transforms the original phase space z into a new phase space Z − via or more explicitly, using Z − . = (Q, P) ⊺ and z . = (q, p) ⊺ , where A, B, C, and D are N × N block matrices that comprise S as follows: The inverse matrix, S −1 , has a similar block decomposition given as which is readily derived using Eqs. (26) and (28). Note that Eq. (27a) preserves the origin of phase space, i.e., z = 0 maps to Z − = 0. Shifting the origin does not affect the projective properties of a plane; hence, for the purposes of developing GO, we can identify all projective planes which differ only by a shift in the origin as equivalent. As a result, when we speak of the 'tangent plane' of the dispersion manifold at z(τ ), we really speak of the plane which is parallel to the tangent plane at z(τ ) and passes through the origin.
In the Hilbert space, linear symplectic transformations are performed using metaplectic operators. A metaplectic operatorM is a unitary operator which, via conjugation, transforms the operatorẑ toẐ − aŝ or in terms ofq andp, Configuration-space basis vectors are transformed as Correspondingly, wavefunctions of q-space are projected onto Q-space as [21] Ψ(Q) = dq Q|q q|ψ Here, σ . = ±1 is the overall sign of the MT, G 1 (q, Q) is the generating function [28] and our branch-cut convention restricts all complex phases onto the interval (−π, π]. Thus, σ determines the sign of √ det B. In writing Eq. (32), we have also assumed that B is invertible for simplicity. The generalization for non-invertible B is straightforward, and involves a δ-function kernel within the nullspace of B [19].
Equation (32) defines Ψ(Q) as the metaplectic transform (MT) of ψ(q). The inverse MT is given as with complex phases restricted to the interval [−π, π). Special cases of the MT include the FT and the fractional FT. When S is near-identity, the integral transform of Eq. (32) is well-approximated by a differential transform. For eikonal ψ(q), the near-identity MT (NIMT) which transforms ψ(q) to Ψ(Q) is given as [21] Ψ where we have dropped the term ∂ 2 qq φ(q) from F(q) because it is higher order in the GO parameter. Now, let us explicitly construct the symplectic matrix that maps q-space to the tangent plane of the dispersion manifold at some z(τ ). Recall that coordinates and coordinate axes transform oppositely (contravariantly versus covariantly). In other words, if the coordinates are transformed by S, then the coordinate axes are transformed by S −1 ; hence, we desire S −1 to map q-space to the local tangent space, rather than S. Let {Ť j (t)} and {Ň j (t)} be a symplectically dual set of N orthonormal tangent vectors and normal vectors to the dispersion manifold at τ = t, respectively. As can be readily verified, the matrix maps q-space to the local tangent space at τ = t. (The arrows emphasize that {Ť j (t)} and {Ň j (t)} form the columns of R t .) Since R t is unitary, we obtain The vectors {Ť j (t)} and {Ň j (t)} can be constructed using a 'symplectic Gram-Schmidt' method as follows. First, designatě Next, the complete set {Ť j (t)} is obtained aš where G s represents the 'modified Gram-Schmidt operator', which returns the orthogonalized version of its final argument with respect to the previous j − 1 arguments [29]. Finally, the complete set Let us confirm that as constructed, S t is indeed both orthogonal and symplectic. By Eq. (39), which immediately implies that Also, since the dispersion manifold is obtained by the gradient lift p = ∂ q θ(q), the dispersion manifold is a Lagrangian manifold [30]. A Lagrangian manifold has the property that any set of tangent vectors satisfy Consequently, Hence, S t is orthogonal per Eqs. (41a), (41b), and (41f), and is symplectic per Eqs. (41c), (41d), and (41e).

B. Geometrical optics in arbitrary projective plane
Let us now develop the GO equations on a projective plane that is obtained from q-space via some symplectic matrix S. We again consider the general wave equation given in Eq. (6), but, rather than introducing an eikonal ansatz on q-space as done in Eq. (7), we now assume the wavefield is eikonal on the desired projective plane. Hence, we perform the unitary transformation such that Eq. (6) becomes In principle, Eq. (43) is sufficient to develop GO; however, the simultaneous presence ofQ andẑ is inconvenient. To rectify this, let us introduce into Eq. (43) the metaplectic operator corresponding to S as where we have used the unitarity ofM . This allows us to rigorously transformẑ toẐ −. As shown in Appendix C, which also demonstrates the well-known 'symplectic covariance' property of the Weyl symbol [19,20]. In other words, the Weyl symbol of an operator at a given phasespace location does not depend on how this location is parameterized, as long as different parameterizations (here, z and Z) are connected via symplectic transformations.
Since symplectic transformations preserve the Poisson bracket, they also preserve the Moyal star product (Appendix A). Hence, the GO limit of Eq. (44) can be obtained using the procedure outlined in Appendix B, but replacing D(z) with D(S −1 Z) and z by Z. This yields where, given recent interest [22,31,32], we have included the quasioptical (QO) term which governs diffraction as We have also defined the following quantities: Dynamical equations that govern Φ(Q) . = Q|Φ are obtained by projecting Eq. (46) onto {|Q }. Neglecting diffraction, the GO equations on this projective plane are As before, Eq. (49a) is solved via ray tracing, while Eq. (49b) can be formally solved as However, let us make an observation regarding the rays generated by Eq. (49a). These rays satisfy Since S is constant, the chain rule yields Moreover, since S is symplectic, multiplication by S −1 from the left yields Finally, a comparison with Eq. (16) yields the relationship between the original and the transformed rays: Thus, one does not need to re-trace rays every time the projection plane is changed, but rather, one can simply perform the same symplectic transformation to the rays that one applies to the ambient phase space. The wavefield on the projective plane is constructed as summed over all branches of Θ(Q) if Θ(Q) is multivalued. The wavefield on q-space can then be obtained by taking the inverse MT of Ψ(Q) using Eq. (34).

C. Sequenced geometrical optics in a piecewise-linear tangent space
We are now equipped to discuss the new method of caustic removal outlined in Fig. 3. In this subsection we shall discuss steps (i) and (ii) of Fig. 3, that is, performing GO in the various tangent planes of the dispersion manifold and linking the obtained GO solutions using the NIMT, while in the following subsection we shall discuss step (iii). We shall assume that the dispersion manifold z(τ ) has already been obtained via Eq. (16). This is not a restrictive assumption, though, because the ray trajectories themselves are unaffected by caustics.
Let us consider some point q in configuration space and attempt to construct ψ(q). To do so, we shall map q to the dispersion manifold using the ray map τ (q), solve for Ψ(Q) in the optimal projection plane, that is, the tangent plane of the dispersion manifold at τ , and then map Ψ(Q) to ψ(q) using an inverse MT. In general, however, τ (q) will be multi-valued, and for the aforementioned scheme to work, the contributions to ψ(q) from each branch must be considered separately.
Let t ∈ τ (q) be a branch of τ (q), and let us construct the GO wavefield in the tangent plane at t. Equations (36) and (37) yield the S t that transforms q-space to the tangent plane. The rays are transformed into the new coordinates using Eq. (54) as Using the rays, P t (Q) can be constructed as where τ t (Q) is the function inverse of Q t (τ ). Then, the phase on the tangent plane is computed as while the envelope is computed using Eq. (49b) or its QO analogue. Importantly, when τ (q) is multi-valued, then P t (Q) should be restricted to the branch containing P t (t). This is because we are only interested in the contribution near t on the dispersion manifold. Let Φ t (Q) be the solution to Eq. (49b) (or its QO analogue) subject to the initial condition Φ t [Q t (t)] = 1. The wavefield on the tangent plane is where is a function of t which is required to make ψ(q) continuous. It can be found from where N S [ψ(q)] denotes the NIMT of ψ(q) with respect to S via Eq. (35). In other words, the arbitrary constant in Ψ t+h (Q) is obtained by projecting the neighboring Ψ t (Q) onto the tangent plane at τ = t + h.
In the continuous limit, Eq. (61) can also be written as a differential equation. Using Eqs. (35) and results presented in Appendix C of Ref. [21], one can show that where we have defined η t as The N × N matrices U, V, and W are obtained through the matrix decomposition which is possible because (∂ h S t ) S −1 t is a Hamiltonian matrix [21]. Consequently, U t and W t are both symmetric. We have also defined the directional derivative as Note that ∂ h is a total (directional) derivative in t, so it acts on both arguments of Q t (t), including the subscript. Then, Eq. (61) yields When α t is evolved along a ray, then ∂ h = ∂ t1 , and Eq. (66) is trivially integrated as The remaining constant function α (0,t ⊥ ) is determined by initial conditions.

D. Projecting tangent-space GO solution to configuration space
Equations (59) and (67), along with initial conditions, completely specify the wavefield in the tangent space of the dispersion manifold as Ψ(τ ) .
. This wavefield is free from caustics, and may be sufficient for certain applications. When ψ(q) is required, then Ψ t (Q) must be projected back onto q-space using Eq. (34) as 3 and the contribution to ψ t (q) from the field near t must be isolated by using the SPA. Indeed, we calculate The roots to Eq. (69) are Q = Q t [τ (q)]. Hence, restricting the integration domain near Q t (t) will isolate the contribution from t. To this end, let us define After performing a change in variables, Eq. (68) becomes where we have defined The overall sign σ t is constant unless det B t crosses the branch cut for the MT. Then, σ t changes sign to ensure that ψ t evolves smoothly in t. (See Sec. 2 of Ref. [21] for details or Sec. IV B for an example.) Invoking the SPA, the integration domain is restricted to the neighborhood of ǫ = 0, denoted δǫ. This yields No further approximations to Υ t can be made without additional knowledge of the caustic structure. We shall discuss some of these additional approximations in Sec. IV. Nevertheless, by properly choosing δǫ, it should be possible to numerically evaluate Eq. (73) to sufficiently high accuracy in the general case. Having isolated the field contribution from a single branch of the dispersion manifold, we sum over all such contributions to obtain E. Curvature-dependent adaptive discretization Since the accuracy of the NIMT depends on the rotation angle between neighboring tangent planes, the discretization of the dispersion manifold would ideally have smaller step sizes in regions with higher curvature. This can be achieved by replacing the ray Hamiltonian D(z) in Eq. (16) with a new function D(z) that generates the same ray trajectories in phase space but with a different 'time' parameterization. In this manner, adaptive integration schemes can be developed which are selfsupervising (do not need to be error-controlled in the traditional sense [33]), and amenable to symplectic methods of numerical integration [34,35].
Let us consider a modified ray Hamiltonian where f (z) is some smooth function. Let z D and z f denote the zero sets of D(z) and f (z) respectively, i.e., Clearly, the zero set of D(z) is where denotes the set union. For D(z) and D(z) to generate the same set of rays, f (z) and D(z) cannot be zero simultaneously. Hence, we require where denotes the set intersection and ∅ is the empty set. Equation (78) is most easily satisfied by requiring f (z) be sign-definite, say, positive (so z f = ∅).  [37]. For visualization purposes, the dispersion manifolds are displaced slightly from the origin, and the τ1 discretization is reflected about the p axis.
Let us now compute the rays generated by D(z). Analogous to Eq. (16), the new rays satisfy for some new parameterization τ 1 . From Eq. (79), rays initialized within z D will always remain in z D , at least in exact arithmetic. For such rays, the second term in Eq. (79) is identically zero, making Eq. (79) simply a reparameterization of Eq. (16) with For inexact arithmetic, however, the second term in Eq. (79) is not exactly zero, and therefore must be retained to preserve the Hamiltonian structure [35,36]. By Eq. (80), τ 1 is non-uniformly discretized when τ 1 is uniformly discretized; hence, a curvature-dependent adaptive discretization can be achieved by properly designing f (z). First, we restrict f (z) to only depend on the local curvature of the dispersion manifold, denoted K(z) 4 . Next, we impose that the uniform and the adaptive discretizations are equivalent when K(z) = 0. Hence, Finally, for the adaptive discretization to congregate in regions of high curvature, we require f (K) to be a strictly decreasing function of K, that is Any function which satisfies Eqs. (78), (81), and (82) will be a suitable choice, for example, with µ ≥ 0 a free parameter. Figure 4 shows an example f (K) which satisfies these three requirements, and shows the adaptive discretization generated by Eq. (83). As a final remark, when reparameterized rays are used to calculate Φ(Q), additional terms related to ∂ z f (z) will arise. These can be interpreted as the 'gravitational' forces associated with the 'time dilation' τ 1 → τ 1 [22,38].

IV. EXAMPLES
We now illustrate our methodology with a pair of examples, performed in 1-D for simplicity.

A. Airy's equation in one dimension
As a first example, let us consider a simple fold caustic in 1-D, which occurs when a wave encounters an isolated cutoff. For slowly varying media, this situation is often modeled with Airy's equation [11], Like with Eq. (6), we can also write Eq. (84) aŝ where we have first multiplied Eq. (84) by minus one for convenience. Using results presented in Appendix A, the Weyl symbol ofD(ẑ) is calculated to be In this case, the dispersion manifold D(z) = 0 is a parabola which opens along the negative q-axis. By Eq. (16), the rays are obtained by integrating The solution to Eq. (87) is where p 0 is a constant determined by initial conditions. From Eqs. (38) and (88), we compute the unit tangent vector to the dispersion manifold at some τ = t aš where we have defined Correspondingly, the normal vector to the dispersion manifold is calculated via Eq. (40) aš Using Eqs. (36) and (37), we therefore construct and identify the 1 × 1 block matrices Since B t never changes sign, then the sign of the inverse MT never changes either, and we can choose Using Eq. (54), the rays are transformed by S t as Equations (95) can be combined to obtain P t (Q) as Hence, P t (Q) is double-valued. As discussed following Eq. (58), we must restrict P t (Q) to the branch where P t [Q t (t)] = P t (t). This is accomplished by choosing the (+) sign in Eq. (96). The wavefield phase near Q t (t) is obtained by integrating Eq. (58); this ultimately yields where as a reminder, ǫ . = Q − Q t (t). It is easier analytically to obtain Φ t (Q) via Eq. (49b) rather than Eq. (50), since using Eq. (50) would require inverting Q t (τ ). (In numerical implementations, though, this may not be an issue). Therefore, we must compute Since  107) and (96), at t = −3, t = −0.5, t = 0, and t = 1. For each value of t, the projective plane Q is the tangent plane at Qt(t), denoted by the black dashed line. The caustic, denoted by the red dashed line, is located at Qc(t) where Pt(Q) has a vertical tangent line; we therefore see that under this rotation scheme, the neighborhood of Qt(t) is always free from caustics.

Re[¡t(Q)] Pt(Q)
we compute Depending on implementation, it may be more convenient to compute V t (Q) by using the chain rule to express where we have imposed that Φ t [Q t (t)] = 1 in accordance with the discussion preceding Eq. (59). As our final step in constructing Ψ t (Q), we must calculate α t . From Eq. (64), we compute We also compute Hence, using Eq. (63) we compute Performing a change in variables allows the integral in Eq. (67) to be evaluated as where we have chosen p(0) = 0, i.e., p 0 = 0. Thus, using Eqs. (59), (67), and (106), we obtain the wavefield on the tangent plane, where Θ t (Q) and Φ t (Q) are provided in Eqs. (97) and (102) respectively. For reference, Ψ t [Q t (t)] is plotted for select values of t in Fig. 5.
Next, we perform the inverse MT via Eq. (71) to obtain ψ t [q(t)]. From Eq. (72a), we compute Hence, Eq. (71) yields where Υ t is given in Eq. (73). We can now make further approximations to Υ t . First, since the integral of Υ t is restricted to δǫ, let us approximate Second, let us expand the exponential term about ǫ = 0. From Eq. (72c), we compute Therefore, and consequently, Since the coefficient of ǫ 2 vanishes when t = p 0 , while that of ǫ 3 remains nonzero for all t, Eq. (113) contains a Steepest descent paths for Airy's equation
pair of coalescing saddle points; isolating the contribution from the saddlepoint at ǫ = 0 analytically is therefore nontrivial. For now, let us suppose there is some 'saddlepoint filter' F which can perform this operation. [We shall soon show that F is related to choosing a steepest descent curve to evaluate Eq. (113).] Then, by definition of F , we can extend the integration bounds to infinity without incurring large errors. This yields By making the substitution ε . = iǫ/ϑ(t)−ip(t)ϑ 2 (t) and using Cauchy's integral theorem, the integral of Eq. (114) is placed into standard form where C 0 is a contour from re −iπ/3 to re iπ/3 , where r → ∞ (see Fig. 6). There are two saddlepoints in Eq. (115), located respectively at ε = ε ± . = ∓i|p(t)ϑ 2 (t)|. When Eq. (115) is evaluated using the method of steepest descent [39], the integration contour must be split into the two pieces denoted C + and C − , each of which encounters only one of the two saddle points (Fig. 6).
It is now clear how the 'saddlepoint filter' F can be designed: when F acts on Eq. (115), rather than evaluating the integral over both C + and C − , the integral is evaluated over either C + or C − . To determine which, note that when p > 0, ǫ = 0 corresponds to ε = −ipϑ 2 = −i|pϑ 2 | = ǫ + , while when p < 0, ǫ = 0 corresponds to ε = −ipϑ 2 = +i|pϑ 2 | = ǫ − . Hence, with As can readily be shown, where Ai(x) and Bi(x) are the Airy functions of the first and second kind, respectively [27]. Thus, using Eqs. (109), (114), and (116), we obtain We next sum over both branches of τ (q), per Eq. (74). Upon setting α 0 = − i/2π, we obtain where we have defined Figure 7 compares Eq. (120) with the exact solution, and with the standard GO approximation [27], As can be seen, the MGO solution is almost indistinguishable from the standard GO approximation far from the caustic at q = 0. However, in contrast to the standard GO solution, our solution remains finite for all q ≤ 0, like the exact solution of Eq. (84).

B. Weber's equation in one dimension
Next, let us consider a bounded wave in a 1-D harmonic potential, which exhibits two adjacent fold caustics. This situation is described mathematically by Weber's equation, which is also the Schrödinger equation for a quantum harmonic oscillator [3]. Equation (124) can be written aŝ where E . = ν + 1/2. As before, we have multiplied Eq. (124) by minus one for convenience. The Weyl symbol is readily computed to be In this case, the dispersion manifold D(z) = 0 is a circle of radius R . = √ 2E. From Eq. (16), the ray equations are Their solutions have the form where we have assumed the initial condition q(0) = R. The unit tangent and normal vectors at τ = t are calculated from Eqs. (38) and (40) aš Using Eqs. (36) and (37), we therefore construct and identify the 1 × 1 block matrices Unlike the previous example, B t can now change sign, and consequently, σ t will change sign as well. Let us choose to have B t cross the branch cut whenever B t changes from positive to negative. This is encapsulated by the phase convention where ⌊ ⌋ denotes the floor operation. Hence, choosing will ensure continuity in t across the branch cut. Using Eq. (54), the rays are transformed by S t as Equations (134) can be combined to obtain P t (Q) as Hence, P t (Q) is double-valued. As before, we restrict P t (Q) to the (+) branch. After integrating Eq. (58), we obtain the wavefield phase where ǫ . = Q−Q t (t) = Q, since per Eqs. (134), Q t (t) = 0. Next, since we compute Thus, using Eq. (49b) we obtain the wavefield envelope where as before, we set Φ t [Q t (t)] = Φ(0) = 1. We now calculate α t . From Eq. (64), we compute We also compute Hence, using Eq. (63) we compute

Qt(t) Qc(t) Qc(t)
FIG. 8: Same as Fig. 5 but for Weber's equation (124) for two different values of t and R. In all cases, Pt(Q) is a circle of radius R centered at the origin.
Thus, we evaluate The wavefield on the tangent plane is therefore Figure 8 shows Ψ t (Q) at select values of t and R. Next, we perform the inverse MT via Eq. (71) to obtain ψ t [q(t)]. From Eq. (72a), we compute Hence, Eq. (71) yields where Υ t is given in Eq. (73). Before proceeding, note that Υ t only involves functions which are either constant (Φ t , Θ t ) or π-periodic in time (S t , q). Thus, where we have used σ t+π = σ t exp (−iπ). For ψ t to be single-valued over the dispersion manifold, R 2 − 1 must be an even integer, which in turn requires ν to be an integer. Since E ≥ 0 is also needed for R to be real, the integer must be nonnegative. All together, this leads to the Bohr-Sommerfeld quantization of Weber's equation, more commonly known as [3] To evaluate Υ t , note that Eq. (72c) leads to We then approximate Υ t in the same manner as in the previous example; namely, we approximate and we approximate Consequently, Equation (152) is of the same basic form as Eq. (113). Therefore, we immediately conclude that where we have defined Hence, Upon setting α 0 = −i(iπ) −1/2 (2R) −5/6 , we sum over both branches of τ (q) to obtain the MGO solution where we have defined where D ν (x) is Whittaker's parabolic cylinder function [27], and with the standard GO approximation [40], ψ GO = 2 1/6 cos q 2 R 2 − q 2 − R 2 2 cos −1 q R + π 4 √ π R 1/3 (R 2 − q 2 ) 1/4 .
(159) Already for ν = 0, our MGO solution generally captures the exact solution behavior, albeit with some small error. Beginning from ν = 1, the agreement with the exact solution becomes even more remarkable, for either parity. Just like the previous example, our solution remains finite everywhere, whereas the standard GO solution becomes singular near the cutoffs at q = ±R.

V. CONCLUSION
In this work, we develop a reformulation of GO, called 'metaplectic GO', that is well-behaved near caustics and can be applied to any linear wave equation. MGO uses sequenced MTs to rotate the phase space continually along a ray such that caustics are never encountered. For each point on the dispersion manifold, (i) the phase space is rotated to align configuration space with the local tangent plane, (ii) GO is applied in the rotated phase space, (iii) the GO solution Ψ(Q) is linked to previous and subsequent GO calculations via NIMTs to ensure continuity, and (iv) Ψ(Q) is transformed back to the original phase space using an MT, summing over distinct branches of the dispersion manifold if applicable. This procedure should also be suitable for quasioptical modeling when generalized to non-Euclidean ray-based coordinates.
Our auxiliary results include: (i) an explicit construction of the rotation matrix for the tangent plane based on Gram-Schmidt orthogonalization with symplectic modifications; (ii) the GO equations for an arbitrarily rotated phase space (and more generally, a phase space obtained by an arbitrary linear symplectic transformation); and (iii) a simplified version of the MT to obtain ψ(q) from Ψ(Q) using the stationary phase approximation. This final result is restricted to det B = 0, which prohibits its use to analyze wave propagation in homogeneous media; however, a generalization that would allow for det B = 0 is straightforward. We also discuss two examples of 1-D linear wave problems and show analytically how MGO successfully approximates the exact wave dynamics.
where ⋆ denotes the Moyal star product. The Moyal star product is an associative non-commutative product rule for phase-space functions given explicitly as with J provided by Eq. (10). Importantly, Eqs. (A5) imply that the Weyl transform preserves both hermiticity and locality: the Weyl symbol of a hermitian operator is a function that is purely real, and the Weyl transforms of two functions which are 'close' to each other in the function space will also be 'close' in the operator space. Both of these properties make the Weyl symbol calculus an attractive means to approximate wave equations.
The following are some relevant Weyl transforms: where : denotes the double contraction.
Appendix B: Approximating the envelope equation in the GO limit via the Weyl symbol calculus Here we derive Eq. (12) from Eq. (8) using the Weyl symbol calculus. (See also Refs. [22,24].) To obtain the GO envelope operator, we shall (i) calculate the Weyl symbol of the envelope operator, (ii) approximate the Weyl symbol in the GO limit, and (iii) take the Weyl transform of the GO Weyl symbol.
Using Eqs. Here, κ is a dimensionless wavenumber that has been formally introduced to elucidate the GO ordering, where κ → ∞. We shall keep terms up to O(κ −2 ), rather than