Quantum Nonadiabatic Cloning of Entangled ... - ACS Publications

Letter pubs.acs.org/JPCL

Quantum Nonadiabatic Cloning of Entangled Coherent States Artur F. Izmaylov* and Loïc Joubert-Doriol Department of Physical and Environmental Sciences, University of Toronto Scarborough, Toronto, Ontario M1C 1A4, Canada Chemical Physics Theory Group, Department of Chemistry, University of Toronto, Toronto, Ontario M5S 3H6, Canada S Supporting Information *

ABSTRACT: We propose a systematic approach to the basis set extension for nonadiabatic dynamics of entangled combination of nuclear coherent states (CSs) evolving according to the time-dependent variational principle (TDVP). The TDVP provides a rigorous framework for fully quantum nonadiabatic dynamics of closed systems; however, the quality of results strongly depends on available basis functions. Starting with a single nuclear CS replicated vertically on all electronic states, our approach clones this function when replicas of the CS on different electronic states experience increasingly different forces. Created clones move away from each other (decohere), extending the basis set. To determine a moment for cloning, we introduce generalized forces based on derivatives that maximally contribute to a variation of the total quantum action and thus account for entanglement of all basis functions.

T

he time-dependent variational principle (TDVP)1−3 provides variationally optimal equations of motion (EOM) for the system wave function specified by a certain ansatz. The TDVP allows one to model the quantum nuclear wave function in a computationally efficient way for both adiabatic and nonadiabatic nuclear dynamics in molecules. There are two main popular forms of the nuclear wave function: (1) originating from the multiconfiguration timedependent Hartree (MCTDH) method4−6 and its multilayer generalizations7,8 and (2) based on frozen-width Gaussians,9−17 which are moving either classically9−11,16,17 or quantum mechanically.12−15 The latter ansatz, due to locality of involved basis functions, is very well suited to be used in conjunction with the on-the-fly solution of the electronic structure problem.10,11,17 The main practical difficulty for any dynamical method based on the TDVP is basis set limitation. If we consider nonadiabatic dynamics using a linear combination of frozen-width Gaussians NG

|Ψ(t )⟩ =

ically.9,10 If a Gaussian arrives at a region of strong coupling between electronic states and there is no Gaussian on the other state to interact with it, the spawning algorithm creates the counterpart needed for population exchange (Figure 1S). This

Ns

∑ ∑ CI(s)(t )|GI(s)(t )⟩|ϕs⟩ I=1 s=1

Figure 1. Illustration of the spawning (S) and cloning (C) procedures. (1)

consideration may seem ad hoc and does not account for the fact that each Gaussian basis function is a part of the total nuclear wave function. However, the spawning approach can be also rigorously introduced using time-dependent perturbation theory15 that takes the total wave function into account and provides a route for dynamical basis set extension. This perturbative spawning has been extended to the fully quantum propagation schemes such as the variational multiconfiguration Gaussian (vMCG) method where the dynamic of Gaussians has a highly entangled quantum character.15

where CI(s) are time-dependent coefficients (amplitudes), indices s and I enumerate electronic states |ϕs⟩ and Gaussians |G(s) I ⟩, respectively, the population transfer between electronic states can only take place when Gaussians located on different electronic states have significant overlap in nuclear degrees of (s ) freedom (DOF), ⟨G(s) I |GJ ′ ⟩ ≫ 0. However, considering the localized nature of Gaussians and that different electronic surfaces provide different forces in the same area of nuclear geometry, these overlaps generally quickly decay along the dynamics. This decoherence process artificially reduces the electronic population transfer. To address this issue, the spawning technique was introduced for a linear combination of frozen-width Gaussians whose parameters evolved classically while the amplitudes were propagated quantum mechan© XXXX American Chemical Society

Received: March 10, 2017 Accepted: April 4, 2017 Published: April 4, 2017 1793

DOI: 10.1021/acs.jpclett.7b00596 J. Phys. Chem. Lett. 2017, 8, 1793−1797

Letter

The Journal of Physical Chemistry Letters ⎡ ω ⎛ ωj ⎞1/4 j ⎜ ⎟ ⟨x|GI (t )⟩ = ∏ exp⎢ − [xj − qjI (t )]2 ⎝ ⎠ ⎢ π 2 j=1 ⎣

