Finite-difference time-domain method for calculating photoexcited molecular vibration wave-packet dynamics

We propose a simple numerical calculation method for solving molecular vibrational dynamics. A finite-difference time-domain (FDTD) method is utilized to solve the Schrödinger equation for wave packet dynamics driven by incident light pulses. Vibrational wave function is discretized and divided into real and imaginary parts, and then the Schrödinger equation is solved by temporally and alternately calculating the real and imaginary parts, in a similar manner to the conventional FDTD method. Taking generation of a hot molecule by pump-dump scheme and multi-photon excitation process as examples, we show the validity and usefulness of the proposed method. The dynamics of the vibrational wave packets in the pump-dump scheme and transient virtual excitations of intermediate vibronic state in the three-photon absorption can be easily visualized, and the intensity dependence of excitation population for ultra-high intensity beyond the perturbation theory can be accurately calculated.


Introduction
Coherent control is a technique for controlling molecular vibrational wave packets using ultrashort light pulses [1] and has been intensively investigated in various fields of physics, chemistry and biology, combined with newgeneration of techniques, such as quantum information technologies. In particular, the coherent control of photochemical reactions has a high potential for directly controlling the elementary processes of the photochemical reactions, such as ionization [2], photodissociation [3] and photosynthetic energy transfer [4]. So far, various methods of coherent control have been proposed [5][6][7] and demonstrated experimentally for simple molecular systems. For a number of complex molecules of interest, however, we have to adequately control irradiation timing and shape of pulses, and therefore the prediction using numerical calculation is highly important.
Basically, for the numerical analysis of the quantum dynamics of molecular wave packets, we only have to solve the Schrödinger equation as an initial-value problem. In fact, some numerical calculation methods such as the Crank-Nicolson method and the split-step (Fourier) method have been proposed [8,9]. However, to implement these methods, expertise of quantum theory and knowledge about computational algorithms, especially for stability and boundary conditions, are required, and these could be a barrier for non-experts. In addition, from the viewpoint of multi-photon excitation processes driven by multiple pulse irradiations, the prediction of proper timing of pulse irradiation will be of increasing significance, combined with femtochemistry [10] and rapidly developing attochemistry [11], and hence the demand for a simple numerical calculation will increase.
One of promising candidates for the simple numerical calculation of wave packet dynamics is a finitedifference time-domain (FDTD) method. The FDTD method is a numerical algorithm introduced originally to calculate the propagation of electromagnetic waves with comparatively high stability and simple implementation, in which Maxwell's equations are solved as an initial-value problem by temporally and alternately calculating the electric and magnetic fields. In fact, using the FDTD algorithm, we can simply implement a computer program by directly utilizing the Maxwell's equation except for spatiotemporal discretization. This is the reason why the FDTD calculation is now one of the most widely used tools in the analysis of electromagnetic waves, and the packaged software for the FDTD calculation is now readily accessible to everyone.
An FDTD method for electronic systems has already been proposed [12][13][14], called the FDTD method for quantum device (FDTD-Q). In the FDTD-Q method, electronic wave function is discretized and divided into real and imaginary parts, ψ(r)=ψ R (r)+iψ I (r), and then the Schrödinger equation is solved by temporally and alternately calculating ψ R (t) and ψ I (r), in a similar manner to the FDTD method. However, the application of the FDTD-Q method is restricted to the electronic quantum devices, such as quantum dots and quantum wells, and to eigenmode analysis for the quantum devices. In addition, in the conventional FDTD-Q method, electronic transitions by light irradiation are not considered. The advantage of the FDTD algorithm is that we can calculate the dynamics of wave propagation, and therefore the above applications do not make full use of the characteristic of the FDTD method. These restrictions reduce the versatility of application of the FDTD-Q method.
In this study, we develop a new FDTD method for simply calculating the wave packet dynamics driven by light pulse excitation (hereinafter referred to as FDTD-WD method) and apply the FDTD-WD method to photoexcited molecular vibrational motions. Taking high-energy vibrational wave packet generation using pump-dump scheme [15] and multi-photon absorption process as examples, we show the validity and usefulness of the FDTD-WD method. Using the FDTD-WD method, we can easily calculate wave packet dynamics only by providing adiabatic potentials and incident light pulses without molecular eigenmode analysis or Franck-Condon analysis.
The rest of this paper is organized as follows. In section 2, the formulation of the FDTD-WD method is given in sequence. In section 3, we apply the FDTD-WD method to two concrete systems, namely, the generation of a hot molecule using a pump-dump scheme and three-photon absorption process, and show the validity and usefulness of the FDTD-WD method. In section 4, we summarize and discuss our results.

