Open access peer-reviewed chapter

Electronic and Molecular Dynamics by the Quantum Wave Packet Method

Written By

Zhigang Sun

Submitted: 23 October 2015 Reviewed: 28 April 2016 Published: 24 August 2016

DOI: 10.5772/63999

Chapter metrics overview

1,575 Chapter Downloads

View Full Metrics


A solution to the time‐dependent Schrödinger equation is required in a variety of problems in physics and chemistry. In this chapter, recent developments of numerical and theoretical techniques for quantum wave packet methods efficiently describe the dynamics of molecular dynamics, and electronic dynamics induced by ultrashort laser pulses in atoms and molecules will be reviewed, particularly on the development of grid methods and time‐propagation or pseudo‐time evolution methods developed recently. Applications of the quantum wave packet for studying the reactive resonances in F + H2/HD and O + O2 reaction, dissociative chemisorption of water on transition‐metal surfaces, state‐to‐state reaction dynamics, state‐to‐state tetra‐atomic reaction dynamics using transition wave packet method and reactant coordinate method, and electronic dynamics in H2+ and H2 molecules will be presented.


  • wave packet method
  • molecular and electronic dynamics
  • grid methods
  • time propagators
  • molecular dynamics on surface

1. Introduction

Accurate and efficient solution to the Schrödinger equation is important in the fields of atomic physics, molecular physics, and chemical dynamics, which include the dynamics of atoms and molecules in time‐dependent electromagnetic fields, a variety of atomic and/or molecular collision problems, molecular photodissociation, molecular dynamics on a surface, and in describing the behavior of materials subjected to internal and external forces, and the production of solitons and quantum vortices in Bose‐Einstein condensates. In many situations, especially at the nanoscale level, the evolution of the physical system occurs over multiple length and timescales, and a computationally efficient approach is a necessity to make progress. Theoretical modeling of all of these phenomena is important to obtain a fundamental understanding of atomic processes at the nanoscale level and to provide a scientific base for the design and development of nanostructured materials and for the optimal control over a chemical reaction.

The numerical solution of the time‐dependent (TD) Schrödinger equation relies heavily on the discretization of the variables (r, t) in coordinate space, or named as grid method [14]. The usage of numerical grid methods, which have numerically spectral convergence, has become the most important tools to solve the Schrödinger equation and has been well developed in the past decades because the grid methods are convenient for evolving the TD Schrödinger equation with an initial value with iterative approaches. The early attempts to demonstrate the computational viability of various grid methods for molecular dynamics have been largely superseded by applications to specific problems and deeper research into more sophisticated molecular systems. The grid methods, however, for efficiently solving the electronic Schrödinger equation are under active development currently, mainly due to the Coulomb‐attractive and repulsive singularities, which are quite different from the potential energy surface (PES) experienced by the nuclei in molecular dynamics calculations.

Once a grid method has been chosen, a suitable iterative method, or an evolving method for the initial wave function, must be determined accordingly. The efficiency of the evolving methods of the Schrödinger equation with time‐independent (TID) or TD Hamiltonian is crucial for computer modeling of the quantum dynamics. They have been well developed before 1990; among them, the most famous two are the Chebyshev polynomial expansion with a TID Hamiltonian and split operator with a TD or TID Hamiltonian [5,6]. Recently, there is particular interest in developing efficient and accurate propagators with a TD Hamiltonian.

An efficient quantum dynamics code comes from an optimal combination of the time propagator and the spatial grid method. Especially, recent rapid developments in hardware capabilities, as exemplified by the Intel's many integrated cores, Nvidia graphics processing unit cards, have added another dimension to increasing the number of floating‐point operations per second that can be executed.

In the final of this part, we would like to note that, in the past, usually we would like to state that the solutions of the Schrödinger equation with the TID method, that is, matrix diagonalization methods, are much more computer consuming than the solutions with the TD method, that is, iterative methods for an initial value problem. Along with the development of the numerical method, such as the filter method or spectral transform method, this is true only conditionally [4,7,8]. Especially, with the intention of the real Chebyshev wave packet method, there is no distinction between the TID and TD method anymore [4].

In this work, the spatial grid methods and time propagators to solve the Schrödinger equation with TD and TID Hamiltonian will be reviewed, and some of their applications, that is, describing the resonances and isotope effects in molecular reaction dynamics, molecular dynamics on surface, and electronic dynamics in atoms and molecules induced by ultrashort laser pulses, will be presented.


2. The Schrödinger equation

The time‐dependent Schrödinger equation usually can be written as


and for the dynamics of a molecule, the Hamiltonian and wave function are functions of R and r, where r is a collection of the coordinates of electrons and R is a collection of the coordinates of nuclei. In the case where we are interested in the dynamics of nuclei, usually, first the electronic structure theory is applied to construct the potential energy surface (PES) using the Born‐Oppenheimer approximation, then the coordinate r disappears in the equation and we can conveniently carry out molecular dynamics simulations. In the case where we are interested in the electronic dynamics processes, especially which are induced by ultrashort laser pulses, the coordinate r has to be considered explicitly, and usually the coordinate R is fixed since the movement of the nuclei is much more slow.

In any case, in a modern solution to the TD Schrödinger equation, the spatial coordinates of the Hamiltonian have to be discretized to facilitate the wave function propagation, where iterative evaluations of the action of Hamiltonian operator on the wave function are required. Sometimes, the time variable also requires discretization; thus, the wave function propagation is accomplished by a series of advancing the wave function from t to t+Δt, where Δt is the time interval or time step. In the following, recent developments concerning these two aspects will be briefly reviewed. The author notes that the selected contents purely come from personal interest and expresses his apologies to the authors whose intelligent works are not included here.


3. The spatial discretization methods

Discretization of the spatial coordinates is very important in an efficient and accurate quantum dynamics calculation. For a special problem, a careful designed grid method would lead to substantial decrease of the computational effort. Basically, there are three approaches to the discretization of the spatial coordinates of the Hamiltonian: (1) local methods, such as the finite difference methods; (2) global methods, which includes all of the spectral methods and the corresponding grid methods; and (3) spectral element methods. Each has its advantages and limitations.

The local method usually applies with low‐order finite difference method and is simple to use, thus leading to very sparse Hamiltonian matrix. However, they usually have limited accuracy since they converge with both the number of the grid points and the spacing of the grid points. When boundary region is important in the model, we need to impose the boundary conditions carefully, which results in an asymmetric Hamiltonian matrix, and unphysical states may arise. The grid points in finite difference calculations connect with each other only “locally,” that is, connect with the nearby grid points; thus, parallel algorithm can be efficiently applied with MPI.

On the contrary, the global methods usually are derived based on the classical orthogonal polynomial. The grid point connects with all of the other points, necessarily or unnecessarily, which leads to the dense Hamiltonian matrix. However, the global method usually has spectral convergence, and results of machinery accuracy can be obtained. Due to limited accuracy of finite element method, spectral element method is much more popular in a quantum molecular and atomic dynamics calculation [911]. The grid points in an element of the spectral element are determined by the Lobatto‐Chebyshev or Lobatto‐Legendre quadrature. To set up the connection between different elements, many proposals have been put forward. The spectral element usually leads to numerical accuracy between the spectral method and low‐order finite difference method, and results in Hamiltonian matrix sparser than that using the global method. Currently, the spectral element method has popular applications in solving electronic Schrödinger equation induced by ultrashort laser pulses. However, due to the denser grid points around both ends of each element, the spectral range of the Hamiltonian matrix usually is necessarily large; especially, high‐order Lagrange polynomials are applied in each element for obtaining results of high accuracy. This unnecessary huge spectral range results in much difficulty in searching for an efficient matched time propagator.

Besides these three spatial discretization methods, we also have higher order finite difference method, spectral difference method, and the distributed approximation functional (DAF) method [1216]. They usually exhibit spectral convergence, exactly same as the global spectral methods. However, they connect only with the nearby “necessary” grid points, thus leading to denser Hamiltonian matrix than that with low‐order finite difference method, but sparser Hamiltonian matrix than that with the global methods. The grid points distribute evenly in the spectral difference method, but the grid points distribute exactly the same as their precedent Gauss quadrature of the kernel in the DAF method. They avoid the discontinuities between elements in the spectral method, thus an optimal choice in a calculation with smooth interaction potential functions. Since the grid points connect with each other semilocally, methods of this kind are suitable in a parallel calculation also, similar to the spectral element method. The higher order difference or spectral difference method and DAF method have been applied successfully in the molecular dynamics field. Recently, Sun also puts forward DAF method, which is capable of accurately describing the Coulomb singularity.

In the following, brief introductions to several spatial discretization methods, which currently are popular in a quantum dynamics calculation, will be given.

3.1. Discrete variable representation (DVR)

The DVR has been used quite successfully for a variety of problems, including the calculation of the ro‐vibrational spectra of complex molecules and the reactive scattering cross sections of atoms and molecules [17,18]. In a usual variational or spectral calculation, one chooses some, in principle, complete, orthonormal basis set, φ, with truncation to a finite number of terms, N, and then evaluates all of the matrix elements of the Hamiltonian analytically, or by a quadrature technique. If the spectral basis is chosen to be one of the classical orthogonal functions, satisfying the three‐term recursion relationship,


then there exists an associated Gauss quadrature rule with the same number of points and weights, which integrates exactly any integrand of degree (2N - 1) or less. By diagonalizing the tridiagonal matrix built from the recursion coefficients α and β, the points and weights may be obtained. Under these conditions, it is possible to transform from the original spectral basis of N functions, φ, to a new basis of DVR or coordinate functions, called O, as follows:


Here, all of the integrals are performed using the Gauss quadrature rule; there are no approximations as a consequence of the fact that the quadrature rule is exact for all integrands (2N - 1) or less. The Fourier functions or particle in a box, which are not even polynomials, can be verified to exactly satisfy Eq. (2) with a suitably chosen quadrature rule. In all of these cases, there exists a unitary transformation between the original spectral and coordinate basis. In the DVR, the potential V(x) can be written as


since V(x) can be written as


Thus, the matrix elements of the potential operator are equal to their values at the Gauss quadrature points and diagonal. In a wide class of problems, this approximation has been proved to be accurate. As a result, there are significant advantages in transforming to the DVR basis over the original spectral basis in that the complex matrix elements of the potential are diagonal and are simply equal to their values as the quadrature points. This is very convenient in a numerical calculation since it leads to a sparse form of the Hamiltonian matrix.