Alternatively, one can approach the problem of population transfer in TDVP-based nonadiabatic dynamics by introducing a nuclear basis with the condition |GI(s)⟩ = |GI(s ′)⟩ = |GI ⟩

s ≠ s′

Ndim

(2)

+ ipjI [xj − qjI (t )] +

This condition will ensure the maximum overlap between (s ) Gaussians on different electronic states ⟨G(s) I |GI ′ ⟩ = 1. The wave function becomes NG

|Ψ(t )⟩ =

I=1 s=1

(3)

or equivalently |Ψ(t )⟩ =

(5)

Here, x are nuclear coordinates, dim(x) = Ndim, and qI(t) = {qjI(t)}j=1,Ndim and pI(t) = {pjI(t)}j=1,Ndim are time-dependent positions and momenta. EOM for all parameters of the wave function in eq 4 can be obtained by finding an extremum of the action S = ∫ ⟨Ψ(t)|Ĥ − i∂t|Ψ(t)⟩ dt, which is equivalent to solving24

Ns

∑ ∑ CI(s)(t )|GI (t )⟩|ϕs⟩

⎤ i pjI qjI ⎥ ⎥⎦ 2

Re⟨δ Ψ(t )|Ĥ − i∂t|Ψ(t )⟩ = 0

NG

Ns

I=1

s=1

∑ |GI (t )⟩[∑ CI(s)(t )|ϕs⟩]

(6)

Here, Ĥ is the system Hamiltonian. For the parametrization of eq 4, it is easy to show that such a version of the TDVP is equivalent to those of Dirac−Frenkel2,3 and McLachlan.25 Introducing the SS variations

(4)

Therefore, this basis, also known as the single-set (SS) basis, can be thought of as consisting of either stacks of identical Gaussians replicated for all electronic states (eq 3) or Gaussians with individual time-dependent electronic functions (eq 4). Although the SS Gaussians always can exchange the population between electronic states, a new problem arises: replicas cannot take individual paths or decohere; instead, each SS Gaussian moves on an average Ehrenfest-like surface. To introduce more freedom, the cloning technique was suggested:16 the algorithm monitors the difference in forces that replicas within a SS stack experience on different electronic states. When the force difference becomes large, the cloning scheme splits the stack of Gaussians in two clones and adds empty replicas of Gaussians for parts of the stack that went to another clone (Figure 1C). As in the case of spawning, cloning has been introduced for frozen-width Gaussians that are evolving classically on Ehrenfest-like surfaces.16,17 The algorithm treats every stack of Gaussians independently, and thus, evaluation of forces is straightforward. However, such cloning has the same drawback as spawning: it does not treat each Gaussian as a part of the total wave function. In order to put the cloning idea on a rigorous quantum basis as well as to extend it to fully quantum treatment of nuclear dynamics, one should consider the case of quantum entangled Gaussians with corresponding quantum forces that originate from the total nuclear wave function. This is exactly the aim of the current Letter, where we propose a cloning algorithm for the fully quantum nonadiabatic dynamics in the basis of SS frozen-width Gaussians. For the sake of simplicity, our technique will be illustrated on a set of two-state low-dimensional diabatic models where the exact quantum results can be easily obtained. However, nothing prevents the use of the approach in the adiabatic representation with the on-the-fly generation of potential electronic surfaces. To treat challenging geometric phase effects arising in the conical intersection case,18−22 one can use a recently introduced scheme evaluating adiabatic electronic functions only at Gaussian centers.16,17,23 Equations of Motion for the SS Representation. We start with the total nonstationary wave function given by eq 4 where Gaussians are taken in the coherent state (CS) form, which improves numerical stability of the vMCG scheme15

∑

|δ Ψ⟩ =

I ,s

∂Ψ ∂CI(s)

δCI(s) +

∑ j,I

∂Ψ ∂qjI

δqjI +

∂Ψ ∂pjI

δpjI (7)

and accounting for independence and arbitrariness of individual variations δC(s) I , δqjI, and δpjI leads to the EOM for all parameters13,26 (s) iCİ =

∑ [S−1(Hss ′ − iτδss′)]IJ C(J s ′) (8)

J ,s′

