Geometric phase-induced nuclear quantum interference is robust against quantum dissipation

TL;DR

This study demonstrates through numerical simulations that geometric phase-induced nuclear quantum interference remains robust even in the presence of complex quantum dissipation, highlighting its potential stability in realistic environments.

Contribution

The paper introduces a numerically exact method combining local diabatic representation and hierarchical equations of motion to analyze geometric phase effects under dissipation.

Findings

Destructive interference persists despite non-Markovian quantum dissipation.

Robustness is confirmed for both vibrational and electronic environments.

Provides an intuitive path integral-like explanation for the observed robustness.

Abstract

One of the intriguing effects due to conical intersections is the geometric phase, manifested as destructive quantum interference in the nuclear probability distribution. However, whether such geometric phaseinduced interference can survive in dissipative environments remains an open question. We demonstrate by numerically exact dissipative conical intersection dynamics simulations that the destructive interference is highly robust against non-Markovian quantum dissipation. To do so, we integrate the recently proposed local diabatic representation to describe vibronic couplings and the hierarchical equations of motion for system-bath interactions. Both vibrational and electronic environments are considered. An intuitive path integral-like picture isprovided to explain the robustness of geometric phase-induced interference.

Equations

$$\Psi (\mathbf{r},\mathbf{R},t)=\sum _{\alpha }{\varphi }_{\alpha }(\mathbf{r};\mathbf{R}){\chi }_{\alpha }(\mathbf{R},t).$$

$$H_{\mathrm{BO}}(\mathbf{R})|{\varphi }_{\alpha }(\mathbf{R})⟩=V_{\alpha }(\mathbf{R})|{\varphi }_{\alpha }(\mathbf{R})⟩$$

$$\Psi (\mathbf{r},\mathbf{R},t)$$

$$\equiv \sum _{\mathbf{n}\alpha }{C}_{\mathbf{n}\alpha }(t)|\mathbf{n}\alpha ⟩$$

$$i{\dot{C}}_{\mathbf{m}\beta }(t)=V_{\mathbf{m}\beta }{C}_{\mathbf{m}\beta }(t)+\sum _{\mathbf{n},\alpha }{T}_{\mathbf{mn}}{A}_{\mathbf{m}\beta ,\mathbf{n}\alpha }{C}_{\mathbf{n}\alpha }(t).$$

$$A_{\mathbf{m}\beta ,\mathbf{n}\alpha }=⟨{\varphi }_{\beta }({\mathbf{R}}_{\mathbf{m}})|{\varphi }_{\alpha }({\mathbf{R}}_{\mathbf{n}}){⟩}_{\mathbf{r}},$$

$$e^{i{\theta }_{\alpha }}=$$

$$\cdots |{\varphi }_{\alpha }({\mathbf{R}}_N){⟩}_{\mathbf{r}}⟨{\varphi }_{\alpha }({\mathbf{R}}_N)|{\varphi }_{\alpha }({\mathbf{R}}_1){⟩}_{\mathbf{r}}$$

$$=$$

$$G({\mathbf{R}}_f,{\alpha }_f,t;{\mathbf{R}}_0,{\alpha }_0,t_0)\equiv ⟨{\mathbf{R}}_f{\alpha }_f|e^{-iHt}|{\mathbf{R}}_0{\alpha }_0⟩$$

$$=\sum _{\{{\alpha }_j,{\mathbf{R}}_j\}}\prod _{j=0}^{L-1}⟨{\mathbf{R}}_{j+1}{\alpha }_{j+1}|e^{-i({\hat{T}}_{\mathrm{N}}+{\hat{H}}_{\mathrm{BO}})\Delta t}|{\mathbf{R}}_j{\alpha }_j⟩$$

$$=\sum _{\{\xi (t)\}}W[\xi (t)]e^{iS[\xi (t)]}$$

$$\displaystyle S[\xi(t)]=-\sum_{j=0}^{N_{t}}\Big{(}K_{{\mathbf{R}}_{j+1},{\mathbf{R}}_{j}}+V_{\alpha_{j}}({\mathbf{R}}_{j})\Big{)}\Delta t,$$

