A grid-based Ehrenfest model to study electron-nuclear processes

TL;DR

This paper introduces a grid-based Ehrenfest model for electron-nuclear processes that simplifies the computational complexity of simulating diatomic molecules by combining quantum and semi-classical dynamics.

Contribution

The work presents a novel application of the Ehrenfest approach to treat nuclear motion in electron-nuclear systems, enabling scalable simulations with multiple nuclear degrees of freedom.

Findings

Ehrenfest model shows good agreement with fully quantum results for H₂⁺.

The approach is computationally efficient and easily extendable to more nuclear degrees of freedom.

Demonstrates potential for modeling complex electron-nuclear interactions in molecules.

Abstract

The two-dimensional electron-nuclear Schr\"odinger equation using soft-core Coulomb potentials has been a cornerstone for modeling and predicting the behavior of one-active-electron diatomic molecules, particularly for processes where both bound and continuum states are important. The model, however, is computationally expensive to extend to more electron or nuclear coordinates. Here we propose to use the Ehrenfest approach to treat the nuclear motion, while the electronic motion is still solved by quantum propagation on a grid. In this work we present results for a one-dimensional treatment of H$_2^+$, where the quantum and semi-classical dynamics can be directly compared, showing remarkably good agreement for a variety of situations. The advantage of the Ehrenfest approach is that it can be easily extended to treat as many nuclear degrees of freedom as needed.

Figures

Equations

$$V_{sc}(z;R(t),E(t))=-\frac{1}{\sqrt{(z-R(t)/2{)}^2+{ϵ}^2}}-\frac{1}{\sqrt{(z+R(t)/2{)}^2+{ϵ}^2}}+\frac{1}{R(t)}-zE(t)$$

$$i\frac{\partial }{\partial t}\psi (z,t)=-\frac{1}{2}\frac{{\partial }^2}{\partial z^2}\psi (z,t)+V_{sc}(z;R(t),E(t))\psi (z,t)$$

$$\frac{d^2}{d{t}^2}R(t)=-\frac{1}{M}⟨\psi (z,t)\left|\frac{\partial }{\partial z}{V}_{sc}(z;R(t),E(t))\right|\psi (z,t)⟩$$

Full text

I Introduction

Full control over chemical reactions and molecular properties requires the manipulation of electronic and nuclear degrees of freedom (dofs). When the system is in the ground electronic state, the energy and time-scales associated to both dofs is often very different, allowing for a separate study of their dynamics, in the spirit of the Born-Oppenheimer approximation. This ceases to be the case when the dynamics occurs in excited electronic states, where non-adiabatic processes (i.e. by definition beyond Born-Oppenheimer) are the usual suspects behind the fate of most molecular events leading to energy deactivation or to a specific channel in a chemical reaction Yarkony (1996); Worth and Cederbaum (2004); Domcke et al. (2011).Because of computational constraints, most non-adiabatic processes have been analyzed expanding the wave function on a small subset of electronic states Persico and Granucci (2014); Beck et al. (2000); Meyer et al. (2009); Tully (1998); Yonehara et al. (2012); Richter et al. (2011); Tully (2012); Mai et al. (2018); Ben-Nun et al. (2000); Shalashilin (2009). However, if the energy in the system is large, so is the density of states, and different types of models must be elaborated to fully account for electron and nuclear correlations and their control, particularly in the presence of strong fields.

Recently, there has been a surge of studies that aim to analyze and control electron-nuclear processesKim et al. (2012); Braun et al. (2014); Calegari et al. (2014); Sola et al. (2015); Wolter et al. (2016); Wang et al. (2017); Attar et al. (2017); Nisoli et al. (2017); Amini et al. (2019) to create novel transient molecular properties in the presence of strong fields, for instance huge electronic dipoles Chang et al. (2013, 2015a, 2015b), to manipulate non-adiabatic transitions Hofmann and de Vivie-Riedle (2001); Gonzalez-Vazquez et al. (2006, 2009); Falge et al. (2012) or to unravel the electron-nuclear dynamical features in conical intersections Hader et al. (2017); Arnold et al. (2018); Csehi et al. (2017) particularly at so-called light-induced conical intersections Moiseyev et al. (2008); Halász et al. (2012); Demekhin and Cederbaum (2013); Halász et al. (2015), Our goal is to design a simple model that can provide qualitative predictions of quantum control in one-active electron systems, treating strong field laser couplings, non-adiabatic couplings and ionization on equal footing, and that can be extended to polyatomic molecules. One possible avenue is to use the Multiconfiguration Time-Dependent Hartree (MCTDH) methodBeck et al. (2000); Meyer et al. (2009), the Ab-Initio Multiple Spawning (AIMS) methodBen-Nun et al. (2000); Shalashilin (2009) or other schemes that incorporate quantum features to the nuclear motionPrezhdo (2006); Zamstein and Tannor (2012a, b); Subotnik (2010); Curchod and Tavernelli (2013).

