Efficient Frozen Gaussian Sampling Algorithms for Nonadiabatic Quantum Dynamics at Metal Surfaces

Zhen Huang [email protected] Limin Xu [email protected] Zhennan Zhou [email protected] Department of Mathematics, University of California, Berkeley, California 94720 USA Department of Mathematical Sciences, Tsinghua University, Beijing 100084, China Beijing International Center for Mathematical Research, Peking University, Beijing, 100871, China

Abstract

In this article, we propose a Frozen Gaussian Sampling (FGS) algorithm for simulating nonadiabatic quantum dynamics at metal surfaces with a continuous spectrum. This method consists of a Monte-Carlo algorithm for sampling the initial wave packets on the phase space and a surface-hopping type stochastic time propagation scheme for the wave packets. We prove that to reach a certain accuracy threshold, the sample size required is independent of both the semiclassical parameter $\varepsilon$ and the number of metal orbitals $N$ , which makes it one of the most promising methods to study the nonadiabatic dynamics. The algorithm and its convergence properties are also validated numerically. Furthermore, we carry out numerical experiments including exploring the nuclei dynamics, electron transfer and finite-temperature effects, and demonstrate that our method captures the physics which can not be captured by classical surface hopping trajectories.

keywords:

metal surfaces , Frozen Gaussian Sampling , nonadiabatic quantum dynamics , Semiclassical Schrödinger equation system

1 Introduction

Quantum dynamics of atomic systems are well known to exhibit multiscale structures, where the nuclei are much heavier than the electrons and therefore move much slower. The widely-used Born-Oppenheimer (BO) approximation [1] is proposed based on such a scale-separation structure, where one regards the nuclei as classical particles, and ignore the couplings between different electronic states. The latter treatment is also termed the adiabatic assumption [2].

On the one hand, the BO approximation of the Schrödinger equation, unfortunately, breaks down dramatically in many scenarios since the adiabatic assumptions are often violated [3, 4, 5]. On the other hand, the Schrödinger equation, as a high-dimensional model for quantum dynamics, can not be solved in a brute-force way because of the curse of dimensionality. Therefore, how to efficiently calculate dynamics in the non-adiabatic regime is a central topic in chemistry and material science. In order to tackle these issues, various methods based on mixed quantum-classical dynamics [6, 7, 8, 9, 10] have been developed to simulate non-adiabatic quantum dynamics.

Molecule-metal interfaces are of particular interest since it is related to many different phenomenon in experimental chemistry, such as chemisorption [11], electrochemistry [12], heterogeneous catalysis [13] and molecular junctions [14]. The breakdown of BO approximations at metal surfaces leads to theoretical study of many interesting physical phenomenon: such as electronic friction [15, 16], electron transfer [17] and energy transfer [18]. The well-known Newns-Anderson model [11], or the Anderson-Holstein model [19], is widely used to model nonadiabatic quantum dynamics at metal surfaces [16, 20]. In the first quantization, this model could be written as the following Schrödinger equation system:

\begin{array}[]{c}{\left\{\begin{aligned} \mathrm{i}\varepsilon\partial_{t}u_{0}(t,x)&=-\frac{\varepsilon^{2}}{2}\Delta u_{0}(t,x)+U_{0}(x)u_{0}(t,x)+\varepsilon\sum_{k=1}^{N}hV\left({\mathcal{E}_{k}},x\right)u_{k}(t,x),\\ \mathrm{i}\varepsilon\partial_{t}u_{k}(t,x)&=-\frac{\varepsilon^{2}}{2}\Delta u_{k}(t,x)+\left(U_{1}(x)+{\mathcal{E}}_{k}\right)u_{k}(t,x)+\varepsilon\overline{V}({\mathcal{E}_{k}},x)u_{0}(t,x).\end{aligned}\right.}\\ x\in\mathbb{R}^{m},\quad t\in\mathbb{R},\quad k=1,\cdots,N.\end{array}

(1.1)

with initial value

\left\{\begin{aligned} u_{0}(0,x)&=u_{0}^{\mathrm{in}}(x),\\ u_{k}\left(0,x\right)&=0,\qquad k=1,\cdots,N.\end{aligned}\right.

(1.2)

We will discuss the derivation and details of 1.1 in section 2. For now, let us emphasize its major difficulties as of numerical simulations. On the one hand, there is a non-dimensional parameter $\varepsilon$ presumed to be very small, therefore this model is exactly in the regime of semiclassical dynamics [21, 22]. The solution is highly oscillatory both in space and time [21]. On the other hand, we need to solve an $(N+1)\times(N+1)$ matrix Schrödinger equation system, where $N$ comes from discretizing the continuum band of metal, and $N$ needs to be large enough in order to capture the physics of continuum band emerging from condensed phases of matter. We need a reasonable algorithm whose cost doesn’t scale dramatically as $\varepsilon\rightarrow 0$ and $N\rightarrow\infty$ . Methods based on direct discretization of spatial freedom, such as time-splitting spectral methods [23], is not ideal since its cost scales (at least) $O(1/\varepsilon^{1+m})$ ( $m$ is the dimension of spatial coordinates) as $\varepsilon\rightarrow 0$ and (at least) $O(N)$ as $N\rightarrow\infty$ .

In scientific literature [20], there are generally two kinds of numerical schemes to calculate nonadiabatic quantum dynamics at metal surfaces: Erhenfest dynamics [24] (possibly with electronic friction [15]) and surface hopping methods [25]. With or without electronic friction, Erhenfest dynamics are in fact a mean-field approximation of the coupled electron–nuclear system, therefore is not valid when strong nonadiabatic dynamics occur [26]. Independent Electron Surface Hopping (IESH) [25], as the surface hopping method specifically designed for metal surfaces, is proposed within the same spirit of Tully’s original surface hopping method [27]. There are many variants of surface hopping methods in literature, most of which lack a rigorous mathematical foundation and numerical analysis. In [28], the authors provide a Frozen Gaussian Approximation with surface hopping (FGA-SH) method along with a rigorous mathematical justification for $O(\varepsilon)$ accuracy. However, its generalization to simulating metal surfaces demands proper treatment of a large number of bands with a particular coupling structure, thus such an investigation poses non-trivial challenges in terms of algorithm design and numerical analysis.

In this article, we propose a Frozen Gaussian Sampling (FGS) algorithm for simulating nonadiabatic quantum dynamics at metal surfaces. Based on a stochastic approximation of the Frozen Gaussian representation of the Schrödinger equation system, the wave packets to be sampled incorporate randomness from both the initial state and from the stochastic evolution. We rigorously prove and numerically validate that the sample size required by the FGS algorithm is independent from both the number of metal orbitals $N$ and the semiclassical parameter $\varepsilon$ . From this sense, FGS as an asymptotic method to the semiclassical Schrödinger equation system is superior to the direct application of other existing asymptotic schemes, such as the WKB methods [29], Landau-Zener transition asymptotic methods [30, 31] and various wave packet based asymptotic methods [32, 33, 34, 35]. A by-product is we are able to prove that the computational cost of FGS algorithm in finite-band systems [28] is also essentially independent of the semiclassical parameter $\varepsilon$ . What’s more, using numerical experiments, we demonstrate that the FGS method can capture information of observables that are intrinsically quantum, in the sense that such information can not be captured by classical surface hopping trajectories.

The article is organized as followed. In section 2, we give an introduction of the Anderson-Holstein model (Newns-Anderson model) and its nondimensionalization, along with its discretization on the continuous spectrum. In section 3, we present our algorithm for simulating nonadiabatic dynamics at metal surfaces. After reviewing Frozen Gaussian representation for one-level system in section 3.1, we derive the integral representation of the approximate solution based on the surface hopping ansatz in section 3.2. Based on its stochastic interpretation in section 3.3, we present the Frozen Gaussian Sampling algorithm in section 3.4. In section 4, we carry out the numerical analysis of our method, indicating that the cost of this stochastic method is independent of both $\varepsilon$ and $N$ (See Theorem 4.1). Finally, we provide extensive numerical experiments to illustrate the numerical advantages of our method and some explorations of this model in section 5.

2 Quantum dynamics of nuclei at metal surfaces

The quantum dynamics of nuclei at metal surfaces are described by the Newns-Anderson model [11], also known as the Anderson-Holstein model [19], which is a special kind of Anderson impurity model: the electronic ground state of the molecule is treated as the system orbital, while the metal electronic orbitals are treated as bath orbitals. The energies of metal electronic orbitals form a continuous spectrum $[\mathcal{E}_{a},\mathcal{E}_{b}]$ . The model Hamiltonian is usually written in the following second quantization form:

	$\displaystyle\hat{H}$	$\displaystyle=\frac{\hat{p}^{2}}{2m_{\text{n}}}+U_{1}(\hat{x})+h(\hat{x})\hat{d}^{\dagger}\hat{d}$		(2.1)
		$\displaystyle+\int_{\mathcal{E}_{a}}^{\mathcal{E}_{b}}(\mathcal{E}-\mu)\hat{c}_{\mathcal{E}}^{\dagger}\hat{c}_{\mathcal{E}}\mathrm{d}\mathcal{E}+\int_{\mathcal{E}_{a}}^{\mathcal{E}_{b}}\left(V(\mathcal{E},\hat{x})\hat{c}_{\mathcal{E}}^{\dagger}\hat{d}+\overline{V}(\mathcal{E},\hat{x})\hat{d}^{\dagger}\hat{c}_{\mathcal{E}}\right)\mathrm{d}\mathcal{E}.$		(2.1)

Here $\hat{p}$ is the momentum operator, $\hat{x}$ is the position operator, $m_{\text{n}}$ is the mass of the nuclei, $U_{1}(\hat{x})$ is the nuclear potential for the neutral molecule, $U_{1}(\hat{x})+h(\hat{x})$ is the nuclear potential for the charged molecule, $\mu$ is the chemical potential. And, $\hat{d}$ and $\hat{d}^{\dagger}$ are the annihilation and creation operators for the electronic ground state of the molecule, $\hat{c}_{\mathcal{E}}$ and $\hat{c}_{\mathcal{E}}^{\dagger}$ are the annihilation and creation operators for metal electronic orbitals with energy level $\mathcal{E}\in[\mathcal{E}_{a},\mathcal{E}_{b}]$ , $V\left(\mathcal{E},x\right)$ describes the coupling between the molecule and metal orbitals, and $\overline{V}\left(\mathcal{E},x\right)$ means the complex conjugate of ${V}\left(\mathcal{E},x\right)$ .

Now let us demonstrate the non-dimensionalization of this model. Let $\ell$ be the characteristic length scale and $E$ be the characteristic energy scale, we introduce the dimensionless parameter $\varepsilon$ , which is referred to as the semiclassical parameter:

\varepsilon=\frac{\hbar}{\ell\sqrt{mE}}.

(2.2)

Let $E_{V}$ be the characteristic scale of molecule-metal coupling energy, we are ready to perform the following rescaling (similar as in [36]):

\begin{array}[]{c}x=\ell\tilde{x},\quad t=\sqrt{\frac{m\ell^{2}}{E}}\tilde{t},\quad\mathcal{E}=E\tilde{\mathcal{E}},\quad\mu=E\tilde{\mu},\\ \tilde{U}_{1}(\tilde{x})=\frac{U_{1}(x)}{E},\quad\tilde{h}(\tilde{x})=\frac{h(x)}{E},\quad,\tilde{V}(\tilde{\mathcal{E}},\tilde{x})=\frac{V(\mathcal{E},x)}{E_{V}},\quad\hat{H}=E{{\hat{H}_{\text{non}}}}.\end{array}

(2.3)

And particularly, we choose the regime $E_{V}=\varepsilon E$ , and the corresponding non-dimensionalized Hamiltonian is

	$\displaystyle\hat{H}_{\text{non}}$	$\displaystyle=-\frac{\varepsilon^{2}}{2}\nabla_{\tilde{x}}^{2}+\tilde{U}_{1}(\tilde{x})+\tilde{h}(\tilde{x})\hat{d}^{\dagger}\hat{d}$		(2.4)
		$\displaystyle+\int_{\tilde{\mathcal{E}}_{a}}^{\tilde{\mathcal{E}}_{b}}(\tilde{\mathcal{E}}-\tilde{\mu})\hat{c}_{\tilde{\mathcal{E}}}^{\dagger}\hat{c}_{\tilde{\mathcal{E}}}\mathrm{d}\mathcal{E}+\varepsilon\int_{\tilde{\mathcal{E}}_{a}}^{\tilde{\mathcal{E}}_{b}}\left(\tilde{V}(\tilde{\mathcal{E}},\tilde{x})\hat{c}_{\mathcal{E}}^{\dagger}\hat{d}+\overline{V}(\mathcal{E},\hat{x})\hat{d}^{\dagger}\hat{c}_{\mathcal{E}}\right)\mathrm{d}\mathcal{E}.$		(2.4)

We can see that the coupling between the molecule and metal electronic orbitals is chosen to be $O(\varepsilon)$ . By the scientific literature [17, 16], this scenario can be intepreted as the weak-coupling regime of nonadiabatic dynamics at metal surfaces. However, even though the coupling is considered to be weak, since the Hamiltonian has a multiscale structure, such an $O(\varepsilon)$ coupling still causes an $O(1)$ change in the density of states. As a matter of fact, in the adiabatic representation, the coupling between different levels is also formally $O(\varepsilon)$ , see [28].

In this article, our starting point is to rewrite $\hat{H}_{\operatorname{non}}$ in the first quantization. Let $x\in\mathbb{R}^{m}$ represents the nuclei degrees of freedom, and let $U_{0}(x)=U_{1}(x)+h(x)$ denote the nuclear potential for the charged molecule, and recall that $U_{1}(x)$ represents the nuclear potential for the neutral molecule. The above model is described by the following Schrödinger equation system:

\left\{\begin{aligned} \mathrm{i}\varepsilon\partial_{t}\psi_{0}(t,x)&=-\frac{\varepsilon^{2}}{2}\Delta\psi_{0}(t,x)+U_{0}(x)\psi_{0}(t,x)+\varepsilon\int_{\mathcal{E}_{a}}^{{\mathcal{E}_{b}}}\mathrm{d}\tilde{\mathcal{E}}V(\tilde{\mathcal{E}},x)\psi_{1}(t,\tilde{\mathcal{E}},x),\\ \mathrm{i}\varepsilon\partial_{t}\psi_{1}(t,\mathcal{E},x)&=-\frac{\varepsilon^{2}}{2}\Delta\psi_{1}(t,\mathcal{E},x)+\left(U_{1}(x)+{\mathcal{E}}\right)\psi_{1}(t,\mathcal{E},x)+\varepsilon\overline{V}\left({\mathcal{E}},x\right)\psi_{0}(t,x).\end{aligned}\right.

(2.5)

We have dropped the “tilde” and without loss of generality we assume $\mu=0$ for simplification. Here $\psi_{0}(t,x)\in L^{2}(\mathbb{R}^{m})$ represents the nuclei wave function when the independent electron propagates into the molecule and therefore the molecule is charged, while $\psi_{1}\left(t,\mathcal{E},x\right)\in L^{2}\left([\mathcal{E}_{a},\mathcal{E}_{b}]\times\mathbb{R}^{m}\right)$ represents the nuclei wave function when the molecule is neutral and the electron is in the metal orbital with energy $\mathcal{E}$ . $\psi_{0}(t,x)$ and $\psi_{1}\left(t,\mathcal{E},x\right)$ satisfy the following normalization condition (also known as, mass conservation):

1=m(t)=\int_{\mathbb{R}^{m}}\left(|\psi_{0}(t,x)|^{2}\mathrm{d}x+\int_{\mathcal{E}_{a}}^{\mathcal{E}_{b}}|\psi_{1}\left(t,\mathcal{E},x\right)|^{2}\mathrm{d}\mathcal{E}\right)\mathrm{d}x,\quad\forall t\geq 0.

(2.6)

We remark that in this work we are considering a closed quantum system. This is exactly the viewpoint of Newn’s pioneering work [11]. In some scientific literature [37, 20], the setup of metal surfaces are open quantum systems, where one focus only on the molecule orbital (the orbital corresponding to $\psi_{0}$ ) and treat the metal orbitals (the orbitals corresponding to $\psi_{1}$ ) as heat baths with Boltzmann distribution. We would like to emphasize that even in the closed-system treatment, the nonadiabatic features of quantum dynamics have already arisen. In other words, one do not need the open quantum system and canonical ensemble assumption to encounter nonadiabatic quantum dynamics at metal surfaces. We leave the open quantum system treatment of metal surfaces to future work.

The energy $E(t)$ of the above system 2.5 is defined as:

$\displaystyle E(t)$	$\displaystyle=\int_{\mathbb{R}^{m}}\left(\frac{\varepsilon^{2}}{2}\left\|\nabla_{x}\psi_{0}(t,x)\right\|^{2}+U_{0}(x)\left\|\psi_{0}(t,x)\right\|^{2}\right)\mathrm{d}x$	(2.7)
	$\displaystyle+\int_{\mathbb{R}^{m}}\int_{\mathcal{E}_{a}}^{\mathcal{E}_{b}}\left(\frac{\varepsilon^{2}}{2}\left\|\nabla_{x}\psi_{1}\left(t,\mathcal{E},x\right)\right\|^{2}+\left(U_{1}(x)+\mathcal{E}\right)\left\|\psi_{1}\left(t,\mathcal{E},x\right)\right\|^{2}\right)\mathrm{d}\mathcal{E}\mathrm{d}x$
	$\displaystyle+\varepsilon\int_{\mathbb{R}^{m}}\int_{\mathcal{E}_{a}}^{\mathcal{E}_{b}}\left(\bar{\psi}_{0}(t,x)V\left(\mathcal{E},x\right)\psi_{1}\left(t,\mathcal{E},x\right)+\psi_{0}(t,x)\overline{V}\left(\mathcal{E},x\right)\bar{\psi}_{1}\left(t,\mathcal{E},x\right)\right)\mathrm{d}\mathcal{E}\mathrm{d}x.$

It’s standard to verify that 2.5 satisfies the energy conservation: $E(t)=E(0)$ for any $t\geq 0$ .

For simplicity, we assume at the initial state the molecule is charged, and thus the electron is in the adsorbate orbital. In other words, we consider the following initial condition:

\psi_{0}(0,x)=\psi_{0}^{\mathrm{in}}(x),\quad\psi_{1}\left(0,\mathcal{E},x\right)=0.

(2.8)

However, since 2.5 is a linear equation system, the results in this article are applicable to any initial conditions.

A typical treatment of 2.5 is to discretize the electronic continuum spectrum $[{\mathcal{E}_{a}},{\mathcal{E}_{b}}]$ into $N$ separate orbitals. Then we obtain the discrete Anderson-Holstein model, which is a commonly used model to investigate surface hopping at molecule-metal surfaces in scientific literature (for example, see [25, 37, 20]). For the sake of simplicity, we use the following equidistant discretization throughout this article:

\mathcal{E}_{k}=\mathcal{E}_{a}+(i-\frac{1}{2})h,\quad h=\frac{1}{N}(\mathcal{E}_{b}-\mathcal{E}_{a}),\quad k=1,\cdots,N.

(2.9)

Note that the generalization to any non-equidistant discretization is trivial. Let

u_{0}(t,x)=\psi_{0}(t,x),\quad u_{k}(t,x)=\psi_{1}(t,\mathcal{E}_{k},x),\quad k=1,\cdots,N.

Then we can write down the Anderson-Holstein model for $\{u_{k}(t,x)\}_{k=0}^{N}$ , which we have already presented in section 1:

\begin{array}[]{c}{\left\{\begin{aligned} \mathrm{i}\varepsilon\partial_{t}u_{0}(t,x)&=-\frac{\varepsilon^{2}}{2}\Delta u_{0}(t,x)+U_{0}(x)u_{0}(t,x)+\varepsilon\sum_{k=1}^{N}hV\left({\mathcal{E}_{k}},x\right)u_{k}(t,x),\\ \mathrm{i}\varepsilon\partial_{t}u_{k}(t,x)&=-\frac{\varepsilon^{2}}{2}\Delta u_{k}(t,x)+\left(U_{1}(x)+{\mathcal{E}}_{k}\right)u_{k}(t,x)+\varepsilon\overline{V}({\mathcal{E}_{k}},x)u_{0}(t,x).\end{aligned}\right.}\\ x\in\mathbb{R}^{m},\quad t\in\mathbb{R},\quad k=1,\cdots,N.\end{array}

with initial value

u_{0}(0,x)=u_{0}^{\mathrm{in}}(x),\quad u_{k}\left(0,x\right)=0,\quad k=1,\cdots,N.

Then for the systems of $\{u_{k}(t,x)\}_{k=0}^{N}$ , the total mass is defined as

m_{N}(t)=\int_{\mathbb{R}^{m}}\left(|u_{0}(t,x)|^{2}\mathrm{d}x+h\sum_{k=1}^{N}|u_{k}(t,x)|^{2}\right)\mathrm{d}x.

(2.10)

Here the index $N$ in $m_{N}(t)$ represents that the metal continuum in this model is discretized into $N$ orbitals. Similarly, we can define the system energy $E_{N}(t)$ :

	$\displaystyle E_{N}(t)$	$\displaystyle=\int_{\mathbb{R}^{m}}\frac{\varepsilon^{2}}{2}\left\|\nabla_{x}u_{0}\right\|^{2}+U_{0}(x)\left\|u_{0}\right\|^{2}+h\sum_{k=1}^{N}\left(\frac{\varepsilon^{2}}{2}\|\nabla_{x}u_{k}\|^{2}+\left(U_{1}(x)+\mathcal{E}\right)\left\|u_{k}\right\|^{2}\right)\mathrm{d}x$		(2.11)
		$\displaystyle+\varepsilon\int_{\mathbb{R}^{m}}h\sum_{k=1}^{N}\left(V(\mathcal{E}_{k},x)\bar{u}_{0}(t,x)u_{k}(t,x)+\overline{V}(\mathcal{E}_{k},x)u_{0}(t,x)\bar{u}_{k}(t,x)\right)\mathrm{d}x.$		(2.11)

It’s also straightforward to verify that solutions to 1.1 satisfy the conservation of mass and energy:

m_{N}(t)=m_{N}(0),\quad E_{N}(t)=E_{N}(0),\quad\forall t\geq 0.

(2.12)

Dynamics 1.1 could be regarded as a $(N+1)-$ level semiclassical Schrödinger equation system. As we have mentioned in section 1, our goal is to design an algorithm for 1.1, whose computational cost does not essentially depend on both the semiclassical parameter $\varepsilon$ and the number of metal orbitals $N$ .

3 Algorithm

In this section, we introduce the Frozen Gaussian Sampling (FGS) algorithm for the system 1.1.

3.1 Review of Frozen Gaussian representation for one-level system

Before dealing with the multi-level system 1.1, we first give a brief review of the Frozen Gaussian representation for a one-level system. Consider the following semiclassical Schrödinger equation:

i\varepsilon\frac{\partial}{\partial t}u(t,x)=-\frac{\varepsilon^{2}}{2}\Delta u(t,x)+U(x)u(t,x),\quad u(0,x)=u_{\mathrm{in}}(x).

(3.1)

It is known that the Frozen Gaussian approximation of 3.1 is an approximation derived from asymptotic analysis with $O(\varepsilon)$ error [38]. It is based on an integral representation of Frozen Gaussian wave packets on phase space:

u_{\mathrm{FG}}(t,x)=\frac{1}{(2\pi\varepsilon)^{\frac{3m}{2}}}\int_{\mathbb{R}^{2m}}A(t,q,p)\mathrm{e}^{\frac{\mathrm{i}\Theta(t,x,q,p)}{\varepsilon}}\mathrm{d}q\mathrm{~{}d}p.

(3.2)

where the phase function $\Theta$ is

\Theta(t,x,q,p)=S(t,q,p)+P(t,q,p)\cdot(x-Q(t,q,p))+\frac{\mathrm{i}}{2}|x-Q(t,q,p)|^{2}.

(3.3)

The evolution of $(Q(t,q,p),P(t,q,p))$ is governed by the Hamiltonian flow with $(q,p)$ as its initial value:

	$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}Q$	$\displaystyle=\partial_{p}h(Q,P),\quad Q(0,q,p)=q.$		(3.4)
	$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}P$	$\displaystyle=-\partial_{q}h(Q,P),\quad P(0,q,p)=p.$		(3.4)

where $h(q,p)=\frac{1}{2}|p|^{2}+U(q)$ is the classical Hamiltonian. The dynamics of the action $S(t,q,p)$ and amplitude $A(t,q,p)$ are

\frac{\mathrm{d}}{\mathrm{d}t}S(t,q,p)=\frac{1}{2}|P|^{2}-U(Q),\quad S(0,q,p)=0.

(3.5)

	$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}A(t,q,p)$	$\displaystyle=\frac{A}{2}\operatorname{tr}\left(Z^{-1}\left(\partial_{z}P-\mathrm{i}\partial_{z}Q\nabla_{Q}^{2}U(Q)\right)\right),$		(3.6)
	$\displaystyle A(0,q,p)$	$\displaystyle=2^{\frac{m}{2}}\int_{\mathbb{R}^{m}}u_{\mathrm{in}}(y)e^{\frac{\mathrm{i}}{\varepsilon}\left(-p\cdot(y-q)+\frac{i}{2}\|y-q\|^{2}\right)}\mathrm{d}y.$		(3.6)

Here we use the short hand notation:

\partial_{z}=\partial_{q}-i\partial_{p}\quad\text{ and }\quad Z=\partial_{z}(Q+iP).

(3.7)

We emphasize here that the initial condition $u_{\text{in}}$ has already been encoded into $A(0,q,p)$ .

The above Frozen Gaussian dynamics are derived using asymptotic matching [38, 33]. The Frozen Gaussian representation for the initial value $u_{\mathrm{in}}$ has the following property, which we quote directly from [38]:

Theorem 3.1.

For $u\in L^{2}\left(\mathbb{R}^{m}\right)$ , we have the following decomposition

u_{\mathrm{in}}(x)=\frac{1}{(2\pi\varepsilon)^{\frac{3m}{2}}}\int_{\mathbb{R}^{2m}}A(0,q,p)\mathrm{e}^{\frac{\mathrm{i}\Theta(0,x,q,p)}{\varepsilon}}\mathrm{d}q\mathrm{~{}d}p,

(3.8)

where $A(0,q,p)$ is given as in 3.6.

And, for the Frozen Gaussian representation 3.2 of the one-level semiclassical dynamics 3.1, we have the following error estimate [38, 33]:

Theorem 3.2.

Let $u(t,x)$ be the exact solution for the one-level semiclassical dynamics 3.1. Assume that $U(x)\in C^{\infty}\left(\mathbb{R}^{m}\right)$ . For any $t>0$ , we have the following $L^{2}$ error estimate for the Frozen Gaussian representation 3.2

\left\|u_{\mathrm{FG}}(t,\cdot)-u(t,\cdot)\right\|_{L^{2}}\leq C_{m}(t)\varepsilon\left\|u_{\mathrm{in}}(\cdot)\right\|_{L^{2}},

(3.9)

where $C_{m}(t)$ is a constant independent of $\varepsilon$ .

3.2 Integral representation based on the surface hopping ansatz

The integral representation of the solution of 1.1 is based on the following surface hopping ansatz: (as illustrated in Figure 1)

1.

A wave packet starts from the molecule orbital.
2.

While the wave packet is on the molecule orbital, it might hop onto one of the metal surfaces; while the wave packet is on one of the metal surfaces, it might hop back onto the molecule orbital. No hopping between different metal surfaces is allowed.
3.

When the wave packet hasn’t hopped out of the current energy surface, it would propagate subject to the dynamics specified by the current energy surface.

Refer to caption — Figure 1: Illustration of the surface hopping ansatz: $k_{1},k_{2},k_{3}\in\{1,2,\cdots,N\}$ represents the indices of the metal orbitals that the wave packet has hopped onto. $u_{0}^{(0)}$ represents the wave packet that hasn’t ever hopped; $u_{k_{1}}^{(1)}$ represents the wave packet that has hopped once and is on the $k_{1}$ -th metal surface; $u_{0,k_{1}}^{(2)}$ represents the wave packet that has hopped twice and was on the $k_{1}$ -th metal surface before the second jump; $u_{k_{2},k_{1}}^{(3)}$ represents the wave packet that has hopped three times, is currently on the $k_{2}$ -th metal surface and was on the $k_{1}$ -th metal surfaces in history. The meaning of $u_{0,k_{2:1}}^{(4)}$ and $u_{k_{3},k_{2:1}}^{(5)}$ could be determined similarly.

The hopping probability and the detailed wave packet dynamics will be specified later. For now, let us introduce some notations. Let $k_{1}\in\{1,2,\cdots,N\}$ be the index of the metal orbital that the wave packet is on when it is the first time for the wave packet to ever hop onto the metal surfaces. In the same way, we can define $k_{2},k_{3},\cdots$ . The wave packet on the molecule orbital is denoted as

u^{(2n)}_{0,k_{n:1}},\quad n=0,1,2,\cdots,\quad k_{1},k_{2},\cdots,k_{n}\in\{1,2,\cdots,N\},

where the superscript $2n$ indicates the times that the wave packet has already jumped (which must be an even integer since the wave packet has to hop out and back to be on the molecule orbital), the subscript $0$ means that the wave packet is currently on the molecule orbital, and $k_{n:1}$ is the short hand notation for $k_{1},k_{2},\cdots,k_{n}$ , which documents the indices of all metal orbitals that this wave packet has ever jumped onto. Similarly, the wave packet on the metal surfaces are denoted as

u^{(2n+1)}_{k,k_{n:1}},\quad n=0,1,2,\cdots,\quad k,k_{1},k_{2},\cdots,k_{n}\in\{1,2,\cdots,N\},

where the superscript is still the hopping times (an odd integer) and the subscript is the current metal orbital index $k$ and the historical metal orbitals indices $k_{n:1}$ . The wave packets could be written in the following unified formulation:

u_{k,k_{[\nu/2]:1}}^{(\nu)},\quad\nu=0,1,2,\cdots,\quad k_{1},k_{2},\cdots,k_{n}\in\{1,2,\cdots,N\},

where $k=0$ when $\nu$ is even, which means the wave packet is on the molecule orbital, and $k\in\{1,\cdots,N\}$ when $\nu$ is odd, which means that the wave packet is on the $k$ -th metal surface. Later, when we introduce other quantities of the wave packet, we will use the same index system.

Now we are ready to present the surface hopping ansatz of the wavefunction:

	$\displaystyle{u}_{0,\text{FG}}(t,x)$	$\displaystyle=u_{0}^{(0)}(t,x)+\sum_{k_{1}=1}^{N}u_{0,k_{1}}^{(2)}(t,x)+\sum_{k_{1},k_{2}=1}^{N}u_{0,k_{2:1}}^{(4)}(t,x)+\cdots,$		(3.10)
	$\displaystyle{u}_{k,\text{FG}}(t,x)$	$\displaystyle=u_{k}^{(1)}(t,x)+\sum_{k_{1}=1}^{N}u_{k,k_{1}}^{(3)}(t,x)+\sum_{k_{1},k_{2}=1}^{N}u_{k,k_{2:1}}^{(5)}(t,x)+\cdots,k=1:N.$		(3.10)

Since $u_{0}^{(0)}(t,x)$ represents the wave packet that only propagates on $U_{0}$ and hasn’t switched to the metal surfaces, it is thus given by the same ansatz as in one-level system:

u_{0}^{(0)}(t,x)=\frac{1}{(2\pi\varepsilon)^{3m/2}}\int A_{0}^{(0)}\left(t,z_{0}\right)\exp\left(\frac{\mathrm{i}}{\varepsilon}\Theta_{0}^{(0)}\left(t,z_{0},x\right)\right)\mathrm{d}z_{0}.

(3.11)

Here $z_{0}=(q_{0},p_{0})$ , and

\Theta_{0}^{(0)}(t,x,q,p)=S_{0}^{(0)}(t,q,p)+P_{0}^{(0)}(t,q,p)\cdot(x-Q_{0}^{(0)}(t,q,p))+\frac{\mathrm{i}}{2}|x-Q_{0}^{(0)}(t,q,p)|^{2},

(3.12)

and the dynamics of $Q_{0}^{(0)}$ , $P_{0}^{(0)}$ , $A_{0}^{(0)}$ , $S_{0}^{(0)}$ will be specified below.

For $\nu>0$ , let us denote $T_{\nu:1}=\left(t_{\nu},\ldots,t_{1}\right)$ as a sequence of times in decreasing order:

t_{0}=0\leq t_{1}\leq t_{2}\leq\cdots\leq t_{\nu}\leq t=t_{\nu+1}.

And now we are ready to introduce our ansatz for $\nu=2n$ and $\nu=2n+1$ :

$\displaystyle u_{0,k_{n:1}}^{(2n)}(t,x)=\frac{1}{(2\pi\varepsilon)^{3m/2}}$	$\displaystyle\int\mathrm{d}z_{0}\int_{0}^{t}\mathrm{~{}d}t_{2n}\int_{0}^{t_{2n}}\mathrm{~{}d}t_{2n-1}\cdots\int_{0}^{t_{2}}\mathrm{d}t_{1}$	(3.13)
	$\displaystyle\left(\prod_{j=1}^{n}\left(h\tau_{k_{j},k_{j-1:1}}^{(2j-1)}\left(T_{2j-1:1},z_{0}\right)\tau_{{0,k_{j:1}}}^{(2j)}\left(T_{2j:1},z_{0}\right)\right)\right)$
	$\displaystyle\times A_{0,k_{n:1}}^{(2n)}\left(t,T_{2n:1},z_{0}\right)\exp\left(\frac{\mathrm{i}}{\varepsilon}\Theta_{0,k_{n:1}}^{(2n)}\left(t,T_{2n:1},z_{0},x\right)\right).$

$\displaystyle u_{k,k_{n:1}}^{(2n+1)}(t,x)=$	$\displaystyle\frac{1}{(2\pi\varepsilon)^{3m/2}}\int\mathrm{d}z_{0}\int_{0}^{t}\mathrm{~{}d}t_{2n+1}\int_{0}^{t_{2n}}\mathrm{~{}d}t_{2n}\cdots\int_{0}^{t_{2}}\mathrm{d}t_{1}$	(3.14)
	$\displaystyle\tau_{k,k_{n:1}}^{(2n+1)}\left(T_{2n+1:1},z_{0}\right)\times\left(\prod_{j=1}^{n}\left(h\tau_{k_{j},k_{j-1:1}}^{(2j-1)}\left(T_{2j-1:1},z_{0}\right)\tau_{0,k_{j:1}}^{(2j)}\left(T_{2j:1},z_{0}\right)\right)\right)$
	$\displaystyle\times A_{k,k_{n:1}}^{(2n+1)}\left(t,T_{2n+1:1},z_{0}\right)\exp\left(\frac{\mathrm{i}}{\varepsilon}\Theta_{k,k_{n:1}}^{(2n+1)}\left(t,T_{2n+1:1},z_{0},x\right)\right).$

Here $h=(\mathcal{E}_{b}-\mathcal{E}_{a})/N$ . For both $\nu=2n$ and $\nu=2n+1$ , we have

	$\displaystyle\Theta_{k,k_{[{\nu}/{2}]:1}}^{(\nu)}\left(t,T_{\nu:1},z_{0},x\right)=$	$\displaystyle S_{k,k_{[{\nu}/{2}]:1}}^{(\nu)}\left(t,T_{\nu:1},z_{0}\right)+\frac{\mathrm{i}}{2}\left\|x-Q_{k,k_{[{\nu}/{2}]:1}}^{(\nu)}\left(t,T_{\nu:1},z_{0}\right)\right\|^{2}$		(3.15)
	$\displaystyle+$	$\displaystyle P_{k,k_{[{\nu}/{2}]:1}}^{(\nu)}\left(t,T_{\nu:1},z_{0}\right)\cdot\left(x-Q_{k,k_{[{\nu}/{2}]:1}}^{(\nu)}\left(t,T_{\nu:1},z_{0}\right)\right).$		(3.15)

Compared with the Frozen Gaussian representation for one-level system 3.2, the above ansatz includes the integration over all possible jumping times, while the parameters $\tau_{k,k_{[\nu/2]:1}}^{(\nu)}$ accounting for these jumps are to be determined.

Based on the ansatz, we can determine the dynamics of $S_{k,k_{[{\nu}/{2}]:1}}^{(\nu)},P_{k,k_{[{\nu}/{2}]:1}}^{(\nu)}$ , $Q_{k,k_{[{\nu}/{2}]:1}}^{(\nu)}$ , $A_{k,k_{[{\nu}/{2}]:1}}^{(\nu)}$ , and the value of $\tau_{k,k_{[\nu/2]:1}}^{(\nu)}$ . The detailed derivation is based on the asymptotic matching, and is stated in A. The dynamics of $S_{k,k_{[{\nu}/{2}]:1}}^{(\nu)},P_{k,k_{[{\nu}/{2}]:1}}^{(\nu)},Q_{k,k_{[{\nu}/{2}]:1}}^{(\nu)}$ and $A_{k,k_{[{\nu}/{2}]:1}}^{(\nu)}$ for $t\in[t_{\nu},t_{\nu+1}]$ are

$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}Q_{k,k_{[{\nu}/{2}]:1}}^{(\nu)}=$	$\displaystyle P_{k,k_{[{\nu}/{2}]:1}}^{(\nu)},$	(3.16)
$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}P_{k,k_{[{\nu}/{2}]:1}}^{(\nu)}=$	$\displaystyle-\nabla_{Q}\widetilde{U}\left(t,Q_{k,k_{[{\nu}/{2}]:1}}^{(\nu)}\right),$
$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}S_{k,k_{[{\nu}/{2}]:1}}^{(\nu)}=$	$\displaystyle\frac{1}{2}\left(P_{k,k_{[{\nu}/{2}]:1}}^{(\nu)}\right)^{2}-\widetilde{U}\left(t,Q_{k,k_{[{\nu}/{2}]:1}}^{(\nu)}\right),$
$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}A_{k,k_{[{\nu}/{2}]:1}}^{(\nu)}=$	$\displaystyle\frac{1}{2}A_{k,k_{[{\nu}/{2}]:1}}^{(\nu)}\operatorname{tr}\left(\left(Z^{(\nu)}\right)^{-1}\left(\partial_{z}P_{k,k_{[{\nu}/{2}]:1}}^{(\nu)}-i\partial_{z}Q_{k,k_{[{\nu}/{2}]:1}}^{(\nu)}\nabla_{Q}^{2}\widetilde{U}\left(t,Q_{k,k_{[{\nu}/{2}]:1}}^{(\nu)}\right)\right)\right),$

