Numerical approaches to solving the time-dependent Schrödinger equation with different potentials

TVE 16 040 juni

Examensarbete 15 hp

Juni 2016

Numerical approaches to solving

the time-dependent Schrödinger

equation with different potentials

Lina Viklund

Teknisk- naturvetenskaplig fakultet UTH-enheten Besöksadress: Ångströmlaboratoriet Lägerhyddsvägen 1 Hus 4, Plan 0 Postadress: Box 536 751 21 Uppsala Telefon: 018 – 471 30 03 Telefax: 018 – 471 30 00 Hemsida: http://www.teknat.uu.se/student

Abstract

Numerical approaches to solving the time-dependent

Schrödinger equation with different potentials

Lina Viklund, Louise Augustsson, Jonas Melander

This project is an immersive study in numerical methods, focusing on quantum molecular dynamics and methods for solving the time-dependent Schrödinger equation. First the Schrödinger equation was solved with finite differences and a basic propagator in time, and it was then concluded that this method is far too slow and compuationally heavy for its use to be justified for this type of problem. Instead pseudo-spectral methods with split-operators were implemented, and this proved to be a far more favourable method for solving, both in regards to time and memory requirements.

Further, the pseudo-spectral methods with split-operators were used to solve the dynamics resulting from the excitation of sodium iodide by an ultra-fast laser pulse. This was modeled as two Schrödinger equations coupled with a potential modeling the laser pulse. The resulting solution made the quantum nature of the system clear, but also the limitations and advantages of different numerical methods.

ISSN: 1401-5757, TVE 16 040 juni Examinator: Martin Sjödin

Ämnesgranskare: Rikard Emanuelsson

Popul¨

arvetenskaplig sammanfattning

Allting vi ser i v˚art dagliga liv, inklusive oss sj¨alva, ¨ar uppbyggt av atomer och molekyler, och molekyl¨ar kvantdynamik beskriver hur atomer och molekyler r¨or sig och interagerar med varandra. Genom kun-skap och f¨orst˚aelse i omr˚adet ¨oppnas nya m¨ojligheter att simulera och f¨orutsp˚a resultatet av kemiska reaktioner, vilket ¨ar mycket f¨ordelaktigt inom vissa forskningsomr˚aden.

Det f¨orsta steget till att simulera molekylernas r¨orelse ¨ar att ta fram matematiska modeller som beskriver denna r¨orelse som st¨ammer v¨al ¨overens med verkligheten. Den inom kvantfysiken my-cket ber¨omda Schr¨odingerekvationen beskriver partiklar som v˚agfunktioner med viss sannolikhet att befinna sig i en viss punkt vid ett visst tillf¨alle. Partiklarna r¨or sig p˚a olika energiniv˚aer (potentialer) som kan vara konstanta eller ¨andras ¨over tiden, och beroende p˚a potential och antal partiklar blir ekvationerna olika sv˚ara att l¨osa. F¨or vissa specialfall finns ’exakta’ l¨osningar, dvs. fall som g˚ar att l¨osa ”med papper och penna”, men f¨or det mesta ¨ar modellerna komplexa och saknar exakta l¨osningar ¨

aven f¨or mycket sm˚a system.

Ett s¨att att g˚a tillv¨aga f¨or att l¨osa dessa komplexa system ¨ar att diskretisera dem och l¨osa dem nu-meriskt. Detta inneb¨ar att tids- och rumsintervallet som systemet l¨oses p˚a delas upp i ett begr¨ansat antal punkter, och problemet l¨oses sedan i varje s˚adan punkt med hj¨alp av datorer. Bland det stora antalet numeriska metoder som kan anv¨andas finns givetvis mer eller mindre noggranna, mer eller mindre tidskr¨avande samt mer eller mindre minneskr¨avande metoder. Det g¨aller allts˚a att hitta en metod som ¨ar tillr¨ackligt noggrann, inte tar f¨or l˚ang tid och inte kr¨aver f¨or mycket minne. Det h¨ar arbetet syftar dels till att l¨osa Schr¨odingerekvationen f¨or ett enpartikelsystem b˚ade i tidsoberoende och tidsberoende potentialer, dels till att l¨osa Schr¨odingerekvationen f¨or ett flerpartikelsystem med en tidsberoende potential. Metoderna som anv¨ands analyseras och j¨amf¨ors med varandra med avseende p˚a noggrannhet, tid och minne.

Acknowledgements

Contents

2.1.4 Potential for the NaI-system . . . 5

2.1.5 Coupling potential . . . 6

2.2 Model problem . . . 6

