# Gaussian wave packet transform based numerical scheme for the semi-classical Schrödinger equation with random inputs
^{1}^{1}1
S. Jin was supported by NSF grants DMS-1522184 and DMS-1107291: RNMS KI-Net, and NSFC grants No. 31571071 and No. 11871297.
L. Liu was partially supported by the funding DOE–Simulation Center for Runaway Electron Avoidance and Mitigation, project No. DE-SC0016283.
G. Russo was partially funded by the ITN-ETN Marie-Curie Horizon 2020 program
ModCompShock, Modeling and computation of shocks and interfaces, project ID: 642768 and NSF grant No. DMS-1115252.
Z. Zhou was partially supported by a start-up fund from Peking University and NSFC grant No. 11801016.

###### Abstract

In this work, we study the semi-classical limit of the Schrödinger equation with random inputs, and show that the semi-classical Schrödinger equation produces oscillations in the random variable space. With the Gaussian wave packet transform, the original Schrödinger equation is mapped to an ODE system for the wave packet parameters coupled with a PDE for the quantity in rescaled variables. Further, we show that the equation does not produce dependent oscillations, and thus it is more amenable for numerical simulations. We propose multi-level sampling strategy in implementing the Gaussian wave packet transform, where in the most costly part, i.e. simulating the equation, it is sufficient to use independent samples. We also provide extensive numerical tests as well as meaningful numerical experiments to justify the properties of the numerical algorithm, and hopefully shed light on possible future directions.

## 1 Introduction

In simulating physical systems, which are often modeled by differential equations, there are inevitably modeling errors, imprecise measurements of the initial data or the background coefficients, which may bring about uncertainties to the equation. There has been a growing interest in analyzing such models to understand the impact of these uncertainties, and thus design efficient numerical methods.

When it comes to quantum dynamics, quantifying the effect of uncertainty is even a trickier task. The solution to the Schrödinger equation is a complex valued wave function, whose nonlinear transforms (e.g. position density, flux density) lead to probabilistic measures of the physical observables. Thus, the uncertainty in the Schrödinger equation may or may not result in changes in measurable quantities from the quantum state.

We consider the following semi-classical Schrödinger equation with random inputs

(1.1) |

(1.2) |

Here, is the semi-classical parameter, which is reminiscent of the scaled Plank constant, and is the scalar potential function, which is slow-varying and often used to model the external field. The initial condition will be assumed to be in the form of a semi-classical wave packet, which is a Gaussian wave packet parameterized by the wave packet position, the wave packet momentum, etc.

The uncertainty is described by the random variable , which lies in the random space with a probability measure . We introduce the notation for the expected value of in the random variable ,

(1.3) |

In this paper, we only consider the uncertainty coming from initial data and potential functions, that is the uncertainty is classical. For example, the external classical field, the wave packet position or the wave packet momentum is uncertain, which reflects on the uncertainty in physical observables.

We do not, however, aim to analyze different types of uncertainties in quantum dynamics, but rather, we study how the uncertainty propagates in the semi-classical Schrödinger equation. When , it is well known that the Schrödinger equation is in the high frequency regime, where the solution generates scaled oscillations in space and time. As we shall show in this paper, the solution generically propagates scaled oscillations in the variable as well even if the random variable obeys an independent probability distribution. The high frequency of the solution in space, time and uncertainty leads to unaffordable computational cost, which makes conventional numerical methods infeasible. Thus, it is of great interest to design efficient numerical method based on the multiscale nature of the analytical solutions.

In the semi-classical regime, due to the scaled oscillation in the solution to the Schrödinger equation, the wave function does not converge in the strong sense as . The high frequency nature of the wave function of the semi-classical Schrödinger equation also causes significant computation burdens. If one aims for direct simulation of the wave function, one of the best choices is the time splitting spectral method, as analyzed in [2] by Bao, Jin and Markowich. See also [10, 12, 21, 3], where the meshing strategy and is sufficient for moderate values of . Another advantage of the time splitting methods is that if one is only interested in the physical observables, the time step size can be relaxed to , in other words, independently of , whereas one still needs to resolve the spatial oscillations. As we shall show in the paper, when the uncertainty is present, the wave function is also highly oscillatory in the variable, which means the size of samples in random variable grows as in order to obtain accurate approximations of the quantum dynamics.

