Computation of the Schrödinger equation in the ... - Semantic Scholar

¨ COMPUTATION OF THE SCHRODINGER EQUATION IN THE SEMICLASSICAL REGIME ON UNBOUNDED DOMAIN∗ XU YANG† AND JIWEI ZHANG‡ Abstract. The study of this paper is two-fold: On the one hand, we generalize the high-order local absorbing boundary conditions (LABCs) proposed in [J. Zhang, Z. Sun, X. Wu and D. Wang, Commun. Comput. Phys., 10 (2011), pp. 742–766] to compute the Schr¨ odinger equation in the semiclassical regime on unbounded domain. We analyze the stability of the equation with LABCs and the convergence of the Crank-Nicolson scheme that discretizes it and we conclude that when the rescaled Planck constant ε gets small, the accuracy deteriorates and the requirements on time step and mesh size get tough. This leads to the second part of our study. We propose an asymptotic method based on the frozen Gaussian approximation. The absorbing boundary condition is dealt by a simple strategy that all the effects of the Gaussian functions which contribute to the outgoing waves will be eliminated by stopping the Hamiltonian flow of their centers when they get out of the domain of interest. We present numerical examples in both one and two dimensions to verify the performance of the proposed numerical methods. Key words. Schr¨ odinger equation, semiclassical regime, absorbing boundary condition, CrankNicolson scheme, frozen Gaussian approximation AMS subject classifications. 65M12, 65M06, 65M15

1. Introduction. We are interested in developing efficient numerical methods for computing the Schr¨ odinger equation on unbounded domain iε∂t ψ ε (x, t) = −

ε2 ∆ψ ε (x, t) + V (x, t)ψ ε , 2

x ∈ Rd ,

0 < t ≤ T,

(1.1)

with the initial condition ψ ε (x, 0) = ψ0ε (x).

(1.2)

Here ψ ε (x, t) is the wave function, V (x, t) represents the external potential, ε (0 < ε ≪ 1) is the rescaled Plank constant that describes the ratio between the quantum and macroscopic time/space scales, and d is the dimensionality. We assume ψ0ε is an L2 function and has compact support. The numerical computation of (1.1)-(1.2) has two difficulties: 1. the unboundedness of the domain; 2. the oscillatory feature of the solution when ε is small. To overcome the first difficulty, one powerful tool is to use artificial boundary conditions to reformulate the problem appropriately on a bounded domain ([45, 14, 15]). The key idea is to develop suitable (ideal) boundary conditions to absorb the waves arriving at artificial boundaries. In the literature, much attention has been given to study the artificial boundary conditions (ABCs) of the Schr¨odinger equation, which includes implicit and explicit, local and global ABCs ([3]). For example, Antoine and Besse considered the no-reflected Dirichlet-to-Neumann (DtN) or Neurman-toDirichlet (NtD) boundary conditions for one-dimensional Schr¨odinger equation ([1]). ∗ X.Y. was partially supported by the Regents Junior Faculty Fellowship of University of California, Santa Barbara. J.Z. was partially supported by NSF grant DMS-1009575. † Department of Mathematics, University of California, Santa Barbara, CA 93106, USA ([email protected]). ‡ Department of Mathematics, Courant Institute of Mathematical Sciences, New York University, New York, NY 10012, USA ([email protected]).

1

2

X. YANG AND J. ZHANG

Later Antoine et al extended their method to the two-dimensional case ([2]). Han and Huang obtained the exact artificial boundary conditions for the Schr¨odinger equation in high dimensions ([16, 17]). Recently, motivated by the idea of Kuska ([26]), Zhang et al constructed a number of high-order local absorbing boundary conditions (LABCs) by introducing a physical parameter k0 to optimize the performance of high-order LABCs ([48]). The parameter k0 is related to the frequency of wave impinged on artificial boundaries, which is different from the ABCs obtained by Di Menza ([7]), Fevens and Jiang ([11]), Szeftel ([42]) and Antoine et al ([4]). The influences and adaptive choices of k0 were discussed in [46]. The second difficulty arises when ε ≪ 1, which is called the semiclassical regime. In this case, the wave function ψ ε becomes oscillatory of the wave length O(ε). This means one has to work on a large computational domain that contains thousands to millions of wavelengths, and each of them needs to be resolved if direct numerical methods are applied. For example, a mesh size of O(ϵ) is required when using the timesplitting spectral method ([5]) to compute (1.1)-(1.2) directly; an even worse mesh size of o(ε) is needed if one uses the Crank-Nicolson scheme ([35]) or the Dufort-Frankel scheme ([36]). Moreover, in these methods, a large domain is demanded in order to avoid the boundary effects. Therefore the total number of grid points is huge, which usually leads to unaffordable computational cost. For the applications of high-order LABCs to wave equations, we refer the interested readers to [8, 13, 6]. Semiclassical approximation provides an efficient tool to reduce the computational cost. It studies the mathematical limit of the solutions when the wavelength O(ε) goes to zero in (1.1)-(1.2). Based on that, asymptotic methods have been successfully applied, among which, the geometric optics (GO) and the Gaussian beam method (GBM) are the most well-known currently ([10, 12, 19, 20, 21, 22, 23, 24, 28, 29, 30, 31, 37, 38, 39, 40, 43, 44]). Recently, the frozen Gaussian approximation (FGA) method was developed for general linear strictly hyperbolic systems with smooth coefficients (e.g. acoustic wave equations) in [32, 33, 34], which was motivated by the Herman-Kluk propagator in quantum chemistry ([18]). The method uses Gaussian functions with fixed widths to approximate high frequency waves. Compared to GO and GBM, the FGA method performs more reliably, i.e. on the one hand, it provides an accurate asymptotic solution in the presence of caustics where GO fails to, and on the other hand, it resolves the issue of losing accuracy as beams in GBM spread in the time evolution. We refer interested readers to [21, 32] for the detailed comparisons of these methods. In this paper, we first generalize the high-order LABCs introduced in [48] to compute (1.1)-(1.2) on a bounded domain. We analyze the stability of the Schr¨odinger equation with LABCs, and design the Crank-Nicolson scheme accordingly. The convergence of the numerical scheme will be proved. In particular, we conclude that when ε ≪ 1, in order to get accurate solution, the time step and mesh size have to be chosen small compared to ε, which makes the computational cost very expensive. This motivates us to develop numerical methods based on semiclassical approximation. Specifically, we propose an asymptotic method based on FGA for the computation of (1.1)-(1.2) on the unbounded domain, where the LABCs are replaced by a simple strategy that all the effects of the Gaussian functions which contribute to the outgoing waves will be eliminated by stopping the Hamiltonian flow of their centers when they get out of the domain of interest. The rest of the paper is organized as follows. In Section 2, we construct the high-order LABCs for computing (1.1)-(1.2) on the unbounded domain. In Section 3,

Computation of the Schr¨ odinger equation on unbounded domain

3

the resulted equation with LABCs will be solved by the Crank-Nicolson method, and the stability and convergence of the numerical scheme will be analyzed. We introduce the asymptotic method based on FGA in Section 4. Two numerical examples in one and two dimensions will be presented in Section 5 to verify the performance of the numerical methods. We make some conclusive remarks in Section 6. 2. High-order absorbing boundary conditions. In this section, we systematically construct the high-order absorbing boundary conditions (ABCs) for the computation of (1.1)-(1.2) in one dimension. We shall omit the superscript ε in ψ ε for simplicity without causing any confusion. Assume ψ0 in (1.2) is compactly supported ¯ i, on Ωi := [xl , xr ], and the potential V (x, t) in (1.1) is constant outside Ωe = R \ Ω ¯ i denotes the closure of the set Ωi . These assumptions imply that the waves where Ω in the computational domain only propagate from the interior domain into exterior domain Ωe through the artificial boundaries Γ = {xl , xr } and there are no waves traveling from the exterior domain to the interior. Note that the (boundaries) should be almost transparent for a plane wave of the form ψ(x, t) = exp i(ξx − ωt) , where ω and ξ denote the frequency and wave number respectively. Substituting the plane wave into (1.1) gives the dispersion relation ε2 2 ξ = εω − V, 2

(2.1)

which implies √ ξ=±

2√ εω − V , ε

(2.2)

where “±” correspond to ABCs at the boundaries xr and xl respectively. By the inverse Fourier transformation, the classical nonlocal boundary conditions (the DtN map) are given by ∂n ψ(x, t) = −

e

√ πi itV (x) 24− ε

√ ε

d dt

∫

t

0

iτ V (x)

ψ(x, τ )e ε √ t−τ

dτ,

x ∈ Γ = {xl , xr },

(2.3)

where ∂n stands for the outwardly directed normal derivative. To obtain the local ABCs, we use the framework of Engquist and Majda ([9]) for hyperbolic wave propagation models. This method has also been applied to the Sch¨odinger √ equation, and we refer to [3] for more details. Denoting Z = εω − V and expanding Z around some point Z0 in Pad´e approximation produce √