2.3 Discretization . . . 6

2.3.1 Spatial methods . . . 7

2.3.2 Methods for time stepping . . . 8

2.4 TDSE for NaI-system . . . 10

1

Introduction

Quantum molecular dynamics (QMD) describes the time evolution of atomic and molecular interac-tions on the quantum scale. Simulating such dynamics gives an insight into the behaviour of molecules and atoms, and this knowledge in turn helps to e.g. understand how to better control the outcome of chemical reactions. Experimental studies of chemical reactions by A. H. Zewail on a type of system similar to the one studied in this paper [1] was rewarded by the Nobel Prize in chemistry in 1999. This project will focus on solving the dynamics resulting from the excitation of sodium iodide (NaI) by an ultra-fast laser pulse, operating on the femto-second scale.

To understand QMD, one needs accurate mathematical models to describe the system of particles. In order to simulate these interactions, efficient and accurate numerical methods are needed, since most systems of particles are too complex to be described by equations with analytic solutions. Numerical models are obtained by using discrete representations of the continuous mathematical models, and this can be done in many ways, some more efficient and accurate than others.

A model problem (presented in section 2.2) with Morse potential (presented in section 2.1.2) will be solved using both basic finite differences which are widely used in numerical calculations, and pseudo-spectral methods combined with the split operator technique, which have successfully been applied to problems of this type. The complete NaI system will be solved with pseudo-spectral meth-ods combined with the split operator technique. This will be followed by a discussion comparing the two methods.

2

Theory

The time-dependent Schr¨odinger equation (TDSE) is an n-dimensional complex-valued partial differ-ential equation (PDE) describing the motion of particles in space and time. In general the TDSE can be formulated as i¯h∂ ∂tΨ(~r, t) = − ¯ h2 2m∇ 2_{Ψ(~r, t) + V (~r, t)Ψ(~r, t), (1)}

where ~r is an n-dimensional vector, t is time, i =√−1, m the mass of the object described, and ¯h is Planck’s constant divided by 2π. Equation (1) may be rewritten on the form

i¯h∂

∂tΨ(~r, t) = ˆH(t)Ψ(~r, t), (2) where ˆH(t) is the time dependent Hamiltonian operator given by

ˆ

H(t) = −¯h

2

2m∇

2_{+ V (~r, t) = ˆT + V (~r, t), (3)}

where ˆT is the kinetic energy operator, and V (~r, t) a time-dependent potential.

A diatomic system with two coupled electronic potentials can be modeled as two PDEs.          ˆ H1Ψ1(~r, t) + Vc(~r, t)Ψ2(~r, t) = i¯h ∂ ∂tΨ1(~r, t) ˆ H2Ψ2(~r, t) + Vc†(~r, t)Ψ1(~r, t) = i¯h ∂ ∂tΨ2(~r, t). (4) ˆ

H1 and ˆH2 are Hamiltonian operators containing V1(~r) and V2(~r) respectively. V1(~r) and V2(~r) are

potential energy surfaces corresponding to different states. Vc(t) is a coupling potential and Vc†(t) is

Further, the inner product of the quantum wave function with its complex conjugate, Ψ Ψ = Z ∞ −∞ Ψ∗Ψdr,

can be interpreted as the total probability density, which should always be preserved. A good numer-ical method conserves the probability density, and this is thus a necessary but not sufficient condition for correctness.

When performing calculations on this scale it is convenient to use atomic units. In this unit sys-tem me, the electronic mass, e, the electron charge, and ¯h are all set to unity, and all other units are

scaled accordingly.

2.1 Potentials

The choice of potential V (r) decides what system that will be modeled. A potential with an easy analytic solution can be used for testing the accuracy, convergence and stability of the numerical methods. A more advanced potential that cannot be solved analytically may be a better approximation of the real problem. Below a number of potentials that will be used in the project are presented.

2.1.1 Harmonic oscillator

The accuracy of the numerical solutions to the Schr¨odinger equation may be difficult to ascertain, not knowing the true solution before hand. The harmonic oscillator

V (r) = 1 2mω

2_r2_{, (5)}

where m is the mass and ω is the eigenfrequency of the particle, has a known analytic solution. This can be used to test aspects of the numerical methods that are used. This is possible because for a time-independent case we have that

ˆ

Hψ(r) = Eψ(r).

This means that when the Hamiltonian operates on a stationary wave function, the result will be a constant E multiplied with the same wave function, where E is the eigenenergy of this stationary state. For a harmonic oscillator potential the formula for these eigenenergies are known,

En= n +