As the other part of the Hamiltonian, the matrix elements of the kinetic energy operator, while not diagonal in the DVR basis, usually can be evaluated simply and exactly or analytically. Since the kinetic operator part of the Hamiltonian matrix is a separable sum over the particle and coordinate variables, a direct product DVR basis results in a very sparse matrix representation in n‐body problems. This has a number of advantages which can be well exploited in the iterative stage, where the action of the Hamiltonian matrix on the wave function needs to be evaluated many times, such as the time propagation step or the iterative computation in Lanczos method for calculating the eigenstates. In a practical calculation, the basis functions are not necessary to be the classical polynomials or the trigonometric functions. Usually, the eigenstates calculated in a reduced dimensional model are taken as the basis function, which often leads to efficient potential optimized DVR. In the past decades, complicated but efficient sequential truncation techniques and phase optimization methods have been developed for solving ro‐vibrational states of polyatomic molecules [8,19]. Calculations of reactive scattering and ro‐vibrational state of molecules with dimensionality of nine using the DVR technique in a direct product form have been reported.

Usually, the DVR invented in the field of molecular dynamics does not work well with electronic dynamics; especially, the Coulomb singularity was rigorously included. The exception is the Lobatto‐DVR method, which based on the Legendre polynomials but with Gauss‐Lobatto quadrature [20]. However, with suitable phase optimization technique, the trigonometric functions can work well with the Coulomb singularities, even for describing the electrons of diatomic molecule in cylindrical coordinates.

The Lagrange meshes (LM) method is popular in the field of atomic physics [21]. Same as the DVR method, the LM method is a global method, using information from the whole domain of definition of the studied problem, an approximate spectral method taking the form of a collocation method. The LM method is a special case of DVR method, and when classical orthogonal polynomials are adopted as the basis, the LM and the DVR methods are identical. However, the LM method is derived from the Lagrange functions, which make it to be convenient for regularization of the singularity, starting from the basis function. However, usually, the quadrature in a DVR calculation is obtained by diagonalizing the coordinate matrix Xij=ϕi|x|ϕj, and the kinetic matrix in DVR is given by


where O is an orthogonal matrix, which diagonalizes the coordinate matrix, X, defined as


This fact makes DVR be convenient for various generalizations and for the optimization of the DVR. The basic idea is to get rid of the constraints on the choice of the mesh points and associated weights. They are close to the collocation spirit, but they cannot be so easily related to the LM method. It has been proven that these generalizations are capable of giving very accurate results.

The quadrature discretization (QD) method is efficient for calculating low‐bound states of molecules [22]. Same as the DVR and the LM methods, the QD method is a global method. The QD method is very similar to DVR method, but using nonclassical polynomials. In the QD method, a differential operator is replaced by a matrix. The mth derivative  dmdxm is represented by the mth power of a matrix representing  ddx in a basis of nonclassical orthogonal polynomials. Like in the DVR method, a matrix transformation leads to a basis where the potential is represented by a diagonal matrix. In a basis of N polynomials, ϕj(x), orthogonal with respect to a weight function, ω(x), a Gauss quadrature is obtained with the zeros xi of ϕN(xi) and weights ωi, and the Hamiltonian matrix can be written as


where matrix D is obtained by an orthogonal transformation based on O from the matrix elements of  ddx in the polynomial basis,


with the auxiliary potential that reads as


When ω(x) is chosen as V(x), it means


where ψ0(x) is the ground‐state wave function of the potential V(x). Thus, this method is particularly useful for calculating low‐lying bound states.

3.2. Multidomain spectral element method

The Coulomb‐attractive singularities between nuclei and electrons are very difficult to describe accurately. At the same time, usually the range of the radial coordinate in a calculation has to be large, which leads to large Hamiltonian matrix using a global spectral method and results in heavy computational effort. The Lobatto‐DVR, using the Lagrange interpolation polynomials on the Gauss‐Legendre‐Lobatto quadrature, is particularly suitable for describing the Coulomb‐attractive singularity, whose quadrature, however, clusters densely around the two grid ends. Thus, average efficiency of the quadrature is not high. To overcome these two difficulties, the multidomain spectral method was purposed [3,4]. In this method, the whole range of the radial coordinate is divided into M segmentations or finite elements with boundaries,


The wave functions are represented in a basis of functions which are local to each finite element. In this spectral element method (which is also named as finite element [FE]‐DVR), a few of the Lobatto‐DVR functions,


which fulfill


are applied for each of the element. The basis function is then defined explicitly with N quadrature as


where the “bridge” function is used to ensure the communication between two adjacent elements and guarantees the continuity of the wave function. Correspondingly, the kinetic matrix is given by


assuming each element has same N Lagrange functions. Further, tmni is given by


where sk=[2xki(xi1+xi)/(xixi1)] and LN1(s) are the N‐order Legendre polynomials. For m=k=2,3,,N1,


thus, the first and last diagonal values can be found to be


Figure 1 illustrates the structure for the kinetic energy matrix of one dimension. Since the many‐particle kinetic energy operators are the sum of one‐particle operators, the generalization to tensor product basis sets is straightforward. Figure 2 shows the gird points of a typical FE‐DVR grid with 11 basis functions per element and the resulting grid spacing.

Figure 1.

Left: FE‐DVR basis functions for three elements, where each element has four basis functions. The solid and dashed lines correspond to the element and bridge functions, respectively. Right: The structure of the 1D Hamiltonian matrix, corresponding to the basis in the left panel.

Figure 2.

Grid points (left panels) and their spacing (right panels) of a typical FE‐DVR grid (upper panels) and optimized FE‐DVR grid (lower panels) with 11 basis functions per finite element. The oscillatory structure in the grid spacing of upper panels is a consequence of the Gauss‐Lobatto quadrature. However, after optimization, the oscillation is much more gentle.

Even the finite‐element DVR has been proven to work well in many cases; the oscillatory distribution of the grid and the bridge function leads to unphysically large spectral range of the Hamiltonian matrix, as shown in the upper panels of Figure 3. This fact not only leads to low grid efficiency but also makes the iterative propagator, such as short‐time iterative Lanczos (SIL) method, converge slowly [23]. To overcome this difficulty, one can use the variable mapping method or prolate spheroidal wave functions to get a more uniform spatial grid far from the origin and dense spatial grid around the origin, to increase the numerical efficiency [24,25]. In the bottom panels of Figure 2, the gird distribution becomes nearly even after optimization, which improves the efficiency of the grid. At the same time, the spectral range of the Hamiltonian matrix reduces several times against that of the original FE‐DVR.

Figure 3.

Left panels: The low eigen energies of the kinetic matrix grow quadratically of a usual FE‐DVR (upper) and optimized FE‐DVR (bottom). Right panels: For higher energies, the spectrum contains “unphysical” steps due to the division into finite elements, which reduces much of optimized FE‐DVR (bottom right panel).

On the other hand, in order to exploit the sparseness of the Hamiltonian matrix in finite element DVR, the powerful split operator propagator cannot be applied, which even is insensitive to the spectral range of the Hamiltonian. The common choice is the short‐time iterative Lanczos method, which requires a significant amount of computer memory. And, due to the fact that the continuity between elements only is guaranteed with the wave function other than the derivatives, which introduce numerical complications, results with spectral accuracy cannot be obtained.

3.3. Distributed approximating functionals (DAF) method

Even though it is known that for well‐conditioned problems, that is, when infinitely differentiable solutions exist almost everywhere, global methods provide a better accuracy than local methods using information from pieces of the domain, such that the Lobatto‐DVR method is capable of giving results of higher accuracy than finite element DVR method, the drawback is that the Hamiltonian matrix is dense and sometimes requires more computational effort. The higher‐order finite difference (FD) method is a good choice for obtaining results of high accuracy, but possibly with sparsest Hamiltonian matrix, especially with spectral difference method where accelerating weights are applied.

However, in the case of Coulomb potential with singularity, the boundary conditions are very important to the convergence of the FD method. Usually with well‐defined boundary conditions, the resulted Hamiltonian matrix is asymmetric, which leads to unphysical states that may be numerically unstable.

The distributed approximating functional (DAF) method presents a good solution to these problems, which reduce the correlated region around the grid point by a controllable Gaussian weight ϖ(x,σ), but retaining the numerical accuracy by monitoring the convergence [12,2630]. This kind of method uses the fact that the grid is unnecessary to correlate with all of the remaining grids, especially in the case where large spatial grid range has to be applied. The basis functions with specified DVR functions ϕi(x) as the kernel of the DAF may be defined as


where ϖ(x,σ) is defined as


with a controllable parameter σ. In the DAF, the momentum and kinetic operator can be derived from those of Pi,j and Ti,j in the corresponding DVR method directly by the following equations:




The potential operator has exactly the same form as that in the precedent DVR method. Thus, the Hamiltonian matrix for one dimension has similar form as that using higher‐order FD method, as shown in Figure 4.

Figure 4.

The matrix in the DVR method (red and blue elements) and the corresponding DAF method (red elements). The matrix using the DAF is much sparser.

As observed from Eq. (3), the DAF concept can be conveniently applied with many DVR kernels. The usage of DAF and the corresponding DVR method can be switched back and forth conveniently for checking convergence and saving computational effort, since the introduction of the DAF only reduces the computational effort, without affecting the numerical efficiency of its precedent DVR. The DAF method has been successfully applied for solving vibrational bound state, reactive scattering process, and calculating electronic states with explicit Coulomb singularity.

3.4. Sparse grid (SG)

With increasing dimensionality, the computational effort for solving Schrödinger equation increases exponentially using the basis space by a direct product expansion of a one‐dimensional (1D) basis, regardless of how efficient the grid is for one dimension. As we understand, the physically required phase space does not increase so fast with increasing dimensionality in a specified energy range. Thus, many grid points (combinations) from different degrees of freedom in a direct product form are actually useless.