iξjİ = [B−1Y]jI

(9)

where ξjI = ωjqjI + ipjI are convenient variables encoding both position and momentum components of CSs. Matrices involved in eqs 8 and 9 are τIJ = ⟨GI |∂tGJ ⟩

SIJ = ⟨GI |GJ ⟩

(10)

Hss ′ , IJ = ⟨GI |⟨ϕs|Ĥ |ϕs ′⟩|GJ ⟩

BIk , Jn =

(11)

∑ CI(s) *C(J s)(S(kn) − S(k0)S−1S(0n))IJ

(12)

s

YIk =

∑ CI(s) *C(J s ′)(H(ssk′0) − S(k0)S−1Hss ′)IJ (13)

ss ′ , J

Hss(k′0) , IJ =

∂GI ⟨ϕ|Ĥ |ϕs ′⟩ GJ ∂ξkI s

SIJ(kn) =

∂GI ∂GJ ∂ξkI ∂ξnJ (14)

SIJ(k 0) =

∂GI GJ ∂ξkI

SIJ(0n) =

GI

∂GJ ∂ξnJ

(15)

Time derivatives of CSs needed in the τ matrix are derived using the chain rule |∂tGK ⟩ =

∂GK ∂qK

q̇K (t ) +

∂GK ∂pK

pK̇ (t )

(16)

Solving eqs 8 and 9 constitutes the vMCG approach within the SS basis set. 1794


Letter

The Journal of Physical Chemistry Letters

Figure 2. Electronic population (P) as a function of time for Ĥ SB and Ĥ CI models for different couplings: (SBw) Ĥ SB weak coupling, V = 0.1ω2; (SBs) Ĥ SB strong coupling, V = 0.5ω2; (CIw) Ĥ CI weak coupling, c2 = 0.1ω2; (CIs) Ĥ CI strong coupling, c2 = 0.5ω2. In QC-vMCG (dashed curves), various ε’s produced different numbers of SS pairs at the end of the propagation; they are given by the NG values. Positions of the initial CS are qc = −3 (1D) and qc = (−3,−1) (2D).

where ε is an accuracy threshold. Interestingly, because we use |ΨI⟩ from the SS simulation, the sum of derivatives over electronic states is always zero

Cloning SS Pairs. If we consider a variation of the total wave function that changes positions and momenta of replicas for an Ith CS on different electronic states independently |δI Ψ⟩ =

∑ s,j

∂Ψ ∂ξjI(s)

δξjI(s)

∑ Re (17)

s

the condition of eq 6 will not be satisfied. Formally, to consider such variation, we need to evaluate it on a wave function obtained from |Ψ⟩ by allowing the Ith CS’s replicas to be different for different electronic states |ΨI (t )⟩ =

J≠I

(18)

To determine when and which of the SS pairs split for cloning, it is instructive to consider the variation Re⟨δI Ψ|I Ĥ − i∂t|Ψ⟩ I ⎡ ∂ΨI = Re⎢∑ δξjI(s) Ĥ − i∂t ΨI ⎢⎣ j , s ∂ξjI(s) ≠0

⎤ ⎥ ⎥⎦ (19)

This quantity contains arbitrary variations δξ(s) jI , which can be removed if one is interested in the effect of splitting of the Ith SS pair on the action. Thus, our criterion for splitting the Ith SS pair is

∑ j,s

Re

∂ΨI ∂ξjI(s)

Ĥ − i∂t ΨI

∂ξjI(s)

Ĥ − i∂t ΨI

=0 (21)

which is consistent with a zero state average value. Therefore, the sum in eq 20 corresponds to the norm of the deviation of generalized quantum state-specific forces acting on an individual CS from the state-averaged counterpart. Equation 20 provides a rigorous estimate of a decoherence strain on variation of the quantum action. Further intuitive simplifications can be detrimental for the cloning process because they may remove important components contributing to the variationally optimal character of eq 20. For instance, a (s ̇ ̇ ) criterion based only on ∑j |ξ(s) jI − ξjI ′ | > ε would ignore a difference between populations of different replicas of the Ith CS. Note that more than one SS pair can be split using the criterion of eq 20 at a time, but for the sake of simplicity of further discussion, we assume that only one pair has been split. Once the decision on splitting is made, to avoid linear dependency between clones, we propagate the split pair treated as independent CSs along with NG − 1 unsplit SS pairs. The EOM for such hybrid evolution are obtained using the TDVP applied to the parametrization ΨI in eq 18 and detailed in the SI. CSs of the split pair move on different potential energy surfaces and necessarily decohere so that the overlap integral (2) ⟨G(1) I |GI ⟩ will decrease, allowing one to create two new SS pairs without introducing linear dependency