1 2
¯hω,

where n is the quantized energy level. Since these can be calculated analytically, they can be compared to numerical calculations of the eigenenergies, and thus be used to investigate the numerical methods. Figure 1 gives a visual representation of the potential.

The initial values for harmonic oscillations can be approximated using a Gaussian wave-packet, ψ(r), which can be expressed as

ψ(r) = µeα2(r−r0)2_{, (6)}

where µ, the normalization constant, is given by

µ = q_Rr 1

1

r0 ψ

∗_(r0_)ψ(r0_)dr0

and α governs the width of the function.

2.1.2 Morse potential

The Morse potential, visualized in figure 2, can be a better approximation of a potential energy surface when working with diatomic molecules. It is given by

V (r) = De(1 − ea(r−r0))2, (7)

where r − r0 is the distance from a potential well at r0, a and De are parameters that relate to the

”width” and ”depth” of the potential. These parameters are molecule specific.

Figure 2: The Morse potential.

2.1.3 Time-dependent potential

When the Hamiltonian is time-dependent, an extra time-dependent term is added to V (r). Here the time dependent part model the systems interaction with a photon. In this case with a laserpulse

V (r, t) = V (r) + Acos(ωpt),

where A is a constant, and ωp is the frequency of the time-dependent perturbation.

2.1.4 Potential for the NaI-system

For the NaI-system described in section 2.4 the ground state potential is given by

and the excited state is given by

Ve(r) = Aee −r

β . (9)

These potentials with the parameter values tabulated in table 11 of Appendix A have previously been shown to give accurate results [8] for this kind of simulation. The potentials can be seen in figure 3.

2.1.5 Coupling potential

The laser pulse’s interactions with the NaI-system can be modeled as a coupling potential. The spatial part of the coupling potential is given by

Vc,r(r) = A12e−β12(r−r0)

2

, (10)

and the temporal part is given by

Vc,t(t) = ηe−αt(t−t0)

2

cos(ω ∗ (t − t0)), (11)

where η is the strength of the laser pulse and

αt=

2 log(2)

ξ2 , (12)

where ξ is the width of the pulse.

2.2 Model problem

The model problem that will be used to investigate the numerical methods is a TDSE with one dimension in space, both time-independent and time-dependent. This is given by

i¯h∂

∂tΨ(r, t) = ˆHΨ(r, t), (13) where r − r0 is the distance to the center of some potential well at r0, and ˆH the Hamiltonian. The

time-independent Hamiltonian is given by

ˆ

H = − ¯h 2m

∂2

∂r2 + V (r). (14)

The time-independent case will be used in order to check properties of the numerical methods, such as order of accuracy and convergence, as explained in section 2.1.1. In other regards than this, the time-independent case is not that interesting, since a model of a chemical reactions needs to involve a time-dependent potential. The time-dependent Hamiltonian is given by

ˆ

H = − ¯h 2m

∂2

∂r2 + V (r, t), (15)

and for this case there is no simple analytic solution. An analytic approximation can be made [4], but for the time-dependent cases it is more efficient to use numerical methods to calculate the solutions. Reducing equation (1) to equation (13) offers an easy way to test numerical methods without demand-ing powerful computers to calculate the solutions.

2.3 Discretization

Applying numerical methods to a PDE problem requires it to be dicretized. This can be done by introducing a uniform grid in both time and space

rj = j∆r j = 0, · · · , N ∆r = 1/N

where rj are the grid points in space, N the number of grid points in space and ∆r the step size in

space. tn are the grid points in time, and ∆t is the step size in time. For each grid point j, n we then

get a solution ψn_j that approximates the true solution, ψ_jn≈ ψ(rj, tn).

Discretized initial conditions (I.C.) are required to get accurate solutions when using other than Harmonic potentials. These are obtained by solving the matrix equation

HΨ0(r) = EΨ0(r) (16)

for the eigenvector Ψ0(r). E is the eigenenergy for the initial state.

2.3.1 Spatial methods

The model problem will be solved using both a simple method, the finite difference method with a centered difference, and a more advanced method, the pseudo-spectral method.

The main idea with finite difference methods is to approximate derivatives and solutions to differ-ential equations [9] with linear combinations of function values in neighboring grid points. One of the most basic ways to approximate a second order derivative with finite differences is the centered differ-ence, which uses grid points centered around a middle point in order to approximate the derivative. How many grid points that are used for the centered difference, and how they are weighted depend on the desired order of accuracy. The TDSE approximated with a second order centered difference in space reads i¯h∂ ∂tψj = − ¯ h2 2m ψ_j+1− 2ψ_j+ ψ_j−1 (∆r)2  + V (rj)ψj = Hψj, (17)

