Frozen Gaussian Sampling: A Mesh-free Monte Carlo Method For Approximating Semiclassical Schr\"odinger Equations

In this paper, we develop a Monte Carlo algorithm named the Frozen Gaussian Sampling (FGS) to solve the semiclassical Schr\"odinger equation based on the frozen Gaussian approximation. Due to the highly oscillatory structure of the wave function, traditional mesh-based algorithms suffer from"the curse of dimensionality", which gives rise to more severe computational burden when the semiclassical parameter \(\ep\) is small. The Frozen Gaussian sampling outperforms the existing algorithms in that it is mesh-free in computing the physical observables and is suitable for high dimensional problems. In this work, we provide detailed procedures to implement the FGS for both Gaussian and WKB initial data cases, where the sampling strategies on the phase space balance the need of variance reduction and sampling convenience. Moreover, we rigorously prove that, to reach a certain accuracy, the number of samples needed for the FGS is independent of the scaling parameter \(\ep\). Furthermore, the complexity of the FGS algorithm is of a sublinear scaling with respect to the microscopic degrees of freedom and, in particular, is insensitive to the dimension number. The performance of the FGS is validated through several typical numerical experiments, including simulating scattering by the barrier potential, formation of the caustics and computing the high-dimensional physical observables without mesh.


Introduction
The Schrödinger equation, which depicts the evolution of the wave function of a quantum system, is given as follows: where m o denotes the mass of the particles and denotes the reduced Planck constant. By adopting the nondimensionalisation techniques (see Appendix A), one obtains the following semiclassical Schrödinger equation in the dimensionless form: where the parameter ε, named the semiclassical parameter, describes the microscopic macroscopic scale ratio and is small within the semiclassical regime (i.e. 0 < ε 1). The quantum mechanical wave function u(t, x), valued in C, often serves as an auxiliary quantity to compute the physical quantities of interest. The physical quantities, including position, momentum and energy, etc., are referred to as "observables", which are of the form of quantum averages with respect to the wave functions of certain selfadjoint operators: Here,Ô denotes the operator related to a physical observable. For examples,Ô q = x and O p = −iε∇ are the operators corresponding to the position and momentum observables, respectively. In the semiclassical regime, namely ε 1, the solution to equation (1.2) is highly oscillatory both in space and time on the scale O(ε) and doesn't converge in strong sense as ε → 0, which gives rise to significant computational burdens [13,17]. Nevertheless, various effective numerical algorithms have been proposed so far, succeeding in capturing the highly oscillatory wave functions accurately with affordable computational complexity. The existing algorithms can be divided into two classes: direct discretization methods and asymptotic methods. Among the direct discretization methods, one of the most prevailing choices is the time splitting spectral method proposed in [1], where the meshing strategy ∆t = O(ε) and ∆x = O(ε) is practically sufficient for accurate approximation of the wave function. The time step size of the time splitting spectral method can be even relaxed to O(1), in other words, independently of ε, if one is only interested in capturing the correct physical observables, whereas one still needs to resolve the spatial oscillations. The asymptotic methods, including the WKB methods [5,6] and various wave packets based methods [29,14,18,19,20,27,12,24], provide accurate numerical approximation to the wave function with a relatively less computational cost with respect to the parameter ε compared to the direct discretization methods. The wave packets based methods have been extended to more complicated quantum dynamics, such as the vector potential cases and systems with random inputs [30,27,12,24], providing broad application prospects in many applications.
Although the existing algorithms have made great progress in efficiently simulating the wave functions and the physical observables in the semiclassical regime, they still suffer from the "the curse of dimensionality" when considering high dimensional cases, since those methods necessarily reply on an Eulerian mesh either in the initial decomposition or in the time evolution. When the dimension number m 1, the computation and memory amounts of the both classes of methods increase exponentially as m increases, which often leads to more infeasible computational cost when ε 1. For example, for the time splitting spectral method, the degree of freedom to approximate the wave function on an m−dimensional mesh is O(ε −m ), not to mention the increased computational complexity to evolve the wave function on this high dimensional mesh. Therefore, a novel numerical method for the Schrödinger equation in semiclassical regime, whose cost is insensitive to ε and the dimension number m, is highly desirable.
To the best of our knowledge, there has been very few existing methods, which are feasible or have the potential to be suitable for high dimensional simulations. From this perspective, the path integral sampling method proposed recently in [22,21] can be viewed as a pioneering work, which is constructed on the so-called Frozen Gaussian approximation (FGA). The FGA represents the wave function with a superposition of Gaussian functions in the phase space with fixed widths and is often wrote in the form of an integral on the phase space: A(t, q, p)e iΘ(t,x,q,p) ε dqdp, (1.4) where the variables A and Θ are can be approximately obtained through evolving a system of Lagrangian ODEs, thus the full quantum dynamics is reduced to an ensemble of wave packets dynamics. There have been quite a few works related to the FGA, including the theoretical studies in the chemistry community [8,9,15,16], the rigorous convergence analysis [28] and the FGA based algorithms for solving hyperbolic systems which reply on the direct phase space discretization [18,19,20]. However, it is not until [22,21] that people find another distinct prospective, giving the FGA ansatz (1.4) a probabilistic interpretation, and show that the Monte-Carlo sampling method based on the FGA gives a feasible algorithm which is mesh-free in computing the physical observables. The main spirit of [22,21] can be summarized as follows: firstly, as an integral on the phase space, the FGA ansatz can be viewed as an expectation of random variables on the phase space, which can be calculated by solving the parameter ODEs; in addition, when the initial data is Gaussian, the module function |A(0, q, p)| is also Gaussian and naturally serves as an ideal choice of a sampling density function on the phase space. However, what remains to be done is that there have been no rigorous numerical analysis on the sampling error nor discussion on how to extend the sampling method to general cases when the initial data is not Gaussian. Our current work aims to fill the vacancy on these issues.
In this work, we develop a Monte Carlo method which relies on the probabilistic interpretation on the FGA ansatz, and we shall name it the Frozen Gaussian Sampling (FGS). The essence of the FGS is to properly choose a reference density function on the phase space, and reformulate the FGA representation (1.4) as an expectation of a complex-valued functional with respect to the Frozen Gaussian trajectories. In this view, the FGS is a path-integral approximation of the wave function, and in the dynamical picture, the evolution of the quantum state is replaced by propagating an ensemble of Frozen Gaussian wave packets, where the initial phase space coordinates are sampled from the reference density function. The design of an efficient FGS algorithm boils down to proposing a suitable reference sampling measure, which is closely related the form of the initial condition of the wave function, and the quality of an FGS algorithm is measured by the variance of the functional to be sampled. Hence, the FGS provides a subtle shift in the computation paradigm and the analysis framework, which, as we shall elaborate, turns to show great promises as a simulation tool for high dimensional problems.
We focus on two of the most important types of initial quantum states: the Gaussian wave packet and the WKB data, and study the respective strategies in designing the density functions that balance the need of variance reduction and sampling convenience. For the Gaussian initial condition scenario, we rigorously prove that, with the feasible choice of the sampling density, the variance of the functional to be sampled is independent of the scaling parameter ε . Furthermore, the number of samples to reach a certain sampling error δ is in the order O c m δ −2 , where c is a O(1) constant and m is the dimension number. In particular, if one aims to halve the sampling error, only four times of the samples are needed, which is independent of the dimension m or the semiclassical parameter ε. Different from all the existing methods, this result shows that the number of wave packets does not to "resolve" any spatial oscillations, and is of a sublinear scaling of the microscopic degrees of freedom. Besides, the computation cost grows rather gently as m increases, and the whole FGS algorithm is parallelizable, which makes high dimensional simulations rather affordable. For WKB initial data cases, we propose a sequential sampling strategy leveraging the mono-kinetic structure of the highly oscillatory wave function. We numerically verify that the sample size to reach a certain error is also insensitive to ε, which still outperforms the existing mesh-based algorithms whose computational cost grows rapidly as ε → 0.
It is worthy emphasizing that, we only require that the initial condition to be in cer-tain forms to carry out efficient initial sampling, but we do not require that the solutions remain in such ansatz forms. This is a natural consequence of the construction of the FGS algorithm, while most asymptotic methods break down when the solution ansatz is no longer valid. And we have carried a set of 2D numerical tests to illustrate such a property, where the initial Gaussian wave packet is scattered into various profiles by different potential barriers, and the FGS method is able to capture the correct evolution behavior with a reasonably small sample size.
Another main highlight of the FGS lies in that it is mesh-free in computing the physical observables. By taking advantage of path integral representation, the FGS succeeds in computing a wide class of physical observables without reconstructing the wave functions. Therefore, if one is only interested in physical observables rather than the detained behavior of the wave function, the FGS serves as an efficient algorithm to simulate high dimensional Schrödinger equation, since it is a mesh-free algorithm of a sublinear scaling complexity and is naturally compatible with parallel computing. To demonstrate the performance of the FGS in computing observables, we implement the FGS to compute the physical observables in high dimensional quantum systems with a harmonic potential and it is shown numerically that the computational cost to reach a certain error is rather affordable in various dimensions. With the sample size quadrupling when the dimension number increases by one, we observe that the errors in the observables is uniformly bounded in all dimensions tested.
The rest of this paper is organized as follows: in Section 2, we review the basic setup the FGA and give the detailed procedures to compute the wave functions and physical observables with the FGS. We also present the error analysis framework of the FGS and elaborate the difficulties to be overcome when choosing the density function for the FGS. In Section 3, we focus on Equation (1.2) with Gaussian initial data and prove that the quantitative estimates of the FGS sampling error . In Section 4, we propose a feasible sampling strategy to apply the FGS on Equation (1.2) with WKB initial data. In Section 5, the numerical performance of the FGS is validated through five carefully designed numerical experiments, including simulating the Gaussian barrier potential and the caustics and computing the high-dimensional position and momentum observables without mesh.
we present the frozen Gaussian sampling (FGS) algorithm as a general computation methodology for semiclassical quantum dynamics and discuss the practical implementation of computing physical observables based on the FGS algorithm.