where

\widetilde{U}(t,Q)=\left\{\begin{aligned} U_{0}(Q),\quad&t\in[t_{2j},t_{2j+1}]\\ U_{1}(Q)+\mathcal{E}_{k_{j}},\quad&t\in[t_{2j-1},t_{2j}]\end{aligned}\right.,

(3.17)

and the initial value of time interval $[0,t_{1}]$ is

Q_{0}^{(0)}(0,z_{0})=q,\quad P_{0}^{(0)}(0,z_{0})=q,\quad S_{0}^{(0)}(0,z_{0})=0,

(3.18)

A_{0}^{(0)}(0,z_{0})=2^{\frac{m}{2}}\int_{\mathbb{R}^{m}}u_{0,\mathrm{in}}(y)e^{\frac{\mathrm{i}}{\varepsilon}\left(-p\cdot(y-q)+\frac{i}{2}|y-q|^{2}\right)}\mathrm{d}y.

(3.19)

and the initial value on time interval $[t_{\nu},t_{\nu+1}]$ is chosen to makes sure all of them are continuous functions.

$\displaystyle Q_{k,k_{\lceil{(\nu+1)}/{2}\rceil:1}}^{(\nu+1)}\left(t_{\nu+1},T_{\nu+1:1},z_{0}\right)$	$\displaystyle=Q_{k,k_{[{\nu}/{2}]:1}}^{(\nu)}\left(t_{\nu+1},T_{\nu:1},z_{0}\right),$	(3.20)
$\displaystyle P_{k,k_{\lceil{(\nu+1)}/{2}\rceil:1}}^{(\nu+1)}\left(t_{\nu+1},T_{\nu+1:1},z_{0}\right)$	$\displaystyle=P_{k,k_{[{\nu}/{2}]:1}}^{(\nu)}\left(t_{\nu+1},T_{\nu:1},z_{0}\right),$
$\displaystyle S_{k,k_{\lceil{(\nu+1)}/{2}\rceil:1}}^{(\nu+1)}\left(t_{\nu+1},T_{\nu+1:1},z_{0}\right)$	$\displaystyle=S_{k,k_{[{\nu}/{2}]:1}}^{(\nu)}\left(t_{\nu+1},T_{\nu:1},z_{0}\right),$
$\displaystyle A_{k,k_{\lceil{(\nu+1)}/{2}\rceil:1}}^{(\nu+1)}\left(t_{\nu+1},T_{\nu+1:1},z_{0}\right)$	$\displaystyle=A_{k,k_{[{\nu}/{2}]:1}}^{(\nu)}\left(t_{\nu+1},T_{\nu:1},z_{0}\right).$

Now the only thing left to be specified is $\tau_{k,k_{[\nu/2]:1}}^{(\nu)}$ . The detailed derivation of $\tau_{k}^{\nu}$ could also be seen in appendix A. We have

	$\displaystyle\tau_{0,k_{n:1}}^{(2n)}\left(t,T_{2n-1:1},z_{0}\right)$	$\displaystyle=-\mathrm{i}V\left(\mathcal{E}_{k_{n}},Q_{k_{n},k_{n-1:1}}^{(2n-1)}\left(t,T_{2n-1:1},z_{0}\right)\right),$		(3.21)
	$\displaystyle\tau_{k,k_{n:1}}^{(2n+1)}\left(t,T_{2n:1},z_{0}\right)$	$\displaystyle=-\mathrm{i}{\left(\overline{V}\left(\mathcal{E}_{k},Q_{0,k_{n:1}}^{(2n)}\left(t,T_{2n:1},z_{0}\right)\right)\right)}.$		(3.21)

We refer to $\tau_{k,k_{[\nu/2]:1}}^{(\nu)}$ as the hopping coefficient associated with the $k$ -th level, and we will see in section 3.3 that $\tau_{k,k_{[\nu/2]:1}}^{(\nu)}$ is used to define the hopping rate onto or away from the $k$ -th level in the jump process that we will construct.

3.3 Stochastic interpretation

In order to implement the Frozen Gaussian representation with an efficient algorithm, we introduce its stochastic interpretation. We will construct a stochastic process $\widetilde{z}_{t}$ , such that the above integral representation could be represented as taking expectations over all the path $\widetilde{z}_{t}$ .

The stochastic process $\widetilde{z}_{t}=(z_{t},l_{t})$ is defined on the extended phase space $\mathbb{R}^{2m}\times\{0,1,\cdots,N\}$ . Here $l_{t}\in\{0,1,\cdots,N\}$ represents the energy surface that the trajectory is currently on at time $t$ . The evolution of $z_{t}$ is determined by the energy surface specified by $l_{t}$ :

\mathrm{d}z_{t}=\left\{\begin{aligned} \left(p_{t},-\nabla_{q}U_{0}\left(q_{t}\right)\right)\mathrm{d}t,&\quad l_{t}=0,\\ \left(p_{t},-\nabla_{q}U_{1}\left(q_{t}\right)\right)\mathrm{d}t,&\quad l_{t}=1,\cdots,N.\end{aligned}\right.

(3.22)

Here we take advantage of the fact that $\nabla_{q}\left(U_{1}(q)+\mathcal{E}_{k}\right)=\nabla_{q}U_{1}(q)$ since $\mathcal{E}_{k}$ is a constant. $z_{t}$ is coupled with $l_{t}$ , which follows a jump process, with the following infinitesimal transition rate:

\mathbb{P}\left(l_{t+\delta t}=k^{\prime}\mid l_{t}=k,z_{t}=z\right)=\delta_{kk^{\prime}}+\lambda_{kk^{\prime}}(z)\delta t+o(\delta t),

(3.23)

where

\lambda_{kk^{\prime}}(z)=\left\{\begin{aligned} h\left|V\left(\mathcal{E}_{k},q_{t}\right)\right|,\quad&k=0,k^{\prime}=1,\cdots,N\\ \left|V\left(\mathcal{E}_{k},q_{t}\right)\right|,\quad&k^{\prime}=0,k=1,\cdots,N\\ -h\sum_{k=1}^{N}\left|V\left(\mathcal{E}_{k},q_{t}\right)\right|,\quad&k=k^{\prime}=0\\ -\left|V\left(\mathcal{E}_{k},q_{t}\right)\right|,\quad&k=k^{\prime}=1,\cdots,N\\ 0,\quad&{\text{ otherwise}}\end{aligned}\right.\quad,

(3.24)

	$\displaystyle\lambda(z)$	$\displaystyle=\left(\begin{array}[]{cccc}\lambda_{00}(z)&\lambda_{01}(z)&\cdots&\lambda_{0N}(z)\\ \lambda_{10}(z)&\lambda_{11}(z)&\cdots&\lambda_{1N}(z)\\ \vdots&\vdots&\ddots&\vdots\\ \lambda_{N0}(z)&\lambda_{N1}(z)&\cdots&\lambda_{NN}(z)\end{array}\right)$		(3.25)
		$\displaystyle=\left(\begin{array}[]{cccc}-h\sum_{k=1}^{N}\left\|V\left(\mathcal{E}_{k},q_{t}\right)\right\|&h\|V(\mathcal{E}_{1},q_{t})\|&\cdots&h\|V(\mathcal{E}_{N},q_{t})\|\\ \left\|V\left(\mathcal{E}_{1},q_{t}\right)\right\|&-\left\|V\left(\mathcal{E}_{1},q_{t}\right)\right\|&\cdots&0\\ \vdots&\vdots&\ddots&\vdots\\ \left\|V\left(\mathcal{E}_{N},q_{t}\right)\right\|&0&\cdots&-\left\|V\left(\mathcal{E}_{N},q_{t}\right)\right\|\end{array}\right).$		(3.25)

With $\lambda(z)$ , we can write down now the forward Kolmogorov equation specified by the above hopping strategy

\frac{\partial}{\partial t}\rho_{k}(t,z)+\left\{h_{k},\rho_{k}(t,z)\right\}=\sum_{m=0}^{N}\lambda_{mk}(z)\rho_{k}(t,z),

(3.26)

where $\rho_{k}(t,z)$ represents the phase space density function of the $k$ -th level, and $\{\cdot\}$ is the Poisson bracket:

\{h,F\}=\partial_{p}h\cdot\partial_{q}F-\partial_{q}h\cdot\partial_{p}F.

Note that we have $\tau_{k,k_{[\nu/2]:1}}^{(\nu)}\left(t,T_{\nu:1},z_{0}\right)=\left|V(\mathcal{E}_{k},Q_{k_{[\nu/2]:1}}^{(\nu-1)})\right|$ . For $\nu=2n+1$ , let us define $\widetilde{\lambda}_{k_{n:1}}^{(2n+1)}\left(t,T_{2n:1},z_{0}\right)$ for future purpose:

\widetilde{\lambda}_{k_{n:1}}^{(2n+1)}\left(t,T_{2n:1},z_{0}\right)=\sum_{k=1}^{N}h\left|\tau_{k,k_{n:1}}^{(2n+1)}\left(t,T_{2n:1},z_{0}\right)\right|.

(3.27)

Let $\mathbb{P}_{2n}(t,z_{0})$ be the probability that there is $2n$ jumps in $[0,t]$ . Let $\mathbb{P}_{2n+1}(t,k,z_{0})$ be the probability that there is $(2n+1)$ jumps in $[0,t]$ and $l_{t}=k$ . By definition, we have

1=\sum_{n=0}^{\infty}\mathbb{P}_{2n}(t,z_{0})+\sum_{n=0}^{\infty}\sum_{k=1}^{N}\mathbb{P}_{2n+1}\left(t,k,z_{0}\right).

(3.28)

By the properties of jumping process, we have

	$\displaystyle\mathbb{P}_{2n}(t,z_{0})=\sum_{k_{1},\cdots,k_{n}=1}^{N}\iint\cdots\int_{[0,t]^{2n}}$	$\displaystyle\mathrm{d}t_{2n}\cdots\mathrm{d}t_{1}$		(3.29)
		$\displaystyle\rho^{(2n)}_{0,k_{n:1}}\left(T_{2n:1},z_{0}\right)\Theta\left(0<t_{1}<\cdots<t_{2n}<t\right),$		(3.29)

	$\displaystyle\mathbb{P}_{2n+1}(t,k,z_{0})=\sum_{k_{1},\cdots,k_{n},k=1}^{N}\iint\cdots\int_{[0,t]^{2n+1}}\mathrm{d}t_{2n+1}$	$\displaystyle\cdots\mathrm{d}t_{1}$		(3.30)
	$\displaystyle\rho^{(2n+1)}_{k,k_{n:1}}\left(T_{2n+1:1},z_{0}\right)$	$\displaystyle\Theta\left(0<t_{1}<\cdots<t_{2n+1}<t\right).$		(3.30)

where

		$\displaystyle\rho^{(2n)}_{0,k_{n:1}}\left(T_{2n:1},z_{0}\right)=\prod_{j=0}^{n}\left(\exp\left(-\int_{t_{2j}}^{t_{2j+1}}\mathrm{d}s_{2j+1}\widetilde{\lambda}_{k_{n:1}}^{(2n+1)}\left(s_{2j+1},T_{2n:1},z_{0}\right)\right)\right)$		(3.31)
		$\displaystyle\times\left(\prod_{j=1}^{n}h\left\|\tau_{k,k_{j-1:1}}^{(2j-1)}\left(T_{2j-1:1}\right)\right\|\mathrm{e}^{\left(-\int_{t_{2j-1}}^{t_{2j}}\left\|\tau_{0,k_{j:1}}^{(2j)}\left(s_{2j},T_{2j-1:1}\right)\right\|\mathrm{d}s_{2j}\right)}\left\|\tau_{0,k_{j:1}}^{(2j)}\left(T_{2j:1}\right)\right\|\right),$		(3.31)

and

		$\displaystyle\rho^{(2n+1)}_{k,k_{n:1}}\left(T_{2n+1:1},z_{0}\right)=h\left\|\tau_{k,k_{j:1}}^{(2n+1)}\left(T_{2n+1:1}\right)\right\|\exp\left(-\int_{t_{2n+1}}^{t}\mathrm{d}s\tau_{k,k_{n:1}}^{(2n+2)}\left(s,T_{2j-1:1},z_{0}\right)\right)$		(3.32)
		$\displaystyle\times\prod_{j=0}^{n}\left(\exp\left(-\int_{t_{2j}}^{t_{2j+1}}\mathrm{d}s_{2j+1}\widetilde{\lambda}_{k_{n:1}}^{(2n+1)}\left(s_{2j+1},T_{2n:1},z_{0}\right)\right)\right)$
		$\displaystyle\times\left(\prod_{j=1}^{n}h\left\|\tau_{k,k_{j-1:1}}^{(2j-1)}\left(T_{2j-1:1}\right)\right\|\mathrm{e}^{\left(-\int_{t_{2j-1}}^{t_{2j}}\left\|\tau_{0,k_{j:1}}^{(2j)}\left(s_{2j},T_{2j-1:1}\right)\right\|\mathrm{d}s_{2j}\right)}\left\|\tau_{0,k_{j:1}}^{(2j)}\left(T_{2j:1}\right)\right\|\right).$

Based on these, we have :

	$\displaystyle u_{0,k_{n:1}}^{(2n)}(t,x)$	$\displaystyle=\frac{1}{(2\pi\varepsilon)^{3m/2}}\int\mathrm{d}z_{0}\int_{0}^{t}\mathrm{~{}d}t_{2n}\int_{0}^{t_{2n}}\mathrm{~{}d}t_{2n-1}\cdots\int_{0}^{t_{2}}\mathrm{d}t_{1}\rho^{(2n)}_{0,k_{n:1}}(T_{2n:1})$		(3.33)
		$\displaystyle\times A_{0,k_{n:1}}^{(2n)}\left(t,T_{2n:1},z_{0}\right)\exp\left(\frac{\mathrm{i}}{\varepsilon}\Theta_{0,k_{n:1}}^{(2n)}\left(t,T_{2n:1},z_{0},x\right)\right)\Gamma_{0,k_{n:1}}^{(2n)}(t,T_{2n:1}),$		(3.33)

		$\displaystyle u_{k,k_{n:1}}^{(2n+1)}(t,x)=\frac{1}{(2\pi\varepsilon)^{3m/2}}\frac{1}{h}\int\mathrm{d}z_{0}\int_{0}^{t}\mathrm{~{}d}t_{2n+1}\int_{0}^{t_{2n}}\mathrm{~{}d}t_{2n}\cdots\int_{0}^{t_{2}}\mathrm{d}t_{1}\rho^{(2n+1)}_{k,k_{n:1}}\left(T_{2n+1:1}\right)$		(3.34)
		$\displaystyle\times A_{k,k_{n:1}}^{(2n+1)}\left(t,T_{2n+1:1},z_{0}\right)\exp\left(\frac{\mathrm{i}}{\varepsilon}\Theta_{k,k_{n:1}}^{(2n+1)}\left(t,T_{2n+1:1},z_{0},x\right)\right)\Gamma_{k,k_{n:1}}^{(2n+1)}(t,T_{2n+1:1}),$		(3.34)

where

		$\displaystyle\Gamma_{0,k_{n:1}}^{(2n)}(t,T_{2n:1},z_{0})=\prod_{j=0}^{n}\left(\exp\left(\int_{t_{2j}}^{t_{2j+1}}\mathrm{d}s_{2j+1}\widetilde{\lambda}_{k_{n:1}}^{(2n+1)}\left(s_{2j+1},T_{2n:1},z_{0}\right)\right)\right)$		(3.35)
		$\displaystyle\times\prod_{j=1}^{n}\frac{\tau_{k_{j},k_{j-1:1}}^{(2j-1)}\left(T_{2j-1:1}\right)\tau_{{0,k_{j:1}}}^{(2j)}\left(T_{2j:1}\right)\exp\left({\left(\int_{t_{2j-1}}^{t_{2j}}\left\|\tau_{0,k_{j:1}}^{(2j)}\left(s_{2j},T_{2j-1:1}\right)\right\|\mathrm{d}s_{2j}\right)}\right)}{\left\|\tau_{k_{j},k_{j-1:1}}^{(2j-1)}\left(T_{2j-1:1}\right)\right\|\left\|\tau_{0,k_{j:1}}^{(2j)}\left(T_{2j:1}\right)\right\|},$		(3.35)

		$\displaystyle\Gamma_{k,k_{n:1}}^{(2n+1)}(t,T_{2n+1:1},z_{0})=\prod_{j=0}^{n}\left(\exp\left(\int_{t_{2j}}^{t_{2j+1}}\mathrm{d}s_{2j+1}\widetilde{\lambda}_{k_{n:1}}^{(2n+1)}\left(s_{2j+1},T_{2n:1},z_{0}\right)\right)\right)$		(3.36)
		$\displaystyle\times\exp\left(\int_{t_{2n+1}}^{t}\tau_{0,k_{n+1:1}}^{(2n+2)}\left(s,T_{2j-1:1}\right)\mathrm{d}s\right)\times\frac{\tau_{k,k_{n:1}}^{(2n+1)}\left(T_{2n+1:1}\right)}{\left\|\tau_{k,k_{n:1}}^{(2n+1)}\left(T_{2n+1:1}\right)\right\|}$
		$\displaystyle\times\prod_{j=1}^{n}\frac{\tau_{k_{j},k_{j-1:1}}^{(2j-1)}\left(T_{2j-1:1}\right)\tau_{{0,k_{j:1}}}^{(2j)}\left(T_{2j:1}\right)\exp\left({\left(\int_{t_{2j-1}}^{t_{2j}}\left\|\tau_{0,k_{j:1}}^{(2j)}\left(s_{2j},T_{2j-1:1}\right)\right\|\mathrm{d}s_{2j}\right)}\right)}{\left\|\tau_{k,k_{j-1:1}}^{(2j-1)}\left(T_{2j-1:1}\right)\right\|\left\|\tau_{0,k_{j:1}}^{(2j)}\left(T_{2j:1}\right)\right\|}.$

Finally, we are ready to present the Frozen Gaussian representation as an average over stochastic path. Let

\mathbf{U}(t,x)=\left(\begin{array}[]{c}u_{0}(t,x)\\ u_{1}(t,x)\\ \cdots\\ u_{N}(t,x)\end{array}\right),

(3.37)

And let $\mathbf{U}_{\text{FG}}(t,x)$ be the Frozen Gaussian representation of $\mathbf{U}(t,x)$ . Rewriting 3.33, 3.34 with 3.35 and 3.36, we have

\mathbf{U}_{\text{FG}}(t,x)=\mathbb{E}_{\widetilde{z}_{t}}\left(\Lambda(t,x;\widetilde{z}_{t})\right),

(3.38)

where

\Lambda(t,x;\widetilde{z}_{t})=W(t,x;\widetilde{z}_{t}){\mathbb{V}}^{(k)},

(3.39)

W(t,x;\widetilde{z}_{t})=\frac{1}{(2\pi\varepsilon)^{3m/2}}\frac{1}{\pi(z_{0})}A_{k,k_{[{\nu}/{2}]:1}}^{(\nu)}\exp\left(\frac{\mathrm{i}}{\varepsilon}\Theta_{k,k_{[{\nu}/{2}]:1}}^{(\nu)}\right)\Gamma_{k,k_{[\nu/2]:1}}^{(\nu)}.

(3.40)

Here $k$ is the index of surface that the wave packet is currently on, $\pi(z_{0})$ is the probability density that $z_{0}$ is sampled from, and $\mathbb{V}^{(\nu)}$ is an $(N+1)\times 1$ vector, its $j$ -th component is ( $j=0,1,\cdots,N$ )

{\mathbb{V}}_{j}^{(\nu)}=\left\{\begin{aligned} \delta_{0j},\quad\nu\text{ is even},\\ \frac{1}{h}\delta_{k_{\lceil{\nu}/{2}\rceil}j},\quad\nu\text{ is odd}.\end{aligned}\right.

(3.41)

Note that $\Lambda(t,x;\widetilde{z}_{t})$ is a $(N+1)$ -dimensional vector-valued function and $W(t,x;\widetilde{z}_{t})$ is a scalar-valued function.

We emphasize that $\Lambda(t,x;\widetilde{z}_{t})$ could be viewed as a functional of the stochastic trajectory $\widetilde{z}_{t}$ . As a result, according to 3.38, the wavefunction $\mathbf{U}_{\text{FG}}$ is in fact the average of $\Lambda(t,x;\widetilde{z}_{t})$ 3.39 over the stochastic path $\widetilde{z}_{t}$ . From an algorithmic perspective, the sampling of the stochastic path $\widetilde{z}_{t}$ is achieved by first sampling its initial value $z_{0}\sim\pi(z)$ and then a stochastic evolution subject to the Hamiltonian flow 3.22 and the jump process 3.23.

Remark 3.3.

Later We’ll use $\Lambda_{k}(t,x;\widetilde{z}_{t})$ to represent the $k$ -th component of $\Lambda$ for $k=0,1,\cdots,N$ . We will also use $\mathbf{U}_{k,\operatorname{FG}}(t,x)$ to represent the $k$ -th component of $\mathbf{U}_{\operatorname{FG}}(t,x)$ for $k=0,1,\cdots,N$ .

3.4 FGS Algorithm

Based on the stochastic interpretation 3.38, we are able to constuct the Frozen Gaussian Sampling (FGS) algorithm, as illustrated in Figure 2.

After drawing $M$ independent initial wave packets by a properly chosen phase space distribution, each wavepacket evolves subject to the Hamiltonian flow 3.22 and hops between energy surfaces following the jump process 3.23. Then at time $t$ , we can reconstruct the solution for system 1.1 with the average of $M$ wave packets:

\mathbf{U}_{\text{FG}}(t,x)\approx\frac{1}{M}\sum_{j=1}^{M}\Lambda^{(j)}\left(t,x;\widetilde{z}_{t}^{(j)}\right).

(3.42)

Here the superscript $(j)$ indicates the $j$ -th wavepacket that we have sampled. We summarize these procedures in the following Algorithm 1.

Algorithm 1 Frozen Gaussian Sampling algorithm for nonadiabatic quantum dynamics at metal surfaces

1: Initial Sampling: generate

M

i.i.d samples

\left\{z_{0}^{(j)}\right\}_{j=1}^{M}

according to the distribution

\pi

, and then calculate the corresponding

A_{0}^{(0)}\left(0,z_{0}^{(j)}\right)

using equation 3.19.

2: Trajectory evolving: for each

j=1,\cdots,M

, with

z_{0}^{(j)}

and

A_{0}^{(0)}\left(0,z_{0}^{(j)}\right)

, calculate

A_{k^{j}_{[{\nu_{j}}/{2}]:1}}^{(\nu_{j})}

P_{k^{j}_{[{\nu_{j}}/{2}]:1}}^{(\nu_{j})}

Q_{k^{j}_{[{\nu_{j}}/{2}]:1}}^{(\nu_{j})}

and

S_{k^{j}_{[{\nu_{j}}/{2}]:1}}^{(\nu_{j})}

using 3.16, and then calculate

\Theta_{k^{j}_{[{\nu_{j}}/{2}]:1}}^{(\nu_{j})}

using 3.15 and

\Gamma_{k^{j},k^{j}_{[\nu_{j}/2]:1}}^{(\nu_{j})}

using 3.35 and 3.36.

3: Wavefunction constructing: calculate

\frac{1}{M}\sum_{j=1}^{M}\Lambda^{(j)}\left(t,x;\widetilde{z}_{t}^{(j)}\right)

, where the wave packet

\Lambda^{(j)}

is constructed using 3.39.

We can also calculate time-dependent physical observables. We would like to emphasize that this could be done without the reconstruction of wavefunctions. Consider a physical observable $\hat{O}$ , which corresponds to a compact operator on $L^{2}(\mathbb{R}^{m},\mathbb{C}^{N+1})$ , we have

\hat{O}(t)=\frac{1}{M^{2}}\sum_{j,j^{\prime}=1}^{M}\left\langle\Lambda^{(j)}\left(t,x;\widetilde{z}_{t}^{(j)}\right)\right|\hat{O}\left|\Lambda^{(j^{\prime})}\left(t,x;\widetilde{z}_{t}^{(j^{\prime})}\right)\right\rangle.

(3.43)

It should be noted that since $\Lambda^{(j)}\left(t,x,z_{0}^{(j)}\right)$ has the frozen Gaussian form, therefore for a wide class of observables, the integral in the double sum can be evaluated analytically or approximately by Gaussian integral formulas.

We often consider the Gaussian type initial conditions, where the initial value $u_{0,\text{in}}$ is a Gaussian wave packet:

u_{0,\mathrm{in}}(x)=\left(\prod_{j=1}^{m}a_{j}\right)^{\frac{1}{4}}(\pi\varepsilon)^{-\frac{m}{4}}\exp\left(\frac{\mathrm{i}}{\varepsilon}\tilde{p}\cdot(x-\tilde{q})\right)\exp\left(-\sum_{j=1}^{m}\frac{a_{j}}{2\varepsilon}\left(x_{j}-\tilde{q}_{j}\right)^{2}\right).

(3.44)

Here $\tilde{p},\tilde{q}\in\mathbb{R}^{m}$ and $a_{j}(1\leq j\leq m)$ are positive. In this case, we have

	$\displaystyle A_{0}^{(0)}(0,q,p)=2^{m}(\pi\varepsilon)^{\frac{m}{4}}\prod_{j=1}^{m}$	$\displaystyle{\left[\left(\frac{\sqrt{a_{j}}}{1+a_{j}}\right)^{\frac{1}{2}}\exp\left(-\frac{\left(\tilde{p}_{j}-p_{j}\right)^{2}+a_{j}\left(\tilde{q}_{j}-q_{j}\right)^{2}}{2\left(1+a_{j}\right)\varepsilon}\right)\right.}$		(3.45)
		$\displaystyle\left.\times\exp\left(\frac{\mathrm{i}\left(a_{j}\tilde{q}_{j}+q_{j}\right)\left(\tilde{p}_{j}-p_{j}\right)}{\left(1+a_{j}\right)\varepsilon}+\frac{\mathrm{i}\left(p_{j}q_{j}-\tilde{p}_{j}\tilde{q}_{j}\right)}{\varepsilon}\right)\right],$		(3.45)

A suitable initial sampling distribution for Gaussian wave packet is [28, 39]

\pi(q,p)=\frac{A_{0}^{(0)}(0,q,p)}{\int_{\mathbb{R}^{2m}}\mathrm{d}q\mathrm{d}p|A_{0}^{(0)}(0,q,p)|}.

(3.46)

therefore we have

\pi(q,p)=(2\pi\varepsilon)^{-m}\prod_{j=1}^{m}\left(\frac{\sqrt{a_{j}}}{1+a_{j}}\right)\exp\left(-\sum_{j=1}^{m}\frac{\left(\tilde{p}_{j}-p_{j}\right)^{2}+a_{j}\left(\tilde{q}_{j}-q_{j}\right)^{2}}{2\left(1+a_{j}\right)\varepsilon}\right).

(3.47)

It is clear that $\pi$ is a Gaussian distribution and is extremely easy to sample. In other words, we have the initial condition for 3.26:

\rho_{0}(0,q,p)=\pi(q,p),\quad\rho_{k}(0,q,p)=0,\quad k=1,\cdots,N.

(3.48)

For other intial condition, we refer to [39] for a broad discussion on proper choices of the initial sampling measure.

4 Analysis

In this section, we present the numerical analysis of FGS. We focus on the case of Gaussian initial data 3.44.

4.1 Main result

After sampling $M$ wave packets to approximate $U_{\operatorname{FG}}$ , the sampling error should be evaluated using the following weighted $L^{2}$ error:

\mathbf{e}_{N}\left(\left\{\widetilde{z}_{t}^{(j)}\right\}_{j=1}^{M}\right)=\left\|R_{N}\left(\frac{1}{M}\sum_{j=1}^{M}\Lambda^{(j)}\left(t,x;\widetilde{z}_{t}^{(j)}\right)-{\mathbf{U}}_{\text{FG}}(t,x)\right)\right\|_{L^{2}\left(\mathbb{R}^{m},\mathbb{C}^{N+1}\right)},

(4.1)

where $R_{N}$ is the following $(N+1)\times(N+1)$ weight matrix:

R_{N}=\operatorname{diag}(1,h,h,\cdots,h).

(4.2)

We adopt the weighted $L^{2}$ error 4.1 due to the fact that the $(N+1)$ -level Schrödinger equation system is obtained by discretizing the continuum spectrum. In other words, the $h$ in the weighted matrix should be understood as a numerical integration along the energy spectrum $[\mathcal{E}_{a},\mathcal{E}_{b}]$ .

Since FGS is a stochastic method, its numerical error should be evaluated as the following mean square error

\displaystyle\mathrm{Err}_{M}(t,N)=

\displaystyle\left[\mathbb{E}_{\left\{\widetilde{z}_{t}^{(j)}\right\}_{j=1}^{M}}\left(\left(\mathbf{e}_{N}\left(\left\{\widetilde{z}_{t}^{(j)}\right\}_{j=1}^{M}\right)\right)^{2}\right)\right]^{\frac{1}{2}},

(4.3)

where $L^{2}\left(\mathbb{R}^{m},\mathbb{C}^{N+1}\right)$ is the $(N+1)$ -dimensional-valued $L^{2}$ -space on $\mathbb{R}^{m}$ . In fact, since all $z_{0}^{(j)}$ are sampled independently, we have

\displaystyle{\mathrm{Err}}_{M}(t,N)

\displaystyle=\frac{1}{\sqrt{M}}{\mathrm{Err}}_{1}(t,N).

(4.4)

To proceed with the error estimate, we make the following assumption on the potential functions:

Assumption 1.

Each energy surface $U_{k}(q)\in C^{\infty}\left(\mathbb{R}^{m}\right),k\in\{0,1\}$ and satisfies the following subquadratic condition:

\sup_{q\in\mathbb{R}^{m}}\left|\partial_{\alpha}U_{k}(q)\right|\leq C_{E},\quad\forall|\alpha|=2.

(4.5)

If $U_{k}(q)$ depends on $\varepsilon$ , we require the inequality to hold uniformly for all $\varepsilon$ .

Our main numerical analysis result is that:

Theorem 4.1.

For system 1.1 with initial value 3.44, let $M$ be the sample size. With 1, the sampling error ${\mathrm{Err}}_{M}(t,N)$ of the FGS algorithm (Algorithm 1) satisfies the following inequality:

{\mathrm{Err}}_{M}(t,N)\leq\frac{1}{\sqrt{M}}K(t)\prod_{j=1}^{m}\left(\frac{1+a_{j}}{\sqrt{a_{j}}}\right)^{\frac{1}{2}},

(4.6)

where $K(t)$ is independent of both the semiclassical parameter $\varepsilon$ and the number of metal orbitals $N$ .

This theorem indicates that the sample size required by the FGS algorithm does not increase when $N$ gets larger. In fact, recall 3.40 and 3.39, we observe that

$\displaystyle\mathrm{Err}_{1}(t,N)$	$\displaystyle=\left(\mathbb{E}_{\widetilde{z}_{t}}\left(\left\\|R_{N}\left(\Lambda(t,x;{\widetilde{z}_{t}})-{\mathbf{U}}_{\text{FG}}(t,x)\right)\right\\|^{2}_{L^{2}\left(\mathbb{R}^{m},\mathbb{C}^{N+1}\right)}\right)\right)^{\frac{1}{2}}$	(4.7)
	$\displaystyle\leq\left(\mathbb{E}_{\widetilde{z}_{t}}\left(\left\\|R_{N}\left(\Lambda(t,x;{\widetilde{z}_{t}})\right)\right\\|^{2}_{L^{2}\left(\mathbb{R}^{m},\mathbb{C}^{N+1}\right)}\right)\right)^{\frac{1}{2}}$
	$\displaystyle=\left(\mathbb{E}_{\widetilde{z}_{t}}\left(\left\\|W(t,x;{\widetilde{z}_{t}})\right\\|^{2}_{L^{2}\left(\mathbb{R}^{m}\right)}\right)\right)^{\frac{1}{2}}.$

Here this inequality uses the fact that the variance of a stochastic variable is not bigger than its secondary moment. The value of $\left\|W(t,x;{\widetilde{z}_{t}})\right\|^{2}_{L^{2}\left(\mathbb{R}^{m}\right)}$ is a function of the stochastic trajectory $\widetilde{z}_{t}$ . Let us introduce the following notation:

f({\widetilde{z}_{t}})=\left\|W(t,x,z_{0})\right\|^{2}_{L^{2}\left(\mathbb{R}^{m}\right)}.

(4.8)

Then along with 4.4, it suffices to only prove the following result:

Theorem 4.2.

With the same assumptions in Theorem 4.1, we have the following inequality:

{\mathbb{E}_{\widetilde{z}_{t}}\left(f(\widetilde{z}_{t})\right)}\leq(K(t))^{2}\prod_{j=1}^{m}\left(\frac{1+a_{j}}{\sqrt{a_{j}}}\right),

(4.9)

where $K(t)$ is independent of both the semiclassical parameter $\varepsilon$ and the number of metal orbitals $N$ .

4.2 Proof of Theorem 4.2

We first introduce a lemma that gives an upper bound for $\left|\Gamma_{k,k_{[\nu/2]:1}}^{(\nu)}(t,T_{\nu:1})\right|$ :

Lemma 4.1 (Upper bound for $\Gamma_{k,k_{[\nu/2]:1}}^{(\nu)}(t,T_{\nu:1},z_{0})$ ).

Let $C=\max(1,\mathcal{E}_{b}-\mathcal{E}_{a})$ , $M_{0}=\max_{\mathcal{E}\in[\mathcal{E}_{a},\mathcal{E}_{b}]}\left|V(\mathcal{E})\right|$ . we have

\left|\Gamma_{k,k_{[\nu/2]:1}}^{(\nu)}(t,T_{\nu:1},z_{0})\right|\leq\exp({CM_{0}t}),

(4.10)

for any $\nu\geq 0$ , any $k,k_{[\nu/2]:1}\in\{0,1,\cdots,N\}$ , any $z_{0}\in\mathbb{R}^{2m}$ and any $0<t_{1}<\cdots<t_{\nu}<t$ .

Proof.

From 3.21, we can see that $\left|\tau_{k,k_{[\nu/2]:1}}^{(\nu)}\right|\leq M_{0}$ and according to 3.27, $\left|\widetilde{\lambda}_{k_{n:1}}^{(2n+1)}\left(t,T_{2n:1},z_{0}\right)\right|\leq(\mathcal{E}_{b}-\mathcal{E}_{a})M_{0}$ . Therefore, $\left|\tau_{k,k_{[\nu/2]:1}}^{(\nu)}\right|\leq CM_{0}$ , $\left|\widetilde{\lambda}_{k_{n:1}}^{(2n+1)}\right|\leq CM_{0}$ . For $\nu$ being odd, according to 3.36, we have

		$\displaystyle\left\|\Gamma_{k,k_{n:1}}^{(2n+1)}(t,T_{2n+1:1},z_{0})\right\|=\prod_{j=0}^{n}\left(\exp\left(\int_{t_{2j}}^{t_{2j+1}}\widetilde{\lambda}_{k_{n:1}}^{(2n+1)}\left(s_{2j+1},T_{2n:1},z_{0}\right)\mathrm{d}s_{2j+1}\right)\right)$		(4.11)
		$\displaystyle\times\exp\left(\int_{t_{2n+1}}^{t}\tau_{0,k_{n+1:1}}^{(2n+2)}\mathrm{d}s\right)\times\prod_{j=1}^{n}\exp\left({\left(\int_{t_{2j-1}}^{t_{2j}}\left\|\tau_{0,k_{j:1}}^{(2j)}\left(s_{2j},T_{2j-1:1}\right)\right\|\mathrm{d}s_{2j}\right)}\right)$
		$\displaystyle\leq\prod_{j=0}^{n}\exp\left(\int_{t_{2j}}^{t_{2j+1}}CM_{0}\mathrm{d}s_{2j+1}\right)\exp\left(\int_{t_{2n+1}}^{t}CM_{0}\mathrm{d}s\right)\prod_{j=1}^{n}\exp{\left(\int_{t_{2j-1}}^{t_{2j}}CM_{0}\mathrm{d}s_{2j}\right)}$
		$\displaystyle\leq\exp(CM_{0}t).$

For $\nu$ being even, based on 3.35, we can prove the same result in a similar way as in 4.11. ∎

Using lemma 4.1, we can give an upper bound of $\mathbb{E}_{\widetilde{z}_{t}}\left(f(\widetilde{z}_{t})\right)$ :

Lemma 4.2.

With the same assumptions in Theorem 4.1, we have the following inequality:

\mathbb{E}_{\widetilde{z}_{t}}\left(f(\widetilde{z}_{t})\right)\leq\frac{\exp(2CM_{0}t)}{(2\pi\varepsilon)^{\frac{5}{2}m}}\int\mathrm{d}z_{0}\frac{1}{\pi({z_{0}})}\mathbb{E}_{\widetilde{z}_{t}|z_{0}}\left(\left|A_{k,k_{[{\nu}/{2}]:1}}^{(\nu)}\left(t,T_{\nu:1},z_{0}\right)\right|^{2}\right),

(4.12)

where $C,M_{0}$ is defined the same way as in lemma 4.1.

Proof.

With 3.39 and lemma 4.1, we have

	$\displaystyle\left\|W(t,x,z_{0})\right\|^{2}$	$\displaystyle=\frac{1}{(2\pi\varepsilon)^{3m}}\frac{\left\|A_{k,k_{[{\nu}/{2}]:1}}^{(\nu)}\right\|^{2}}{\pi(z_{0})^{2}}\exp\left(-\frac{\left\|x-Q_{k_{[{\nu}/{2}]:1}}^{(\nu)}\right\|^{2}}{\varepsilon}\right)\left\|\Gamma_{k,k_{[\nu/2]:1}}^{(\nu)}\right\|^{2}$		(4.13)
		$\displaystyle\leq\frac{1}{(2\pi\varepsilon)^{3m}}\frac{\left\|A_{k,k_{[{\nu}/{2}]:1}}^{(\nu)}\right\|^{2}}{\pi(z_{0})^{2}}\exp\left(-\frac{\left\|x-Q_{k_{[{\nu}/{2}]:1}}^{(\nu)}\right\|^{2}}{\varepsilon}\right)\exp(2CM_{0}t).$		(4.13)

Since $\int_{\mathbb{R}^{m}}\exp(-{|x-Q_{k_{[{\nu}/{2}]:1}}^{(\nu)}|^{2}}/{\varepsilon})\mathrm{d}x=(2\pi\varepsilon)^{m/2}$ , we have

$\displaystyle\mathbb{E}_{\widetilde{z}_{t}}\left(f(\widetilde{z}_{t})\right)$	$\displaystyle=\mathbb{E}_{\widetilde{z}_{t}}\left(\int_{\mathbb{R}^{m}}\left\|W(t,x,z_{0})\right\|^{2}\mathrm{d}x\right)$	(4.14)
	$\displaystyle\leq\frac{\exp(2CM_{0}t)}{(2\pi\varepsilon)^{5m/2}}\mathbb{E}_{\widetilde{z}_{t}}\left(\frac{\left\|A_{k,k_{[{\nu}/{2}]:1}}^{(\nu)}\right\|^{2}}{\pi(z_{0})^{2}}\right)$
	$\displaystyle=\frac{\exp(2CM_{0}t)}{(2\pi\varepsilon)^{5m/2}}\mathbb{E}_{z_{0}\sim\pi}\left(\frac{1}{\pi(z_{0})^{2}}\mathbb{E}_{\widetilde{z}_{t}\|z_{0}}\left({\left\|A_{k,k_{[{\nu}/{2}]:1}}^{(\nu)}\right\|^{2}}\right)\right)$
	$\displaystyle=\frac{\exp(2CM_{0}t)}{(2\pi\varepsilon)^{5m/2}}\int\mathrm{d}z_{0}\frac{1}{\pi(z_{0})}\mathbb{E}_{\widetilde{z}_{t}\|z_{0}}\left({\left\|A_{k,k_{[{\nu}/{2}]:1}}^{(\nu)}\right\|^{2}}\right).$

∎

Now, with the following lemma, we can finally prove Theorem 4.2.

Lemma 4.3.

Given time $t$ , with 1, for any $\nu\geq 0$ , any $k,k_{[\nu/2]:1}\in\{0,1,\cdots,N\}$ , any $z_{0}\in\mathbb{R}^{2m}$ and any $0<t_{1}<\cdots<t_{\nu}<t$ , there exist $G(t)$ such that

\frac{\left|A_{k,k_{[{\nu}/{2}]:1}}^{(\nu)}\left(t,T_{n:1},z_{0}\right)\right|}{\left|A_{0}^{(0)}(0,z_{0})\right|}\leq G(t),

(4.15)

where $G(t)$ is independent of $N$ and $\varepsilon$ .

The proof of lemma 4.3 is described later. Let us first demonstrate how to use lemma 4.3 to prove Theorem 4.2.

Proof of Theorem 4.2 using lemma 4.3.

With lemma 4.2 and lemma 4.3, we have:

	$\displaystyle\mathbb{E}_{\widetilde{z}_{t}}\left(f(\widetilde{z}_{t})\right)$	$\displaystyle\leq\frac{\exp(2CM_{0}t)}{(2\pi\varepsilon)^{5m/2}}\int\mathrm{d}z_{0}\frac{1}{\pi(z_{0})}\mathbb{E}_{\widetilde{z}_{t}\|z_{0}}\left({\left\|A_{k,k_{[{\nu}/{2}]:1}}^{(\nu)}\right\|^{2}}\right)$		(4.16)
		$\displaystyle\leq\frac{\exp(2CM_{0}t)G(t)^{2}}{(2\pi\varepsilon)^{5m/2}}\int\mathrm{d}z_{0}\frac{\left\|A_{0}^{(0)}(0,z_{0})\right\|^{2}}{\pi(z_{0})}.$		(4.16)

Recall that $\pi\left(z_{0}\right)=\frac{\left|A_{0}^{(0)}\left(0,z_{0}\right)\right|}{\int\left|A_{0}^{(0)}\left(0,z_{0}\right)\right|\mathrm{d}z_{0}}$ , let $K(t)=\exp(CM_{0}t)G(t)$ , then

\displaystyle\mathbb{E}_{\widetilde{z}_{t}}\left(f(\widetilde{z}_{t})\right)

\displaystyle\leq\frac{K(t)^{2}}{(2\pi\varepsilon)^{5m/2}}\int\mathrm{d}z_{0}\frac{\left|A_{0}^{(0)}(0,z_{0})\right|^{2}}{\pi(z_{0})}=\frac{\left(\int\left|A_{0}^{(0)}\left(0,z_{0}\right)\right|\mathrm{d}z_{0}\right)^{2}}{(2\pi\varepsilon)^{5m/2}}.

(4.17)

Since

\int_{\mathbb{R}^{2m}}\left|A_{0}^{(0)}\left(0,z_{0}\right)\right|\mathrm{d}z_{0}=2^{2m}(\pi\varepsilon)^{\frac{5m}{4}}\prod_{j=1}^{m}\left(\frac{1+a_{j}}{\sqrt{a_{j}}}\right)^{\frac{1}{2}},

then we have

\mathbb{E}_{\widetilde{z}_{t}}\left(f(\widetilde{z}_{t})\right)\leq 2^{\frac{3m}{2}}\left(K(t)\right)^{2}\prod_{j=1}^{m}\left(\frac{1+a_{j}}{\sqrt{a_{j}}}\right).

(4.18)

∎

Note that using 4.4 and 4.7, Theorem 4.2 immediately yields Theorem 4.1.

4.3 Proof of lemma 4.3

Now we only need to provide a proof for lemma 4.3.

Note that the dynamics of $Q_{k,k_{[{\nu}/{2}]:1}}^{(\nu)}\left(t,T_{\nu:1},z_{0}\right)$ and $P_{k,k_{[{\nu}/{2}]:1}}^{(\nu)}\left(t,T_{\nu:1},z_{0}\right)$ doesn’t depend on the value of $k$ and $k_{[{\nu}/{2}]:1}$ . This can be seen in 3.22 where we take advantage of the fact that $U_{k}(x)$ only differs from each other by a constant for $k=1,\cdots,N$ . In other words, the Hamiltonian flow only depends on the hopping time $T_{\nu:1}=\left\{t_{\nu},\ldots,t_{1}\right\}$ . We denote the map on the phase space from initial time 0 to time $s$ by $\kappa_{s,T_{\nu:1}}$ ( $s\in[0,t]$ ):

	$\displaystyle\kappa_{s,T_{\nu:1}}:$	$\displaystyle\mathbb{R}^{2m}\rightarrow\mathbb{R}^{2m}$		(4.19)
		$\displaystyle(q,p)\longmapsto\left(Q^{\kappa_{s,T_{\nu:1}}}(q,p),P^{\kappa_{s,T_{\nu:1}}}(q,p)\right),$		(4.19)

such that

\begin{array}[]{l}\left(Q^{\kappa_{s,T_{j:1}}}(q,p),P^{\kappa_{s,T_{j:1}}(q,p)}\right)\\ \quad=\left\{\begin{array}[]{ll}\left(Q^{(0)}(t,q,p),P^{(0)}(t,q,p)\right),&t\leq t_{1}\\ \left(Q^{(j)}\left(t,T_{j:1},q,p\right),P^{(j)}\left(t,T_{i:1},q,p\right)\right),&t\in\left[t_{j},t_{j+1}\right],j\in\{1,\ldots,\nu\}\\ \left(Q^{(\nu)}\left(t,T_{\nu:1},q,p\right),P^{(\nu)}\left(t,T_{\nu:1},q,p\right)\right),&t\geq t_{\nu}\end{array}\right..\end{array}

(4.20)

Here we omit the subscript $k,k_{[{\nu}/{2}]:1}$ for $Q^{(\nu)}$ and $P^{(\nu)}$ .

For a transformation $\kappa$ on the phase space $\mathbb{R}^{2m}$ , we can define its Jacobian matrix as

J^{\kappa}(q,p)=\left(\begin{array}[]{ll}\left(\partial_{q}Q^{\kappa}\right)^{T}(q,p)&\left(\partial_{p}Q^{\kappa}\right)^{T}(q,p)\\ \left(\partial_{q}P^{\kappa}\right)^{T}(q,p)&\left(\partial_{p}P^{\kappa}\right)^{T}(q,p)\end{array}\right).

(4.21)

We also define

Z^{\kappa}(q,p)=\partial_{z}\left(Q^{\kappa}(q,p)+iP^{\kappa}(q,p)\right),\quad\partial_{z}=\partial_{q}-i\partial_{p}.

(4.22)

Let us quote the following lemma from [28, 33] that states the property of the Hamiltonian flow:

Lemma 4.4.

Given $t>0$ , with 1, for any $\nu$ and any $T_{\nu:1}$ , the asscociated map $\kappa_{t,T_{\nu:1}}$ and $J^{\kappa_{t,T_{\nu:1}}}(q,p)$ , $Z^{\kappa_{t,T_{\nu:1}}}(q,p)$ has the follow properties:

$\kappa_{t,T_{\nu:1}}$ is a canonical transformation, i.e.

\left(J^{\kappa_{t,T_{\nu:1}}}\right)^{T}\left(\begin{array}[]{cc}0&I_{m}\\ -I_{m}&0\end{array}\right)J^{\kappa_{t,T_{\nu:1}}}=\left(\begin{array}[]{cc}0&I_{m}\\ -I_{m}&0\end{array}\right),

(4.23)

Here $I_{m}$ is the $m\times m$ identity matrix.

For any $k\in\mathbb{N}$ , there exists a constant $C_{k}$ such that

\sup_{(q,p)\in\mathbb{R}^{2m}}\max_{\left|\alpha_{p}\right|+\left|\alpha_{q}\right|\leq k}\left|\partial_{q}^{\alpha_{q}}\partial_{p}^{\alpha_{p}}\left[J^{\kappa_{t,T_{\nu:1}}}(q,p)\right]\right|\leq C_{k}(t),

(4.24)

$C_{k}(t)$ is dependent on $U_{0}$ and $U_{1}$ .

$Z^{\kappa_{t,T_{\nu:1}}}$ is an invertible $m\times m$ matrix, and for any $k\in\mathbb{N}$ , there exists a constant $C_{k}$ such that

\sup_{(q,p)\in K}\max_{\left|\alpha_{p}\right|+\left|\alpha_{q}\right|\leq k}\left|\partial_{q}^{\alpha_{q}}\partial_{p}^{\alpha_{p}}\left[\left(Z^{\kappa_{t,T_{j:1}}}(q,p)\right)^{-1}\right]\right|\leq C_{k}(t),

(4.25)

$C_{k}(t)$ is dependent on $U_{0}$ and $U_{1}$ .

Since $\kappa$ doesn’t depend on the value of $k_{[{\nu}/{2}]:1}$ , therefore the value of $C_{k}(t)$ in lemma 4.4 doesn’t depend on $N$ at all.

Based on this lemma we are ready to give a proof for lemma 4.3.

Proof of lemma 4.3.

Recall that on $[t_{\nu},t_{\nu+1}]$ , we have

\frac{\mathrm{d}}{\mathrm{d}t}A_{k,k_{[{\nu}/{2}]:1}}^{(\nu)}=\frac{1}{2}A_{k,k_{[{\nu}/{2}]:1}}^{(\nu)}\operatorname{tr}\left(\left(Z^{(\nu)}\right)^{-1}\left(\partial_{z}P_{k_{[{\nu}/{2}]:1}}^{(\nu)}-i\partial_{z}Q_{k_{[{\nu}/{2}]:1}}^{(\nu)}\nabla_{Q}^{2}\widetilde{U}\right)\right),

where

\widetilde{U}(t,Q)=\left\{\begin{aligned} U_{0}(Q),\quad&t\in[t_{2j},t_{2j+1}]\\ U_{1}(Q)+\mathcal{E}_{k_{j}},\quad&t\in[t_{2j-1},t_{2j}]\end{aligned}\right..

Therefore on $[t_{\nu},t_{\nu+1}]$ , we have

\frac{\mathrm{d}}{\mathrm{d}t}\left|A_{k,k_{[{\nu}/{2}]:1}}^{(\nu)}\right|\leq\frac{1}{2}\left|A_{k,k_{[{\nu}/{2}]:1}}^{(\nu)}\right|\left|\operatorname{tr}\left(\left(Z^{(\nu)}\right)^{-1}\left(\partial_{z}P_{k_{[{\nu}/{2}]:1}}^{(\nu)}-i\partial_{z}Q_{k_{[{\nu}/{2}]:1}}^{(\nu)}\nabla_{Q}^{2}\widetilde{U}\right)\right)\right|.

(4.26)

According to 1 and lemma 4.4, the right hand side could be uniformly bounded for $[0,t]$ by a constant $\gamma(t)$ which doesn’t depend on $N$ and $\varepsilon$ . Then according to Gronwall’s inequality, let $G(t)=\exp(\gamma(t))$ , we have

\frac{\mathrm{d}}{\mathrm{d}t}\left|A_{k,k_{[{\nu}/{2}]:1}}^{(\nu)}\right|\leq\frac{1}{2}\left|A_{k,k_{[{\nu}/{2}]:1}}^{(\nu)}\right|\gamma(t)\Rightarrow\frac{\left|A_{k,k_{[{\nu}/{2}]:1}}^{(\nu)}\left(t,T_{n:1},z_{0}\right)\right|}{\left|A_{0}^{(0)}(0,z_{0})\right|}\leq G(t).

∎

4.4 Discussions

With the above analysis, we have established that to reach a certain accuracy, the sample size required by the FGS algorithm depends neither on the number of metal orbitals $N$ nor on the semiclassical parameter $\varepsilon$ . This is the most important advantage of the FGS algorithm. We will also numerically validate this property in section 5.

The numerical analysis of the surface hopping ansatz is inspired by the numerical analysis of one-level Frozen Gaussian Sampling [39]. However, in one-level system, there is no surface hopping, therefore the trajectory becomes deterministic if its initial value has been specified. As a result, the error control of the propagating scheme is much simpler.

Our result, however, indicates that the hopping between energy surfaces doesn’t induce additional sampling errors. To the best of our knowledge, this is the first rigorous mathematical analysis that prove this particular feature of surface hopping algorithms, which gives an explanation to why surface hopping methods have been found to be efficient even in the semi-classical limit.

Furthermore, we have proved the above result in the context of surface hopping at molecule-metal surface. A similar analysis could be applied to surface hopping on two band systems (or any finite band systems) [28, 40], where the authors have already observed numerically that the computational cost is essentially independent of the semiclassical parameter $\varepsilon$ .

5 Numerical Results

In this section, we present numerical experiments to verify the accuracy and convergence properties of our method. We also numerically explore the Anderson-Holstein model to show that our method can capture the information of the quantum observables that classical trajectories fail to capture.

In the following experiments, we solve equation system 1.1 with the following initial value:

\left\{\begin{aligned} &u_{0}(0,x)=(\pi\varepsilon)^{-1/4}\exp\left(-\frac{(x-q_{0})^{2}}{2\varepsilon}\right)\exp\left(\frac{\mathrm{i}p_{0}(x-q_{0})}{\varepsilon}\right),\\ &u_{k}(0,x)=0,\quad k=1,2,\cdots,N.\end{aligned}\right.

(5.1)

The initial condition corresponds to a wave packet localized at $q=q_{0}$ with momentum $p=p_{0}$ in the molecule orbital. We numerically solve all the time-evolution ODEs 3.16 using the fourth order Runge-Kutta scheme with the time step $\Delta t=0.01$ . Without loss of generality, we take $\mathcal{E}_{a}=0,\mathcal{E}_{b}=1$ unless specified otherwise.

5.1 Accuracy and convergence test

The goals of this subsection are two-fold: on the one hand, we aim to verify the accuracy of the FGS algorithm, using the numerical solution obtained by time-splitting spectral method [23] as the reference solution. On the other hand, we demonstrate the efficiency of the FGS algorithm, by showing that the sample size required to reach a certain error threshold is independent of both the number of metal orbitals $N$ and the semi-classical parameter $\varepsilon$ . Besides, the half order convergence of the FGS algorithm with respect to the sample size is observed.

Let us first consider the following example where $U_{0}$ and $U_{1}$ are harmonic potentials. This is the scenario that are widely used in the modeling of nonadiabatc dynamics at metal surfaces [36, 20, 16].

Example 1 (harmonic oscillators). Let

		$\displaystyle U_{0}(x)=\frac{x^{2}}{2},\quad U_{1}(x)=\frac{x^{2}}{2}+0.1x+0.1,$
		$\displaystyle\quad q_{0}=-1.5,\quad p_{0}=2,\quad V(\mathcal{E},x)=1.$

Let us first demonstrate the accuracy of the FGS algorithm. In Figure 4, with $M=20000$ sampled wave packets, we present the real part of the reconstructed wave functions $u_{0}(t,x)$ at the molecule orbital, and $u_{1}(t,x)$ at one of the metal orbitals, with $\varepsilon=\frac{1}{64}$ , $N=16$ and $t=1$ . We can see that the reconstructed wavefunction has a very high precision.

Now let us further investigate the sampling error of the FGS algorithm and numerically validate that it does not depend on $\varepsilon$ and $N$ . We compare the $L^{2}$ sampling errors of FGS algorithm under two sets of scenarios: (1) fix a large enough number of orbitals $N=256$ , let $\varepsilon=1/32,1/64,1/128$ ; (2) fix a small enough semiclassical parameter $\varepsilon=\frac{1}{64}$ , let $N=16,64,256$ . The sampling errors of these two comparisons are listed in table 1. We can see that to reach the same precision, the number of wave packets that we need in the FGS algorithm is independent of both $\varepsilon$ and $N$ . Furthermore, a log-scale plot of the sampling error is available in Figure 5, which verifies that the convergence order of our stochastic method is indeed $\frac{1}{2}$ .

Example 1 ( $N=256$ )	$\varepsilon=1/32$	$\varepsilon=1/64$	$\varepsilon=1/128$
$M=3200$	1.64e-01	1.52e-01	1.55e-01
$M=6400$	1.11e-01	1.15e-01	1.11e-01
$M=12800$	8.00e-02	7.75e-02	7.59e-02
$M=25600$	5.55e-02	5.52e-02	5.38e-02
$M=51200$	4.08e-02	4.05e-02	4.00e-02
Example 1 ( $\varepsilon=1/64$ )	$N=16$	$N=64$	$N=256$
$M=3200$	1.53e-01	1.59e-01	1.52e-01
$M=6400$	1.10e-01	1.10e-01	1.15e-01
$M=12800$	8.39e-02	8.11e-02	7.75e-02
$M=25600$	5.39e-02	5.54e-02	5.52e-02
$M=51200$	4.01e-02	4.00e-02	4.05e-02

Table 1: (Example 1) The

L^{2}

errors of wave functions (at

t=1.5

) for

\varepsilon=\frac{1}{32},\frac{1}{64},\frac{1}{128}

versus various sample size

M

with fixed

N=256

(top) and for

N=16,64,256

versus various sample size

M

with fixed

\varepsilon=\frac{1}{64}

(bottom).

Next, we move on to more complicated examples. Let us consider the extended coupling with reflection (example 2) and the inhomogeneous transition (example 3). The extended coupling with reflection has been one of the most challenging test cases for surface hopping algorithm since Tully’s pioneering work [27], while the inhomogeneous transition is widely used to incorporate generic coupling between molecule and metal orbitals [16, 20].

Example 2 (extended coupling with reflection) Let

	$\displaystyle U_{0}(x)=\arctan(10x)+\frac{\pi}{2},\quad U_{1}(x)=-U_{0}(x),$
	$\displaystyle q_{0}=-1.5,\quad p_{0}=2,\quad V(\mathcal{E},x)=1.$

Example 3 (inhomogeneous transition) Let

	$\displaystyle U_{0}(x)=0,\quad U_{1}(x)=-0.1\exp\left(-0.28x^{2}\right)+0.05,$
	$\displaystyle q_{0}=-2.5,\quad p_{0}=2,\quad V(\mathcal{E},x)=\exp\left(-0.06x^{2}\right).$

The real part of the reconstructed wave function $u_{0}(t,x)$ of these two examples is presented in Figure 6. In example 2, the wave packet is initialized with the center placed at $q_{0}=-1.5$ and is set to travel to the right, and gets reflected back to the left. In example 3, the jumping intensity peaks around the origin and decays rapidly away from the center. These features make the computation much more difficult than that in example 1. Nonetheless, we can still see that the reconstructed wave functions in Figure 6 have pretty high precision.

Let us examine the convergence property of FGS algorithm once again using example 2, similar as what we did in example 1. The comparison of sampling errors is listed in table 2. We can still see that the sampling error is independent of both $N$ and $\varepsilon$ . A log-scale plot of the sampling error Figure 7 demonstrates its half-order convergence with respect to the number of sampled wavepackets $M$ .

Example 2 ( $N=256$ )	$\varepsilon=1/32$	$\varepsilon=1/64$	$\varepsilon=1/128$
$M=3200$	2.70e00	7.50e-01	9.82e-01
$M=6400$	2.20e-01	1.96e-01	1.87e-01
$M=12800$	1.38e-01	1.35e-01	1.27e-01
$M=25600$	1.00e-01	9.79e-02	8.92e-02
$M=51200$	7.21e-02	6.81e-02	6.41e-02
Example 2 ( $\varepsilon=1/64$ )	$N=16$	$N=64$	$N=256$
$M=3200$	2.66e-01	2.68e-01	7.50e-01
$M=6400$	1.92e-01	1.93e-01	1.96e-01
$M=12800$	1.39e-01	1.37e-01	1.34e-01
$M=25600$	9.46e-02	9.47e-02	9.79e-02
$M=51200$	6.73e-02	6.72e-02	6.81e-02

Table 2: (Example 2) The

L^{2}

error (at

t=1.5

) for various

\varepsilon=\frac{1}{32},\frac{1}{64},\frac{1}{128}

versus various sample size

M

with fixed

N=256

(top) and for various

N=16,64,256

versus various sample size

M

with fixed

\varepsilon=\frac{1}{64}

(bottom) .

5.2 Numerical explorations

With the FGS algorithm, we are now ready to explore further on the nonadiabatic dynamics at metal surfaces. On the one hand, we would like to emphasize that our method is superior to the classical trajectories surface hopping methods by showing that the quantum observables of the nonadiabatic dynamics can not be fully captured by the ensemble average of classical trajectories. On the other hand, we use FGS algorithm to explore the physics of metal surfaces, such as the transition rates with different energy gaps and the finite temperature effect.

5.2.1 Quantum observables v.s. classical ensemble average

Let $\hat{O}\in L^{2}(\mathbb{R}^{m},\mathbb{C}^{N+1})$ be the quantum observable associated with the classical quantity $O(q,p)$ . Based on the FGS algorithm, the expectation value of the quantum observable is

\langle\hat{O}\rangle_{\text{q}}(t;\varepsilon)=\int\overline{u}_{0}(t,x)\hat{O}u_{0}(t,x)+\frac{1}{N}\left(\sum_{k=1}^{N}\overline{u}_{k}(t,x)\hat{O}u_{k}(t,x)\right)\mathrm{d}x,

(5.2)

However, using the classical trajectories sampled from $\widetilde{z}_{t}$ 3.22, we can also calculate the classical ensemble average of $O(q,p)$ :

\langle O\rangle_{\text{c}}(t)=\int O(q,p)\pi_{t}(q,p)\mathrm{d}q\mathrm{d}p,

(5.3)

where $\pi_{t}(q,p)$ is the distribution at time $t$ determined by the initial condition $\pi(q,p)$ and the stochastic evolution 3.22. In practice, we have

\langle O\rangle_{\text{c}}(t)\approx\frac{1}{M}\sum_{i=1}^{M}O(q_{i},p_{i}),

(5.4)

Note that in the FGS algorithm, the classical trajectory does not depend on $\varepsilon$ . In other words, the classical ensemble average is independent of $\varepsilon$ , while the expectation valquantum observable definitely depends on $\varepsilon$ . One might expect that when $\varepsilon\rightarrow 0$ , $\langle\hat{O}\rangle_{\text{q}}(t;\varepsilon)$ would converge to $\langle O\rangle_{\text{c}}(t)$ , which is the motivation of many surface hopping methods. However, we demonstrate below that this is not true.

Let us consider $\hat{O}$ to be the nuclear position $\hat{x}$ . Following 5.2 and 5.3, we can calculate the quantum expectation $\langle\hat{x}\rangle_{\text{q}}(t;\varepsilon)$ ( $\varepsilon$ is taken to be $1/2$ , $1/4$ , $1/8$ , $1/16$ , $1/32$ and $1/64$ ) and the classical average $\langle x\rangle_{\text{c}}(t)$ in example 1. The result is shown in Figure 8.

One can clearly see that as $\varepsilon\rightarrow 0$ , $\langle\hat{x}\rangle_{\text{q}}(t;\varepsilon)$ does not converge to $\langle x\rangle_{\text{c}}(t)$ . This means that the classical trajectories can not capture all the physics of the system and one must rely on a quantum-level method, such as FGS algorithm to capture this information. Formally speaking, the FGS method tracks the phase information of each wave packet, and the superposition of those wave packets yields a subtle cancellation which the ensemble of classical trajectories can not provide.

5.2.2 Electron transfer and finite temperature effect

Now, let us numerically explore the physics of Anderson-Holstein model with more specified setting. Let us first simulate the electron transfer at metal surfaces using example $3$ . The transition rate at $t=0.5$ , calculated by both FGS and the time-splitting scheme, are shown in Figure 9, with various $\mathcal{E}_{a}=0,\frac{1}{32},\frac{1}{16},\frac{1}{8},\frac{1}{4}$ but fixing $\mathcal{E}_{b}-\mathcal{E}_{a}=1$ . We remark that $\mathcal{E}_{a}$ can be interpreted as the energy gap between the molecule orbital and the metal orbitals. As $\mathcal{E}_{a}$ decreases, the energy gap gets smaller, and then it’s easier for electrons to transfer between energy surfaces, and therefore the system has a higher hopping rate. This is indeed reflected in Figure 9.

However, according to 3.22, we can see that different metal orbital energies do not affect the classical trajectories, but only contribute to the quantum phases in the wave function. In other words, the change of transition rate with different energy band gaps are fundamentally rooted in the cancellation of the phase of the wave packet. This is also a truly quantum effect that can not be captured by classical trajectories.

We are also able to efficiently simulate the finite temperature effect at metal surfaces [16], by letting the coupling potential $V(\mathcal{E},x)$ be a Fermi-Dirac distribution:

Example 4.(Finite temperature effect)

		$\displaystyle U_{0}(x)=0,\quad U_{1}(x)=-0.1\exp\left(-0.28x^{2}\right)+0.05,$
		$\displaystyle\quad p_{0}=2,\quad q_{0}=-2.5,\quad V(\mathcal{E},x)=\frac{8}{1+\exp\left(\beta\mathcal{E}\right)}.$

Here $\beta$ is the inverse temperature. For different values of $\beta$ ( $\beta=1,16$ ), we plot the real part of the reconstructed wave functions $u_{0}(t,x)$ of the molecule orbital at $t=0.5$ (with $\varepsilon=1/16$ , $N=16$ ) in Figure 10. We can see that the FGS algorithm can efficiently simulate both low and high temperature scenarios at metal surfaces.

From a physical point of view, a higher temperature (i.e. a smaller $\beta$ ) induces more electron transfer. This is indeed observed in Figure 10: a smaller $\beta$ corresponds to a larger amplitude of the wave function at the molecule orbital, which indicates that there is less electron transfer.

6 Conclusions and Discussions

In this work, we have developed a Frozen Gaussian Sampling (FGS) method for efficiently simulating nonadiabatic quantum dynamics at metal surfaces. The key advantage is that its computational cost is essentially independent of both the semi-classical parameter $\varepsilon$ and the number of metal orbitals $N$ . To our best knowledge, we have provided the first rigorous mathematical proof that surface hopping strategies do not induce additional sampling errors.

The FGS algorithm provides various new possibilities for studying nonadiabatic dynamics at metal surfaces. For example, we are modeling the metal surfaces as a closed quantum system in this work. However, the metal surfaces are often viewed as an open quantum system coupled with multiple baths subject to the Boltzmann distribution [20, 16]. The strategy of the FGS could, in theory, be generalized towards the open system setting. What’s more, nonadiabatic dynamics in the strong-coupling regime is very challenging to simulate. This has been discussed in [41] for the two-level spin-boson model, but it is still unclear whether there exists a generic approach for all strong-coupling models. These are left for future research.

Appendix A Derivation of the method

Recall that

\left\{\begin{aligned} \mathrm{i}\varepsilon\partial_{t}u_{0}(t,x)&=-\frac{\varepsilon^{2}}{2}\Delta u_{0}(t,x)+U_{0}(x)u_{0}(t,x)+\varepsilon\sum_{k=1}^{N}hV\left({\mathcal{E}_{k}},x\right)u_{k}(t,x),\\ \mathrm{i}\varepsilon\partial_{t}u_{k}(t,x)&=-\frac{\varepsilon^{2}}{2}\Delta u_{k}(t,x)+\left(U_{1}(x)+{\mathcal{E}}_{k}\right)u_{k}(t,x)+\varepsilon\overline{V}({\mathcal{E}_{k}},x)u_{0}(t,x),\text{ }k=1:N.\end{aligned}\right.

and our ansatz

	$\displaystyle{u}_{0,\text{FG}}(t,x)$	$\displaystyle=u_{0}^{(0)}(t,x)+\sum_{k_{1}=1}^{N}u_{0,k_{1}}^{(2)}(t,x)+\cdots,$
	$\displaystyle{u}_{k,\text{FG}}(t,x)$	$\displaystyle=\sum_{k=1}^{N}u_{k}^{(1)}(t,x)+\sum_{k_{1}=1}^{N}u_{k,k_{1}}^{(3)}(t,x)+\cdots,\quad k=1,\cdots,N.$

The ansatz for $u_{0}^{(0)}(t,x)$ is

u_{0}^{(0)}(t,x)=\frac{1}{(2\pi\varepsilon)^{3m/2}}\int A_{0}^{(0)}\left(t,z_{0}\right)\exp\left(\frac{\mathrm{i}}{\varepsilon}\Theta_{0}^{(0)}\left(t,z_{0},x\right)\right)\mathrm{d}z_{0}.

(A.1)

The ansatz for $u_{k}^{(1)}(t,x)$ is

	$\displaystyle u_{k}^{(1)}(t,x)=$	$\displaystyle\frac{1}{(2\pi\varepsilon)^{3m/2}}\int\mathrm{d}z_{0}\int_{0}^{t}\mathrm{d}t_{1}$		(A.2)
		$\displaystyle\tau_{k}^{(1)}\left(t_{1},z_{0}\right)\times A_{k}^{(1)}\left(t,t_{1},z_{0}\right)\exp\left(\frac{\mathrm{i}}{\varepsilon}\Theta_{k}^{(1)}\left(t,t_{1},z_{0},x\right)\right).$		(A.2)

We have

		$\displaystyle\mathrm{i}\varepsilon\partial_{t}u_{k}^{(1)}(t,x)+\frac{\varepsilon^{2}}{2}\Delta u_{k}^{(1)}(t,x)-\left(U_{1}(x)+{\mathcal{E}}_{k}\right)u_{k}^{(1)}(t,x)$		(A.3)
	$\displaystyle=$	$\displaystyle\frac{\mathrm{i}\varepsilon}{(2\pi\varepsilon)^{3m/2}}\int\mathrm{d}z_{0}\tau_{k}^{(1)}\left(t,z_{0}\right)\times A_{k}^{(1)}\left(t,t,z_{0}\right)\exp\left(\frac{\mathrm{i}}{\varepsilon}\Theta_{k}^{(1)}\left(t,t,z_{0},x\right)\right).$		(A.3)

With comparison to the equation for $u_{k}$ , we actually want to match the following term:

\displaystyle\mathrm{i}\varepsilon\partial_{t}u_{k}^{(1)}(t,x)+\frac{\varepsilon^{2}}{2}\Delta u_{k}^{(1)}(t,x)-\left(U_{1}(x)+{\mathcal{E}}_{k}\right)u_{k}^{(1)}(t,x)=\varepsilon\overline{V}({\mathcal{E}_{k}},x)u_{0}^{(0)}(t,x),

(A.4)

i.e.

		$\displaystyle\mathrm{i}\int\mathrm{d}z_{0}\tau_{k}^{(1)}\left(t,z_{0}\right)\times A_{k}^{(1)}\left(t,t,z_{0}\right)\exp\left(\frac{\mathrm{i}}{\varepsilon}\Theta_{k}^{(1)}\left(t,t,z_{0},x\right)\right)$		(A.5)
	$\displaystyle=$	$\displaystyle\overline{V}({\mathcal{E}_{k}},x)\int A_{0}^{(0)}\left(t,z_{0}\right)\exp\left(\frac{\mathrm{i}}{\varepsilon}\Theta_{0}^{(0)}\left(t,z_{0},x\right)\right)\mathrm{d}z_{0}.$		(A.5)

Therefore, we let

A_{k}^{(1)}\left(t,t,z_{0}\right)=A_{0}^{(0)}\left(t,z_{0}\right),\quad\Theta_{k}^{(1)}\left(t,t,z_{0}\right)=\Theta_{0}^{(0)}\left(t,z_{0}\right),

and we let

\tau_{k}^{(1)}\left(t,z_{0}\right)=-\mathrm{i}\overline{V}({\mathcal{E}_{k}},Q_{0}^{(0)}),

where we use $\overline{V}({\mathcal{E}_{k}},Q_{0}^{(0)})$ to approximate $\overline{V}({\mathcal{E}_{k}},x)$ . This approximation could be justified using Taylor’s expansion, and it is rigorously proved in [28] that it only introduces an $O(\varepsilon)$ error.

Now let’s turn to the ansatz for

	$\displaystyle u_{0,k_{1}}^{(2)}$	$\displaystyle(t,x)=\frac{1}{(2\pi\varepsilon)^{3m/2}}\int\mathrm{d}z_{0}\int_{0}^{t}\mathrm{~{}d}t_{2}\int_{0}^{t_{2}}\mathrm{~{}d}t_{1}$		(A.6)
		$\displaystyle h\tau_{k_{1}}^{(1)}\left(t_{1},z_{0}\right)\tau_{{k_{1}}}^{(2)}\left(t_{2},t_{1},z_{0}\right)A_{0,k_{1}}^{(2)}\left(t,t_{2},t_{1},z_{0}\right)\exp\left(\frac{\mathrm{i}}{\varepsilon}\Theta_{0,k_{1}}^{(2)}\left(t,t_{2},t_{1},z_{0},x\right).\right)$		(A.6)

Plugging in this ansatz into the equation for $u_{0}$ , we identify the term needs to be matched:

		$\displaystyle\mathrm{i}\varepsilon\int\mathrm{d}z_{0}\int_{0}^{t}\mathrm{~{}d}t_{1}h\tau_{k_{1}}^{(1)}\left(t_{1},z_{0}\right)\tau_{{k_{1}}}^{(2)}\left(t,t_{1},z_{0}\right)A_{0,k_{1}}^{(2)}\left(t,t,t_{1},z_{0}\right)\exp\left(\frac{\mathrm{i}}{\varepsilon}\Theta_{0,k_{1}}^{(2)}\left(t,t,t_{1},z_{0},x\right)\right)$		(A.7)
		$\displaystyle=\varepsilon hV\left({\mathcal{E}_{k_{1}}},x\right)\int\mathrm{d}z_{0}\int_{0}^{t}\mathrm{d}t_{1}\tau_{{k_{1}}}^{(1)}\left(t_{1},z_{0}\right)\times A_{k_{1}}^{(1)}\left(t,t_{1},z_{0}\right)\exp\left(\frac{\mathrm{i}}{\varepsilon}\Theta_{k_{1}}^{(1)}\left(t,t_{1},z_{0},x\right)\right).$		(A.7)

Therefore, we let

A_{0,k_{1}}^{(2)}\left(t,t,t_{1},z_{0}\right)=A_{k_{1}}^{(1)}\left(t,t_{1},z_{0}\right),\quad\Theta_{0,k_{1}}^{(2)}\left(t,t,t_{1},z_{0}\right)=\Theta_{k_{1}}^{(1)}\left(t,t_{1},z_{0}\right),

and we let

\tau_{k_{1}}^{(2)}\left(t,z_{0}\right)=-\mathrm{i}{V}({\mathcal{E}_{k}},Q_{k_{1}}^{(1)}).

where we use ${V}({\mathcal{E}_{k}},Q_{k_{1}}^{(1)})$ to approximate ${V}({\mathcal{E}_{k}},x)$ . This approximation again only introduces an $O(\varepsilon)$ error.

Carrying on this process, we arrive at the dynamics 3.16, the initial value 3.18, 3.19 and 3.20, and the hopping coefficient 3.21. Its mathematical justification could be done by the exact same semiclassical analysis in [28], of which we omit the details here.

Acknowledgement

The work of Z.Z. was supported by the National Key R&D Program of China, Project Number 2020YFA0712000, 2021YFA1001200 and NSFC grant Number 12031013, 12171013. The work of Z.H. was partially supported by the National Science Foundation under Grant No. DMS-1652330. Z.Z. thanks Dr. Wenjie Dou for helpful discussions. L.X. thanks Dr. Hao Wu for his support and encouragement.

References

[1] M. Born, R. Oppenheimer, Zur quantentheorie der molekeln, Annalen der Physik 389 (20) (1927) 457–484.
[2] M. Born, K. Huang, M. Lax, Dynamical theory of crystal lattices, American Journal of Physics 23 (7) (1955) 474–474.
[3] P. Bunker, R. Moss, The breakdown of the Born-Oppenheimer approximation: the effective vibration-rotation Hamiltonian for a diatomic molecule, Molecular Physics 33 (2) (1977) 417–424.
[4] S. Pisana, M. Lazzeri, C. Casiraghi, K. S. Novoselov, A. K. Geim, A. C. Ferrari, F. Mauri, Breakdown of the adiabatic Born–Oppenheimer approximation in graphene, Nature Materials 6 (3) (2007) 198–201.
[5] I. Rahinov, R. Cooper, D. Matsiev, C. Bartels, D. J. Auerbach, A. M. Wodtke, Quantifying the breakdown of the Born-Oppenheimer approximation in surface chemistry, Physical Chemistry Chemical Physics 13 (28) (2011) 12680–12692.
[6] J. C. Tully, Mixed quantum–classical dynamics, Faraday Discussions 110 (1998) 407–419.
[7] R. Kapral, G. Ciccotti, Mixed quantum-classical dynamics, The Journal of Chemical Physics 110 (18) (1999) 8919–8929.
[8] R. Kapral, Progress in the theory of mixed quantum-classical dynamics, Annual Review of Physical Chemistry 57 (2006) 129–157.
[9] A. Abedi, F. Agostini, E. Gross, Mixed quantum-classical dynamics from the exact decomposition of electron-nuclear motion, Europhysics Letters 106 (3) (2014) 33001.
[10] R. Crespo-Otero, M. Barbatti, Recent advances and perspectives on nonadiabatic mixed quantum–classical dynamics, Chemical Reviews 118 (15) (2018) 7026–7068.
[11] D. M. Newns, Self-consistent model of hydrogen chemisorption, Physical Review 178 (3) (1969) 1123–1135.
[12] B. Persson, Applications of surface resistivity to atomic scale friction, to the migration of ”hot” adatoms, and to electrochemistry, The Journal of Chemical Physics 98 (2) (1993) 1659–1672.
[13] X. Luo, B. Jiang, J. I. Juaristi, M. Alducin, H. Guo, Electron-hole pair effects in methane dissociative chemisorption on Ni (111), The Journal of Chemical Physics 145 (4) (2016) 044704.
[14] A. Nitzan, M. A. Ratner, Electron transport in molecular wire junctions, Science 300 (5624) (2003) 1384–1389.
[15] M. Head-Gordon, J. C. Tully, Molecular dynamics with electronic frictions, The Journal of Chemical Physics 103 (23) (1995) 10137–10145.
[16] W. Dou, J. E. Subotnik, Perspective: How to understand electronic friction, The Journal of Chemical Physics 148 (23) (2018) 230901.
[17] C. Lindstrom, X.-Y. Zhu, Photoinduced electron transfer at molecule- metal interfaces, Chemical Reviews 106 (10) (2006) 4281–4300.
[18] P. Whitmore, H. Robota, C. Harris, Mechanisms for electronic energy transfer between molecules and metal surfaces: A comparison of silver and nickel, The Journal of Chemical Physics 77 (3) (1982) 1560–1568.
[19] T. Holstein, Studies of polaron motion: Part I. The molecular-crystal model, Annals of Physics 8 (3) (1959) 325–342.
[20] W. Dou, A. Nitzan, J. E. Subotnik, Surface hopping with a manifold of electronic states. II. Application to the many-body Anderson-Holstein model, The Journal of Chemical Physics 142 (8) (2015) 084110.
[21] S. Jin, P. Markowich, C. Sparber, Mathematical and computational methods for semiclassical Schrödinger equations, Acta Numerica 20 (2011) 121–209.
[22] C. Lasser, C. Lubich, Computing quantum dynamics in the semiclassical regime, Acta Numerica 29 (2020) 229–401.
[23] W. Bao, S. Jin, P. A. Markowich, On time-splitting spectral approximations for the Schrödinger equation in the semiclassical regime, Journal of Computational Physics 175 (2) (2002) 487–524.
[24] C. L. Moss, C. M. Isborn, X. Li, Ehrenfest dynamics with a time-dependent density-functional-theory calculation of lifetimes and resonant widths of charge-transfer states of Li⁺ near an aluminum cluster surface, Physical Review A 80 (2) (2009) 024503.
[25] N. Shenvi, S. Roy, J. C. Tully, Nonadiabatic dynamics at metal surfaces: Independent-electron surface hopping, The Journal of Chemical Physics 130 (17) (2009) 174107.
[26] C. Bartels, R. Cooper, D. J. Auerbach, A. M. Wodtke, Energy transfer at metal surfaces: the need to go beyond the electronic friction picture, Chemical Science 2 (9) (2011) 1647–1655.
[27] J. C. Tully, Molecular dynamics with electronic transitions, The Journal of Chemical Physics 93 (2) (1990) 1061–1071.
[28] J. Lu, Z. Zhou, Frozen Gaussian approximation with surface hopping for mixed quantum-classical dynamics: A mathematical justification of fewest switches surface hopping algorithms, Mathematics of Computation 87 (313) (2018) 2189–2232.
[29] B. Engquist, O. Runborg, Computational high frequency wave propagation, Acta Numerica 12 (2003) 181–266.
[30] G. A. Hagedorn, A. Joye, Landau-Zener transitions through small electronic eigenvalue gaps in the Born-Oppenheimer approximation, Annales de l’I.H.P. Physique Théorique 68 (1) (1998) 85–134.
[31] S. Jin, P. Qi, Z. Zhang, An Eulerian surface hopping method for the Schrödinger equation with conical crossings, Multiscale Modeling & Simulation 9 (1) (2011) 258–281.
[32] S. Jin, H. Wu, X. Yang, Gaussian beam methods for the Schrödinger equation in the semi-classical regime: Lagrangian and Eulerian formulations, Communications in Mathematical Sciences 6 (4) (2008) 995–1020.
[33] J. Lu, X. Yang, Frozen Gaussian approximation for high frequency wave propagation, Communications in Mathematical Sciences 9 (3) (2011) 663–683.
[34] Z. Zhou, G. Russo, The Gaussian wave packets transform for the semi-classical Schrödinger equation with vector potentials, Communications in Computational Physics 26 (2019) 469–505.
[35] B. Miao, G. Russo, Z. Zhou, A novel spectral method for the semi-classical Schrödinger equation based on the Gaussian wave-packet transform, IMA Journal of Numerical Analysis (2022).
[36] Y. Cao, J. Lu, Lindblad equation and its semiclassical limit of the Anderson-Holstein model, Journal of Mathematical Physics 58 (12) (2017) 122105.
[37] Z. Jin, J. E. Subotnik, Nonadiabatic dynamics at metal surfaces: fewest switches surface hopping with electronic relaxation, Journal of Chemical Theory and Computation 17 (2) (2021) 614–626.
[38] T. Swart, V. Rousse, A mathematical justification for the Herman-Kluk propagator, Communications in Mathematical Physics 286 (2) (2009) 725–750.
[39] Y. Xie, Z. Zhou, Frozen Gaussian Sampling: a mesh-free Monte Carlo method for approximating semiclassical Schrödinger equations, arXiv preprint arXiv:2112.05405 (2021).
[40] J. Lu, Z. Zhou, Improved sampling and validation of frozen Gaussian approximation with surface hopping algorithm for nonadiabatic dynamics, The Journal of Chemical Physics 145 (12) (2016) 124109.
[41] Z. Cai, D. Fang, J. Lu, Asymptotic analysis of diabatic surface hopping algorithm in the adiabatic and non-adiabatic limits, arXiv preprint arXiv:2205.02312 (2022).

$\displaystyle\mathrm{Err}_{1}(t,N)$	$\displaystyle=\left(\mathbb{E}_{\widetilde{z}_{t}}\left(\left\\|R_{N}\left(\Lambda(t,x;{\widetilde{z}_{t}})-{\mathbf{U}}_{\text{FG}}(t,x)\right)\right\\|^{2}_{L^{2}\left(\mathbb{R}^{m},\mathbb{C}^{N+1}\right)}\right)\right)^{\frac{1}{2}}$	(4.7)
	$\displaystyle\leq\left(\mathbb{E}_{\widetilde{z}_{t}}\left(\left\\|R_{N}\left(\Lambda(t,x;{\widetilde{z}_{t}})\right)\right\\|^{2}_{L^{2}\left(\mathbb{R}^{m},\mathbb{C}^{N+1}\right)}\right)\right)^{\frac{1}{2}}$
	$\displaystyle=\left(\mathbb{E}_{\widetilde{z}_{t}}\left(\left\\|W(t,x;{\widetilde{z}_{t}})\right\\|^{2}_{L^{2}\left(\mathbb{R}^{m}\right)}\right)\right)^{\frac{1}{2}}.$

	$\displaystyle\left\|W(t,x,z_{0})\right\|^{2}$	$\displaystyle=\frac{1}{(2\pi\varepsilon)^{3m}}\frac{\left\|A_{k,k_{[{\nu}/{2}]:1}}^{(\nu)}\right\|^{2}}{\pi(z_{0})^{2}}\exp\left(-\frac{\left\|x-Q_{k_{[{\nu}/{2}]:1}}^{(\nu)}\right\|^{2}}{\varepsilon}\right)\left\|\Gamma_{k,k_{[\nu/2]:1}}^{(\nu)}\right\|^{2}$		(4.13)
		$\displaystyle\leq\frac{1}{(2\pi\varepsilon)^{3m}}\frac{\left\|A_{k,k_{[{\nu}/{2}]:1}}^{(\nu)}\right\|^{2}}{\pi(z_{0})^{2}}\exp\left(-\frac{\left\|x-Q_{k_{[{\nu}/{2}]:1}}^{(\nu)}\right\|^{2}}{\varepsilon}\right)\exp(2CM_{0}t).$		(4.13)

$\displaystyle\mathbb{E}_{\widetilde{z}_{t}}\left(f(\widetilde{z}_{t})\right)$	$\displaystyle=\mathbb{E}_{\widetilde{z}_{t}}\left(\int_{\mathbb{R}^{m}}\left\|W(t,x,z_{0})\right\|^{2}\mathrm{d}x\right)$	(4.14)
	$\displaystyle\leq\frac{\exp(2CM_{0}t)}{(2\pi\varepsilon)^{5m/2}}\mathbb{E}_{\widetilde{z}_{t}}\left(\frac{\left\|A_{k,k_{[{\nu}/{2}]:1}}^{(\nu)}\right\|^{2}}{\pi(z_{0})^{2}}\right)$
	$\displaystyle=\frac{\exp(2CM_{0}t)}{(2\pi\varepsilon)^{5m/2}}\mathbb{E}_{z_{0}\sim\pi}\left(\frac{1}{\pi(z_{0})^{2}}\mathbb{E}_{\widetilde{z}_{t}\|z_{0}}\left({\left\|A_{k,k_{[{\nu}/{2}]:1}}^{(\nu)}\right\|^{2}}\right)\right)$
	$\displaystyle=\frac{\exp(2CM_{0}t)}{(2\pi\varepsilon)^{5m/2}}\int\mathrm{d}z_{0}\frac{1}{\pi(z_{0})}\mathbb{E}_{\widetilde{z}_{t}\|z_{0}}\left({\left\|A_{k,k_{[{\nu}/{2}]:1}}^{(\nu)}\right\|^{2}}\right).$

Efficient Frozen Gaussian Sampling Algorithms for Nonadiabatic Quantum Dynamics at Metal Surfaces

Abstract

keywords:

1 Introduction

2 Quantum dynamics of nuclei at metal surfaces

3 Algorithm

3.1 Review of Frozen Gaussian representation for one-level system

Theorem 3.1.

Theorem 3.2.

3.2 Integral representation based on the surface hopping ansatz

3.3 Stochastic interpretation

Remark 3.3.

3.4 FGS Algorithm

4 Analysis

4.1 Main result

Assumption 1.

Theorem 4.1.

Theorem 4.2.

4.2 Proof of Theorem 4.2

Lemma 4.1 (Upper bound for Γk,k[ν/2]:1(ν)​(t,Tν:1,z0)\Gamma_{k,k_{[\nu/2]:1}}^{(\nu)}(t,T_{\nu:1},z_{0})).

Proof.

Lemma 4.2.

Proof.

Lemma 4.3.

Proof of Theorem 4.2 using lemma 4.3.

4.3 Proof of lemma 4.3

Lemma 4.4.

Proof of lemma 4.3.

4.4 Discussions

5 Numerical Results

5.1 Accuracy and convergence test

5.2 Numerical explorations

5.2.1 Quantum observables v.s. classical ensemble average

5.2.2 Electron transfer and finite temperature effect

6 Conclusions and Discussions

Appendix A Derivation of the method

Acknowledgement

References

Lemma 4.1 (Upper bound for $\Gamma_{k,k_{[\nu/2]:1}}^{(\nu)}(t,T_{\nu:1},z_{0})$ ).