$$W[\xi (t)]=\prod _j{A}_{{\xi }_{j+1},{\xi }_j}.$$

$$\mathcal{A}({\mathbf{R}}_f,{\mathbf{R}}_0)=e^{iS}(W[{\Gamma }_{-}]+W[{\Gamma }_{+}])=e^{iS}W[{\Gamma }_{-}](1+e^{i{\theta }_{\Gamma }})=0.$$

$$H_{\mathrm{I}}=Q\otimes X.$$

$$C(t)=⟨X(t)X(0){⟩}_{\text{B}},$$

$$C(t)=\frac{1}{\pi }{\int }_{-\infty }^{\infty }\mathrm{d}\omega \frac{e^{-i\omega t}J(\omega )}{1-e^{-\beta \omega }}.$$

$$J(\omega )\equiv \frac{1}{2}{\int }_{-\infty }^{\infty }\mathrm{d}t{e}^{i\omega t}⟨[X(t),X(0)]{⟩}_{\text{mbox}B}.$$

$$\hat{\rho }(t)=\sum _{\mathbf{n}{\mathbf{n}}^{\prime },\alpha {\alpha }^{\prime }}{\rho }_{\mathbf{n}\alpha ,{\mathbf{n}}^{\prime }{\alpha }^{\prime }}(t)|\mathbf{n}\alpha ⟩⟨{\mathbf{n}}^{\prime }{\alpha }^{\prime }|=\sum _{\xi ,{\xi }^{\prime }}{\rho }_{\text{mbox}M}(\xi ,{\xi }^{\prime },t)|\xi ⟩⟨{\xi }^{\prime }|$$

$$\rho (\xi ,{\xi }^{\prime },t)=\sum _{{\xi }_0,{\xi }_0^{\prime }}\mathcal{G}(\xi ,{\xi }^{\prime },t;{\xi }_0,{\xi }_0^{\prime },t_0){\rho }_{\text{mbox}M}({\xi }_0,{\xi }_0^{\prime },t_0)$$

$$\mathcal{G}(\xi ,{\xi }^{\prime },t;{\xi }_0,{\xi }_0^{\prime },t_0)=$$

$$\sum _{\{\xi (t),{\xi }^{\prime }(t)\}}W[\xi (t)]e^{iS[\xi (t)]}\mathcal{F}[Q(\xi (t)),Q({\xi }^{\prime }(t))]\overline{W}[{\xi }^{\prime }(t)]e^{-iS[{\xi }^{\prime }(t)]}$$

$$\mathcal{F}[Q(\xi (t)),Q({\xi }^{\prime }(t))]=\exp \{$$

$$\displaystyle\sum^{j}_{l}\big{[}C(t_{j}-t_{l})Q[\xi_{l}]-C^{*}(t_{j}-t_{l})Q[\xi^{\prime}_{l}]\big{]}\Delta t\Big{\}}$$

$$⟨X^{\text{mbox}B}(t)X^{\text{mbox}B}(0){⟩}_{\text{mbox}B}=\sum _{k=1}^K{\eta }_k{e}^{-{\gamma }_{k}t}.$$

$$⟨X^{\text{mbox}B}(0)X^{\text{mbox}B}(t){⟩}_{\text{mbox}B}=\sum _{k=1}^K{\eta }_k^{\ast }{e}^{-{\gamma }_k^{\ast }t}\equiv \sum _{k=1}^K{\eta }_{\bar{k}}^{\ast }{e}^{-{\gamma }_{k}t}.$$

$${\dot{\rho }}_{\mathbf{L}}^{(L)}=$$

$$\displaystyle-i\sum_{k}L_{k}\big{(}\eta_{k}{Q}\rho^{(L-1)}_{{\mathbf{L}}_{k}^{-}}-\eta_{\bar{k}}^{\ast}\rho^{(L-1)}_{{\mathbf{L}}_{k}^{-}}{Q}\big{)}.$$