The frozen Gaussian representation for the semiclassical Schrödinger equation
The frozen Gaussian approximation is an asymptotic approximation to the solution of Equation (1.2) with O(ε) error [28]. It is based on an integral representation on the phase space: Here the phase function Θ is given by For a given initial condition q and p, the evolution of Q and P is governed by the classical Hamitonian H(q, p) = E(q) + 1 2 |p| 2 and defines a trajectory on phase space R 2m called an FGA tracjectory, which also determines the value of the amplitude A and action S. We name the variables Q, P, S, A as FGA variables. The evolution of the FGA variables is determined by the ODE system as follows: where we use the short hand notations Moreover, the initial condition is given by In addition, we remark that the phase space integral (2.1) can also be viewed as an integral of the FGA trajectory. Denote the FGA tracjectory by z t = (Q(t), P (t)), where we hide the dependence of the FGA variables on (q, p) for simplicity (we use the similar short hand notation in the rest of the paper). Then z t satisfies the Hamiltonian flow and defines a one-to-one map z 0 → z t on the phase space. In other words, the FGA trajectory z t is specified once the initial condition (q, p) is fixed and the FGA variables S and A are also specified with a given trajectory z t . The idea of viewing the FGA ansatz (2.1) through FGA trajectory is important in the introduction of the frozen Gaussian sampling.
To conclude this subsection, we directly quote two theorems from [18,28]. The first theorem shows the FGA ansatz integration (2.1) reproduces the initial condition: where the initial amplitude A(0, q, p) given in Equation (2.10) is determined by u in .
The second theorem shows that the L 2 error of u FGA with respect to the exact solution is O(ε): Theorem 2.2. Consider the semiclassical Schrödinger equation (1.2) with the exact solution u(t, x). Assume that the potential function E(x) ∈ C ∞ (R m ). For any t > 0, we have the following L 2 error estimate result for the FGA ansatz (2.1): where C A (t, m) is a parameter independent of the scaling parameter ε.