which is a system of ordinary differential equations, where

H = −i¯h∆t 2m(∆r)2        2 −1 −1 2 −1 . .. ... ... 2 −1 −1 2        | {z } T +i∆t ¯ h          V (r) . .. . .. . .. V (r)          | {z } V(r) . (18)

T is the differentiation matrix, here for a centered finite difference with order of accuracy 2, and V(r) is the potential matrix. With H as in equation 18, we have Dirichlet boundary conditions, set to zero. Depending on the order of accuracy of the centered difference T will have more diagonals filled in, corresponding to the coefficients in table 1.

Table 1: Coefficients of the central difference for second order derivative

Accuracy -4 -3 -2 -1 0 1 2 3 4

2 1 -2 1

4 -1/12 4/3 -5/2 4/3 -1/12

6 1/90 -3/20 3/2 -49/18 3/2 -3/20 1/90

8 -1/560 8/315 -1/5 8/5 -205/72 8/5 -1/5 8/315 -1/560

order methods can then be compared with the reference solution in order to determine their order of accuracy. The order of accuracy tells how fast the solutions converge. When an n-tupled grid is used the relative error decreases by a factor of np where p is the order of accuracy of the method.

If a system of PDE:s with dirichlet boundary conditions (B.C.) can be approximated with a system of identical PDE:s with periodic B.C., a pseudo-spectral method (PSM) may be applied. The PSM’s spatial accuracy is formally infinite [5], the main error in the PSM is the inherited error from the model approximation. Because of the accuracy of the PSM, equally accurate initial conditions are needed, these are obtained by solving equation (16), where

H = Tt+ V, (19)

Ttis a circulant matrix that replaces the standard finite difference differentiation matrix for the kinetic

energy and V is a potential. Since all potentials in our case are analytic, the use of Tt in H is enough

to get equally accurate I.C.:s. The PSM take advantage of the kinetic energy operator being local in the momentum domain, or k-space, making it a diagonal matrix when discretized, with p

2

2m on the diagonal. Table 2 shows the position, momentum and kinetic energy in the position and momentum domain. Equation (13) can be reformulated with the help of the discrete Fourier transform (DFT)

i¯h∂ ∂tΨ(rj, t) = ˆT ψ(rj, t) + V (rj)Ψ(rj, t) ⇒ i¯h ∂ ∂tΨ(rj, t) = F −1 DkF Ψ(rj, t) + V (rj)Ψ(rj, t), (20) where Dk = diag( p2

2m) is the DFT of the kinetic energy operator. These operations are preferably done with the fast Fourier transform fft in MATLAB.

Table 2: Position, momentum and kinetic energy in the position and momentum domain (k − space).

Operator Position domain Momentum domain ˆ x x i¯h_∂p∂ ˆ p2 2m −i¯h ∂2 ∂x2 ˆ p2 2m

2.3.2 Methods for time stepping

For a time-independent Hamiltonian the analytic solution to the TDSE is given by

ψ(t) = e−¯hiHtψ₀, (21)

where ψ0 is the initial state and ˆH is the Hamiltonian matrix. A simple way of time stepping in this

case is then to use MATLAB’s expm() function, which is used as a propagator. The function of a propagator is the propagation (in time) of the function it is applied to. In equation (21) the expo-nential term is the propagator. In this case the solution is analytic, and thus the order of accuracy is infinity.

However, if the Hamiltonian is time-dependent, which is the case with a time-dependent potential, it will include a matrix for the kinetic energy and a matrix for the time-dependent perturbation, which do not necessarily commute. In this case the analytic solution is not known. Instead the Magnus expansion can be used to approximate the time-dependence effect on the solution[6]. If the expansion is truncated after the first term it can be written as

First-order Magnus expansion is second order accurate in time [2]. Using the first term only will not be a problem since the Magnus expansion is unitary, meaning that the probability is conserved even if it is truncated. As in the time-independent case, expm() will be used for propagating, but since the integral in the exponential is approximated by ( ˆH(t) − ˆH(t0))∆t which has order of accuracy 2, as

may be shown by Taylor-expansion. Thus the total order of accuracy in time will be 2.