There are quite a few approximate methods other than directly simulating the semi-classical Schrödinger equation, which are valid in the limit , such as the level set method and the moment closure method based on the WKB analysis and the Wigner transform, see, for example, [10] for a general discussion. In the past few years, many wave packets based methods have been introduced, which reduce the full quantum dynamics to Gaussian wave packets dynamics [9, 16, 8], and thus gain significant savings in computation cost, such as the Gaussian beam method [25, 11, 24], the Hagedorn wave packet approach [5, 31] and the Frozen Gaussian beam method [13, 14, 18, 20, 19]. When random inputs are considered, in theory, one can design numerical methods based on those approximation tools with independent samples in , however, the approximation errors persist in spite of other potential challenges.

A related work to our current subject is [4], where the authors developed the generalized polynomial chaos (gPC)-based stochastic Galerkin method for a class of highly oscillatory transport equations containing uncertainties that arise in semi-classical modeling of non-adiabatic quantum dynamics. Built upon and modified from the nonlinear geometrical optics based method, this scheme can capture oscillations with frequency-independent time step, mesh size as well as degree of the polynomial.

Clearly, an exact reformulation of the semi-classical Schrödinger equation that separates the multiscales in dynamics is desired for efficient simulation. Very recently, Russo and Smereka proposed a new method based on the so-called Gaussian wave packet transform [26, 27, 32], which reduces the quantum dynamics to Gaussian wave packet dynamics together with the time evolution of a rescaled quantity , which satisfies another equation of the Schrödinger type, in which the modified potential becomes time dependent. We emphasize that the Gaussian wave packet transform is an equivalent reformulation of the full quantum dynamics, and there are no more dependent oscillations in the equation. This motivates us to investigate whether this transform facilitates design of efficient numerical methods for semi-classical Schrödinger equations with random inputs.

In this work, we study the semi-classical limit of the Schrödinger equation with random inputs, which is the Liouville equation with a random force field, and show that the semi-classical Schrödinger equation produces oscillations in the variable in general. However, with the Gaussian wave packet transform, the original Schrödinger equation is mapped to an ODE system for the wave packet parameters coupled with a PDE for the quantity in rescaled variables, where the ODE system and the equation also depend on the random variable. Further, we show that the equation does not produce dependent oscillations in the rescaled spatial variable, thus it is more amenable for numerical simulations. We propose multi-level sampling strategy in implementing the Gaussian wave packet transform, where in the most costly part, simulating the equation, it is sufficient to use independent samples, thus the complexity of the whole algorithm has satisfactory scaling behavior as goes to . We also provide extensive numerical tests as well as meaningful numerical experiments to justify the properties of the numerical algorithm as well as hopefully shed light on possible future directions.

The rest of the paper is outlined as follows. Section 2 discusses the semi-classical limit of the Schrödinger equation and analyzes the regularity of in the random space. In Section 3, we introduce the Gaussian wave packet transform and prove that the equation is not oscillatory in the random space. Section 5 shows extensive numerical tests by using the stochastic collocation method to demonstrate the efficiency and accuracy of our proposed scheme. Relations of different numbers of the collocation points needed in each step of the implementation will be explained and studied numerically. Conclusion and future work are given in Section 6.

## 2 The semi-classical Schrödinger equation with random inputs

### 2.1 The semi-classical limit with random inputs

In this part, we investigate the semi-classical limit of the Schrödinger equation with random inputs by the Wigner transform [6, 10, 1, 17]. Obviously, the potential function can be decomposed as

(2.4) |

such that

For , the Wigner transform is defined as a phase-space function

(2.5) |

Recall that is the exact solution of equation (1.1). Denote it is possible to prove that satisfies the Wigner equation [6, 2]

(2.6) |

in which is a pseudo-differential operator acting on defined by

(2.7) |

By Weyl’s Calculus (see [6]) as , the Wigner measure satisfies the classical Liouville equation

(2.8) |

with

(2.9) |

All the limits above are defined in an appropriate weak sense (see [6, 17]).

With the decomposition of the potential as in (2.4), the classical Liouville equation becomes

(2.10) |

This clearly shows, due to the random potential, the bi-characteristics of the Liouville equation contains the random force term with , and the characteristic equations are

(2.11) |

Also, by definition, it is easy to check that , is real-valued. To sum up, in the semi-classical limit, the Wigner measure picks up the dependence of the random variable though the initial condition and the vector field .

Next, we discuss if one can derive the averaged equation to integrate out the random variable . In equation (2.10), by taking the average with respect to the random variable , one gets

Notice that, in the last term, both and depend on , thus it cannot be directly written as a term involving , rather it connects with the covariance of and . In fact,

since .

To illustrate the effect of the random potential, we consider the following special case

