Toward fully quantum modelling of ultrafast photodissociation imaging experiments. Treating tunnelling in the ab initio multiple cloning approach

We present an account of our recent e ﬀ ort to improve simulation of the photodissociation of small heteroaromatic molecules using the Ab Initio Multiple Cloning (AIMC) algorithm. The ultimate goal is to create a quantitative and converged technique for fully quantum simulations which treats both electrons and nuclei on a fully quantum level. We calculate and analyse the total kinetic energy release (TKER) spectra and Velocity Map Images (VMI), and compare the results directly with experimental measurements. In this work, we perform new extensive calculations using an improved AIMC algorithm that now takes into account the tunnelling of hydrogen atoms. This can play an extremely important role in photodissociation dynamics.


I. Introduction
Quantum non-adiabatic molecular dynamics is a powerful tool for understanding the details of the mechanisms of important photo-induced processes, such as the photodissociation of pyrrole and other heteroaromatic molecules.In these processes, quantum effects such as electronically non-adiabatic transitions and tunnelling are important, and an approach that goes beyond surface hopping, such as multicongurational time dependent Hartree (MCTDH), 1 for example, is oen required.MCTDH can be very accurate, and was recently used to simulate the dissociation of pyrrole. 2 However it needs a parameterized potential energy surface as a starting point, which signicantly restricts its practicality.A good alternative is represented by a variety of methods [3][4][5][6][7][8][9][10][11] based on trajectory-guided Gaussian basis functions (TBF).Despite the fact that such approaches use classical trajectories, they are still fully quantum mechanical because these trajectories are employed only for propagating the basis, while the evolution of their amplitudes and, thus, of the total nuclear wave-function is determined by the time-dependent Schrödinger equation.An important advantage of trajectoryguided quantum dynamics methods is that they are fully compatible with direct or ab initio molecular dynamics where excited state energies, gradients, and nonadiabatic coupling terms are evaluated on the y simultaneously with the nuclear propagation.The disadvantage is that trajectory based direct dynamics is very expensive due to the high cost of electronic structure calculations and typically can afford only a limited number of trajectories, which can be an obstacle to full convergence.
Recently, we introduced the ab initio multiple cloning (AIMC) 10 method, where TBFs are moving along Ehrenfest trajectories, as in the multicongurational Ehrenfest (MCE) 8,9 approach, with bifurcation of the wave-functions taken into account via basis function cloning.While leading to the growth of the number of trajectories, the use of cloning helps to adopt the basis set to quantum dynamics signicantly better than in the classical MCE approach.AIMC also uses a number of tricks to efficiently sample the trajectory basis and to use the information obtained on the y: (1) similar to the previously developed trajectory based methods AIMC relies on importance sampling of initial conditions.(2) AIMC uses the so called time displaced or train basis sets, 10,12,13 which increase the basis set size almost without any additional extra cost by reusing the ab initio data which has already been obtained.(3) The method calculates quantum amplitudes in a "post-processing technique" aer the trajectories of the basis set functions have been found.As a result, the trajectories can be calculated one by one in parallel and good statistics can be accumulated.
In this work, we present a new implementation of the AIMC approach that is improved to take into account the tunnelling of hydrogen atoms by identifying possible tunnelling points and placing additional TBFs of the other side of the barrier.We use this new implementation to simulate the dynamics of the photodissociation of pyrrole, a process where tunnelling can play a very important role.We calculate the TKER spectrum and velocity map image (VMI), and directly compare the results of our calculations with experimental observations. 14he paper is organized as follows.In Section II we describe the proposed implementation of the AIMC approach.Section III contains the computational details of our simulations.In Section IV, we present and discuss the results.Conclusions are given in Section V.