The frozen Gaussian sampling algorithm
We have seen in Section 2.1 that the FGA ansatz (2.1) is an integral on the phase space. With the FGA, the task of solving a PDE is converted into solving an ensemble of ODEs (2.3)-(2.6). However, difficulties still remain in computing the integral (2.1), which can be of high dimensions. Direct mesh method of numerical integration may suffer from "the curse of dimensionality" and does not lead to a practical algorithm. To deal with the high-dimensional integral (2.1), we seek a stochastic method for calculation. To this end, we write the FGA ansatz in the form of an expectation and approximate the integral (2.1) numerically with a Monte Carlo method. The resulting stochastic algorithm is called the frozen Gaussian sampling (FGS), which also gives a meshfree implement to compute physical observables. In this subsection, we elaborate the FGS algorithm and discuss its advantages and difficulties to be further investigated.
To write the FGA ansatz (2.1) into an expectation, we introduce a real-valued probability density function π(·) on the phase space R 2m . The density π then defines a probability measure P 0 on R 2m : where P 0 (R 2m ) = R 2m π(z 0 )dz 0 = 1. With the probability measure, we can rewrite the FGA ansatz (2.1) as (2.14) Thus, we may use a Monte Carlo sampling for u FGA (t, x): are M independent identically distributed samples with the probability density function π. We refer to u FGS as the frozen Gaussian sampling (FGS) wave function. We also note that Λ (t, x, z 0 ) can be treated as a random variable related to z 0 with the expectation Algorithmically, once we have sampled z 0 , the FGA trajectory and FGA variables are available through evolving the ODE system Equation (2.3)-(2.6) up to time t. With enough samples of the FGA trajectories, we approximates the FGS wave function u FGS numerically by Equation (2.15). From the discussion above, we give the pseudo-code of the frozen Gaussian sampling (FGS) in Algorithm 1. For simplicity, in the rest of this paper, we use the short hand notations to the FGA trajectories and FGA variables: We note that analogous short hand notation may be also applied to other variables without confusion. However, there are certain difficulties when actually executing the FGS algorithm, especially in high dimensional cases. First, to compute the initial amplitude, one needs to calculate the integral (2.10). For a Gaussian wave packet initial condition u in , the initial amplitude A(0, q, p) can be calculated explicitly in Equation (3.2), which has been pointed out in [22]. However, for a more complex initial condition u in such as a WKB initial condition, the explicit expression of A(0, q, p) is unavailable. To compute the initial amplitude (2.10), the computational cost of numerical integration is huge due to the high dimension and high frequency oscillation structure of the integrand.
Next, let us discuss the error analysis framework of the FGS wave function u FGS . Obviously, the L 2 error function E 0 = u FGS − u L 2 of the FGS wave function with respect to the exact solution can be divided it into two parts: Here E FGA = u FGA − u L 2 denotes the asymptotic error of the FGA ansatz and E S = u FGS − u FGA L 2 denotes the Monte Carlo sampling error. Thanks to Throrem 2.2, we arrive at where C A (t, m) is a constant depending on time t and dimension m and independent of the rescaled parameter ε. However, with unreasonable choice of density function π, the Monte Carlo sampling error E S may increase dramatically as the dimension number m increases, thus formidably many samples are needed. Therefore, to apply the FGS in high dimensional cases, variance reduction techniques such as importance sampling are necessary to reduce the number of samples. Ideally speaking, we aim to choose a reference density π which not only leads to a minimized sampling variance, but also easy to sample numerically. If we solely consider the variance reduction prospective, an obvious choice of sampling density can be given as follows: For a Gaussian wave packet initial condition, the density function given by Equation (2.20) actually defines a multivariate normal distribution and the resulting FGS algorithm is proven numerically practical and effective [22,21]. For a more complex WKB initial condition when the explicit expression of |A(0, z 0 )| in Equation (2.20) is unavailable, how to choose π(z 0 ) is still unclear. As far as we know, no rigorous numerical analysis on sampling error E S has been carried out prior to this work. The aim of the rest of this paper is to investigate and overcome the above difficulties of the FGS. In Section 3, we prove rigorously that the sampling error E S of FGS with sampling density (2.20) is independent of ε with a Gaussian initial condition. In Section 4, we introduce an asymptotic expression approximating the initial amplitude (2.10) when the initial data is of a WKB form and present a new initial sampling method of the FGS algorithm based on this approximation. In Section 5, we show various numerical examples to test the practical performance of the FGS with these new techniques, especially in high dimensional cases.

Computing observables without mesh
For many applications, the goal is not to approximate the wave function in semiclassical Schrödinger equation (1.2) itself, but rather computing certain observables. To calculate an observable, reconstructing the wave function on a mesh is often neither feasible nor necessary, especially for high dimensional cases. In this subsection, we introduce a mesh-free method to compute a wide class of observables, which only use the FGA trajectories sampled by the FGS and doesn't require the reconstruction of the FGS wave function. The method exhibits great advantages on computing high dimensional problems.
Let O = u|Ô|u whereÔ denotes a self-adjoint operator asscoiated with a physical observable. Then we compute O based on the FGS wave function u FGS as follows: We observe that the only dependence on x in Equation (2.21) is in the term g j,k , which is an integration over the Gaussian wave packets. For a wide class of observablesÔ, the term g j,k (t, x) can be carried out analytically through Gaussian integration formula, hence the observables O = u|Ô|u can be computed without mesh through computing the double summation in Equation (2.21). This class of observables include functions of the position variable and the momentum variable, which is in a polynomial form or can be approximated by a polynomial. Now we take the position observable and the momentum observable as examples to show the analytical expression of g j,k . For the position observableÔ = x, we obtain For the momentum observableÔ = −iε∇ x , we obtain To sum up, given the FGA trajectories and FGA variables, we can reconstruct the observables in a mesh-free fashion without computing the FGS wave function. Clearly, the complexity for the physical observable evaluation (2.21) is O(M 2 ). In the following sections, we shall further explore the dependence of the ensemble size M on the scaling parameter ε and the dimension number m and show great potential of the FGS in computing physical observables in high-dimensional cases.

Frozen Gaussian sampling with Gaussian initial conditions
In Section 3, we consider a m−dimension normalized semiclassical Gaussian wave packet centering atq with momentump: The semiclassical Gaussian wave packet is one of the simplest yet most important wave function in quantum dynamics, whose semiclassical limit arrives at a delta function on the phase space. With the Gaussian integration formula, we are able to calculate the initial amplitude A(0, q, p) (see Equation (2.10)) explicitly and the resulting density function π as shown in Equation (2.20) is Gaussian. Finally, we rigorously prove the sampling error E S of the FGS with density function is independent of ε and discuss its dependence on the dimension number m.