In this case, the Liouville equation simplifies to

(2.12) |

Then, if one considers the following change of variables

then, equation (2.12) becomes

(2.13) |

And correspondingly, the initial condition becomes

Thus, the average with respect to the random variable can be taken and one gets a closed equation for ,

(2.14) |

and

This example shows that the randomness in the slow-varying potential changes the transport part of the Liouville equation, although in the averaged equation the transport structure may be even unchanged.

When the potential is random and with a general form of dependence, the two processes of a) first pushing then taking expected value in of the classical limit (the Liouville equation); and b) first taking expected value in on the Wigner equation then letting do not commute. Though numerically, our simulation results seem to suggest the commutation of these two iterated processes. In later section we will make a numerical comparison to show that the expected values of the position density and flux do converge in the limit.

To conclude this part, we remark that in [1, 23] and subsequent works, the authors have considered apparently a related but fundamentally different model, where the unperturbed system is the semi-classical Schrödinger with fast-varying smaller magnitude and random perturbation in the potential is also fast-varying. And they show that the randomness in that scaling introduces additional scattering terms in the limit equations, while in our case the randomness only persist in the initial data and the force field in the limit equation.

### 2.2 Regularity of in the variable

The semi-classical Schrödinger equation is a family of dispersive wave equations parameterized by , and it is well known that the wave equation propagates scaled oscillations in space and time. However, it is not clear yet whether the small parameter induces oscillations in the random variable .

Here and in subsection 3.2, we will conduct a regularity analysis of in the random space, which enables us to study the oscillatory behavior of solutions in the random space, which gives guidance on how many collocation points needed in each step of the collocation method should depend on the scaled constant .

To investigate the regularity of the wave function in the variable, we check the following averaged norm

(2.15) |

To be more precise, (2.15) denotes the square root of the expected value in of the square of the norm of . We name it -norm for short.

One first observes that ,

thus

which means the -norm of the wave function is conserved in time,

However, we show in the following that has -scaled oscillations in even if and do not have -dependent oscillations in . We first examine the first-order partial derivative of in , and denote and , then by differentiating the semi-classical schrödinger equation (1.1) with respect to , one gets

By direct calculation,

where the Cauchy-Schwarz inequality and the Jensen inequality are used in the last step, more specifically,

Thus

For , the pessimistic estimate implies

Moreover, for , denote , we can similarly conclude that

(2.16) |

Although the estimates above are apparently pessimistic, we would like to show that the high frequency oscillations in can be seen in the following example. For simplicity, we consider be one dimensional. If the potential is quadratic in , it has been shown by Heller in [9] that,

(2.17) |

is an exact solution to the semi-classical Schrödinger equation, provided that, , , and satisfy the following system of equations

(2.18) |

Due to the same reason,

is an exact solution to equation (1.1), when potential is quadratic in . Clearly, this specific solution saturates the estimate (2.16), which implies, even if initially is smooth in , it will pick up -dependent oscillations in the variables. In Section 5, we will also show numerically the -scaled oscillations of such wave functions in the variable.

To conclude this section, we emphasize the numerical challenges with respect to the random variable . Due to the oscillatory behavior in , if one applies the generalized polynomial chaos (gPC)-based stochastic methods directly to the semi-classical Schrödinger equation, one needs at least -dependent basis functions or quadrature points to get an accurate approximation. The stochastic collocation method will be discussed in detail in subsection 5.1.

## 3 The semi-classical Schrödinger equation and the Gaussian wave packet transformation

To overcome the numerical burdens in sampling the random variable , we introduce the Gaussian wave packet transformation (abbreviated by GWPT), which has been proven to be a very efficient tool for computing the Schrödinger equation in the high frequency regime. In essence, the GWPT equivalently transforms the highly oscillatory wave equation to , thus facilitate the design of efficient numerical methods. The main goal of this section is to study the uniform regularity of in the variable, and conclude that the number of basis functions or quadrature points is independent of , if applying the SG or SC method in our GWPT framework.

### 3.1 Review of the Gaussian wave packet transformation

First, we briefly summarize the Gaussian wave packet transformation applied to the semi-classical Schrödinger equation with random inputs, which is a natural extension of the GWPT method for the deterministic problem [26]. Consider the semi-classical Schrödinger equation given by (1.1)–(1.2). Note that in this work, random inputs are assumed to be classical, thus we only consider the cases when the wave packet position and momentum in the GWPT parameters depend on the random variable .

We start by the following ansatz

(3.19) |