√ Z0

N √ √ ∑ Z bm (Z0 − Z) ≈ Z0 − Z0 , Z0 Z − am (Z0 − Z) m=1 0

for

1 − Z < 1, Z0

(2.4)

where am = cos2 (

mπ ), 2N + 1

bm =

2 mπ sin2 ( ), 2N + 1 2N + 1

m = 1, 2, · · · , N.

Substituting (2.4) into (2.2) yields ( ) √ N ∑ bm (k02 − Z) 2k0 1− ξ=± , ε k 2 − am (k02 − Z) m=1 0

(2.5)

4

X. YANG AND J. ZHANG

√ where k0 = Z0 . The equation (2.5) is called the paraxial approximation in the literature, and will cause difficulties in the computation when N is large. To decrease the order of partial derivatives, we use the Lindmann’s trick ([27]) with some auxiliary variables defined by 1 ˆ ψ, (1 − am )k02 + am Z

ϕˆm,l =

where l = 1, 2 representing the left and right boundary auxiliaries respectively, and ϕˆm,l ψˆ are the Fourier transform of ϕm,l and ψ. The correspondences ξ ⇔ −i∂x and ω ⇔ i∂t give the following boundary conditions  ( ) √ N ( 2 ) ∑  0 i∂x ψ ± 2k ψ − b k ϕ − iε∂ ϕ + V ϕ = 0, m t m,l m,l 0 m,l ε m=1  (1 − am ) k02 ϕm,l − V am ϕm,l + iεam ∂t ϕm,l = ψ, m = 1, 2, · · · , N, (2.6) ¯ i with the LABCs on Γ, which produce the Schr¨ odinger equation on Ω ε2 iε∂t ψ(x, t) = − ∂x2 ψ(x, t) + V (x, t)ψ(x, t), xl < x < xr , 0 < t ≤ T, 2 [( ) ] √ N N ∑ ∑ 2k0 bm bm 2 i∂x ψ(xl , t) − 1+ ψ(xl , t) − k0 ϕm,1 (t) = 0, ε a a m=1 m m=1 m √ i∂x ψ(xr , t) +

[(

2k0 ε

1+

N ∑ bm a m=1 m

)

0 < t ≤ T, (2.8)

]

ψ(xr , t) − k02

(2.7)

N ∑ bm ϕm,2 (t) = 0, a m=1 m

0 < t ≤ T, (2.9)

− V (xl )am ϕm,1 (t) + iεam ∂t ϕm,1 (t) = ψ(xl , t),

1 ≤ m ≤ N, (2.10)

(1 − am ) k02 ϕm,2 (t) − V (xr )am ϕm,2 (t) + iεam ∂t ϕm,2 (t) = ψ(xr , t),

1 ≤ m ≤ N, (2.11)

(1 −

am ) k02 ϕm,1 (t)

ψ(x, 0) = ψ0 (x),

xl ≤ x ≤ xr ,

ϕm,1 (0) = 0,

ϕm,2 (0) = 0,

1 ≤ m ≤ N. (2.12)

The following theorem analyzes the stability of the above equations. Theorem 2.1. Let ψ(x, t) be the solution to (2.7)-(2.12) and denote

c=

N ∑ bm 2 a m=1 m N ∑ bm 1+ a m=1 m

,

(2.13)

then we have, for 0 < t ≤ T , ∫

xr xl

∫ xr N ( ) 2 εk 3 ∑ 2 2 |ψ0 (x)|2 dx. |ψ(x, t)|2 dx + √ 0 bm |ϕm,1 (t)| + |ϕm,2 (t)| ≤ eck0 t/ε 2 m=1 xl

The proof can be found in [48], and we omit the details here.

Computation of the Schr¨ odinger equation on unbounded domain

5

3. Crank-Nicolson method. In this section, we describe the Crank-Nicolson method to solve (2.7)-(2.12), and analyze the stability and convergence correspondingly. 3.1. Numerical scheme. We take the mesh size h = (xr −xl )/M and time step τ = T /K where M, K are positive integers. The domain [xl , xr ] × [0, T ] is discretized as Ωh × Ωτ : Ωh = {xj | xj = xl + jh, 0 ≤ j ≤ M },

Ωτ = {tn | tn = nτ, 0 ≤ n ≤ K}.

Denote {ψjn | 0 ≤ j ≤ M, 0 ≤ n ≤ K} to be the grid function on Ωh × Ωτ , and introduce the following notations 1 n 1 n 1 n− 1 n n n (ψj + ψj−1 ), δx ψj− (ψj − ψj−1 ), ψj 2 = (ψjn + ψjn−1 ), 1 = 2 2 2 h 2 1 1 n n− 21 n n δ t ψj = (ψj − ψjn−1 ), δx2 ψjn = 2 (ψj+1 − 2ψjn + ψj−1 ), h vτ   u M −1 u ∑ 2 1 1 u 2 n ψ n + |ψ n |2 . ∥ψ ∥ = th  |ψ0n | + (3.1) j M 2 2 j=1 n ψj− 1 =

Define the midpoints of the grids xj− 12 = (xj + xj−1 )/2,

n− 12

tn− 12 = (tn + tn−1 )/2,

Vj

= V (xj , tn− 12 ),

and two fictitious points x−1 = x0 − h, xM +1 = xM + h. Then (2.7)-(2.12) can be discretized as, for 1 ≤ m ≤ N , 1 ≤ n ≤ K, n− 12

iεδt ψj

ε2 2 n− 12 n− 1 n− 1 + Vj 2 ψj 2 , 0 ≤ j ≤ M, (3.2) δx ψj 2 [( ) ] √ N N ) ∑ ∑ bm n− 21 2k0 bm n− 1 n− 1 1+ ψ0 2 − k02 ϕm,1 = 0, + δx ψ 1 2 − 2 ε a a m=1 m m=1 m

=−