Initial sampling and initial amplitude for Gaussian initial condition
Plug the Gaussian initial condition (3.1) into Equation (2.10), we calculate A(0, q, p) analytically with the Gaussian integration formula: and As we have stated, the probability density function π(z 0 ) given by Equation (2.20) is a feasible choice of initial sampling. For Gaussian initial condition (3.1), the resulting probability distribution for the initial position and momentum z 0 = (q, p) is a multivariate normal distribution centered at q =q and p =p with probability density function: Algorithmically, we can directly sample z 0 = (q, p) from Gaussian density function in Equation (3.4), then the initial amplitude A(0, q, p) can be computed analytically.

Main result of sampling error
In this subsection, we aim to establish a rigorous error estimate for the FGS algorithm. We justify that the Monte Carlo sampling error E S = u FGS − u FGA L 2 is independent of the scaling parameter ε and study the relationship between E S and the dimension number m. With the error estimate, we are able to discuss the proper sample size for the FGS for given ε and m and illustrates the advantages of the FGS algorithm, especially in high-dimensional cases.
Recall that u FGA denotes the FGA ansatz defined in Equation (2.1) and u FGS denotes the FGS wave function given in Equation (2.15). Now that u FGS is determined by a series of random variables z (j) 0 on the phase space, we define the quadratic expectation of sampling error E S and total error E 0 as follows: Then plug Equation (2.15) into the expression of E S : where is regarded as the sampling variance of the FGS wave function. The derivation of Equation (3.7) uses the fact that z (j) 0 are independent identically distributed with density π. We aim to estimate an upper bound for the integral I FGA .
Before we present the main result of error estimate, we introduce the definition of the subquadratic condition. We say the potential function E(x) satisfies the subquadratic condition if there exists a constant C E such that The subquadratic condition is commonly used in analysis of hyperbolic systems [28,19,20,21] and we also assume the subquadratic condition in our theorem. With all of the preparation, let us show the main result of this section: 2) up to time t with the normalized Gaussian initial condition (3.1). Suppose that the potential function E(x) ∈ C ∞ (R m ) and satisfies the subquadratic condition (3.9). When applying the frozen Gaussian sampling (FGS) as in Equation (2.15) with the initial sampling density π as in Equation (2.20), there exists a constant C 0 (m, t) independent of the scaling parameter ε and the initial data u in such that the sampling variance given in Equation Moreover, if the assumptions in Lemma 3.6 are valid, the constant C 0 (t, m) in Equation where γ 0 (t) > 1 is a constant independent of ε and m.
We leave the proof of Theorem 3.1 to Section 3.3 and discuss on the theorem first. As a corollary, the following estimate on the quadratic expectation E S and E 0 defined in Equation (3.12) and (3.6) is obvious: Under the same assumptions in Theorem 3.1, we have the following estimates to E S and E 0 : and To conclude this subsection, let us demonstrate the advantages of the FGS over the time-splitting spectral methods (see [1]) through the error estimate in Corollary 3.2. For simplicity, we assume a j = 1 for j = 1, 2, · · · , m and ignore the asymptotic error (3.14) where is another constant independent of ε and m. In order to guarantee that E 0 ≤ δ, the sample size of the FGS should be more than In contrast, if we use a mesh-based algorithm such as time-splitting spectral methods in [1] to numerically approximate the semiclassical Schrödinger equation (1.2), the spatial step length of the mesh should be at least of order O(ε) to capture the highly oscillatory wave function. For an m−dimensional semiclassical Schrödinger equation, the spatial degree of freedom of the mesh-based method should be at least Recall that the sample size N FGA is the m−th power of an O(1) parameter. Therefore, when scaling parameter ε 1 or the dimension number m is large, the computational cost to solve a given semiclassical Schrödinger equation with the FGS is much smaller than those with the grid based methods, i.e. N FGA N TSP .
Remark 3.3. We remark that though the sampling size N FGS increases exponentially as m increases, we examine from numerical examples that the increasing rate of N FGS is much smaller compared to ε −1 . For more information of this issue, see Example 5 in Section 5.3.

Proof of Theorem 3.1
Before we prove Theorem 3.1, let us first collect some properties of the Hamiltonian flow on the phase space. In this subsection, we assume the potential function E(x) ∈ C ∞ (R m ) and satisfies the subquadratic condition (3.9) without special instructions.
For simplicity, we denote the map on the phase space from initial time 0 to time t as: t, q, p), P (t, q, p)) =: (Q κt (q, p), P κt (q, p)) . For the transformation κ on the phase space, we denote its Jacobian matrix as Note that the transformation κ is canonical, i.e. its Jacobian matrix J κt satisfies for any (q, p) ∈ R 2m . Here I m denotes m × m identity matrix.
The following propositions give uniform bounds of the Jacobian matrix J κt (q, p) and its derivatives. where | · | denotes the matrix norm.
Note that Proposition 3.4 was proven in Proposition 1 of [28]. For the canonical transformation κ, we define where ∂ z = ∂ q − i∂ p . The invertibility and boundness estimate of Z κt can be proven by mimicking the proof of Lemma 5.1 in [20]. We summarize them into the following proposition.
Proposition 3.5. For given time t > 0 and k ∈ N, there exists a constant K 0 (t, k) such that sup
With the properties of the matrix J κt and (Z κt ) −1 , we finally establish the estimate of the FGA variable A(t, q, p): Lemma 3.6. Consider the FGA variable A(t, q, p) that solves Equation (2.6) on the phase space. For given time t > 0, there exists a constant K 1 (t, m) such that Moreover, if we assume the constant K 0 (t, 1) in Proposition 3.4 and 3.5 is independent of the dimension number m, there exists a parameter γ * (t) > 1 such that The proof is straightforward through Gronwall's inequality.
Proof. Based on Equation (2.6), we arrive at Note that the constants in Proposition 3.4 and 3.5 are independent of the dimension number m, we use the constant K * (t) = K 0 (t, 1) to estimate Equation (3.25): Then apply Gronwall's inequality to Equation (3.25): where γ * (t) = e 2(K * (t)) 2 is independent of ε and m.
Remark 3.7. We remark that the importance of the three propositions and lemmas lie in the fact that the estimate to the FGA trajectories and FGA variables is uniform with respect to the initial phase z 0 = (q, p). Hence we can make global estimate on the integration I FGA . Validity of the propositions and lemmas rely on the subquadratic condition of the potential function E(x).
After the preparation above, we finally move to the proof of Theorem 3.1: Proof. According to Equation (3.8), we have (3.28) and it suffices that we focus on the following secondary moment: Now recall that the sampling density π is given by Equation (2.20), then we can simplify Equation (3.30): Next we integrate M (2) (x) in Equation (3.31) with respect to the position variable x: Now we apply Lemma 3.6 in Equation (3.32): (3.33) Here K 1 (t, m) is a constant independent to ε. According to Equation (3.3), we can calculate the integration of |A(0, z 0 )| on the phase space: (3.34) Finally we use Equation (3.34) in Equation (3.33): Hence the right hand side of Equation (3.35) is a constant independent to ε. Let C 0 (t, m) = 2 3m 2 (K 1 (t, m)) 2 , we finally arrive at Equation (3.10). Through the expression of C 0 (t, m), it is obvious that Equation (3.11) is valid provided that Equation (3.24) in Lemma 3.6 is satisfied.

