Do you want to **read the rest** of this article?

# The backward phase flow and FBI-transform-based Eulerian Gaussian beams for the Schrödinger equation

**Article**

*in*Journal of Computational Physics 229(23):8888-8917 · November 2010

*with*93 Reads

**How we measure 'reads'**

A 'read' is counted each time someone views a publication summary (such as the title, abstract, and list of authors), clicks on a figure, or views or downloads the full-text. Learn more

DOI: 10.1016/j.jcp.2010.08.015 · Source: DBLP

Cite this publicationAbstract

We propose the backward phase flow method to implement the Fourier–Bros–Iagolnitzer (FBI)-transform-based Eulerian Gaussian beam method for solving the Schrödinger equation in the semi-classical regime. The idea of Eulerian Gaussian beams has been first proposed in [12]. In this paper we aim at two crucial computational issues of the Eulerian Gaussian beam method: how to carry out long-time beam propagation and how to compute beam ingredients rapidly in phase space. By virtue of the FBI transform, we address the first issue by introducing the reinitialization strategy into the Eulerian Gaussian beam framework. Essentially we reinitialize beam propagation by applying the FBI transform to wavefields at intermediate time steps when the beams become too wide. To address the second issue, inspired by the original phase flow method, we propose the backward phase flow method which allows us to compute beam ingredients rapidly. Numerical examples demonstrate the efficiency and accuracy of the proposed algorithms.

- ... In this paper, motivated by the work in [18] which designed fast multiscale Gaussian wavepacket transforms and multiscale Gaussian beam methods for pure initial value problems of wave equations, we propose to develop fast multiscale Gaussian beam methods for wave equations in bounded convex domains. Although Gaussian-beam based numerical methods are well developed for pure initial value problems of Schrödinger equation and wave equations [11,22,9,21,10,17,18,14,23,20,12], it seems that no efficient numerical Gaussian beam method has been developed for wave equations and Schrödinger equations in bounded domains. In the case of wave equations in bounded domains, some essential difficulties arise in developing numerical Gaussian beams. ...... An integral superposition of these solutions can be used to define a more general solution that is not necessarily concentrated on a single curve. Gaussian beams can be used to treat pseudo-differential equations in a natural way, including Helmholtz and Schrödinger equations [11,22,9,21,10,17,18, 14,20,12]. Gaussian beam superpositions have been used in geophysical applications for seismic wave modeling [4] and migra- tion [6] . ...... The numerical implementations in these areas are based on ray-centered coordinates which prove to be computationally inefficient [4,6]. More recently, based on [19,22] the first Eulerian Gaussian beam method was proposed in [11] which overcomes some of these difficulties; it can be easily applied to both high frequency waves and semi-classical quantum mechanics [9,10]. In [22] Lagrangian Gaussian beams are successfully constructed to simulate mountain waves, a kind of stationary gravity wave forming over mountain peaks and interfering with aviation. ...Article
- Mar 2014
- J COMPUT PHYS

Motivated by fast multiscale Gaussian wavepacket transforms and multiscale Gaussian beam methods which were originally designed for pure initial-value problems of wave equations, we develop fast multiscale Gaussian beam methods for initial boundary value problems of wave equations in bounded convex domains in the high frequency regime. To compute the wave propagation in bounded convex domains, we have to take into account reflecting multiscale Gaussian beams, which are accomplished by enforcing reflecting boundary conditions during beam propagation and carrying out suitable reflecting beam summation. To propagate multiscale beams efficiently, we prove that the ratio of the squared magnitude of beam amplitude and the beam width is roughly conserved, and accordingly we propose an effective indicator to identify significant beams. We also prove that the resulting multiscale Gaussian beam methods converge asymptotically. Numerical examples demonstrate the accuracy and efficiency of the method. - ... Sparse approximations. Motivated by recent works in [14, 15, 24, 25, 27], we develop sparse evolution approaches to further speed up the above algorithms. To develop a sparse evolution approach, one needs to view the solution in the transformed domain, where the original solution has a sparse representation in the sense that only a few significant coefficients suffice to capture the overall behavior of the solution. ...... To develop a sparse evolution approach, one needs to view the solution in the transformed domain, where the original solution has a sparse representation in the sense that only a few significant coefficients suffice to capture the overall behavior of the solution. In [14, 15, 24, 25] the high frequency wave solutions are first represented in phase space by using a phase space transform, such as the wavepacket transform, then a hard-thresholding process is applied to the phase space representation so that a few modes in the phase space suffice to capture the solution; this strategy results in efficient high-frequency wave propagation. In [27], the authors have proposed a new framework to efficiently approximate solutions to PDEs by sparse dynamics. ...... Because the Schrödinger equation generates new scales during the evolution process , the methodology of [27] cannot be applied here. On the other hand, we may still replace the hard-thresholding process as used in [14, 15, 24, 25] with the softthresholding in the evolution process; consequently, although we are using the softthresholding operator, our methodology is different from the one in [27]. In the following, we adopt the soft-thresholding process to both the Taylorexpansion based spectral methods and the Taylor-expansion based convolution methods . ...ArticleFull-text available
- Jan 2014