In 1963, Smolyak first introduced a multivariable integral method to overcome this difficulty [31]. This approach constructs a multidimensional multilevel basis by a special truncation of the tensor product expansion of a one‐dimensional multilevel basis, which holds the promise of alleviating the restriction of excessively large numbers of unknowns and solves a set of differential equations using significantly fewer unknown. In the limit of high accuracy, the number of unknowns required by the sparse grid combination technique is independent of the spatial dimensionality of the problem. For example, asymptotically, a spatial 3D problem requires the same order of unknowns as a spatial 1D problem.

Sparse grids for solving a discretized partial differential equation (PDE) were introduced in 1990 by Zenger [32], in order to significantly reduce the number of degrees of freedom, while causing only a marginal increase in the representation error relative to the standard discretization. In 1992, Griebel et al. [33] showed that, for two and three dimensions, the sparse grid complexity and representation error can also be achieved by the so‐called combination technique which is suitable for a parallel computation.

The sparse grid combination technique can be understood as a multivariate extrapolation technique. Instead of solving a set of differential equations on a single grid, solutions are obtained on a number of semi‐coarsened grids. After solving these semi‐coarsened problems, the solutions are combined to obtain a single, more accurate solution. Currently, there is focused interest for solving various PDE using sparse grid technique with different quadratures and in combination with multidomain techniques.

There are many studies on the introduction of the sparse grid method in the literature, and their theoretical details are not presented here. With multidimensional sparse grid and hierarchical bases, vibrational states of a molecule with dimensionality of 64 have been obtained [34], electronic structures of small molecules have been solved [3540], and solution of the time‐dependent Schrödinger equation of model potentials with split operator method, with hierarchical Fourier basis, have been reported [41].

In passing, we note that the coordinate choice is important in an efficient quantum dynamics calculation, especially for a calculation using sparse grid technique where the combination technique can be effectively realized. In spirit, the usage of hierarchical basis is very similar to the high‐dimensional representations (HDMR) and multiconfiguration time‐dependent Hartree (MCTDH) method [42], and general applications of the multidimensional sparse grid and hierarchical basis in the high‐dimensional atomic and molecular dynamics are expected in the near future [43]. Inspecting the prospective development of multidomain technique for solving PDE and the coordinate problem in a state‐to‐state reactive scattering calculation, we might also expect that an exploration of multidomain technique in the field of molecular reaction dynamics would be prospective.


4. Time discretization methods

Besides the importance of the coordinate discretization methods for solving the TD Schrödinger equation, the efficiency of a time propagator is also vital. Especially, it is not surprised that the efficiency of a time propagator is only high for special spatial discretization methods. During the past decades, various time propagators have been developed, which roughly can be classified into two classes: One is short‐time propagator which often is of low order and the other one is long‐time propagator by a polynomial expansion which often is accurate. In the following, only the most useful ones are briefly reviewed.

4.1. Split operator method (second and higher order)

The second‐order split operator method is the most popular one in a quantum dynamics calculation, especially in a molecular quantum dynamics calculation in combination with the DVR method [1]. The second‐order split operator method is realized as




for the TID Hamiltonian or TD Hamiltonian, ignoring the time ordering. With the DVR method or the Fourier basis, where fast Fourier transform (FFT) can be applied, the split operator can be conveniently implemented. Higher order split operator can be numerically much more efficient than the second‐order split operator for evolving the wave function of a diatomic molecule or triatomic reactive scattering [44,45].

4.2. Crank‐Nicolson and spectral transformed Crank‐Nicolson propagator

The Crank‐Nicolson is particularly popular for solving a TD Schrödinger equation in the field of atomic physics. One reason, which is well known, is that the Coulomb singularity is difficult for a spectral method, but the finite difference is good for describing the Coulomb singularity thus often adopted. The finite difference in combination with the Crank‐Nicolson propagator is very efficient as a propagator of TD Schrödinger equation by solving a linear equation. The other one reason is not so explicit and needs a bit more explanation, which comes from the concept of spectral transformed Hamiltonian [46]. As we know, with a TID Hamiltonian H^,






By retaining the first Taylor expansion term of tan(x), we obtain the usual Crank‐Nicolson propagator


On the other hand, we can define a spectral transformed Hamiltonian Θ^ as


Hamiltonian Θ^ has the exactly same eigenfunctions as the original Hamiltonian H^, but with different eigenvalues. The TD Schrödinger equation with Hamiltonian Θ^ can be written as


For this equation, we conveniently have an exact short‐time propagator


Therefore, in the traditional usage of the Crank‐Nicolson propagator, the wave function actually is propagated with the spectral transformed Hamiltonian Θ^. In the case where the time step is small, the energy of a wave function of the Hamiltonian H^ can be effectively approximated by the one of Θ^.

In the calculation involving Coulomb singularity, the spectral range of the Hamiltonian usually is quite huge, which imposes the upper limit for the time step by


As seen from Eq. (4), this spectral transformation reduces the original huge spectral range of Hamiltonian H^ into a very limited one, [π/2,π/2], which alleviates the upper limit of the time step much. This fact explains why the split operator in the exponential form is less accurate than that in the Crank‐Nicolson form, contrary to intuition.

4.3. Short‐time iterative Lanczos method

Another popular short‐time propagator is the short‐time iterative Lanczos (SIL) method [47]. The Lanczos method initially was proposed to calculate the eigenvalues and eigenfunctions of a symmetric Hamiltonian matrix by solving a tridiagonal matrix, which was given by a similar transformation of the original Hamiltonian matrix through the Krylov vectors. Therefore, the fundamental idea of the Lanczos method is to iteratively build a small, orthornormal set of vectors using the Krylov subspace generated from the solution to the time‐dependent Schrödinger equation at the previous time step as the initial vector. This small set of vectors are used to diagonalize the Hamiltonian, which enables the exponentiation of the matrix to be carried out. The SIL method is realized by initializing


and by expressing generally


where the coefficients are given by


The projected subspace representation of the Hamiltonian operator has the tridiagonal form


Then, the propagation can be accomplished by diagonalizing a small matrix Z:


The SIL method requires memory usually one order larger than that required in a calculation using the split operator or Crank‐Nicolson method. However, since SIL does not need to switch between different representations, unlike that using the split operator method, the SIL method can explore the block sparseness of the kinetic operator matrix for one dimension. Several works have been reported on the combination usage of the SIL method and the FE‐DVR method. On the other hand, the iteration number for a convergent result using the SIL has direct ratio to the spectral range of the Hamiltonian. Therefore, even the SIL method takes well the advantage of the sparseness of the Hamiltonian matrix using the FE‐DVR method; the unphysically large spectral range resulted from the finite element method reduces the convergence of the SIL method, which requires more iterations, and thus more memory and more computational effort.

4.4. Symplectic propagator

Symplectic time propagators for time‐dependent Schrödinger equation are of particular interest, especially higher order ones. They can be very accurate and implemented simply [48]. For a typical symplectic propagator, the implementation can be done by separately expressing the wave function into real and imaginary parts,


where both φ and γ are real. Then, the time‐dependent Schrödinger equation can be written as


which can be rewritten as




The evolution operator of Eq. (5) can be written as


which is an orthogonal and symplectic matrix. Direct evaluation of Eq. (6) is expensive. By using composition methods, evaluation of Eq. (6) can be written as




Various coefficients determined by the composition methods have been reported, which have been proven highly accurate for solving a time‐dependent Schrödinger equation to molecular dynamics.

Even symplectic propagator implemented in this form is accurate and straightforward; it is unstable and gives diverged results when the adopted time step is larger than an upper limit. Usually, this upper limit is not large, which is determined by the spectral range of the Hamiltonian. Thus, the numerical efficiency of symplectic propagator, usually, is not high. This is quite different from the split operator in exponential form whose time step is only determined by the numerical error. This fact explains its much less popularities in a practical quantum dynamics calculation.

4.5. Chebyshev and real Chebyshev polynomial expansion with a TID Hamiltonian

The most popular long‐time propagator is Chebyshev polynomial propagator, introduced by Taz and Kosloff [2]. It is accomplished as follows:


where Pn is the n‐order of Chebyshev polynomial. Emax and Emin are the maximal and minimal energies of the Hamiltonian. The expansion coefficient an is given by


where Jn is the n‐order Bessel function of first class. Since the Chebyshev polynomials can be obtained by the following iteration,


the propagation of the wave function by the Chebyshev polynomial expansion can be implemented by an iterative procedure.

For the initial wave function which has only the real part, the time evolution of the wave function can be achieved in the order domain of Chebyshev polynomial, instead of the time domain [49,50]. This can be understood as follows:


With Ψ=φ+iγ, we have




Therefore, the real and imaginary parts of the wave function can be propagated separately. By introducing the spectral transformed Hamiltonian,


Then, in the space of Ξ, the evolution of the real‐wave function φ can be written as


Here, the time step Δt is only an arbitrary parameter; thus, we simply set it as 1.0. Equation (7) then becomes


with φ(1)=H^φ(0); the evolution of the wave function can then be carried out simply by the iterative action of the Hamiltonian on the wave function, according to Eq. (8). It has been proven that the real Chebyshev wave packet method is about two times faster than the Chebyshev polynomial expansion method, but the DVR method for direct reaction is less efficient than the second‐order split operator [51,52].

4.6. Propagator with explicitly TD Hamiltonian

When the Hamiltonian is time‐dependent, such as in a case where there is interaction between laser pulse and molecule or atom, the time ordering in a calculation using short‐time propagator for numerical results of high accuracy has to be considered sometimes. For short‐time propagators, one can supplement with a Euler‐MacLaurin expansion for the time‐integrated nonhomogeneous term or with a Magnus expansion for considering the time ordering [53,54]. Tal‐Ezer and his coworkers recently introduced a new Chebyshev polynomial expansion method for nonhomogenous time‐dependent Schrödinger equation. For more details, one may refer Ref. [55].


5. Boundary conditions

In a solution to TD Schrödinger equation using spatial discretization method, the spatial grid range has to be limited. However, during the propagation time, the wave function usually can extend out of the grid. We have to artificially impose the boundary condition, which is capable of damping the wave function before reaching the end of the grid, but without affecting the interior wave function. A usual choice is by including the absorbing potential in the potential operator or by imposing damping function to the wave function. Smooth exterior complex scaling method was reported to be a more accurate way to impose the boundary conditions in a numerical solution to TD Schrödinger equation, which has been applied in a calculation of electronic dynamics induced by ultrashort laser pulses.


6. Applications