II.1 Working equations
The AIMC method 10 is based on the same ansatz as the multicongurational Ehrenfest (MCE) approach, 8,9 in which the total wave-function |J(t)i is represented in a trajectory-guided basis |j n (t)i: The basis functions |j n (t)i are composed of nuclear and electronic parts: The nuclear part |c n (t)i is a Gaussian coherent state moving along an Ehrenfest trajectory: where R n (t) and P n (t) are the phase space coordinate and momentum vectors of the basis function centre, g n (t) is a phase, and the parameter a determines the width of the Gaussians.The electronic part of basis functions |j n (t)i is represented as a superposition of several adiabatic eigenstates |f I i with quantum amplitudes a (n) I .The time dependence of the Ehrenfest amplitudes a (n) I is given by the equations where the matrix elements of electronic Hamiltonian H el(n) IJ are expressed as: here The motion of the centres of the Gaussians follows the standard Newton's equations: where the force F n is an Ehrenfest force that includes both the usual gradient term and the additional term related to the change of quantum amplitudes as a result of non-adiabatic coupling: Finally, the phase g n evolves as: Eqn ( 3)-( 8) form a complete set, determining the basis and its time evolution.The evolution of the total wave-function |J(t)i (eqn (1)) is dened by both the evolution of the basis functions |j n (t)i and the evolution of the relevant amplitudes c n (t).The time dependence of the amplitudes c n (t) is given by the equation Assuming that the second derivative of the electronic wave-function |f I i with respect to R can be disregarded, we get: The matrix elements of the kinetic energy operator be calculated analytically.For potential energy and non-adiabatic coupling matrix elements, we use a simple approximation: The approximation ( 12) represents a linear interpolation of the potential energy between the two points and can be improved further at the cost of calculating higher derivatives of the potential energy along the trajectories.It has been tested previously, 10 and no visible change of the results was found when this approximation was applied compared to the saddle point approximation which expands around a distinct centroid for each pair of TBFs. 4 The term hj m ðtÞj d dt jj n ðtÞi in eqn ( 9), which originates from the time dependence of the basis, can be expressed as: where Notice that in the AIMC approach, all off-diagonal matrix elements entering eqn (9) are calculated from the electronic structure data at the TBF centres, which is needed for the propagation of the basis.Thus, quantum coupling between the congurations comes at almost no extra cost.Moreover, eqn (9) can be solved aer the trajectories have been calculated, provided the appropriate electronic structure information has been saved.
The detailed derivation of MCE equations together with the expressions for relevant matrix elements can be found in our previous works. 10,11

II.2 Basis set sampling and cloning
The Ehrenfest basis set is guided by an average potential, which can be advantageous when quantum transitions are frequent.However, it becomes unphysical in regions of low non-adiabatic coupling when two or more electronic states have signicant amplitudes: in this case, the difference of the shapes of potential energy surfaces for different electronic states should lead to branching of the wavepacket.
In order to reproduce the bifurcation of the wave-function aer leaving the non-adiabatic coupling region, AIMC methods adopt the cloning procedure, 10 where the appropriate basis function is replaced by two basis functions, each guided (mostly) by a single potential energy surface.Aer the cloning event, an Ehrenfest conguration and j The rst clone conguration has non-zero amplitudes for only one electronic state, and the second clone contains contributions of all other electronic states.The amplitudes of the two new congurations become: so that the contribution of the two clones |j 0 n i and |j 0 n 0 i to the whole wavefunction (1) remains the same as the contribution of original function: We apply the cloning procedure shortly aer a trajectory passes near a conical intersection, when the non-adiabatic coupling is lower than a threshold, and, at the same time, the so-called breaking force which is the force pulling the Ith state away from the remaining states, is sufficiently strong.The cloning procedure is very much in spirit of the spawning, used in the Ab Initio Multiple Spawning approach (AIMS).Cloning does not require any backpropagation of spawned/cloned basis functions, unlike many 4 (but not all 15,16 ) implementations of spawning.
As has been described in our previous work, 7 we rely on importance sampling when generating the initial conditions.Using the linearity of the Schrödinger equation, we rst represent the initial wave-function as a superposition of Gaussians and then propagate each of them independently, "bit-by-bit". 7We use a time-displaced basis set (coherent state trains), where several Gaussian basis functions are moving along the same trajectory but with a time-shi Dt, allowing us to reuse the same electronic structure data for each of the basis functions in the "train."Fig. 1 shows a time displaced basis guided by a trajectory and its bifurcation via cloning.The best possible result with AIMC can be achieved when a swarm of trains is used to propagate each "bit" of the initial wave-function.