Schrödinger equation for vibrational motion driven by light pulses
In this section, we derive in sequence the Schrödinger equation describing molecular vibration driven by incident light pulses. Generally, electronic transition by light-molecule interaction takes place much faster than the timescale of nuclear motion, and therefore the motions of nuclei and electrons can be separated on the assumption of the so-called Born-Oppenheimer approximation. The Hamiltonian for the nuclei is then given by where M i and X i are the mass and position of the i-th nucleus, respectively. ÿ=h/2π is the Dirac constant, ε n is the adiabatic potential of the n-th electronic level and { } X denotes  { } X X , , 1 2 . In order to focus on molecular vibration, we rewrite equation ( whereL i 2 is the square angular momentum operator, related to rotational motion, and is ignored in this study.
Consequently, we can obtain the Hamiltonian for the molecular vibration as follows, In general, the spatial scale of molecules is quite small in comparison with that of the wavelength of light and therefore for light-molecule interactions we can approximate the spatial distribution of the optical electric field as a constant, which is the so-called long-wave approximation. In addition, in this study, we consider only electric dipole transitions, for simplicity. After the second quantization of electronic states, the dipole operator describing the transition between the electronic states of ñ |n and + ñ |n 1 is given by ñá + + + ñá (| | | |) d n n n n 1 1 n , where d denotes the electric dipole moment. Using the optical electric field ( ) E t at the position of the molecule, the Hamiltonian for the light-molecule interaction can be described as The vibrational wave packet in the n-th adiabatic potential can be described as a superposition state of the direct product of ñ |n and the ν-th vibrational mode nñ | n , given by where n c n, is the probability amplitude. The Schrödinger equation for y ñ | ( ) t n at time t can then be given by The wavefunction for vibrational motion is obtained by operating á á | { }| n r in equation (8). The Schrödinger equation for a vibrational wave function in the n-th adiabatic potential can be described as 1 , related to the Rabi frequency, and we use y In accordance with the standard method in quantum mechanics, we replace y to simplify the description of equation (9). With this replacement, the Schrödinger equation for a vibrational wave packet Ψ n is given by The normalization is taken as . The last two terms in equation (10) indicate the electronic transitions between the adiabatic potentials induced by light irradiation.

Formulation of the FDTD-WD method
In this section, we transform equation (10) into the FDTD formulation. First, we discretize the time t and relative distance r i . Using the time division of Δ t and the spatial division of Δ i , the vibrational wavefunction can be rewritten as where the same spatial division Δ is used for all relative distances r i , for simplicity. The time derivative in equation (10) can be approximated by forward difference, as The spatial second derivative can be rewritten using the central difference, as where derivative is performed only for the i-th relative distance.
We can now formulate the FDTD algorithm. According to [12], the FDTD algorithm for Schrödinger equation is solved by separating Ψ into the real and imaginary parts as Y = Y + Y i R I . The essential point is to alternately calculate Ψ R (t) and Ψ I (r), by introducing the intermediate . The FDTD algorithm for the vibrational dynamics in n-th adiabatic potential is then given by The key is to calculate Y + ℓ n I , , , and hence we can program the FDTD-WD algorithm using equations (14) and (15) as is. The steps for a concrete numerical loop iteration are as follows:

Increment ℓ by one and go back to 2
Thus, we can calculate vibrational wave packet dynamics using the above FDTD-WD method. We only have to provide the initial condition of Ψ, ε n and Ω. In particular, for e n , we can utilize ab initio quantum chemical calculation and the approximate solution, e.g., the Morse potential.
In general, numerical analysis of spatial regions requires the boundary conditions at the ends of the analytical domain. The most famous one is the absorbing boundary condition as typified by the perfectly matched layer (PML). However, for the molecular vibration considered in this study, we need not consider the complicated PML when widening the analytical domain, because the molecular vibration is bound and usually localized in the adiabatic potential.
However, in order to stably execute the FDTD-WD calculation, the Courant condition is required for Δ t and Δ r . The Courant condition for equations (14) and (15) is given by [16]   When we choose small Δ, spatial resolution becomes high. However, Δ t must also be small to satisfy the Courant condition, leading to vast loop iteration and calculation time. Conversely, large Δ t can reduce the calculation time, but then the spatial resolution becomes low. Thus, the trade-off between calculation time and spatial resolution always exists and we have to choose proper values of Δ t and Δ. In the next section, we will analyse a method to select Δ.

Adiabatic potentials: Morse potential approximation for diatomic molecule
In this study, we adopt an diatomic molecule and the Morse potential [17] as an adiabatic potential, for simplicity, because a simple two-body system is fundamental for understanding anharmonic molecular vibration and furthermore we can analytically solve the eigenvalues and eigenfunctions from the Morse potential.
For the two-body system, it is convenient to utilize the relative coordinate =r X X where L represents the generalized Laguerre polynomials, and N is the normalized coefficient, given by where Γ is the Gamma function. In order to introduce the Morse potential into the FDTD-WD method, we have to discretize the spatial variable r using Δ. In addition, in actual calculations, we also have to choose the value of Δ, the range of analytical domain, and the cut-off energy of V(r) for small r to avoid the divergence of  . The eigenfunctions for ν=0, 10, 20, 30, 40 and 50 are plotted together at the corresponding eigenenergies in the potential. Even for ν=50, the vibrational mode is well localized and the probability amplitude is zero at the ends of the analytical domain. Since a vibrational wave packet in the potential is formed as a superposition state of the eigen modes, the absorbing boundary condition is not needed for the vibrational FDTD-WD calculation having wide analytical domain. Figure 1  To appropriately approximate the wave properties using the FDTD method, it is said that at least 20 points per one oscillation are needed. The wave function for ν=50 has 50 nodes, oscillating 25.5 times spatially, and is spatially distributed from a 6 B to a 13 B . Therefore, from rough estimation, we can obtain B , and this is in good agreement with our numerical result shown in figure 1(b). Though we use the Morse potential as an adiabatic potential in this study, the Rydberg-Klein-Rees (RKR) method [20] or ab-initio quantum chemical calculations can also be used. In such cases, the above rough estimation is useful to determine the value of Δ.

Numerical results
In this section, we apply the FDTD-WD method to two concrete systems, namely, vibrational wave packet generation using the pump-dump scheme [15] and three-photon absorption, and visualize their transient processes to validate the FDTD-WD method.

Vibrational wave packet generation using pump-dump scheme
The pump-dump scheme is one of the techniques to generate a high-vibrational wave packet in the electronic ground state. The concrete process is described as follows: we first generate a vibrational wave packet in the excited state using a pump pulse, and then dump the wave packet to high vibrational levels in the ground state using a dump pulse with a higher intensity and an energy lower than the pump pulse, as illustrated in figure 2(a). As an analytical mode, we adopt a diatomic molecule K 2 , similarly in figure 1, where S + X g 1 and P B u 1 are used as a ground state and an excited state, respectively. The molecular parameters are the same as in figure 1 for the S + X g . In this study we use Gaussian pulses, given by , where ξ[y(r)] is the eigenfunction for the ground state, defined by equation (19). The results show that two local maxima exist at τ=190 and 300 fs. This indicates that proper adjustment of τ is needed to achieve an efficient pump-dump scheme. Thus, the FDTD-WD method is useful for visualizing the wave packet dynamics and for estimating properτ.

Three-photon absorption
In this section, we apply the FDTD-WD method to multi photon excitation. Multi photon excitation is a highorder nonlinear optical process. In particular for more than three photons, theoretical analysis becomes complicated and numerical calculation is highly important. In this study, we calculate a three-photon absorption process in a diatomic molecule with four sets of vibronic states, as depicted in figure 4(a). For the vibronic states, S + X g by an intense single pulse. The central wavelength of the pulse is set to l = 885.6 nm, and other parameters except for λ are the same as the pump pulse in figure 2. The incident pulse is linearly polarized and has a Gaussian shape, given by . The molecular parameters for S + X g  , however, the molecule exhibits rapid Rabi oscillation during interacting with incident intense pulse as shown in figure 5(b) and ρ e derails from the line predicted by the perturbation theory. Thus, the FDTD-WD method accurately reproduces the perturbation theory for low intensity and can calculate the Rabi oscillation in high intensity, which cannot be predicted using the perturbation theory.
It should be mentioned that actual molecules undergo structural changes for ultra-high intensity (∼TW cm −2 ) considered in figure 5. Therefore, we might not give a practical meaning to the calculation result of Rabi oscillation. However, figure 5 is given to show that we can easily calculate the optical response of molecules  to ultra-high intensity beyond the perturbation theory, using the FDTD-WD method. If we consider a molecule with large dipole moment (10 Debye), the Rabi oscillation can occur for lower intensity (∼100 GW cm −2 ) and could be actually observed. In fact, for molecules confined in a microcavity, the vacuum Rabi oscillation has already been observed experimentally utilizing the cavity QED effect [21]. Finally, we briefly refer to the advantage of the FDTD-WD, compared with available standard methods. The important key in numerically calculating the Schrödinger equation is to conserve the energy and probability amplitude. In the Crank-Nicolson method, they are kept by combining the implicit and explicit methods of a half time step Δ t /2, e.g., . This equation can be reduced to a system of equations and we can obtain molecular wave packet dynamics by solving the system of equations at each step. The Crank-Nicolson method is stable and very useful for low dimension and small time steps, however the calculation becomes complicated and heavy for high dimension and large time steps because we must solve the system of equations at each step. On the other hand, in the split-step (Fourier) method, the time evolution operator is split into two parts of kinetic and potential evolutions as , where + =ˆT U H, andT and U are kinetic and potential evolutions, respectively. In calculating the equation, we have to Fourier transform the part of kinetic evolution to wave-number domain in order to convert the partial derivative of ∂/∂x inT to ik. Therefore, we have to compute fast Fourier transform at each step, and hence much computation time becomes required for high dimension and large time steps for the same reason in the Crank-Nicolson method. However, the proposed FDTD-WD method does not require other calculation algorithms to compute the time evolution, in contrast to the above available methods. We only have to compute equations (14) and (15) even for higher dimension. In addition, the leap-frog alternating calculation of the real and imaginary parts of wave function ensures the energy and probability conservations. Thus, the FDTD-WD method has advantages in simplicity and convenience of program implementation, which just meet the original purpose of this study.

Summary and discussions
We have proposed the FDTD-WD method for numerically solving Schrödinger equation for molecular vibrations driven by light pulse irradiation. Taking the generation of a hot molecule using pump-dump scheme and three-photon absorption in a diatomic molecule as examples, we have demonstrated the validity of the FDTD-WD calculation. The optical transition of vibrational wave packets in the pump-dump scheme and virtual excitations of intermediate vibronic state in the three-photon absorption process can be easily visualized. In addition, the intensity dependence of population for ultra-high intensity beyond the perturbation theory is accurately calculated. The FDTD-WD method is thus a very useful tool for numerical analysing and visualizing the transient processes in photoinduced chemical reaction.
Throughout this study, we focus on the excitation of relatively low-energy vibrations and hence the absorbing boundary condition can be ignored. However, if we target high-energy vibrations around the highestenergy eigenmode, the absorbing boundary condition must be considered. One possible way is to adopt a scheme of perfectly matched layer for Schrödinger equation, proposed in [22]. Though this scheme is originally developed for quantum device application, it is possible to apply it to the FDTD-WD method.
Finally, we refer to the introduction of relaxation effects in the FDTD-WD method. For general molecular systems, the relaxation due to photon emission can be ignored because the time scale of the emission process is much longer than that of ultra-short pulse excitation. For phase and thermal relaxation, however, the time scale is comparable to that of the vibrational process. Therefore, except for low-temperature environment, we have to consider the relaxation effects. A straightforward way of including the phase and thermal relaxation effects is to calculate the dynamics using quantum master equation. However, to calculate the quantum master equation, some expertise of quantum mechanics is required, and this does not meet the original purpose of this study. One promising candidate is the stochastic Schrödinger equation [23] because the scheme can directly applied to the FDTD-WD method proposed in this study. This formulation will be our next step.