With the pioneer work in the solution to the Schrödinger equation in the field of molecular quantum dynamics, the developments of quantum wave packet now have a general application in the field of atomic physics, molecular physics, molecular reaction dynamics, and the structure simulations of semiconductors or nanomaterials. In the following sections, only those chosen from the author's personal interest are briefly reviewed, which may not be the most representative works.

6.1. Feschbash resonance in the reaction of F + H2/HD (v0 = 0, j0 = 0) and Cl + HD (v0 = 1, j0 = 0)

The reactive resonances discover the quasi‐bound states of the reaction complex with unique clarity, and they do exist. Recognition of the reactive resonances is helpful for understanding how elementary chemical processes at a single quantum‐state level take place. F + H2 and its isotopic analogs are the most classical prototypes [56].

Interest in the F + H2 reaction largely results from Lee's benchmark molecular beam studies [57] and early chemical laser work [58]. In 2000, the work on the F + HD reaction of Dong et al. [59] discovered the existence of reactive resonances in crossed beam experiments, where a resonance‐enhanced step in the excitation function was observed.

Figure 5.

Schematic diagram indicating the resonance‐mediated reaction mechanism for the F + H2 reaction, which has two resonance states trapped in the HF(v’ = 3)‐H vibrational adiabatic potential (VAP)well. The 1D wave functions of the two resonance states are also given. The (003) state is the lowest ground resonance state; the (103) resonance is the first excited resonance state. Calculated van der Waals states for the lower VAPs are also given. OP, overtone pumping; Eb, barrier height; Ec, collision energy.

In 2006, with both theory and experiment exhibiting consistent behavior on XXZ PES, the work of Qiu et al. [60] reported the evidence for the resonances in the F + H2 reaction, as shown in Figure 5. The appearance of this report is enabled by the crossed molecular beam combined with high‐resolution Rydberg state tagging technique, the developments on quantum scattering method development, and ab initio method. The sharp forward peak arising at the collision energy of 0.52 kcal/mol is actually due to the interference between the first two Feshbach resonances. With total angular momentum as zero, there are two resonance states at collision energies of 0.26 and 0.46 kcal/mol. Along with increasing J, the resonance energy will shift to higher collision energy. The three‐dimensional (3D) scattering wave function at the collision energy of 0.26 kcal/mol suggests the existence of three nodes along the H‐F coordinate (which correlates to the HF product) in the HF‐H’ complex with zero node along the reaction coordinate. The projection of the J = 0 scattering wave function at 0.26 kcal/mol to the HF vibrational states suggests that the dominant character of this wave function is HF (v’ = 3), with the outgoing waves mostly on HF (v’ = 2). This implies that the resonance state at 0.26 kcal/mol is the lowest resonance state, (003), which is trapped in the HF (v’ = 3)‐H’ vibrational adiabatic potential (VAP) well.

The 3D scattering wave function for J = 0 at the collision energy of 0.46 kcal/mol suggests the existence of three nodes along the HF coordinate (which correlates to the HF product) in the HF‐H’ complex, with a single node along the reaction coordinate. The projection of the J = 0 scattering wave function at 0.46 kcal/mol to the HF vibrational states indicates that the dominant character in this wave function is HF (v’ = 3), but with the outgoing waves also mostly on HF (v’ = 2). This indicates that the resonance state at 0.46 kcal/mol is the excited reaction resonance state trapped in the HF (v’ = 3)‐H’ VAP well. This resonance state can be assigned to the (103) resonance state, which has one‐quantum vibration along the reaction coordinate, zero‐quantum vibration on the bending motion (or hindered rotation), and three‐quanta vibration along the HF stretching. The resonance schemes were shown in Figure 5.

Their next work at higher collision energy (0.94 kcal/mol) revealed the tunneling and shape resonance effects [61], other than Feschbash resonance leading to the forward scattering in the reaction of F + H2, which indicated that the reactive resonances played quite distinguished roles in the same reaction but at different collision energies.

Figure 6.

Theoretical and experimental 3D plots for the product of translational energy and angular distributions for the F(2P3/2) + HD(j0 = 0) → HF + D reaction at different collision energies: 0.43 kcal/mol (a); 0.48 kcal/mol (b); 0.52 kcal/mol; and (c) 0.71 kcal/mol (d).

In 2008, Ren et al. [62] measured the DCSs at several different collision energies, which suggested strong variation as a function of collision energies, resulting from the existence of strong reactive resonance states. The theoretical DCSs on the new version PES developed by Zhang and his coworkers (FXZ PES) showed good agreement with the experimental observation, as shown in Figure 6, which demonstrated that the F + H2 reaction is the first reaction which can be investigated with spectroscopy accuracy, besides the H + H2 reaction.

To have a better understanding about the reactive resonance state, the ground resonance state wave function of the F + HD → HF + D, along with the two‐dimensional (2D) minimal potential, which was optimized along the angle degree of freedom, is presented in Figure 7. It is seen there that the wave function shows features of a semibound state. The outgoing part, corresponding to the HF(v’ = 2) product, has two nodes of structure, but the inside peak has three nodes, corresponding to an excited vibrational state of v’ = 3.

Figure 7.

The 3D ground reaction resonance wave function of the F(2P3/2) + HD(j0 = 0) → HF + D reaction, along with the potential, which is optimized along angle degree of freedom.

In 2015, the resonances in the Cl +HD (v0 = 1, j0 = 0) reaction were experimentally observed, a first example besides the resonances in F + H2 and F + HD [63]. The quantum dynamics revealed that the resonances in the Cl + HD reaction come from the bond‐softening due to the strong anharmonicity over the transition barrier. Reactive resonances of this kind are expected to arise in many reactions; thus, the existence of reactive resonances is not so rare as we had thought.

6.2. Isotope effects in the O + O2 reaction

The O + O2 isotope exchange reactions play an important role in measuring the oxygen isotopic composition of a number of trace gases in the atmosphere. Their temperature dependence and kinetic isotope effects (KIEs) provide important constraints on our understanding of the origin and mechanism of these and other uncommon oxygen KIEs, which are important in the atmosphere [64,65]. Since there is a deep potential well in the process of the reaction and relative large mass of O atom, the rigorous state‐to‐state quantum dynamics calculation requires much computation. Thanks to the development of the numerical methods and computer techniques, on the most recent Dawes‐Lolur‐Li‐Jiang‐Guo (DLLJG) potential energy surface [66], calculation of the state‐to‐state differential cross sections and thermal reaction rate constant is possible [6770].

For studying the KIEs in the O + O2 reaction, using the reactant coordinate‐based time‐dependent wave packet method, where DVR method was applied for each degree of freedom, the integral cross sections (ICSs) of the 18O + 32O2 and 16O + 36O2 reactions with the initial states of (v0, j0) = (0, 1), (0, 5), (0, 9), and (0, 21) have been calculated explicitly on the newly constructed DLLJG PES. The ICSs for other j0 ≤ 21 values were estimated by a j0‐interpolation method and those for j0 > 21 were estimated by extrapolation. These ICSs yield the corresponding initial state‐specified reaction rate coefficients. Thermal reaction rate coefficients with Boltzmann, averaging over all relevant initial states were thus calculated approximately, and they exhibit clear negative temperature dependences for both reactions.

The calculated thermal reaction rate coefficients of both the 18O + 32O2 and 16O + 36O2 reactions on the DLLJG PES show a clear negative temperature dependence, in sharp contrast with the positive temperature dependence obtained on the earlier modified Siebert‐Schinke‐Bittererova (mSSB) PES [71]. In addition, there is an improved agreement between the calculated KIE and the experiment. These results validated the absence of the “reef” structure in the entrance/exit channels of the DLLJG PES, which is present in the mSSB PES.

The calculated total reaction probabilities of the 18O + 32O2 and 16O + 36O2 reactions on the DLLJG PES as a function of the total angular momentum (J) suggest that the O + O2 exchange reactions are dominated by resonances at very low collision energies (<0.2 eV) immediately above the reaction threshold, as shown in Figure 8. These resonances depend strongly on the masses of the oxygen atoms involved and/or the zero‐point energy difference between the reactant and product diatoms. Though it appears that the isotopic effects in the exchange reactions come from the ZPE difference, the underlying physical mechanism for the isotope effects is shown here to result from strong near‐threshold reactive resonances which mediate the reactions. Thus, it was first time pointed out that the accurate characterization of the reactivity for these near‐thermoneutral reactions immediately above the reaction threshold is important for correct characterization of the thermal reaction rate coefficients [47].

Figure 8.

The 2D plot of smoothened total reaction probabilities as a function of collision energy and total angular momentum J of 16O + 36O2 (6 + 88) and 18O + 32O2 (8 + 66), with the initial state of j0 = 1 on the mSSB PES. The J‐shifting rule is clearly observed in both plots, and the reactivity of 8 + 66 is larger at collision energies below 0.15 eV, resulting from the resonance enhancement.

6.3. State‐to‐state cross sections calculation of a chemical reaction

Theoretical method for a state‐to‐state calculation is of vital importance in the field of molecular reaction dynamics, which enables a detailed comparison between theoretical results and the experimental differential cross sections by crossed molecular beams [72]. Traditionally, there is only product coordinate‐based method to achieve this goal using time‐dependent quantum wave packet method in Jacobi coordinates. In 1996, John Zhang and coworkers put forward an efficient reactant‐product decoupling method for a direct reaction. Recently, reactant coordinate‐based method and transition state wave packet method were also proposed for a state‐to‐state reaction calculation.

The transition state wave packet (TSWP) method is a relatively new method for computing the reactive S‐matrix elements, and has recently been demonstrated to be accurate and efficient for computing J = 0 state‐to‐state reaction probabilities for three, four, and also six atoms [7379]. The TSWP method is based on the quantum transition state theory of Miller [80,81], who provided a direct way to calculate the cumulative reaction probability (CRP), N(E), and thermal rate coefficient, k(T), from a flux correlation function on the dividing surfaces located near the transition state [82]. Recently, Manthe and coworkers extended this formulism to compute the S‐matrix elements from generalized flux correlation functions of TSWPs [49,83]. Specifically, a set of initial TSWPs are calculated as the eigenstates of the thermal flux operator, which is defined in the transition state region as [84]