Since we are interested in few active electrons (one, in this work) under a strong field that can lead to a large deformation of the molecule and its charge distribution over large distances, the typical orbital basis that are used to describe the electron density in Quantum Chemistry are not well developed for this purpose. Hence, most methodologies, which are based on expanding few electronic states on a basis, will not perform adequately. In this work we will follow a different methodology. We solve the time-dependent Schrödinger equation for the electron on a gridKosloff (1994), incorporating the vibrational motions in a semi-classical manner, via a mean-field Ehrenfest approachTully (1998, 2012)

In order to compare the results of the Ehrenfest approach to the full quantum calculation, in this paper we will focus on the simplest molecule under strong fields, treated in low-dimensionality by assuming that all particles are aligned and interact via soft-core Coulomb potentials. Although implying strong approximations, the quantum one-electron plus one-nuclear dimension Hamiltonian has been shown to give results in qualitative agreement with experiments, and as a model it has provided invaluable guidance to find and to understand new processes of molecules in strong fieldsSola et al. (2018).

We will consider three different situations to test the performance of the Ehrenfest approach: Laser-free dynamics of a coherent superposition of electronic states, strong-field laser-driven dynamics in the ground state, and the dynamics in the excited dissociative state under a strong static field.

The first case is interesting because it is well known that the dynamics in a superposition state cannot be well reproduced by an Ehrenfest approach, since each quantum wave packet in the superposition will feel a different potential than the mean-field potential governing the evolution of the nuclear trajectory. Because of decoherence, however, the impact of this difference over the averaged results may be much smaller than anticipated. The opposite situation occurs when the dynamics is adiabatic (i.e. slow) under a strong field. Then the nuclei experience an average, so-called light-induced potentialYuan and George (1978); Bandrauk and Sink (1981); Chang et al. (2016) (LIP), which can be perfectly reproduced by the Ehrenfest methodBajo et al. (2012). In the second case we will use optical non-resonant fields where in principle non-adiabatic effects might be important. The difference between the quantum and semiclassical results will quantify to a certain extent the importance of non-adiabatic effects. In the third case we use strong static fields, where the dynamics is expected to be fully adiabatic (that is, to occur in a single LIP). Here, we expect that the difference between the quantum and semi-classical results measure to a certain degree the impact of ionization and dissociation, which cannot be well reproduced in a mean-field theory.

II Model

A simple yet powerful enough model often used to integrate the dynamics of a system composed of two nuclei and one electron under a strong field is provided by the Hamiltonian of the aligned particles moving under the effect of the soft-core Coulomb potentialJavanainen et al. (1988); Su and Eberly (1991); Kulander et al. (1996). Within the range of validity of the model one can treat quantum mechanically both the electron and nuclear motion and study strong field effects beyond the Born-Oppenheimer approximation. In this work we will treat the nuclear motion classically under the Ehrenfest approximation testing how the results compared with those obtained in a fully quantum $2$ dimensional model, under different conditions, so that one can caliber the quality of the approximation before working in models with larger degrees of freedom where the quantum results can not be directly obtained. In our semiclassical approximation, both the field $E(t)$ and the bond distance $R(t)$ are treated as classical variables that enter as parameters in the potential. For H${}_{2}^{+}$,

[TABLE]

The electron motion in the $z$ coordinate obeys the time-dependent Schrödinger equation (in a.u.)

[TABLE]

