Semiclassical analysis of the electron-nuclear coupling in electronic non-adiabatic processes

0.55 MB·English

by  Federica Agostini

Semiclassicalanalysisoftheelectron-nuclearcouplinginelectronicnon-adiabaticprocesses
Federica Agostinia, Seung Kyu Mina and E. K. U. Grossa
aMax-Planck Institute of Microstructure Physics, Weinberg 2, D-06120 Halle, Germany
(Dated:March31,2015)
In the context of the exact factorization of the electron-nuclear wave function, the coupling be-
tweenelectronsandnucleibeyondtheadiabaticregimeisencoded(i)inthetime-dependentvector
and scalar potentials and (ii) in the electron-nuclear coupling operator. The former appear in the
5 Schro¨dinger-likeequationthatdrivestheevolutionofthenucleardegreesoffreedom, whereasthe
1
latter is responsible for inducing non-adiabatic effects in the electronic evolution equation. As we
0
havedevotedpreviousstudiestotheanalysisofthevectorandscalarpotentials,inthispaperwefo-
2
cusonthepropertiesoftheelectron-nuclearcouplingoperator,withtheaimofdescribinganumerical

r proceduretoapproximateitwithinasemiclassicaltreatmentofthenucleardynamics.
a
M
I. INTRODUCTION as a single product of a nuclear wave function and an
0
electronicfactor,thatparametricallydependsonthenu-
3
clearconfiguration. Severaladvantagesofthisreformu-
Modelling the dynamical coupling of electrons and
lation have been pointed out. First of all, it has been
]
nuclei beyond the Born-Oppenheimer (BO), or adia-
h shown [42, 43] that the nuclear wave function evolves
p batic, regime is currently among the most challenging
accordingtoamodifiedTDSEwhereatime-dependent
- problems in the fields of Theoretical Chemistry and
m vector potential and a time-dependent scalar potential
Condensed Matter Physics. Within the BO framework,
represent the effect of the electrons on the nuclei, be-
e molecularsystemsarevisualizedasasetofnucleimov-
yond the adiabatic regime. When a classical treatment
h ing on a single potential energy surface that represents
c ofthenucleardegreesoffreedomisintroduced,thecou-
the effect of the electrons in a given eigenstate. Many
. plingtotheelectronsisexactlyrepresentedbytheforce
s interesting phenomena, however, such as vision [1, 2],
c determinedfromthegradientofthetime-dependentpo-
charge separation in organic photovoltaic materials [3,
i tentials[43,44]. Ontheotherhand,theelectronicfactor
s 4] or Joule heating in molecular junctions [5, 6], occur
y evolves according to a (less standard) evolution equa-
in non-adiabatic conditions. In these situations, solv-
h tion, coupled to the nuclear TDSE, where the effect of
p ing exactly the time-dependent Schro¨dinger equation thenucleiisrepresentedbyanelectron-nuclearcoupling
[ (TDSE)forthecoupledsystemofelectronsandnuclear
operatorexplicitlydependingonthenuclearwavefunc-
is not feasible, as the computational cost scales expo-
2 tion.
nentiallywiththenumberofdegreesoffreedom. How-
v
ever,sincethefullquantumtreatmentrequirestorepre- Theanalysisofthetime-dependentpotentialsandof
7
1 sent the problem in terms of adiabatic states and tran- the classical nuclear force has been the subject of pre-
1 sitions among them in regions of strong non-adiabatic vious work. We have been able to analyze a simple
0 coupling, wave packet propagation techniques have model system to pinpoint some relevant features [42–
0 beendeveloped,retainingthequantumcharacterofthe 44]ofthepotentialsthatshouldbeaccountedforwhen
2. nuclear dynamics [7–20]. And these techniques are developing approximations. Based on these observa-
0 presentlythestateoftheartinquantumdynamicscom- tions,startingfromtheexactformulation,wehavepro-
5 putationalmethods,provingthebenchmarkforapprox- posed [45, 46] a novel mixed quantum-classical algo-
1 imate methods. In fact, for large systems, the dimen- rithm to solve the coupled electronic and nuclear evo-
:
v sionality of the problem does not allow to employ a lutionequationsinafullyapproximatedway. Theclas-
i quantummechanicaldescription,thustheonlyfeasible sical limit is considered as the lowest-order, in a (cid:126)-
X
approachistocombineaclassicaldescriptionforthenu- expansion,ofthenuclearwavefunctioninthecomplex-
r clei(orions)withaquantumtreatmentofafewotheres- phaserepresentation[47].
a
sentialdegreesoffreedom,e.g. electronsorprotons. In
In the present paper, the focus is directed towards
thiscontext,thequestionofhowtomodelthecoupling
the analysis of the electron-nuclear coupling operator
between the quantum and classical subsystems still re-
that, in the electronic equation, mediates the coupling
mains open, despite the fact that several schemes [20–
to the nuclei and that depends explicitly on (the gra-
39]havebeenproposedintheliteraturetryingtosettle
dient of) the nuclear wave function. It is fundamen-
thisissue.
tal to be able to correctly approximate such term, as it
In recent work, we have addressed this problem in isresponsibletoinduceelectronicnon-adiabatictransi-
the context of the exact factorization of the electron- tions[45,46]anddecoherence[48]. Whenaclassicalde-
nuclear wave function [40, 41]. In such a treatment of scription of the nuclei is adopted, the concept of wave
quantumdynamics,thesolutionoftheTDSEiswritten function is somehow lost and problems arise when ap-
proximating this operator. If a distribution of trajecto- canbewrittenastheproduct
ries [44] is used to mimic the evolution of the nuclear
wavefunction,itsmodulusandphasecannotbesmooth Ψ(r,R,t)=ΦR(r,t)χ(R,t), (2)
functionsofspace. Thenumericalerrorthusintroduced
affectsthecalculations,butitcanbecuredifrefinedap- ofthenuclearwavefunction,χ(R,t),andtheelectronic
proximations are considered. The goal of this paper is wave function, ΦR(r,t), which parametrically depends
todescribeaproceduretoavoidtheaboveissueandto ofthenuclearconfiguration[59,60].Throughoutthepa-
testitsefficiencyagainstexactcalculationsforasimple per the symbols r,R indicate the coordinates of the Ne
modelsystem. Therefore,weproposehere(i)toemploy electrons and Nn nuclei, respectively. Eq. (2) is unique
arepresentationofthenucleardensityintermsofevolv- underthepartialnormalizationcondition(PNC)
ingfrozengaussians(FGs)[49],ratherthantrajectories,
(cid:90) (cid:12) (cid:12)2
following the scheme presented in Ref. [44] and (ii) to dr(cid:12)Φ (r,t)(cid:12) =1 ∀R,t (3)
(cid:12) R (cid:12)
estimatethephaseofthenuclearwavefunctionadopt-
ingsuchaFGspicture,basedonasimplifiedformofthe
up to within a gauge-like phase transformation. The
semiclassicalHerman-Kluk[50–53]propagationscheme
evolutionequationsforΦ (r,t)andχ(R,t),
withintheinitialvaluerepresentation(IVR)theory[54– R
57]. The model system for non-adiabatic charge trans-
(cid:16) (cid:17)
ferofRef.[58]allowsforanexactnumericalsolutionof Hˆ −(cid:15)(R,t) Φ (r,t)=i(cid:126)∂ Φ (r,t) (4)
el R t R
the full quantum mechanical problem, thus providing
Hˆ χ(R,t)=i(cid:126)∂ χ(R,t), (5)
abenchmarktoanyapproximationthatwillbeconsid- n t
ered. Starting from these results, we will compute the
are derived by applying Frenkel’s action principle [61,
exact time-dependent scalar potential, also referred to
62] with respect to the two wave functions and are ex-
astime-dependentpotentialenergysurface(TDPES),in
actlyequivalent[40,41]totheTDSE(1). Eqs.(4)and(5)
a gauge where the vector potential can be set to zero.
areobtainedbyimposingthePNC[63,64]bymeansof
Theeffectoftheelectronsisthenfullyaccountedforby
Lagrangemultipliers.
thisTDPES,thatisadoptedtoevolvetheFGs. Sincethe
The electronic equation (4) contains the electronic
electronicpartoftheproblemissolvedexactly,theonly
Hamiltonian
source of error will be in this semiclassical approxima-
tion. It is worth stressing that this procedure does not Hˆ =Hˆ (r,R)+Uˆcoup[Φ ,χ], (6)
result in the development of a new algorithm, but is a el BO en R
testoftheperformanceoftheFGsapproximationofthe
which is the sum of the BO Hamiltonian and the ENC
nuclearmotion. operatorUˆcoup[Φ ,χ],
en R
The paper is organized as follows. In Section II we
bonrietflhye arencaallylstihseoffatchteoreilzeacttiroonn-fnourmclaealirsmcouanpdlinwge(EfoNcCus) Uˆcoup[Φ ,χ]=(cid:88)Nn 1 (cid:34)(cid:2)−i(cid:126)∇ν −Aν(R,t)(cid:3)2 (7)
en R M 2
operatorintheelectronicevolutionequation. SectionIII ν
ν=1
ipslodyeevdotiendthtoeacadlcisucluatsisoinosn.oNnutmheeraipcarolnrexsimulatstioanrespemre-- +(cid:18)−i(cid:126)χ∇νχ +Aν(R,t)(cid:19)(cid:0)−i(cid:126)∇ν −Aν(R,t)(cid:1)(cid:21).
sented in Section IV, by comparing situations with dif-
ferent non-adiabatic coupling strengths and testing the InEq.(4),(cid:15)(R,t)istheTDPES,definedas
approximations employed to evaluate the ENC term.
ConclusionsarestatedinSectionV. (cid:15)(R,t)=(cid:68)Φ (t)(cid:12)(cid:12)Hˆ −i(cid:126)∂ (cid:12)(cid:12)Φ (t)(cid:69) . (8)
R (cid:12) el t(cid:12) R
r
Uˆcoup and (cid:15)(R,t), along with the vector potential
en
II. THEEXACTFACTORIZATIONFRAMEWORK A(R,t),
(cid:68) (cid:12) (cid:69)
Intheabsenceofanexternalfield,thenon-relativistic A(R,t)= ΦR(t)(cid:12)(cid:12)−i(cid:126)∇νΦR(t) , (9)
Hamiltonian Hˆ = Tˆ +Hˆ describes a system of in- r
n BO
teractingnucleiandelectrons. Here,Tˆn denotesthenu- mediatethecouplingbetweenelectronsandnucleiina
clearkineticenergyandHˆBO(r,R)=Tˆe(r)+Vˆe,n(r,R) formally exact way. Here, the symbol (cid:104)·|·(cid:105)r stands for
istheBOHamiltonian,containingtheelectronickinetic anintegrationoverelectroniccoordinates.
energyTˆe(r)andallinteractionsVˆe,n(r,R). Asrecently The nuclear evolution is generated by the Hamilto-
proven [40, 41], the full wave function, Ψ(r,R,t), solu- nian
tionoftheTDSE
Hˆ (R,t)=(cid:88)Nn (cid:2)−i(cid:126)∇ν +Aν(R,t)(cid:3)2 +(cid:15)(R,t), (10)
HˆΨ(r,R,t)=i(cid:126)∂tΨ(r,R,t), (1) n ν=1 2Mν
2
accordingtotheTDSE(5). thesecondtermontheRHS,as|χ|isasmoothfunction
TheTDPESandthevectorpotentialareuniquelyde- ofthenuclearcoordinates. Thephaseinformationisin-
termined up to within gauge-like transformations [40, stead encoded in the classical action accumulated over
41]. Theuniquenesscanbestraightforwardlyprovenby timeandassociatedtoeachFG,aswewillshowinSec-
following the steps of the current density version [65] tionIII.
of the Runge-Gross theorem [66]. In this paper, as a Henceforth,wewilldropthebold-doubleunderlined
choiceofgauge, weintroducetheadditionalconstraint notationforelectronicandnuclearpositionsaswewill
(cid:0) (cid:1)
A R,t = 0(seeRef.[44]foradetaileddiscussionon dealwithone-dimensional(1D)quantities.
ν
howthisconditioncanbeimposed)[67].
Asdiscussedintheintroduction, westudytheprop-
erties of the ENC term −i(cid:126)∇ χ/χ that in the expres- III. SEMICLASSICALAPPROXIMATION
ν
sion of the operator Uˆcoup[Φ ,χ] explicitly depends on
en R
the nuclear wave function. This analysis is based on A. Nucleardensity
the interest in developing a procedure to approximate
itwhenaclassicalorsemiclassicaltreatmentofthenu- According to the procedure presented in Ref. [44], a
clear motion is adopted. For instance, in the classical setofindependentclassicaltrajectoriesevolvingonthe
limit,wehavederived[45,46]itsexpressionintermsof exact TDPES are able to reproduce the nuclear density
thenuclearmomentum. Thishasbeendonebywriting in almost perfect agreement with quantum results. As
the nuclear wave function as χ(R,t) = exp(cid:2)iS(R,t)(cid:3), pointedoutintheIntroduction,however,constructinga
(cid:126)
with S(R,t) a complex function [47]. If we now sup- histogram from the distribution of the trajectories does
pose [45, 46] that this function can be expanded as an not allow to compute the second term on the RHS of
asymptotic series in powers of (cid:126), namely S(R,t) = Eq.(13),thatinvolvesthegradientof|χ|,withoutalarge
(cid:80) (cid:126)αS (R,t),theENCtermbecomes numerical error. The reason is that the “classical” den-
α α
sityisnotasmoothfunctionofspace. Thesolutionpro-
−i(cid:126)∇ χ(R,t) posedhereistoimprovethepreviousapproximationof
ν
χ(R,t) =∇νS0(R,t), (11) nuclear dynamics, by propagating the mean positions
andmomentaofasetofFGsontheexactTDPES,rather
atthelowest-orderin(cid:126). Ontheright-hand-side(RHS), thanclassicaltrajectories.
thefunction∇ S (R,t)istheclassicalnuclearmomen- Given a set of Ntraj initial positions and momenta,
ν 0
tumevaluatedalongthetrajectory,since[45,46]S sat- R0,P0, sampled as described in Section IV, complex
0
gaussians, also referred to as coherent states, are con-
isfiesaHamilton-JacobiequationwithHamiltonian
structedas
Hn =(cid:88)Nn (cid:2)∇νS0(R,t2)M+Aν(R,t)(cid:3)2 +(cid:15)(R,t). (12) g(R;Rl,0Pl,0,γ)=(cid:16)πγ(cid:17)14 e−γ2(R−Rl,0)2+(cid:126)iPl,0(R−Rl,0),
ν
ν=1
(14)
Thevectorpotentialappearsintheaboveexpressionof
with width γ to be determined below. Each FG is also
theclassicalHamiltonianbecausethisresulthasgeneral
associateda“weight”,
validity,notonlyinthegaugeadoptedinthefollowing
calculations. (cid:90)
Alternatively, if the nuclear wave function is written wl = dRg∗(R;Rl,0Pl,0,γ)χ0(R), (15)
intermsofitsmodulusandphase,χ=|χ|eiS/(cid:126),the(ex-
act)expressionoftheENCtermbecomes corresponding to the projection of the initial nuclear
wave function χ (R) on the FGs. The nuclear density
−i(cid:126)∇νχ(R,t) −(cid:126)∇ν(cid:12)(cid:12)χ(R,t)(cid:12)(cid:12) ateachtimeisth0enobtainedas
χ(R,t) =∇νS(R,t)+i (cid:12)(cid:12)χ(R,t)(cid:12)(cid:12) . (13) Ntraj
|χ(R,t)|2 (cid:39) (cid:88) |w |2|g(R;R (t)P (t),γ)|2, (16)
l l l
It is clear at this point that a good estimate of the ENC
l=1
term is only possible when both the modulus and the
phase of the nuclear wave function are correctly de- where R (t),P (t) are the time-evolved positions and
l l
scribed. A classical treatment, as in Eq. (11), only pro- momentaofR ,P . Incomparisontothepurelyclas-
l,0 l,0
vides an approximation to the real part of the ENC sical approximation, Eq. (16) allows not only to repro-
term,whiletheinformationabouttheimaginarypartis duce the nuclear density in very good agreement with
lost[68]. quantum results, as will be shown in Section IV, but
In the following, we will introduce a FG-based ap- alsotocalculate(analytically)thegradientofthenuclear
proach to determine an approximation to Eq. (13). density(orofthemodulus,asitappearsinEq.(13)).
FGs [49] evolving on the exact TDPES are used to re- It is important to notice that Eq. (16) is an approxi-
constructthenucleardensity,thusallowingtocalculate mation to the nuclear density when the nuclear wave
3
function is represented as a superposition of coherent use a FGs approximation rather than a rigorous semi-
states. Infact,coherentstatesformanovercompleteba- classicalapproximation.
sis and, in writing Eq. (16), we neglect the overlaps of Intheprocedureemployedhere,thesemiclassicalnu-
coherent states. The reason for this further approxima- clearwavefunctionisestimatedasasumovertrajecto-
tionisrelatedtothechoiceoftheinitialsetofpositions riesoftime-evolvedcoherentstates,namely
and momenta, the mean positions and mean momenta
oftheFGs. Ononehand,wewanttomaintainthesame N(cid:88)traj e(cid:126)iSl
choice of initial conditions done for the classical prop- χFG(R,t)= wl 2π(cid:126)g(R;Rl(t),Pl(t),γ), (19)
agation (see Section IV), in order to be able to directly l=1
compare classical and FG results. On the other hand,
wherethel-thclassicalactioniscalculatedas
we want to obtain an initial nuclear density as close as
possible to the exact density. We have thus computed (cid:90) t (cid:18)P2(τ) (cid:19)
S = dτ l −(cid:15)(R (τ),τ) . (20)
therootmeansquaredeviation(RMSD)forthenuclear l 2M l
0
density when either (1) Eq. (16) is employed or (2) the
full expression is considered, i.e. with overlap terms. Here, the (exact) TDPES is evaluated, at time τ, at the
Thebestagreementachievedincase(1)isforγ =7.0a−02 classicalpositionRl(τ)anditsexpressionisobtainedby
(usedinSectionIV),withRMSD=0.020,0.017,0.006for solving the TDSE for the full wave function Ψ and by
N =2000,5000,10000,respectively,whereasincase directly calculating Eq. (8) once the factorization (2) is
traj
(2) is for γ = 3.0 a−2 with RMSD = 0.037, 0.024 0.0013. applied.
0
It is evident that a better agreement at the initial time
is obtained for case (1), namely when the approxima-
tioninEq.(16)isused,forallthevaluesofN . Once IV. NUMERICALRESULTS
traj
theparametersdefiningthecoherentstatesareselected,
namely R ,P ,γ, the weights w associated to each FG The expression of the TDPES according to Eq. (8) is
l l l
areautomaticallydeterminedbyEq.(15)andkeptcon- determinedbycalculatingtheelectronicwavefunction
stantthroughoutthepropagation. ΦR(r,t)fromthefullwavefunctionΨ(r,R,t), whichis
knownatalltimesbysolvingtheTDSE(1)forthemodel
Hamiltonian
B. Nuclearphase
1 ∂2 1 ∂2 1 1
Hˆ(r,R)=− − + +
Thephaseofthenuclearwavefunctionwillbedeter- 2∂r2 2M ∂R2 (cid:12)(cid:12)L2 −R(cid:12)(cid:12) (cid:12)(cid:12)L2 +R(cid:12)(cid:12)
minedaccordingto (21)
(cid:18) (cid:19) (cid:18) (cid:19)
(cid:16) (cid:17) |r−L| |r+L|
S(R,t)(cid:39)arctan(cid:61)(cid:60)[[χχSC((RR,,tt))]], (17) −erf |RR−fr| − erf Rr2 − erf Rl2 .
SC |R−r| (cid:12)(cid:12)r− L(cid:12)(cid:12) (cid:12)(cid:12)r+ L(cid:12)(cid:12)
2 2
where χ (R,t) is a semiclassical approximation to
SC
the exact χ(R,t). Following the Herman-Kluk proce- ThissystemhasbeenintroducebyShinandMetiu[58]
dure [50–53] to approximate the quantum propagator, as a prototype for non-adiabatic charge transfer. The
the expression of the time-evolved wave function at systemis1Dandconsistsofthreeionsandasingleelec-
timetis tron, as depicted in Fig. 1. Two ions are fixed at a dis-
(cid:90) dR dP
χSC(R,t)= 20π(cid:126) 0Ct(R0,P0)e(cid:126)iS(R0,P0;t) ❘
r
(cid:104)R|R P ,γ(cid:105)(cid:104)R P ,γ|χ (cid:105).
t t 0 0 0
(18)
Here, (cid:104)R P ,γ|χ (cid:105) denotes the projection of the initial (cid:0)☎✆ ✵ ✂❡✂✝✞✟☎✆
0 0 0 ❢(cid:0)✁✂✄(cid:0)☎✆ ❢(cid:0)✁✂✄(cid:0)☎✆
nuclear wave function on the coherent states, similarly
to Eq. (15), while (cid:104)R|R P ,γ(cid:105) is an alternative expres-
t t
sion for the coherent states (see Eq. (14)). R ,P are ▲
t t
the (classically) evolved positions and momenta corre- FIG. 1. Schematic representation of the model system de-
sponding to the initial conditions R ,P . S(R ,P ;t) scribedbytheHamiltonian(21).
0 0 0 0
is the classical action accumulated up to time t along
the trajectory whose initial conditions are R ,P . The tance of L = 19.0 a , the third ion and the electron are
0 0 0
Herman-Klukpre-factorisindicatedherewiththesym- free to move in 1D along the line joining the two fixed
bol C (R ,P ) [50–53], but it will be set equal to unity ions. Here, the symbols r and R are the coordinates
t 0 0
throughout the calculations. Therefore, the symbol oftheelectronandthemovableionmeasuredfromthe
χ (R,t) will be replaced by χ (R,t), since we will centerofthetwofixedions. Theionicmassischosenas
SC FG
4
M = 1836,theprotonmass,whereastheotherparame- Hamilton’sequations
ters are tuned in order to make the system essentially
a two-electronic-state model. We present here the re- R˙l(t)=Pl(t)/M; P˙l(t)=−∂R(cid:15)(R,t)(cid:12)(cid:12)Rl(t), (22)
sultsobtainedbychoosingtwodifferentsetsofparam-
eters, producing strong and weak non-adiabatic cou- integrated by using the velocity-Verlet algorithm with
plings, between the first, (cid:15)(1), and the second BOPES, a time-step of 12.0 × 10−3 fs (or 0.5 a.u.). The initial
BO
(cid:15)(2), around the avoided crossing at R = −1.90 a . conditions are sampled from the Wigner phase-space
BO ac 0 distribution corresponding to the initial nuclear wave
The values of the parameters in the Hamiltonian (21)
function. The initial coherent states used in Eq. (15)
are: R = 5.0 a , R = 3.1 a and R = 4.0 a , for
f 0 l 0 r 0 arethusconstructedtodeterminethe(complex-valued)
the strong coupling case; R = 3.8 a , R = 2.9 a and
f 0 l 0 weights w . Numerical results have been obtained for
R =5.5a ,fortheweakcouplingcase.TheBOsurfaces l
r 0 sets of N = 2000,5000,10000 trajectories, with the
areshowninFig.2(upperpanels). traj
value [70] γ = 7.0 a−2 for the width of the FGs. The
0
coherentstatesusedtorepresentthenuclearwavefunc-
strong coupling weak coupling
tion form an overcomplete basis, therefore the sum in
Eq. (16) has to be truncated. In order to choose the ad-
-0.1
equatenumberofbasisfunctionstoincludeinthesum,
()εεBOh --00..32 εε|BB χOO |(( 122)) teatwhnlleeevrceaganlysuetechNsoenotresfaxeNjarvct=rtaat5nji.0ou0Tnc0hleh.eaWarrseedsbhueeanletvnssiewttyceioslatmlenbdpdeuapittnreseddsaectphonpentrefiRodrxMmoimnSeDladytibffoooenrr-
-6 -3 0 3 6 -6 -3 0 3 6 in Eq. (16) for a set of values of γ and we have chosen
R (a ) R (a ) thevalueofthisparameterforwhichtheRMSDismini-
0 0
1 mum(seealsothediscussionattheendofSectionIIIA),
n butclearlyothervaluescanbeselected.
o
ati ρ1 We will confirm below that the semiclassical FG
pul 0.5 ρ2 scheme, adapted to the factorization approach, is an
o
p
O accurate and efficient way to approximate the nuclear
B
wave function. Before presenting the FG results, let us
0
first show different snapshots taken along the dynam-
0 10 20 30 0 10 20 30
ics, showing the TDPES. Its feature have been exten-
time (fs) time (fs)
sively discussed in previous work [42–44], but for the
FIG.2.Upperpanels:LowestfourBOsurfaces,asfunctionsof sake of completeness, we report here a few configura-
thenuclearcoordinate,forstrong(left)andweak(right)non- tions.Moreover,wewouldliketounderlinethatherere-
adiabatic coupling strengths. The first (red line) and second sultsfordifferentnon-adiabaticcouplingstrengthswill
(greenline)surfacescorrespondtotheadiabaticstatesthatare bepresented,inordertotesttheefficiencyofthemethod
populatedduringthedynamics,whereasthethirdandfourth
alsointheweakcouplingregime.
surfaces (dashed black lines) are shown for reference. The
Fig. 3 shows the gauge-invariant (GI) and gauge-
squared modulus (reduced by ten times and rigidly shifted
dependent (GD) components of the TDPES, namely
in order to superimpose it on the energy curves) of the ini-
the two terms that can be identified in Eq. (8) as
tialnuclearwavepacketisalsoshown(thinblackline).Lower
panels: Populationsofthetwolowestadiabaticstates(ρ1and (cid:15)GI(cid:0)R,t(cid:1) = (cid:104)ΦR(t)|Hˆel|ΦR(t)(cid:105)r and (cid:15)GD(cid:0)R,t(cid:1) =
ρ2)asfunctionsoftime. Thearrowsrepresentthetime-steps (cid:104)ΦR(t)|−i(cid:126)∂t|ΦR(t)(cid:105)r. The snapshots shown in Fig. 3
showninthefollowingfigures.
are taken along the evolution at the times indicated by
the arrows in Fig. 2 (lower panels) and the two lowest
adiabatic surfaces are shown for reference. As previ-
Westudythetimeevolutionofthissystembychoos- ously discussed [42–44], before the splitting of the nu-
ing the initial wave function as the product of a real- clearwavepacket,theGIpartoftheTDPESisdiabatic,
valued normalized Gaussian wave pack√et, centered at whereas the GD part is constant, and after the split-
Rc = −4.0 a0 with variance σ = 1/ 2.85 a0 (thin ting the GI component develops steps that bridge be-
black line in Fig. 2, upper panels), and the second BO tweendifferentadiabaticsurfaces,whereastheGDpart
electronic state. To calculate the TDPES, we first solve is piecewise constant. The TDPES is only calculated in
the TDSE (1) for the complete system, with Hamilto- the regions where the nuclear density is (numerically)
nian (21), and obtain the full wave function, Ψ(r,R,t). notzero. Thelackofreliableinformationbeyondthere-
This is done by numerical integration of the TDSE us- gionsshowninthefiguresisnotanissuewhenclassical
ing the split-operator-technique [69], with time-step of trajectories or FGs are employed to mimic the nuclear
2.4×10−3fs(or0.1a.u.). density, as the regions where the exact density is expo-
The mean positions and momenta of the FGs evolve nentiallysmallarenot,orarepoorly,sampled.
along classical trajectories, generated according to ItisevidentfromFig.3thatwewilldiscussresultsfor
5
strong coupling weak coupling strong coupling weak coupling
0
t = 4.84 fs |χ|2
ε)h-0.1 GD 1 |χFG|2
y ( histogram
erg-0.2 GI 0.5
n
e-0.3 |χ|2 εBO(1 ,2) t = 4.84 fs
0
0
t = 14.52 fs
1
ε)h-0.1
gy ( 0.5
er-0.2
n
e
-0.3 t = 14.52 fs 0
0 0.6 t = 30.25 fs
ε)h-0.1
y ( 0.3
g
er-0.2
n
e 0
-0.3 t = 30.25 fs -6 -3 0 3 6 -6 -3 0 3 6
-6 -3 0 3 6 -6 -3 0 3 6 R (a0) R (a0)
R (a ) R (a )
0 0
FIG.4. Nucleardensityforstrong(left)andweak(right)non-
FIG.3. TDPESatdifferenttime-steps,asindicatedbythear- adiabatic coupling strengths. Exact results (cyan) are com-
rowsinFig.2,forstrong(left)andweak(right)non-adiabatic pared with the semiclassical density (red), expressed as sum
coupling strengths. The two lowest BO surfaces (black) are ofFGs, andwiththehistogram(green)constructedfromthe
plottedforreference. TheGIpartoftheTDPES(red)presents distributionofclassicalpositionsalongthetrajectories.
the well-studied [42, 43] dynamical steps that bridge piece-
wise adiabatic shapes potential energy surface, whereas the
GIpart[44](cyan)ispiecewiseconstant. Thenucleardensity
(blue)isshownforreference.
tions (4) and (5). They only represent a benchmark for
anyquantum-classicalalgorithm, sincetheeffectofthe
electrons, via the TDPES, is treated exactly. The semi-
classical density, constructed as the weighted sum of
shortdynamics,limitedtothefirsthalfoftheoscillation FGsgiveninEq.(16),isalsoaccurate. Incomparisonto
periodofthenuclearwavepacketinthepotentialwell. theclassicalhistogram,thegainhereisthesmoothness
Interestingdynamicsmayariseatlatertimes,whenfor of the density, not achievable with purely classical tra-
instancethenuclearwavepacketcrossesasecondtime jectories. Thisfeatureisextremelyimportantforthecal-
the non-adiabatic coupling region. However, here we culationoftheENCtermcontainingthegradientofthe
focusontheinitialnon-adiabaticeventandtesthowthe nucleardensity,viatheterm−i(cid:126)∇ |χ|/|χ|inEq.(13).
ν
FGapproximationcapturethisprocess. Duetothefact Thesecondimportantcharacteristicofthesemiclassi-
that the analysis reported below is the first attempt to cal approach is that each FG contributes a phase factor
incorporate semiclassically nuclear quantum effects in tothefullnuclearwavefunction,thusallowingtodeter-
theexactfactorizationformalism,westudyasimplesit- minealsothefirsttermontheRHSofEq.(13). Weshow
uation that nonetheless captures the main features of thisterm,i.e. ∇ S(R,t),inFig.5,comparingonceagain
ν
a non-adiabatic event. Also, situations where, for in- exact results with the corresponding FG and classical
stance, reproducing tunneling dynamics might repre- approximations. The semiclassical value of ∇ S(R,t)
ν
sentaproblemforclassicallyevolvingFGs[71–74]will is determined as the gradient of the phase in Eq. (17),
bethesubjectoffurtherstudyandnotaddressedhere. whereasthefunction∇ S(R,t)isclassicallyinterpreted
ν
A first set of results is shown in Fig. 4, where we as the nuclear momentum evaluated along each trajec-
compare the nuclear density, |χ(R,t)|2, from three dif- tory. WhiletheagreementbetweenquantumandFGre-
ferentcalculations: quantum(cyan),employingthefull sultsisremarkable,thephase-spacepointscorrespond-
electron-nuclearwavefunction;classical(green),where ing to the classical trajectories do not allow to recon-
the histogram is constructed from the distribution of struct a smooth function of R. Even if the number of
classical trajectories evolving on the TDPES according classicaltrajectoriesisincreased,thephase-spacepoints
to Eqs. (22) (as in Ref. [44]); FG (red), with the nu- are too “noisy” to allow for reconstructing a smooth
clear density given in Eq. (16). As expected from pre- function. A smoothing algorithm should then be em-
vious calculations [44], the use of classical trajectories ployed,butthenumericalefficiencyofthewholeproce-
seems to be enough accurate to reproduce the nuclear duremightbecomequestionable.
density. It is however important to stress again, that Fig.6showstheimaginarypartoftheENCtermfrom
theseresultsarenotobtainedbysolvingafullyapprox- Eq. (13) at different time-steps during the dynamics, as
imateformofthecoupledelectronicandnuclearequa- in previous figures. The exact results (cyan) are com-
6
strong coupling weak coupling
simplifiedformoftheHerman-Klukpropagator,where
20
thepre-factorissetto1,semiclassicalresultsareinsat-
FG
C] 15 t = 4.84 fs trajeecxtaocrties isfactorygoodagreementwithexactresults.
EN 10 One could consider improving the approach devel-
Re [ 5 oped here by, for instance, explicitly computing the
|χ|2 Herman-Kluk pre-factor, at the expenses of increasing
0
thecomputationalcost.
-5 -4 -3 -2 -5 -4 -3 -2
20
15
C]
N V. CONCLUSION
E 10 t = 14.52 fs
e [
R 5
We have reported our first semiclassical procedure
0 adaptedtotheformalismoftheexactfactorizationofthe
-2 -1 0 -2 -1 0
electron-nuclear wave function [40, 41]. The approach
20
has been used to estimate the ENC term that explic-
NC] 15 t = 30.25 fs itly depends on the nuclear wave function (gradient of
e [E 10 its modulus and phase) in the electronic equation. In
R
5 previous work [45, 46] on the development of a mixed
0 quantum-classical algorithm in the context of the exact
2 4 6 2 4 6 factorization,suchtermhasbeentreatedfullyclassically
R (a) R (a)
0 0 andidentifiedasthenuclearmomentum. However,we
FIG.5. RealpartoftheENCtermforstrong(left)andweak observedthat,despitethefactthatthisapproximationis
(right)non-adiabaticcouplingstrengths. Thenucleardensity widelyusedintheliterature[23,26,27],correctionterms
(black) is shown for reference. Exact results (cyan) are com- naturally arise in the factorization framework. The ex-
paredwithsemiclassicalcalculations(red)andwiththeclassi- tentandroleofthecorrectionsmaydependonthepro-
calphase-spacepoints(green). cess that we intend to study, therefore we have shown
in the present paper how to evaluate such corrections
basedonthesemiclassicalpropagationofFGs,forasim-
ple model case of non-adiabatic charge transfer. The
strong coupling weak coupling
10 maingaininusingFGs,ratherthanapurelyclassicalap-
t = 4.84 fs FG
exact proach,liesinthepossibilityofobtainingsmoothfunc-
C] |χ|2 tionswhosegradientscanbeeasilydeterminedwithout
N
m [E 0 introducinglargenumericalerrors.
I The semiclassical results shown in the paper have
been obtained by evolving FGs on the exact TDPES,
-10
-5 -4 -3 -2 -5 -4 -3 -2 thatisknownforthesimplemodelsystemstudiedhere
10
as the outcome of exact calculations based on the nu-
t = 14.52 fs
C] merical solution of the full TDSE. This procedure is a
N
E 0 test, since the only approximation is the semiclassical
m [ treatment of the nuclear dynamics, whereas the elec-
I
trons are treated exactly, via the information encoded
-10
-2 -1 0 -2 -1 0 in the TDPES. Moreover, this study provides a bench-
10 mark for future development, aiming at improving the
t = 30.25 fs
initial,andlowest-order,mixedquantum-classicalalgo-
C]
N rithmderivedfromthefactorization[45,46]. Theproce-
E 0
m [ duredescribedherecanbeeasilyimplementedinanal-
I
gorithm,aswewillpresentelsewhere[48],resultingina
-10 novelmixedquantum-semiclassicalschemeforsolving
2 4 6 2 4 6
coupledelectron-nucleardynamics.
R (a) R (a)
0 0
FIG.6. ImaginarypartoftheENCtermforstrong(left)and
weak(right)non-adiabaticcouplingstrengths.Thecolorcode ACKNOWLEDGEMENTS
isthesameasinFig.5.
The authors would like to thank Neepa T. Maitra for
her help in improving the presentation of the results.
Partial support from the Deutsche Forschungsgemein-
pared only with the approximation (red) based on the schaft (SFB 762) and from the European Commission
semiclassical propagation of FGs. Even if we employ a (FP7-NMP-CRONOS)isgratefullyacknowledged.
7
[1] D. Polli, P. Altoe`, O. Weingart, K. M. Spillane, C. Man- [32] B.F.E.Curchod,I.Tavernelli, andU.Rothlisberger,Phys.
zoni,D.Brida,G.Tomasello,G.Orlandi,P.Kukura,R.A. Chem.Chem.Phys.13,3231(2011).
Mathies, M. Garavelli, and G. Cerullo, Nature 467, 440 [33] R.E.Wyatt,C.L.Lopreore, andG.Parlant,J.Chem.Phys.
(2010). 114,5113(2001).
[2] W.C.Chung,S.Nanbu, andT.Ishida,J.Phys.Chem.B [34] I.Burghardt,J.Chem.Phys.122,094103(2005).
116,8009(2012). [35] O.V.PrezhdoandC.Brooksby,Phys.Rev.Lett.86,3215
[3] C.A.Rozzi,S.M.Falke,N.Spallanzani,A.Rubio,E.Moli- (2001).
nari, D. Brida, M. Maiuri, G. Cerullo, H. Schramm, [36] N.ZamsteinandD.J.Tannor,J.Chem.Phys.137,22A518
J. Christoffers, and C. Lienau, Nat. Communic. 4, 1602 (2012).
(2013). [37] D. Bousquet, K. H. Hughes, D. A. Micha, and
[4] A.E.Jailaubekov,A.P.Willard,J.R.Tritsch,W.-L.Chan, I.Burghardt,J.Chem.Phys.134,064116(2011).
N.Sai,R.Gearba,L.G.Kaake,K.J.Williams,K.Leung, [38] J.O.RichardsonandM.Thoss,J.Chem.Phys.139,031102
P.J.Rossky, andX.-Y.Zhu,Nat.Mater.12,66(2013). (2013).
[5] E. J. McEniry, D. R. Bowler, D. Dundas, A. P. Horsfield, [39] N.Ananth,J.Chem.Phys.139,124102(2013).
C.G.Sa´nchez, andT.N.Todorov,J.Phys.:Condens.Mat- [40] A. Abedi, N. T. Maitra, and E. K. U. Gross, Phys. Rev.
ter19,196201(2007). Lett.105,123002(2010).
[6] A. P. Horsfield, D. R. Bowler, A. J. Fisher, and T. N. [41] A.Abedi,N.T.Maitra, andE.K.U.Gross,J.Chem.Phys.
Todorov,J.Phys.:Condens.Matter16,3609(2004). 137,22A530(2012).
[7] H.-D. Meyer, U. Manthe, and L. S. Cederbaum, Chem. [42] A.Abedi,F.Agostini,Y.Suzuki, andE.K.U.Gross,Phys.
Phys.Lett.165,73(1990). Rev.Lett110,263001(2013).
[8] U.Manthe,H.-D.Meyer, andL.S.Cederbaum,J.Chem. [43] F.Agostini,A.Abedi,Y.Suzuki, andE.K.U.Gross,Mol.
Phys.97,3199(1992). Phys.111,3625(2013).
[9] I.Burghardt,H.-D.Meyer, andL.S.Cederbaum,J.Chem. [44] F.Agostini,A.Abedi,Y.Suzuki,S.K.Min,N.T.Maitra,
Phys.111,2927(1999). and E. K. U. Gross, arXiv:1406.4667 [physics.chem-ph]
[10] H.-D.MeyerandG.A.Worth,Theor.Chim.Acta109,251 (2014).
(2003). [45] A.Abedi,F.Agostini, andE.K.U.Gross,Europhys.Lett.
[11] M. Thoss, W. Domcke, and H. Wang, Chem. Phys. 296, 106,33001(2014).
217 (2004). [46] F.Agostini,A.Abedi, andE.K.U.Gross,J.Chem.Phys.
[12] H.-D.M.L.WangandV.May,J.Chem.Phys.125,014102 141,214101(2014).
(2003). [47] J.H.V.Vleck,Proc.Nat.Ac.Sci.14,178(1928).
[13] T. Westermann, R. Brodbeck, W. S. Alexander B. [48] S.K.Min,F.Agostini, andE.K.U.Gross, (tobesubmit-
Rozhenko and, and U. Manthe, J. Chem. Phys. 135, ted).
184102(2011). [49] E.J.Heller,J.Chem.Phys.75,2923(1981).
[14] H.WangandM.Thoss,J.Chem.Phys.119,1289(2003). [50] M.F.HermanandE.Kluk,Chem.Phys.91,27(1984).
[15] J.Li,I.Kondov,H.Wang, andM.Thoss,J.Phys.Chem. [51] E.Kluk,M.F.Herman, andH.L.Davis,J.Chem.Phys.
C114,18481(2010). 84,326(1986).
[16] Q.Meng,S.Faraji,O.Vendrell, andH.-D.Meyer,J.Chem. [52] K.G.Kay,Chem.Phys.322,3(2006).
Phys.137,134302(2012). [53] W.H.Miller,Mol.Phys.100,397(2002).
[17] M. Schro¨der, J.-L. C. n Macedo, and A. Brown*, Phys. [54] W.H.Miller,J.Chem.Phys.53,3578(1970).
Chem.Chem.Phys.10,850(2008). [55] M. Thoss and H. Wang, Ann. Rev. Phys. Chem. 55, 299
[18] T.J.Mart´ınezandR.D.Levine,Chem.Phys.Lett.259,252 (2004).
(1996). [56] K.G.Kay,Ann.Rev.Phys.Chem.56,255(2005).
[19] T. J. Mart´ınez, M. Ben-Nun, and R. D. Levine, J. Phys. [57] W.H.Miller,J.Phys.ChemA105,2942(2001).
Chem.100,7884(1996). [58] S.ShinandH.Metiu,J.Chem.Phys.102,23(1995).
[20] T.J.Mart`ınez,Acc.Chem.Res.39,119(2006). [59] G.Hunter,Int.J.QuantumChem9,237(1975).
[21] P.Ehrenfest,Zeitschriftfu¨rPhysik45,455(1927). [60] N.I.GidopoulosandE.K.U.Gross,Phil.Trans.R.Soc.A
[22] P.Pechukas,Phys.Rev.181,174(1969). 372,20130059(2014).
[23] J.C.Tully,J.Chem.Phys.93,1061(1990). [61] J.Frenkel,Wavemechanics,Clarendon,Oxforded.(1934).
[24] G.StockandM.Thoss,Phys.Rev.Lett.78,578(1997). [62] A.D.McLachlan,Mol.Phys.8,39(1964).
[25] X.SunandW.H.Miller,J.Chem.Phys.106,6346(1997). [63] J. L. Alonso, J. Clemente-Gallardo, P. Echeniche-Robba,
[26] A.W.Jasper,C.Zhu,S.Nangia, andD.G.Truhlar,Fara- andJ.A.Jover-Galtier,J.Chem.Phys.139,087101(2013).
dayDiscuss.127,1(2004). [64] A.Abedi,N.T.Maitra, andE.K.U.Gross,J.Chem.Phys.
[27] S.Nielsen,R.Kapral, andG.Ciccotti,J.Chem.Phys.112, 139,087102(2013).
6543(2000). [65] S.K.GhoshandA.K.Dhara,Phys.Rev.A38,1149(1988).
[28] C.-Y. Hsieh and R. Kapral, J. Chem. Phys. 138, 134110 [66] E. Runge and E. K. U. Gross, Phys. Rev. Lett. 52, 997
(2013). (1984).
[29] S. Bonella and D. F. Coker, J. Chem. Phys. 122, 194102 [67] S.K.Min,A.Abedi,K.S.Kim, andE.K.U.Gross,Phys.
(2005). Rev.Lett.113,263004(2014).
[30] P.HuoandD.F.Coker,J.Chem.Phys.137,22A535(2012). [68] Tobeprecise,thenucleardensityisapproximated,within
[31] N. L. Doltsinis and D. Marx, Phys. Rev. lett. 88, 166402 theclassicaltreatment,asδ-function,centeredatalltimes
(2002). at the classical trajectory. But this contribution is totally
8
omittedintheexpressionoftheENCterm. [72] J.Ankerhold,Quantumtunnellingincomplexsystems.The
[69] M.D.Feit,F.A.FleckJr., andA.Steiger,J.Comput.Phys. semiclassicalapproach,SpringerTractsinModernPhysics,
47,412(1982). Vol.224(Springer-VerlagBerlinHeidelberg,2007).
[70] Thevalueofγischosentohavethebestoverlapbetween [73] M. Saltzer and J. Ankerhold, Phys. Rev. A 68, 042108
theinitialsemiclassicaldensityandtheexactdensity, in (2003).
ordertoimposeinitialconditionsascloseaspossibleto [74] K.G.Kay,Phys.Rev.A88,012122(2013).
theexactones.
[71] D.GelmanandS.D.Schwartz,J.Chem.Phys.129,024504
(2008).
9