Frozen Gaussian sampling with WKB initial conditions
In this section, we study the semiclassical Schrödinger equation (1.2) with the following normalized WKB initial condition: where a in and S in are two real-valued functions on R m , which are both independent of ε. The initial amplitude a in follows With the celebrated WKB approximation (see Appendix B for a brief introduction), the exact solution to Equation (1.2) has an asymptotic description of WKB type before caustics onset, provided that a in is smooth and decays rapidly [13]. Practically, we mainly consider the cases with the following Gaussian initial amplitude: where a j > 0 are sufficiently big andx j ∈ R (j = 1, 2, ·, m) are given parameters. We denote the mean of a in asx = (x 1 , x 2 , · · · , x n ) T .

(4.4)
The cases with general initial amplitude a in may be treated in a similar fashion and we refer to Remark 4.4 for discussion on this issue. Compared to Gaussian wave packet initial data cases where the explicit expression of A(0, q, p) is available, the difficulty of applying the FGS in WKB initial data cases lies in that there is no natural choice for the reference density function π(z 0 ) which achieves the balance between variance reduction and convenience of sampling. Moreover, the WKB approximation to Equation (1.2) may break down after caustics onset, which brings additional difficulties to the numerical approximation of Equation (1.2). In this section, we put forward a new initial sampling algorithm based on the stationary phase approximation to the initial amplitude (2.10). Then we validate the numerical performance of the FGS associating with the new initial sampling strategy through various typical examples in Section 5.