where T is a reference temperature in Kelvin, and kB denotes the Boltzmann constant. These TSWPs are then propagated independently into both the reactant and product arrangement channels, and the S‐matrix elements are obtained from the resulting cross‐correlation functions:


where the energy‐resolved projection amplitudes in the product (p) and reactant (r) channels, Aυpjpn(E) and Anυrjr(E), are given as Fourier transforms of the appropriate cross‐correlation functions:


where {fTn, |fTn} are the thermal flux eigenpairs at the reference temperature of T, and ηυpjp+(E) and ηυrjr(E) are energy‐normalizing factors of the asymptotic wave functions (|Φυpjp+ and |Φυrjr) for the two channels.

In addition to its conceptual clarity, the TSWP method has several numerical advantages. First, the propagation of each TSWP is essentially an inelastic calculation in the appropriate Jacobi coordinates, in which much smaller bases/grids are needed. One of course needs to transform the initial TSWPs prepared in one coordinate into the other, but this transformation is performed only once in a small region near the transition state. Second, the existence of either a pre‐reaction or a post‐reaction well has a limited impact on the computational costs, as the analysis plane can be placed deep into the asymptotic region without significantly impeding the numerical efficiency. Third, in contrast to the initial state‐specific wave packet (ISSWP) approach, the entire S‐matrix can be obtained at all energies in one complete calculation. Finally, the multiple propagations of TSWPs are “embarrassingly” parallel, adding further to its numerical efficiency. However, with the increase of energy, the number of involved TSWPs grows rapidly, especially for those systems with small vibrational frequencies of the activated complex.

The reactant coordinate‐based (RCB) method for computing the reactive S‐matrix elements has been widely used for the triatomic reaction systems [40,8589], but it is only recently that the RCB method is extended to tetra‐atomic systems [90,91]. The key idea of the RCB method is to propagate the time‐dependent initial wave packet deep into the product asymptotic region, where the product states are well defined on an analysis plane. There are generally two schemes to define the product states on the analysis plane: one evaluates product states on the grids of reactant coordinates with an interpolation scheme, while the other projects product states onto a set of intermediate coordinates. The interpolation scheme saves grid points in expressing the product wave functions, thus saving computer memory, and the intermediate coordinate method is numerically more efficient.

In the intermediate coordinate scheme, the projection plane for a particular product channel is designed to include the corresponding scattering coordinate that is defined as the separation of the two products. This enables the expression of the product wave packet in a product form, in which the product wave packet in the scattering coordinate is simply a delta. The remaining degrees of freedom are chosen from reactant coordinate, which renders the projection of time‐dependent wave packet on the intermediate coordinate very efficient. The calculation of overlaps between the time‐dependent wave packets and the product states is carried out on the intermediate coordinate by transforming them from the respective reactant and product Jacobi coordinates into a set of intermediate coordinates.

Compared to the TSWP method, the RCB method has the advantage of propagating the wave packet in only the reactant Jacobi coordinates, which allows the simultaneous analysis of the two product channels with only a single propagation. However, RCB has a genesis of ISSWP method, and it only calculates a column of the S‐matrix.

6.4. Dissociative chemisorption of water on transition‐metal surfaces

The interaction of molecular species with metal surfaces is of great importance to many industrial applications. The dissociative chemisorption of water on transition‐metal surfaces plays a pivotal role in understanding many industrial heterogeneous processes such as steam reforming, and represents the rate‐limiting step in low‐temperature WGS reaction on copper catalysts [92]. Tremendous progress has been made for the dissociative chemisorption dynamics.

Figure 9.

Schematic of the TSWP method. The TSWPs are prepared near the transition state region and then propagated to both the reactant and product sides to resolve the state‐to‐state information (Adopted from Ref. [4]).

Figure 10.

Two‐dimensional PES plot in the RCB method for both the H’ + H2O → H’H + OH abstraction (Abs.) and H’ + H2O → H + H’OH exchange (Exc.) channels in the reactant Jacobi coordinates R and r1. Other degrees of freedom are optimized. The initial wave packet is prepared in the reaction‐asymptotic region, and two projection planes are located in the asymptotic regions of the corresponding product channels (adopted from Ref. [67]).

A total of 9 degrees of freedom should be considered for the dynamics of H2O molecule on a rigid surface (as shown in Figure 11), rendering it formidable to carry out a fully coupled nine‐dimensional (9D) quantum dynamics calculations. Guo and coworkers performed six‐dimensional (6D) quantum dynamics calculations to study the mode specificity of H2O and bond selectivity of HOD on a rigid flat Cu(111) surface [9395], employing their 6D PESs developed by permutationally invariant polynomials (PIP). The 6D model neglects the effects of impact sites and surface corrugation, because the surface‐lateral coordinates (X and Y) and the azimuthal angle were fixed at the values of transition state. The first quantum‐state resolved molecular beam experiment on the dissociative chemisorptions of D2O on Ni(111) demonstrated a large enhancement in reactivity upon excitation of the asymmetric stretching mode of D2O [96]. The observed mode specificity was semiquantitatively understood by Guo and coworkers using the 6D quantum dynamical approach on a 6D PES. Recently, Jiang and Guo investigated 6D site‐specific dissociation probabilities of this system, but based on a new 9D PES developed by permutation‐invariant polynomial‐neural network (PIP‐NN) method [97].

Figure 11.

The 9D dissociation probability and 7D site‐specific (top, bridge, fcc, hcp, and TS sites) probabilities with H2O initially in the ground ro‐vibrational state (000). The 9D molecular coordinates of the H2O/Cu(111) system are shown in the inset.

Very recently, we reported the first seven‐dimensional (7D) quantum dynamics study for the dissociative chemisorption of H2O on Cu(111), based on an accurate 9D PES developed by neural network approach [98,99]. The dissociation probabilities exhibit strong azimuthal angle dependence, and large differences were seen between 7D and 6D results, indicating that the 6D quantum model, neglecting the azimuthal angle, can introduce substantial errors. A significant step forward in simulating gas‐surface reactions from the first‐principles was recently achieved, where we reported the first full‐dimensional quantum dynamics study of chemisorption process of H2O on rigid Cu(111), with all 9 degrees of freedom (9D) fully coupled on a global potential energy surface [100]. We found that the full‐dimensional quantum‐mechanical reactivity is quite different from the corresponding results obtained by previous reduced‐dimensional models, indicating the importance of the challenging full‐dimensional quantum‐mechanical calculations to achieve a quantitatively accurate understanding of this reaction. The excitations in vibrational modes of H2O, in particular the stretching modes, are more efficacious than increasing the translational energy in promoting the reaction, much stronger than observed in previous reduced‐dimensionality quantum studies. The full‐dimensional quantum‐mechanical calculations not only offer the dynamic features with the highest accuracy but also allow a conclusive examination of previous dynamical approximations.

6.5. Electronic states of H2+ and H2 by mapped DVR using cylindrical coordinate

For describing the low‐lying electronic state of a diatomic molecule, the prolate spheroidal coordinate is a natural choice. However, it may be not an optimal choice for describing high‐lying state. The Hamiltonian operator of electrons in a diatomic molecule is very simple and easy for numerical evaluation. It is well known that the Coulomb singularity numerically is very difficult to deal with, and in cylindrical coordinate, it is even more difficult.

To overcome this difficulty, a numerical scheme by a combination of mapped sinc DVR and cosine DVR method was introduced in the work by Lin and Sun [101]. The mapped function for variables is as follows:


Figure 12.

A typical 2D distribution of the grid points of H2+ in cylindrical coordinate after variable's mapping. In the bottom, an enlarged part of the distribution is shown.

The parameter A0 is applied to determine the grid range of ρ, but the parameter β is to adjust the relative distribution of the grid points of ρ degree of freedom. The parameter E0, which approximately corresponds to the energy of the reference state, following a phase optimization philosophy, is used to adjust the distribution of the grid points of z degree of freedom, along with Z0 parameter. Z0 approximately corresponds with the nuclear charges. The parameter z0 makes function g(0)=0, and ρ1 is the starting grid point of ρ degree of freedom. The evenly distributed grid points in the coordinate x1 and x2 are shown in Figure 12 as a function of ρ and z, where the grid points cluster around the nuclear positions.

With this mapped DVR method, the energy of the lowest electronic state, which is most sensitive to the accuracy of the description of the singularity, can be given with 10 significant numbers only with 65 and 100 grid points for ρ and z. The method also was applied to H2 molecule, including correlation of electrons. It was found that a straightforward extension of the method for H2+ only is capable of giving the energy of ground electronic state with three or four significant digits, with the parameters as follows:


The nuclear distance is 1.4 Bohr (a.u.). The grid range for ρ1 and ρ2 is [0, 20>] a.u., for z1 and z2, it is [-30, 30] a.u., and for θ, it is [0, π].

Figure 13.

The ground state and low‐lying excited states potential energy curves of H2, calculated by the mapped DVR method and the MRCI/AV6Z method by the Molpro package.

Figure 13 presents the potential energy curves of the ground electronic state and low‐lying excited states of H2 as a function of nuclear distance, calculated by the mapped DVR method (symbols) with the same parameters listed above, along with those (lines) obtained with the MRCI/AV6Z method in the Molpro package [44]. The potential energy curves of the excited electronic states obtained by the mapped DVR method are always lower than those calculated by the MRCI/AV6Z method; however, the curve of the ground electronic state by the mapped DVR method is a little higher than the one obtained by the MRCI/AV6Z method in the dissociation region, because of the strong electronic correlation energy of the ground electronic state. Further numerical studies indicate that by making the grid point distribute more densely around the middle position of the two nuclei, the accuracy of description of the mapped DVR method for the electronic correlation effects can be improved by two orders, which may be enough for a quantum calculation of the electronic dynamics induced by ultrashort laser pulses [102].

6.6. Electronic dynamics in bichromatic circularly polarized attosecond laser fields