i( n− 1 δ x ψ− 1 2 2 2

(3.3) [( ) ] √ N N ( ) ∑ bm n− 1 ∑ bm i 2k0 n− 1 n− 1 n− 1 δx ψM −21 + δx ψM +21 + 1+ ψM 2 − k02 ϕm,22 = 0, 2 2 2 ε a a m m m=1 m=1 (3.4) (1 −

n− 1 am ) k02 ϕm,12 n− 21

−

n− 1 n− 1 V0 2 am ϕm,12 n− 12

(1 − am ) k02 ϕm,2 − VM ψj0

= ψ0 (xj ),

+

n− 1 iεam δt ϕm,12

=

n− 12

n− 12

n− 1 ψ0 2 , n− 12

am ϕm,2 + iεam δt ϕm,2 = ψM

0 ≤ j ≤ M,

ϕ0m,1

= 0,

ϕ0m,2

= 0.

,

(3.5) (3.6) (3.7)

At the n-th time level, (3.2)-(3.7) is a system of linear algebraic equations on ψjn and ϕnm,1 , ϕnm,2 . This Crank-Nicolson method is unconditionally stable and of the second-order accuracy in space and time. In the LABCs (3.3) - (3.7), the choice of parameter k0 is related to the frequency of the wave impinging on the artificial boundaries. When the order N of Pad´e expansion is large enough, the effectiveness of LABCs is not sensitive to the choice of k0 ([48]). We refer to [46] for the adaptive method of finding the value of k0 .

6

X. YANG AND J. ZHANG n− 1

Eliminating ψ−1 2 in (3.4) by (3.2) at j = 0 gives n− 12

iεδt ψ0

n− 12

= V0

ε2 n− 1 − δx ψ 1 2 2 2h [( ) ] √ N N ∑ ∑ 1 bm i 2εk0 b n− 12 n− m 1+ ψ0 ϕm,12 . − − k02 h a a m m m=1 m=1 n− 21

ψ0

(3.8)

Similarly, n− 21

iεδt ψM

n− 12

= VM

ε2 n− 1 + δx ψM −21 2 2h [( ) ] √ N N ∑ ∑ bm i 2εk0 bm n− 12 n− 12 2 1+ ψM − k0 ϕm,2 . − h a a m=1 m m=1 m n− 21

ψM

(3.9)

Then (3.2)-(3.7) become, for 1 ≤ m ≤ N , 1 ≤ n ≤ K, n− 21

iεδt ψj

n− 12

iεδt ψ0

ε2 2 n− 12 n− 1 n− 1 δ ψ + Vj 2 ψj 2 , 1 ≤ j ≤ M − 1, (3.10) 2 x j [( { ) ]} N N ∑ ∑ √ bm n− 12 1 ε2 bm n− 12 n− 12 2 =− 1+ − k0 δx ψ 1 + i 2εk0 ψ0 ϕm,1 2 h 2 a a m=1 m m=1 m =−

n− 12

n− 21

iεδt ψM

=−

1 h

{

√ i 2εk0

[( 1+

+ V0 )

N ∑

bm a m=1 m

n− 12

n− 12

(1 −

−

n− 1 n− 1 V0 2 am ϕm,12 n− 12

n− 12

(1 − am ) k02 ϕm,2 − VM ψj0

= ψ0 (xj ),

+

n− 12

ψM

n− 1 iεam δt ϕm,12

n− 12

,

− k02

ψM

+ VM n− 1 am ) k02 ϕm,12

n− 12

ψ0

=

n− 12

bm n− 12 ε2 n− 1 ϕm,2 − δx ψM −21 2 a 2 m=1 m ,

0 ≤ j ≤ M,

= 0,

ϕ0m,2

(3.12)

n− 1 ψ0 2 , n− 12

am ϕm,2 + iεam δt ϕm,2 = ψM ϕ0m,1

(3.11) }

]

N ∑

(3.13)

,

(3.14)

= 0.

(3.15)

3.2. Stability and solvability. Theorem 3.1. Denote {ψjn | 0 ≤ j ≤ M, 0 ≤ √ n ≤ K} as the solution to (3.2)-(3.7). When τ ≤ ε(6 2ck02 )−1 , we have 2

∥ψ n ∥ +

√ N ( 2T ∑ √ 6 2ck0 2 2 )

ψ 0 2 , 2εk03 bm ϕnm,1 + ϕnm,2 ≤ e ε

1 ≤ n ≤ K, (3.16)

m=1

where c is defined in (2.13). n− 1

n− 12

Proof. Multiplying (3.10) by hψ j 2 , (3.11) by 12 hψ 0 them, and taking the imaginary part producel

n− 1

, (3.12) by 21 hψ M 2 , adding

( )( N ) ∑ √

2 ) b 1 ( n 2 n− 21 2 n− 12 2 m n−1

= − 2k0 1 + ∥ψ ∥ − ψ ψ0 + ψM 2τ a m=1 m } { N N √ 3 n− 12 ∑ bm n− 12 ∑ bm n− 21 n− 12 ϕm,1 + ψ M ϕm,2 . (3.17) + 2k0 Re ψ 0 a a m=1 m m=1 m

7

Computation of the Schr¨ odinger equation on unbounded domain n− 1

Multiplying (3.13) by ϕm,12 gives n− 1 n− 1 n− 1 n− 1 n− 1 n− 1 2 n− 1 2 (1 − am ) k02 ϕm,12 − V0 2 am ϕm,12 + iεam ϕm,12 δt ϕm,12 = ϕm,12 ψ0 2 , the imaginary part of which yields { n− 1 } εam ( n 2 n−1 2 ) n− 1 ϕm,1 − ϕm,1 = Im ϕm,12 ψ0 2 . 2τ √ Multiplying the above equation by 2k03 bm /am and summing up m from 1 to N give } { √ N N ( 1 ∑ b 2 n−1 2 ) √ 3 2εk03 ∑ n− 12 n− m n 2 . (3.18) bm ϕm,1 − ϕm,1 = 2k0 Im ψ0 ϕ 2τ m=1 a m,1 m=1 m Similarly, from (3.14), { } √ N N ( ∑ n−1 2 ) √ 3 2εk03 ∑ bm n− 12 n− 12 n 2 bm ϕm,2 − ϕm,2 = 2k0 Im ψM ϕ . 2τ m=1 a m,2 m=1 m

(3.19)

Define N ( ∑ √ 2 2 ) 3 E = ∥ψ ∥ + 2εk0 bm ϕnm,1 + ϕnm,2 . n 2

n

m=1

Adding (3.17), (3.18) and (3.19) brings ( )( N ) ∑ ) bm 1 ( n n− 12 2 n− 21 2 n−1 E −E ≤ −k0 1 + ψ0 + ψM + B, 2τ a m=1 m

(3.20)

where ( ) N { } { } √ 3∑ bm n− 12 n− 12 n− 12 n− 12 n− 12 n− 12 n− 21 n− 21 B ≡ 2k0 Re ψ0 ϕm,1 + ψM ϕm,2 + Im ψ0 ϕm,1 + ψM ϕm,2 . a m=1 m √ 2 |Z| and the θ inequality imply ( N ) N ( ( ) k6 ∑ ) ∑ b n− 21 2 n− 12 2 n− 12 2 n− 21 2 m 0 B ≤ θ ψ0 + ψM + · bm ϕm,1 + ϕm,2 . θ m=1 a2m m=1

|Re(Z)| + |Im(Z)| ≤

( ∑N Taking θ = k0 1 + m=1 ( B ≤ k0

N ∑ bm 1+ a m=1 m

bm am

) gives

)( ( N ) ) ∑ n− 21 2 n− 12 2 n− 1 2 n− 1 2 bm ϕm,12 + ϕm,22 . ψ0 + ψM + ck05 m=1

(3.21) Inserting the inequality (3.21) into Eq. (3.20), one has ( N ) ∑ ) 1( n n− 1 2 n− 1 2 E − E n−1 ≤ 4ck05 bm ϕm,12 + ϕm,22 , τ m=1

8

X. YANG AND J. ZHANG

which produces √ 1 n 2 2ck02 n n−1 (E − E )≤ (E + E n−1 ). τ ε

√ When τ ≤ ε(6 2ck02 )−1 , the above equation gives √ 6 2ck02 τ n−1 n E ≤ (1 + )E , ε which implies En ≤ e

√ 2 nτ 6 2ck0 ε

E0.

This is equivalent to (3.16), which completes the proof. The above stability result implies the following theorem. Theorem 3.2. The difference scheme (3.2)-(3.7), or (3.10)-(3.15) is uniquely solvable. 3.3. Convergence. Define the function values on the grids as Ψnj = ψ(xj , tn ), Φnm,1 = ϕm,1 (tn ), Φnm,2 = ϕm,2 (tn ), 0 ≤ j ≤ M, 0 ≤ n ≤ K, 0 ≤ n ≤ K. Denote n− 12

pj

n− 1 p0 2

n− 12

= iεδt Ψj

=

n− 1 iεδt Ψ0 2

{ 2 } ε n− 1 n− 1 n− 1 − − δx2 Ψj 2 + Vj 2 Ψj 2 , 2

1 + h

{

√ ε n− 1 δx Ψ 1 2 + i 2εk0 2 2 2

[(

n− 12

n− 1 pM 2

=

n− 1 iεδt ΨM 2

1 + h

{

√ i 2εk0

[(

− V0 N ∑

bm 1+ a m=1 m

N ∑

bm 1+ a m=1 m n− 12

ψ0 )

n− 12

− VM

1 ≤ j ≤ M − 1, 1 ≤ n ≤ K, (3.22)

) n− 1 Ψ0 2

,

n− 1 ΨM 2

−

k02

−

k02

N ∑

bm n− 12 Φm,1 a m=1 m ]

N ∑

(3.23)

bm n− 12 ε2 n− 1 Φm,2 − δx ΨM −21 2 a 2 m=1 m

n− 1

ΨM 2 ,

(3.24)

n− 1

n− 1

n− 12

am Φm,1 + iεam δt Φm,12 − Ψ0

,

1 ≤ m ≤ N, (3.25)

n− 1

n− 1

n− 12

am Φm,2 + iεam δt Φm,22 − ΨM 2 ,

1 ≤ m ≤ N. (3.26)

rm,12 = (1 − am ) k02 Φm,12 − V0

rm,22 = (1 − am ) k02 Φm,22 − VM

]}

n− 1

n− 12

n− 1

n− 1

Then (2.7)-(2.11) and Taylor expansion imply there exists a constant c1 such that n− 21 (3.27) pj ≤ c1 (ετ 2 + ε2 h2 ), 1 ≤ j ≤ M − 1, 1 ≤ n ≤ K, n− 12 n− 1 (3.28) p0 ≤ c1 (ετ 2 + ε2 h), pM 2 ≤ c1 (ετ 2 + ε2 h), 1 ≤ n ≤ K, 1 1 n− n− 2 (3.29) rm,1 ≤ c1 ετ 2 , rm,22 ≤ c1 ετ 2 , 1 ≤ n ≤ K.

}

9

Computation of the Schr¨ odinger equation on unbounded domain

In the derivation of (3.28) we have used the facts that ε2 iε∂t ψ(x0 , t) = − ∂x2 ψ(x0 , t) + V (x0 , t)ψ(x0 , t) 2 { } ε2 1 [ψ(x1 , t) − ψ(x0 , t)] − ∂x ψ(x0 , t) + V (x0 , t)ψ(x0 , t) + O(εh) =− h h { [( ) ]} N N ∑ ∑ √ 1 ε2 bm bm 2 =− ψ(x0 , t) − k0 [ψ(x1 , t) − ψ(x0 , t)] + i 2εk0 1+ ϕm,1 (t) h h a a m=1 m m=1 m 0 < t ≤ T.

+ V (x0 , t)ψ(x0 , t) + O(εh),

Moreover, (3.22)-(3.26) and (2.12) yield n− 12

iεδt Ψj

n− 12

iεδt Ψ0

=−

=−

ε2 2 n− 12 n− 1 n− 1 n− 1 δx Ψj + Vj 2 Ψj 2 + pj 2 , 2

1 h

{

√ ε2 n− 1 δx Ψ 1 2 + i 2εk0 2 2

[( 1+

N ∑ bm a m=1 m

n− 12

n− 12

iεδt ΨM

=−

1 h

{

√ i 2εk0

[( 1+

N ∑

+ V0 )

bm a m=1 m

n− 12

n− 12

n− 1

n− 12

n− 1

n− 12

(1 − am ) k02 Φm,12 − V0

(1 − am ) k02 Φm,22 − VM Ψ0j = ψ0 (xj ),

n− 12

ψ0

ΨM

+ VM

1 ≤ j ≤ M − 1, 1 ≤ n ≤ K,

− k02 n− 12

ΨM

(3.30) ]}

) n− 12

Ψ0

n− 12

+ p0

(3.31) }

]

ε2 bm n− 21 n− 1 Φm,2 − δx ΨM −21 2 a 2 m=1 m n− 1

+ pM 2 ,

n− 1

n− 12

n− 1

n− 12

am Φm,2 + iεam δt Φm,22 = ΨM Φ0m,1 = 0,

,

N ∑

am Φm,1 + iεam δt Φm,12 = Ψ0

0 ≤ j ≤ M,

N ∑ bm n− 21 Φm,1 a m=1 m

− k02

(3.32)

n− 1

1 ≤ m ≤ N, (3.33)

n− 1

1 ≤ m ≤ N, (3.34)

+ rm,12 , + rm,22 ,

Φ0m,2 = 0.

(3.35)

Define ψ˜jn = Ψnj − ψjn , 0 ≤ j ≤ M, 0 ≤ n ≤ K; ϕ˜nm,1 = Φnm,1 − ϕnm,1 , ϕ˜nm,2 = Φnm,2 − ϕnm,2 ,

1 ≤ m ≤ N,

0 ≤ n ≤ K,

then we have the following theorem. Theorem 3.3. Denote {ψjn | 0 ≤ j ≤ M, 0 ≤ n ≤ K} as the solution to the difference scheme (3.2)-(3.7), or (3.10)-(3.15), then √ ( N ) ∑ √ 3(1 + 2ck02 )T ˜n 2 ˜n 2 n 2 3 ˜ ∥ψ ∥ + 2εk0 }ε(τ 2 +εh2 )2 , bm ϕm,1 + ϕm,2 ≤ c2 exp{ ε m=1 [√ where c2 =

2 2k0

√ ∑N + (xr − xl ) + 2 2ε2 k03 m=1

bm a2m

]

c21 .

1 ≤ n ≤ K,

10

X. YANG AND J. ZHANG

Proof. Subtracting (3.10)-(3.15) from (3.30)-(3.35) produces 2

ε n− n− n− n− n− iεδt ψ˜j 2 = − δx2 ψ˜j 2 + Vj 2 ψ˜j 2 + pj 2 , 2 1

1 n− iεδt ψ˜0 2 = − h

{

1

1

1

1

1

√ ε2 ˜n− 12 + i 2εk0 δx ψ 1 2 2

[( 1+

1 ≤ j ≤ M − 1, 1 ≤ n ≤ K,

N ∑

bm a m=1 m

n− 1

1

n− iεδt ψ˜M 2

n− 12

n− 1 am ) k02 ϕ˜m,12

−

n− 1 V0 2 am ϕ˜m,1

n− 1

n− 12

(1 − am ) k02 ϕ˜m,22 − VM ψ˜j0 = 0,

1

n− ψ˜0 2 − k02

n− 1

n− 1

1

1

N ∑

(3.36) ]}

bm ˜n− 12 ϕm,1 a m=1 m

+ V0 2 ψ˜0 2 + p0 2 , 1 ≤ n ≤ K, (3.37) { [( ) ] } N N ∑ ∑ bm ˜n− 21 1 √ bm ˜n− 12 ε2 ˜n− 21 2 =− i 2εk0 1+ ψM − k 0 ϕm,2 − δx ψM − 1 2 h a a 2 m=1 m m=1 m + VM

(1 −

)

0 ≤ j ≤ M,

+

n− n− ψ˜M 2 + pM 2 ,

n− 1 iεam δt ϕ˜m,12 1

=

n− 1 ψ˜0 2

+

1 ≤ n ≤ K, (3.38)

n− 1 rm,12 ,

1

1

n− n− n− am ϕ˜m,2 + iεam δt ϕ˜m,22 = ψ˜M 2 + rm,22 ,

ϕ˜0m,1 = 0,

ϕ˜0m,2 = 0,

n− 1

n− 1

1 ≤ m ≤ N, (3.39) 1 ≤ m ≤ N, (3.40)

1 ≤ m ≤ N.

(3.41) n− 1

2 2 2 Multiplying (3.36) by hψ˜j , (3.37) by 12 hψ˜0 , (3.38) by 21 hψ˜M , adding them together, and taking the imaginary part of the results bring )( ( ( N

)

) ∑ √ 1 bm

˜n 2 ˜n−1 2 ˜n− 21 2 ˜n− 21 2

ψ − ψ ψ0 + ψM

≤ − 2k0 1 + 2τ a m=1 m { } N N n− 12 ∑ b n− 12 ∑ b √ 3 1 1 n− n− m m ϕ˜m,12 + ψ˜M ϕ˜m,22 + 2k0 Re ψ˜0 a a m=1 m m=1 m    M −1  1 1 n− 1 ∑ n− 21 n− 1 n− 12 n− 1  1 2 n− 1  ψ˜ + Im p0 2 + ψ˜j pj 2 + ψ˜M pM 2  . (3.42)  εh 2 0  2

j=1

n− 12

Multiplying (3.39) by ϕ˜m,1 , and taking the imaginary part of the result yields } ( ) { n− 12 n− 12 n− 1 εam ˜n 2 ˜n−1 2 n− 1 ϕm,1 − ϕm,1 = Im ϕ˜m,1 ψ˜0 2 + ϕ˜m,1 rm,12 . 2τ √ Multiplying the above equation by 2k03 bm /am and summing m from 1 to N produce √ } { ( N N 2 2 ) √ ∑ n− 1 n− 1 1 1 2εk03 ∑ bm 3 ˜n− 2 ϕ˜ 2 + ϕ˜ 2 rn− 2 . 2k Im ψ bm ϕ˜nm,1 − ϕ˜n−1 = m,1 m,1 0 0 m,1 m,1 2τ m=1 a m=1 m (3.43) Similarly, from (3.40), √ } { ( N N 1 ) √ ∑ n− 12 n− 1 2εk03 ∑ bm n− 12 ˜n− 2 ˜n 2 ˜n−1 2 3 2 ˜ ˜ Im ψM ϕm,2 + ϕm,2 rm,2 . bm ϕm,2 − ϕm,2 = 2k0 2τ m=1 a m=1 m (3.44)

11

Computation of the Schr¨ odinger equation on unbounded domain

Denote ( N

2 √ 2 2 ) ∑

F n = ψ˜n + 2εk03 bm ϕ˜nm,1 + ϕ˜nm,2 . m=1

Adding (3.42), (3.43) and (3.44) yields ( )( N ) ∑ √ ) bm 1 ( n ˜n− 12 2 ˜n− 12 2 n−1 F −F ≤ − 2k0 1 + ψ0 + ψM + B 2τ a m=1 m { } N n− 12 n− 1 n− 21 n− 1 √ 3∑ bm 2 2 ˜ ˜ + 2k0 Im ϕm,1 rm,1 + ϕm,2 rm,2 a m=1 m    M −1  h 1 n− 1 ∑ n− 12 n− 1 n− 12 n− 1  1 2 n− 1  ψ˜ + Im p0 2 + ψ˜j pj 2 + ψ˜M pM 2  ε 2 0  2 j=1 (3.45) with ( { } { }) N 1 1 1 1 √ 3∑ bm n− 12 ˜n− 2 n− 12 ˜n− 2 n− 21 ˜n− 2 n− 12 ˜n− 2 ˜ ˜ ˜ ˜ B ≡ 2k0 Re ψ0 ϕm,1 + ψM ϕm,2 + Im ψ0 ϕm,1 + ψM ϕm,2 . a m=1 m By |Re(Z)|+|Im(Z)| ≤ gives

√ 2 |Z|, the θ inequality and taking θ =

√ 2 2 k0

( 1+

∑N

bm m=1 am

)

( N ) N ( ( ) ) k6 ∑ ∑ b ˜n− 12 2 ˜n− 12 2 ˜n− 12 2 ˜n− 21 2 m 0 B ≤ θ ψ0 + ψM + · bm ϕm,1 + ϕm,2 θ m=1 a2m m=1 ( )( √ ( N N ) ) √ ∑ ∑ 2k0 bm ˜n− 12 2 ˜n− 12 2 ˜n− 12 2 ˜n− 12 2 5 ≤ + 2ck b ϕ + 1+ ψ ψ + ϕ 0 M m,2 . m m,1 0 2 a m=1 m=1 m (3.46) Moreover, one has { } N N ∑ ∑ √ 3 bm ˜n− 12 n− 12 bm ˜n− 12 n− 21 + ϕ ϕ 2k0 Im r r a m,1 m,1 a m,2 m,2 m=1 m m=1 m [ N ( ( )] N ) √2 ∑ √ 3 ∑ bm n− 12 2 n− 12 2 ˜n− 12 2 ˜n− 21 2 ≤ 2k0 bm ϕm,1 + ϕm,2 + rm,1 + rm,2 4 m=1 a2m m=1 (3.47) and   √ ( ) −1 1 ∑ n− 12 n− 1 n− 1 h 1 n− 21 n− 1 M 2k0 ˜n− 12 2 ˜n− 21 2 1 ˜ 2 n− 2  2 ˜  ψ˜0 p 2 + ≤ ψ p + ψ p ψ + ψ 0 M j 0 j ε 2 2 M M 2 j=1 1 +∥ψ˜n− 2 ∥2 +

√ 2 ( M −1 ) 2h h ∑ n− 12 2 n− 12 2 n− 12 2 + + p . p p 0 M 16ε2 k0 4ε2 j=1 j

(3.48)

12

X. YANG AND J. ZHANG

Inserting (3.46)-(3.48) into (3.45) and by (3.27)-(3.29) yield ) 1 ( n F − F n−1 2τ √ ( N )

∑

˜n− 12 2 1 + 2ck02 √ n− 1 2 n− 1 2 bm ϕ˜m,12 + ϕ˜m,22 ≤ ψ 2εk03

+ ε m=1 √ 2 ( √ 3 N ) ( ) M −1 2 2 2h h ∑ n− 12 2 2k0 ∑ bm n− 12 2 n− 12 2 n− 12 n− 12 + p0 + pM + 2 p + rm,1 + rm,2 16ε2 k0 4ε j=1 j 4 m=1 a2m √ ) 1 ( )2 1 + 2ck02 ( n F + F n−1 + c2 τ 2 + εh2 . ≤ 2ε 4 √ 2 When τ ≤ ε/[3(1 + 2ck0 )], it holds that [ √ 2 ] )2 3(1 + 2ck0 ) 3 ( τ F n−1 + c2 τ 2 + εh2 . Fn ≤ 1 + ε 4 The Gronwall inequality yields F ≤ exp{ n

3(1 +

√

} { c2 ε(τ 2 + εh2 )2 2ck02 )nτ 0 √ . }· F + ε 4(1 + 2ck02 )

Noticing that F 0 = 0 completes the proof. 4. Frozen Gaussian approximation. The conclusion of Theorem 3.3 suggests that, in the semiclassical regime ε ≪ 1, one needs τ, h ≪ ε to obtain accurate results, 1 since lim ϵ2 e ϵ = ∞. In this case, we resort to the semiclassical approximation. The ϵ→0

idea is, instead of computing ψ ε directly in (1.1), one uses asymptotic analysis to derive and solve the equations of the macroscopic quantities (e.g. amplitude and action functions) associated with the wave function. Specifically, we propose the frozen Gaussian approximation (FGA) with appropriate boundary conditions to compute (1.1)-(1.2) on unbounded domain. FGA was originally proposed by Herman and Kluk ([18]) in quantum chemistry, and was recently generalized to strictly hyperbolic systems in [32, 33, 34]. 4.1. Formulation. We first state the formulation of FGA for the global space in this subsection, followed by the discussion of the boundary strategy in the next. FGA approximates the solution to the Schr¨odinger equation (1.1) by ∫ 1 ε ψFGA (t, x) = a(t, q, p)eiΦ(t,x,y,q,p)/ε ψ0ε (y)dydpdq, (4.1) (2πε)3d/2 where the phase function Φ is given by i i Φ(t, x, y, q, p) = S(t, q, p) + |x − Q|2 + P · (x − Q) + |y − q|2 − p · (y − q). (4.2) 2 2 The essence of FGA is to use Gaussian functions with fixed widths, instead of using those that might spread over time (e.g. [21, 29]), to approximate the solution. Specifically, as one can see from (4.1)-(4.2), the solution to the Schr¨odinger equation is approximated by a superposition of Gaussian functions living in the phase space.

Computation of the Schr¨ odinger equation on unbounded domain

13

Compared to the Gaussian beam method (GBM), each Gaussian function is not necessarily an asymptotic solution (e.g. [32, 33, 34]), and hence the method relies more on the superposition property of the Schr¨odinger equation than GBM. In (4.2), the i’s in front of |x − Q|2 and |y − q|2 make the function eiΦ/ε have a Gaussian profile with centers at Q and q in x and y spaces respectively. Q(t, q, p) and P (t, q, p) are defined as the center and momentum functions of each Gaussian function, which satisfy the Hamiltonian flow  dQ   = P,  dt (4.3)    dP = −∂Q V (Q, t), dt with the initial conditions Q(0, q, p) = q,

and P (0, q, p) = p.

(4.4)

There are also the classic trajectories given by the Hamiltonian flow H(t, Q, P ) = |P |2 /2 + V (Q, t). The action function S(t, q, p) satisfies dS |P |2 = − V (Q, t), dt 2

(4.5)

with the initial condition S(0, q, p) = 0. The right-hand side of (4.5) is the Lagrangian associated with the Hamiltonian H(t, Q, P ). The evolution of the amplitude function a(t, q, p) satisfies ( ) ( ) 1 da −1 2 = a tr Z ∂z P − i∂z Q∂Q V , (4.6) dt 2 with the initial condition a(0, q, p) = 2d/2 . In (4.6), we have introduced the short hand notations ∂z = ∂q − i∂p ,

Z = ∂z (Q + iP ).

(4.7)

Note that the amplitude function a(t, q, p) describes the weight of each Gaussian function in the summation integral (4.1). Remark 4.1. We refer interested readers to [18, 25] for the details of the derivation of (4.3)-(4.6) for the Schr¨ odinger equation (1.1). 4.2. Boundary strategy. In this subsection, we describe the strategy of implementing absorbing boundary conditions on FGA in details. Firstly, we rewrite (4.1) as ∫ ( ) i 2 1 1 ε ε ε S(t,q,p)+P ·(x−Q) e− 2ε |x−Q| dpdq, a(t, q, p)w (q, p)e ψFGA (t, x) = (2πε)3d/2 (4.8) where the weight function wε (q, p) is defined as ) ( ∫ 1 i 2 ε (4.9) w (q, p) = exp − p · (y − q) − |y − q| ψ0ε (y)dy. ε 2ε

14

X. YANG AND J. ZHANG

Secondly, we observe that (4.8) approximates the solution by an√integral (on q and p) of Gaussian functions centered at x = Q with the width of O( ε). Therefore ε ψFGA (t, x) only depends on the √ value of Gaussian functions whose centers are in the neighborhood of x of size O( ε). This implies that the waves going out of the domain are approximated by these Gaussian functions whose centers are at least √ O( ε) away from the boundary of domain. So the absorbing boundary conditions can be understood as the elimination of √ such Gaussian functions that are away from the boundary at a distance of at least O( ε). An intuitive explanation of the strategy is given in Fig. 4.1.

Q(t,q,p) Particle Trajectory

Initial Position q

Ω Ωi

ε i

Fig. 4.1. A cartoon illustration of the strategy of implementing the artificial boundary conditions for FGA: The particle trajectory starting at q which the center of a Gaussian function Q(t, q, p) follows will be ceased when it goes out the domain Ωεi , i.e. when Q(t, q, p) ∈ / Ωεi

Now we implement the absorbing boundary conditions: Denote Ωi as the domain of interest, and Ωεi ⊃ Ωi such that √ for any x ∈ ∂Ωεi , dist(x, ∂Ωi ) ∼ O( ε). The centers of Gaussian functions Q follow the Hamiltonian flow (4.3), and when Q∈ / Ωεi meaning the Gaussian centered at Q has no contribution to the solution in Ωi , we stop computing the dynamics of Q and ignore the Gaussian function centered at Q. 4.3. Algorithm. We first give a description of the overall algorithm. To construct the frozen Gaussian approximation on a mesh of x, one needs to compute the integral (4.8) numerically with a mesh on (q, p). This will relate to the numerical computation of (4.9) with a mesh of y. Hence three different meshes are needed in the algorithm. To better illustrate the idea of constructing the meshes, we consider the following WKB initial condition for (1.1), ψ0ε (x) = A(x) exp(iS0 (x)/ε),

(4.10)

then the stationary phase approximation implies that wε in (4.9) is localized around the submanifold p = ∇q S0 (q) on the phase plane when ε is small. This means we only need to put the mesh grids of p around ∇q S0 (q) initially to get a good approximation of (4.10). A one-dimensional example is given to illustrate this localization property of wε in Figure 4.2 (left). The associated mesh grids are shown in Figure 4.2 (right). Next we describe in details all the meshes used in the algorithm.

15

Computation of the Schr¨ odinger equation on unbounded domain 1 0.2

0.5

0.5

0

p

p

0.15

0

0.1 −0.5 0.05 −1 0

0.5

1

−0.5

−1 0

1.5

0.5

1

1.5

2

q

q

ε ε ) Left: an illustration of the localization of w on (q, p) domain for ψ0 (y) = (Fig. 4.2. sin(6y) , ε = 1/128; the black solid curve is p = cos(6q)/2. Right: the corresponding mesh exp i 12ε

grids of (q, p).

1. Discrete mesh of (q, p) for initializing Q, P . Denote δq = (δq1 , · · · , δqd ) and δp = (δp1 , · · · , δpd ) as the mesh size. Suppose q 0 = (q10 , · · · , qd0 ) is the starting point, then the mesh grids q k , k = (k1 , · · · , kd ), are defined as ( ) q k = q10 + (k1 − 1)δq1 , · · · , qd0 + (kd − 1)δqd , where kj = 1, · · · , Nq for each j ∈ {1, · · · , d}. The mesh grids pk,ℓ , ℓ = (ℓ1 , · · · , ℓd ), are defined associated with the mesh grids q k , ( ) pk,ℓ = ∂q1 S0 (q k ) + ℓ1 δp1 , · · · , ∂qd S0 (q k ) + ℓd δpd , where ℓj = −Np , · · · , Np for each j ∈ {1, · · · , d}. 2. Discrete mesh of y for evaluating wε in (4.9). δy = (δy1 , · · · , δyd ) is the mesh size. Denote y 0 = (y10 , · · · , yd0 ) as the starting point. The mesh grids y m are, m = (m1 , · · · , md ), ( ) y m = y10 + (m1 − 1)δy1 , · · · , yd0 + (md − 1)δyd , where mj = 1, · · · , Ny for each j ∈ {1, · · · , d}. 3. Discrete mesh of x for reconstructing the final solution. δx = (δx1 , · · · , δxd ) is the mesh size. Denote x0 = (x01 , · · · , x0d ) as the starting point. The mesh grids xn are, n = (n1 , · · · , nd ), ) ( xn = x01 + (n1 − 1)δx1 , · · · , x0d + (nd − 1)δxd , where nj = 1, · · · , Nx for each j ∈ {1, · · · , d}. With the preparation of the meshes, we introduce the algorithm as follows. Step 1. Compute the weight function wε by (4.9) for (Q, P ) initialized at (q k , pk,ℓ ). wε (q k , pk,ℓ ) =

∑

i

e ε (−p

k,ℓ

·(y m −q k )+ 2i |y m −q k |2 )

ψ0 (y m )rθ (|y m −q k |)δy1 · · · δyd ,

m

(4.11) where rθ is a cutoff function such that rθ = 1 in the ball of radius θ > 0 centered at origin and rθ = 0 outside the ball.

16

X. YANG AND J. ZHANG

Step 2. Solve the evolution equations (4.3) for (Q, P ), (4.5) for the action function S k,ℓ , and (4.6) for the amplitude function a by standard numerical integrator for ODE, for example the fourth-order Runge-Kutta scheme. Denote the numerical solutions as (Qk,ℓ , P k,ℓ ) and ak,ℓ . / Ωεi , Step 3. At each time step, check the position of Gaussian function Qk,ℓ . If Qk,ℓ ∈ stop computing the dynamics of Q and eliminate the corresponding index (k, ℓ). Step 4. Reconstruct the solution on the mesh of x by summing up the values of √ Gaussian functions whose centers Q are within an O( ε) distance of x, ) ( ) ∑( ak,ℓ rθ 1 |xn −Qk,ℓ |2 ε ε k k,ℓ εi S k,ℓ +P k,ℓ ·(xn −Qk,ℓ ) − 2ε ψFGA (t, x) = w (q , p )e . (2πε)3d/2 k,ℓ (4.12) Remark 4.2. 1. In setting up the meshes, we assume that the initial condition (4.10) either has compact support or decays sufficiently fast to zero as x → ∞ so that we only need finite number of mesh points in physical space. 2. The role of the truncation function rθ is to save computational cost, since although a Gaussian function √ is not localized, it decays quickly away from the center. In practice we take θ = O( ε), the same order as the width of each Gaussian, when we evaluate (4.11) and (4.12) numerically. 3. There are two types of errors present in the method. The first type comes from the asymptotic approximation to the wave equation ψ ε , which is O(ε), with the detailed analysis given in [41, 33]. This error can not be reduced unless one includes higher order corrections. The other type is the numerical error which comes from two sources: one is from the ODE numerical integrator; the other is from the discrete approximations of integrals (4.8) and (4.9). It can be reduced by either taking small mesh size and time step or using higher order numerical methods. 4. The boundary strategy of FGA works for any dimension, while the LABCs introduced in Section 2 has difficulties in being generalized to high dimensional cases when the boundary of Ωi has singularities (e.g. Ωi is a square.) 5. Step 2 and 4 can be expedited by making use of the discrete fast Gaussian transform, as in [38, 39]. 5. Numerical examples. Two numerical examples are presented to demonstrate not only the validity of the proposed LABCs, but also the accuracy of FGA with the boundary strategy. Example 5.1 focuses on the effectiveness of LABCs and investigate the dependence on the computational domain length, on the values of different ε, and on the spatial mesh size. When ε is very small (e.g. ε = 1/2048), we need to use tough mesh size to resolute the solutions. Finally, FGA method is given to demonstrate its effectiveness for small ε = 1/2048. Examples 5.2 only shows that FGA algorithm can be easily extended to two-dimension. In this example, we do not give the simulation by artificial boundary method since it is firstly hard to design suitable high-order boundary conditions at corners if the square computational domain is used, secondly, the computation cost will be expensive if low-order boundary conditions are implemented for small ε. In fact, the requirement of the spatial and temporal mesh sizes is quite tough in the semiclassical regime when the traditional methods are used. For example, a mesh size of O(ϵ) is required when using the time-splitting spectral method ([5]) to compute (1.1)-(1.2) directly; an even worse mesh size of o(ε) is needed if one uses the Crank-Nicolson schemes ([35]) or the Dufort-Frankel scheme

17

Computation of the Schr¨ odinger equation on unbounded domain Table 5.1 Example 5.1, the dependence of errors on the computational domain lengths.

[xl , xr ] E(h, τ )

[−2, 2] 6.949e-4

[−2.5, 2.5] 3.941e-5

[−3, 3] 9.505.0e-7

[−3.5, 3.5] 9.505e-7

[−4, 4] 9.505e-7

Table 5.2 Example 5.1, the errors corresponding to different ε’s for domain [-2,2] and τ = 10−5 .

ε E(h, τ )

1/32 4.648e-5

1/64 2.301e-5

1/128 2.047e-5

1/256 1.300e-4

1/512 1.027e-3

1/1024 7.723e-3

1/2048 2.730e-1

([36]). Example 5.1. (1d Schr¨ odinger) We choose the initial condition and potential function as ψ0 (x) = exp(−25x2 + i(x2 + 3x)/ε),

V (x) = exp(−0.2x2 ).

(5.1)

We shall fix N = 40 and k0 = 10 in the investigation of LABCs, including the dependence of errors on the domain length, the value of ε, and the spatial mesh size. We will also test the performance of FGA when ε ≪ 1. The dependence of error on domain length. In the calculation, we use the fixed spatial mesh size h = 1.0e − 4, and the time step τ = 1.0e − 5, ε = 1/512, and the final time T = 1. We compare the numerical solution ψnum to the reference solution ψex (x, t). The reference solution ψex (x, t) is obtained by computing the numerical solution in the large interval [−40, 40] with the time-splitting spectral method ([5]). For the different lengths of the computational domains, we list the errors between the reference solution and numerical solution in L2 -norm in Table 5.1. The error in L2 -norm is defined by E(h, τ ) = ∥ψexact − ψnum ∥ and ∥ · ∥ is given by Eq. (3.1). One can see that large computational interval gives small error in L2 -norm. This is actually due to the fact that, the larger the length of the computation domain is, the smaller the wave is reflected since the wave spends more time to arrive at the artificial boundary and crosses it. We remark that the error is calculated on the interval [−2, 2] only. The dependence of error on ε. In the calculation, we fixed h = 2.5e − 5, τ = 1.0e − 5, T = 1, [xl , xr ] = [−2, 2]. For the different choices of ε, the errors are showed in Table 5.2. One can see that the smaller the parameter ε is, the larger the L2 -error is. The error is caused by two factors: One is the reflected wave from boundary condition; the other is the small value of ε. Since the wave function becomes oscillatory of the wave length O(ε) when ε ≪ 1, we need to refine the mesh size to resolve the numerical solution. Fig. 5.1 shows the errors for the real and imaginary parts of the numerical ¯ e = [−2, 2] and changing ε solutions by fixing M = 1.6e5, K = 1.0e5, T = 1 and Ω from 1/256 to 1/2048. It also shows that smaller parameter ε gives larger error. The dependence of error on spatial mesh size. We compute the error by refining the mesh grid size and fixing ε = 1/512, τ = 1.0e − 5, T = 1, [xl , xr ] = [−2, 2]. From Table 5.3, one can see that the error becomes smaller when the spacial mesh size is refined, and the convergence rate is about 1.5. To further test the effectiveness of

18

X. YANG AND J. ZHANG Table 5.3 Example 5.1, the errors for different mesh sizes.

h E(h, τ ) order

1.0e-4 1.104e-2 –

5.0e-5 3.098e-3 1.746

2.5e-5 1.027e-3 1.593

Table 5.4 Example 5.1, the errors corresponding to different ε’s for domain [-1,1] and τ = 1.0e − 5.

ε E(h=5e-5,τ ) E(h=2.5e-5,τ ) E(h=1.25e-5,τ ) E(h=6.25e-6,τ )

1/128 1.803e-5 1.276e-5 1.223e-5 1.219e-5

1/256 5.492e-5 1.818e-5 1.112e-5 1.012e-5

1/512 2.495e-4 7.742e-5 3.824e-5 3.042e-5

1/1024 1.321e-3 4.493e-4 2.231e-4 1.708e-4

LABCs, we use a smaller computational interval [−1, 1] with the errors given in Table 5.4. In this case, the wave passes through the boundary quickly, hence it is reasonable to take smaller final time T = 0.7. In the computation, we fix τ = 1.0e − 5. From Table 5.4, the error becomes smaller when the spatial mesh size is refined for fixed ε. This implies that we can capture the accuracy of numerical solution with smaller mesh sizes for small ε. The error becomes larger when ε gets smaller in Fig. 5.2, which shows the errors of the real and imaginary parts of the numerical solutions by ¯ e = [−1, 1] and changing ε from 1/256 fixing M = 1.6e5, K = 1.0e5, T = 0.7 and Ω to 1/2048. Finally, to show the efficiency of our high-order ABCs, the evolution of the numerical solution and reference solution are plotted in Fig. 5.3 using parameters ¯ e = [−1, 1]. One can see that the main M = 8e4, K = 7e4, T = 0.7, ε = 1/512 and Ω part of wave is absorbed by ABCs and there is no obvious reflection. The performance of FGA for ε = 1/2048. We take δy = δq = δp = 1/210 , and the time step is δt = 1/211 . Note that wε in (4.9) has exponential decay in p-direction, hence one does not need to put mesh grids of p in the regions where wε is negligible; see Fig. 2.1 in [32] for example. The comparison of the real and imaginary parts of the FGA solution ψF GA with the reference solution ψex is shown in Fig. 5.4. Compared to the results in Figs. 5.1 and 5.2, which are obtained by high-order LABCs method, one can see that FGA gives better numerical results in the semiclassical regime ε ≪ 1. Example 5.2. (2d Schr¨ odinger) iε

∂ψ ε ε2 ω2 2 ε = − ∆ψ ε + |x| ψ , ∂t 2 2

x ∈ R2 ,

and the initial condition (1.2) as ) ( |x − x0 |2 , ψ0ε (x) = exp − 2ε where x = (x1 , x2 ) and x0 = (0.5, 0.5). In this example, we investigate the performance of FGA in two dimensions. The

19

Computation of the Schr¨ odinger equation on unbounded domain

%

$

ex

−ψ )|)