$$H_{\text{mbox}M}={\hat{T}}_{\text{N}}\mathbf{I}+\mathbf{V}$$

$$\begin{array}{rl}{V}_{11}& =\frac{{\omega }_1^2}{2}\text{quantity}(x+\frac{a}{2}{)}^2+\frac{{\omega }_2^2}{2}{y}^2\\ V_{22}& =A{e}^{-(x+b)/D}+\frac{{\omega }_2^2}{2}{y}^2-\Delta \end{array}$$

$$J(\omega )=\frac{2\lambda \gamma \omega }{{\omega }^2+{\gamma }^2},$$

$${\chi }_0(x,y)=\frac{1}{\sqrt{\pi }}\exp -\frac{1}{2}\text{quantity}((x-x_0{)}^2+(y-y_0{)}^2)e^{i{p}_{x}x}$$

$$\mathcal{F}[x^{+}(t),x^{\prime }(t)]=\mathcal{F}[x^{-}(t),x^{\prime }(t)]$$

Full text

I Introduction

Conical intersections associated with electronic degeneracy play critical roles in photochemistry and photophysics [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 3, 11, 12, 13, 14, 15, 16, 17, 18, 19] . One of the intriguing effects due to conical intersections (CIs) is the geometric phase. For time-reversal-invariant electronic Hamiltonian, the adiabatic electronic wavefunction will acquire a phase of $\pi$ upon traversing a closed loop around a CI [20, 21]. As the total molecular wavefunction is single-valued, the multivaluedness of the adiabatic wavefunction is compensated by the same phase in the nuclear wavefunction. The geometric phase is originated from the nontrivial geometry of the adiabatic electronic state space [6, 22] and often manifests as a destructive interference pattern as e.g. a nodal line in the nuclear wave packet dynamics [23]. The geometric phase can have a significant influence on reaction pathways, product distributions, and vibrational spectrum [24, 25, 26]. It should be emphasized that the geometric phase is only meaningful using the adiabatic electronic states (or more precisely, a geometrically nontrival molecular fiber bundle). If a crude adiabatic representation is employed, there will be no geometric phase.

Chemical reactions rarely occur in isolation but are subject to dissipative interactions with environments such as solvents or interfaces. Environmental dissipation can strongly alter the reaction dynamics by quantum dissipation and decoherence [27, 28]. Exploring how the environment influce nonadiabatic dynamics near CIs is important for understanding reactions in complex envrionments. Previous studies have focused on the reaction rates and population dynamics [29, 30, 31, 32, 33], whether the geometric phase-induced interference in the wavepacket dynamics can survive in dissipative environments remains an open question.

In this article, we investigate nuclear wave packet dynamics near CIs under non-Markovian quantum dissipation by an numerically exact modeling of dissipative conical intersection dynamics method. We focus on how quantum decoherence and dissipation influence geometric phase, or more precisely, the geometric phase-induced quantum interference. Theoretically, we employ Feynman path integrals and influence functional theory [34] to dissect the impact of dissipation on geometric-phase-induced wavepacket interference. We combine the recently proposed local diabatic representation (LDR) to describe vibronic couplings [35, 36] and the hierarchical equations of motion (HEOM) method to describe the system-bath coupling [37, 38, 39, 40, 41, 42, 43] for the numerically exact modeling of dissipative conical intersection dynamics.

With a two-dimensional vibronic coupling model [44] coupled to either an electronic or a vibrational bath, we demonstrate that despite dissipation alters the nonadiabatic transitions of the electrons and the density distribution of the nuclears, the destructive interference pattern induced by geometric phase remains robust against strong dissipation.

This paper is organized as follows: In Sec.II, we present the theories of LDR, Feynman-Vernon path integral and HEOM. In Sec. III, we investigate the dynamics of nuclear wavepacket and electron population with the vibrational bath and the electronic bath via the LDR–HEOM method. The following discussion about the effect of quantum dissipation is presented based on the Feynman–Vernon path integral theory. We summarize the paper in Sec. IV.