Bichromatic circularly polarized attosecond laser pulses are new tools for investigating electron dynamics in strong field ionization of atoms, molecules, and even surfaces, following early experiments on polarization properties of high‐intensity high‐order harmonic generation [103]. Photoionization in atoms and molecules by various frequency combinations of co‐rotating and counter‐rotating bichromatic circularly polarized attosecond ultraviolet pulses exhibits signature of spiral interference patterns in photoelectron momentum distributions. The generation of electronic currents such as vortices [104106] has been proposed as sources of attosecond magnetic field pulses [107,108]. It is found that the helicity of circularly polarized pulses is the essence for the spiral photoelectron distributions. Starace and coworkers [82] have shown that for the cases with two identical frequency circularly polarized attosecond ultraviolet pulses, circularly symmetric patterns independent of the electron ejection angle are produced in photoelectron momentum distributions in atomic He. The corresponding momentum rings are also insensitive to the pulse phases. Electronic vortices can only be obtained by opposite helicity laser pulses.

We numerically solve the three‐dimensional time‐dependent Schrödinger equation of the aligned molecular ion H2+ at equilibrium. A five‐point finite difference method combined with fast Fourier transform technique is adopted to describe the molecular Hamiltonian in spatial space [85]. The electron wave function is propagated in time domain with high‐order split‐operator methods [109,110]. We use an efficient method by calculating a radial flux (electron current density at an asymptotic point) to simulate the ionization spectra [111]. At large asymptotic point, the angular flux distributions can be ignored. As a result, we only need to consider the radial part of the electronic flux along the radial direction. Results show that in bichromatic circular polarization ionization, spiral interference patterns can be observed for both cases with co‐rotating and counter‐rotating components, which are sensitive to the pulse frequencies, phase, and time delays, as shown in Figure 14. We analyze these results by an attosecond perturbation ionization model. Coherent electron wave packets with same kinetic energies created respectively by the two color pulses in the continuum interfere with each other. Photoionization distributions are functions of the photoelectron momentum and the ejection angle, thus giving rise to spiral distributions. Such sensitivity of the spiral electron distributions to the laser parameters offers an ideal means of characterizing these pulses and timing ultrafast photoionization processes by such circular attosecond pulses.

Figure 14.

Photoelectron momentum distributions of z‐aligned H2+ in bichromatic (x, y) circularly polarized attosecond ultraviolet laser pulses with (upper row) co‐rotating and (bottom row) counter‐rotating components at different frequencies.


7. Conclusion

Recent developments of quantum wave packet in atomic and molecular dynamics and their applications are reviewed. The continuous expansion of computing power and numerical techniques in quantum wave packet method leads to solve the molecular dynamics with increasingly complicated problems. Foreseeably, more general applications of the quantum wave packet method will result in more exciting discoveries in exploration of the quantum dynamics of molecules and atoms.