We propose fast Huygens sweeping methods for Schrödinger equations in the semi-classical regime by incorporating short-time Wentzel-Kramers-Brillouin-Jeffreys (WKBJ) propaga-tors into Huygens' principle. Even though the WKBJ solution is valid only for a short time period due to the occurrence of caustics, Huygens' principle allows us to construct the global-in-time semi-classical solution. To improve the computational efficiency, we develop analytic approximation formulas for the short-time WKBJ propagator by using the Taylor expansion in time. These analytic formulas allow us to develop two classes of fast Huygens sweeping methods, among which one is posed in the momentum space, and the other is posed in the position space, and both of these methods are of computational complexity O(N log N) for each time step, where N is the total number of sampling points in the d-dimensional position space. To further speed up these methods, we also incorporate the soft-thresholding sparsification strategy into our new algorithms so that the computational cost can be further reduced. The methodology can also be extended to nonlinear Schrödinger equations. One, two, and three dimensional examples demonstrate the performance of the new algorithms. - ... The sum can be discrete [38] [47] [2] or continuous [30] [31], the selection of the beams to be superposed can be done according to several criteria. In [4], a weighted integral of Gaussian beams was designed to build an approximate solution of the IBVP (1.1) under an additional assumption (H5) on the initial data (see p.9). See also [27] [28] [43] for recent numerical implementations related to this method. Gaussian beams seem to be very well suited for the study of Wigner measures. ...... The sum can be discrete [38, 47, 2] or continuous [30, 31], the selection of the beams to be superposed can be done according to several criteria. In [4], a weighted integral of Gaussian beams was designed to build an approximate solution of the IBVP (1.1) under an additional assumption (H5) on the initial data (see p.9). See also [27, 28, 43] for recent numerical implementations related to this method. Gaussian beams seem to be very well suited for the study of Wigner measures. ...A Gaussian beam method is presented for the analysis of the energy of the high frequency solution to the mixed problem of the scalar wave equation in an open and convex subset, with initial conditions compactly supported in this set, and Dirichlet or Neumann type boundary condition. The transport of the microlocal energy density along the broken bicharacteristic flow at the high frequency limit is proved through the use of Wigner measures. Our approach consists first in computing explicitly the Wigner measures under an additional control of the initial data allowing to approach the solution by a superposition of first order Gaussian beams. The results are then generalized to standard initial conditions.
- ... Consequently, Operto et al. (2000) extend the asymptotic ray- Born migration/inversion originally designed to process one single arrival to the case of multiple arrivals using ray theory-based dynamically correct multiarrival Green's functions, so as to account for the cross contributions of all the source-receiver raypaths; naturally , it is nontrivial to construct such ray-theory Green's functions because they have to explicitly keep track of caustics, ray branches, and the Keller-Maslov index, so that multivalued phases and amplitudes are correctly taken into consideration (Červený et al., 1977; Chapman, 1985). As an alternative, Gaussian beams are able to take care of multiple arrivals and caustics automatically in ray theory asymptotics for Green's functions at the cost of expensive beam summations (Červený et al., 1982; Popov, 1982; Ralston, 1983; Hill, 1990; Albertin et al., 2004; Gray, 2005; Leung et al., 2007; Tanushev et al., 2007; Hu and Stoffa, 2009; Leung and Qian, 2010; Qian and Ying, 2010a, 2010b; Bao et al., 2014). On the other hand, Bevc (1997) proposes a semirecursive Kirchhoff depth migration method to image complex structures by combining wave equation datuming (Berryhill, 1979) and first arrival-based Kirchhoff migration (Geoltrain and Brac, 1993). ...Multiarrival Green's functions are essential in seismic modeling, migration, and inversion. Huygens-Kirchhoff (HK) integrals provide a bridge to integrate locally valid first-arrival Green's functions into a globally valid multiarrival Green's function. We have designed robust and accurate finite-difference methods to compute first-arrival traveltimes and amplitudes, so that first-arrival Green's functions can be constructed rapidly. We adapted a fast butterfly algorithm to evaluate discretized HK integrals. The resulting fast Huygens' sweeping method has the following unique features: (1) it precomputes a set of local traveltime and amplitude tables, (2) it automatically takes care of caustics, (3) it constructs Green's functions of the Helmholtz equation for arbitrary frequencies and for many point sources, and (4) for a fixed number of points per wavelength, it constructs each Green's function in nearly optimal complexity O{eth}N log N{Thorn} in terms of the total number of mesh points N, where the prefactor of the complexity only depends on the specified accuracy, and is independent of the frequency. The 2D and 3D examples revealed the performance of the method.
- ... The third approach is based on Gaussian beam methods [12,40,44,52,25, 48,47]. Although Gaussian beam methods can treat caustics automatically along a central ray, the method itself suffers from expensive beam summation and exponential growth of beam width as analyzed and illustrated in [42,43,23,24,36,47], and such shortcomings sometimes have hindered applications of Gaussian beam methods to complicated inhomogeneous media. Our proposed new method is different from the above three approaches. ...
- ... In order to apply the Gaussian beam methods efficiently and accurately, one has to generate a beam decomposition for a general initial condition. There have been some recent advances in beam decomposition for the wave equation [19] [26] [27] and the Schrödinger equation [17] [18] [22]. In particular, [22] introduced single-scale Gaussian wavepacket transforms and developed on top of them a highly efficient initialization algorithm for the Schrödinger equation. ...The Gaussian beam method is an asymptotic method for wave equations with highly oscillatory data. In a recent published paper by two of the authors, a multiscale Gaussian beam method was first proposed for wave equations by utilizing the parabolic scaling principle and multiscale Gaussian wavepacket transforms, and numerical examples there demonstrated excellent performance of the multiscale Gaussian beam method. This article is concerned with the important convergence properties of this multiscale method. Specifically, the following results are established. If the Cauchy data are in the form of non-truncated multiscale Gaussian wavepackets, the multiscale Gaussian beam method provides a convergent parametrix for the wave equation with highly oscillatory data, and the convergence rate is , where λ is the smallest frequency contained in the highly oscillatory data. If the highly oscillatory Cauchy data are in the form of truncated multiscale Gaussian wavepackets, the multiscale Gaussian beam method converges with a rate controlled by , where ε is the error from initializing the Gaussian beam method by multiscale Gaussian wavepacket transforms. To prove these convergence results, it is essential to characterize multiscale properties of wavepacket interaction and beam decaying by carrying out some highly-oscillatory integrals of Fourier-integral-operator type, so that those multiscale properties lead to precise convergence orders for the multiscale Gaussian beam method.
- Article
- Jan 2014
- SIAM J NUMER ANAL