∑ [CI(s)(t )|GI(s)(t )⟩ + ∑ C(J s)(t )|GJ (t )⟩]|ϕs⟩ s

∂ΨI

>ε (20) 1795


Letter

The Journal of Physical Chemistry Letters ⎛C (1)|G (1)⟩⎞ ⎛ 0 ⎞ ⎛ (1)⎞ I (2) ⎜ I ⎟ → ⎜CI ⎟|G (1)⟩, ⎜ I ⎜ (2)⎟⎟|GI ⟩ ⎟ ⎜ ⎜ (2) (2) ⎟ C ⎠ ⎝ ⎠ ⎝ 0 I ⎝CI |GI ⟩⎠

but also derivatives of unsplit CSs are reduced. The latter is the effect that comes from quantum entanglement of CSs in the total nuclear wave function. In conclusion, we have introduced a novel general algorithm to extend the basis set when needed in quantum dynamical simulations based on the magnitude of the quantum forces that have maximal effect on the action variation. These derivatives become large when pairs of nuclear CSs located on different potential energy surfaces experience very different forces. Similar developments were done for classically moving CSs, where ad hoc criteria of pair separation were introduced based on force differences. We rigorously extended these intuitive techniques to full quantum dynamics of CSs, with account for entanglement between different CSs in the total wave function. Our approach can be easily extended to more than two electronic states, the adiabatic representation, and on-the-fly generation of potential energy surfaces. Another useful extension can be a formulation of a spawning technique, which will use similar derivatives to determine a spawning event variationally. The work on this approach is underway and will be reported elsewhere.

(22)

where vectors are written in the basis of electronic functions {ϕs}s=1,2. Once the split CSs are cloned into two new SS pairs, the regular EOM (eqs 8 and 9) for the NG + 1 SS pairs are employed. We will refer to this algorithm as the quantum cloning vMCG (QC-vMCG) approach. Numerical Examples. We illustrate the performance of QCvMCG in modeling nuclear dynamics of one- and twodimensional two-state diabatic models Ĥ =

Ndim

∑ j=1

⎛[p ̂ 2 + ω 2x 2]/2 ⎞ cjxj j j ⎜ j ⎟ ⎛0 V⎞ ⎟ ⎜⎜ ⎟+⎜ 2 2 2 cjxj [pĵ + ωj (xj − dj) ]/2 ⎟⎠ ⎝ V Δ ⎠ ⎝

(23)