whereas the internuclear motion follows Newton equations with the mean-field electronic potential in the Hellmann-Feynman approximation Tully (1998).

[TABLE]

where $M$, the reduced mass of the molecular frame, is approximately the mass of the proton.

The initial conditions for the dynamical equations can be obtained under different procedures. In the fully quantum TDSE for both electronic and nuclear dof the initial wave function, $\Psi(z,R,0)$, will typically be a product state of the energy eigenstate of the electronic Hamiltonian $\psi_{j}(z;R)$ (or a superposition of different electronic states) times a Gaussian nuclear wave function, $\chi(R)$, centered at different internuclear distances (for instance, the equilibrium bond distance). We will reproduce laser-free or laser-driven dynamics following the examples of Chang et al.Chang et al. (2015b)

In the semiclassical approach the initial internuclear position and momentum are obtained by a Monte Carlo sampling from the Wigner distribution of $\chi(R)$, $\rho_{W}(R,P)$. For each trajectory $k$, one obtains a different set of phase space coordinates $(R_{k}(0),P_{k}(0))$. However, we sometimes change the distribution to test the sensitivity of the scheme to the initial nuclear coordinates. In addition, since the $\psi_{j}(z;R)$ are obtained discretized on a grid, by the Fourier grid Hamiltonian methodMarston and Balint-Kurti (1989), while the phase space points are continuous, we choose the electronic wave function $\psi_{j}(z;\tilde{R}_{k})$ at the grid point $\tilde{R}_{k}$ closer to $R_{k}$. In order to compare the results obtained by the semiclassical method with those obtained solving the TDSE for the coupled electron-nuclear motion, one needs to average the results of different trajectories. Eq.(2) is solved using the split-operator methodFeit et al. (1982) while Eq.(3) is solved using the forth-order Runge-Kutta method.

III Results

III.1 Field-free Dynamics

A stern test for a mean-field theory to pass is to compare the results of its dynamics to the quantum case when the initial wave function is a superposition of different electronic states, in the absence of a field that couples (i.e. mixes) the potentials. Although no single Ehrenfest trajectory can reproduce the average observables in this situation, the ensemble average of the trajectories may give reasonable results due to the decoherence that the vibrational motion induces on the electronic dynamics. Focusing only on the electronic dof, the decoherence provokes dumping in the oscillation of the average electron position $\langle z(t)\rangle$, because the interfering pathways that create the oscillations between the electronic states become dephased due to the different vibrational periods in each stateChang et al. (2015a).

We first consider the dynamics starting in a symmetric superposition state $\psi(z;R(0))=\frac{1}{\sqrt{2}}\left(\psi_{1}(z;R(0)))+\psi_{2}(z;R(0))\right)$, where $\psi_{j}(z;R)$ ($j=1,2$) are the ground/first excited electronic states of H${}_{2}^{+}$ in the soft-core Coulomb potential. In the full quantum results we use for the initial nuclear wave function a Gaussian wave packet $\chi(R)$ with $0.5$ $a_{0}$ width, centered at $R_{0}=2$ $a_{0}$, that approximately represents the ground state of the parent $H_{2}$ molecule before ionization.

In  1 we compare the results of the quantum simulation with those of the Ehrenfest approach, after averaging $N=10000$ trajectories. To test the sensibility of the Ehrenfest approach to the characteristics of the nuclear ensemble of trajectories, the initial nuclear coordinates and momentum are taken from random sampling of three different Wigner distributions. One corresponds to the nuclear wavefunction of the quantum calculation ($\tilde{\rho}(R)=\int dP\rho_{W}(R,P)=|\chi(R)|^{2}$, $\tilde{\rho}(P)=\int dR\rho_{W}(R,P)=|\chi(P)|^{2}$). The others are squeezed distributions in the coordinate or momentum representations, as shown in 2, as if we have doubled/halved the width of the wave packet.

The results of the simulation show good qualitative agreement between the quantum and the Ehrenfest average observables. Quantum electronic coherence survives several oscillations until it is completely suppressed by decoherence. In our model this occurs when the nuclear wave packets on the two electronic states cease to overlap each other. Obviously, because we are using a one-dimensional model and the wave packets remain bounded, revivals would be observed at a later time.