The calculations are performed in Matlab. For a time-independent Hamiltonian the kinetic energy matrix, the Hamiltonian matrix and the matrix exponential are computed once. They are then used in a for-loop to do the time stepping. When the Hamiltonian is time-dependent it changes for every time step. Therefore both the Hamiltonian and the matrix exponential must be computed anew for every time step. Matrix-vector multiplication takes O(N2) operations and the matrix exponential operator takes approximately O(N3) operations. This gives that each time step in the solution has computational complexity O(N2) for time-independent potential, and O(N3+N2) for time-dependent. Concerning the memory requirements it is the same at the end stage with both time independent and time dependent potential. However the time-dependent solution will pre-allocate more memory be-cause it needs to compute the matrix exponential in every loop and a full matrix, O(N2), is needed to store the matrix exponential. Computations of the matrix exponential are computational expensive and takes a lot of memory, O(N2).

A more advanced method that addresses this concern is the split-operator method [7]. This method is based on splitting the Hamiltonian into a kinetic energy part and a potential energy part. This means that the solution to the TDSE can be written as

ψ(r, t) = ei(T +V )∆tψ0. (24)

Since T and V are matrices that do not necessarily commute, ei(T +V )∆t6= eiT ∆t_eiV ∆t _{in this case, but}

equation 24 can be approximated as

ψ(r, t + ∆t) ≈ eiT ∆teiV ∆tψ(r, t) + O((∆t)2). (25) In order to achieve a better approximation [2] the kinetic energy matrix can be split into

ψ(r, t + ∆t) ≈ eiT ∆t/2eiV ∆teiT ∆t/2ψ(r, t) + O((∆t)3). (26) Pairing the split-operator method with the pseudo-spectral method described in section 2.3.1 makes for an n-efficient method, since the splitting makes it easy to diagonalize the kinetic energy matrix. It takes O(N log N ) operations and O(N ) in allocated memory.

2.4 TDSE for NaI-system

The complete NaI system can be modeled as

i¯h∂ ∂t ψ1(r, t) ψ2(r, t)  =T + Vg(r) V † c(r, t) Vc(r, t) T + Ve(r)  ψ1(r, t) ψ2(r, t)  , (27)

where Vg is the ground state potential matrix, Ve is the first excited state potential matrix, Vc is the

coupling potential matrix, and Vc† is the complex conjugate of Vc. Ψ1(r, t) and Ψ2(r, t) are vectors

corresponding to Vg and Ve respectively.

Applying PSM and split operators to equation (27) yields Ψ(r, tn+1) = F−1e− i 2¯hΞ∆tF e− i ¯ hW ∆tF−1e− i 2¯hΞ∆tF Ψ(r, t), (28)

Which when applied from left to right will result in stepping the TDSE one time step forward. Here we have that Ξ = diag(Dk), with Dk= diag(

p2

2m). W is a matrix of potential matrices defined as W = Vg(r) V † c(r, t) Vc(r, t) Ve(r)  (29)

which is not diagonal. Because Vc(r, t) is not a complex-valued function in the case treated here, we

get Vc(r, t) = Vc†(r, t). The pre calculation of e−

i ¯

hW ∆t is desired for increased computational speed

and reduced memory requirements. The calculation of the matrix exponential can be done by

e−¯hiW ∆t = exp − i ¯ h  Vg(r) Vc(r, t) Vc(r, t) Ve(r)  ∆t ! = exp − i ¯ h Vg(r) + Ve(r)
∆t 2 ! cospQ∆t 2¯h 1 0 0 1  + (30) isin √ Q∆t_2¯h
√ Q Ve(r) − Vg(r) −2Vc(r, t) −2Vc(r, t) Vg(r) − Ve(r) ! = exp − i ¯ h Vg(r) + Ve(r)
∆t 2 ! γ

which is shown in [8] and [2], and Q = (Ve(r) − Vg(r))2+ 4Vc2(r, t). This results in

Ψ(r, tn+1) = F−1e− i 2¯hΞ∆tF e− i ¯ h(Vg(r)+Ve(r))∆tγF−1e− i 2¯hΞ∆tF Ψ(r, t) (31)

which is the final form representing the time stepping of the TDSE when applied from right to left.

Solving the TDSE on coupled potential surfaces, as is the case for the NaI-system, give rise to the phenomenon where the probability density function splits over the two surfaces shown in figure 3. The probability density on a surface is called a population and, as always when dealing with probabilty densities, the sum of the populations are unity. When the laser pulse excites the molecule a population will appear on the excited state and then move along this surface until it reaches the crossing of the excited state and the ground state. Here the population will split, so there is a probability both of the probability density continuing along the excited state and of it moving back down on the ground state. If the density continues on its excited state or drift further apart along the ground state, the molecule will separate in to its individual atoms. Depending on which potential surface the density traveled along it will either be an ion separation or a separation of neutral atoms. If the molecule separates into two atoms the models described will no longer be applicable on the system.