where xj and p̂j are nuclear coordinates and associated momenta, and V, Δ, dj, cj, and ωj are constants. In the onedimensional (1D) model (Ndim = 1), which is also known as spin-boson, Ĥ SB = Ĥ , where cj = 0. In the two-dimensional (2D) model (Ndim = 2), Ĥ CI = Ĥ , where V = d2 = c1 = 0; this setup gives rise to the conical intersection of potential energies if transformed to the adiabatic representation. Other parameters in both cases are ω1 = 0.89, ω2 = 0.9, d1 = 5, and Δ = −ω1. The last condition ensures resonance between vibrational levels of diabats coupled with linear potential coupling. Such resonances are unavoidable in large dimensional problems with conical intersections but can be missing in 2D models. We consider systems with strong and weak couplings, which are characterized by V and c2 for Ĥ SB and Ĥ CI, respectively. Weak couplings simulate diabatically trapped systems,27 while strong couplings bring systems closer to the adiabatic limit. However, strong nonadiabatic couplings are present in both cases. Also, large reorganization energy is maintained in all systems (d1 = 5) to make them challenging for the SS basis. Note that the dj = cj = 0 case can be solved with a single SS pair because diabatic states have identical nuclear dependence in this limit. We simulate nuclear dynamics starting with an initial wave function constituting a single SS pair |Ψ(t = 0)⟩ = |G1⟩[1 · |ϕ1⟩ + 0 · |ϕ2⟩] with zero initial momentum and centered at point qc. For each Hamiltonian, we simulate time-dependent wave functions and monitor the population of the first electronic state, P(t) = Trn[|⟨ϕ1|Ψ(t)⟩|2], where Trn is the trace over the nuclear coordinates (Figure 2). In all QC-vMCG calculations, lowering ε allowed us to converge to the exact dynamics generated by the split operator approach.28 It may seem that lower couplings require lower thresholds, but it partly comes from the scale of the plots. Stronger couplings make initially empty replicas of CSs be populated faster and generate a force difference for faster decoherence. Compared to 1D, in 2D, there are more ways for CSs to avoid each other and to lower the overlap between different pairs. Mutual help of CSs is weaker in 2D, and thus, more CSs are required in 2D for convergence. Besides cases in Figure 2, decoherence forces of eq 20 for extreme limits of the spin-boson model have been considered: (1) Ĥ SB with d1 = 0 produces zero derivatives in eq 20 because diabatic surfaces are parallel; (2) Ĥ SB with V → 0 also produces zero derivatives in eq 20 because the population transfer is negligible, and the initial CSs evolves on a single harmonic oscillator. Also, in a general case, it was confirmed that upon splitting not only do derivatives in eq 20 of the split pair vanish

■

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jpclett.7b00596. Explicit form of equations of motion for the case when two coherent states are allowed to move on their electronic states and the rest of coherent states form the single-set basis and details of the split operator numerical integration scheme (PDF)

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. ORCID

Artur F. Izmaylov: 0000-0001-8035-6020 Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS A.F.I. thanks I. G. Ryabinkin for critical reading of the manuscript and acknowledges funding from a Sloan Research Fellowship and the Natural Sciences and Engineering Research Council of Canada (NSERC) through the Discovery Grants Program.

■

REFERENCES

(1) Kramer, P.; Saraceno, M. Geometry of the Time-Dependent Variational Principle in Quantum Mechanics; Springer: New York, 1981. (2) Dirac, P. A. M. The Principles of Quantum Mechanics, 4th ed.; Clarendon Press: Oxford, U.K., 1958. (3) Frenkel, J. Wave Mechanics; Clarendon Press: Oxford, U.K., 1934. (4) Meyer, H.-D.; Manthe, U.; Cederbaum, L. S. The MultiConfigurational Time-Dependent Hartree Approach. Chem. Phys. Lett. 1990, 165, 73−78. (5) Wang, H.; Thoss, M. Multilayer Formulation of the Multiconfiguration Time-Dependent Hartree Theory. J. Chem. Phys. 2003, 119, 1289−1299. (6) Worth, G. A., J, A.; Beck, M. H.; Meyer, H.-D The MCTDH Package, development version 9.0; University of Heidelberg: Heidelberg, Germany, 2009.

1796


Letter

The Journal of Physical Chemistry Letters

Electronic Transitions through Conical Intersections. J. Chem. Phys. 2011, 135, 234106. (28) Tannor, D. J. Introduction to Quantum Mechanics: A TimeDependent Perspective; University Science Books: Sausalito, CA, 2007; p 214.