Gaussian beams are generally local asymptotic solutions to the linear wave equations in the high-frequency regime. Each Gaussian beam is concentrated around a specific ray path determined by the underlying Hamiltonian system. Expressed as some superposition of Gaussian beams, Gaussian beam approximation is expected to be a high-frequency asymptotic solution which remains globally valid even around caustics. We derive optimal first-order error estimates for first-order Gaussian beam approximations, in both continuous and discrete levels, to the Schrodinger equation equipped with a WKB initial data. Our error estimates are valid for any spatial dimension and unaffected by the presence of caustics. Some numerical tests are presented to validate the theoretical results. - Article
- Jul 2015
- COMMUN COMPUT PHYS

We propose a new semi-implicit level set approach to a class of curvature dependent flows. The method generalizes a recent algorithm proposed for the motion by mean curvature where the interface is updated by solving the Rudin-Osher-Fatemi (ROF) model for image regularization. Our proposal is general enough so that one can easily extend and apply the method to other curvature dependent motions. Since the derivation is based on a semi-implicit time discretization, this suggests that the numerical scheme is stable even using a time-step significantly larger than that of the corresponding explicit method. As an interesting application of the numerical approach, we propose a new variational approach for extracting limit cycles in dynamical systems. The resulting algorithm can automatically detect multiple limit cycles staying inside the initial guess with no condition imposed on the number nor the location of the limit cycles. Further, we also propose in this work an Eulerian approach based on the level set method to test if the limit cycles are stable or unstable. - Article
- Jul 2017
- J COMPUT PHYS

We propose a numerically efficient algorithm for simulating the multi-color optical self-focusing phenomena in nematic liquid crystals. The propagation of the nematicon is modeled by a parabolic wave equation coupled with a nonlinear elliptic partial differential equation governing the angle between the crystal and the direction of propagation. Numerically, the paraxial parabolic wave equation is solved by a fast Huygens sweeping method, while the nonlinear elliptic PDE is handled by the alternating direction explicit (ADE) method. The overall algorithm is shown to be numerically efficient for computing high frequency beam propagations. - Article
- May 2016
- J COMPUT PHYS