In this article, we use atomic units with $k_{\text{B}}=\hbar=1$.

II Formalism

II.1 Local diabatic representation of conical intersection dynamics

The most widely used framework for nonadiabatic molecular dynamics is based on the Born-Huang expansion,

[TABLE]

where $\phi_{\alpha}({\mathbf{r}};{\mathbf{R}})$ is the adiabatic electronic states which depends parametrically on the nuclear configuration $\mathbf{R}$, i.e.

[TABLE]

where $H_{\rm BO}({\mathbf{R}})=H_{\text{M}}-\hat{T}_{\text{N}}$ is the electronic Hamiltonian, the full molecular Hamiltonian subtracting the nuclear kinetic energy operator, $V_{\alpha}({\mathbf{R}})$ is the $\alpha$th adiabatic PES, and $\chi_{\alpha}(R,t)$ is the nuclear wave packet evolving on the $\alpha$th adiabatic PES. However, in this Born-Huang representation, the non-Born-Oppenheimer effects including geometric phase and nonadiabatic transitions, critical for conical intersection dynamics, are accounted for by singular terms such as the Berry connection and derivative couplings. To remove the divergences, we recently proposed a local diabatic representation (LDR) [35, 36]. In it, the ansatz for the full molecular wavefunction is given by

[TABLE]

where $\phi_{\alpha}(\mathbf{R}_{\mathbf{n}})$ is the $\alpha$th adiabatic electronic eigenstate of the electronic Born-Oppenheimer Hamiltonian at the nuclear geometry $\mathbf{R}_{\mathbf{n}}$ with energy $V_{\alpha}(\mathbf{R}_{\mathbf{n}})$, $\chi_{\mathbf{n}}({\mathbf{R}})$ is the orthonormal discrete variable representation nuclear basis for the nuclear wavefunction, localized at ${\mathbf{R}}_{\mathbf{n}}$.

Inserting Eq. (II.1) into the molecular time-dependent Schrödinger equation $i\partial{\Psi(\mathbf{r},\mathbf{R},t)}/\partial t=H\Psi(\mathbf{r},\mathbf{R},t)$ with the molecular Hamiltonian $H=\hat{T}_{\text{N}}+H_{\text{BO}}(\mathbf{r};\mathbf{R}),$ and left multiply $\bra{\mathbf{m}\beta}$ yields the equation of motion for the expansion coefficients

[TABLE]

Here $T_{\mathbf{m}\mathbf{n}}=\langle{\chi_{\mathbf{m}}|\hat{T}_{\text{N}}|\chi_{\mathbf{n}}}\rangle_{\mathbf{R}}$ is the kinetic energy operator matrix elements and the electronic overlap matrix

[TABLE]

where $\langle{\cdots}\rangle_{\mathbf{r}}$ ($\langle{\cdots}\rangle_{\mathbf{R}}$) denotes the integration over electronic (nuclear) degrees of freedom. In deriving Eq. (4), we have made use of $H_{\text{BO}}(\mathbf{R})|{{\mathbf{n}}\alpha}\rangle\approx V_{\alpha}(\mathbf{R}_{\mathbf{n}})|{{\mathbf{n}}\alpha}\rangle$ as the nuclear state $|{\mathbf{n}}\rangle$ is an eigenstate of all position operators. Although we have employed the adiabatic electronic states, Eq. (4) does not contain any singularities because the nuclear kinetic energy operator does not operate on the electronic states. Therefore, in contrast to wavepacket dynamics in the Born-Huang representation, the singularity of the derivative couplings at the CI will not affect the simulations based on the LDR.

The geometrical information of the vibronic Hilbert space is encoded in the global overlap matrix $A_{{\mathbf{m}}\beta,{\mathbf{n}}\alpha}$. Specifically, the geometric phase around the CI can be obtained from the product of overlap matrices along a closed loop

[TABLE]