The Ehrenfest model matches the results during the first oscillations. However, the classical averaging induces a faster decay than the decoherence rate at the beginning, although after a few femtoseconds the killing of coherence is also more complete in the quantum case than in the semiclassical approach. Regarding the nuclear motion, the average bond distance is also very well described at short times. Although we are following the dynamics for the first three femtoseconds, the duration is already significant since the bond stretches one Bohr. Of course, in the quantum case one can separately follow the average distance in each electronic state, which are quite different than the mean value.

Comparing the results obtained for different Wigner distributions, $\rho_{W}(R,P)$, we find that the induced dephasing is stronger both when we start with a wider and a narrower distribution in the internuclear distance $\tilde{\rho}(R)$. The former case is expected, since we sample a wider range of initial positions, each originating a different trajectory, which induce the dephasing after averaging. To understand the latter case one should view the distribution in the momentum representation, $\tilde{\rho}(P)$, which by the Fourier transform is wider the narrower $\tilde{\rho}(R)$ is. For each system there will be an “optimal” initial Wigner distribution that minimizes the dephasing rate. Just by coincidence, in our example this distribution is approximately the one obtained from the chosen initial nuclear Gaussian wave function.

In the semiclassical approach the results will always be sensitive to the number of trajectories that are included in the ensemble averageNelson et al. (2012). This is particularly important when the mean-field departs clearly from the different electronic potentials that act on the wave packets. Obviously, a single trajectory under this mean field will evolve in a fully coherent way: the period of the oscillations in $\langle z(t)\rangle$ will depend on the energy difference between the potentials $V_{2}(R(t))-V_{1}(R(t))$, while the amplitude of the oscillation will depend on the distance between the turning points in the mean field potential. In  3 we give an example of a typical trajectory ($N=1$) and then show the results of the average $\langle z(t)\rangle$ as a function of the number of trajectories $N$ included in the ensemble. Even a small number of trajectories can suffice to mimic the results of the full distribution during the first oscillations. However, the dephasing is not complete and one observes immediate revivals. To obtain fully converged results for larger times one needs to sum over a very large number of trajectories, which in this case is around $N\sim 5000$.

III.2 Field-driven dynamics

We now consider what happens when the dynamics is driven by a strong nonresonant laser pulse in the ground state. The quantum results show that the field $E(t)$ moves the electron (the electron position is actually anti-correlated with the field) generating a small dipole synchronized with the external frequency proportional to the bond distance. But the proportionality constant is small unless the fields are very intense. On the other hand, for too intense fields the molecule ionizesChang et al. (2015b).

In 4 we show the average electron position (the dipole) and bond distance using a $3$-cycle $1.8$ $\mu$m pulse of $100$ TW/cm2 peak intensity. The Ehrenfest results match closely the quantum dynamics for the electron, even if the deviation of the average bond distance is noticeable. For these pulse parameters, the Keldysh adiabatic parameterSeideman et al. (1995) $\gamma=\omega\sqrt{I_{p}}/E_{0}$ ($I_{p}$ is the ionization potential from the initial state, $E_{0}$ the peak field amplitude and $\omega$ the field frequency) is $\gamma\sim 0.5$, which indicates that there is some tunneling ionization. It is interesting to observe that Ehrenfest reproduces all the important features of the process, including the exact form of the wiggles in $\langle z(t)\rangle$. Ionization in the absence of dissociation does not deteriorate the quality of the results. Small deviations in the amplitudes occur. They are expected whenever nonadiabatic effects may be important, as induced by the shortness of the pulse.

When the pulses are too short, it might be expected that the effects of the nuclear motion on the electron dynamics are small, and hence any failure in its description will not modify the electronic observables. In 5 we follow the dynamics driven by a $30$-cycle $800$ nm pulse of $25$ TW/cm2 peak intensity. In this case the Keldysh parameter is $\gamma\sim 2$ and given the pulse amplitude, multi-photon ionization is not expected (it is smaller than 0.1%).