In some applications, it is reasonable to assume that geodesics (rays) have a consistent orientation so that the Helmholtz equation can be viewed as an evolution equation in one of the spatial directions. With such applications in mind, starting from Babich's expansion, we develop a new high-order asymptotic method, which we dub the fast Huygens sweeping method, for solving point-source Helmholtz equations in inhomogeneous media in the high-frequency regime and in the presence of caustics. The first novelty of this method is that we develop a new Eulerian approach to compute the asymptotics, i.e. the traveltime function and amplitude coefficients that arise in Babich's expansion, yielding a locally valid solution, which is accurate close enough to the source. The second novelty is that we utilize the Huygens–Kirchhoff integral to integrate many locally valid wavefields to construct globally valid wavefields. This automatically treats caustics and yields uniformly accurate solutions both near the source and remote from it. The third novelty is that the butterfly algorithm is adapted to accelerate the Huygens–Kirchhoff summation, achieving nearly optimal complexity , where N is the number of mesh points; the complexity prefactor depends on the desired accuracy and is independent of the frequency. To reduce the storage of the resulting tables of asymptotics in Babich's expansion, we use the multivariable Chebyshev series expansion to compress each table by encoding the information into a small number of coefficients. - Article
- Jun 2014
- J COMPUT PHYS

We propose a new Eulerian tool to study complicated dynamical systems based on the average growth in the surface area of a family of level surfaces represented implicitly by a level set function. Since this proposed quantity determines the temporal variation of the averaged surface area of all level surfaces, we name the quantity the Variation of the Integral over Area of Level Surfaces (VIALS). Numerically, all these infinitely many level surfaces are advected according to the given dynamics by solving one single linear advection equation. To develop a computationally efficient approach, we apply the coarea formula and rewrite the surface area integral as a simple integral relating the total variation (TV) of the level set function. The proposed method can be easily incorporated with a recent Eulerian algorithm for efficient computation of flow maps to speed up our approach. We will also prove that the proposed VIALS is closely related to the computation of the finite time Lyapunov exponent (FTLE) in the Lagrangian coherent structure (LCS) extraction. This connects our proposed Eulerian approach to widely used Lagrangian techniques for understanding complicated dynamical systems. - Article
- May 2014
- J COMPUT PHYS

We develop an efficient Eulerian numerical approach to extract invariant sets in a continuous dynamical system in the extended phase space (the x–t space). We extend the idea of ergodic partition and propose a concept called coherent ergodic partition for visualizing ergodic components in a continuous flow. Numerically, we first apply the level set method [33] and extend the backward phase flow method [25] to determine the long time flow map. To compute all required long time averages of observables along particle trajectories, we propose an Eulerian approach by simply incorporating flow maps to iteratively interpolate those short time averages. Numerical experiments will demonstrate the effectiveness of the approach. As an application of the method, we apply the approach to the field of geometrical optics for high frequency wave propagation and propose to use the result from the coherent ergodic partition as a criteria for adaptivity in typical Lagrangian ray tracing methods. - Article
- Mar 2017
- J SCI COMPUT

We propose and analyze a new class of Eulerian methods for constructing both the forward and the backward flow maps of sufficiently smooth dynamical systems. These methods improve previous Eulerian approaches so that the computations of the forward flow map can be done on the fly as one imports or measures the velocity field forward in time. Similar to typical Lagrangian or semi-Lagrangian methods, the proposed methods require an interpolation at each step. Having said that, the Eulerian method interpolates d components of the flow maps in the d dimensional space but does not require any \((d+1)\)-dimensional spatial-temporal interpolation as in the Lagrangian approaches. We will also extend these Eulerian methods to compute line integrals along any Lagrangian particle. The paper gives a computational complexity analysis and an error estimate of these Eulerian methods. The method can be applied to a wide range of applications for flow map constructions including the finite time Lyapunov exponent computations, the coherent ergodic partition, and high frequency wave propagations using geometric optic. - Article
- Dec 2013
- Chaos

We propose a simple Eulerian approach to compute the moderate to long time flow map for approximating the Lyapunov exponent of a (periodic or aperiodic) dynamical system. The idea is to generalize a recently proposed backward phase flow method which is specially designed for long time level set propagation. Unlike the original phase flow method or the backward phase flow method, which is applicable only to autonomous systems, the current approach can also be applied to any time-dependent (periodic or aperiodic) flow. We will discuss the stability of the proposed method. Numerical examples will be given to demonstrate the effectiveness of the algorithm. - Article
- Sep 2013
- WAVE MOTION