0

−10

log10(|Im( ψ

log10(|Re( ψ

−5

num

−5

num

ex

−ψ )|)

0

ε = 1/256 ε = 1/512 ε = 1/1024 ε = 1/2048

−15 −20 −2

−1

X 0

1

−10

−15 −20 −2

2

−1

X0

1

2

¯ e = [−2, 2], we change ε Fig. 5.1. Example 5.1, fixing M = 1.6e5, K = 1.0e5, T = 1 and Ω from 1/256 to 1/2048 to show the absolute errors of the real (in panel A) and imaginary (in panel B) parts between the reference solutions and numerical solutions in log10.

%

$

−6 ε = 1/256 ε = 1/512 ε = 1/1024 ε = 1/2048

−8 −10 −1

−0.5

0

ex

−2

num

−ψ )|)

0

−4

log10(|Im( ψ

ex

−4

log10(|Re( ψ

−2

num

−ψ )|)

0

0.5

1

−6 −8 −10 −1

−0.5

0

0.5

¯ e = [−1, 1], we change Fig. 5.2. Example 5.1, fixing M = 1.6e5, K = 1.0e5, T = 0.7 and Ω ε from 1/256 to 1/2048 to compare the absolute errors of the real (in panel A) and imaginary (in panel B) parts between the reference solutions and numerical solutions in log10.

exact solution of Schr¨ odinger equation is given by ( )−1 i ψ ε (x, t) = cos(ωt) + sin(ωt) ω { ( 1 i − ω sin(2ωt)|x0 |2 − ω sin(ωt)x0 · (x − x0 cos(ωt)) × exp ε 4 )} −ω sin(ωt) + i cos ωt 2 ) . + ( |x − x cos(ωt)| 0 2 cos(ωt) + ωi sin ωt In this test, we simulate the numerical solution in a relatively short time such that the wave is out-going and reach the artificial boundaries. We take the artificial domain Ωi to be the square [0, 1] × [0, 1]. We choose ε = 1/128, ω = 1/4, and the mesh size δq1 = δq2 = δp1 = δp2 = δy1 = δy2 = 1/32. The numerical results are shown in Fig. 5.5, from which one can see the solution by FGA agrees with the reference

1

20

X. YANG AND J. ZHANG

%

A$

1

1

|ψnu |

|ψ | ex

0 0.7

0 0.7 1

0.35 Time

1

0.35 0 −1

0 X

Time

0 0 −1

X

¯ e = [−1, 1], ε = 1/512. The Fig. 5.3. Example 5.1, fixing M = 8e4, K = 1.0e5, T = 0.7 and Ω numerical solutions (in panel A) and the reference solutions (in panel B) are evolved in successive times versus x.

solution well, and the error is in the order of O(ε). 6. Conclusion. In this paper, we studied the computation of the Schr¨odinger equation on unbounded domain in the semiclassical regime. First, we generalized the high-order local absorbing boundary conditions (LABCs) proposed in [48] to compute the Schr¨odinger equation, and analyzed the stability of the equation with LABCs and the convergence of the corresponding Crank-Nicolson scheme. The conclusion shows that when the rescaled Planck constant ε gets small, the accuracy deteriorates and the requirements on time step and mesh size become restricted. Then we propose the asymptotic method based on frozen Gaussian approximation for this case of ε ≪ 1, with the absorbing boundary condition dealt by a simple strategy that all the effects of the Gaussian functions which contribute to the outgoing waves will be eliminated by stopping the Hamiltonian flow of their centers when they get out of the domain of interest. We present numerical examples in both one and two dimensions to verify the performance of the proposed numerical methods. The extension of the method proposed in this paper to nonlinear Schr¨ odinger equation in the semiclassical regime on unbounded domain will be discussed in the future. Acknowledgement. X.Y. would like to thank Prof. Xiaonan Wu for his valuable discussions. REFERENCES [1] X. Antoine, C. Besse. Unconditional stable discretization schemes of non-reflected boundary conditions for the one-dimensional Schr¨ odinger equation. J. Comput. Phys., 188 (2003) 157–175. [2] X. Antoine, C. Besse, V. Mouysset. Numerical schemes for simulation of the two-dimensional Schr¨ odinger equation using non-reflected boundary conditions. Math. Comput., 73 (2004) 1779–1799. adle, A review of transparent and [3] X. Antoine, A. Arnold, C. Besse, M. Ehrhardt, and A. Sch¨ artificial boundary conditions techniques for linear and nonlinear Schr¨ odinger equations, Commun. Comput. Phys., 4 (4) (2008), 729–796. odinger equation with [4] X. Antoine, C. Besse, P. Klein, Absorbing boundary conditions for Schr¨ an exterior repulsive potential, J. Comput. Phys., 228 (2009), 312–335.

21

Computation of the Schr¨ odinger equation on unbounded domain

ε=1/2048

−3

Re(ΨFGA−Ψex)

5

x 10

0 −5 −2

−1

0 x

1

2

−1

0 x

1

2

−3

Im(ΨFGA−Ψex)

5

x 10

0

−5 −2

ε=1/2048

−3

Re(ΨFGA−Ψex)

2

x 10

0 −2 −1

−0.5

0 x

0.5

1

−0.5

0 x

0.5

1

−3

Im(ΨFGA−Ψex)

2

x 10

0

−2 −1

Fig. 5.4. Example 5.1, the comparison of the FGA solution to the reference solution. Left: Final time t = 1 and domain [−2, 2]; right: Final time t = 0.7 and domain [−1, 1]. Re(ΨFGA−Ψex) at T= 6

Im(ΨFGA−Ψex) at T= 6

−3

x 10

0.9

0.9 2

0.8

1

1

0

2

0.5

0

0.6 x

0.6 2

0.8 0.7

0.7

x

−3

x 10

0.5 −1

0.4

0.4 −1

0.3 0.2

−2

0.3 −2

0.2 0.1

0.1 0.2

0.4

0.6 x1

0.8

0.2

0.4

0.6

0.8

−3

x1

Fig. 5.5. Example 5.2, the comparisons of the errors of the real (left figure) and imagine (right figure) parts at the final time T = 6 between the reference solution and the numerical solution by FGA. The parameters are ε = 1/128, ω = 1/4, ∆q1 = ∆q2 = ∆p1 = ∆p2 = ∆y1 = ∆y2 = 1/32.

22

X. YANG AND J. ZHANG

[5] W. Bao, S. Jin and P.A. Markowich, On time-splitting spectral approximations for the Schr¨ odinger equation in the semiclassical regime, J. Comput. Phys. 175 (2002), 2, 487– 524. [6] E. B´ ecache, D. Givoli, T. Hagstrom. High-order absorbing boundary conditions for anisotropic and convective wave equations, J. Comput. Phys. 229 (2010), 1099–1129. odinger equation [7] L. Di Menza, Absorbing boundary conditions on a hypersurface for the Schr¨ in a half space, Appl. Math. Lett., 9 (1996), 55–59. [8] B. Enguist, A. Majda, Absorbing boundary conditions for the numerical simulation of waves, Math. Comp. 31 (1977), 629–651. [9] B. Engquist and A. Majda, Radiation boundary conditions for acoustic and elastic wave calculations, Comm. Pure Appl. Math., 32 (1979) 314–358. [10] B. Engquist and O. Runborg, Computational high frequency wave propagation, Acta Numer., 12 (2003), 181–266. [11] T. Fevens and H. Jiang, Absorbing boundary conditions for the Schr¨ odinger equation, SIAM J. Sci. Comput., 21(1) (1999), 255–282. [12] G. Hagedorn, Semiclassical quantum mechanics I: ~ → 0 limit for coherent states. Commun. Math. Phys., 71 (1980), 77–93. [13] T. Hagstrom, A. Mar-Or, D. Givoli, High-order local absorbing conditions for the wave equation: extensions and improvements, J. Comput. Phys., 227 (2008), 3322–3357. [14] H. Han, The artificial boundary method - numerical solutions of partial differential equations in unbounded domains, Frontiers and Prospents of Contemporary Applied Mathematics, Higher Education Press, World Scientific, 2005. [15] H. Han and X. Wu, Artificial boundary method, Springer-Verlag, 2013. [16] H. Han, Z. Huang, Exact artificial boundary conditions for the Schr¨ odinger equation in R2 , Comm. Math. Sci., 2(1) (2004), 79–94. [17] H. Han,D. Yin, Z. Huang, Numerical solution of Schr¨ odinger equation in R3 , Numer. Methods Partial Differ. Equ., 23 (2007), 511–533. [18] M. Herman and E. Kluk, A semiclassical justification for the use of non-spreading wavepackets in dynamics calculations, Chem. Phys., 91 (1984), 27–34. [19] S. Jin, P.A. Markowich and C. Sparber, Mathematical and computational methods for semiclassical Schr¨ odinger equations, Acta Numer., 20 (2011), 211–289. [20] S. Jin, H. Wu and Z. Huang, A Hybrid Phase-Flow Method for Hamiltonian Systems with Discontinuous Hamiltonians, SIAM J. Sci. Comput., 31 (2008), 1303–1321. [21] S. Jin, H. Wu and X. Yang, Gaussian beam methods for the Schr¨ odinger equation in the semiclassical regime: Lagrangian and Eulerian Formulations, Commun. Math. Sci., 6 (2008), 995–1020. [22] S. Jin, H. Wu, X. Yang, and Z. Huang, Bloch decomposition-based Gaussian beam method for the Schr¨ odinger equation with periodic potentials, J. Comput. Phys., 229 (2010), 4869– 4883. [23] S. Jin, H. Wu, and X. Yang, A numerical study of the Gaussian beam methods for onedimensional Schr¨ odinger-Poisson equations, J. Comput. Math., 28 (2010), 261–272. [24] S. Jin, H. Wu, and X. Yang, Semi-Eulerian and high order Gaussian beam methods for the Schr¨ odinger equation in the semiclassical regime, Commun. Comput. Phys., 9 (2011), 668– 687. [25] K. Kay, The Herman-Kluk approximation: Derivation and semiclassical corrections, Chem. Phys., 322 (2006), 3–12. [26] J. Kuska, Absorbing boundary conditions for the Schr¨ odinger equation on finite intervals, Phys. Rev. B, 46(8) (1999), 5000–5003. [27] E. Lindmann, Free space boundary conditions for the time dependent wave equation, J. Comput. Phys., 18 (1985), 16–78. [28] S. Leung, J. Qian and R. Burridge, Eulerian Gaussian beams for high-frequency wave propagation, Geophysics, 72 (2007), 61–76. [29] S. Leung and J. Qian, Eulerian Gaussian beams for Schr¨ odinger equations in the semiclassical regime, J. Comput. Phys., 228 (2009), 2951–2977. [30] H. Liu and J. Ralston, Recovery of high frequency wave fields for the acoustic wave equation, Multiscale Model. Simul., 8 (2009), 428–444. [31] H. Liu and J. Ralston, Recovery of high frequency wave fields from phase space based measurements, Multiscale Model. Simul., 8 (2010), 622–644. [32] J. Lu and X. Yang, Frozen Gaussian approximation for high frequency wave propagation, Commun. Math. Sci., 9 (2011), 663–683. [33] J. Lu and X. Yang, Convergence of frozen Gaussian approximation for high frequency wave propagation, Comm. Pure Appl. Math., 65 (2012), 759–789.

Computation of the Schr¨ odinger equation on unbounded domain

23

[34] J. Lu and X. Yang, Frozen Gaussian approximation for general linear strictly hyperbolic system: Formulation and Eulerian methods, Multiscale Model. Simul., 10 (2012), 451–472. [35] P. Markowich, P. Pietra and C. Pohl, Numerical approximation of quadratic observable of Schr¨ odinger equation-type equations in the semiclassical limit. Numer. Math., 81(1999), 595-630. [36] P. Markowich, P. Pietra, C. Pohl and H. Stimming, A wigner-measure analysis of the DufortFrankel scheme for the Schr¨ odinger equation, SIAM J. Numer. Anal., 40 (2002),1281–1310. [37] M. Motamed and O. Runborg, Taylor expansion and discretization errors in Gaussian beam superposition, Wave Motion, 47 (2010), 421–439. [38] J. Qian and L. Ying, Fast Gaussian wavepacket transforms and Gaussian beams for the Schr¨ odinger equation, J. Comput. Phys., 229 (2010), 7848–7873. [39] J. Qian and L. Ying, Fast multiscale Gaussian wavepacket transforms and multiscale Gaussian beams for the wave equation, Multiscale Model. Simul., 8 (2010), 1803–1837. [40] J. Ralston, Gaussian beams and the propagation of singularities, Studies in PDEs, MAA Stud. Math., 23 (1982), 206–248. [41] T. Swart and V. Rousse, A mathematical justification of the Herman-Kluk propagator, Commun. Math. Phys., 286 (2009), 725–750. [42] J. Szeftel, Design of absorbing boundary conditions for Schr¨ odinger equations in Rd , SIAM J. Numer. Anal., 42(4) (2004), 1527–1551. [43] N. Tanushev, J. Qian, and J. Ralston, Mountain waves and Gaussian beams, Multiscale Model. Simul., 6 (2007), 688–709. [44] N. Tanushev, Superpositions and higher order Gaussian beams, Commun. Math. Sci., 6 (2008), 449–475. [45] S.V. Tsynkov, Numerical solution of problems on unbounded domain: a review. Appl. Numer. Math., 27 (1998), 465–532. [46] Z. Xu, H. Han, X. Wu, Adaptive absorbing boundary conditions for Schr¨ odinger-type equations: application to nonlinear and multi-dimensional problems, J. Comput. Phys. 225 (2007) 1577–1589. [47] J. Zhang, Z. Xu, X. Wu and D. Wang. Stability analysis of local absorbing boundary conditions for general nonlinear Schr¨ odinger equations in one and two dimensions. To appear. [48] J. Zhang, Z. Sun, X. Wu and D. Wang, Analysis of high-order absorbing boundary conditions for the Schr¨ odinger equation. Commun. Comput. Phys. 10(2011), 742–766.