The PSM with split operators assumes a spatial periodicity in the system i.e. that the boundary conditions can be set as periodic. This is not actually true for this system, an approximation of periodic boundary conditions is made at such a distance between the atoms as to make the problem uninteresting and the TDSE sufficiently small.

3

Method

The calculations are set up so that they will be working upwards, beginning with more basic methods such as finite differences, to the more advanced pseudo-spectral methods with split operators. The reason for this is that in this way a thorough comparison between the two methods can be made, and it will be made apparent why pseudo-spectral methods are more suitable for this type of problem, compared with finite differences.

First the Schr¨odinger equation with a known analytic solution, i.e with harmonic oscillator poten-tial, is solved with finite difference methods. Centered finite differences of different order of accuracy are used. The analytical eigenenergies are compared with the numerically computed eigenenergies. By comparing the change in relative error for different number of grid points, the order of accuracy and convergence can be evaluated. The order of accuracy is analyzed a second time by computing a reference solution and comparing the numerical solutions with the reference solution. When the order of accuracy and convergence are confirmed, the method is used to solve the Schr¨odinger equation with Morse potential.

The methods for solving the TDSE with a time-dependent potential are then analyzed. A centered finite difference of order eight is used for solving in space. Here there is no known analytic solution, so the method’s convergence is assured by solving the problem with a decreased step size in space and a small time step. If the solution converges the method is assumed to converge. In order to analyze the order of accuracy a reference solution is computed and the numerical solutions are compared with the reference solution. The computational time, computational work and memory requirements are noted for future comparison.

4

Results

4.1 Model problem

Using the harmonic oscillator potential for the time-independent model problem the potential’s eigenen-ergies were computed both analytically and numerically with finite differences in order to investigate order of accuracy of the centered difference. The results and the relative errors are shown in Table 3 and 4 for a second and fourth order method respectively. It can be seen that the relative error decreases with a factor 22 _{when doubling the number of grid points for the second order method, and}

with a factor 24for the fourth order method. It can also be seen that for all methods, the relative error increases for higher eigenenergies for the same number of gridpoints. When doing the same analysis for a sixth and eighth order method it is found that the error decreases with a factor 26 _{and 2}8_.

Table 3: Comparison of eigenenergies, second order method. n Analytic Number of gridpoints Second order

Eigenenergy Eigenenergy Relative error %

0 0.5 25 0.4807 3.8 50 0.4951 0.98 100 0.4988 0.24 1 1.5 25 1.4001 6.7 50 1.4755 1.6 100 1.4938 0.41 2 2.5 25 2.2193 11.2 50 2.4358 2.6 100 2.4840 0.64

Table 4: Comparison of eigenenergies, fourth order method. n Analytic Number of gridpoints Fourth order

Eigenenergy Eigenenergy Relative error %

0 0.5 25 0.48968 0.64 50 0.4998 0.04 100 0.5000 0 1 1.5 25 1.4781 1.46 50 1.4984 0.11 100 1.4999 0.007 2 2.5 25 2.4208 3.2 50 2.4942 0.23 100 2.4996 0.02

Table 5: The slope of the curves in figure 4 corresponding to the order of accuracy of the method. Order of accuracy Slope

2 1.9

4 3.8

6 4.8

8 7.3

Figure 4: Discretisation error in the numerical methods depending on the size of the time step.

The solution to the TDSE with harmonic oscillator potential solved with an eighth order centered finite difference after 0, 1, 3.4 and 4.5 atomic units of time are shown in Figure 5. The black curve is the initial value, and the other curves show the solution oscillating back and forth. Using the analytic solution to the ground state as initial value causes the solution to oscillates vertically.

The solution to the TDSE with Morse potential solved with an eighth order centered difference is shown in Figure 6. The solution oscillates in this case as well, but it displays a more oscillatory behavior the longer the simulation progress.

Figure 6: Numerical solution to the TDSE using a time-independent Morse potential.

The TDSE with a time-dependent perturbation to the harmonic oscillator potential was solved with an eighth order centered finite difference and simple time-stepping. Figure 7 shows the discretisation error as a function of step size. It can be seen that the error decreases with decreased step size. The curve has a slope of 7 in the steepest part and a slope of 2 for smaller steps.

Figure 7: Error depending on step size when the TDSE with a time-dependent harmonic oscillator potential was solved with an eighth order centered difference.

Figure 8: The movement of a particle in a harmonic oscillator potential. The upper graph shows the time-independent case and the lower the time-dependent case.