If the loop encircles a CI, $\theta_{\alpha}=\pi$. For time-reversal symmetric Hamiltonian, the geometric phase $\theta_{\alpha}\in\set{0,\pi}$ is topological as it is invariant under any local changes of the loop.

In the LDR, the time-evolution operator is given by

[TABLE]

where $\xi(t)=\set{\xi_{0},\xi_{1},\cdots,\xi_{L}}$ represents all discrete paths in the joint electronic-nuclear space satisfying the boundary conditions $\xi_{0}=(\mathbf{R}_{0},\alpha_{0}),\xi_{L}=(\mathbf{R}_{f},\alpha_{f})$. Each $\xi_{j}=(\mathbf{R}_{j},\alpha_{j})$ labels both the nuclear configuration and electronic state at time $t_{j}=j\Delta t$. Here the action

[TABLE]

with $K_{{\mathbf{R}}_{j+1},{\mathbf{R}}_{j}}=\ln(\langle{\mathbf{R}}_{j+1}|e^{-i\hat{T}_{\rm N}\Delta t}|{\mathbf{R}}_{j}\rangle)/\Delta t$ and the geometric factor

[TABLE]

accounts for the geometric phase effects and nonadiabatic transitions.

The path integral-like picture is particularly convenient to understand how geometric phase can induce destructive nuclear quantum interference around a conical intersection. Consider two adiabatic paths $\Gamma_{-}$ and $\Gamma_{+}$ surrounding a CI. When the dynamical actions of two paths are equal ${S}[\Gamma_{-}]={S}[\Gamma_{+}]=S$, the sum of the transition amplitudes is

[TABLE]

where $\Gamma$ refers to the closed loop traversing along $\Gamma_{+}$ first and then along $\Gamma_{-}$ backwards.

II.2 Influence functional for dissipative conical intersection dynamics

We employ the influence functional to develop an intuitive picture of how a dissipative environment can influence the conical intersection dynamics. We consider a bosonic Gaussian bath, which is widely used for modeling the solvent environments. It contains noninteracting bosons and have linear coupling with the system. The total Hamiltonian is $H=H_{\text{M}}+H_{\text{I}}+H_{\text{B}},$ where $H_{\text{B}}$ is the bath Hamiltonian. The system-bath coupling is assumed to be

[TABLE]

where $Q$ is an arbitrary Hermitian operator acting on either nuclear or electronic space, $X$ is a collective bath coordinate, which is a summation of microscopic degrees of freedom. The bath influence on the system dynamics is completely characterized by the bath correlation function

[TABLE]

with $X(t)=e^{-iH_{\mbox{\tiny B}}t}Xe^{iH_{\mbox{\tiny B}}t}$ and $\langle\cdots\rangle_{\mbox{\tiny B}}={\rm tr}_{\mbox{\tiny B}}\{\cdots\rho^{\rm eq}_{\mbox{\tiny B}}\}$ for the canonical state $\rho^{\rm eq}_{\text{B}}=e^{-\beta H_{\text{B}}}/Z$. The bath correlation function can be obtained from the spectral density $J(\omega)$ via the fluctuation-dissipation theorem,

[TABLE]

with inverse temperature $\beta=1/T$. Here the spectral density, reads

[TABLE]

For Gaussian bath, $J(\omega)$ is a coupling strength-weighted density of states, independent with the state of bath.

In the LDR the vibronic density operator is represented by

[TABLE]

To lighten the notation, we introduce the composite index $\xi=\set{\mathbf{n},\alpha}$. The evolution of reduced density operator can be represented as

[TABLE]

Here the propagator $\mathcal{G}(\xi,\xi^{\prime},t;\xi_{0},\xi_{0}^{\prime},t_{0})$ can be obtained from the path integral.

[TABLE]

$Q[\xi_{j}]$ is the trajectory of $Q$, depending on the composite path $\xi_{j}$. $\mathcal{F}[Q(\xi(t)),Q(\xi^{\prime}(t))]$ is the influence functional which represents the bath’s influence on the evolution of the molecular system [34]. For Gaussian bath, the influence functional is given by