II.3 Tunnelling
The tunnelling of hydrogen atoms can play an important role in photodissociation processes.As mentioned above, MCE, AIMC and AIMS are fully quantum methods because classical trajectories are used only to propagate the basis, while the amplitudes c n (t) are found by solving the time dependent Schrödinger equation.When Gaussian basis functions are present on two sides of the potential barrier, the interaction between them can provide quantum tunnelling through the barrier.However, in the case of direct ab initio dynamics, the basis is usually very small, far from being complete.As a result, no basis functions normally would be present on the other side, and they must be placed there by hand in order to take tunnelling into account.
In this paper we adopt the ideas 17,18 previously used in the AIMS method to describe tunnelling for use with the AIMC technique.Fig. 2 illustrates the algorithm that we apply.First, we calculate the usual AIMC trajectories and nd turning points, where the distance between the hydrogen atom and the radical reaches a local maximum.Then, for each of these turning points, we calculate the shape of the potential barrier: we increase manually the length of the N-H bond keeping all other degrees of freedom frozen, calculate potential energies, and nd the point on the other side of the barrier with the same energy as in the turning point.If this point lies further than a set threshold from the turning point, we assume that tunnelling is not possible here, as the potential barrier is too wide.Otherwise, we use it as a starting point for an additional AIMC trajectory.The new trajectory is calculated both forward and backward in time, and the initial momenta are taken as the same as in the turning point, ensuring that new trajectories have the same total classical energies as their parent trajectories.This is exactly the procedure used in the multiple spawning approach, thus our method combines cloning for non-adiabatic events and spawning for tunnelling events.The forward propagation of new trajectories oen involves branching as a result of cloning; backward propagation is performed without cloning and for a sufficiently short time, until new and parent trajectories separate in phase space.When all the trajectories are calculated, we solve eqn (9) for quantum amplitudes c n (t) in a time-displaced basis set (coherent state trains).This is similar to our previous approach 10,11 but with the difference that now the basis is better adapted to treat tunnelling.The train basis on the new trajectory is placed in such way that it reaches the tunnelling point at the same time as the train basis on the parent trajectory.Because the new trajectory differs from its parent by only one coordinate at a tunnelling point, namely by the length of the N-H bond, there is a signicant overlap between Gaussian basis functions belonging to these two trajectories.This interaction is retained for a signicant time while the coherent state trains are passing the tunnelling point, ensuring the transfer of quantum amplitude across the barrier.