where , is a real-valued symmetric matrix and is a complex-valued scalar.

Denote , and . Insert the ansatz (3.19) into (1.1), then satisfies

provided , , satisfy the following equations

(3.20) |

while

At last, introduce the change of variables , where

(3.21) |

with

then and

(3.22) |

Note that in the variable,

so the equation (3.22) is not oscillatory in nor in . Furthermore, if one drops those terms, one expects to recover the leading order Gaussian beam method [27].

In our numerical tests, we will only consider the initial data given by a Gaussian wave packet, i.e.

(3.23) |

More general initial conditions, (with a numerical support proportional to ) can be approximated with the desired accuracy as a superposition of relatively small number of Gaussian wave packets. See [26, Section 2.6] for more details.

### 3.2 Regularity of in the variable

It is well understood that the equation no longer propagates -dependent oscillations in space or in time. We show in the following that the equation is not oscillatory in the variable either. In this section, we assume the spatial variable is one dimensional to simplify the analysis. Now the equation (3.22) reduces to

(3.24) |

where

and . Observe that is a function of independent variables , but it obtains the dependence of the parameter through the coefficients. We emphasize that, although the change of variable (3.21) is dependent, in the equation, and are independent variables. This is due to the fact that it is , not . With the random inputs, the Gaussian wave packet transform is straightforward for all . Besides the time dependence, the Gaussian wave packet parameters , , , and also depend on . The smooth components and depend on . Then, it is not yet clear whether has dependent derivatives.

We make the assumption that the potential is infinitely smooth with bounded derivatives in both and , namely, for and , there exists a constant such that,

(3.25) |

where . We also assume that the equation (3.24) is equipped with an initial condition, , which has a sized support, satisfying the following assumption: for and , there exists an -independent constant such that,

(3.26) |

Our goal is to show that the smoothness in the variable will be preserved in time. Before calculating the regularity of in the variable, we first show Lemma 3.1 and Lemma 3.2, which will be used for the main result of this section namely Theorem 3.1.

###### Lemma 3.1.

Assume boundedness condition (3.25), and initial data for wave packet parameters satisfying the following: there exists an independent constant , such that

Then, there exists an independent constant , such that for ,

The proof follows standard estimations of the ODE system of the parameters, which we shall omit here. We also remark that we have only listed the parameters needed for showing the regularity properties of the equation (3.24), but this argument clearly works for other wave packets parameters as well.

###### Lemma 3.2.

With the boundedness assumptions (3.25), for all , it is

###### Proof.

Recall that

(3.27) |

We observe that, assumption (3.25) together with the Taylor’s Theorem implies, for and ,

(3.28) |

Thus, it is clear that .

Next, we examine the first order derivative in . By direct calculation,

We give the definition of the averaged norm of ,

which is analogous to the –norm defined in (2.15) for while using variable in the norm here. We now present the main theorem of this section:

###### Theorem 3.1.

###### Proof.

The equation can be written as

(3.30) |

where

and given in (3.27). We first look at the -norm of ,

which implies the averaged norm of is preserved, .

Differentiating (3.30) with respect to , by the chain rule, one has

By direct calculation (omit the -subscript in the norm for notation simplicity),

where integration by parts and the Cauchy-Schwarz inequality are used. Thus

(3.31) |

Clearly, to prove the boundedness of , it suffices to show the boundedness of the right hand side of (3.31). Fortunately, the right hand side of (3.31) does not involve the derivative of . The estimates of the derivatives of are standard, which can be carried out in the following deductive way.

We calculate that

(3.32) |

Since , and by the chain rule,

where are constants. Thus

(3.33) |

For , the boundedness of follows from the Grönwall’s inequality and assumption (3.26). Similarly, we get

where the Cauchy-Schwarz inequality is used again just as in (3.32). By the chain rule, such that

thus

(3.34) |

The boundedness of follows from Grönwall’s inequality and assumption (3.26), and similarly for . , are also bounded. By Lemma 3.2, . Using Grönwall’s inequality on (3.31), one gets

where is a constant, independent of . By induction, we obtain

where . Therefore, we have shown Theorem 3.1. ∎

## 4 Quantum and classical uncertainty

In this section we briefly discuss about quantum and classical uncertainty, and about the comparison between quantum and classical systems, for small values of the rescaled Planck’s constant. For simplicity, we first consider the case with one degree of freedom, , and scalar random variable .

### 4.1 Moments and expectations

A quantum system is completely determined by its wave function . For each realisation of the random variable , the quantum system is described by .

The primary physical quantities of interest include the position density,