Analysis of geometric phase effects in the quantum-classical Liouville formalism
Ilya G. Ryabinkin, Chang-Yu Hsieh, Raymond Kapral, Artur F. Izmaylov
aa r X i v : . [ qu a n t - ph ] D ec Analysis of geometric phase effects in the quantum-classical Liouvilleformalism
Ilya G. Ryabinkin,
Chang-Yu Hsieh, Raymond Kapral, and Artur F. Izmaylov
1, 2 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 (Dated: 12 September 2018)
We analyze two approaches to the quantum-classical Liouville (QCL) formalism that differ in the order of twooperations: Wigner transformation and projection onto adiabatic electronic states. The analysis is carriedout on a two-dimensional linear vibronic model where geometric phase (GP) effects arising from a conicalintersection profoundly affect nuclear dynamics. We find that the Wigner-then-Adiabatic (WA) QCL approachcaptures GP effects, whereas the Adiabatic-then-Wigner (AW) QCL approach does not. Moreover, the Wignertransform in AW-QCL leads to an ill-defined Fourier transform of double-valued functions. The double-valuedcharacter of these functions stems from the nontrivial GP of adiabatic electronic states in the presence of aconical intersection. In contrast, WA-QCL avoids this issue by starting with the Wigner transform of single-valued quantities of the full problem. Since the WA-QCL approach uses solely the adiabatic potentials andnon-adiabatic derivative couplings as an input, our results indicate that WA-QCL can capture GP effectsin general two-state crossing problems using first-principles electronic structure calculations without priordiabatization or introduction of explicit phase factors.
I. INTRODUCTION
Molecular electronic adiabatic surfaces often crossforming degenerate manifolds of nuclear configurationswith the topology of conical intersections (CIs).
CIsare the most common triggers of radiationless transi-tions that drive photo-induced chemistry, and trans-fers of electronic energy and charge. Besides facil-itating electronic transitions CIs change the topologyof the nuclear subspace: if adiabatic electronic wavefunctions are infinitely slowly (adiabatically) transportedaround a closed loop that encircles the CI seam, they ac-quire an extra ( −
1) phase. This is the geometric phase(GP) that makes the electronic wave functions double-valued functions of nuclear coordinates.
In order tohave a single-valued total electron-nuclear wave function,the nuclear wave functions have to compensate for signchanges in their electronic counterparts; hence, a nuclearSchr¨odinger equation must be solved with double-valuedboundary conditions. The extra phase causes interfer-ence between parts of the nuclear wave function, whichresults in different vibronic spectra and nuclear dynamicsas compared to the case without explicit account for theGP.
The double-valued character of the nuclear wavefunctions in the adiabatic representation is challengingfor numerical simulations. Switching to the diabaticrepresentation makes the nuclear diabatic wave func-tion single-valued. However, the diabatic representationwithin a finite electronic subspace cannot be obtained for a general polyatomic system ( N atoms > To remain in the unambiguous adiabatic representation,Mead and Truhlar proposed to compensate for the GPof individual nuclear states by attaching an extra phase factor e iλ ( R ) , where λ ( R ) is a function that increases by π on encircling a closed path around the CI seam. Thisapproach was successfully implemented and appliedto many real molecules by Kendrick. However, thedefinition of λ ( R ) is based on diabatic model considera-tions, and its application to a general molecular systemassumes some diabatic model potential. An alternativeapproach to the definition of λ ( R ) was suggested by Baerand coworkers. They related ∇ λ ( R ) to the matrixelements of derivative couplings d (1)12 ( R ), thus removingany arbitrariness in the definition of λ ( R ). However, thissuggestion was challenged by Kendrick, Mead, and Truh-lar who found that Baer’s approach inconsistently ne-glects some terms, and this inconsistency leads one toquestion the validity of the final result. Therefore, ac-counting for GP effects for general systems using onlyresults of electronic structure calculations seems to bechallenging within a fully quantum approach.In this context, it is interesting that the quantum-classical Liouville (QCL) formalism is able to captureGP effects using only adiabatic input: energies andnon-adiabatic couplings. The QCL framework uses theWigner transform (WT) of the nuclear degrees of free-dom (DOF) to arrive at a mixed quantum-classical de-scription. Derivations of the QCL equations for non-adiabatic problems using the adiabatic representation forelectronic DOF can proceed along two paths that differin the order in which the WT and the projection to theadiabatic basis are applied. We will denote by WA-QCLthe approach where the WT is done first, and by AW-QCL the approach where the adiabatic projection pre-cedes the WT. It is not a priori obvious which ofthe two approaches is better. Previously, the WA-QCLand AW-QCL approaches were used to simulate the spin-boson model and no significant differences were found. In this work we analyze the efficacy of these approachesto describe the dynamics in a two-dimensional (2D) lin-ear vibronic coupling (LVC) model. As has been shownrecently , the 2D LVC model exhibits prominent GP ef-fects, such as a distinct interference pattern due to theGP as well as a strong impact of the GP on nuclear dy-namics. This model allows us to assess both QCL ap-proaches with respect to CI-introduced topological fea-tures in non-adiabatic dynamics.The paper is organized as follows: first, we briefly illus-trate the emergence of double-valued adiabatic electronicwave functions in the 2D LVC model. Then, we outlinethe derivations of the matrix form of the WA-QCL equa-tion starting from the full electron-nuclear density ma-trix and the AW-QCL equation starting from a projec-tion of the Liouville equation onto an adiabatic electronicbasis. Comparing these two approaches we identifya term which is responsible for GP effects. We concludeour paper with numerical results and their analysis.
II. GEOMETRIC PHASE IN A TWO-DIMENSIONALLINEAR VIBRONIC COUPLING MODEL
The two-dimensional LVC is a prototype for systemscontaining a CI and exhibiting nontrivial GP effects. Theelectron-nuclear Hamiltonian of the model isˆ H = ˆ T N + ˆ V ˆ V ˆ V ˆ V ! , (1)where ˆ T N = − ~ ( ∂ /∂x + ∂ /∂y ) is the nuclear kineticenergy operator, ˆ V and ˆ V are the diabatic potentialsrepresented by identical 2D parabolas shifted in the x -direction by a and coupled by the ˆ V potential:ˆ V = ω (cid:20)(cid:16) x + a (cid:17) + y (cid:21) , ˆ V = cy, (2)ˆ V = ω (cid:20)(cid:16) x − a (cid:17) + y (cid:21) . (3)Electronic DOF are represented as position-independentdiabatic states | i and | i in a two-dimensional electronicsubspace.Transformation to the adiabatic representation is madeby diagonalizing the potential matrix in Eq. (1) by meansof the unitary transformation U U ( θ ) = cos θ − sin θ sin θ cos θ ! , (4)where θ is a mixing angle between the diabatic states, andthe factor is introduced for convenience. The adiabaticelectronic wave functions are related to the columns of U as | ψ adi1 i = cos θ | i + sin θ | i (5) | ψ adi2 i = − sin θ | i + cos θ | i . (6) The angle θ is determined by the matrix elements of thediabatic Hamiltonian (1) as θ = arctan 2 V V − V , (7)and changes by 2 π for any closed path in the nuclear( x, y ) subspace that encircles the CI point. Consideringthat U (2 π ) = − , both adiabatic wave functions | ψ adi i i acquire an extra ( − ) sign, which is a consequence of theGP. To have a single-valued total wave function | Ψ i = | χ adi1 i | ψ adi1 i + | χ adi2 i | ψ adi2 i , (8)a sign change must be also imposed on the adiabatic nu-clear wave functions | χ adi1 , i .This consideration can be extended to the total densityoperator ˆ ρ = | Ψ i h Ψ | , which is also a single-valued func-tion of the nuclear ( x, y ) coordinates. In contrast, thenuclear density matrix in the adiabatic representationˆ X αα ′ = (cid:10) ψ adi α (cid:12)(cid:12) ˆ ρ (cid:12)(cid:12) ψ adi α ′ (cid:11) = | χ adi α i h χ adi α ′ | , (9)is a complicated object due to the double-value boundaryconditions imposed on | χ adi α i and h χ adi α ′ | .To model dynamics of the nuclear adiabatic densitiesˆ X we solve the corresponding Liouville equation i ~ ∂ ˆ X ∂t = [ ˆ H adi , ˆ X ] . (10)Here, the Hamiltonian ˆ H adi is obtained by projectingthe electronic DOF in Hamiltonian (1) onto the | ψ adi1 , i subspace:ˆ H adi = ˆ T N + ˆ τ + ˆ W − ˆ τ ˆ τ ˆ T N + ˆ τ + ˆ W + ! , (11)whereˆ τ ij = − ~ [ h ψ adi i |∇ ψ adi j i · ∇ + 12 h ψ adi i |∇ ψ adi j i ] ≡ − ~ [ d (1) ij · ∇ + 12 d (2) ij ] . (12)with ∇ = ( ∂/∂x, ∂/∂y ). The non-adiabatic couplings, d (1) ij and d (2) ij are the vector and scalar derivative couplingmatrix elements, and ˆ W ± are the matrix elements of theadiabatic potential ˆW = U ˆVU † :ˆ W ± = ˆ V + ˆ V ± q ( ˆ V − ˆ V ) + 4 ˆ V . (13)Equation (10) has to be solved by imposing double-value boundary conditions for ˆ X αα ′ due to the pres-ence of the nontrivial GP. Single-valued boundary con-ditions lead to a completely different dynamics as shownin Fig. 1.To avoid GP complications in this model one can sim-ulate dynamics for the diabatic nuclear densityˆ ρ αα ′ = h α | ˆ ρ | α ′ i = | χ dia α i h χ dia α ′ | , (14) P t, a.uDouble−valued BCSingle−valued BC FIG. 1. The solutions of Eq. (10) with double- and single-valued boundary conditions. The initial wavepacket is a sim-ple Gaussian located in the left well; P is the fraction of thewavepacket located in the left well. The model parametersare ω = 2 , a = 3 , c = 4. by solving the corresponding Liouville equation i ~ ∂ ˆ ρ ∂t = [ ˆ H , ˆ ρ ] . (15)Of course, this diabatic path cannot be strictly followedfor real molecules where the diabatic representation isnot rigorously defined. But in this work we will use thediabatic representation to generate the exact quantumdynamics for our model problem. III. QUANTUM-CLASSICAL LIOUVILLE FORMALISM
In what follows we review the main steps of the WA-QCL and AW-QCL derivations to see the conse-quences of the double-valued character of the adiabaticelectronic states in both approaches. We consider only asingle set of nuclear coordinates associated with a parti-cle of mass M to simplify our discussion and to remainclose to that in Sec. II; generalization to include addi-tional nuclear coordinates is straightforward. A. Wigner-then-adiabatic path
In line with Ref. 40 we start with the WT of the nu-clear DOF in the Liouville equation for the total densitymatrix in the coordinate representation ˆ ρ ( r , R , r ′ , R ′ ) =Ψ ∗ ( r , R )Ψ( r ′ , R ′ ), i ~ ∂ ˆ ρ W ∂t = ( ˆ H ˆ ρ ) W − (ˆ ρ ˆ H ) W . (16)Here, we used the Wigner representation for operatorsˆ A W ( R c , P ) = Z D R c + s (cid:12)(cid:12)(cid:12) ˆ A (cid:12)(cid:12)(cid:12) R c − s E e − i ~ s · P d s , (17) where R c = ( R + R ′ ) / s = R − R ′ , and P is a parameter that can beassociated with the classical momentum in the classicallimit. Operator products (cid:16) ˆ A ˆ B (cid:17) W can be transformedfurther using the Wigner-Moyal operator e − i ~ ˆΛ / as (cid:16) ˆ A ˆ B (cid:17) W = ˆ A W e − i ~ ˆΛ / ˆ B W , (18)where ˆΛ is the Poisson bracket operator ˆΛ = ←−∇ P · −→∇ R c −←−∇ R c ·−→∇ P , and the arrows indicate the directions in whichthe differential operators act. Using identity (18) theWigner-transformed Liouville equation becomes i ~ ∂ ˆ ρ W ∂t = ˆ H W e − i ~ ˆΛ / ˆ ρ W − ˆ ρ W e − i ~ ˆΛ / ˆ H W . (19)The ˆΛ operator acts only on the nuclear DOF, thus, tointroduce a semi-classical description of nuclear dynamicswe expand the exponent in a Taylor series and keeponly terms linear in ~ : i ~ ∂ ˆ ρ W ∂t = ˆ H W (cid:16) − i ~ ˆΛ / (cid:17) ˆ ρ W − ˆ ρ W (cid:16) − i ~ ˆΛ / (cid:17) ˆ H W . (20)The partially Wigner-transformed H W derived fromEq. (1) is quadratic in R and P , hence, the Wigner-transformed Liouville equation (20) is exact for our modelcase.Both the Wigner-transformed density matrix ˆ ρ W andHamiltonian ˆ H W in Eq. (20) are still quantum operatorsin the electronic subspace. Let us project the electronicDOF on the adiabatic electronic eigenfuctions | ψ adi α ( R c ) i i ~ ∂ ˆ ρ αα ′ W ∂t = h ψ adi α | ˆ H W (cid:18) − i ~ (cid:19) ˆ ρ W | ψ adi α ′ i− h ψ adi α | ˆ ρ W (cid:18) − i ~ (cid:19) ˆ H W | ψ adi α ′ i . (21)Inserting the resolution-of-the-identity operator in theelectronic subspace ˆ I = P β | ψ adi β ( R c ) i h ψ adi β ( R c ) | be-tween the operator products, we express Eq. (21) in termsof operator matrices defined as A W ≡ h ψ adi α | ˆ A W | ψ adi α ′ i ∂ ρ W ∂t = − i ~ [ H W , ρ W ] − (cid:16) H W ˆΛ ρ W − ρ W ˆΛ H W (cid:17) + 12 (cid:16) [ D , H W ] · ∇ P ρ W + ∇ P ρ W · [ D , H W ] (cid:17) − P M · [ D , ρ W ] , (22)where D is the matrix of the vector derivative cou-plings d (1) αβ . Equation (22) accounts for mutual inter-actions between the electronic and nuclear DOF. Notethat all Wigner-transformed operators depend on a sin-gle nuclear coordinate R c , and even though the adia-batic functions | ψ adi α ( R c ) i are double-valued, the diago-nal density matrix elements h ψ adi α | ˆ ρ W | ψ adi α i are all single-valued functions of R c . As for the off-diagonal elements h ψ adi α | ˆ ρ W | ψ adi β i (coherences), in general, they are double-valued functions, because the GP related sign flip of thebra and ket components may not be simultaneous forany closed contour in the R c subspace. An example ofsuch situation is a three-state model where two electronicstates | ψ adi1 , i have a CI, and the third state | ψ adi3 i has anavoided crossing with the first two states. In this exam-ple, the matrix elements h ψ adi α | ˆ ρ W | ψ adi3 i α =1 , are double-valued functions because the sign change in h ψ adi1 , | is notcompensated by that in | ψ adi3 i on encircling the CI. Thus,generally, Eq. (22) must be simulated with double-valuedboundary conditions for the off-diagonal elements. How-ever, for problems where the parametric dependence ofall adiabatic states leads to the same GP change after en-circling any closed contour all matrix elements in Eq. (22)are single-valued functions, and Eq. (22) does not needdouble-valued boundary conditions. Our 2D LVC modelis an example where the latter condition is satisfied be-cause both functions | ψ adi1 ( R c ) i and | ψ adi2 ( R c ) i eitherchange their sign or do not depending on whether theCI point is enclosed by the contour. At any case, allquantities in Eq. (22) are well defined mathematically,and their physical interpretation is well-known. B. Adiabatic-then-Wigner path
We start by projecting the Liouville equation onto theadiabatic electronic basis | ψ adi α i to obtain i ~ ∂ ˆ X αα ′ ∂t = X β (cid:16) ˆ H adi αβ ˆ X βα ′ − ˆ X αβ ˆ H adi βα ′ (cid:17) , (23)where ˆ H adi αβ and ˆ X αβ are multi-state generalizations ofthe two-state adiabatic Hamiltonian (11) and nucleardensity (9). Applying the WT to Eq. (23) we obtain i ~ ∂ X W ∂t = H adi W e − i ~ ˆΛ / X W − X W e − i ~ ˆΛ / H adi W . (24)For the double-valued adiabatic densities ˆ X αα ′ the WT X αα ′ W = Z D R c + s (cid:12)(cid:12)(cid:12) ˆ X αα ′ (cid:12)(cid:12)(cid:12) R c − s E e − i ~ s · P d s (25)is not well defined, because, as is shown in the Appendix,this operation is equivalent to the Fourier transform ofa double-valued function. Unless special care is taken atthis step the double-valued character of the ˆ X αα ′ densitywill be lost. Here, we follow previous derivations ofthe AW-QCL equation, noting that it is inconsistent withsubtleties arising from a nontrivial GP.The Wigner-transformed adiabatic Hamiltonian H adi W in Eq. (24) differs from the H W matrix only by the WT ofthe non-adiabatic coupling operators ˆ τ αβ [see Eq. (12)], H adi W = H W + τ W , (26) where the τ W matrix has elements τ αα ′ W = − i ~ P M · d (1) αα ′ + ~ M (cid:16) ∇ R c · d (1) αα ′ − d (2) αα ′ (cid:17) . (27)The last term in the parentheses defines a positive-definite matrix K with elements K αα ′ = ~ M (cid:16) ∇ · d (1) αα ′ − d (2) αα ′ (cid:17) = ~ M h∇ ψ adi α |∇ ψ adi α ′ i . (28)The explicit form of the adiabatic Hamiltonian obtainedfrom Eqs. (26–28) is H adi W = H W − i ~ P M · D + ~ K . (29)Expanding the Wigner-Moyal exponent in the Taylor se-ries and keeping only first order terms in ~ in Eq. (24)we obtain the AW-QCL equation, ∂ X W ∂t = − i ~ [ H W , X W ] − (cid:16) H W ˆΛ X W − X W ˆΛ H W (cid:17) − P M · [ D , X W ] . (30)It is important to note that Eq. (30) is no longer an ex-act evolution equation for the quantum dynamics of our2D LVC model. This is in contrast to Eq. (20) whichis fully equivalent to the quantum Liouville equation forthis model. The difference stems from two sources: First,the adiabatic Hamiltonian [Eq. (29)] is not a quadraticfunction of the nuclear coordinates, and consequently thetruncation of the Wigner–Moyal exponent is no longer ex-act. Second, the AW-QCL derivation ignores the double-valued character of the nuclear density, which resultsin Eq. (30) that corresponds to the evolution of a nu-clear density with the incorrect (single-valued) boundaryconditions. Comparing Eqs. (22) and (30) reveals thatboth sources of the difference between the AW-QCL andWA-QCL approaches can be related to the third term inEq. (22). Since this term accounts for direct changes inthe nuclear momenta as a result of coupling to the elec-tronic DOF, one might expect significant influence of thisterm on the nuclear dynamics. Moreover, as we shall seebelow, neglecting this term leads to loss of GP relatedfeatures in the nuclear dynamics. IV. NUMERICAL SIMULATIONS
To quantify the difference between the WA-QCL andAW-QCL approaches and to assess the relative impor-tance of GP effects we simulate non-adiabatic dynamicsfor the 2D LVC model. Our choice of initial conditionsand the model parameters aims to minimize differencesbetween the WA-QCL and AW-QCL approaches that arerelated to the omission of the higher order quantum con-tributions in AW-QCL. We employ the following set ofparameters for the 2D LVC model [Eqs. (1–3)]: ω = 2, P t, a.u DiabaticAdiabatic (no GP)WA−QCLAW−QCL FIG. 2. The fraction of the probability density remaining inthe region x < P ) for the 2D LVC model simulated withdifferent methods. a = 1, c = 4, M = 1. The initial density is taken as aproduct of two Gaussian wave packets centred at the min-imum of the | i diabatic surface ( − a/ ,
0) with widths σ = p /ω and an initial momentum (5 / ,
0) (motiontowards the positive x values). This setup results in arelatively high initial energy and reduction of quantumtunnelling. We consider four different approaches to non-adiabatic dynamics: a) quantum diabatic [Eq. (15)], b)quantum adiabatic without GP [Eq. (10)], c) WA-QCL,d) AW-QCL. We investigate two properties: the timeevolution of the fraction of the nuclear density remainingin the region x < and use of the “momentum jump”approximation for the last two terms in Eq. (22). Bothapproaches predict almost full transfer of the densityfrom the initial x < x > . t = 1 . x >
0. This prominent feature is completely ab-sent in the quantum adiabatic dynamics [Fig. 3(b)] and,thus, can be seen as a manifestation of the GP.
Moreover, the adiabatic model predicts completely spu-rious constructive interference between parts of the wavepacket that skirt the CI point from opposite sides, giv-
FIG. 3. Snapshots of the square root of the probability densityat t = 1.6 a.u. for the 2D LVC model dynamics: a) fullyquantum (diabatic) , b) quantum adiabatic (no GP), c) WA-QCL, d) AW-QCL. x g (x) g (x) ( g + g ) g + g ( g − g ) FIG. 4. Total density of two Gaussian wave packets g and g in cases of (gold) constructive interference, (red) destructiveinterference, (green) no interference. ing rise to a peak in the nuclear density at x ≈ . y = 0 line, however, this dip is notas pronounced as for the models accounting for the GP.The AW-QCL density profile suggests that this approachdoes not produce any nuclear interference. To supportthis assertion we considered different relations betweentwo spatially separated Gaussian wave packets. Figure 4shows that in the absence of interference the total den-sity may have a central minimum, which, however, is notas deep as in the case of destructive interference. Thisresult also indicates that electronic transitions constitutethe dominant quantum effect in the AW-QCL dynamics,while quantum nuclear interference within each electronicstate is negligible. V. CONCLUSIONS
The analysis of the QCL dynamics for the 2D LVCmodel presented above showed that the WA-QCL equa-tion is the only approach where complications arisingfrom a non-trivial GP in the adiabatic representationdo not make formalism ill-defined. In situations whenthe GP behavior of the electronic adiabatic states in-volved in the dynamics is similar (e.g., CI of two elec-tronic states), the WA-QCL approach can treat GP ef-fects exactly without imposing double-valued boundaryconditions. In contrast, the AW-QCL approach involvesthe WT integral of the adiabatic density operator witha double-valued kernel. Practically, this transformationis equivalent to the Fourier transformation of a double-valued function, which is not a well-defined procedure.Yet, if one proceeds with the Fourier transformation ig-noring the double-valued character of the quantum den-sity, the resulting AW-QCL equation describes the quan-tum density with the incorrect single-valued boundaryconditions. Thus, the AW-QCL approach is plaguedby the same problems as the quantum dynamics in theadiabatic representation without including the GP (seeFig. 1).Interestingly, the AW-QCL and WA-QCL equationsdiffer only by one term. This term is responsible fordifferences associated with the Moyal–Wigner exponenttruncation in different representations and GP effects.Further separation of the GP-related terms from thosecorresponding to the Moyal–Wigner exponent truncationdoes not seem to be feasible. However, this analysis sug-gests why the previous numerical assessment of the twoapproaches did not find a significant difference betweenresults for the spin-boson model where GP effects areabsent.It is worth noting that in contrast to the quantumadiabatic picture, where the GP appears in the form ofnon-trivial boundary conditions, the WA-QCL approachconfronts the problem of the GP in much milder fashion.The WT maps similar GP behavior of different electronicstates into the third term of the WA-QCL dynamicalequation (22). Thus, the only components that are leftwith the double-valued boundary conditions are coher-ences of electronic states that have different GP behav-iors. An interesting question for future investigation iswhether the trajectory based techniques used to simulateEq. (22) will be able to handle the double-valued charac-ter of the coherences in general multi-level systems owingto the space locality of a propagated object?It should be also emphasized that the WA-QCL ap-proach, though accounting for the GP in our model, can-not be reduced to the Mead and Truhlar treatment of the GP. This can be easily seen by considering a singleelectronic state version of Eq. (22) that describes classicaldynamics of nuclear DOF on a given electronic potentialwithout any GP terms, which inevitably appear in theMead and Truhlar treatment.All quantities that are required to perform WA-QCLdynamics, such as the adiabatic potentials and thederivative couplings, are readily available from first-principles quantum chemistry calculations. Thus, theWA-QCL formalism is probably the only formalism thatcan naturally account for GP effects in two electronicstate crossing problems with on-the-fly quantum chem-istry calculations. ACKNOWLEDGMENTS
RK and AFI acknowledge support from the NaturalScience and Engineering Research Council (NSERC) ofCanada through the Discovery Grants Program.
Appendix A: Wigner transformation of adiabatic nucleardensities
Here we consider the WT of the adiabatic densities, X αα ′ W ( P ; R c ) = Z D R c + s (cid:12)(cid:12)(cid:12) ˆ X αα ′ (cid:12)(cid:12)(cid:12) R c − s E e − i ~ s · P d s = Z χ α (cid:16) R c + s (cid:17) ∗ χ α ′ (cid:16) R c − s (cid:17) e − i ~ s · P d s = Z f αα ′ ( s ; R c ) e − i ~ s · P d s . (A1)Equation (A1) illustrates that the WT can be formu-lated as the Fourier transform (FT) of the s variablefor functions f αα ′ ( s ; R c ). In what follows we will showthat for the 2D LVC model the f αα ′ ( s ; R c ) are double-valued functions of the s variable, and therefore, the FTin Eq. (A1) is not well defined.Even for our model we do not have an exactanalytical representation of the f αα ′ ( s ; R c ) func-tions, but we know that the total density functionˆ ρ = P αα ′ | ψ adi α ′ ( R c − s / i f αα ′ ( s ; R c ) h ψ adi α ( R c + s / | is always a single-valued function. Thus, if f αα ′ ( s ; R c )are double-valued then | ψ adi α ′ ( R c − s / i h ψ adi α ( R c + s / | must be too, and vice versa . Hence, we will show that theelectronic components | ψ adi α ′ ( R c − s / i h ψ adi α ( R c + s / | are double-valued functions of the s parameter. Thisis not a trivial check because the double-valued char-acters of the bra and ket components can potentiallycompensate each other. We will demonstrate by explicitconstruction that there exist at least some closed con-tours for the s parameter which, for a fixed R c , causesonly | ψ adi α ′ ( R c − s / i or h ψ adi α ( R c + s / | to change itssign.For the 2D LVC model the functions | ψ adi α ′ ( R c − s / i and h ψ adi α ( R c + s / | are defined by Eqs. (5) and (6). To FIG. 5. Two contours z ± s of Eq. (A2): the red contour ( z − s )encircles the CI point, and the blue contour ( z + s ) does not.Parameters were set to z c = x c + iy c = 1, d = 10, r = 9. f / pq / p FIG. 6. θ - φ relation for the contours on Fig. 5: the red (blue)dots follow the contour that encircles (does not encircle) theCI point. Parameters of the Hamiltonian (1) are chosen sothat θ = arctan( y/x ) construct desired contours we fixed R c = ( x c , y c ) andobtained the ( x, y ) coordinates of the R c − s / R c + s / z ± s = z c ± ( d + re iφ ) / . (A2)Figure 5 shows that this choice generates two contourswhich are topologically different with respect to the CIpoint. Thus, the contour that encompasses the CI pro-duces a change in the mixing angle θ [Eq. (7)] by 2 π ,while moving along the other contour returns θ to its ini-tial value (see Fig. 6). Adiabatic transport along thesecontours will lead to a sign change in all products of theelectronic components | ψ adi α ′ ( R c − s / i h ψ adi α ( R c + s / | and, thus, this concludes our illustration of the double-valued character of the f αα ′ ( s ; R c ) functions. A. Migani and M. Olivucci, in
Conical Intersection ElectronicStructure, Dynamics and Spectroscopy , edited by W. Domcke,D. R. Yarkony, and H. K¨oppel (World Scientific, Singapore,2004) p. 271. D. R. Yarkony, Rev. Mod. Phys. , 985 (1996). L. Blancafort, P. Hunt, and M. A. Robb, J. Am. Chem. Soc. , 3391 (2005). L. Blancafort, F. Jolibois, M. Olivucci, and M. A. Robb, J. Am.Chem. Soc. , 722 (2001). A. F. Izmaylov, D. Mendive-Tapia, M. J. Bearpark, M. A. Robb,J. C. Tully, and M. J. Frisch, J. Chem. Phys. , 234106 (2011). H. C. Longuet-Higgins, U. Opik, M. H. L. Pryce, and R. A.Sack, Proc. R. Soc. A , 1 (1958). M. V. Berry, Proc. R. Soc. A , 45 (1984). B. Simon, Phys. Rev. Lett. , 2167 (1983). F. S. Ham, Phys. Rev. Lett. , 725 (1987). M. V. Berry, Proc. R. Soc. A , 31 (1987). I. J. R. Aitchison, Phys. Scr. , 12 (1988). J. Sch¨on and H. K¨oppel, Chem. Phys. Lett. , 55 (1994). J. Sch¨on and H. K¨oppel, J. Chem. Phys. , 9292 (1995). I. G. Ryabinkin and A. F. Izmaylov, Phys. Rev. Lett. , 220406(2013). L. Joubert-Doriol, I. G. Ryabinkin, and A. F. Izmaylov, J. Chem.Phys. (2013), in press, arXiv:1310.2929. L. S. Cederbaum, in
Conical Intersections , edited by W. Dom-cke, D. R. Yarkony, and H. K¨oppel (World Scientific, Singapore,2004) pp. 3–40. M. Baer, Chem. Phys. Lett. , 112 (1975). H. K¨oppel and B. Schubert, Mol. Phys. , 1069 (2006). T. Van Voorhis, T. Kowalczyk, B. Kaduk, L.-P. Wang, C.-L.Cheng, and Q. Wu, Ann. Rev. Phys. Chem. , 149 (2010). A. Sirjoosingh and S. Hammes-Schiffer, J. Chem. Theory Com-put. , 2831 (2011). D. R. Yarkony, Chem. Rev. , 481 (2012). J. E. Subotnik, R. J. Cave, R. P. Steele, and N. Shenvi, J. Chem.Phys. , 234102 (2009). D. Opalka and W. Domcke, J. Chem. Phys. , 154108 (2010). A. Troisi and G. Orlandi, J. Chem. Phys. , 5356 (2003). B. N. Papas, M. S. Schuurman, and D. R. Yarkony, J. Chem.Phys. , 124104 (2008). H. Nakamura and D. G. Truhlar, J. Chem. Phys. , 10353(2001). C. A. Mead and D. G. Truhlar, J. Chem. Phys. , 2284 (1979). B. Kendrick and C. A. Mead, J. Chem. Phys. , 4160 (1995). B. K. Kendrick, C. Alden Mead, and D. G. Truhlar, Chem. Phys. , 31 (2002). B. K. Kendrick, in
Conical Intersections. Electronic Structure,Dynamics and Spectroscopy , Advanced Series in Physical Chem-istry, Vol. 15, edited by W. Domcke, D. R. Yarkony, andH. K¨oppel (World Scientific, Singapore, 2003) Chap. 12, pp. 521–553. B. Kendrick and R. T. Pack, J. Chem. Phys. , 7475 (1996). B. Kendrick and R. T. Pack, J. Chem. Phys. , 7502 (1996). B. Kendrick, Phys. Rev. Lett. , 2431 (1997). B. K. Kendrick, J. Phys. Chem. A , 6739 (2003). M. Baer, J. Chem. Phys. , 2694 (1997). Z. Xu, M. Baer, and A. J. C. Varandas, J. Chem. Phys. ,2746 (2000). B. K. Kendrick, C. A. Mead, and D. G. Truhlar, J. Chem. Phys. , 7594 (1999). R. Kapral, Annual Review of Physical Chemistry , 129 (2006). A. Kelly and R. Kapral, J. Chem. Phys. , 084502 (2010). R. Kapral and G. Ciccotti, J. Chem. Phys. , 8919 (1999). I. Horenko, C. Salzmann, B. Schmidt, and C. Sch¨utte, J. Chem.Phys. , 11075 (2002). K. Ando and M. Santer, J. Chem. Phys. , 10399 (2003). K. ˙Imre, E. ¨Ozizmir, M. Rosenbaum, and P. F. Zweifel, J. Math.Phys. , 1097 (1967). A more physically appealing derivation proceeds by scaling theequation so that an expansion of the Liouville operator in smallparameter ( m/M ) / , where m and M are the characteristicmasses of the quantum subsystem and bath DOF, can be car-ried out ( see Ref 40). D. M. Kernan, G. Ciccotti, and R. Kapral, The Journal of Phys- ical Chemistry B , 424 (2008). A. Ferretti, G. Granucci, A. Lami, M. Persico, and G. Villani, J. Chem. Phys.104