This work was supported by the National Basic Research Program of China (973 program, No. 2013CB922200), the National Natural Science Foundation of China (Grant Nos. 21222308, 21103187, and 21133006), the Chinese Academy of Sciences, and the Key Research Program of the Chinese Academy of Sciences.


  1. 1. Kosloff R, Time‐dependent quantum‐mechanical methods for molecular dynamics. Journal of Physical Chemistry. 1988;92:2087–2100. DOI:10.1021/j100319a003
  2. 2. Kosloff R, Propagation methods for quantum molecular dynamics. Annual Review of Physical Chemistry. 1994;45:145–178. DOI:10.1146/annurev.pc.45.100194.001045
  3. 3. Leforestier C, Bisseling RH, Cerjan C, Feit MD, Friesner R, Guldberg A, Hammerich A, Jolicard G, Karrlein W, Meyer HD, Lipkin N, Roncero O, and Kosloff R, A comparison of different propagation schemes for the time dependent Schrödinger equation. Journal of Computational Physics. 1991;94:59–80. DOI:10.1016/0021‐9991(91)90137‐A
  4. 4. Guo H, Recursive solutions to large eigenproblems in molecular spectroscopy and reaction dynamics. Reviews in Computational Chemistry. 2007;25:285–347. DOI:10.1002/9780470189078.ch7
  5. 5. Feit MD, Fleck JA, and Steiger A, Solution of the Schrödinger equation by a spectral method. Journal of Computational Physics. 1982;47:412–433. DOI:10.1016/0021‐9991(82)90091‐2
  6. 6. Kosloff R and Tal‐Ezer H, A direct relaxation method for calculating eigenfunctions and eigenvalues of the schrödinger equation on a grid. Chemical Physics Letters. 1986;127:223–230. DOI:10.1016/0009‐2614(86)80262‐7
  7. 7. Yao GH and Wyatt RE, Stationary approaches for solving the Schrödinger equation with time-dependent Hamiltonians. Journal of Chemical Physics. 1994;101:1904–1913. DOI:10.1063/1.467700
  8. 8. Peskin U, Kosloff R, and Moiseye N, The solution of the time dependent Schrödinger equation by the (t,t’) method: The use of global polynomial propagators for time dependent Hamiltonians. Journal Chemical Physics. 1994;100:8849–8855. DOI:10.1063/1.466739
  9. 9. Luo M, Liu QH, and Li ZB, Spectral element method for band structures of two‐dimensional anisotropic photonic crystals. Physical Review E. 2009;79:026705. DOI:10.1103/PhysRevE.79.026705
  10. 10. Rescigno TN and McCurdy CW, Numerical grid methods for quantum‐mechanical scattering problems. Physical Review A. 2000;62:032706. DOI:10.1103/PhysRevA.62.032706
  11. 11. Schneider BI and Collins LA, The discrete variable method for the solution of the time‐dependent Schrödinger equation. Journal of Non‐Crystalline Solids. 2005;351:1551–1558. DOI:10.1016/j.jnoncrysol.2005.03.028
  12. 12. Mazziotti DA, Spectral difference methods for solving differential equations. Chemical Physics Letters. 1999;299:473–480. DOI:10.1016/S0009‐2614(98)01324‐4
  13. 13. Boyd JP, Sum‐accelerated pseudospectral methods: The Euler‐accelerated sinc algorithm. Applied Numerical Mathematics. 1991;7:287–296. DOI:10.1016/0168‐9274(91)90065‐8
  14. 14. Gray SK and Goldfield EM, Dispersion fitted finite difference method with applications to molecular quantum mechanics. Journal of Chemical Physics. 2001;115:8331–8344. DOI:10.1063/1.1408285
  15. 15. Mazziotti DA, Spectral difference methods for solving the differential equations of chemical physics. Journal of Chemical Physics. 2002;117:2455–2468. DOI:10.1063/1.1490344
  16. 16. Wei GW, Zhang DS, Kouri DJ, and Hoffman DK, Lagrange distributed approximating functionals. Physical Review Letters. 1997;79:775–778. DOI:10.1103/PhysRevLett.79.775
  17. 17. Light JC, Hamilton IP, and Lill JV, Generalized discrete variable approximation in quantum mechanics. Journal of Chemical Physics. 1985;82:1400–1408. DOI:10.1063/1.448462
  18. 18. Light JC, Discrete‐variable representations and their utilization. Advances in Chemical Physics. 2000;114:263–310. DOI:10.1002/9780470141731.ch4
  19. 19. Wang XG and Carrington T, A contracted basis‐Lanczos calculation of vibrational levels of methane: Solving the Schrödinger equation in nine dimensions. Journal of Chemical Physics. 2003;119:101–117. DOI:10.1063/1.1574016
  20. 20. Manolopoulos DE and Wyatt RE, Quantum scattering via the log derivative version of the Kohn variational principle. Chemical Physics Letters. 1988;152:23–32. DOI:10.1016/0009‐2614(88)87322‐6
  21. 21. Baye D, The Lagrange‐mesh method. Physics Report. 2015;565:1–107. DOI:10.1016/j.physrep.2014.11.006
  22. 22. Shizgal BD and Chen HL, The quadrature discretization method (QDM) in the solution of the Schrödinger equation with nonclassical basis functions. Journal of Chemical Physics. 1996;104:4137–4150. DOI:10.1063/1.471225
  23. 23. Chen R and Guo H, Effect of spectral range on the convergence in Lanczos algorithm, a numerical study. Chemical Physical Letters. 2003;369:650–655. DOI:10.1016/S0009‐2614(02)02040‐7
  24. 24. Lin WB, KovvaIi N, and Carin L, Pseudospectral method based on prolate spheroidal wave functions for semiconductor nanodevice simulation. Computer Physics Communications. 2006;175:78–85. DOI:10.1016/j.cpc.2006.02.006
  25. 25. Yu DQ, Cong SL, and Sun ZG, An improved Lobatto discrete variable representation by a phase optimisation and variable mapping method. Chemical Physics. 2015;458:41–51. DOI:10.1016/j.chemphys.2015.07.009
  26. 26. Huang Y, Zhu W, Kouri DJ, and Hoffman DK, Distributed approximating function approach to atom‐diatom reactive scattering: Time‐dependent and time‐independent wavepacket treatments. Journal of Physical Chemistry. 1994;98:1868–1874. DOI:10.1021/j100058a025
  27. 27. Hoffman DK, Wei GW, Zhang DS, and Kouri DJ, Shannon–Gabor wavelet distributed approximating functional. Chemical Physics Letter. 1998;287:119–124. DOI:10.1016/S0009‐2614(98)00130‐4
  28. 28. Hoffman DK, Arnold M, and Kouri DJ, Properties of the optimum distributed approximating function class propagator for discretized and continuous wave packet propagations. Journal of Physical Chemistry. 1992;96:6539–6545. DOI:10.1021/j100195a007
  29. 29. Hoffman DK and Kouri DJ, Distributed approximating function theory: A general, fully quantal approach to wave propagation. Journal of Physical Chemistry. 1992;96:1179–1184. DOI:10.1021/j100182a030
  30. 30. Wei GW, Althorpe SC, Zhang DS, Kouri DJ, and Hoffman DK, Lagrange‐distributed approximating‐functional approach to wave‐packet propagation: Application to the time‐independent wave‐packet reactant‐product decoupling method. Physical Review A. 1998;57:3309–3316. DOI:10.1103/PhysRevA.57.3309
  31. 31. Smolyak SA, Quadrature and interpolation formulas for tensor products of certain classes of functions. Soviet Mathematics Doklady. 1963;4:240–243.
  32. 32. Zenger C, Sparse grids. Technical Report 18–90 SFB342, TU Mönchen, 1990.
  33. 33. Griebel M, Schneider M, and Zenger C, A combination technique for the solution of sparse grid problems, in de Groen P, Beauwens R, editors. Iterative Methods in Linear Algebra. North Holland: IMACS, Elsevier; 1992, pp. 263–281.
  34. 34. Thomas PS and Carrington T, Using nested contractions and a hierarchical tensor format to compute vibrational spectra of molecules with seven atoms. Journal of Physical Chemistry A. 2015;119:13074–13091. DOI:10.1021/acs.jpca.5b10015
  35. 35. Shen J, Wang YW, and Yu HJ, Efficient spectral‐element methods for the electronic schrödinger equation. Sparse Grids and Applications – Stuttgart. 2014;2014:265–289.
  36. 36. Carrington JRT, Two new methods for computing vibrational energy levels. Canadian Journal of Chemistry. 2015;93:589–593. DOI:10.1139/cjc‐2014‐0590
  37. 37. Manzhos S, Dawes R, and Carrington JRT, Neural network‐based approaches for building high dimensional and quantum dynamics‐friendly potential energy surfaces. International Journal of Quantum Chemistry. 2014;115:1012–1020. DOI:10.1002/qua.24795
  38. 38. Avila G and Carrington JRT, Solving the Schroedinger equation using Smolyak interpolants. Journal of Chemical Physics. 2013;139:134114. DOI:10.1063/1.4821348
  39. 39. Avila G and Carrington JRT, Using multi‐dimensional Smolyak interpolation to make a sum‐of‐products potential. Journal of Chemical Physics. 2015;143:044106. DOI:10.1063/1.4926651
  40. 40. Avila G and Carrington JRT, Solving the vibrational Schrödinger equation using bases pruned to include strongly coupled functions and compatible quadratures. Journal of Chemical Physics. 2012;137:174108. DOI:10.1063/1.4764099
  41. 41. Gradinaru V, Fourier transform on sparse grids: Code design and the time dependent Schrödinger equation. Computing. 2007;80:1–22. DOI:10.1007/s00607‐007‐0225‐3
  42. 42. Beck MH, Jäckle A, Worth GA, and Meyer HD, The multiconfiguration time‐dependent Hartree (MCTDH) method: A highly efficient algorithm for propagating wavepackets. Physics Reports. 2000;324:1–105. DOI:10.1016/S0370‐1573(99)00047‐2
  43. 43. Li GY, Rosenthal C, and Rabitz H, High dimensional model representations. Journal of Physical Chemistry A. 2001;105:7765–7777. DOI:10.1021/jp010450t
  44. 44. Sun ZG, Yang WT, and Zhang DH, Higher‐order split operator schemes for solving the Schrödinger equation in the time‐dependent wave packet method: Applications to triatomic reactive scattering calculations. Physical Chemistry Chemical Physics. 2012;14:1827–1845. DOI:10.1039/C1CP22790D
  45. 45. Liu WT, Zhang DH, and Sun ZG, Efficient fourth‐order split operator for solving the triatomic reactive schrödinger equation in the time‐dependent wavepacket approach. Journal of Physical Chemistry A. 2014;118:9801–9810. DOI:10.1021/jp5074158
  46. 46. Sun ZG and Yang WT, An exact short‐time solver for the time‐dependent Schrödinger equation. Journal of Chemical Physics. 2011;134:041101. DOI:10.1063/1.3549570
  47. 47. Park TJ and Light JC, Unitary quantum time evolution by iterative Lanczos reduction. Journal of Chemical Physics. 1986;85:5870–5878. DOI:10.1063/1.451548
  48. 48. Blanes S, Casas F, and Murua A, Symplectic splitting operator methods for the time‐dependent Schrödinger equation. Journal of Chemical Physics. 2006;124:234105. DOI:10.1063/1.2203609
  49. 49. Chen RQ and Guo H, Evolution of quantum system in order domain of Chebyshev operator. Journal of Chemical Physics. 1996;105:3569–3581. DOI:10.1063/1.472228
  50. 50. Gray SK and Balint‐Kurti GG, Quantum dynamics with real wave packets, including application to three‐dimensional (J=0) D+H2→DH+H reactive scattering. Journal of Chemical Physics. 1998;108:950–962. DOI:10.1063/1.475495
  51. 51. Kroes GJ and Neuhauser D, Performance of a time‐independent scattering wave packet technique using real operators and wave functions. Journal of Chemical Physics. 1996;105:8690–8701. DOI:10.1063/1.472650
  52. 52. Sun ZG, Guo H, and Zhang DH, Comparison of second‐order split operator and Chebyshev propagator in wave packet based state‐to‐state reactive scattering calculations. Journal of Chemical Physics. 2009;130:174102. DOI:10.1063/1.3126363
  53. 53. van Dijk W, Brown J, and Spyksma K, Efficiency and accuracy of numerical solutions to the time‐dependent Schrodinger equation. Physical Review E. 2011;84:056703. DOI:10.1103/PhysRevE.84.056703
  54. 54. Blanes S, Casas F, Oteo JA, and Ros J, The Magnus expansion and some of its applications. Physics Reports. 2009;470:151–238. DOI:10.1016/j.physrep.2008.11.001
  55. 55. Tal‐Ezer H, Kosloff R, and Schaefer I, New, highly accurate propagator for the linear and nonlinear Schrodinger equation. Journal of Scientific Computation. 2012;53:211–221. DOI:10.1007/s10915‐012‐9583‐x
  56. 56. Sun Z, Zhao B, Liu S, and Zhang DH, Reactive scattering and resonance, in Gatti F, editor. Molecular Quantum Dynamics: From Theory to Applications. Berlin: Springer; 2014, pp. 81–116. DOI:10.1007/978‐3‐642‐45290‐1_4
  57. 57. Schaefer TP, Siska PE, Parson JM, Tully FP, Wong YC, and Lee Y, Crossed molecular beam study of F + D2. Journal of Chemical Physics. 1970;53:3385–3387. DOI:10.1063/1.1674500
  58. 58. Polanyi JC and Schreiber JL, The reaction of F + H2 →HF + H. A case study in reaction dynamics. Faraday Discussions of the Chemical Society. 1977;62:267–290. DOI:10.1039/DC9776200267
  59. 59. Dong F, Lee SH, and Liu K, Reactive excitation functions for F+p‐H2/n‐H2/D2 and the vibrational branching for F + HD. Journal of Chemical Physics. 2000;113:3633–3640. DOI:10.1063/1.1287840
  60. 60. Qiu M, Ren Z, Che L, Dai D, Harich SA, Wang X, Yang X, Xu C, Xie D, Gustafsson M, Skoedje RT, Sun Z, and Zhang D‐H, Observation of feshbach resonances in the F + H2 → HF + H reaction. Science. 2006;311:1440–1443. DOI:10.1126/science.1123452
  61. 61. Wang XA, Dong WR, Qiu MH, Ren ZF, Che L, Dai DX, Wang XY, Yang XM, Sun ZG, Fu B, Lee SY, Xu X, and Zhang DH, HF(v’ = 3) forward scattering in the F → H2 reaction: Shape resonance and slow‐down mechanism. Proceedings of the National Academy of Sciences of the United States of America. 2008;105:6227–6231. DOI:10.1073/pnas.0710840105
  62. 62. Ren ZF, Che L, Qiu MH, Wang XA, Dong WR, Dai DX, Wang XY, Yang XM, Sun ZG, Fu B, Lee SY, Xu X, and Zhang DH, Probing the resonance potential in the f atom reaction with hydrogen deuteride with spectroscopic accuracy. Proceedings of the National Academy of Sciences of the United States of America. 2008;105:12662–12666. DOI:10.1073/pnas.0709974105
  63. 63. Yang TG, Chen J, Huang L, Wang T, Xiao CL, Sun ZG, Dai DX, Yang XM, and Zhang DH, Extremely short‐lived reaction resonances in Cl + HD (v = 1) → DCl + H due to chemical bond softening. Science. 2015;347:60–63. DOI:10.1126/science.1260527
  64. 64. Fleurat‐Lessard P, Grebenshchikov SY, Schinke R, Janssen C, and Krankowsky D, Isotope dependence of the O+O2 exchange reaction: Experiment and theory. Journal of Chemical Physics. 2003;119:4700–4712. DOI:10.1063/1.1595091
  65. 65. Gao YQ and Marcus RA, Strange and unconventional isotope effects in ozone formation. Science. 2001;293:259–263. DOI:10.1126/science.1058528
  66. 66. Dawes R, Lolur P, Li A, Jiang B, and Guo H, Communication: An accurate global potential energy surface for the ground electronic state of ozone. Journal of Chemical Physics. 2013;139:201103. DOI:10.1063/1.4837175
  67. 67. Li Y, Sun Z, Jiang B, Xie D, Dawes R, and Guo H, Communication: Rigorous quantum dynamics of O + O2 exchange reactions on an ab initio potential energy surface substantiate the negative temperature dependence of rate coefficients. Journal of Chemical Physics. 2014;141:081102. DOI:10.1063/1.4894069
  68. 68. Sun Z, Liu L, Lin SY, Schinke R, Guo H, and Zhang DH, State‐to‐state quantum dynamics of O + O2 isotope exchange reactions reveals nonstatistical behavior at atmospheric conditions. Proceedings of the National Academy of Sciences of the United States of America. 2010;107:555–558. DOI:10.1073/pnas.0911356107
  69. 69. Xie W, Liu L, Sun Z, Guo H, and Dawes R, State‐to‐state reaction dynamics of 18O+32O2 studied by a time‐dependent quantum wavepacket method. Journal of Chemical Physics. 2015;142:064308. DOI:10.1063/1.4907229
  70. 70. Sun Z, Yu DQ, Xie WB, Hou JY, Dawes R, and Guo H, Kinetic isotope effect of the 16O + 36O2 and 18O + 32O2 isotope exchange reactions: Dominant role of reactive resonances revealed by an accurate time‐dependent quantum wave packet study. Journal of Chemical Physics. 2015;142:174312. DOI:10.1063/1.4919861
  71. 71. Siebert R, Schinke R, and Bittererova M, Spectroscopy of ozone at the dissociation threshold: Quantum calculations of bound and resonance states on a new global potential energy surface. Physical Chemistry Chemical Physics. 2001;3:1795–1798. DOI:10.1039/B102830H
  72. 72. Xiao CL, Xu X, Liu S, Wang T, Dong WR, Yang TG, Sun Z, Dai DX, Xu X, Zhang DH, and Yang XM, Experimental and theoretical differential cross sections for a four‐atom reaction: HD+OH→H2O+D. Science. 2011;333:440–442. DOI:10.1126/science.1205770
  73. 73. Welsch R, Huarte‐Larrañaga F, and Manthe U, State‐to‐state reaction probabilities within the quantum transition state framework. Journal of Chemical Physics. 2012;136:064117. DOI:10.1063/1.3684631
  74. 74. Zhao B, Sun Z, and Guo H, Calculation of state‐to‐state differential and integral cross sections for atom‐diatom reactions with transition‐state wave packets. Journal of Chemical Physics. 2014;140:234110. DOI:10.1063/1.4883615
  75. 75. Zhao B, Sun Z, and Guo H, Calculation of the state‐to‐state S‐matrix for tetra‐atomic reactions with transition‐state wave packets: H2/D2 + OH → H/D + H2O/HOD. Journal of Chemical Physics. 2014;141:154112. DOI:10.1063/1.4898100
  76. 76. Zhao B and Guo H, Modulations of transition‐state control of state‐to‐state dynamics in the F + H2O → HF + OH reaction. Journal of Physical Chemistry Letters. 2015;6:676–680. DOI:10.1021/acs.jpclett.5b00071
  77. 77. Zhao B, Sun Z, and Guo H, Communication: State‐to‐state dynamics of the Cl + H2O → HCl + OH reaction: Energy flow into reaction coordinate and transition‐state control of product energy disposal. Journal of Chemical Physics. 2015;142:241101. DOI:10.1063/1.4922650
  78. 78. Welsch R and Manthe U, Loss of memory in H + CH4 → H2 + CH3 state‐to‐state reactive scattering. Journal of Physical Chemistry Letters. 2015;6:338–342. DOI:10.1021/jz502525p
  79. 79. Zhao B, Sun Z, and Guo H, State‐to‐state mode specificity: Energy sequestration and flow gated by transition state. Journal of the American Chemical Society. 2015;137:15964–15970. DOI:10.1021/jacs.5b11404
  80. 80. Miller WH, Quantum‐mechanical transition‐state theory and a new semiclassical model for reaction‐rate constants. Journal of Chemical Physics. 1974;61:1823–1834. DOI:10.1063/1.1682181
  81. 81. Miller WH, Schwartz SD, and Tromp JW, Quantum‐mechanical rate constants for bimolecular reactions. Journal of Chemical Physics. 1983;79:4889–4899. DOI:10.1063/1.445581
  82. 82. Miller WH, “Direct” and “correct” calculation of canonical and microcanonical rate constants for chemical reactions. Journal of Physical Chemistry A. 1998;102:793–806. DOI:10.1021/jp973208o
  83. 83. Manthe U and Welsch R, Correlation functions for fully or partially state‐resolved reactive scattering calculations. Journal of Chemical Physics. 2014;140:244113. DOI:10.1063/1.4884716
  84. 84. Park TJ and Light JC, Quantum flux operators and thermal rate constant: Collinear H+H2. Journal of Chemical Physics. 1988;88:4897–4912. DOI:10.1063/1.454702
  85. 85. Wang T, Chen J, Yang T, Xiao C, Sun Z, Huang L, Dai D, Yang X, and Zhang DH, Dynamical resonances accessible only by reagent vibrational excitation in the F + HD → HF + D reaction. Science. 2013;342:1499–1502. DOI:10.1126/science.1246546
  86. 86. Gómez‐Carrasco S and Roncero O, Coordinate transformation methods to calculate state‐to‐state reaction probabilities with wave packet treatments. Journal of Chemical Physics. 2006;125:054102. DOI:10.1063/1.2218337
  87. 87. Sun Z, Guo H, and Zhang DH, Extraction of state‐to‐state reactive scattering attributes from wave packet in reactant Jacobi coordinates. Journal of Chemical Physics. 2010;132:084112. DOI:10.1063/1.3328109
  88. 88. Dai JQ and Zhang JZH, Time‐dependent wave packet approach to state‐to‐state reactive scattering and application to H+O2 reaction. Journal of Physical Chemistry. 1996;100:6898. DOI:10.1021/jp9536662
  89. 89. Sun Z, Lin X, Lee SY, and Zhang DH, A reactant‐coordinate‐based time‐dependent wave packet method for triatomic state‐to‐state reaction dynamics: Application to the H + O2 reaction. Journal of Physical Chemistry A. 2009;113:4145–4154. DOI:10.1021/jp810512j
  90. 90. Zhao B, Sun Z, and Guo H, A reactant‐coordinate‐based wave packet method for full‐dimensional state‐to‐state quantum dynamics of tetra‐atomic reactions: Application to both the abstraction and exchange channels in the H + H2O reaction. Journal of Chemical Physics. 2016;144:064104. DOI:10.1063/1.4941671
  91. 91. Zhao B, Sun Z, and Guo H, Reactant‐coordinate‐based wave packet method for full‐dimensional state‐to‐state quantum dynamics of tetra‐atomic reactions: Application to both the abstraction and exchange channels in the H + H2O reaction, Journal Chemical Physics. 2016;144:064104. DOI:10.1063/1.4941671
  92. 92. Chorkendorff I and Niemantsverdriet JW, Concepts of Modern Catalysis and Kinetics. Weinheim: Wiley‐VCH; 2003.
  93. 93. Jiang B, Xie D, and Guo H, Vibrationally mediated bond selective dissociative chemisorption of HOD on Cu(111). Chemical Science. 2013;4:503–508. DOI:10.1039/C2SC21393A
  94. 94. Jiang B, Li J, Xie D, and Guo H, Effects of reactant internal excitation and orientation on dissociative chemisorption of H2O on Cu(111): Quasi‐seven‐dimensional quantum dynamics on a refined potential energy surface. Journal of Chemical Physics. 2013;138:044704. DOI:10.1063/1.4776770
  95. 95. Jiang B, Ren X, Xie D, and Guo H, Enhancing dissociative chemisorption of H2O on Cu(111) via vibrational excitation. Proceedings of the National Academy of Sciences of the United States of America. 2012;109:10224–10227. doi:10.1073/pnas.1203895109
  96. 96. Hundt PM, Jiang B, van Reijzen, Maarten E, Guo H, and Beck RD, Vibrationally promoted dissociation of water on Ni(111). Science. 2014;344:504–507. DOI:10.1126/science.1251277
  97. 97. Jiang B and Guo H, Quantum and classical dynamics of water dissociation on Ni(111): A test of the site‐averaging model in dissociative chemisorption of polyatomic molecules. Journal of Chemical Physics. 2015;143:164705. DOI:10.1063/1.4934357
  98. 98. Liu T, Zhang ZJ, Fu B, Yang X, and Zhang DH, A seven‐dimensional quantum dynamicsstudy of the dissociative chemisorption of H2O on Cu(111): Effects of azimuthal angles and azimuthal angle‐averaging. Chemical Science. 2016;7:1840–1845. DOI:10.1039/C5SC03689E
  99. 99. Liu T, Zhang ZJ, Fu B, Yang X, and Zhang DH, Mode specificity for the dissociative chemisorption of H2O on Cu(111): A quantum dynamics study on an accurately fitted potential energy surface. Physical Chemistry Chemical Physics. 2016;18:8537–8544. DOI:10.1039/C6CP00034G
  100. 100. Zhang ZJ, Liu T, Fu B, Yang X, and Zhang DH, First‐principles quantum dynamical theory for the dissociative chemisorption of H2O on rigid Cu(111). Nature Communation. Under revision, 2016.
  101. 101. Lin XS and Sun ZG, Accurate mapped trigonometric discrete variable representations for Coulomb singularities in molecules: Applications with H2+ and H2 in cylindrical coordinates. Chemical Physics Letters. 2015;621:35–40. DOI:10.1016/j.cplett.2014.12.043
  102. 102. Lin XS and Sun ZG, To be published.
  103. 103. Weihe FA, Dutta SK, Korn G, Du D, Bucksbaum PH, and Shkolnikov PL, Polarization of high‐intensity high‐harmonic generation. Physical Review A. 1995;51:R3433. DOI:10.1103/PhysRevA.51.R3433
  104. 104. Barth I, Manz J, Shigeta Y, and Yagi, K, Unidirectional electronic ring current driven by a few cycle circularly polarized laser pulse: Quantum model simulations for Mg‐porphyrin. Journal of the American Chemical Society. 2006;128:7043. DOI:10.1021/ja057197l
  105. 105. Djiokap JMN, Hu SX, Madsen LB, Manakov NL, Meremianin AV, and Starace AF, Electron vortices in photoionization by circularly polarized attosecond pulses. Physical Review Letters. 2015;115:113004. DOI:10.1103/PhysRevLett.115.113004
  106. 106. Yuan KJ, Lu H, and Bandrauk AD, High‐order‐harmonic generation in molecular sequential double ionization by intense circularly polarized laser pulses. Physical Review A. 2015;92:023415. DOI:10.1103/PhysRevA.92.023415
  107. 107. Yuan KJ and Bandrauk AD, Attosecond‐magnetic‐field‐pulse generation by electronic currents in bichromatic circularly polarized UV laser fields. Physical Review A. 2015;92:023415. DOI:10.1103/PhysRevA.92.063401
  108. 108. Bandrauk AD and Lu HZ, High‐Dimensional Partial Differential Equations in Science and Engineering, CRM Lecture Series, Vol. 41. Philadelphia: American Math. Soc.; 2007, pp. 1–15.
  109. 109. Bandrauk AD and Shen H, Exponential split operator methods for solving coupled time-dependent Schrödinger equations. Journal of Chemical Physics. 993;99:1185. DOI:10.1063/1.465362
  110. 110. Bandrauk AD and Lu HZ, Exponential propagators (integrators) for the time‐dependent Schrödinger equation. Journal of Theoretical and Computational Chemistry. 2013;12:1340001. DOI:10.1142/S0219633613400014
  111. 111. Zhang DH and Zhang JZH, Quantum reactive scattering with a deep well: Time-dependent calculation for H+O2 reaction and bound state characterization for HO2. Journal of Chemical Physics. 1994;101:3671. DOI:10.1063/1.467551

Written By

Zhigang Sun

Submitted: 23 October 2015 Reviewed: 28 April 2016 Published: 24 August 2016