[TABLE]

Although Markovian quantum master equations such as Lindblad and Redfield equations are widely used for open quantum dynamics [45], the exact reduced system dynamics is non-Markovian, meaning that the time-evolution is not only decided by the current state $\rho(t)$, but also by its history $\rho(t^{\prime}),t^{\prime}<t$. The memory length is characterized by the bath correlation function $C(t)$.

II.3 Hierarchical equations of motion

We employ the hierarchical equations of motion (HEOM), an exact and non-perturbative method for the simulation of non-Markovian open quantum dynamics [37, 38, 39, 40, 41, 42, 43, 46, 47]. In HEOM, the bath correlation function is decomposed by a sum of exponentials

[TABLE]

with the time reversal symmetry of the correlation function $C(-t)=C^{*}(t)$,

[TABLE]

Such expansion can be achieved by a sum–over–poles expression for the Fourier integrand on the right–hand–side of Eq. (13), followed by the Cauchy’s contour integration [48, 49, 50, 51, 52] , or using the time–domain fitting decomposition scheme [53, 54]. The second equality of Eq. (20) is due to the fact that the exponents in Eq. (19) must be either real or complex conjugate paired [53], and thus we may determine $\bar{k}$ in the index set $\{k=1,2,...,K\}$ by the pairwise equality $\gamma_{\bar{k}}=\gamma_{k}^{\ast}$ [55]. Via above exponential decomposition we can obtain the HEOM from the Feynman-Vernon influence functional:

[TABLE]

Here we denote ${\bm{\ell}}=[\ell_{1},\cdots,\ell_{K}]$ as the levels of the hierarchy and $L=\sum_{k}\ell_{k}$ with $\ell_{k}=0,1,2,\cdots$. ${\bm{\ell}}^{\pm}_{k}$ means the $[L_{1},\cdots,L_{k}\pm 1,\cdots,L_{K}]$. The time-dependence of the density operators is suppressed in section II.3. The zeroth-tier $\rho^{(L)}_{\mathbf{L}}$ is the molecular reduced density matrix $\rho^{(0)}_{\mathbf{0}}={\rho}$. The higher-order $\{\rho^{(L)}_{\mathbf{L}}\}$ are auxiliary density operators (ADOs), which contain the information of the history of evolution and the system-bath correlation [56, 57, 47, 58]. Each $\rho^{(L)}_{\mathbf{L}}$ connects with higher-tier ADOs $\rho^{(L+1)}_{{\mathbf{L}}_{k}^{+}}$ and lower-tier ADOs $\rho^{(L-1)}_{{\mathbf{L}}_{k}^{-}}$, forming the hierarchical structure. As $\{\ell_{k}\}$ increase, the damping terms $\{-\sum_{k}\ell_{k}\gamma_{k}\rho^{(L)}_{\bm{\ell}}\}$ will suppress the high-order ADOs, which allows us to truncate the HEOM at an enough high tier of the hierarchy [59, 60].

III Results and Discussion

We apply the LDR-HEOM method to simulate the dissipative conical intersection dynamics of

a two-state, two-dimensional vibronic coupling model [44], resembling the photodissociation of phenol, with the Hamiltonian given by

[TABLE]

where $\mathbf{I}$ is the identity matrix in the electronic space, and the nuclear kinetic energy operator $\hat{T}_{\text{N}}=-\frac{1}{2}({{\partial}^{2}/{\partial x^{2}}+{\partial}^{2}/{\partial y^{2}}})$. The diabatic potential energy matrix $\mathbf{V}$ consists of diabatic potential energy surfaces [44]

[TABLE]