The Eulerian Gaussian beam method is an efficient way to compute high frequency wave propagation, which was originally studied in Leung et al. (2007) [17]. Later Jin, Wu and Yang developed a new way of computing the Hessian functions from the derivatives of the level set functions in Jin et al. (2008) [19], which greatly reduced the number of equations in the Eulerian Gaussian beam method.In this paper, we generalize this new method (JWY-EGBM) to compute high frequency wave propagation in the reduced momentum space. The difficulty lies in that, the dimensionality of momentum space is one less than that of configuration space, while JWY-EGBM requires configuration and momentum spaces be equally dimensional. We present two numerical examples to show the performance of the proposed method. - Article
- Jan 2013
- J COMPUT PHYS

In this paper an efficient and accurate method for simulating the propagation of a localized solution of the Schrödinger equation with a smooth potential near the semiclassical limit is presented. We are interested computing arbitrarily accurate solutions when the non dimensional Planck’s constant, εε, is small, but not negligible. The method is based on a time dependent transformation of the independent variables, closely related to Gaussian wave packets. A rescaled wave function, w, satisfies a new Schrödinger equation with a time dependent potential which is a perturbation of the harmonic oscillator, the perturbation being O(ε), so that all stiffness (in space and time) is greatly reduced. In fact, for integration in a fixed time interval, the number of modes required to fully resolve the problem decreases when εε is decreased. The original wave function may be reconstructed by Fourier interpolation, although expectation values of the observables can be computed directly from the function w itself. If the initial condition is a Gaussian wave packet, very few modes are necessary to fully resolve the w variable, so for short time very accurate solutions can be obtained at low computational cost. Initial conditions other than Gaussians wave packets can also be used. In this paper, the Gaussian wave packet transform is carefully outlined and applied to the Schrödinger equation in one dimension. Detailed numerical tests show the efficiency and accuracy of the approach. In the sequel of this paper, this approach has been extended to the multidimensional case. - Article
- Nov 2017
- J SCI COMPUT

We propose a new class of Lagrangian approaches for constructing the flow maps of given dynamical systems. In the case when only discrete velocity data at mesh points is available, an interpolation step will be required. However, in our proposed approaches all particle trajectories share a common global interpolation at each time step and therefore interpolation operations will not increase the overall computational complexity. The old Lagrangian approaches propose to solve the corresponding ordinary differential equations (ODEs) backward in time to obtain the backward flow map. It is inconvenient and not natural, especially when incorporated with certain computational fluid dynamic solvers, because the velocity field needs to be loaded from the terminal time backward to the initial time. In contrast, our proposed approaches for computing the backward flow map propose to solve the corresponding ODEs forward in time which is more practical. We will also extend the proposed approaches to compute line integrals along any particle trajectory. This paper gives a detailed analysis on the computational complexity and error estimate of the proposed Lagrangian approach. Finally, a wide range of applications of our approaches will be given, including the so-called coherent ergodic partition and the high frequency wave propagations based on geometric optic.

- We develop an efficient numerical method to compute single slit or dou-ble slit diffraction patterns from high frequency wave in inhomogeneous media. We approximate the high frequency asymptotic solution to the Helmholtz equation using the Eulerian Gaussian beam summation proposed in [20, 21]. The emitted rays from a slit are embedded in the phase space using an open segment. The evolution of this open curve is accurately computed using the recently developed Grid Based Particle Method [24] which results in a very efficient computational algorithm. Following the grid based particle method we proposed in [23, 24], we represent the open curve or the open surface by meshless Lagrangian particles sampled according to an under-lying fixed Eulerian mesh. The end-points of the open curve are tracked explicitly and consistently with interior particles. To construct the overall wavefield, each of these sampling particles also carry necessary quantities that are obtained by solving advection-reaction equations. Numerical experiments show that the resulting method can model diffraction patterns in inhomogeneous media accurately, even in the occur-rence of caustics.
- Book
- Jan 2002