Since the laser period, $2\pi/\omega$, is much shorter than the time-scale of the molecular motion, the latter will be mostly decoupled from the electron motion. As observed in 5 the vibrational motion in the ground state induces a modulation on the dipole, more noticeable during the first oscillations which occur before the nuclear wave packet spreads. This effect is perfectly described by the Ehrenfest model, despite the fact that the average bond distance stretches in time, an indication that the mean-field potential shows bond-softening, which is not playing an important role in the quantum dynamics. In the quantum case, nonadiabatic effects lead to some dissociation (smaller than $10$%), but the dissociating packet, after some threshold value, is not included in the average bond distance. On the other hand, it always contributes to the average potential, leading to the excess in bond softening simulated by the Ehrenfest approach. Still, the effect over the electron observables is small.

The dynamics in the excited state in the presence of a constant field is a fine example of fully correlated electron-nuclear motion, which can easily be explained as the adiabatic motion of a wave packet in the excited LIP, that shows bond hardening. In perfect conditions, for large bond distances, the dipole is (in atomic units) $\langle z(t)\rangle=\langle R(t)\rangle/2$Chang et al. (2015a). The dressed electronic state implies the electron moving with the left proton for a positive field, and the latter moving back and forth as it reaches the classical turning points of the LIP. This creates the fully correlated oscillations in both observables. As 6 shows, both the dipole and the bond oscillations are dumped. Behind this dumping is the nuclear wave packet spreading due to the anharmonicity of the potential, as can be easily verified looking at revivals at later times. As already shown in Sec.III.1, the effects of the spreading of the packet can be well accounted for by averaging trajectories. In this example convergent results can be obtained with $N\sim 1000$ or less.

Small nonadiabatic effects occur at early times if the initial electronic wave function is the excited molecular electronic state, and not the dressed state. The electron then departs from $\langle z(0)\rangle=0$ and not from $\langle R(0)\rangle/2$ and some population can leak to the lowest-energy LIP where the packet dissociates due to bond softening. In addition, this also explains why $\langle z(t)\rangle$ is smaller than $\langle R(t)\rangle/2$. Unlike in weaker fields, under strong fields like $E_{0}=0.04$ a.u. the LIP is strongly bounded and the electron simply does not have enough time to overcome the initial separation and reach half the bond distance before the nuclei get to the turning point of the potential and both the bond and the dipole shrink.

The effect of the dissociation is always reflected on a slightly weaker mean-field potential which leads to a small relaxation in the Ehrenfest values of $\langle R(t)\rangle$. Due to the electron-nuclear correlation this is reflected also in $\langle z(t)\rangle$ in this case. Interestingly, this effect can be partly compensated for by changing the initial Wigner distribution. Sampling a wider initial distribution of internuclear bond distances the average values of the Ehrenfest method get closer to the quantum values, although at the expense of having (as expected from the results of Sec.III.1) over-dumped oscillations due to the larger decoherence.

IV Conclusions

In this work we have developed and test an algorithm that mixes grid propagation for the electron (active) coordinates with an Ehrenfest approach to the nuclear coordinates. We have simulated the dynamics of a collinear model for H${}_{2}^{+}$ with or without an external strong field, starting in the electronic ground state, excited state, or a superposition of states and compared the results with those obtained by solving the full electronic plus nuclear TDSE of the system. The results show quantitatively good agreement between the methods, at least for short enough times when the dissociation is not prevalent. For sufficient sampling of trajectories, in general we observe that the semiclassical results tend to over-characterize the initial dephasing and the bond stretching. The first problem can be at least partially handled by modifying the initial Wigner distribution or biasing the sampling. The mean field results cannot reproduce the dynamics when the wave function branches over electronic states with very different forces, so that one cannot expect good agreement for the long-time behavior of the nuclear dynamics in chemical reactions. However, the results prove that the Ehrenfest approach can provide a feasible alternative to test strong-field laser-control processes over transition states of molecules with several bonds, characterizing well ionization processes and severe bond geometrical distortions, which will be the subject of further work.

Acknowledgments

This work is supported by the Korean government through the Basic Science Research program (2017R1A2B1010215) and the National Creative Research Initiative Grant (NRF-2014R1A3A2030423), by the Spanish government through the MINECO Project No. CTQ2015-65033-P and by the Comunidad de Madrid through project Y2018/NMT-5028. We thank the Army Research Laboratory for the hospitality during which part of this work was created.