A novel spectral method for the semi-classical Schrödinger equation based on the Gaussian wave-packet transform
aa r X i v : . [ m a t h . NA ] F e b IMA Journal of Numerical Analysis (2021) Page 1 of 31doi:10.1093/imanum/drnxxx
A novel spectral method for the semi-classical Schr¨odinger equationbased on the Gaussian wave-packet transform B ORUI M IAO † Department of Mathematical Sciences, Tsinghua University, Beijing, 100084, China
AND G IOVANNI R USSO ‡ Department of Mathematics and Computer Science,University of Catania, Catania, 95125, Italy
AND Z HENNAN Z HOU § Beijing International Center for Mathematical Research,Peking University, Beijing, 100871, China [Received on ; revised on ]
In this article, we develop and analyse a new spectral method to solve the semi-classical Schr¨odingerequation based on the Gaussian wave-packet transform (GWPT) and Hagedorn’s semi-classical wave-packets (HWP). The GWPT equivalently recasts the highly oscillatory wave equation as a much lessoscillatory one (the w equation) coupled with a set of ordinary differential equations governing the dy-namics of the so-called GWPT parameters. The Hamiltonian of the w equation consists of a quadratic partand a small non-quadratic perturbation, which is of order O ( √ ε ) , where ε ≪ w equation as a superposition of Hagedorn’s wave-packets, weconstruct a spectral method while the O ( √ ε ) perturbation part is treated by the Galerkin approximation.This numerical implementation of the GWPT avoids imposing artificial boundary conditions and facili-tates rigorous numerical analysis. For arbitrary dimensional cases, we establish how the error of solvingthe semi-classical Schr¨odinger equation with the GWPT is determined by the errors of solving the w equation and the GWPT parameters. We prove that this scheme has the spectral convergence with respectto the number of Hagedorn’s wave-packets in one dimension. Extensive numerical tests are provided todemonstrate the properties of the proposed method. Keywords : semi-classical Schr¨odiner equation; Gaussian wave-packet transform; Hagedrn’s wave-packets; spectral method. † Email: [email protected] ‡ Email: [email protected] § Email: [email protected] © The author 2021. Published by Oxford University Press on behalf of the Institute of Mathematics and its Applications. All rights reserved. of 31
1. Introduction
We consider the time-dependent semi-classical Schr¨odinger equation that governs the dynamics of aquantum particle in an external potential:i ε∂ t ψ ( x , t ) = ˆ H ( x , − i ε∇ x ) ψ ( x , t ) = − ε ∆ x ψ ( x , t ) + V ( x ) ψ ( x , t ) , x ∈ R d , t ∈ R + ; (1.1) ψ ( x , ) = ψ ( x ) , x ∈ R d . (1.2)Here, ε ≪ H denotes the quantum Hamiltonian operator obtained from the classical Hamiltonian symbol H ( x , p ) by the Weyl quantization. V ( x ) in (1.1) is a real-valued scalar potential function and ψ ( x , t ) isa complex-valued wave function. This is one of the fundamental equations in quantum mechanics, andbridges the gap between realms of quantum mechanics and classical mechanics as ε →
0. It also givesrise to many interesting problems in analysis and numerical computation, see e.g. Heller (2018); Jin et al. (2011a); Lasser & Lubich (2020); Zworski (2012).Due to the smallness of the semi-classical parameter ε , the solutions ψ ( x , t ) is highly oscillatory inboth space and time on a scale O ( ε ) (see Jin et al. , 2011a). Also, many physical applications demandaccurate simulations of the semi-classical Schr¨odinger equation (1.1) in high dimensions. Both factslead to significant computational burdens if conventional algorithms are directly applied to equation(1.1)-(1.2). Comparatively speaking, the Fourier spectral method studied in Bao et al. (2002). In theirpaper, it is most efficient among grid-based numerical methods, where the meshing strategy ∆ t = O ( ε ) and ∆ x = O ( ε ) is needed to approximate the wave function. If one only aims to capture the correctphysical observables, the time step size can be relaxed to O ( ) , while the spatial mesh constraint ∆ x = O ( ε ) is still necessary (see Golse et al. (2020) for a rigorous justification). In recent years, this methodhas been extended to the semi-classical Schr¨odinger equations with vector potentials and time-dependentself-consistent field systems, see Jin & Zhou (2013); Ma et al. (2017); Jin et al. (2017).However, to implement the grid-based methods in high dimensions, one has to at least put a fewcollocation points for each wavelength on a scale O ( ε ) in every space dimension. This will causeunaffordable costs in numerical implementation. To overcome the curse of dimensionality, one canapply the splitting methods on a sparse grid (see Gradinaru, 2008) when ε ∼ O ( ) . Unfortunately,the authors pointed out in Lasser & Lubich (2020) that more than O ( ε − d ) points are needed in thesemi-classical regime, making the Fourier spectral method on a sparse grid infeasible.For simulating the semi-classical Sch¨odinger equation in high dimensions, many asymptotic ap-proaches have been developed to avoid resolving the highly oscillatory wave function. For example,the Gaussian beam method (see e.g. Liu et al. (2012); Jin et al. (2008)) reduces the quantum dynamicsto Gaussian beam dynamics that is determined by a set of ordinary differential equations with O ( md ) variables. Here m is the number of Gaussian beams which is determined by the initial condition. Thismethod also brings an approximation error of order O ( ε / ) , which is not satisfactory unless ε is ex-tremely small. High order Gaussian beams (see Liu et al. (2012); Jin et al. (2011b)) can be used toreduce the error of order C k ( T ) ε k / , but this may not lead to higher accuracy for moderate ε either.In particular, for a given ε , it is not guaranteed that increasing k provides a decrease in the error. Forother numerical methods to the semi-classical Schr¨odinger equation, the readers may refer to the recentreviews Jin et al. (2011a); Lasser & Lubich (2020).In Hagedorn (1998), Hagedorn constructed a set of semi-classical wave-packets that constitute abasis of L ( R d ) . When the Hamiltonian operator is quadratic (i.e. the Hamiltonian symbol H ( x , p ) isa second order polynomial in x , p ), the full quantum dynamics is equivalent to the dynamics of the NOVEL SPECTRAL METHOD BASED ON GWPT et al. , 2009) first turned them intoa computational tool (later abbreviated as the HWP method). They use the Galerkin approximation todeal with the general potential, making it feasible for numerical implementation. They also observedthat sparse grid techniques can be applied when choosing the superposition of wave-packets in highdimensions. Compared with Gaussian beam method, the HWP method improved the approximationerror to C ( T ) ε k / (see, e.g., Lasser & Lubich, 2020). For the validity of the Galerkin approximation andthe generalization of HWP to the semi-classical Schr¨odinger equations with vector potentials, see Zhou(2014) for details.Russo and Smereka proposed a new computational framework called the Gaussian wave-packettransform in Russo & Smereka (2013, 2014) (abbreviated as GWPT in later text) to solve the semi-classical Schr¨odinger equation in arbitrary spatial dimensions. It is a reversible transformation andreduces the semi-classical Schr¨odinger equation to the wave-packet dynamics together with the timeevolution of a rescaled wave-function w . The GWPT parameters that describe the wave-packet dynamicssatisfy a set of ordinary differential equations. The rescaled wave function w satisfies a Schr¨odinger-typeequation with a time-dependent Hamiltonian operator, and is less oscillatory than the original solutionto the Schr¨odinger equation (1.1). Another fact about the w equation is that its quantum Hamiltonianconsists of a quadratic part and a O ( √ ε ) perturbation part. This fact motivated us to treat the non-quadratic part as a perturbation. For extension of the GWPT to semi-classical Schr¨odinger equation withvector potential, see Zhou & Russo (2019). In all these works (Russo & Smereka, 2013, 2014; Zhou &Russo, 2019) the GWPT formulation has been implemented as follows. For the w equation, a second orfourth order splitting method is used in time and a Fourier-spectral discretization in a truncated spatialdomain has been adopted, with periodic boundary conditions, which may be a reasonable approximationif the wave function decays fast enough. The GWPT parameters are numerically solved by a fourth orderRunge-Kutta scheme. Thus, the whole scheme is spectrally accurate in space and second or fourth orderaccurate in time, given that the domain truncation is proper.Now we are concerned with the following aspects of the GWPT. First, although this transformationcan effectively reduce the computational burden, its related numerical analysis remains undone. Hereby,we ask how the error in solving the Gaussian wave-packet dynamics and the w equation affects the errorof the reconstructed solution ψ . Second, in Russo & Smereka (2013, 2014) the authors proposed tosolve the w equation with a Fourier spectral method. As we mentioned before, an artificial periodicboundary condition may cause difficulties for the numerical analysis of the method.In this paper, we propose a spectral method based on the GWPT and the HWP method. Recallthat the quantum Hamiltonian operator of the w equation is divided into a quadratic part and a O ( √ ε ) non-quadratic perturbation. We can expand the solution to the w equation as superpositions of Hagedornwave-packets, and thus the time evolution of the quadratic Hamiltonian is recast in terms of wave packet-dynamics. Then, we use the Galerkin approximation to deal with the O ( √ ε ) non-quadratic part. Thisnumerical implementation of the GWPT avoids imposing artificial boundary conditions since it solvesthe w equation on the whole space.For the numerical analysis, we first develop a framework to estimate the error bound in the wavefunction ψ for generic numerical methods which are based on the GWPT. Here, we only need to knowthe L ∞ error when solving the GWPT parameters and discrete L error when solving the w equation to of 31give such an error estimate. Second, we prove an asymptotic error bound with respect to the number ofwave-packets rigorously in one dimension, showing the spectral accuracy. Compared with the previousresults, which are at most algebraic accurate with respect to the number of wave-packets, our proposedscheme is certainly superior and appears promising in high dimensional simulations, since it greatlysaves the computational cost. Extensive numerical tests are provided to demonstrate the properties ofthe method and verify the above arguments.This article is organized as follows: In Section 2 we will briefly review the GWPT and the HWP forthe Schr¨odinger equation and propose a novel spectral method in solving the semi-classical Schr¨odingerequation based on the GWPT. In Section 3 we carry out the related numerical analysis. We providesystematic numerical tests in Section 4, demonstrating the spectral convergence and other propertiesof the proposed scheme. We give some concluding remarks in Section 5. The appendix is devoted toexplain some notation and technical details.
2. Description of the Numerical Method
In this section, we will present a novel spectral method for the semi-classical Schr¨odinger equationbased on the GWPT and the HWP. In particular, the efficiency of the method is not affected by thesmallness of the parameter ε , and the approximation error can be rigorously quantified. To describe themethod in detail, we will first review some properties of the GWPT and HWP.2.1 A Brief Review of the GWPT
For heuristic purposes, let’s first consider the semi-classical Schr¨odinger equation (1.1) with a normal-ized Gaussian initial value ψ ( x , ) = d / ( πε ) − d / exp n ( i / ε ) h(cid:0) x − q (cid:1) T α (cid:0) x − q (cid:1) + p (cid:0) x − q (cid:1) + γ io . (2.1)where q , p are vectors in R d × , γ is a complex scalar and α is a complex-valued d × d matrixwhose imaginary part is positive-definite. Here γ must be chosen such that k ψ ( · , ) k L x =
1. Suchnormalization implies that exp [ − γ , I / ε ] = { det [ α , I ] } / . (2.2)where α I is the imaginary part of the complex matrix α . This Cauchy problem has three typical scales:the O ( ε ) oscillations in the phase, the O ( √ ε ) beam width, and its dynamics over a O ( ) space-and-timedomain. The external scalar potential varies in a O ( ) scale as well. Many asymptotic approaches havebeen extensively studied, while the introduction of the asymptotic approximation error was inevitable.In Russo & Smereka (2013), the authors proposed the Gaussian wave-packet transform (GWPT) toseparate the three scales in the time evolution of highly oscillatory wave function, and the resultingsystem is equivalent to the original semi-classical Schr¨odinger equation.We now investigate this method in detail. Suppose the solution is in the form of ψ ( x , t ) = W ( ξ , t ) exp (cid:2) ( i / ε ) (cid:0) ξ T α R ξ + p T ξ + γ (cid:1)(cid:3) , (2.3)where ξ = x − q ∈ R d × and α R ∈ R d × d is the real part of a complex d × d matrix α . Furthermore sup-pose q , p , γ , α are functions of time whose dynamics are governed by the following ordinary differential NOVEL SPECTRAL METHOD BASED ON GWPT q = p , ˙ p = − ( ∇ V ( q )) T , ˙ γ = p T p − V ( q ) + i ε tr ( α R ) , ˙ α = d ( α R + i α I ) dt = − αα − ∇∇ V ( q ) . (2.4)The initial values for q , p , γ , α can be determined from the initial value (2.1), and are, respectively, q , p , γ , α . Here, q and p are interpreted as the position and momentum center of the Gaussian wavepacket and they follow the Hamiltonian’s equations of classical dynamics. Then the equation for W ( ξ , t ) becomes ∂ t W = i ε ∆ ξ W − ξ T α R ( ∇ ξ W ) T − i ε F ( ξ ) W , (2.5)where F = V ( ξ ; q ) + ξ T α I α I ξ , V ( ξ ; q ) = V ( q + ξ ) − V ( q ) − ∇ V ( q ) ξ − ξ T ∇∇ V ( q ) ξ . Notice that V is the difference between the potential and its second order Taylor expansion around q .The ansatz (2.3) and the ODE system (2.4) separate the O ( ε ) oscillations and ensure that the W functionis concentrated around ξ =
0, i.e. where the position variable x is close to the wave packet center q .Finally, we take the change of coordinates η = B ξ / √ ε , (2.6)where B is a function of time that satisfies the following ODE˙ B = − B α R . (2.7) B = √ α , I . (2.8)It can be shown that α I = B T B for all t >
0, thus the matrix B is invertible at any time. By the change ofcoordinates (2.6), the wave-packets has width O ( ) on η space, thus resolving the O ( √ ε ) scale. Nowwe obtain the equation for w ( η , t ) = W ( ξ ( η , t ) , t ) and its initial valuei ∂ t w = −
12 tr ( B T ∇ η ∇ η wB ) + η T BB T η w + ε U ( η ; q ) w , ˆ H ( t ) w , (2.9) w ( η , ) = d / ( πε ) d / exp ( − η T η ) , (2.10)where U ( η ; q ) = V ( √ ε B − η ; q ) is given by U ( η ; q ) = V ( q + √ ε B − η ) − V ( q ) − √ ε∇ V ( q )( B − η ) − ε ( B − η ) T ∇∇ V ( q ) B − η . (2.11)This Schr¨odinger-typed equation (2.9) and the above ordinary differential equations (2.4) are not oscilla-tory, making them easier to solve than the original Cauchy problem (1.1)-(1.2). Another thing to notice of 31here is that the L η norm of the initial value (2.10) is ε − d / . Besides, this transform can be extended todeal with more general initial conditions, see Russo & Smereka (2013, 2014) for more details.To propose a numerical method based on the GWPT, we have to solve a set of ordinary differentialequations (2.4) and the w equation (2.9) efficiently. In Russo & Smereka (2013), the authors proposeda numerical method based on the GWPT. They implement a Fourier spectral method to solve the w equation (2.9) on some artificial computational domain [ − σ / , σ / ] , which will bring in numericalerrors that are difficult to quantify.In the rest of the article, the parameters q , p , γ , α R , α I , B will be called GWPT parameters . And wewill define the quadratic part ˆ H q ( t ) of the Hamiltonian operator ˆ H ( t ) asˆ H q ( t ) w , −
12 tr ( B T ∇ η ∇ η wB ) + η T BB T η w . (2.12)It is easy to see that both operators ˆ H q ( t ) , ˆ H ( t ) are self-adjoint.2.2 Numerical Method Description
In devising a numerical method based on the GWPT, the main issue concerns the construction of anefficient algorithm for the w equation (2.9). Notice first that the residue term U / ε given by (2.11) is order O ( ε / ) , which can be treated as a perturbation. Therefore, one may expect to use superpositions ofeigenfunctions of ˆ H q ( t ) to represent the solution, and the perturbation term determines the evolution ofthe expansion coefficients, which can be dealt with in a variational way. In Hagedorn (1980), Hagedornproposed a constructive method to identify the eigenfunctions of time-dependent quadratic quantumHamiltonians, which can be viewed as a generalization of the quantum harmonic oscillator theory. Theeigenfunctions can referred to as the semi-classical wave-packets, or Hagedorn’s wave packets. We willfirst review some basic properties of the semi-classical wave-packets.2.2.1 Review of Hagedorn’s Wave-packets (HWP).
To construct such semi-classical wave-packetsrecursively with fixed real vectors q h , p h ∈ R d × and invertible matrices Q h , P h ∈ C d × d that satisfy Q Th P h − P Th Q h = , (2.13) Q ∗ h P h − P ∗ h Q h = E d , × diag ( , , · · · , ) . (2.14)Here Q Th , P Th denote the transpose of Q h , P h and Q ∗ h , P ∗ h denote the conjugate transpose of Q h , P h . Vector q h , p h can be interpreted as the center of the wave-packet on position and momentum space, while thematrices Q h , P h are related to the width of the wave-packets. We add subscript “ h ” to denote the HWPparameters, which are to be distinguished from the GWPT parameters. Based on them, Hagedorn con-structed a set of orthonormal basis functions, also known as the semi-classical wave packets { ϕ δ k } k ∈ N d in L x ( R d ) (see Hagedorn, 1980, 1998) based on a generic quadratic Hamiltonian operator ˆ H . Here δ denotes the semi-classical scale of the HWP, while √ δ is related to the typical width of the HWP. Inparticular, in 1D, those basis functions coincide with the eigenfunction of the Hamiltonian operator. Itis worth noting that H q as in (2.12) is indeed one example of such Hamiltonians. The construction of thebasis functions is carried out by the raising and lowering operators, which are given by A † = ( A † j ) dj = NOVEL SPECTRAL METHOD BASED ON GWPT A = ( A j ) dj = A [ q h , p h , Q h , P h ] , √ δ (cid:2) i Q Th ( ˆ p − p h ) − i P Th ( x − q h ) (cid:3) , (2.15) A † [ q h , p h , Q h , P h ] , √ δ [ − i Q ∗ h ( ˆ p − p h ) + i P ∗ h ( x − q h )] , (2.16)where ˆ p , − i δ∇ x is the momentum operator. The first wave-packet is given by ϕ δ [ q h , p h , Q h , P h ]( x ) = ( πδ ) − d / ( det Q h ) − / exp (cid:18) i2 δ ( x − q h ) T P h Q − h ( x − q h ) + i δ p Th ( x − q h ) (cid:19) , (2.17)while other wave-packets, labeled by multi-index k = ( k , k , · · · , k d ) with non-negative integer entriesare given by ϕ δ k + h j i = p k j + A † j ϕ δ k , (2.18) ϕ δ k −h j i = p k j A j ϕ δ k . (2.19)Here h j i denotes the j -th unit vector in d -dimension and j takes its value in 1 , , · · · , d . ϕ k −h j i is set tobe zero if one of the coordinates in multi-index k − h j i is negative. And it can be proved that { ϕ δ k } k ∈ N d constitute an orthonormal basis of L x ( R d ) , which can be also related to the quantum Hamiltonian oper-ator defined in (2.12). From the above two operators, we can give the following three term recurrencerelations (see, for example, formula (3.28) in Hagedorn (1998)) x − q h = r δ (cid:0) Q h A † + Q h A (cid:1) . (2.20)This formula can be more explicitly written as (see Lasser & Lubich (2020) for detail) Q h (cid:16)p k j + ϕ δ k + h j i ( x ) (cid:17) dj = = r δ ( x − q h ) ϕ δ k ( x ) − Q h (cid:16)p k j ϕ δ k −h j i ( x ) (cid:17) dj = . (2.21)Here (cid:16)p k j + ϕ δ k + h j i ( x ) (cid:17) dj = and q δ ( x − q h ) ϕ δ k ( x ) are regarded as column vectors whose j -th com-ponent are p k j + ϕ δ k + h j i ( x ) and q δ ( x j − q h , j ) ϕ δ k ( x ) respectively. ϕ k −h j i is set to be zero if one of thecoordinates in multi-index k − h j i is negative. Q h denotes the element-wise conjugation of matrix Q h .The recurrence relation provides a practical method to compute the basis ϕ δ k for any multi-index k withfixed δ .In Faou et al. (2009), the authors use finite superposition of { ϕ ε k } k ∈ N d to represent the solution ofthe original equation (1.1) directly, while there is no rigorous analysis in the article. In Zhou (2014), it isshown that the error is order O ( ε k / ) , while k is related to the number of wave-packets. In the followingsections, we will propose the GWPT+HWP method based on the aforementioned methods. of 312.2.2 Formulation of the GWPT+HWP Method–quadratic case.
To explain how the method works,let’s start from a simple case when V ( x ) is a second degree polynomial, so that the residue potential U as in (2.11) vanishes. Then the Hamiltonian operator ˆ H ( t ) in (2.9) can be written in the form ofˆ H ( t ) = ˆ H q ( t ) = (cid:18) − i ∇ x x (cid:19) T (cid:18) BB T
00 4 BB T (cid:19) (cid:18) − i ∇ x x (cid:19) . (2.22)It is proved in Hagedorn (1998) that if we let q h , p h , Q h , P h be time-dependent parameters and introduceanother time-dependent parameter S h ∈ R such that˙ q h = BB T p h , ˙ p h = − BB T q h , ˙ S h = p Th BB T p h − q Th BB T q h , ˙ Q h = BB T P h , ˙ P h = − BB T Q h , (2.23)and suppose further that the initial value can be represented by a convergent series w ( η , ) = e i S h ( ) ∑ k ∈ N d c k ( ) ϕ k [ q h ( ) , p h ( ) , Q h ( ) , P h ( )]( η ) , { c k } k ∈ N d ⊂ C , (2.24)then, the solution to equation (2.9) is given by w ( η , t ) = e i S h ( t ) ∑ k ∈ N d c k ( ) ϕ k [ q h ( t ) , p h ( t ) , Q h ( t ) , P h ( t )]( η ) . (2.25)It can be proved that Q h ( t ) , P h ( t ) are invertible and relation (2.13)-(2.14) holds for all t >
0, whichguarantees the existence of the first wave-packet. In the case when the initial value is in the form of(2.10), we can derive the initial values for (2.23) q h ( ) , p h ( ) = , Q h ( ) = √ E d , P h ( ) = √ E d , S h ( ) = . c ( ) = ε − d / , c k ( ) = , k = . (2.26)which implies that the sum in (2.25) is finite. With a further inspection, we see that only Q h , P h dependon time, while q h , p h , S h , c k can be shown to be constant in time.Therefore, when the external potential V is quadratic, the full semi-classical quantum dynamics areequivalently recast as the dynamics of the GWPT parameters and the HWP parameters, while signifi-cantly less computational cost is needed for simulating the latter systems.2.2.3 Formulation of the GWPT+HWP method–General Case.
To approximate the w equation (2.9)with general potential function, the residue term U in (2.9) won’t vanish. But recalling that semi-classical wave-packets form a basis of L x ( R d ) , we can represent the exact solution w ( η , t ) by an infinitesuperposition of wave-packets with time-dependent expansion coefficients w ( η , t ) = e i S h ( t ) ∑ k ∈ N d c k ( t ) ϕ k [ q h ( t ) , p h ( t ) , Q h ( t ) , P h ( t )]( η ) . (2.27)Here q h , p h , Q h , P h , S h are determined by (2.23) and { c k } k ∈ N d are time-dependent coefficients. To de-termine the dynamics for the coefficients { c k } k ∈ N d from (2.9), we solve an infinite system of ordinarydifferential equations i ε ˙ c k = ∑ l ∈ N d f kl c l , k ∈ N d . (2.28) NOVEL SPECTRAL METHOD BASED ON GWPT { f kl } k , l ∈ N d is an infinite dimensional Hermitian matrix whose elements are determined by f kl ( t ) = Z R d ϕ k ( η , t ) U ( η , t ) ϕ l ( η , t ) d η , ∀ k , l ∈ N d . (2.29)It is clear that solving (2.28)-(2.29) and (2.23) is equivalent to solving the Cauchy problem for (2.9).To construct a practical algorithm, we can truncate the series representation in (2.27), obtaining theapproximation ˜ w ( η , t ) = e i S h ∑ k ∈ K c k ( t ) ϕ k [ q h , p h , Q h , P h ]( η ) , ( K ) < + ∞ , (2.30)and the dynamics of the coefficients are given byi ε ˙ c k = ∑ l ∈ K f kl c l , k ∈ K . (2.31)where F , { f kl } k , l ∈ K is determined in (2.29) and can be approximated by certain quadrature rule.R EMARK F = { f kl } k , l ∈ K can be evaluated by a recurrence-freemethod if the residue U ( η , t ) can be approximated by polynomial functions. Equation (2.20) tells usthat the polynomials can be represented by a linear combination and composition of raising and loweringoperators, thus forming the F matrix can be achieved with far less computations than direct computationof the integral in Eq. (2.29). But we will leave this topic to our next project.R EMARK { c k } k ∈ N d , so q h , p h , S h are still constants. Specifically, q h = p h = ∆ t c , one needs GWPT parameters q , p , γ , α R , α I , B , pa-rameters Q h , P h and evaluations of recurrence relations (2.20) on some finer time grid with step ∆ t gt withhigher order accuracy. Therefore, to update the expansion coefficients { c k } k ∈ K with time step ∆ t c , thethree term recurrence relations (2.21) must be estimated on time step of O ( ∆ t c ) . From the argumentsabove, we can give our algorithm: Algorithm 1
GWPT+HWP method to advance by ∆ t c Update the GWPT parameter q , p , γ , α R , α I , B and parameters Q h , P h by fourth order RK4 with timestep ∆ t gt , from current ti me t to time t + ∆ t c . Apply the three-term recurrence relation (2.20) to calculate { ϕ k ( · , t ) } k ∈ K on quadrature points ontime step of order O ( ∆ t c ) . Update the coefficient vector c = { c k } k ∈ K by solving the ODE:i ε ˙ c = F ( t ) c , by RK4 with time step ∆ t c . F is an Hermitian matrix whose elements are given by (2.29)To conclude, this method uses superposition of semi-classical wave-packets to represent the solution,which allows us to deal with the quadratic and perturbation part respectively. This method can also becarried out with sparse grid techniques in several dimensions, which presents a possibility to solve (2.9)0 of 31in moderately high dimension. What is more, the method only uses wave-packets with fixed width overthe computational domain when implemented with GWPT. This fact facilitates the numerical analysis.We will show that this method can resolve the three typical scales, as mentioned in Section 2.1.R EMARK ∆ t c can be chosen independent of ε , while ∆ t gt needs to be of order O ( ε m ) , where m is the order of the numerical ODE solver. However, since the cost of solving the ODEsystem is negligible compared with the cost of solving the w equation, we conclude that this method canachieve high accuracy while saving computational cost, when ε ≪ Numerical Complexity of GWPT+HWP.
Now we analyze the numerical complexity of Algo-rithm 1, using quadrature rules to evaluate the integral in (2.29). Here we suppose the total number ofquadrature points is N Q . First, we present two possible choices of wave-packets in multiple dimension,both of which can be characterized with a parameter n : K n ∞ = { ϕ k : | k | ∞ n , k ∈ N d } , (2.32) K n = { ϕ k : | k | n , k ∈ N d } . (2.33)We will show later in the article that K n allows higher accuracy for the same computational cost, seeTable 2. So the number of wave-packets is C nd + n in the case of K n , where C mn denotes the binomialcoefficient representing picking m unordered outcomes from n possibilities. We will also show how thischoice of wave-packets affects the complexity.To analyse the numerical complexity, we first recover what sort of operations are involved in eachstep in Algorithm 1. In the first step, the GWPT parameters are calculated. Matrix-vector multiplicationtakes most of the effort in solving the ODE system (2.23) and (2.4). So the overall cost will be O ( N gt d ) ,where N gt = t f / ∆ t gt is the total time steps in solving the parameters q , p , γ , α R , α I , B and Q h , P h .In the second step, we shall assume that the inverse operation when evaluating wave-packets onquadrature points: (cid:16)p k j + ϕ ε k + h j i (cid:17) dj = = Q − h r ε ( x − q h ) ϕ ε k − Q (cid:0)p k j ϕ ε k −h j i (cid:1) dj = ! . will cost O ( d N Q ) . About C nd + n three-term recurrence relations are needed to recover the full wave-packets in K n , so the total calculation will cost O (cid:0) C nd + n d N Q (cid:1) operations.In the last step in each time loop, we need to calculate the Galerkin matrix { f kl } k , l ∈ K given by(2.29). To this purpose, one has to calculate the potential function (2.11) and the wave-packets on atotal number of N Q quadrature points. O ( d N Q ) operations are needed to evaluate the potential oneach quadrature points. Then from the same procedure to estimate the complexity of reconstructingwave-packets on quadrature points, we have to use O (cid:2) C nd + n N Q + ( C nd + n ) N Q (cid:3) operations in the last step.Combining the above arguments, the total computational cost is: O n N c C nd + n d N Q + N c ( C nd + n ) N Q + N gt d o , (2.34)where N c = t f / ∆ t c is the total time step in solving the coefficients { c k } k ∈ K .In the above analysis, the number of C nd + n is reduced largely compared with n d if d ≫ n and N Q is at the same scale as C nd + n . This fact follows from the asymptotes for binomial number by Stirling’s NOVEL SPECTRAL METHOD BASED ON GWPT
11 of 31formula C nd + n ∼ (cid:18) + nd (cid:19) d + ( n + d ) n n !e n ∼ ( n + d ) n n ! , (2.35)which grows polynomially with respect to d . Now supposing N Q grows polynomially in d , then from(2.34), the overall complexity grows polynomially with respect to d . We will show later that this methodwill converge faster with respect to the number of wave-packets, and thus requires far less number ofwave-packets compared with the Fourier spectral method.R EMARK F matrix with improved efficiency. In that scenario, we expect that the complexity can be further reduced.Whereas, how the complexity depends on the dimensionality is still yet to be explored, and we shallwork on that in the future as well.
3. Convergence Analysis of the GWPT+HWP Method
Since there is no rigorous convergence analysis for GWPT-based numerical methods prior to this work,we shall estimate the error of the GWPT+HWP method in two steps. First, we establish an error estimateframework in Section 3.1 for generic GWPT-based methods in assembling the errors in approximatingthe GWPT parameters (2.4) and the w equation (2.9) respectively. Considering the convergence analysisfor the numerical approximation of the GWPT parameters is rather standard, the error estimate boilsdown to quantifying the accuracy of approximating the w equation by the Hagedorn’s wave-packetswhen δ =
1, which we will show to be spectrally accurate with respect to the number of wave-packetsin one dimension in Theorem 3.5.3.1
Error Estimate of the GWPT
We first study how the errors in solving the w equation (2.9)-(2.10) and the GWPT parameters (2.4)affect the error in solving the semi-classical Schr¨odinger equation (1.1)-(1.2).We begin by showing the following identity.L EMMA t >
0, the following identity holds ( det ( α I )) / = exp ( − γ I / ε ) . Proof of Lemma 3.1.
One only has to check that ( det ( α I )) / and exp ( − γ I / ε ) share the same equationand initial value. Their initial values coincide by (2.2).And by chain rule and (2.4), one can derive the equation for exp ( − γ I / ε ) . ddt exp ( − γ I / ε ) = − ε exp ( − γ I / ε )( ε tr ( α R )) = − tr ( α R ) exp ( − γ I / ε ) . Then another direct calculation leads to ddt ( det ( α I )) / =
14 det ( α I ) / ddt ( det ( α I )) = − tr ( α R )( det ( α I )) / . The uniqueness of the solution to an ordinary differential equation gives the desired result. (cid:3) w and q , p , γ , α , B respectively. Furthermore, let˜ w and ˜ q , ˜ p , ˜ γ , ˜ α , ˜ B denote the approximate solutions with error bounds given bymax t ∈ [ , T ] k w − ˜ w k L η E ε − d / , (3.1)max t ∈ [ , T ] max (cid:0) k q − ˜ q k ∞ , k p − ˜ p k ∞ , | γ − ˜ γ | , k α − ˜ α k ∞ , k B − ˜ B k ∞ (cid:1) E . (3.2)For a complex d × d matrix A the standard natural norm is adopted: k A k ∞ = max k v k = k Av k ∞ , v ∈ C d . Here, E denotes the error in solving the GWPT parameters, and E ε − d / is the error in solving the w function, whose L norm is O (cid:0) ε − d / (cid:1) manifesting the magnitude of the initial condition for w as in(2.10).To investigate where the error comes from when rebuilding the wave function by approximate so-lutions ˜ w and ˜ q , ˜ p , ˜ γ , ˜ α , ˜ B , let’s first recall some facts about GWPT. Notice the wave functions ψ , ˜ ψ arereconstructed by ψ ( x , t ) = w ( η ( x ) , t ) exp { ( i / ε ) θ ( x , t ) } , (3.3)˜ ψ ( x , t ) = ˜ w ( ˜ η ( x ) , t ) exp { ( i / ε ) ˜ θ ( x , t ) } . (3.4)where ˜ η ( x ) = ε − / ˜ B ( x − ˜ q ) , (3.5) η ( x ) = ε − / B ( x − q ) , (3.6) θ ( x , t ) = ( x − q ( t )) T α R ( x − q ( t )) + p T ( t )( x − q ( t )) + γ R ( t ) + i γ I ( t ) , (3.7) = θ R + i θ I , ˜ θ ( x , t ) = ( x − ˜ q ( t )) T ˜ α R ( x − ˜ q ( t ) + ˜ p T ( t )( x − ˜ q ( t )) + ˜ γ R ( t ) + i ˜ γ I ( t ) , (3.8) = ˜ θ R + i ˜ θ I . By triangle inequality k ψ − ˜ ψ k L x k w ( η ( x )) exp { ( i / ε ) θ ( x ) } − w ( η ( x )) exp { ( i / ε ) ˜ θ ( x ) }k L x + k w ( η ( x )) exp { ( i / ε ) ˜ θ ( x ) } − w ( ˜ η ( x )) exp { ( i / ε ) ˜ θ ( x ) }k L x + k w ( ˜ η ( x )) exp { ( i / ε ) ˜ θ ( x ) } − ˜ w ( ˜ η ( x )) exp { ( i / ε ) ˜ θ ( x ) }k L x . (3.9)Here t is omitted for simplicity.To deal with the third part of RHS in (3.9), where the error in w equation is dominant, one firstnotices that, by definition, k w ( ˜ η ( x )) exp { ( i / ε ) ˜ θ } − ˜ w ( ˜ η ( x )) exp { ( i / ε ) ˜ θ }k L x = exp ( − ˜ γ I / ε ) k w ( ˜ η ( x )) − ˜ w ( ˜ η ( x )) k L x . Then a change of coordinates leads toexp ( − ˜ γ I / ε ) k w ( ˜ η ( x )) − ˜ w ( ˜ η ( x )) k L x = s ε d det ( ˜ B ) exp ( − ˜ γ I / ε ) r Z R d | w − ˜ w | d ˜ η . (3.10) NOVEL SPECTRAL METHOD BASED ON GWPT
13 of 31Combining (3.5) and (3.6), we have ˜ η = ˜ BB − η + ε − / ˜ B ( q − ˜ q ) . (3.11)And it follows that s ε d det ( B ) exp ( − ˜ γ I / ε ) r Z R d | w − ˜ w | d η exp ( − ˜ γ I / ε )[ det ( α I )] / E . (3.12)Then, we investigate the second part of the RHS in (3.9), where the error is mainly due to thecoordinate transformation (3.5). From (3.10), we have k w ( ˜ η ( x )) exp { ( i / ε ) ˜ θ ( x ) } − w ( η ( x )) exp { ( i / ε ) ˜ θ ( x ) }k L x = s ε d det ( α I ) exp ( − ˜ γ I / ε ) k w ( η ) − w ( ˜ BB − η ) k L η . (3.13)By triangle inequality, we have C ε k w ( η ) − w ( ˜ BB − η ) k L η C ε h k w ( η ) − w ( η + ε − / ˜ B ( q − ˜ q )) k L η + k w ( η + ε − / ˜ B ( q − ˜ q )) − w ( ˜ BB − η + ε − / ˜ B ( q − ˜ q )) k L η i . (3.14)Here C ε equals to q ε d det ( α I ) exp ( − ˜ γ I / ε ) . Strictly speaking, it is not a constant, although it will be veryclose to ε d / if the ODE solver is accurate enough. Provided w ∈ H ( R d ) and k ∇ w k L η = C ( T , d ) ε − d / , ∀ t ∈ [ , T ] , (3.15)we can estimate the first part of RHS in (3.14) as follows C ε k w ( η ) − w ( η + ε − / ˜ B ( q − ˜ q )) k L η C ε k ∇ w ( η ) k L η k ε − / ˜ B ( q − ˜ q ) k C ( T , d ) exp ( − ˜ γ I / ε ) det ( α I ) ( ε − + ) E , (3.16)as we are going to prove below. The first inequality follows from Proposition 9.3 in Brezis (2010).L EMMA u ∈ H x ( R d ) , then there exist a constant C = k ∇ u k L x such that for all h ∈ R d we have k u ( x + h ) − u ( x ) k L x C k h k To estimate the second part of RHS in (3.14) we will first give a lemmaL
EMMA u ∈ S ( R d ) , we have k u ( η + ε − / ˜ B ( q − ˜ q )) − u ( ˜ BB − η + ε − / ˜ B ( q − ˜ q )) k L η C ( T , d ) k η∇ u k L η E . (3.17)Here k η∇ u k L η denotes the quantity k η∇ u k L η , (cid:18) Z R d | ∇ u | | η | d η (cid:19) / . (3.18)4 of 31 Proof of Lemma 3.3.
By Taylor’s expansion with integral reminder, and a simple identity˜ BB − = E d + ( ˜ B − B ) B − , (3.19)we have u ( ˜ BB − η + ε − / ˜ B ( q − ˜ q )) − u ( η + ε − / ˜ B ( q − ˜ q ))= Z h ∇ u ( η + λ ( ˜ B − B ) B − η ) · ( ˜ B − B ) B − η i d λ . Then by Cauchy’s inequality k u ( ˜ BB − η + ε − / ˜ B ( q − ˜ q )) − u ( η + ε − / ˜ B ( q − ˜ q )) k L η (cid:26) Z R d Z (cid:12)(cid:12)(cid:12)h ∇ u ( η + λ ( ˜ B − B ) B − η ) · ( ˜ B − B ) B − η i(cid:12)(cid:12)(cid:12) d λ d η (cid:27) / C ( T , d ) k η∇ u k L η E . Now we have finished the proof of the lemma. (cid:3)
Then, supposing k η∇ w k L η = C ( T , d ) ε − d / , ∀ t ∈ [ , T ] , (3.20)the norm appeared in second part in (3.14) can be estimated as follows by Lemma 3.3 and argument ofdensity, C ε k w ( η + ε − / ˜ B ( q − ˜ q )) − w ( ˜ BB − η + ε − / ˜ B ( q − ˜ q )) k L η C ( T , d ) exp ( − ˜ γ I / ε ) det ( α I ) E . Combining the above arguments, we have k w ( η ( x )) exp { ( i / ε ) ˜ θ ( x ) } − w ( ˜ η ( x )) exp { ( i / ε ) ˜ θ ( x ) }k L x C ( T , d ) exp ( − ˜ γ I / ε ) det ( α I ) ( ε − + ) E . (3.21)Now we investigate the first part of the inequality in (3.9), where phase error is dominant. First,notice that the phase error can be expanded as follows: | exp { ( i / ε ) θ } − exp { ( i / ε ) ˜ θ }| = ( exp ( − γ I / ε ) − exp ( − ˜ γ I / ε )) + exp [ − ( γ I + ˜ γ I ) / ε ] { − exp [( i / ε )( θ R − ˜ θ R )] − exp [( i / ε )( ˜ θ R − θ R )] } . (3.22)where exp [( i / ε )( ˜ θ R − θ R )] is given by a very long formula:exp [( i / ε )( ˜ θ R − θ R )] = exp (cid:26) i ε h ( x − q ) T α R ( q − ˜ q ) + ( q − ˜ q ) T α R ( x − q ) + ( q − ˜ q ) T α R ( q − ˜ q )+ ( x − q ) T ( ˜ α R − α R )( x − q ) + ( q − ˜ q ) T ( ˜ α R − α R )( x − q )+ ( x − q ) T ( ˜ α R − α R )( q − ˜ q ) + ( q − ˜ q ) T ( ˜ α R − α R )( q − ˜ q )+ ( ˜ p − p ) T ( x − q ) + p T ( q − ˜ q ) + ( ˜ p − p ) T ( q − ˜ q ) + ( γ R − ˜ γ R ) i(cid:27) . (3.23) NOVEL SPECTRAL METHOD BASED ON GWPT
15 of 31Since one can write exp [ − ( γ I + ˜ γ I ) / ε ] { − exp [( i / ε )( θ − ˜ θ )] − exp [( i / ε )( ˜ θ − θ )] } = exp ( − ( γ I + ˜ γ I ) / ε )( − [( i / ε )( θ R − ˜ θ R ))] ε exp ( − ( γ I + ˜ γ I ) / ε ) | θ R − ˜ θ R | . By (3.22) and triangle inequality, one has k w ( η ( x )) exp { ( i / ε ) θ } − w ( η ( x )) exp { ( i / ε ) ˜ θ }k L x k w | exp ( − γ I / ε ) − exp ( − ˜ γ I / ε ) |k L x + ε k w exp ( − ( γ I + ˜ γ I ) / ε ) | θ R − ˜ θ R |k L x . Combining (3.23) and supposing E is sufficiently small, we have, | θ R − ˜ θ R | C ( α R , q , p ) (cid:0) E ( x − q ) T ( x − q ) + E | x − q | + E (cid:1) . (3.24)Here C ( α R , q , p ) is a function of the GWPT parameters α R , q , p and in turn is a function of time.This estimate follows from (3.23) by investigating the terms and omitting the higher order term in E .Therefore, by formal calculation, we have k w ( η ( x )) exp { ( i / ε ) θ } − w ( η ( x )) exp { ( i / ε ) ˜ θ }k L x (cid:12)(cid:12) exp ( − γ I / ε ) − exp ( − ˜ γ I / ε ) (cid:12)(cid:12) s ε d det ( α I ) k w k L η + C ( α R , q , p ) exp [ − ( ˜ γ I + γ I ) / ( ε )]( E / ε ) (cid:16) k| x − q | w ( η ( x )) k L x + k| x − q | w ( η ( x )) k L x + k w ( η ( x )) k L x (cid:17) . where C ( α R , q , p ) is a time-dependent constant. Let M be the fourth moment of ψ defined by M , max t ∈ [ , T ] h ψ ( t ) , | x − q ( t ) | ψ ( t ) i . (3.25)Then we have, by Cauchy’s inequality: k w ( η ( x )) exp { ( i / ε ) θ } − w ( η ( x )) exp { ( i / ε ) ˜ θ }k L x (cid:12)(cid:12) exp ( − γ I / ε ) − exp ( − ˜ γ I / ε ) (cid:12)(cid:12) s ε d det ( α I ) k w k L η + C ( α R , q , p ) exp [ − ( ˜ γ I − γ I ) / ( ε )]( E / ε )( √ M + √ M + ) The L ( R d ) norm of w is a constant ε − d / supposing the initial wave function is normalized. Noticeanother trivial inequality 2 √ M √ M +
1, the estimate becomes: k w ( η ( x )) exp { ( i / ε ) θ } − w ( η ( x )) exp { ( i / ε ) ˜ θ }k L x (cid:12)(cid:12) exp ( − γ I / ε ) − exp ( − ˜ γ I / ε ) (cid:12)(cid:12) det ( α I ) − / + C ( α R , q , p ) exp ( − ( ˜ γ I − γ I ) / ε )( E / ε )( √ M + ) . E = O ( ε ) , it can be guaranteed that: | exp ( − ( γ I − ˜ γ I ) / ε ) − | A E / ε , (3.26)exp ( − ( ˜ γ I − γ I ) / ε ) A . (3.27)where A , A are constants independent of ε , but may depend on T . Now combining this statement andLemma 3.1, we haveT HEOREM q , p , γ , α , B be the exact solution to the ordinary differential equations (2.4) and w be the exact solution to the w equation (2.9). Furthermore suppose the w equation and the GWPTparameters are approximated by ˜ w and ˜ q , ˜ p , ˜ γ , ˜ α , ˜ B with error as follows k w − ˜ w k L η E ε − d / , (3.28)max (cid:0) k q − ˜ q k ∞ , k p − ˜ p k ∞ , | γ − ˜ γ | , k α − ˜ α k ∞ , k B − ˜ B k ∞ (cid:1) E ∼ O ( ε ) . (3.29)If for t ∈ [ , T ] the fourth moment is finite M = h ψ , | x − q | ψ i < + ∞ , t ∈ [ , T ] . and condition (3.1), (3.2), (3.15), (3.20), (3.26) and (3.27) hold. Then the following estimate holds k ψ − ˜ ψ k L x A E + A E / ε + A ( α R , q , p )( E / ε ) √ M . (3.30)where A ( α R , q , p ) is a function of α R , q , p , and therefore depends on time. A , A are positive constantswhich are only dependent on T and dimensionality d .3.2 The Error of Solving (2.9)
Using GWPT+HWP, 1d case
Theorem 3.1 tells us that we can control the error when solving the semi-classical Schr¨odinger equation(1.1) with the GWPT. The error of ODE solvers can be given by standard analysis, while the error ofsolving (2.9), (2.11) with HWP requires some analysis, which will be performed in this section.An available framework to prove the spectral convergence is given in Pasciak (1980). To apply thisframework, we have to give approximating property w.r.t. the number of some certain basis, e.g., Fouriermodes in Pasciak (1980). We will present an approximation theorem of semi-classical wave-packets insection 3.2.2.Then, combining the above arguments, we will prove the a priori error bound in solving the w equation (2.9) by the HWP method.3.2.1 General formula of P h , Q h in 1d. As we mentioned in Section 2.2, HWP parameters Q h , P h satisfy the following sets of differential equations in 1d, where scalars are commutative,˙ Q h = α I P h , ˙ P h = − α I Q h . (3.31)Furthermore, one notice that Q h , P h satisfy the relation (2.14), which in 1d becomes Q h P h − P h Q h = . (3.32) NOVEL SPECTRAL METHOD BASED ON GWPT
17 of 31From (2.26), one can derive the initial value for (3.31), Q h ( ) = √ , P h ( ) = i √ . (3.33)Now it is well known that Q h = √ (cid:18) Z t α I ( s ) ds (cid:19) , P h = √
2i exp (cid:18) Z t α I ( s ) ds (cid:19) . (3.34)This turns out to be a very good property, which guarantees the semi-classical wave-packets given by(2.17) and recurrence relations have fixed width.3.2.2 Approximation with wave-packets after the GWPT.
Now we can state and prove the approxi-mating properties of semi-classical wave-packets after the GWPT. The key observation is that Hagedornwave-packets in 1d are rescaled Hermite functions.(see Hagedorn, 1998)And it is indeed proved in Hagedorn (1998) that for all k ∈ N , Hagedorn’s wave-packet can beexplicitly represented by ϕ δ n [ , , Q h , P h ]( x ) = − n / ( n ! ) − / F n ( Q h , δ , x ) ϕ δ [ , , Q h , P h ]( x ) , (3.35)where F n ( Q h , δ , x ) = Q − n / h ( Q h ) n / H n ( δ − / | Q h | − x ) , (3.36) H n ( x ) = e x (cid:18) − ∂∂ x (cid:19) n e − x , (3.37) ϕ δ [ , , Q h , P h ]( x ) = ( πδ ) − / Q − / h exp (cid:18) i2 δ P h Q − h x (cid:19) . (3.38)So we can expect some approximations similar to the ones appearing in the theorems given in Shen et al. (2011). The final goal is to prove similar properties of Hagedorn’s wave-packets. In the followingpart, we will take the weight to be ω = exp ( −| Q h | − x ) and F n ( x ) = F n ( Q h , , x ) . We first review somebasic properties of Hermite functions and Hermite polynomials, whose proof can be found in, e.g.,section 7.2.1 of Shen et al. (2011).L EMMA ∂ kx F n ( x ) = k n ! Q kh ( n − k ) ! F n − k , (3.39) Z R H m H n exp ( − x ) dx = √ π n n ! δ mn , (3.40) Z R ∂ kx F m ( x ) ∂ kx F n ( x ) exp ( −| Q h | − x ) dx = h n , k δ mn , (3.41)where δ mn is the Kronecker symbol, and h n , k = k n ! | Q h | k − ( n − k ) ! √ π n n !8 of 31To state the first theorem, we will present the definition of weight L ω spaces. Definition 1
The space of weighted L x functions, denoted by L x , ω is given by L x , ω , (cid:26) f ∈ L loc : Z R | f | ω dx < ∞ (cid:27) , (3.42)where ω > x denotes the integral variable. It is a Hilbert space whose normis given by k f k L x , ω , Z R | f | ω ( x ) dx . (3.43)Similarly, we can define the weighted Sobolev spaces by H mx , ω , n f ∈ L loc : ∂ β x f ∈ L x , ω , ∀| β | m , β ∈ N d o . (3.44)And its corresponding norm k · k H mx , ω is given by k f k H mx , ω , ∑ | β | m k ∂ β x f k L x , ω . (3.45)T HEOREM u ∈ S ( R ) with 0 l m N +
1. Then we have: k ∂ lx ( Π N u − u ) k L x , ω | Q h | m − l ( m − l ) / s ( N − m + ) ! ( N − l + ) ! k ∂ mx u k L x , ω . (3.46)Here Π N is the projection onto the space spanned by the first N + Z R ( u − Π N u ) v N ω dx = , ∀ v N ∈ P N = span (cid:8) H i ( | Q h | − x ) (cid:9) Ni = . (3.47)and the weighted L x norm is given by (3.42). S ( R ) denotes the Schwartz space. Proof.
It is obvious that for u ∈ L x , ω , u ( x ) = ∞ ∑ n = ˆ u n F n ( x ) , with ˆ u n = h u , F n i x , ω h F n , F n i x , ω . Then one has k ∂ lx ( Π N u − u ) k L x , ω = ∞ ∑ n = N + h n , l | ˆ u n | max n > N + (cid:26) h n , l h n , m (cid:27) ∞ ∑ n = N + h n , m | ˆ u n | . where ˆ u n are the expansion coefficients. By the following fact,max n > N + (cid:26) | Q h | m − l m − l ( n − m ) ! ( n − l ) ! (cid:27) = | Q h | m − l m − l ( N − l + )( N − l ) · · · ( N − m + ) . (3.48) NOVEL SPECTRAL METHOD BASED ON GWPT
19 of 31We have, k ∂ lx ( Π N u − u ) k L x , ω h N + , l h N + , m k ∂ mx u k L x , ω . (3.49) (cid:3) From the estimates for the approximation error in weighted norm, we can give the error from theraising operator. First we will give a useful lemma representing the lowering operator in 1d.L
EMMA u ∈ S ( R ) and fixed positive integer l . Suppose further that p h = q h =
0, we haveexp (cid:18) i2 P h Q − h x (cid:19) ∂ lx (cid:20) u exp (cid:18) − i2 P h Q − h x (cid:19)(cid:21) = l / Q − l / h A l u . (3.50)where A is the lowering operator given in (2.15) when δ = Proof.
The proof is given by induction on l . For l = l = n −
1. By substituting u for 2 ( n − ) / Q − ( n − ) / h A n − u and again apply the above equality when l =
1, we have:exp (cid:18) i2 P h Q − h x (cid:19) ∂ x (cid:26) exp (cid:18) − i2 P h Q − h x (cid:19) exp (cid:18) i2 P h Q − h x (cid:19) ∂ n − x (cid:20) u exp (cid:18) − i2 P h Q − h x (cid:19)(cid:21)(cid:27) = exp (cid:18) i2 P h Q − h x (cid:19) ∂ nx (cid:20) u exp (cid:18) − i2 P h Q − h x (cid:19)(cid:21) = n / Q − n / h A (cid:0) A n − u (cid:1) . we can obtain the result for any positive integer l . (cid:3) Note that (2.26) and Remark 2.2 justify the special choice of p h , q h . Now we can state and prove anapproximation result related to the lowering operators:T HEOREM u ∈ S ( R ) , with 0 l m N +
1, we have: (cid:13)(cid:13)(cid:13) A l ( ˆ Π N u − u ) (cid:13)(cid:13)(cid:13) L x s ( N − m + ) ! ( N − l + ) ! k A m u k L x . (3.51)Here ˆ Π N is the projection operator onto the space spanned by the first N + Π N u , exp (cid:18) i P h x Q h (cid:19) Π N (cid:18) u exp (cid:18) − iP h x Q h (cid:19)(cid:19) , ∀ u ∈ S ( R ) , (3.52)and A is the lowering operator given by (2.15) when ε = Proof.
From the above lemma, we have k A l ( ˆ Π N u − u ) k L x = (cid:18) | Q h |√ (cid:19) l (cid:13)(cid:13)(cid:13)(cid:13) exp (cid:18) i2 P h Q − h x (cid:19) ∂ lx (cid:20) ( ˆ Π N u − u ) exp (cid:18) − i2 P h Q − h x (cid:19)(cid:21)(cid:13)(cid:13)(cid:13)(cid:13) L x . (3.53)Furthermore, notice another equality, which can be found in Faou et al. (2009). ℑ ( P h Q − h ) = ( Q h Q ∗ h ) − = | Q h | − , (3.54)0 of 31that leads to ω = exp ( −| Q h | − x ) . (3.55)Together with (3.53), which tells us k A l ( ˆ Π N u − u ) k L x = (cid:18) | Q h |√ (cid:19) l (cid:13)(cid:13)(cid:13)(cid:13) ∂ lx (cid:20) ( ˆ Π N u − u ) exp (cid:18) − i2 P h Q − h x (cid:19)(cid:21)(cid:13)(cid:13)(cid:13)(cid:13) L x , ω = (cid:18) | Q h |√ (cid:19) l (cid:13)(cid:13)(cid:13)(cid:13) ∂ lx (cid:26) ( Π N − Id ) (cid:20) u exp (cid:18) − i P h x Q h (cid:19)(cid:21)(cid:27)(cid:13)(cid:13)(cid:13)(cid:13) L x , ω . Now using Theorem 3.2, we have k A l ( ˆ Π N u − u ) k L x (cid:18) | Q h |√ (cid:19) m s ( N − m + ) ! ( N − l + ) ! (cid:13)(cid:13)(cid:13)(cid:13) ∂ mx (cid:20) u exp (cid:18) − i P h x Q h (cid:19)(cid:21)(cid:13)(cid:13)(cid:13)(cid:13) L x , ω . Again by Lemma 3.5 k A l ( ˆ Π N u − u ) k L x s ( N − m + ) ! ( N − l + ) ! k A m u k L x . Now the conclusion immediately follows. (cid:3)
From Theorem 3.2 and Theorem 3.3, we can state and prove the approximation error bound afterthe GWPTT
HEOREM A be the lowering operator defined by (2.15). If u ∈ S ( R ) , k = , , · · · , m and3 m N +
1, then k ∂ lx ( ˆ Π N u − u ) k L x C s ( N − m + ) ! ( N − l + ) ! k A m u k L x , l = , , , . (3.56)where C is a positive constant independent of m , N and u , but dependent on l . Here ˆ Π N is the projectionoperator defined by (3.52). Proof.
We will prove this theorem when Q h , P h satisfy the general formula (3.34). Here we notice auseful fact that ω = exp ( − x ) .When l =
0, this reduces to the previous theorem.When l =
1, we can prove a lemma similar to Formula (B.36b) in Shen et al. (2011) using integrationby parts k xu k L x | Q h | − | Q h | k u k H x , ω = k u k H x , ω . (3.57)And by a direct calculation ∂ x ( ˆ Π N u − u ) = exp (cid:0) − x (cid:1) ∂ x (cid:2) Π N (cid:0) u exp (cid:0) x (cid:1)(cid:1) − u exp (cid:0) x (cid:1)(cid:3) − x exp (cid:0) − x (cid:1) (cid:2) Π N (cid:0) u exp (cid:0) x (cid:1)(cid:1) − u exp (cid:0) x (cid:1)(cid:3) . NOVEL SPECTRAL METHOD BASED ON GWPT
21 of 31Hence by (3.57), we have k ∂ x ( ˆ Π N u − u ) k L x (cid:13)(cid:13) ∂ x (cid:2) Π N ( u exp ( x )) − u exp ( x ) (cid:3)(cid:13)(cid:13) L x , ω + (cid:13)(cid:13) x (cid:2) Π N ( u exp ( x )) − u exp ( x ) (cid:3)(cid:13)(cid:13) L x , ω C k Π N ( u exp ( x )) − u exp ( x ) k H x , ω = C ( (cid:13)(cid:13) ∂ x (cid:2) Π N ( u exp ( x )) − u exp ( x ) (cid:3)(cid:13)(cid:13) L x , ω + (cid:13)(cid:13)(cid:2) Π N ( u exp ( x )) − u exp ( x ) (cid:3)(cid:13)(cid:13) L x , ω ) . Where C is some positive constant independent of N .By Theorem 3.2 and the explicit formula of Q h , we can estimate the two parts in RHS respectively. k ∂ x ( ˆ Π N u − u ) k L x C ( m − ) r ( N − m + ) ! N ! k ∂ mx ( u exp ( x )) k L x , ω + C m s ( N − m + ) ! ( N + ) ! k ∂ mx ( u exp ( x )) k L x , ω C m r ( N − m + ) ! N ! k ∂ mx ( u exp ( x )) k L x , ω = C m r ( N − m + ) ! N ! k exp ( − x ) ∂ mx ( u exp ( x )) k L x . The third equality holds from definition. Then we apply (3.50) that statesexp ( − x ) ∂ mx (cid:2) u exp ( x ) (cid:3) = m exp (cid:18) − m Z t α I ( s ) ds (cid:19) A m u . This completes the proof when l = l = , (cid:3) Considering the argument of density, the error estimates of the spectral interpolations as in The-orem 3.4 also hold for functions in certain Sobolev spaces. Such estimates are crucial for the erroranalysis for solving the w equation to be presented in the next section.3.2.3 Dynamics for the truncated GWPT+HWP method.
To present the final theorem, we first intro-duce some notation. Let w N ( η , t ) be the solution of the equation: i ∂ t w N = ˆ Π N H ( t ) w N w N ( η , ) = ˆ Π N w where H ( t ) is defined in the w equation. We have to give another lemma concerning the dynamics ofHWP method.L EMMA K = { , , , · · · , N } and suppose the GWPT parameters ( q , p , γ , α R , α I , B ) , HWPparameters ( Q h , P h ) and the coefficients ( { c k } k ∈ K ) are solved exactly. Then the approximate solution˜ w of w equation, where ˜ w is given by˜ w ( η , t ) = ∑ k ∈ K c k ( t ) ϕ k [ , , Q h , P h ]( η ) , satisfies i ∂ t ˜ w = ˆ Π N ˆ H ( t ) ˜ w . (3.58)2 of 31Here ˆ Π N is the projection operator onto the space spanned by the first N + H ( t ) is the quantum Hamiltonian operator defined in (2.9). Proof.
We do this by direct calculation.i ∂ t ˜ w = N ∑ k = i ( ∂ t c k ) ϕ k + N ∑ k = c k i ∂ t ϕ k = N ∑ k = N ∑ l = f kl c l ! ϕ k + N ∑ k = c k ˆ H q ϕ k = N ∑ k = N ∑ l = f kl c l ! ϕ k + N ∑ k = c k ( k + ) α I ϕ k . The last equality holds since Hagedorn wave-packets are the eigenfunctions of the quadratic Hamilto-nian ˆ H q whose corresponding eigenvalues are α I ( k + ) (see Remark 2.4 in Hagedorn (1998)). (cid:3) Now we can prove the spectral convergence of the GWPT+HWP method at given time T :T HEOREM w = w ( t ) be the solution of the equation (2.9) and ˆ Π N w ( t ) − w ( t ) ∈ H x ( R ) up totime T . Suppose A m w ( t ) ∈ L x ( R ) , V ( ) ∈ L ∞ x ( R ) and the stability condition: k A m w ( t ) k L η B ( m , T ) k A m w ( ) k L η for t ∈ [ , T ] . (3.59)Then there exist constants B ( m , T ) , B ( m , T ) independent of w such that k w ( t ) − w N ( t ) k L η B ( m , T ) s ( N − m + ) ! ( N + ) ! k A m w ( ) k L η + B ( m , T ) T ε / ˜ α / I s ( N − m + ) ! ( N − ) ! k A m w ( ) k L η . (3.60)Here ˜ α I = min t T α I and V ( ) is the third derivative of the potential function. Proof.
The idea is simply an extension of a technique introduced in the proof in Theorem 2 in Pasciak(1980). Let X ( t ) = ˆ Π N w ( t ) and v = w ( t ) − X ( t ) , where ˆ Π N is a projection operator that depends on time.We first calculate the time derivative of X ( t ) . Here we abbreviate the operator ˆ H ( t ) , ˆ H q ( t ) as ˆ H , ˆ H q . ∂ t ( ˆ Π N ( t ) w ( t )) = ∂ t N ∑ k = h w , ϕ k i ϕ k ! = N ∑ k = h
1i ˆ Hw , ϕ k i ϕ k + N ∑ k = h w ,
1i ˆ H q ϕ k i ϕ k + N ∑ k = h w , ϕ k i
1i ˆ H q ϕ k = N ∑ k = h
1i ˆ Hw , ϕ k i ϕ k − N ∑ k = h
1i ˆ H q w , ϕ k i ϕ k + N ∑ k = h w , ϕ k i
1i ˆ H q ϕ k . If θ ∈ S N ( t ) , where S N ( t ) is the linear space spanned by the first N + { ϕ k } Nk = , we then have h ∂ t ( ˆ Π N w ( t )) , θ i = h
1i ˆ Π N ˆ Hw , θ i + h ( ˆ H q ˆ Π N − ˆ Π N ˆ H q ) w , θ i , that leads to h ∂ t ( ˆ Π N w ( t )) , θ i − h
1i ˆ H ˆ Π N w ( t ) , θ i = h ε (cid:0) ˆ Π N U − U ˆ Π N (cid:1) w ( t ) , θ i . (3.61) NOVEL SPECTRAL METHOD BASED ON GWPT
23 of 31Here U denotes the non-quadratic part of the potential function defined in (2.11). Further by Lemma3.6, we have h ∂ t w N ( t ) , θ i − h
1i ˆ Hw N ( t ) , θ i = . (3.62)Let e = w N ( t ) − X ( t ) ∈ S N ( t ) , we have, by subtracting (3.61) from (3.62) and setting θ = e : h ∂ t e , e i − h
1i ˆ He , e i = h ε (cid:0) − ˆ Π N U + U ˆ Π N (cid:1) w , e i = h ε UX − Uw , e i . (3.63)The last equality holds by the self-adjointness of the operator ˆ Π N . In the same manner, we have: h e , ∂ t e i − h e ,
1i ˆ He i = h e , ε UX − Uw i . (3.64)Thus we have, adding the above two equalities (3.63), (3.64) and notice that − i ˆ H is skew-symmetric,we have ∂ t k e k L η = ℜ h ε ( UX − Uw ) , e i ε k e k L η k U ( ˆ Π N w − w ) k L η . We will give a lemma to prove the boundedness of U ,L EMMA Π N w − w ∈ H ( R ) up to time T . And we suppose V ( ) ∈ L ∞ ( R ) , then we have theestimate: (cid:13)(cid:13) U ( ˆ Π N w − w ) (cid:13)(cid:13) L η C (cid:18) εα I (cid:19) / (cid:13)(cid:13) Π N w exp ( η ) − w exp ( η ) (cid:13)(cid:13) H η , ω . (3.65) Proof.
We notice the Taylor expansion with integral remainder U = (cid:18) εα I (cid:19) / η Z ( − θ ) V ( ) (cid:18) q + θη r εα I (cid:19) d θ . (3.66)So we have, by the assumption that V ( ) ∈ L ∞ k U ( ˆ Π N w − w ) k L η C (cid:18) εα I (cid:19) / k η ( ˆ Π N w − w ) k L η = C (cid:18) εα I (cid:19) / k η ( Π N − Id ) w exp ( η ) k L η , ω . where C is a constant that only depends on the potential function V . To give a bound for the RHS, wecan apply the technique applied in Theorem 3.4. Thus we have k η ( Π N w exp ( η ) − w exp ( η )) k L η , ω C k Π N w exp ( η ) − w exp ( η ) k H η , ω . So we have finished the proof of the lemma. (cid:3) ∂ t k e k L η C ε / α / I k Π N w exp ( η ) − w exp ( η ) k H η , ω C ε / α / I ∑ l = k A l ( ˆ Π N w − w ) k L η C ε / α / I ∑ l = s ( N − m + ) ! ( N − l + ) ! k A m w k L η . By the stability condition (3.59), we have ∂ t k e k L η C ε / α / I s ( N − m + ) ! ( N − ) ! k A m w k L η C ( m , T ) ε / ˆ α / I s ( N − m + ) ! ( N − ) ! k A m w ( ) k L η . Here ˆ α I = min t T α I and C ( m , T ) is a constant that depends on m and T . The above inequality infers k e k L η C ( m , T ) T ε / ˆ α / I s ( N − m + ) ! ( N − ) ! k A m w ( ) k L η . (3.67)From Theorem 3.4, we have: k w − w N k L η k w − ˆ Π N w k L η + k ˆ Π N w − w N k L η C s ( N − m + ) ! ( N + ) ! k A m w ( t ) k L η + e ( t ) . Again by the stability condition (3.59) and (3.67), k w − w N k L η C ( m , T ) s ( N − m + ) ! ( N + ) ! k A m w ( ) k L η + C ( m , T ) T ε / ˜ α / I s ( N − m + ) ! ( N − ) ! k A m w ( ) k L η . Now we have finished the proof of Theorem 3.5. (cid:3)
4. Numerical Tests
In the following section, we will perform several numerical tests for 1d and 2d cases to provide imme-diate evidence of the effectiveness of the method, and to verify the numerical analysis.4.1
Example 1, 1d case
In this example (see Russo & Smereka, 2013), we choose d = V ( x ) = − cos ( x ) . Theinitial value is chosen to be ψ d ( x , ) = (cid:18) πε (cid:19) / exp (cid:20) − ( x − π / ) ε (cid:21) . NOVEL SPECTRAL METHOD BASED ON GWPT
25 of 31 -8 -6 -4 -2 0 2 4 6 8 η -1-0.500.511.522.53 GWPT+HWP, w equation -0.5 0 0.5 1 1.5 x -2-1.5-1-0.500.511.522.53 Reference solutionGWPT+HWPGWPT, W(x) F IG . 1. We choose ε = / ∆ t c = / ∆ t gt = / t f = . w equation (2.9) on quadrature points is plotted in η space at final time. The cross represent the computed solutionon the grid point, while the dotted line represents the computed solution on {− π + π ( k − ) / } k = . Right panel: thereal part of the reconstructed wave function and the real part of the reference solution are plotted. The W function is alsoplotted on quadrature points transformed to the real space. The dark green dotted line represent the computed solution of W on {− π + π ( k − ) / } k = . The wave function on {− π + π ( k − ) / } k = is also reconstructed and is represented by bluedots. It is clear that our method requires fewer points in real space. Reference solution for ψ : SP4 on ( − π , π ) with ∆ t REF = ε / ∆ x REF = πε / where q = π / , p = , α = i , γ = , B = , Q h , = / √ , P h , = i √
2. Unless otherwise specified,the reference solution is obtained by the Fourier spectral method applied to the Schr¨odinger equation(1.1), using the SP4 method proposed in Chin & Chen (2001), with uniform grid in [ − π , π ) with spacestep ∆ x REF = πε /
64 and time step ∆ t REF = ε / ψ , the wave function w is notoscillatory, neither in space nor in time. Indeed we find out that it requires much fewer grid points onthe computational domain, and a much larger time step ∆ t c , independently on the small scale ε , stillreaching high accuracy.Second, to determine how one chooses time step ∆ t c and ∆ t gt to achieve best approximation, we plotthe L error for various time step ∆ t c and ∆ t gt with 30 wave-packets. The result is shown in Figure 2. Inleft panel, for fixed ∆ t c = / ∆ t gt at final time t f = . ∆ t gt = / L error for different values of ∆ t c at final time t f = . ∆ t gt is small enough, ∆ t c can be chosen independent of ε . Third, we will verify the spectral convergence with respect to the number of wave-packets n . Wetake various values of n and let the quadrature points be N Q =
30 to achieve the best approximation ofGalerkin matrix. The time step is chosen to be ∆ t c = / ∆ t gt = / L error at t f = . ∆ t c = / ∆ t gt = / n =
10 wave-packets.The result is shown in the right panel of Figure 3. This shows the Galerkin matrix is solved with spectral6 of 31 -4 -3 ∆ t gt -14 -13 -12 -11 -10 -9 -8 × ( ∆ t gt ) -3 -2 -1 ∆ t c -14 -12 -10 -8 -6 -4 ( ∆ t c ) /100 F IG . 2. Left panel: We fix ∆ t c = /
128 and compare the L error for different values of ∆ t gt at final time t f = . n =
30 wave-packets in each test. N Q = ∆ t gt = / L error for different values of ∆ t c at final time t f = . ψ : SP4 on ( − π , π ) with ∆ t REF = ε /
64 and spatial grid size ∆ x REF = πε / n: number of wave-packets -14 -12 -10 -8 -6 -4 -2 N Q -12 -10 -8 -6 -4 -2 F IG . 3. Left panel: For fixed time step size ∆ t c = / ∆ t gt = / t f = . L error ofGWPT+HWP in 1d with different values of wave-packets n and ε . 30-point Gauss-Hermite quadrature rule is used to evaluate theGalerkin matrix. Right panel: For fixed time step satisfying ∆ t c = / ∆ t gt = / t f = . L error of GWPT+HWP in 1d with different values of quadrature points N Q and ε . We choose 10 wave-packets in each test.Reference solution for ψ : SP4 on ( − π , π ) with ∆ t REF = ε /
64 and spatial grid size ∆ x REF = πε /
64. NOVEL SPECTRAL METHOD BASED ON GWPT
27 of 31 n CPU time ε Table 1. For different values of ε , we compare the CPU time for different numbers of wave-packets n with fixed time step size ∆ t c = / ∆ t gt = / t f = n + accuracy with respect to the number of quadrature points N Q . The results also tell us how to choose thenumber of quadrature points properly to reduce computational cost.In the following part, we test the relation between α I and the error of GWPT method. We will fix n =
30 and time step ∆ t c = / ∆ t gt = / N Q =
50 whenevaluating the Galerkin matrix. The result is shown in Figure 4. As we can see from the above figure, t -15 -10 -5 × α I F IG . 4. For fixed time step ∆ t c = / ∆ t gt = / L error with α I from t = t =
18. 100-pointquadrature rule is used to evaluate the Galerkin matrix. Reference solution for ψ : SP4 on ( − π , π ) with ∆ t REF = ε /
64 and spatialgrid size ∆ x REF = πε / the error has a sudden increase when α I is small. This behavior is typical in the GWPT, as mentioned inRusso & Smereka (2013, 2014). This phenomenon also verifies the factor appeared in the error estimate.Last but not least, for fixed final time t f = ∆ t c = / ∆ t gt = / n , the number of wave-packets, we choose n + Example 2, 2d Case
Besides testing the properties of GWPT+HWP in 2d, we will explore the possibility of applying themethod to multi-dimensional cases. Note that we do not have rigorous numerical analysis for multi-dimensional cases. In this example, we choose the potential function as V ( x ) = − cos ( x ) − cos ( y ) .8 of 31 n K = K n ∞ ( K n ∞ ) K = K n ( K n ) Table 2. For ε = / ∆ t c = / ∆ t gt = / t f =
2, we compare the L error of K = K n ∞ with the errorof K n . The tensor product of two 50-point Gauss-Hermite rule is used to evaluate the Galerkin matrix. Reference solution for ψ :SP4 on ( − π , π ) with ∆ t REF = ε /
32 and spatial grid size ∆ x REF = πε / The initial value is chosen to be: ψ d ( x , y , ) = (cid:18) πε (cid:19) / exp (cid:20) − ( x − π / ) + ( y − π / ) ε (cid:21) . The reference solution is computed by the SP4 method with uniform spatial collocation point on [ − π , π ) and with spatial grid size ∆ x REF = πε /
16 and time step ∆ t REF = ε / K = K n ∞ and K = K n , as defined in (2.32)-(2.33), is negligible. We choose ε = / ∆ t c = / ∆ t gt = / K = K n ∞ and K = K n is quite small.Second, we will numerically verify the spectral convergence w.r.t. n in K T . We choose ∆ t c = / ∆ t gt = / t f = . L error. Wechoose K = K and ∆ t c = / ∆ t gt = / L error at t f = .
125 for variousnumber of quadrature points N Q . The quadrature rule is chosen to be the tensor product of two p N Q -point Gauss-Hermite rules to evaluate the Galerkin matrix. The result is shown in the right panel ofFigure 5.Last but mot least, we will also compare CPU time for different values of n in K n . We fixed thefinal time t f = ∆ t c = / ∆ t gt = / n + ε .
5. Conclusions
In this article, we propose a novel spectral method based on the GWPT and Hagedorn’s wave-packets. Itcircumvents the problem of artificial boundary condition of the original GWPT method, and facilitatesnumerical analysis. We provide a general framework for the error analysis of GWPT based numericalmethods, and prove that our proposed method has spectral convergence with respect to the number ofthe wave-packets in one dimension. Various numerical tests are presented to show different propertiesof the novel spectral method. The numerical results also validates our numerical analysis. This methoddoes not use fast Fourier transform, which shows the possibility of parallel algorithm in moderate high
NOVEL SPECTRAL METHOD BASED ON GWPT
29 of 31 n in K n1 -14 -12 -10 -8 -6 -4 -2 square root of N Q -12 -10 -8 -6 -4 -2 F IG . 5. Left panel: For fixed time step size ∆ t c = / ∆ t gt = / t f = . L error ofGWPT+HWP in 2d with different values of n in K = K n . The tensor product of two 50-point Gauss-Hermite rule is used toevaluate the Galerkin matrix. Right panel: For fixed time step size ∆ t c = / ∆ t gt = / t f = . L error of GWPT+HWP in 2d. The wave-packets is K = K in each test. Reference solution for ψ : SP4 on ( − π , π ) with ∆ t REF = ε /
32 and spatial grid size ∆ x REF = πε / n CPU time ε Table 3. For different values of ε , we compare the CPU time for different numbers of wave-packets n with fixed time step size ∆ t c = / ∆ t gt = / t f =
4. The tensor product of two n + Acknowledgements
Z. Zhou is supported by the National Key R&D Program of China, Project Number 2020YFA0712000and NSFC grant No. 11801016, No. 12031013. Z. Zhou is also partially supported by Beijing Academyof Artificial Intelligence (BAAI). G. Russo acknowledges partial support from ITN-ETN Horizon 2020Project ModCompShock, “Modeling and Computation of Shocks and Interfaces”, Project Reference642768, and from the Italian Ministry of University and Research (MIUR), PRIN Project 2017 (No.2017KKJP4X entitled “Innovative numerical methods for evolutionary partial differential equations andapplications”. R
EFERENCES B AO , W., J IN , S. & M ARKOWICH , P. A. (2002) On time-splitting spectral approximations for the Schr¨odingerequation in the semiclassical regime.
Journal of Computational Physics , , 487–524.B REZIS , H. (2010)
Functional Analysis, Sobolev Spaces and Partial Differential Equations . Springer New York.C
HIN , S. A. & C
HEN , C. R. (2001) Fourth order gradient symplectic integrator methods for solving the time-dependent Schr¨odinger equation.
The Journal of Chemical Physics , , 7338–7341.F AOU , E., G
RADINARU , V. & L
UBICH , C. (2009) Computing semiclassical quantum dynamics with Hagedornwavepackets.
SIAM Journal on Scientific Computing , , 3027–3041.G OLSE , F., J IN , S. & P AUL , T. (2020) On the convergence of time splitting methods for quantum dynamics in thesemiclassical regime.
Foundations of Computational Mathematics .G RADINARU , V. (2008) Strang splitting for the time-dependent Schr¨odinger equation on sparse grids.
SIAMJournal on Numerical Analysis , , 103–123.H AGEDORN , G. A. (1980) Semiclassical quantum mechanics.
Communications in Mathematical Physics , ,77–93.H AGEDORN , G. A. (1998) Raising and lowering operators for semiclassical wave packets.
Annals of Physics , ,77–104.H ELLER , E. (2018)
The Semiclassical Way to Dynamics and Spectroscopy . Princeton University Press.J IN , S., W U , H. & Y ANG , X. (2008) Gaussian beam methods for the Schr¨odinger equation in the semi-classicalregime: Lagrangian and Eulerian formulations.
Communications in Mathematical Sciences , , 995–1020.J IN , S., M ARKOWICH , P. & S
PARBER , C. (2011a) Mathematical and computational methods for semiclassicalSchr¨odinger equations.
Acta Numerica , , 121–209.J IN , S., W U , H. & Y ANG , X. (2011b) Semi-Eulerian and high order Gaussian beam methods for the Schr¨odingerequation in the semiclassical regime.
Communications in Computational Physics , , 668–687.J IN , S., , S PARBER , C. & Z
HOU , Z. (2017) On the classical limit of a time-dependent self-consistent field system:Analysis and computation.
Kinetic & Related Models , , 263–298.J IN , S. & Z HOU , Z. (2013) A semi-Lagrangian time splitting method for the Schr¨odinger equation with vectorpotentials.
Communications in Information and Systems , , 247–289.L ASSER , C. & L
UBICH , C. (2020) Computing quantum dynamics in the semiclassical regime.
Acta Numerica , , 229–401.L IU , H., R UNBORG , O. & T
ANUSHEV , N. M. (2012) Error estimates for Gaussian beam superpositions.
Mathe-matics of Computation , , 919–952. NOVEL SPECTRAL METHOD BASED ON GWPT
31 of 31 M A , Z., Z HANG , Y. & Z
HOU , Z. (2017) An improved semi-Lagrangian time splitting spectral method for thesemi-classical Schr¨odinger equation with vector potentials using NUFFT.
Applied Numerical Mathematics , , 144–159.P ASCIAK , J. E. (1980) Spectral and pseudo spectral methods for advection equations.
Mathematics of Computa-tion , , 1081.R USSO , G. & S
MEREKA , P. (2013) The gaussian wave packet transform: Efficient computation of the semi-classical limit of the schr¨odinger equation. part 1 – formulation and the one dimensional case.
Journal ofComputational Physics , , 192–209.R USSO , G. & S
MEREKA , P. (2014) The Gaussian wave packet transform: Efficient computation of the semi-classical limit of the Schr¨odinger equation. part 2. multidimensional case.
Journal of Computational Physics , , 1022–1038.S HEN , J., T
ANG , T. & W
ANG , L.-L. (2011)
Spectral Methods . Springer Berlin Heidelberg.Z
HOU , Z. (2014) Numerical approximation of the Schr¨odinger equation with the electromagnetic field by theHagedorn wave packets.
Journal of Computational Physics , , 386–407.Z HOU , Z. & R
USSO , G. (2019) The Gaussian wave packet transform for the semi-classical Schr¨odinger equationwith vector potentials.
Communications in Computational Physics , , 469–505.Z WORSKI , M. (2012)