Introduction * Semiclassical Pseudodifferential Calculus * Microlocalization * Applications to the Solutions of Analytic Linear PDEs * Complements: Symplectic Aspects * Appendix: List of Formulae * Bibliography * Index * List of Notations - We consider the scalar wave equation in a bounded convex domain of ℝ n . The boundary condition is of Dirichlet or Neumann type and the initial conditions have a compact support in the considered doamin. We construct a family of approximate high frequency solutions by a Gaussian beams summation. We give a rigorous justification of the asymptotics in the sense of an energy estimate and show that the error can be reduced to any arbitrary power of ε, which is the high frequency parameter.
- We design an Eulerian Gaussian beam summation method for solving Helmholtz equations in the high-frequency regime. The traditional Gaussian beam summation method is based on Lagrangian ray tracing and local ray-centered coordinates. We propose a new Eulerian formulation of Gaussian beam theory which adopts global Cartesian coordinates, level sets, and Lionville equations, yielding uniformly distributed Eulerian traveltimes and amplitudes in phase space simultaneously for multiple sources. The time harmonic wavefield can be constructed by summing up Gaussian beams with ingredients provided by the new Eulerian formulation. The conventional Gaussian beam summation method can be derived from the proposed method. There are three advantages of the new method: (1) We have uniform resolution of ray distribution. (2) We can obtain wavefields excited at different sources by varying only source locations in the summation formula. (3) We can obtain wavefields excited at different frequencies by varying only frequencies in the summation formula. Numerical experiments indicate that the Gaussian beam summation method yields accurate asymptotic wavefields even at caustics. The new method may be used for seismic modeling and migration.
- Gaussian beams are approximate solutions to hyperbolic partial differential equations that are concentrated on a curve in space-time. In this paper, we present a method for computing the stationary in time wave field that results from steady air flow over topography as a superposition of Gaussian beams. We derive the system of equations that governs these mountain waves as a linearization of the basic equations of fluid dynamics and show that this system is well-posed. Furthermore, we show that the approximate Gaussian beam stationary solution is close to a true time-dependent solution of the linearized system.
- A new semiclassical approach that constructs the full semiclassical Green’s function propagation of any initial wave function directly from an ensemble of real trajectories, without root searching, is presented. Each trajectory controls a cell of initial conditions in phase space, but the cell area is not constrained by Planck’s constant. The method is shown to be accurate for rather long times in anharmonic oscillators, indicating the semiclassical time‐dependent Green’s function is clearly worthy of more study. The evolution of wave functions in anharmonic potentials is examined and a spectrum from the semiclassical correlation function is calculated, comparing with exact fast Fourier transform results.
- Article
- Dec 2006
- J COMPUT PHYS

This paper introduces the phase flow method, a novel, accurate and fast approach for constructing phase maps for nonlinear autonomous ordinary differential equations. The method operates by initially constructing the phase map for small times using a standard ODE integration rule and builds up the phase map for larger times with the help of a local interpolation scheme together with the group property of the phase flow. The computational complexity of building up the complete phase map is usually that of tracing a few rays. In addition, the phase flow method is provably and empirically very accurate. Once the phase map is available, integrating the ODE for initial conditions on the invariant manifold only makes use of local interpolation, thus having constant complexity.The paper develops applications in the field of high frequency wave propagation, and shows how to use the phase flow method to (1) rapidly propagate wave fronts, (2) rapidly calculate wave amplitudes along these wave fronts, and (3) rapidly evaluate multiple wave arrival times at arbitrary locations. - Article
- Jan 1986
- J CHEM PHYS

The propagation of an initially highly excited localized wave packet in an anharmonic oscillator potential is studied within the frozen Gaussian approximation. Comparison is made to quantum mechanical basis set calculations. The frozen Gaussian approximation involves the expansion of the initial wave function in terms of an overcomplete Gaussian basis set. The wave function evolution is evaluated by allowing each Gaussian to travel along a classical trajectory with its shape held rigid. A Monte Carlo algorithm is employed in the selection of the initial Gaussian basis functions. The frozen Gaussian results are very good for times on the order of a few vibrational periods of the oscillator and remain qualitatively correct for the entire length of the calculations which is 12 vibrational periods. The dependence of the calculations on the width of the Gaussian basis functions is investigated and the effect of a simplifying approximation for the prefactor of the Gaussians is tested. - Article
- Jan 1985
- J CHEM PHYS

The propagation of an initially highly excited localized wave packet in an anharmonic oscillator potential is studied within the frozen Gaussian approximation. Comparison is made to quantum mechanical basis set calculations. The frozen Gaussian approximation involves the expansion of the initial wave function in terms of an overcomplete Gaussian basis set. The wave function evolution is evaluated by allowing each Gaussian to travel along a classical trajectory with its shape held rigid. A Monte Carlo algorithm is employed in the selection of the initial Gaussian basis functions. The frozen Gaussian results are very good for times on the order of a few vibrational periods of the oscillator and remain qualitatively correct for the entire length of the calculations which is 12 vibrational periods. The dependence of the calculations on the width of the Gaussian basis functions is investigated and the effect of a simplifying approximation for the prefactor of the Gaussians is tested. - Article
- Nov 1990
- GEOPHYSICS