Description of the algorithm
As we have stated, one of the difficulties of the FGS algorithm with a WKB initial condition is that A(0, q, p) (given in Equation (2.10)) can not be calculated explicitly, hence there is no obvious way to do initial sampling. In this section, we apply the stationary phase method to approximate the initial amplitude A(0, q, p) in Equation (2.10) and give a practical sampling method based on the analytical approximation to A(0, q, p). Before we start, let us summarize the stationary phase method into the following lemma: Lemma 4.1. Let α(x) and Ψ(x) be two smooth functions on R m . Suppose that α(x) ∈ C ∞ 0 (R m ) and Ψ(x) ∈ C ∞ (R m ). Assume that ∇Ψ(x) has a finite number of zero points {y k } N k=1 on the support of α and det ∇ 2 Ψ(y k ) = 0 for each zero point y k . We consider the highly-oscillatory integration As ε → 0, I ε has the following asymptotic behavior: where the term sgn ∇ 2 Ψ (y k ) denotes the number of positive eigenvalues of the matrix ∇ 2 Ψ (y k ) minus the number of negative eigenvalues.
A proof to Lemma 4.1 can be found in [3].
According to Equation (4.6), the initial amplitude A(0, q, p) in Equation (2.10) can be calculated asymptotically as follows: In order to illustrate the spirit of sampling strategy, we make some reasonable assumptions which simplify the asymptotic description of A(0, q, p) in Equation (4.7). The cases when the assumptions break down will be discussed in Remark 4.3. We will show the Assumption 1. For any p ∈ R m , there exists at most one y ∈ R m such that ∇S in (y) = p, i.e. ∇S in is an injective mapping. The inverse mapping of ∇S in is given by T : p → y. The definition domain of T may not be R n and we denote it by D(T ), i.e.
T : D(T ) → R m , p → y. (4.8) Recall that we assume the amplitude of the WKB initial data a in (x) decays exponentially as x → ∞, the support of a in can be viewed as a small interval. Hence it is reasonable to assume that for any p ∈ R m , there exists at most one y ∈ R m such that ∇S in (y) = p and a in (y) = 0. With this assumption, the summation in Equation (4.7) can be reduced to only one term. Besides, we will show later in Remark 4.3 that the difficulties caused by multivalueness of the mapping T are not essential.
For simplicity, we use the variable y to denote the unique image of the mapping T defined in Equation (4.8) at p ∈ D(T ) in the rest of this section without special instructions.
The next assumption is a technique assumption without lose of generality: Assumption 2. The inverse mapping T defined in Equation (4.8) is smooth and satisfies det ∇ 2 S in (y) = 0, ∀p ∈ D(T ). (4.9) Based on the two assumptions, the asymptotic approximation (4.7) of the initial amplitude A(0, q, p) can be simplified as follows: where χ D(T ) (p) stands for the indicative function on D(T ). Then the modulus of Equation (4.10) is given by For Gaussian initial condition cases, we directly use the real-valued Gaussian function |A(0, q, p)| on the phase space as the sampling density π in Equation (2.14), which results in the minimum sampling variance. However, for the WKB initial condition cases, it is impossible to sample towards |A sp (q, p)| given by Equation (4.11) directly since the non-Gaussian function |A sp (q, p)| is hard to sample. To balance between variance reduction and convenience of sampling, we choose the following sampling density function: where Z is a normalization parameter to ensure π a probability density, i.e. R 2m π(q, p)dqdp = 1.
We note that the only difference between |A sp (q, p)| in Equation (4.11) and the probability density π in Equation (4.12) lies in the exponent of the term det ∇ 2 S in (y) in the denominator. Hence the choice of π in Equation (4.12) is rather acceptable from the prospective of variance reduction. The following proposition shows that the probability density π in Equation (4.12) is also easy to sample, provided that the the amplitude a in is given by Equation (4.3): Proposition 4.2. Consider the WKB initial data (4.1) where the amplitude function a in of is given by Equation (4.3). Assume that Assumption 1 and Assumption 2 are valid. Suppose that the random variables Q and P are sampled with respect to the probability density π(q, p) defined in Equation (4.12), then Q and P follows 1. The marginal distribution of the random variable Q is a normal distribution with expectation E π q =x, wherex as in Equation (4.4) denotes the mean of a in .
2. Denote Y = T (P), then the conditional distribution of the random variable Y with respect to Q is a normal distribution.
3. The normalization parameter Z in Equation (4.12) is given as follows: (4.13) Proof. Recall that the mapping T as in Equation (4.8) is the inverse mapping of the gradient function ∇S in . We refer to J T (p) as the Jacobian matrix of the mapping T at p. Based on Assumption 2, the determinant of J T (p) follows: , ∀p ∈ D(T ). (4.14) Then we apply Equation (4.14) and the integral substitution rule to compute the marginal density π Q (q) of the random variable Q: Then the marginal distribution of the random variable Q is a normal distribution N (x, Σ 1 ), where the covariance matrix Σ is a diagonal matrix as follows: , · · · , ε + 1 a n . (4.16) To calculate the expression of the normalization parameter Z given in Equation (4.13), we just need to integrate π Q (q) by q to ensure that π Q (q)dq = 1, i.e.
For fixed Q = q, we denote the conditional density of the random variables P and Y with respect to Q as π P|Q (p) and π Y|Q (y). According to Equation (4.15), we compute π P|Q (p) as follows: Then we can compute π Y|Q (y): Note that the exponent term of the conditional density function π Y|Q (y) is quadratic of y, therefore the conditional distribution of Y with fixed Q is a normal distribution. For fixed Q = q, the conditional distribution of Y follows N (µ 2 , Σ 2 ), where the mean vector and covariance matrix is given by and , · · · , ε εa n + 1 . Motivated by Proposition 4.2, one can sample the random variables (q, p) obeying the probability density π defined in Equation (4.12) on the phase space by sampling two Algorithm 2 Initial sampling of (q, p) on the phase space for WKB initial data (4.1) with amplitude (4.3) 1: Initial position sampling: Sample q ∼ N (x, Σ 1 ), wherex and Σ 1 is given by Equation (4.4) and (4.16).
normally distributed random variable q and y, where p = ∇S in (y). The initial sampling method bases on Proposition 4.2 can be summarized into the following algorithm: To conclude this subsection, we make two remarks on the initial sampling algorithm 2. The first remark elaborate what happens when Assumption 1 is not satisfied. The second remark discuss the cases of WKB initial data (4.1) with general amplitude a in and generalize Proposition 4.2 to these cases. Remark 4.3. Now We discuss the possible treatments when Assumption 1 is violated. If there exists more than one y ∈ R n such that ∇S in (y) = p for some p ∈ R n , the asymptotic approximation A sp (q, p) of the initial amplitude A(0, q, p)(see Equation (4.7)) can't be reduced to only one term as in Equation (4.11), i.e.
The sampling density associating with A sp (q, p) in Equation (4.22) can be given through mimicking the density function (4.12): Note that the density function π sp determines a mixture distribution and each component of the summation in the density function π sp can be sampled through algorithm 2. Then it is natural to adopt the method of sampling the mixture model (see Section 4.2 in [4]) to sample π sp . Therefore, the difficulties resulted from the invalidity of Assumption 1 is not essential.
Remark 4.4. We discuss the possible treatments when a in is not Gaussian. Since a in is of a general profile, Proposition 4.2 no longer holds and there is no direct way to sample the probability density π given in Equation (4.12). To overcome these problems, we may approximate a in with a Gaussian amplitude a S in and sample the following probability with Algorithm 2: where we have replaced a in in Equation (4.12) with the approximated Gaussian amplitude a S in . Ideally, the density π S balances the need of variance reduction and convenience of sampling such that it serves as the sampling density of the FGS when a in is not Gaussian.

Numerical examples
In the previous sections, we have elaborated the FGS algorithm with Gaussian initial data and WKB initial data and studied the sampling error with respect to the sampling size M , the scaled Planck parameter ε and the dimension number m. In this section, we present five typical examples using the FGS to numerically approximate the semiclassical Schrödinger equation (1.2) to illustrate the theoretical results on the sampling error and manifest the practical potential of the FGS in solving Equation (1.2), especially in high-dimensional cases with small ε. The numerical examples include exploring the relationship between the sampling error and parameter ε, applying the FGS in challenging cases such as caustics and computing the physical observable with the FGS in high-dimensional cases. The rest of this section is organized as follows. In Section 5.1 and 5.2, we investigate the cases with Gaussian initial data and WKB initial data respectively. In section 5.3, we compute the physical observables for a high-dimensional problem with the FGS and compare the numerical results with the exact value obtained through the Hagedorn wave packet method.