and the diabatic coupling $V_{12}=V_{21}=cye^{-(x-x_{\text{CI}})^{2}/{2\sigma^{2}_{x}}-y^{2}/2\sigma^{2}_{y}}$. Here, $x$ resembles the stretching of the O-H bond and $y$ represents the coupling mode. The diabatic coupling is linear around the CI and damped by a Gaussian function away from it. The model parameters are (in a.u.) $\omega_{1}=\omega_{2}=1$, $a=4$, $b=-11$, $c=2$, $A=5$, $\Delta=12$, $x_{\text{CI}}=0$, $D=10$, $\sigma_{x}=1.699$, and $\sigma_{y}=0.849$.

The adiabatic potential energy surfaces (Fig. 3) ${\rm S}_{1}$ and ${\rm S}_{2}$, obtained by diagonalizing the diabatic potential energy matrix, show an energetically inaccessible conical intersection flanked by two energetically lower saddle points. The CI is located at $(x,y)=(0.275,0)$ with energy $E_{\text{CI}}=2.586$. The energy of the two equivalent saddle points is 1.854, forming a potential barrier along the tuning mode.

We consider two types of dissipation, as illustrated in Fig. 2: (i) a vibrational bath coupled to the reaction coordinate, $Q=x$; (ii) an electronic bath introducing energy gap fluctuations, $Q=\Pi_{1}=\sum_{\mathbf{n}}|{\mathbf{R}_{\mathbf{n}}},\phi_{1}({\mathbf{r}};{\mathbf{R}_{\mathbf{n}}})\rangle\langle\phi_{1}({\mathbf{r}};{\mathbf{R}_{\mathbf{n}}}),{\mathbf{R}_{\mathbf{n}}}|$ with $\Pi_{1}$ being the projection operator of the lower adiabatic electronic state. We employ the Drude bath spectral density

[TABLE]

where $\lambda$ is the coupling strength, $\gamma$ is the decay rate. For the electronic bath, the solvent will reorganize when the electronic state changes character. The molecular Hamiltonian is modified by the reorganization of solvent $H_{\mbox{\tiny S}}\rightarrow H_{\mbox{\tiny S}}+\lambda\sum_{\mathbf{n}}|{\mathbf{R}_{\mathbf{n}}},\phi_{1}({\mathbf{R}_{\mathbf{n}}})\rangle\langle\phi_{1}({\mathbf{R}_{\mathbf{n}}}),{\mathbf{R}_{\mathbf{n}}}|$ for the Drude spectrum, reducing the energy gap between ${\rm S}_{1}$ and ${\rm S}_{2}$. (Details in Appendix A)

The initial state is set as a Gaussian wavepacket,

[TABLE]

centered at $(x_{0},y_{0})=(-4,0)$ on the lower state ${\rm S}_{1}$ with initial momentum $p_{x}=2$. Fig. 4 depicts the electronic population and nuclear wave packet dynamics $p_{\alpha}(\mathbf{R},t)=\rho_{\alpha\alpha}(\mathbf{R},\mathbf{R},t)$ with vibrational relaxation at different system-bath coupling strength $\lambda$. Without the bath, as the nuclear wavepacket reaches the CI on ${\rm S}_{1}$, there is significant nonadiabatic transitions and a nodel line along $y=0$ in the probability distribution, which is a hallmark of the geometric phase. With the vibrational bath, the diffusion of the wavepacket along $x$-direction is enhanced, and the nonadiabatic transitions is weakened. Fig. 5 depicts the electron population and wavepacket dynamics with vibrational bath coupling at different temperature $T$. As the temperature increase, we can see the diffusion of the wavepacket along $x$-direction, and transition of electron is enhanced. However, the nodal line in the interference pattern induced by the geometric phase remains intact.

Fig. 6 depicts the electron population and nuclear wavepacket dynamics with an electronic bath at different coupling strength $\lambda$. When electronic bath coupling $\lambda$ increasing, the growth of the excited electron population on ${\rm S}_{2}$ is suppressed. The electron population gets redistributed due to the reorganization of the solvent which modifies the gap between PESs. These changes can also be observed when temperature $T$ increases, as Fig. 7 shows. However, as coupling strength $\lambda$ and temperature $T$ increase, the characteristic interference pattern of geometric phase not only survives, but is further enhanced.