III. Computational details
Using our AIMC approach, we have simulated the dynamics of pyrrole following excitation to the rst excited state.Trajectories were calculated using the AIMS-MOLPRO 19 computational package, which has been modied to incorporate Ehrenfest dynamics.Electronic structure calculations were performed with the complete active space self-consistent eld (CASSCF) method using the cc-pVDZ basis set.As in our previous works, 9,11 we used an active space of eight electrons in seven orbitals (three ring p orbitals and two corresponding p* orbitals, one s orbital and a corresponding s* orbital).State averaging was performed over four singlet states using equal weights, i.e. the electronic wave-function is SA4-CAS(8,7)/cc-pVDZ.The width of Gaussian functions a was taken as 4.7 bohr À2 for hydrogen, 22.7 bohr À2 for carbon, and 19.0 bohr À2 for nitrogen atoms, as suggested in ref. 20.Three electronic states were taken into consideration during the dynamicsthe ground state and the two lowest singlet excited states.
The initial positions and momenta were randomly sampled from the ground state vibrational Wigner distribution in the harmonic approximation using vibrational frequencies and normal modes were calculated at the same CASSCF level of theory.We approximate the photoexcitation by simply liing the ground state wavepacket to the excited state, as would be appropriate for an instantaneous excitation pulse within the Condon approximation.Of course, the ne details of the initial photoexcited wavepacket are lost in this approximation, however, we do not expect these details to have much effect on the observables shown in this paper.
We have run 900 initial Ehrenfest trajectories, each propagated with a timestep of $0.06 fs (2.5 a.u.) for 200 fs or until the dissociation occurred, dened as an N-H distance exceeding 4.0 Å.For a small number of trajectories, simulations exhibiting N-H dissociation were carried out to the full 200 fs in order to investigate the dynamics of the radical.Cloning was applied to TBFs when the breaking acceleration of eqn (20) exceeded a threshold of 5 Â 10 À6 a.u. and the norm of the non-adiabatic coupling vector was simultaneously less than 2 Â 10 À3 a.u.For all initial trajectories, as well as for their branches resulting from cloning, we identied turning points for the N-H bond length and calculated the width of the potential barrier.Additional trajectories on the other side of the barrier were placed if the width of the barrier did not exceed 0.5 bohr, which corresponds to an overlap of $0.3 between Gaussian basis functions.The new trajectories were propagated backward for 20 fs to accommodate the train basis set, and forward until dissociation or until the trajectory time exceeds 200 fs.
For each initial trajectory with all its branches and tunnelling sub-trajectories, we solved eqn (9) using a train basis set of N ¼ 21 Gaussians per branch, separated by 10 time steps, which corresponds to an average overlap of $0.6 between the nearest Gaussians in the train.The total size of the basis is constantly changing because of the inclusion of new branches.The nal amplitudes c n give statistical weights for each of the branches, which are used in the analysis that follows.

IV. Results
As a result of cloning, 900 initial congurations give rise to 1131 trajectory branches.This corresponds to an average of $0.25 cloning events per initial trajectory.For these branches, we have found 7702 local maxima of N-H bond length, of which 2376 have been identied as possible tunnelling points.For all these points, we run sub-trajectories, which nally gives 3203 additional branches, 4334 branches in total.The majority of these branches undergo N-H dissociation within our computational time of 200 fs: the total statistical weight of dissociative trajectories is 92%, of which 53% is the contribution of tunnelling sub-trajectories.The kinetic energy distribution of the ejected hydrogen atom is presented in Fig. 3 together with the experimental TKER spectrum. 14Both distributions clearly exhibit two contributions: a large peak at higher energies, and a small contribution at lower energies.It is important to note that adopting the basis set to tunnelling shis the high-energy peak of TKER spectrum toward the lower energies by about $1000 cm À1 and makes the low-energy peak slightly more pronounced.While the calculated energies are still on average about 1.5 times higher than experimental values, this difference can be ascribed to the lack of dynamic electron correlation in the CASSCF potential energy surfaces.We previously showed 11 that a more accurate MS-CASPT2 PES would lead to a shi in the kinetic energy peak of approximately $1800-1900 cm À1 towards lower energies, signicantly improving the agreement with experiment.
Analysis of the electronic state amplitudes in the Ehrenfest congurations (eqn ( 2)) shows that the bifurcation of the wave-function while passing through a conical intersection plays an important role in the formation of a two-peak spectrum: the high kinetic energy product is predominantly in the ground state, while the low energy peak is formed by mostly low-weight branches with substantial contribution from excited electronic states.Fig. 4 presents an example of such a bifurcating trajectory.At about 55 fs aer photoexcitation, this trajectory reaches an intersection for the rst time.Aer passing the intersection, the ground and rst excited states of the original TBF are approximately equally populated, so the cloning procedure is applied creating instead two TBF, one in the ground state and one in the excited state.At this point, the potential energy surfaces for the ground and excited states have opposite gradients.This leads to the acceleration of the hydrogen atom for the TBF associated with the ground state and, at the same time, slows it down for the exited state TBF.As a result, although both branches are leading to dissociation, the kinetic energies of the ejected atoms are signicantly different: the ground state branch contributes to the high energy peak of the distribution in Fig. 3, while the excited state branch contributes to the low energy peak.For the ground state branch, the remaining vibrational energy of the radical is low, so it remains in the ground state for the rest of the run and does not reach the intersection again.For the excited state branch, the energy taken away by the hydrogen atom is lower, leaving the pyrrolyl radical with sufficient energy to pass through numerous intersections with population transfer between the ground and both excited states.Naturally, quenching to the ground state will happen eventually for this branch but the time scale of this process is much longer than that for the dissociation, while the TKER spectrum is only affected by the radical dynamics until the H atom is lost.
In order to calculate the velocity map image with respect to the laser pulse polarization, we must average the velocity distribution of hydrogen atoms relative to the axes of the molecule, given by calculations, over all possible orientations of the molecule: Iðr; 4Þ $ dðr À jvjÞ∭da db dg cos 2 ðxða; b; gÞÞdð4 À fðq; a; b; gÞÞ; where a, b and g are Euler angles, q is the angle between the atom velocity vector v and the transition dipole of the molecule, x(a,b,g) is the angle between the transition dipole and light polarization vectors, and f(q,a,b,g) is the angle between the light polarization vector and atom velocity.Here we take into account that the probability of excitation is proportional to cos 2 (x).Integrating over Euler angles and replacing, as usual, the d-function for |v| with a narrow Gaussian function, we obtain Fig. 5 shows the simulated velocity map with respect to the laser pulse polarization assuming that the transition dipole is normal to the molecular plane.The simulations reproduce well the main feature of the velocity map image, which is the anisotropy of the intense high energy part.Our results are also consistent with experiment 14 in the low energy region showing an isotropic distribution, although admittedly the statistics of both experiment and simulation are poorer in the region of low energy.