Example 1: 1D Gaussian initial data examples
In this example, we aim to justify the theoretical results on the sampling error E S in Section 3 through numerical examples. We implement the FGS to approximate Equation (1.2) at time t = 0.5 and study the relationship among the sampling error E S , the scaled Planck parameter ε and ensemble size M .
In the Gaussian initial condition (3.1), let m = 1, a 1 = 2,p = − 1 2 andq = 1 2 . Then we obtain the following Gaussian wave packet initial data: Here we consider two different potential functions respectively: the harmonic potential E 1 (x) = x 2 2 and the torsion potential E 2 (x) = 1 − cos x. Recall that the L 2 sampling error E S is defined as follows: where u FGA denotes the FGA ansatz (2.1), which deviates from the exact solution to Equation (1.2) by the asymptotic error of the FGA ansatz. Note that u FGA is numerically unavailable, we choose the FGS wave function with a sufficiently large ensemble size M 0 as an accurate approximation of the FGA ansatz where M 0 = 1.5 × 10 5 is chosen large enough.
To compute the FGS wave function u FGS and u FGA , we use the fourth order Runge-Kutta method to solve the ODE system (2.3)-(2.6) in the time interval [0, 0.5] with time step length ∆t = 0.01 and reconstruct the wave function u FGS at t = 0.5. We emphasize that ∆t is chosen sufficiently small, hence the error due to the approximation of the ODE system is negligible. The wave function reconstruction is carried on the spatial interval [−π, π]. Also note that the FGS wave function u FGS with a fixed ensemble size M may vary with different initial samples z (j) 0 M j=1 , the sampling error E S is actually a random variable. Therefore, when approximating the FGS wave function u FGS , we average over numerical results obtained through 30 sets of initial samples for each simulation of E S (M, ε).
In Table 1 and Figure 1, we present the sampling error E S (M, ε) of the FGS with for different scaling Planck parameter ε = 0.0256, 0.0064, 0.0016, 0.0004 and two potential functions E 1 and E 2 respectively. One can easily notice that E S is independent of ε for both potential functions and the convergence rate of the sampling error E S with respect to the ensemble size M is of order O(M − 1 2 ), which is in accordance with the theoretical results in Theorem 3.1 and Corollary 3.2. We also underline that it only takes 400 samples for the sampling error E S to reach O(10 −2 ). As a reference, we have stated in Section 3 that the computational cost of the mesh-based algorithms to reach certain error should be at least O(ε −1 ), which indicates that the FGS is a practical and accurate algorithm to numerically approximate Equation (1.2) with 1D Gaussian initial data, especially for the cases when ε is small. Example 2: 2D Gaussian initial data example As we have stated, the FGS wave functions is constructed by an ensemble of Gaussian wave packets Λ z (j) 0 (see Equation (2.16) for definition) related to an ensemble of FGA trajectories with a certain probability distribution. In Example 1, the potential functions are so simple that the exact solutions to Equation (1.2) are still Gaussian wave packets at time t = 0.5, which lowers the difficulty of constructing the wave function from the ensemble of Gaussian wave packets Λ z (j) 0 . In this example, we implement the FGS algorithm on 2D semiclassical Schrödinger equation with more complicated potential functions, where the exact solutions to Equation (1.2) are not Gaussian wave packets as time evolves. Though it is much more difficult to approximate a wave function that is not Gaussian from an ensembles of Gaussian wave packets, the FGS method is able to capture the correct evolution behavior with a reasonably small simple size.
In the Gaussian initial condition (3.1), let m = 2, a 1 = a 2 = 2,p = (−1, −1) and q = (1, 1). The scaling parameter is fixed at ε = 1 96 . Then we obtain the following 2D Gaussian wave packet: 2 + 10e −5|x| 2 . We illustrate of the potential E 3 and E 4 in Figure 2. Intuitively, due to the effect of the exponential term, the shapes of the potential functions E 3 and E 4 show peaks at the center. The main difference between the two potentials lies in that the peak of E 4 is much higher than E 3 , which leads to different behaviors of the exact solutions to Equation (1.2).
In this example, we apply the Strang Splitting Spectral Method (SP2) [1] with sufficiently small step lengths to compute the exact solution to Equation (1.2). We reconstruct the FGS wave function on spatial region [−π, π]×[−π, π] through Algorithm 1 and compare the resulting FGS density function with the exact density function ρ(t, x) approximated by the SP2. Similar to the previous example, we use the fourth order Runge-Kutta method to solve the ODE system (2.3)-(2.6) in the time interval [0, 0.5] with time step length ∆t = 0.01 and reconstruct the FGS wave function at different times. We choose the ensemble size M = 2000. Note that M = 2000 is of lower order than O(ε −2 ), which is the degree of freedom of the mesh-based method SP2.
In Figure 3 and 4, we show the evolution of the FGS density functions with two different potential functions E 3 and E 4 and compare them with the density functions computed via the SP2 algorithm. The density functions under both potential functions are no longer Gaussian after certain time. For E 3 in Figure 3, the wave packets pass through the potential peak at the center and disperse into different directions, which makes the wave packet crescent. In contrast, the potential peak of E 4 is so high and steep that the wave packets in Figure 4 rebound against the potential peak and disperse into the opposite directions against the wave functions in Figure 3. For both cases, the profile of the FGS density functions conform with the exact solutions when the wave packets are no longer Gaussian, which ensure the accuracy of the FGS method in computing non-Gaussian wave packets. In future work, we are ready to apply the FGS to approximating more complicated wave packets arising in scientific problems.

Examples with WKB initial condition
In this subsection, we apply the FGS with the initial sampling algorithm put forward in Section 4 to two examples of Equation (1.2) with WKB initial data. Caustics forms in the former example and does not form in the latter, which makes the former example more challenging for the FGS. To avoid confusion, we would like to point out that caustics become relevant to Equation (1.2) only in the semiclassical limit, ε → 0, but it also causes numerical challenges for moderately small ε.
Example 3: WKB initial data example with caustics.
In the WKB initial condition We choose the constant potential E 5 (x) = 10. Due to the compressive initial velocity d dx S in , caustics will form. This example was used in [23,1,14] as a reference to test whether certain algorithms can capture caustics.
In this example, we plot the FGS wave function and the resulting FGS current density function and compared them with the reference solutions computed via the SP2 method with sufficiently small step lengths. The computation is carried out on spatial interval [−π, π]. Similarly, we use the fourth order Runge-Kutta method to solve the ODE system (2.3)-(2.6) in the time interval [0, 0.54] with time step length ∆t = 0.01 and construct the wave function and current density at t = 0.54. We fix the scaling Planck parameter as ε = 0.0016 and plot the FGS wave functions and current density functions with ensemble sizes M = 400, 800, 1600 in Figure 5. The wave functions and current density functions computed via the FGS can reproduce the same profile as the reference solutions, which ensures the sampling strategy proposed in Section 4 is reliable and of low computational cost to approximate Equation (1.2) with the WKB initial conditions even in the presence of caustics.