(7) Wang, H.; Thoss, M. Multilayer Formulation of the Multiconfiguration Time-Dependent Hartree Theory. J. Chem. Phys. 2003, 119, 1289. (8) Manthe, U. A Multilayer Multiconfigurational Time-Dependent Hartree Approach for Quantum Dynamics on General Potential Energy Surfaces. J. Chem. Phys. 2008, 128, 164116. (9) Yang, S.; Coe, J. D.; Kaduk, B.; Martínez, T. J. An “Optimal” Spawning Algorithm for Adaptive Basis Set Expansion in Nonadiabatic Dynamics. J. Chem. Phys. 2009, 130, 134113. (10) Ben-Nun, M.; Martinez, T. J. Ab Initio Quantum Molecular Dynamics. Adv. Chem. Phys. 2002, 121, 439−512. (11) Shalashilin, D. V. Quantum Mechanics with the Basis Set Guided by Ehrenfest Trajectories: Theory and Application to SpinBoson Model. J. Chem. Phys. 2009, 130, 244101. (12) Burghardt, I.; Giri, K.; Worth, G. A. Multimode Quantum Dynamics Using Gaussian Wavepackets: The Gaussian-Based Multiconfiguration Time-Dependent Hartree (G-MCTDH) Method Applied to the Absorption Spectrum of Pyrazine. J. Chem. Phys. 2008, 129, 174104. (13) Worth, G. A.; Robb, M. A.; Lasorne, B. Solving the TimeDependent Schrödinger Equation for Nuclear Motion in One Step: Direct Dynamics of Non-Adiabatic Systems. Mol. Phys. 2008, 106, 2077−2091. (14) Worth, G. A.; Robb, M. A.; Burghardt, I. A Novel Algorithm for Non-Adiabatic Direct Dynamics Using Variational Gaussian Wavepackets. Faraday Discuss. 2004, 127, 307−323. (15) Izmaylov, A. F. Perturbative Wave-Packet Spawning Procedure for Non-Adiabatic Dynamics in Diabatic Representation. J. Chem. Phys. 2013, 138, 104115. (16) Makhov, D. V.; Glover, W. J.; Martinez, T. J.; Shalashilin, D. V. Ab Initio Multiple Cloning Algorithm for Quantum Nonadiabatic Molecular Dynamics. J. Chem. Phys. 2014, 141, 054110. (17) Fernandez-Alberti, S.; Makhov, D. V.; Tretiak, S.; Shalashilin, D. V. Non-Adiabatic Excited State Molecular Dynamics of Phenylene Ethynylene Dendrimer Using a Multiconfigurational Ehrenfest Approach. Phys. Chem. Chem. Phys. 2016, 18, 10028−10040. (18) Mead, C. A.; Truhlar, D. G. On the Determination of BornOppenheimer Nuclear Motion Wave Functions Including Complications Due to Conical Intersections and Identical Nuclei. J. Chem. Phys. 1979, 70, 2284−2296. (19) Ryabinkin, I. G.; Izmaylov, A. F. Geometric Phase Effects in Dynamics Near Conical Intersections: Symmetry Breaking and Spatial Localization. Phys. Rev. Lett. 2013, 111, 220406. (20) Joubert-Doriol, L.; Ryabinkin, I. G.; Izmaylov, A. F. Geometric Phase Effects in Low-Energy Dynamics Near Conical Intersections: A Study of the Multidimensional Linear Vibronic Coupling Model. J. Chem. Phys. 2013, 139, 234103. (21) Ryabinkin, I. G.; Joubert-Doriol, L.; Izmaylov, A. F. When Do We Need to Account for the Geometric Phase in Excited State Dynamics? J. Chem. Phys. 2014, 140, 214116. (22) Xie, C.; Ma, J.; Zhu, X.; Yarkony, D. R.; Xie, D.; Guo, H. Nonadiabatic Tunneling in Photodissociation of Phenol. J. Am. Chem. Soc. 2016, 138, 7828−7831. (23) Joubert-Doriol, L.; Sivasubramanium, J.; Ryabinkin, I. G.; Izmaylov, A. F. Topologically Correct Quantum Nonadiabatic Formalism for On-The-Fly Dynamics. J. Phys. Chem. Lett. 2017, 8, 452−456. (24) Broeckhove, J.; Lathouwers, L.; Kesteloot, E.; Van Leuven, P. On the Equivalence of Time-Dependent Variational Principles. Chem. Phys. Lett. 1988, 149, 547−550. (25) McLachlan, A. D. A Variational Solution of the TimeDependent Schrödinger Equation. Mol. Phys. 1964, 8, 39−44. (26) Richings, G. W.; Polyak, I.; Spinlove, K. E.; Worth, G. A.; Burghardt, I.; Lasorne, B. Quantum Dynamics Simulations Using Gaussian Wavepackets: the vMCG Method. Int. Rev. Phys. Chem. 2015, 34, 269. (27) Izmaylov, A. F.; Mendive Tapia, D.; Bearpark, M. J.; Robb, M. A.; Tully, J. C.; Frisch, M. J. Nonequilibrium Fermi Golden Rule for 1797


Quantum Nonadiabatic Cloning of Entangled ... - ACS Publications

Recommend Documents