We choose an asymmetrical initial nuclear wavepacket by shifting the center of the Gaussian wavepacket to $(-4,0.5)$. The interference pattern is still clearly visible despite strong dissipation for both electronic and vibrational baths (Fig. 8). We have used a stronger molecule-solvent interaction than typical conditions. This indicates that the robustness of the interference pattern is not due to symmetry.

The robustness of the geometric phase-induced quantum interference can be understood using the path integral picture in Sec. II.1. Consider the two-dimensional system with $y$-reflection symmetry as in Fig. 1. For APES $V(x,y)$ with $V(x,y)=V(x,-y)$, when adiabatic path pair $\Gamma_{+}$ and $\Gamma_{-}$ surrounding CI and connecting ${\mathbf{R}}_{0}$ and ${\mathbf{R}}_{f}$ are symmetric about $y=0$, their dynamical actions are equal ${S}[\Gamma_{+}(t)]={S}[\Gamma_{-}(t)]$. With the phase difference between $\Gamma_{+}$ and $\Gamma_{-}$, $W[\Gamma_{+}]+W[\Gamma_{-}]=W[\Gamma_{-}](1+e^{i\pi})=0$, as this applies to all pairs of paths, there will always be a destructive interference leading to a nodal line in the nuclear distribution along $y=0$ even in the presence of dissipation.

For the vibrational bath $Q=x$, the influence functionals of the path pair $\Gamma_{+}$ and $\Gamma_{-}$ are equal

[TABLE]

the propagator connecting ${\mathbf{R}}_{0}$ and ${\mathbf{R}}_{f}$ can be presented as a sum of all the path pairs $\{\Gamma_{+},\Gamma_{-}\}$

[TABLE]

As the influence functional of the two paths are equal, it can be factored out such that the transition amplitudes interfere destructively.

The same argument can be applied to the electronic bath $Q=\Pi_{1}$, as the influence functional remain equal for the pair of paths with reflection symmetry. Therefore, we have demonstrated that, both the vibrational and electronic bath do not eliminate the destructive interference of different path pairs $\{\Gamma_{-},\Gamma_{+}\}$. Thus as the dissipation strengthening, the geometric phase-induced quantum interference can remain robust. Although this analysis is only rigorously valid in the close vicinity of a conical intersection, it provides useful insights to understand the computational results.

IV Conclusion

Through numerically exact dissipative conical intersection dynamics modeling, we have demonstrated that quantum decoherence does not eliminate the quantum interference originated from the geometric phase. Although the population dynamics can be influenced by the bath, the geometric phase effects is highly robust against both vibrational relaxation and electronic dephasing. Further simulations with asymmetric initial conditions show that the interference pattern surviving from the dissipation is protected by the topology of the CI, not only by the space reflex symmetry.

Our results suggest that geometric phase effects, neglected in most mixed quantum-classical methods may be important even in condensed phase environments.

V Acknowledgement

This work is supported by the National Natural Science Foundation of China (Grant Nos. 22473090 and 92356310). We appreciate the help and advice from Y. Su, Y. Wang, Y.J. Xie, L.Z. Ye and X.T. Zhu.

Appendix A Details about the electronic bath coupling

The total Hamiltonian of bath-electronic state coupling reads

[TABLE]

Here $\{x_{j}\}$ represent the microscopic degrees of freedoms of the bosonic bath and $\{d_{j}({\mathbf{R}_{\mathbf{n}}})\}$ represent the displacement of solvent.

For simplification, we make the approximation that the displacement of solvent $d_{j}({\mathbf{R}_{\mathbf{n}}})$ is independent of the nuclear configuration $d_{j}({\mathbf{R}_{\mathbf{n}}})\approx d_{j}$. Under the approximation, the terms in Eq. (A) can be rearranged as

[TABLE]

[TABLE]

and

[TABLE]

for Drude spectrum.

Thus the total Hamiltonian of bath-electronic state coupling can be decomposed as

[TABLE]