This paper describes a zero-offset depth migration method that employs Gaussian beam downward continuation of the recorded wave field. The Gaussian-beam migration method has advantages for imaging complex structures. Like finite-difference migration, it is especially compatible with lateral variations in velocity, but Gaussian beam migration can image steeply dipping reflectors and will not produce unwanted reflections from structure in the velocity model. Unlike other raypath methods. Gaussian beam migration has guaranteed regular behavior at caustics and shadows. In addition, the method determines the beam spacing that ensures efficient, accurate calculations. The images produced by Gaussian beam migration are usually stable with respect to changes in beam parameters. -from Author - Article
- Jul 1982
- Geophys J Int

An asymptotic procedure for the computation of wave fields in two-dimensional laterally inhomogeneous media is proposed. It is based on the simulation of the wave field by a system of Gaussian beams. Each beam is continued independently through an arbitrary inhomogeneous structure. The complete wave field at a receiver is then obtained as an integral superposition of all Gaussian beams arriving in some neighbourhood of the receiver. The corresponding integral formula is valid even in various singular regions where the ray method fails (the vicinity of caustic, critical point, etc.). Numerical examples are given. - We review some of the issues facing semiclassical methods in classically chaotic systems, then demonstrate the long-time accuracy of semiclassical propagation of a nonstationary wave packet using the quantum baker's map of Balazs and Voros. We show why some of the standard arguments against the efficacy of semiclassical dynamics for long-time chaotic motion are incorrect.
- Article
- Oct 2010
- J COMPUT PHYS

This paper introduces a wavepacket-transform-based Gaussian beam method for solving the Schrödinger equation. We focus on addressing two computational issues of the Gaussian beam method: how to generate a Gaussian beam representation for general initial conditions and how to perform long time propagation for any finite period of time. To address the first question, we introduce fast Gaussian wavepacket transforms and develop on top of them an efficient initialization algorithm for general initial conditions. Based on this new initialization algorithm, we address the second question by reinitializing the beam representation when the beams become too wide. Numerical examples in one, two, and three dimensions demonstrate the efficiency and accuracy of the proposed algorithms. The methodology can be readily generalized to deal with other semi-classical quantum mechanical problems. - Article
- Mar 1990
- J COMPUT PHYS

We compare four pseudo-spectral split-step methods for solving a class of nonlinear Schrödinger (NLS) equations. The importance of observing the L2 invariance of the continuous problem is demonstrated through numerical experiments. The best performance is obtained by transforming the given equation to an NLS equation where two of the coefficients satisfy a simple algebraic relationship. The problem can be solved efficiently in terms of the new variables, and the l2 norm of the computed solution is time-invariant. - In this paper we study time-splitting spectral approximations for the linear Schrödinger equation in the semiclassical regime, where the Planck constant ϵ is small. In this regime, the equation propagates oscillations with a wavelength of O(ϵ), and finite difference approximations require the spatial mesh size h=o(ϵ) and the time step k=o(ϵ) in order to obtain physically correct observables. Much sharper mesh-size constraints are necessary for a uniform L2-approximation of the wave function. The spectral time-splitting approximation under study will be proved to be unconditionally stable, time reversible, and gauge invariant. It conserves the position density and gives uniform L2-approximation of the wave function for k=o(ϵ) and h=O(ϵ). Extensive numerical examples in both one and two space dimensions and analytical considerations based on the Wigner transform e ven show that weaker constraints (e.g., k independent of ϵ, and h=O(ϵ)) are admissible for obtaining “correct” observables. Finally, we address the application to nonlinear Schrödinger equations and conduct some numerical experiments to predict the corresponding admissible meshing strategies.
- We introduce a new multiscale Gaussian beam method for the numerical solution of the wave equation with smooth variable coefficients. The first computational question addressed in this paper is how to generate a Gaussian beam representation from general initial conditions for the wave equation. We propose fast multiscale Gaussian wavepacket transforms and introduce a highly efficient algorithm for generating the multiscale beam representation for a general initial condition. Starting from this multiscale decomposition of initial data, we propose the multiscale Gaussian beam method for solving the wave equation. The second question is how to perform long time propagation. Based on this new initialization algorithm, we utilize a simple reinitialization procedure that regenerates the beam representation when the beams become too wide. Numerical results in one, two, and three dimensions illustrate the properties of the proposed algorithm. The methodology can be readily generalized to treat other wave propagation problems.
- Article
- May 2009
- J COMPUT PHYS