Example 4: WKB initial data example without caustics
Now that we have shown the reliable results of the FGS in computing complicated WKB initial data problems with caustics, we further quantitatively study the sampling we obtain the following WKB initial data: We consider the null potential function E 6 (x) = 0 and the quadratic potential function These examples were used in [14] and caustics does not form. Similar to Example 1 in Section 5.1, we aim to compute the sampling error E S defined in Equation (5.2). The FGA ansatz u FGA in Equation (5.2) is approximated by a FGS wave function with sufficiently large ensemble size M 0 = 1.5 × 10 5 as shown in Equation (5.3). We use the fourth order Runge-Kutta method to solve the ODE system (2.3)-(2.6) in the time interval [0, 0.5] and reconstruct the wave function u FGS at t = 0.5. The computation is carried out on spatial interval [−π, π].
In Table 2 and Figure 6, we present the sampling errors E S (M, ε) of the FGS with for different scaling Planck parameter ε = 0.0256, 0.0064, 0.0016, 0.0004 and two potential functions E 6 and E 7 respectively. Similar to Example 1, the convergence of the sampling error E S is of order O(M − 1 2 ) as M → +∞. Moreover, we observe that the FGS sample errors of the WKB initial data cases mildly increase as the scaling parameter ε → 0 for both potential functions. Nevertheless, we underline that the increase of E S is insensitive to ε, which is manifested by the results in Table 2: when ε becomes 1 64 times smaller, the errors grow at most a couple times larger. Compared to the mesh-based algorithm whose computational cost is of order O(ε −1 ), the increase of the sampling error of the FGS is rather acceptable.

Computing physical observables in high dimensions
Example 5: High dimensional example In this subsection, we aim to compute the position observable and momentum observable as follows: where u denotes the exact solution to semiclassical Schrödinger equation (1.2). According to Section 2.3, the FGS provides a mesh-free algorithm to compute the position observable and momentum observable without reconstructing the FGS wave function u FGS . The physical observables computed via the FGS are denoted as follows: where we refer to Equation (2.21) for the computing expressions of O q FGA and O p FGA . We compute the mean-square error of the quantities of the physical observables: where we view the observables O q FGS and O p FGS computed via the FGS as a random variables, whose probabilistic distribution is determined by the FGS random samples z (j) 0 m j=1 on the phase space. We aim to discuss the mean-square errors E q and E p with respect to the semiclassical parameter ε and the ensemble size M .
To approximate the mean-square errors, we compute FGS are sampled by independent phase space samples. We choose the ensemble size M = 60 in our tests. To compute the exact quantities O q and O p of the position and momentum observables which are analytically unavailable, we consider a special case when the quantities of O q and O p can be given explicitly by a semiclassical Gaussian wave packet named the Hagedorn wave packet [10,11,30,17], which can be approximated numerically through evolving ODEs. For this special case, the position and momentum observables involved to the wave packet can be derived explicitly through the wave packet parameters. We directly quote Proposition 3.18 in [17] as the theoretical foundation of the Hagedorn wave packet method: Theorem 5.1. Consider the semiclassical Schrödinger equation (1.2) with a quadratic potential E(x) = |x| 2 2 . Let vectors q H (t), p H (t) ∈ R m and matrices Q H (t), P H (t) ∈ C m×m be the solution to the following ODE system: ) ds be the corresponding action function. Then the following complex Gaussian wave packet Moreover, the initial amplitude (a 1 , a 2 , · · · , a m ) ∈ R m is given as follows a j = 1 + {0.2j}, j = 1, 2, · · · , m, (5.27) where {x} denotes the fractional part of a rel number x.
Remark 5.2. As shown in Equation (3.12), the sampling error of the FGS is relevant to the term m j=1 1+a j √ a j 1 2 , which is related to the initial amplitude. To avoid the influence resulting from the initial conditions, we choose a j to be restricted in the interval [1,2).
We aim to compute the mean-square error E q and E p given by Equation ( The computation is carried out on spatial interval [−π, π] and the scaling parameter is fixed to ε = 0.0256. In figure 7, we plot the position and momentum observable errors E q and E p with different ensemble sizes M and dimension number m = 3, 4, · · · , 7. For fixed m, the convergence rate of E q and E p is roughly O(M −2 ), whereas the convergence slows down when M is large for m = 3. This is due to the fact that, when M is sufficiently large, the Monte Carlo sampling errors becomes much smaller than the asymptotic errors of the FGA ansatz that dominate the observable errors E q and E p . Moreover, we observe that, with fixed ensemble size M , the errors of the physical observables mildly increase as the dimension number m increases. As is manifested in , in order to reached a uniformly bounded mean-square error, the sample size should quadruple as the dimension number increases by one. Therefore, it is recommended that the sample size to compute physical observable in a m-dimensional problem should be of the order O(4 m ).
Furthermore, we record in Table 4 the system time to implement the FGS to compute an m−dimensional physical observables O q FGA or O p FGA . Since the complexity to compute the position and momentum observables are almost the same, we only present the system time to compute O q FGA in Table 4. To reached a uniformly bounded mean-square error, the system time increases at five times as m increases by one.

Conclusion and discussion
In this work, we have developed a novel stochastic simulation methodology for the semiclassical Schrödinger equation that shows great promises for high dimensional simulations. Its numerical cost is insensitive to the semiclassical parameter, grows gently with respect to the dimension number and it provides an efficient way to compute the physical observables without reconstructing the wave function. The core component of the algorithm is to properly choose a density function, that achieves a balance between variance reduction and sampling convenience, and should be specified for different types of initial wave functions. Although we have already obtained satisfactory results for the Gaussian initial conditions, and the WKB initial conditions, more advanced sampling techniques, in theory, can be utilized in the framework of the FGS methodology, and hence further improve its performance. Besides, the idea of turning computing highly oscillatory waves into sampling and propagating wave packets may also be applied to other wave equations in the high frequency regime. Thus, this work initiates many directions that we shall explore in the future.
With the characteristic flow X t : y → x W (r, y), the ODE system (6.7) yields a phase function as follows: Note that the characteristic flow X t : y → x W (r, y) may not be one-by-one, hence the phase function S W given in Equation (6.8) may be smooth only in a small time interval when characteristic flow X t : y → x W (r, y) doesn't cross. The cases when X t ceases to be a diffeomorphism are called caustics phenomena. See Figure 2 in [7] for an illustration of the caustics characteristic. Before caustics onset, the WKB ansatz in Equation (6.4) is an asymptotic approximation to the exact solution to Equation (1.2) of order O(ε), provided that the initial considition is exact and the amplitude a in decays rapidly (see Theorem 2.3 in [13]). After caustics forms, the exact solution may no longer be approximated by the WKB ansatz (6.4), hence the WKB approximation breaks down. The readers may refer to [2] for detailed discussion on the WKB approximation and analytic techniques after caustics onset.