Tabulated in table 6 is the runtime required for solving the Schr¨odinger equation with both a time-independent and time-dependent Morse potential, i.e. the more realistic model problem, with finite differences and the simple method for time-stepping. As can be seen, this takes a very long time, over 2 hours in total. Table 6 also shows the allocated memory required to run the simulation.

Table 6: Runtime and memory required to solve the Schr¨odinger equation with both a time-independent and a time-dependent Morse potential, using an 8th order finite difference method with N = 1024 and ∆t = 0.1.

Potential Runtime [s] Memory required [kb] Time-independent 19.7463 45140.00

Time-dependent 7.53 · 103 4.95 · 108

Table 7 shows the data for solving the Schr¨odinger equation with both a independent and time-dependent Morse potential using the PSM with split-operators. Compared with the data presented in table 6 these methods are evidently much faster than using finite differences with simple time-stepping.

Table 7: Runtime and memory required to solve the Schr¨odinger equation with both a time-independent and a time-dependent Morse potential, using the PSM with split-operators, with N = 1024 and ∆t = 0.1.

Potential Runtime [s] Memory required [kb] Time-independent 0.473658 5304.00

4.2 NaI-system

Figure 9 and 10 shows how the populations behave for two different sized grids when the wavelength of the laser pulse is 360 nm. For N = 512 the fluctuations of the population on the excited state cannot be seen.

Figure 9: The populations when solving with a grid of N = 512.

Figure 10: The populations when solving with a grid of N = 2048.

Figure 11: The intensity of the populations over time, where the left graph shows the ground state and the right the excited state, for λ = 360 nm.

Figure 12 and 13 show the same calculations as in figures 10 and 11, but with a laser pulse of higher energy. Here the change in the populations is more distinct. It can be seen that the probability that the molecule will be moving on the excited potential surface is much higher than in the case with λ = 360 nm. In the right graph in Figure 13 it can be seen that there is a probability that the molecule dissociates into two ions. This can bee seen since there is a probability for an increased inter molecular distance for the excited state.

Figure 13: The intensity of the populations over time, where the left graph shows the ground state and the right the excited state, for λ = 335 nm.

Figures 14 and 15 once again shows the behaviour of the populations, but this time for a wavelength of 310 nm. Here the population on the excited state is less intense than in the case with λ = 335, even though λ = 310 translates to a higher energy perturbation. There is still a probability that the molecule dissociates into two ions, even if it is smaller than the previous case.

Figure 15: The intensity of the populations over time, where the left graph shows the ground state and the right the excited state, for λ = 310 nm.

5

Discussion

The three first eigenenergies for a harmonic oscillator were computed numerically with centered finite differences of different order of accuracy. Since the eigenenergies converged to the analytic values for all methods, all methods were assumed to converge. The relative error in the eigenenergies confirmed the order of accuracy for a second, fourth, sixth and eighth order method. The error decreased with a factor 22_{, 2}4_{, 2}6 _{and 2}8 _{respectively when doubling the number of grid points, which was expected.}

The error decreased with the same factor for all three eigenenergies, but the errors were larger for the third and second eigenenergies than for the first. This is due to that the eigenvectors contain higher frequencies that are badly represented on the discretisation grid. Looking at the slopes of the log-log plot over discretisation error as a function of step size, tabulated in table 5, is another way to determine the order of accuracy. The slopes almost correspond to the theoretic order of accuracy. It can be seen in the plot that the error goes to zero for smaller time steps for all methods but at a certain point the slope is no longer as steep as before. This especially occurs for the eighth order method for small steps. This is because the errors are of the same order of magnitude as the truncation errors.

The wave function in an harmonic oscillator potential oscillated as expected. However, when in-serting the analytic ground state for a harmonic oscillator potential as initial value, the numerical solution oscillated. This happens because the exact analytic solution is not the exact solution to the discretized model, due to the error in the discretized model. According to the numerical model the initial value is not the ground state and the solution will therefore oscillate. To avoid this, the numer-ical solution to the ground state should be used instead. In the Morse potential the wave moved more irregularly. This is not due to discretization errors, instead it is a real property of the anharmonicity of the Morse potential. If one would simulate long enough the solution would once again return to it’s smooth appearance.

smaller steps. Since the slope corresponds to the order of accuracy it is probable that the accuracy in time affects the accuracy of the whole solution for small steps. The wave function did harmonic oscillations in the case with time dependent potential as well as shown in figure 8. The difference was that the oscillations moved in a cosine pattern, just like the time dependent perturbation. It was expected that the time dependent potential should affect the wave to move in the same pattern as the perturbation.