We propose Gaussian-beam based Eulerian methods to compute semi-classical solutions of the Schrödinger equation. Traditional Gaussian beam type methods for the Schrödinger equation are based on the Lagrangian ray tracing. Based on the first Eulerian Gaussian beam framework proposed in Leung et al. [S. Leung, J. Qian, R. Burridge, Eulerian Gaussian beams for high frequency wave propagation, Geophysics 72 (2007) SM61–SM76], we develop a new Eulerian Gaussian beam method which uses global Cartesian coordinates, level-set based implicit representation and Liouville equations. The resulting method gives uniformly distributed phases and amplitudes in phase space simultaneously. To obtain semi-classical solutions to the Schrödinger equation with different initial wave functions, we only need to slightly modify the summation formula. This yields a very efficient method for computing semi-classical solutions to the Schrödinger equation. For instance, in the one-dimensional case the proposed algorithm requires only O(sNm2) operations to compute s different solutions with s different initial wave functions under the influence of the same potential, where N=O(1/ℏ),ℏ is the Planck constant, and m≪N is the number of computed beams which depends on ℏ weakly. Numerical experiments indicate that this Eulerian Gaussian beam approach yields accurate semi-classical solutions even at caustics. - Article
- Dec 2009
- J COMPUT PHYS

In this paper, we present a method of decomposing a highly oscillatory wave field into a sparse superposition of Gaussian beams. The goal is to extract the necessary parameters for a Gaussian beam superposition from this wave field, so that further evolution of the high frequency waves can be computed by the method of Gaussian beams. The methodology is described for Rd with numerical examples for d=2. In the first example, a field generated by an interface reflection of Gaussian beams is decomposed into a superposition of Gaussian beams. The beam parameters are reconstructed to a very high accuracy. The data in the second example is not a superposition of a finite number of Gaussian beams. The wave field to be approximated is generated by a finite difference method for a geometry with two slits. The accuracy in the decomposition increases monotonically with the number of beams. - Article
- Aug 2009
- WAVE MOTION

The Gaussian beam superposition method is an asymptotic method for computing high frequency wave fields in smoothly varying inhomogeneous media. In this paper we study the accuracy of the Gaussian beam superposition method and derive error estimates related to the discretization of the superposition integral and the Taylor expansion of the phase and amplitude off the center of the beam. We show that in the case of odd order beams, the error is smaller than a simple analysis would indicate because of error cancellation effects between the beams. Since the cancellation happens only when odd order beams are used, there is no remarkable gain in using even order beams. Moreover, applying the error estimate to the problem with constant speed of propagation, we show that in this case the local beam width is not a good indicator of accuracy, and there is no direct relation between the error and the beam width. We present numerical examples to verify the error estimates. - Article
- Feb 1992
- Chaos

The presence of chaos in the classical dynamics is not necessarily the destructive element it was thought to be for semiclassical approximations in the time domain. The method of calculating the semiclassical propagation of initial states and correlation functions for nonlinear and chaotic dynamics is shown, and the excellent accuracy is noted for rather long times. The breakdown timescale is much longer than the infamous "log time" for the cases investigated here. - Article
- Mar 2006
- ACCOUNTS CHEM RES

We discuss several semiclassical Gaussian wave packet approaches with emphasis on one that is not very well-known, that is, the off-center guiding approach. Off-center guiding of (thawed) Gaussian wave packets uses the same information as other Gaussian propagation schemes, that is, the full van Vleck determinant and multiple trajectories. It retains the well-known caustic smoothing property of Gaussians. The off-center guiding of Gaussians can handle hard chaos and other highly nonlinear situations, where the van Vleck propagator, for example, has over 30,000 separate branches. - We apply the level set method to compute the three dimensional multivalued geometrical optics term in the paraxial formulation. The paraxial formulation is obtained from the 3-D stationary eikonal equation by using one of the spatial directions as the artificial evolution direction. The advection velocity field used to move level sets is obtained by the method of characteristics; therefore the motion of level sets is defined in a phase space. The multivalued traveltime and amplitude-related quantity are obtained from solving advection equations with source terms. We derive an amplitude formula in the reduced phase space which is very convenient to use in the level set framework. By using a semi-Lagrangian method in the paraxial formulation, the method has O(N rather than O(N ) memory storage requirement in the five dimensional phase space, where N is the number of mesh points along one direction. Although the computational complexity is still O(MN ), where M is the number of steps in the ODE solver for the semi-Lagrangian scheme, this disadvantage is largely overcome by the fact that up to ) multiple sources can be treated simultaneously. Three dimensional numerical examples demonstrate the e#ciency and accuracy of the method. 1