V. Conclusion
We simulated the photodissociation dynamics of pyrrole excited to the lowest singlet excited state ( 1 A 1 / 1 A 2 ) using a new implementation of the AIMC approach, which now is modied to take into account the tunnelling of hydrogen atoms more accurately.AIMC is a fully quantum technique but its computational cost in our implementation is compatible with classical "on the y" molecular dynamics, which allows the accumulation of sufficient statistics to clarify the details of photo-induced processes in pyrrole.The treatment of tunnelling in our implementation provides a promising starting point for the further development of fully quantum methods for non-adiabatic dynamics and tunnelling with the ultimate goal of reaching well converged quantitative results.The current version of AIMC is already accurate enough to reproduce features of the experimentally observed TKER spectrum and velocity map images.

Fig. 1 A
Fig.1A sketch of the AIMC propagation scheme.The wave-function is represented as a superposition of Gaussian coherent states, which form a train moving along the trajectory.After passing the intersection, the train branches in the process of cloning.The figure shows a single train with cloning.In the most detailed AIMC calculation, a basis of several cloning trains interacting with each other is used.

Fig. 2
Fig. 2 Illustration of the algorithm used to treat tunnelling in our approach.(A) Identify turning point; (B) find a point with the same potential energy on the opposite side of the barrier; (C) run an additional trajectory through this point; (D) solve time-dependent Schr ödinger equation in the basis of a coherent state trains 10 moving along the trajectories on both sides of the barrier.

Fig. 3
Fig.3Total kinetic energy release (TKER) spectrum of hydrogen atoms after dissociation calculated with (solid) and without (dash) taking tunnelling into account.Both spectra are averaged over the same ensemble of initial configuration.The curves are smoothed by replacing delta-functions with Gaussian functions (s ¼ 200 cm À1 ).The inset shows the experimentally measured spectrum.14

Fig. 4 2 ! cos 2 ðqÞcos 2 ð4Þ þ 1
Fig. 4 An example of trajectory bifurcation on conical intersection.Electronic state populations (a), the kinetic energy of the H atom (b) and the N-H distance (c) as a function of time.Fast and slow branches are referred as (1) and (2) respectively.The black vertical line indicates the moment when cloning was applied.

Fig. 5
Fig. 5 Simulated velocity map image with respect to the laser pulse polarization assuming that the transition dipole moment is normal to the molecule plane.The experimental VMI 14 is shown in the inset.