In table 8 the data from tables 6 and 7 have been compiled to show a comparison of run time, memory requirements and number of calculations for the methods. It is evident that PSM with split-operators are superior to finite differences with simple time-stepping in every respect. Solving using PSM with split-operators is more than ten thousand times faster than performing the same simulation with finite differences and simple time-stepping, and requires about a thousand times less memory.

Table 8: Comparison of runtime, memory requirements and number of operations for finite differences with simple time-stepping and PSM with split-operators, when solving the Schr¨odinger equation with a Morse potential, both time-independent and time-dependent. N = 1024, ∆t = 0.1.

Finite differences Pseudo spectral Time-indep. Time-dep. Time-indep Time-dep. Runtime [s] 19.7463 7.53 · 103 0.473658 0.59669 Number of operations/time step O(N2_{) O(N}3_{) O(N logN ) O(N logN )}

Memory usage [kb] 45140.00 4.95 · 108 5304.00 1.1 · 105

Comparing figure 9 and 10, it can be seen that to be able to see the behaviour of the population on the excited potential surface, a large enough grid is needed. 512 gridpoints is too small, but 1024 gridpoints works well enough. However, for all other simulations of the TDSE with coupled potentials a grid of 2048 points were used.

A laser pulse with a wavelength corresponding to the difference in energy between the ground- and excited states can be seen to produce a larger population transfer. Thus the molecule can with a higher probability be excited with a laser pulse with the correct wavelength, duration and intensity. It can also be seen that there is a probability that the molecule dissociates into two ions for certain wavelengths. This can be seen by comparing figures 11, 13 and 15. A laser pulse with either λ = 360 nm or λ = 310 nm give roughly the same population intensity on the excited state, while λ = 335 nm gives a noticeably more intense population on the excited state. For λ = 360 nm there is almost no probability that the molecule dissociates, while for λ = 310 nm there is a small probability and for λ = 335 nm there is a greater probability. This is because of the quantized levels of energy in the molecule making λ = 335 nm simply fit better.

6

Conclusions

Finite differences paired with a simple time-stepping method such as expm() can be used to solve simpler models of the TDSE, with adequate results. However, when working with more complex mod-els that include a time-dependent Hamiltonian these simple numerical methods become hopelessly demanding in regards to compuational time, if one wants accurate results. For these problems it is much more suitable to use the pseudo-spectral method with split-operators, which is remarkably faster when comparing with the simpler methods.

References

[1] T. S. Rose, M. J. Rosker, A. H. Zewail. Femtosecond realtime probing of reaction. IV. The reactions of alkali halides. J. Chem. 1989 Dec; 91(12): 7415-7436. DOI:10.1063/1.457266

[2] Tannor. Introduction to Quantum Mechanics: A Time-Dependent Perspective. Sausalito: Uni-versity Science Books; 2007.

[3] Fornberg B. Generation of Finite Difference Formulas on Arbitrarily Spaced Grids. Math Comput. 1988 Oct; 51(184): 699–706. DOI:10.1090/S0025-5718-1988-0935077-0

[4] Gasiorowicz S. Quantumm Physics. Third edition. Hoboken, NJ: John Wiley & Sons, Inc.; 2003 [5] Kormann K, Holmgren S, Karlsson H. Accurate time propagation for the Schr¨odinger equa-tion with an explicitly time-dependent Hamiltonian. J Chem Phys. 2008 May; 128(18). DOI:10.1063/1.2916581.

[6] Magnus W. On the Exponential Solution of Differential Equations for a Linear Operator. Commun Pure Appl Math. 1954 Nov; 7(4): 649-673. DOI:10.1002/cpa.3160070404

[7] Strang, G. On the Construction and Comparison of Difference Schemes. SIAM J Numer Anal. 1968 Sep; 5(3): 506-517. DOI:10.1137/0705041

[8] Schwendner P. Seyl F. Schinke R. Photodissociation of Ar₂+in strong laser fields. J. Chem. Phys. (1997); 217: 233-247.

A

Parameters

Here the values for the parameters that have been used for the simulations are tabulated, scaled according to atomic units.

Table 9: Values for the parameters used for the Harmonic oscillator, equation 5. Value [au]

ω 0.05 r0 2

α 2

Table 10: Values for the parameters used for the Morse potential, equation 7. Value [au]

De 12440

a 1.875 r0 2.656

Table 11: Values for the parameters used for the potentials used for the NaI-system, described in equations 8, 10 and 12.