State Dependent Ring Polymer Molecular Dynamics for Investigating Excited Nonadiabatic Dynamics

TL;DR

This paper rigorously derives and tests the non-adiabatic ring polymer molecular dynamics (NRPMD) method, demonstrating its accuracy in simulating excited state non-adiabatic dynamics and its ability to address zero-point energy leakage.

Contribution

The paper provides a rigorous derivation of the NRPMD Hamiltonian using MMST mapping and path-integral nuclear description, validating its effectiveness for excited state dynamics.

Findings

NRPMD alleviates zero-point energy leakage.

NRPMD accurately simulates excited state non-adiabatic dynamics.

Theoretical foundation established for state-dependent RPMD methods.

Abstract

Recently proposed non-adiabatic ring polymer molecular dynamics (NRPMD) approach has shown to provide accurate quantum dynamics by incorporating explicit electronic state descriptions and nuclear quantizations. Here, we present a rigorous derivation of the NRPMD Hamiltonian and investigate its performance on simulating excited state non-adiabatic dynamics. Our derivation is based on the Meyer-Miller-Stock-Thoss (MMST) mapping representation for electronic states and the ring-polymer path-integral description for nuclei, resulting in the same Hamiltonian proposed in the original NRPMD approach. In addition, we investigate the accuracy of using NRPMD to simulate photoinduced non-adiabatic dynamics in simple model systems. These model calculations suggest that NRPMD can alleviate the zero-point energy leakage problem that is commonly encountered in the classical Wigner dynamics, and provide accurate excited states non-adiabatic dynamics. This work provides a solid theoretical foundation of the promising NRPMD Hamiltonian and demonstrates the possibility of using state-dependent RPMD approach to accurately simulate electronic non-adiabatic dynamics while explicitly quantize nuclei.

Figures

Tables

Table 1. Table 1: Parameters for Model IA-IC (in atomic units).

	Model IA			Model IB			Model IC
$i$	1	2	3	1	2	3	1	2	3
$D_{ii}$	0.02	0.02	0.003	0.02	0.01	0.003	0.003	0.004	0.003
${\alpha }_{ii}$	0.4	0.65	0.65	0.65	0.4	0.65	0.65	0.6	0.65
$R_{ii}$	4.0	4.5	6.0	4.5	4.0	4.4	5.0	4.0	6.0
$c_{ii}$	0.02	0.0	0.02	0.0	0.01	0.02	0.0	0.01	0.006
$ij$	12	13	23	12	13	23	12	13	23
$A_{ij}$	0.005	0.005	0.0	0.005	0.005	0.0	0.002	0.0	0.002
${\alpha }_{ij}$	32.0	32.0	0.0	32.0	32.0	0.0	16.0	0.0	16.0
$R_{ij}$	3.40	4.97	0.0	3.66	3.34	0.0	3.40	0.0	4.8

Equations

$$\hat{H}=\hat{T}+{\hat{V}}_0+{\hat{H}}_e=\frac{{\hat{\mathbf{P}}}^2}{2\mathbf{M}}+V_0(\hat{\mathbf{R}})+\sum _{n,m=1}^{\mathcal{K}}{V}_{nm}(\hat{\mathbf{R}})|n⟩⟨m|,$$

$$\mathcal{Z}=\underset{N\to \infty }{\lim }\int d\{{\mathbf{P}}_{\alpha }\}d\{{\mathbf{R}}_{\alpha }\}{e}^{-{\beta }_N{H}_N^0}{\text{Tr}}_{\mathrm{e}}\left[\prod _{\alpha =1}^N{e}^{-{\beta }_N{\hat{H}}_e({\mathbf{R}}_{\alpha })}\right],$$

$$H_N^0=\sum _{\alpha =1}^N\frac{{{\mathbf{P}}_{\alpha }}^2}{2\mathbf{M}}+\frac{\mathbf{M}}{2{\beta }_N^2{\hslash }^2}({\mathbf{R}}_{\alpha }-{\mathbf{R}}_{\alpha -1}{)}^2+V_0({\mathbf{R}}_{\alpha }),$$

$$|n⟩\to |0_1...1_n...0_{\mathcal{K}}⟩={\hat{a}}_n^{†}|0_1...0_n...0_{\mathcal{K}}⟩.$$

$$|n⟩⟨m|\to {\hat{a}}_n^{†}{\hat{a}}_m,$$

$${\hat{H}}_{\mathrm{e}}({\mathbf{R}}_{\alpha })\equiv \sum _{n,m}{V}_{nm}({\mathbf{R}}_{\alpha })|n⟩⟨m|\to \sum _{n,m}{V}_{nm}({\mathbf{R}}_{\alpha }){\hat{a}}_n^{†}{\hat{a}}_m.$$

$$\displaystyle\text{Tr}_{\mathrm{e}}\big{[}\prod_{\alpha=1}^{N}e^{-\beta_{N}\hat{H}_{e}({\bf R}_{\alpha})}\big{]}$$

$$\displaystyle=\frac{1}{(2\pi\hbar)^{\mathcal{K}}}\int d{\bf q}_{1}d{\bf p}_{1}\big{[}e^{-\beta_{N}\hat{H}_{e}({\bf R}_{1})}\prod_{\alpha=2}^{N}e^{-\beta_{N}\hat{H}_{e}({\bf R}_{\alpha})}\mathcal{P}\big{]}_{\mathrm{w}_{1}}.$$

$$[O{]}_{{\mathrm{w}}_{\alpha }}=\int d{\mathbf{\Delta }}_{\alpha }{e}^{i{\mathbf{p}}_{\alpha }^{\mathrm{T}}{\mathbf{\Delta }}_{\alpha }/\hslash }⟨{\mathbf{q}}_{\alpha }-\frac{{\mathbf{\Delta }}_{\alpha }}{2}\left|\hat{\mathrm{O}}\right|{\mathbf{q}}_{\alpha }+\frac{{\mathbf{\Delta }}_{\alpha }}{2}⟩.$$

$$\int d\mathbf{q}d\mathbf{p}[AB{]}_{\mathrm{w}}=\int d\mathbf{q}d\mathbf{p}[A{]}_{\mathrm{w}}[B{]}_{\mathrm{w}},$$

$$\displaystyle\text{Tr}_{\mathrm{e}}\bigg{[}\prod_{\alpha=1}^{N}e^{-\beta_{N}\hat{H}_{\mathrm{e}}({\bf R}_{\alpha})}\bigg{]}$$

$$\displaystyle=\frac{1}{(2\pi\hbar)^{\mathcal{K}}}\int d{\bf q}_{1}d{\bf p}_{1}\bigg{[}e^{-\beta_{N}\hat{H}_{\mathrm{e}}({\bf R}_{1})}\prod_{\alpha=2}^{N}e^{-\beta_{N}\hat{H}_{\mathrm{e}}({\bf R}_{\alpha})}\mathcal{P}\bigg{]}_{\mathrm{w}_{1}}$$

$$\displaystyle=\frac{1}{(2\pi\hbar)^{\mathcal{K}}}\int d{\bf q}_{1}d{\bf p}_{1}\bigg{[}e^{-\beta_{N}\hat{H}_{\mathrm{e}}({\bf R}_{1})}\bigg{]}_{\mathrm{w}_{1}}\bigg{[}\prod_{\alpha=2}^{N}e^{-\beta_{N}\hat{H}_{e}({\bf R}_{\alpha})}\mathcal{P}\bigg{]}_{\mathrm{w}_{1}}$$

$$\displaystyle=\frac{1}{(2\pi\hbar)^{\mathcal{K}}}\int d{\bf q}_{1}d{\bf p}_{1}\bigg{[}1-\beta_{N}\hat{H}_{\mathrm{e}}({\bf R}_{1})\bigg{]}_{\mathrm{w}_{1}}\bigg{[}\prod_{\alpha=2}^{N}e^{-\beta_{N}\hat{H}_{e}({\bf R}_{\alpha})}\mathcal{P}\bigg{]}_{\mathrm{w}_{1}}$$

$$=\frac{1}{(2\pi \hslash {)}^{\mathcal{K}}}\int d{\mathbf{q}}_{1}d{\mathbf{p}}_1{e}^{-{\beta }_N[{\hat{H}}_e({\mathbf{R}}_1){]}_{{\mathrm{w}}_1}}{\left[\prod _{\alpha =2}^N{e}^{-{\beta }_N{\hat{H}}_{\mathrm{e}}({\mathbf{R}}_{\alpha })}\mathcal{P}\right]}_{{\mathrm{w}}_1}$$

$$[{\hat{H}}_{\mathrm{e}}({\mathbf{R}}_1){]}_{{\mathrm{w}}_1}$$

$$[{\hat{H}}_e({\mathbf{R}}_{\alpha }){]}_{{\mathrm{w}}_{\alpha }}$$

$$\displaystyle=\frac{1}{2\hbar}\sum_{nm}V_{nm}({\bf R}_{\alpha})\bigg{(}[{\bf q}_{\alpha}]_{n}[{\bf q}_{\alpha}]_{m}+[{\bf p}_{\alpha}]_{n}[{\bf p}_{\alpha}]_{m}-\delta_{nm}\hbar\bigg{)}.$$

$${\text{Tr}}_{\mathrm{e}}\left[\prod _{\alpha =1}^N{e}^{-{\beta }_N{\hat{H}}_{\mathrm{e}}({\mathbf{R}}_{\alpha })}\right]$$

$$=\frac{1}{(2\pi \hslash {)}^{N\mathcal{K}}}\int d\{{\mathbf{q}}_{\alpha }\}d\{{\mathbf{p}}_{\alpha }\}{e}^{-{\beta }_N{\sum }_{\alpha }[{\hat{H}}_e({\mathbf{R}}_{\alpha }){]}_{{\mathrm{w}}_{\alpha }}}$$

$$\times \int d\{{\mathbf{\Delta }}_{\alpha }\}\prod _{\alpha =1}^N{e}^{i{\mathbf{p}}_{\alpha }^{\mathrm{T}}{\mathbf{\Delta }}_{\alpha }/\hslash }⟨{\mathbf{q}}_{\alpha }-\frac{{\mathbf{\Delta }}_{\alpha }}{2}|\mathcal{P}|{\mathbf{q}}_{\alpha +1}+\frac{{\mathbf{\Delta }}_{\alpha +1}}{2}⟩.$$

$${\mathcal{Z}}_N=\underset{N\to \infty }{\lim }\frac{1}{(2\pi \hslash {)}^{N\mathcal{K}}}\int d\{{\mathbf{P}}_{\alpha }\}d\{{\mathbf{R}}_{\alpha }\}d\{{\mathbf{p}}_{\alpha }\}d\{{\mathbf{q}}_{\alpha }\}\mathbf{\Gamma }{e}^{-{\beta }_N{H}_N},$$

$$\mathbf{\Gamma }$$

$$\displaystyle{H}_{N}=\sum_{\alpha=1}^{N}\bigg{[}{{{\bf P_{\alpha}}^{2}}\over{2M}}+V_{0}({\bf R_{\alpha}})+{{\bf M}\over{2\beta^{2}_{N}\hbar^{2}}}({\bf R}_{\alpha}-{\bf R}_{\alpha-1})^{2}$$

$$\displaystyle+\frac{1}{2\hbar}\sum_{nm}V_{nm}({\bf R}_{\alpha})\big{(}[{\bf q}_{\alpha}]_{n}[{\bf q}_{\alpha}]_{m}+[{\bf p}_{\alpha}]_{n}[{\bf p}_{\alpha}]_{m}-\delta_{nm}\hbar\big{)}\bigg{]}.$$

$$\displaystyle\dot{\bf P}_{\alpha}=-{{\bf M}\over{\beta^{2}_{N}\hbar^{2}}}\big{(}2{\bf R}_{\alpha}-{\bf R}_{\alpha+1}-{\bf R}_{\alpha-1}\big{)}-\nabla_{{\bf R}_{\alpha}}V_{0}({\bf R}_{\alpha})$$

$$\displaystyle~{}~{}-\frac{1}{2\hbar}\sum_{nm}\nabla_{{\bf R}_{\alpha}}V_{nm}({\bf R}_{\alpha})\big{(}[{\bf q}_{\alpha}]_{n}[{\bf q}_{\alpha}]_{m}+[{\bf p}_{\alpha}]_{n}[{\bf p}_{\alpha}]_{m}-\delta_{nm}\hbar\big{)},$$

$$[{\dot{\mathbf{q}}}_{\alpha }{]}_n=\frac{\partial H_N}{\partial [{\mathbf{p}}_{\alpha }{]}_n}=\frac{1}{\hslash }\sum _m{V}_{nm}({\mathbf{R}}_{\alpha })[{\mathbf{p}}_{\alpha }{]}_m$$

$$[{\dot{\mathbf{p}}}_{\alpha }{]}_n=-\frac{\partial H_N}{\partial [{\mathbf{q}}_{\alpha }{]}_n}=-\frac{1}{\hslash }\sum _m{V}_{nm}({\mathbf{R}}_{\alpha })[{\mathbf{q}}_{\alpha }{]}_m.$$

$${\rho }_{jj}(t)={\mathrm{Tr}}_{\mathrm{n}}[{\hat{\rho }}_0{e}^{i\hat{H}t/\hslash }{\hat{\mathcal{P}}}_j{e}^{-i\hat{H}t/\hslash }],$$

$${\rho }_{jj}(t)\approx \int d\mathbf{\tau }{\mathcal{P}}^0(\{{\mathbf{q}}_{\alpha }(0),{\mathbf{p}}_{\alpha }(0)\}){\rho }_{\mathrm{rp}}(\{{\mathbf{R}}_{\alpha }(0),{\mathbf{P}}_{\alpha }(0)\})\cdot {\bar{\mathcal{P}}}_j(t).$$

$${\bar{\mathcal{P}}}_j=\frac{1}{N}\sum _{\alpha }{\mathcal{P}}_j(\alpha )=\frac{1}{N}\sum _{\alpha =1}^N\frac{1}{2}([{\mathbf{q}}_{\alpha }{]}_j^2+[{\mathbf{p}}_{\alpha }{]}_j^2-1),$$

$${\mathcal{P}}^0(\{{\mathbf{q}}_{\alpha }(0),{\mathbf{p}}_{\alpha }(0)\})=\prod _{\alpha =1}^N\prod _{j=1}^{\mathcal{K}}\delta ({\mathcal{P}}_j(\alpha )-{\rho }_{jj}(0)).$$

$${\mathcal{P}}_j(\alpha )=\frac{1}{2}([{\mathbf{q}}_{\alpha }{]}_j^2+[{\mathbf{p}}_{\alpha }{]}_j^2-1)={\delta }_{ij}.$$

$$V_{ii}$$

$$V_{ij}$$

Full text

I Introduction

Accurately and efficiently simulating electronic non-adiabatic transitions and nuclear quantum effects remain one of the central challenges in theoretical chemistry.J.C.Tully (2012); Althorpe et al. (2016) Directly performing exact quantum dynamics simulations in large systems are computationally demanding, despite exciting recent progress.Meyer, Manthe, and Cederbaum (1990); Beck et al. (2000); Wang and Thoss (2003); Eisenbrandt et al. (2018); Richings and Habershon (2018); Makhov et al. (2014); Curchod and Martinez (2018); Ishizaki and Tanimura (2005); Tanimura (2006); Chen et al. (2009); Strumpfer and Schulten (2012); Makri and Makarov (1995a, b); Strathearn et al. (2017); Walters and Makri (2016, 2015); Greene and Batista (2017); Ren, Shuai, and Chan (2018) It is thus ideal to develop trajectory-based approximate methods that can accurately describe the electronic and nuclear quantum effects. To this end, a large number of these approaches are developed, including the popular trajectory surface-hopping method (mixed quantum-classical approach),Tully (1990); Subotnik et al. (2016); Wang, Akimov, and Prezhdo (2016); Crespo-Otero and Barbatti (2018) the semi-classical path-integral approaches,Miller (2001); W.H.Miller (2009); Sun, Wang, and Miller (1998); Shi and Geva (2004); Bonella and Coker (2005); Makri (2011); Huo and Coker (2011); Lee, Huo, and Coker (2016) the mixed quantum-classical Liouville equation, Kernan, Ciccotti, and Kapral (2008); Kim, Nassimi, and Kapral (2008); Hsieh and Kapral (2012, 2013); Kapral (2015) the symmetrical quasi-classical (SQC) approach,Miller and Cotton (2016); Tao (2016); Kananenka et al. (2018); Sandoval, Mandal, and Huo (2018); Cotton and Miller (2019) and the generalized quantum master equation (GQME).Zhang et al. (2006); Pfalzgraff, Kelly, and Markland (2015); Kelly, Brackbill, and Markland (2015); Kelly et al. (2016); Mulvihill et al. (2019) Despite that they are promising to accurately describe electronic non-adiabatic dynamics, these approaches generally do not preserve quantum Boltzmann distribution (QBD)Parandekar and Tully (2006); Schmidt, Parandekar, and Tully (2008) or zero point energy (ZPE) associated with the nuclear degrees of freedom (DOF), and often suffer from numerical issues such as ZPE leakage, Habershon and Manolopoulos (2009); Muller and Stock (1999) although significant improvements are accomplished through the SQCMiller and Cotton (2016) and GQME approaches.Kelly, Brackbill, and Markland (2015); Kelly et al. (2016); Mulvihill et al. (2019)

Imaginary-time path-integral approaches,Berne and Thirumalai (1986); Ceperley (1995); Chandler and Wolynes (1981) such as the centroid molecular dynamics (CMD)Cao and Voth (1994); Jang and Voth (1999) and the ring-polymer molecular dynamics (RPMD),Habershon et al. (2013); Craig and Manolopoulos (2004) resemble classical MD in an extended phase space and provide a convenient way to compute approximate quantum time-correlation functions.Habershon et al. (2013) The classical evolution of RPMD preserves its initial quantum distribution captured by the ring-polymer Hamiltonian, and it is free of the zero-point energy leaking problem.Habershon and Manolopoulos (2009); Habershon et al. (2013) Despite its success on describing quantum effects in the condensed phase, RPMD is limited to one-electron non-adiabatic dynamicsMenzeleev, Ananth, and Miller (2010, 2011); Kretchmer and Miller (2013, 2016); Kenion and Ananth (2016) or nuclear quantization,Habershon et al. (2013); N. Boekelheide and Miller (2011); L.Wang and Markland (2014); Marsalek and Markland (2017); Markland and Ceriotti (2018) as well as the lack of real-time electronic coherence effects.Menzeleev, Ananth, and Miller (2010, 2011) Recently emerged state-dependent RPMD approaches, such as non-adiabatic RPMD (NRPMD),Richardson and Thoss (2013); Hele and Ananth (2016); Richardson et al. (2017) mapping variable RPMD (MV-RPMD),Ananth (2013); J.R.Duke and Ananth (2015); Pierre et al. (2017) Kinetically-constrained RPMD (KC-RPMD),Menzeleev, Bell, and Miller (2014); Kretchmer and Miller (2016); Kretchmer et al. (2018) coherent state RPMD (CS-RPMD),S.N.Chowdhury and Huo (2017) and ring-polymer surface hopping (RPSH) Shushkov, Li, and Tully (2012); Shakib and Huo (2017); Tao, Shushkov, and Miller (2018, 2019) are promising to provide accurate non-adiabatic dynamics with an explicit description of electronic states, in addition to the reliable treatment of nuclear quantum effects through ring polymer quantization.

Among these state-dependent RPMD approaches, NRPMDRichardson and Thoss (2013); Richardson et al. (2017) and CS-RPMDS.N.Chowdhury and Huo (2017) have shown to accurately describe both the electronic dynamics and nuclear quantum effects. The NRPMDRichardson and Thoss (2013); Richardson et al. (2017) and CS-RPMDS.N.Chowdhury and Huo (2017) Hamiltonian can be viewed as the generalization of the Meyer-Miller-Stock-Thoss (MMST) mapping HamiltonianMeyer and Miller (1979); Stock and Thoss (1997, 1999) with the ring-polymer description of the nuclei. Both Hamiltonians have many desired properties, such as a clear adiabatic limit (that returns to the original RPMD Hamiltonian) and one bead limit (that returns to the original MMST Hamiltonian), and recovering the correct electronic Rabi oscillations when the electronic states and nuclei are decoupled.Richardson and Thoss (2013); Richardson et al. (2017); S.N.Chowdhury and Huo (2017) Nevertheless, the promising NRPMD Hamiltonian is proposed through a physically motivated but ad hoc fashion.Richardson and Thoss (2013)

In this paper, we provide a rigorous derivation of the NRPMD Hamiltonian, which is based on the MMST mapping formalismMeyer and Miller (1979); Stock and Thoss (1997, 1999) for the electronic DOF and the ring polymer path-integral representationFeynman and Hibbs (1965); Berne and Thirumalai (1986); Ceperley (1995); Chandler and Wolynes (1981) for the nuclear DOF, leading to the same Hamiltonian that has been previously proposed in the NRPMD approach.Richardson and Thoss (2013); Richardson et al. (2017) The NRPMD Hamiltonian and the previously derived CS-RPMD HamiltonianS.N.Chowdhury and Huo (2017) can be viewed as a unified classical theory for electronic states (through the MMST mapping formalism) and nuclei (through the ring polymer quantization).

In addition, we perform numerical simulations to investigate the accuracy of NRPMD Hamiltonian for excited states non-adiabatic dynamics. Despite that RPMD-based approaches are initially developed for investigating quantum dynamics under thermal-equilibrium conditions,Habershon et al. (2013) recent workR.Welsch et al. (2016) based on the Matsubara dynamics frameworkHele et al. (2015); T. J. H. Hele and Althorpe (2015) (which does not subject to the restriction of equilibrium conditionsR.Welsch et al. (2016)) has demonstrated that RPMD yields reliable non-equilibrium time-correlation function. Further, RPMD based approachesMenzeleev, Ananth, and Miller (2010); J.R.Duke and Ananth (2015) have already been used to simulate non-adiabatic dynamics under non-equilibrium initial conditions. These studies inspire us to investigate the performance of NRPMD for simulating excited states non-adiabatic dynamics.

Our numerical results with model calculations suggest that NRPMD can provide an accurate short-time non-adiabatic branching dynamics among many coupled electronic states, as well as the longer time dynamics such as the plateau value or recurrence of the oscillations of the electronic population. Quantizing nuclei through ring polymer instead of using the Wigner distribution alleviates zero-point energy leakage problems that classical Wigner dynamics encountered, leading to more accurate electronic quantum dynamics. These encouraging numerical results, together with our rigorous derivation of the NRPMD Hamiltonian open up new possibilities to accurately simulate non-adiabatic quantum dynamics while explicitly quantize nuclei.

II Theory

In this section, we provide a rigorous derivation of the NRPMD HamiltonianRichardson and Thoss (2013); Richardson et al. (2017) through the imaginary-time formalism. Based upon that, we propose a numerical procedure to compute the electronic population for excited state non-adiabatic dynamics. We start by expressing the total Hamiltonian operator as

[TABLE]

where $\{|n\rangle\}$ is the diabatic basis, $\hat{T}$ is the nuclear kinetic energy operator, $\hat{\bf R}\equiv\{\hat{R}_{1},...,\hat{R}_{\mathcal{F}}\}$ is the nuclear position operator associated with $\mathcal{F}$ different nuclear DOF, with the corresponding conjugate momentum operator $\hat{\bf P}\equiv\{\hat{P}_{1},...,\hat{P}_{\mathcal{F}}$} and the nuclear mass ${\bf M}\equiv\{M_{1},M_{2},...,M_{\mathcal{F}}\}$. In addition, $V_{0}({\bf\hat{R}})$ is the state-independent potential operator, whereas $\hat{H}_{\mathrm{e}}=\sum_{nm}V_{nm}({\bf\hat{R}})|n\rangle\langle m|$ is the state-dependent potential operator (electronic part of the Hamiltonian) with $\mathcal{K}$ total diabatic electronic states.

The canonical partition function of the system is defined as $\mathcal{Z}=\mathrm{Tr_{en}}[e^{-\beta\hat{H}}]$, where $\mathrm{Tr_{en}=Tr_{e}Tr_{n}}$ represents the trace over both the electronic and nuclear DOFs, with $\mathrm{Tr_{e}}[...]\equiv\sum_{n=1}^{\mathcal{K}}\langle n|...|n\rangle$, and $\mathrm{Tr_{n}}[...]\equiv\int(...)d{\bf R}$. Further, $\beta=1/k_{\mathrm{B}}T$ is the reciprocal temperature, and $\hat{H}$ is the total Hamiltonian operator defined in Eqn. 1. The partition function can be exactly evaluated as $\mathcal{Z}=\mathrm{Tr_{en}}\left[\prod_{\alpha=1}^{N}e^{-\beta_{N}\hat{H}}\right]$, with a higher effective temperature defined as $\beta_{N}=\beta/N$. Further splitting the Boltzmann operator by trotter expansion under the infinite bead limit $N\rightarrow\infty$ gives $\mathcal{Z}=\lim_{N\to\infty}\mathrm{Tr_{en}}\left[\prod_{\alpha=1}^{N}e^{-\beta_{N}(\hat{T}+\hat{V}_{0})}e^{-\beta_{N}\hat{H}_{e}}\right]$. Inserting $N$ copies of the resolution of identity $I_{\bf R}=\int{\bf dR_{\alpha}}|\bf{R_{\alpha}}\rangle\langle{\bf R_{\alpha}}|$ and $I_{\bf P}=\int{\bf dP_{\alpha}}|\bf{P_{\alpha}}\rangle\langle{\bf P_{\alpha}}|$, and explicitly performing the trace over the nuclear DOF based on the standard path-integral technique, Feynman and Hibbs (1965); Berne and Thirumalai (1986); Ceperley (1995); Chandler and Wolynes (1981) we have

[TABLE]

with $\int d\{{\bf P_{\alpha}}\}d\{{\bf R_{\alpha}}\}=\int\prod_{\alpha=1}^{N}d{\bf P_{\alpha}}d{\bf R_{\alpha}}$. The trace over electronic DOF can be performed with diabatic states as $\mathrm{Tr}_{\mathrm{e}}[...]=\sum_{n=1}^{\mathcal{K}}\langle n|...|n\rangle$. The state-independent Hamiltonian $H^{0}_{N}$ is expressed as follows

[TABLE]

where $H_{\mathrm{rp}}=H^{0}_{N}-\sum_{\alpha}V_{0}({\bf R_{\alpha}})$ is so-called the free ring-polymer Hamiltonian.Craig and Manolopoulos (2004); M.Ceriotti et al. (2010) The above partition function is a common expression for all state-dependent RPMD approachesMenzeleev, Bell, and Miller (2014); Ananth (2013); Richardson and Thoss (2013); Hele (2011); Duke and Ananth (2016) and PIMC methods.Alexander (2001); Schmidt and Tully (2007); Ananth and Miller (2010); Lu and Zhou (2017) The only difference among these approaches arises from the treatment of the electronic term $\text{Tr}_{\mathrm{e}}[\prod_{\alpha=1}^{N}e^{-\beta_{N}\hat{H}_{e}({\bf R_{\alpha}})}]$. For example, in the mean-field RPMD approach,Hele (2011); Duke and Ananth (2016) the electronic potential is obtained from a weighted average of ring-polymer in different electronic configurations; in the MV-RPMD approach,Ananth (2013) the electronic states are explicitly described with mapping variables in the Wigner representation.

In this work, we adapt the MMST mapping representationMeyer and Miller (1979); Stock and Thoss (1997, 1999) to transform discrete electronic states into continuous variables. Based on this representation, $\mathcal{K}$ diabatic electronic states are mapped onto $\mathcal{K}$ harmonic oscillators’ ground and first excited states through the following relation

[TABLE]

Here, $|0_{1}...1_{n}...0_{\mathcal{K}}\rangle$ is the singly excited oscillator (SEO) state with $(\mathcal{K}-1)$ oscillators in their ground states and the $n_{\mathrm{th}}$ oscillator in its first excited state. Thus, the MMST formalism providesStock and Thoss (1997, 1999) the following mapping relation

[TABLE]

with $\hat{a}^{\dagger}_{n}={1/{\sqrt{2\hbar}}}\left(\hat{q}_{n}-\mathrm{i}\hat{p}_{n}\right)$ and $\hat{a}_{m}={1/{\sqrt{2\hbar}}}\left(\hat{q}_{m}+\mathrm{i}\hat{p}_{m}\right)$ as the creation and annihilation operators of the harmonic oscillator.

With the MMST representation, we express the state-dependent potential operator $\hat{H}_{\mathrm{e}}({\bf R_{\alpha}})$ in Eqn. 2 as follows

[TABLE]

This is commonly referred as the MMST Hamiltonian.Meyer and Miller (1979); Stock and Thoss (1997) Using the above MMST representation, we can derive a partition function that contains the NRPMD Hamiltonian.

II.1 Derivation of the NRPMD Hamiltonian

We begin by replacing the trace over the electronic DOF in Eqn. 2 with the phase space integral of mapping variables in the Wigner representationCase (2008); Hele and Ananth (2016); Ananth (2013); J.R.Duke and Ananth (2015) as follows

[TABLE]

In the above equation, the Wigner-Weyl transformE.Wigner (1932); Hillery et al. (1984); Case (2008) of the $\alpha_{\mathrm{th}}$ bead’s mapping DOF is defined as

[TABLE]

We use the notation ${\bf q}_{\alpha}\equiv\{[{\bf q}_{\alpha}]_{1},...[{\bf q}_{\alpha}]_{n},...[{\bf q}_{\alpha}]_{\mathcal{K}}\}$ and ${\bf p}_{\alpha}\equiv\{[{\bf p}_{\alpha}]_{1},...[{\bf p}_{\alpha}]_{n},...[{\bf p}_{\alpha}]_{\mathcal{K}}\}$ to represent $\mathcal{K}$ mapping variables for $\mathcal{K}$ electronic states associated with the $\alpha_{\mathrm{th}}$ imaginary-time slice. Further, we have inserted an electronic projection operator $\mathcal{P}=\sum_{n}|n\rangle\langle n|$ to restrain the mapping variables within the SEO subspace.Ananth (2013); Ananth and Miller (2010)

Recall that the Wigner transform has the following propertyCase (2008)

[TABLE]

where $[AB]_{\mathrm{w}}$,$[A]_{\mathrm{w}}$, and $[B]_{\mathrm{w}}$ are Wigner transforms of operator $\hat{A}\hat{B}$, $\hat{A}$, and $\hat{B}$, respectively, defined in Eqn. 8. With the above equality, we can rewrite Eqn. 7 as follows

[TABLE]

To derive the last two lines of the above equation, we use the fact that under the limit $\beta_{N}\rightarrow 0$, Boltzmann operator becomes $e^{-\beta_{N}\hat{H}_{e}({\bf R}_{1})}=1-\beta_{N}\hat{H}_{e}({\bf R}_{1})+\mathcal{O}(\beta_{N}^{2})$. Under this limit, $[e^{-\beta_{N}\hat{H}_{e}({\bf R}_{1})}]_{\mathrm{w}_{1}}=[1-\beta_{N}\hat{H}_{e}({\bf R}_{1})]_{\mathrm{w}_{1}}=1-\beta_{N}[\hat{H}_{e}({\bf R}_{1})]_{\mathrm{w}_{1}}$. Expressing this term back to the full exponential factor, we arrived at the last line of Eqn. 10.

To evaluate $[\hat{H}_{\mathrm{e}}({\bf R}_{1})]_{\mathrm{w}_{1}}$, we use the mapping Hamiltonian relation described in Eqn. 6 and the Wigner transform defined in Eqn. 8, leading to the following expression

[TABLE]

To obtain an explicit expression, we use $\hat{a}^{\dagger}_{n}={1/{\sqrt{2\hbar}}}\left(\hat{q}_{n}-\mathrm{i}\hat{p}_{n}\right)$ and $\hat{a}_{m}={1/{\sqrt{2\hbar}}}\left(\hat{q}_{m}+\mathrm{i}\hat{p}_{m}\right)$, and evaluate these Wigner integrals. This derivation is provided in Appendix A, giving the Wigner transformed MMST mapping Hamiltonian of the $\alpha_{\mathrm{th}}$ bead as follows

[TABLE]

To evaluate $\left[\prod_{\alpha=2}^{N}e^{-\beta_{N}\hat{H}_{e}({\bf R}_{\alpha})}\mathcal{P}\right]_{\mathrm{w}_{1}}$ in Eqn 10, we repeat the same procedure used in Eqn. 10 for the rest of the imaginary-time slices. The details of this derivation are provided in Appendix A, and here, we briefly summarize the procedure. First, we insert the resolution of identity $\int d{\bf q}_{2}\mathcal{P}|{\bf q}_{2}\rangle\langle{\bf q}_{2}|$ into $\left[\prod_{\alpha=2}^{N}e^{-\beta_{N}\hat{H}_{e}({\bf R}_{\alpha})}\mathcal{P}\right]_{\mathrm{w}_{1}}$. Second, we rearrange the order of terms in the integral and reexpress it as a trace over ${\bf q}_{2}$. Finally, replacing the trace over ${\bf q}_{2}$ by a Wigner transform, we can use the property in Eqn. 9 again and repeat the same procedure outlined in Eqn. 10 to factorize the total Wigner function as a product of two Wigner functions inside the $\int d{\bf q}_{2}d{\bf p}_{2}$ integral. Repeating the above outlined process for every single bead, we ended up with the following expression of the electronic trace

[TABLE]

Analytically evaluate the last term of the above expressionAnanth (2013) (with details provided in Appendix A), and plugging the result into Eqn. 2, we arrive at the final expression of the partition function

[TABLE]

where ${\bf\Gamma}$ has the following expression

[TABLE]

The derived NRPMD Hamiltonian in Eqn. 14, which is one of the central results in this paper, has the following expression

[TABLE]

The above Hamiltonian has been proposed in the original NRPMD approach,Richardson and Thoss (2013); Richardson et al. (2017) and the electronic part of this Hamiltonian (i.e., the second line of Eqn. 16) has also been rigorously derived through a mapping Liouvillian in the state-dependent generalized Kubo-transformed time-correlation function formalism.Hele and Ananth (2016) Here, we provide a rigorous derivation of the full Hamiltonian.

In the NRPMD approach,Richardson and Thoss (2013); Richardson et al. (2017) classical trajectories are propagated according to the Hamilton’s equation of motion associated with ${H}_{N}$ in Eqn. 16. The motion of the nuclei is governed by $\dot{\bf P}_{\alpha}=-\nabla_{{\bf R}_{\alpha}}H_{N}$ as follows

[TABLE]

whereas the bead-specific mapping variables are propagated based on the following Hamilton’s equation of motion

[TABLE]

The above NRPMD equations of motionRichardson and Thoss (2013); Richardson et al. (2017) has been used to compute approximate Kubo-transformed time-correlation functions under thermal conditions, such as the position auto-correlation function,Richardson and Thoss (2013) population auto-correlation function,Richardson and Thoss (2013) and absorption spectraRichardson et al. (2017) based on the Fourier transformed dipole auto-correlation function.

II.2 Excited States Non-adiabatic Dynamics with NRPMD

With the derived NRPMD Hamiltonian, we propose to investigate the excited states non-adiabatic dynamics. The central quantity that we aim to compute is the following reduced density matrix

[TABLE]

where, ${\rho}_{jj}(t)$ is the time-dependent population of state $|j\rangle$, $\mathcal{\hat{P}}_{j}=|j\rangle\langle j|$ is the projection operator associated with state $|j\rangle$, $\mathrm{Tr}_{\mathrm{n}}$ represents the trace over the nuclear DOF, and the initial density operator for the entire system is $\hat{\rho}_{0}=|i\rangle\langle i|\otimes\hat{\rho}_{\mathrm{n}}$, which is a direct product of the initial electronic state $|i\rangle$ and the initial nuclear density operator $\hat{\rho}_{\mathrm{n}}$.

Despite that the RPMD-based approaches are originally developed for investigating quantum dynamics under thermal-equilibrium conditions,Habershon et al. (2013) recent workR.Welsch et al. (2016) has shown that RPMD yields the exact non-equilibrium time-correlation function under high temperatures, short-time, and harmonic potential limits through rigorous derivations with the Matsubara dynamics frameworkHele et al. (2015); T. J. H. Hele and Althorpe (2015) which does not subject to any restriction to equilibrium conditions.R.Welsch et al. (2016) In that work, RPMD is used to investigate the non-equilibrium (photoinduced) adiabatic dynamics on a single electronic state, and has shown to accurately describe quantum dynamics compared to the numerically exact results.R.Welsch et al. (2016) Further, the original RPMD method,Menzeleev, Ananth, and Miller (2010) the MV-RPMD approachJ.R.Duke and Ananth (2015) and the state-dependent centroid molecular dynamics approachLiao and Voth (2002) have already been used to simulate non-adiabatic dynamics under non-equilibrium initial conditions. These early studies inspire us to investigate the numerical performance of NRPMD for simulating photoinduced non-equilibrium dynamics.

To compute the time-dependent reduced density matrix ${\rho}_{jj}(t)$ in Eqn. 19 we propose the following NRPMD population expression

[TABLE]

Here, $d{\boldsymbol{\tau}}\equiv\int d\{{\bf R}_{\alpha}\}d\{{\bf P}_{\alpha}\}d\{{\bf q}_{\alpha}\}d\{{\bf p}_{\alpha}\}$, with shorthand notation $d\{{\boldsymbol{\chi}}_{\alpha}\}=\prod_{\alpha=1}^{N}d{\boldsymbol{\chi}}_{\alpha}$. In addition, $\mathcal{P}^{0}(\{{\bf q}_{\alpha}(0),{\bf p}_{\alpha}(0)\})$ represents the distribution of the initial electronic mapping variables, $\rho_{\mathrm{rp}}(\{{\bf R}_{\alpha}(0),{\bf P}_{\alpha}(0)\})$ is the ring-polymer density for the initial nuclear density operator $\hat{\rho}_{\mathrm{n}}$, and $\mathcal{\bar{P}}_{j}(t)$ is the time-dependent population estimator. The above proposed expression is akin to the excited states population used in the MV-RPMD approach.J.R.Duke and Ananth (2015)

The electronic population estimator $\mathcal{\bar{P}}$ has many possible choices.J.R.Duke and Ananth (2015); Richardson and Thoss (2013); Richardson et al. (2017); Hele and Ananth (2016) Here, we use the following estimator

[TABLE]

which was originally proposed in the NRPMD approachRichardson and Thoss (2013) and recently derived in various state-dependent RPMD methods,Hele and Ananth (2016); Richardson et al. (2017) and it is similar to the original MMST population expression.Meyer and Miller (1979); Stock and Thoss (1997)

The initial mapping density $\mathcal{P}^{0}({\bf q}_{\alpha}(0),{\bf p}_{\alpha}(0))$ is required to describe the initial electronic population $\rho_{jj}(0)=\delta_{ij}$, which corresponds to the initial electronic state $|i\rangle$. Here, we adapt the focused initial condition used in the MV-RPMD approachJ.R.Duke and Ananth (2015); Duke and Ananth (2016) and linearized path-integral methodsBonella and Coker (2003, 2005); Dunkel, Bonella, and Coker (2008) to represent a constrained mapping density

[TABLE]

The above expression requires the mapping variables to satisfy the following relation

[TABLE]

This can be viewed as the “Bohr-Sommerfeld” quantization relationMiller and S.J.Cotton (2016) that has been used to initialize mapping variables.Meyer and Miller (1979); Muller and Stock (1999); W.H.Miller (2009); Miller and Cotton (2016) Solving the above equation provides the values of the action variables $[{\bf q}_{\alpha}]_{j}^{2}+[{\bf p}_{\alpha}]_{j}^{2}$ with a value 3 for the occupied state $|i\rangle$, or 1 for unoccupied other states $|j\rangle$, whereas the angle variables $[{\boldsymbol{\theta}}_{\alpha}]_{j}=-\tan^{-1}([{\bf p}_{\alpha}]_{j}/[{\bf q}_{\alpha}]_{j})$ are randomly sampledJ.R.Duke and Ananth (2015); Miller and Cotton (2016) within the range of $[0,2\pi]$. Other choices, such as Window estimators,Cotton and Miller (2016, 2013); Miller and Cotton (2016); Cotton and Miller (2019) Wigner transformed projection operators,Saller, Kelly, and Richardson (2019) or MV-RPMD based estimatorsAnanth (2013); J.R.Duke and Ananth (2015); Duke and Ananth (2016); Pierre et al. (2017) are possible and subject to future investigations.

To compute non-adiabatic dynamics, real-time trajectories are propagated based on Eqn. 17-18, with the mapping and nuclear initial distributions sampled from $\mathcal{P}^{0}(\{{\bf q}_{\alpha}(0),{\bf p}_{\alpha}(0)\})$ and $\rho_{\mathrm{rp}}(\{{\bf R}_{\alpha}(0),{\bf P}_{\alpha}(0)\})$, respectively. The time-dependent population is computed from the ensemble average described in Eqn. 20.

II.3 Computational Details

Model Systems. In this paper, we adapt two widely used model systems to investigate the performance of NRPMD for simulating excited state non-adiabatic dynamics. Model I is a widely used three-state Morse potential for photo-dissociation dynamics.Coronado, Xing, and Miller (2001) The Hamiltonian operator $\hat{H}=\hat{P}^{2}/2M+\hat{V}$ of Model I has the following potential

[TABLE]

Here, $V_{ii}$ and $V_{ij}$ are diabatic potentials and couplings, respectively. Parameters of Model I are provided in Table 1. These potentials and couplings are visualized in Fig. 1(a)-(c). The nuclear mass is $M=20,000$ a.u.

Model II is a one dimensional spin-boson systemSato, Kelly, and Rubio (2018); Richardson and Thoss (2013); Ananth (2013) with the following Hamiltonian

[TABLE]

In the above expression, $\hat{\sigma}_{z}$ and $\hat{\sigma}_{x}$ are Pauli spin matrices in the electronic subspace $\{|1\rangle,|2\rangle\}$, $\Delta=1$ a.u., $\hat{R}$ and $\hat{P}$ are the position and momentum operators of a harmonic boson mode with frequency $\omega=1$ a.u., and the nuclear mass is $M$=1 a.u. The two-level system and the boson mode interact with each other through a bi-linear coupling potential with a constant coupling strength $\gamma$.

Initial Conditions. For all results presented in this work, the initial photo-excitation is modeled by the following density operator

[TABLE]

where $|1\rangle$ is the first diabatic state in both models, and the initial nuclear density operator is

[TABLE]

The above initial nuclear density is chosen as the canonical thermal density associated with the ground state Hamiltonian

[TABLE]

where $M$ is the nuclear mass, and $R_{0}$ is the position of the Franck-Condon vertical excitation. Here, we assume that the ground state and the excited states are electronically decoupled; there is no communication between them except during the initial Franck-Condon excitation. For Model I, the initial excitation is indicated with black arrows in Fig. 1a-c, with $R_{0}$=2.1, 3.3, 2.9 for Model IA, IB, and IC, respectively. The frequency of the ground state is chosenJ.R.Duke and Ananth (2015) as $\omega_{0}=0.005$ a.u., and $\beta$ is the inverse temperature that corresponds to 300 K. For Model II, $R_{0}=0$, the frequency of the ground state is $\omega_{0}=1$ a.u., and the inverse temperature is $\beta=16$ a.u., such that the initial nuclear quantum distribution is significantly different compared to the classical one, and the dynamics cannot be treated accurately through the classical Wigner model (that samples the initial Wigner distribution and propagate the trajectories classically).

Simulation Details. For all NRPMD results presented in this paper, a total of $10^{4}$ trajectories are used to generate the converged population dynamics. The initial thermal nuclear density $\hat{\rho}_{n}$ is sampled by the normal-mode path-integral Monte-Carlo (PIMC)Tuckerman et al. (1993) in the ground state $\hat{H}_{\mathrm{g}}$ (Eqn. 28), which generates the ring-polymer initial density

[TABLE]

In the above equation, $\beta_{N}=\beta/N$ with $N$ as the number of beads (imaginary-time slices), and the ring-polymer Hamiltonian associated with the ground state $\hat{H}_{\mathrm{g}}$ is expressed as follows

[TABLE]

In this study, we follow the recent works of state-dependent RPMD that treat $N$ as a convergence parameter.Richardson et al. (2017); J.R.Duke and Ananth (2015) For the sampling of the nuclear initial condition (Eqn. 29), a large enough $N$ is used to ensure a converged $\rho_{\mathrm{rp}}({\bf R}_{\alpha},{\bf P}_{\alpha})$. This requires $N=4$ for Model I and $N=16$ for Model II, and $N$ remains fixed for the dynamics propagation. We have also performed convergence tests for population dynamics with a higher number of beads, suggesting that these choices of $N$ are sufficient to provide converged results. The initial conditions for the mapping variables are sampled based on Eqn. 23. Each configuration is then propagated with Eqn. 17-18. A symplectic integration scheme is used to numerically propagate the dynamics,Kelly et al. (2012); Church et al. (2018) although other simpler schemeRichardson et al. (2017) generates the same numerical results for the model calculations studied here.

The real-time NRPMD dynamics governed by $H_{N}$ (Eqn. 16) requires a fictitious temperature $\beta$ as the parameter of the dynamics.de la Pe$\tilde{\mathrm{n}}$a (2014); J.R.Duke and Ananth (2015); Shakib and Huo (2017) For Model I, we follow the previous MV-RPMD workJ.R.Duke and Ananth (2015) on choosing $\beta^{-1}$, which is the energy gap between the zero-point energy (ZPE) of the ground state plus the potential energy gap between the lowest excited state and the initially occupied excited state at $R_{0}$. This provides the fictitious temperatures of 15288 K, 9605 K, and 8843 K for Model IA, IB, and IC, respectively.J.R.Duke and Ananth (2015) For Model II, we directly use $\beta=16$ (associated with the initial distribution) during the NRPMD simulation.

Numerically exact results for Model I are obtained from the discrete variable representation (DVR) calculations,Colbert and Miller (1992) with a grid spacing of 0.009 a.u. in the range of $R\in[0.5,20]$ a.u. to ensure convergence. For Model II, exact results are directly obtained by computing the reduced density matrix $\rho_{jj}(t)=\mathrm{Tr}_{\mathrm{n}}[\hat{\rho}_{0}e^{i\hat{H}t/\hbar}|j\rangle\langle j|e^{-i\hat{H}t/\hbar}]$ evaluated with the basis $|i\rangle\otimes|n\rangle$, where $|i\rangle$ is the diabatic basis in Model II and $|n\rangle$ is the eigenstate of the harmonic oscillator centered at $R=0$.

III Results and Discussions

Fig. 1 presents the population dynamics of Model IA-IC for photoinduced dissociation dynamics. Panels (a)-(c) present the diabatic potentials. Black arrows indicate the Franck-Condon vertical excitations. Panels (d)-(f) provide the diabatic populations with the same color coding used in (a)-(c) for the diabatic states, obtained from NRPMD simulations (open circles) as well as numerically exact results (solid lines). For all three cases, NRPMD provides a reasonable agreement with the exact results for both the short-time relaxation and non-adiabatic branching dynamics, as well as the longer time asymptotic populations. We emphasize that Model I is a challenging test case for many approximate non-adiabatic dynamics approachesCoronado, Xing, and Miller (2001); Huo and Coker (2012a); J.R.Duke and Ananth (2015); Tao (2016) due to its highly anharmonic potential and non-linear diabatic couplings. Nevertheless, NRPMD provides accurate predictions for the key features of these non-adiabatic events associated with multiple curve crossings.

Fig. 2 presents the state-dependent population dynamics of Model IA computed from various recently developed dynamics method, including (a) NRPMD (same as the result presented in Fig. 1d), (b) partial linearized density matrix (PLDM) path-integral approach,Huo and Coker (2012a, 2011) (c) MV-RPMD approach,J.R.Duke and Ananth (2015); Ananth (2013) and (d) symmetrical quasi-classical (SQC) approach.Tao (2016); Miller and S.J.Cotton (2016); Cotton and Miller (2013) In panel (b), PLDM generates accurate short-time non-adiabatic branching dynamics, but starts to deviate from the exact result at a longer time, potentially due to the less accurate partial linearization approximation at a longer time.Huo and Coker (2012a, 2011) This is confirmed by using an iterative version of PLDM,Huo and Coker (2012b); Lee, Huo, and Coker (2016) which only requires linearization approximation for a short time propagator, and then concatenate these short-time PLDM propagators together by important sampling technique.Huo and Coker (2012b); Lee, Huo, and Coker (2016); Dunkel, Bonella, and Coker (2008) The iterative-PLDM provides accurate dynamics for Model I, although a large number of trajectories are required to converge the results.Huo and Coker (2012b) Further, the PLDM Hamiltonian can be viewed as a particular limit of the NRPMD Hamiltonian (or the CS-RPMD HamiltonianS.N.Chowdhury and Huo (2017)) with one nuclear bead, and two mapping beads for describing the forward and backward propagations.Huo and Coker (2011) NRPMD which uses multiple beads for all DOFs seems to provide more accurate short-time branching dynamics as well as long-time populations.

In Fig. 2c, MV-RPMD approachJ.R.Duke and Ananth (2015) provides less accurate non-adiabatic branching dynamics, probably due to the presence of the inter-bead coupling of the mapping variables which contaminates the electronic dynamics.Richardson and Thoss (2013); Ananth (2013) Further, its inability to correctly capture electronic coherenceAlthorpe et al. (2016) could also mitigate its accuracy of describing the quantum branching dynamics. The NRPMD Hamiltonian, on the other hand, does not contain any inter-bead coupling for mapping variables, thus can reliably capture electronic coherence dynamicsRichardson and Thoss (2013); Richardson et al. (2017) and provide accurate non-adiabatic population transfer as shown in (a).

In Fig. 2d, the recently developed SQCTao (2016); Cotton and Miller (2013); Miller and Cotton (2016) approach provides less accurate results for this model, compared to the other three approaches. Despite that SQC provides accurate non-adiabatic dynamics for many model systems,Miller and Cotton (2016) the population dynamics for Model IA starts to deviate from the exact results even at a very short time. That is also the case for Model IB and IC.Tao (2016) Further, the closely related LSC-IVR approach generates very similar results compared to SQC, except some negative populations.Coronado, Xing, and Miller (2001); Huo and Coker (2012a); J.R.Duke and Ananth (2015) These results suggest that the Ehrenfest type of the nuclear force, together with the nuclear Wigner distribution (assumed by both SQC and LSC-IVR) might be the cause for this less accurate dynamics. A recently proposed coherence-controlled SQC (cc-SQC) approachTao (2016) has significantly improved the accuracy by using different nuclear forces based on the time-dependent action variables.Tao (2016)

Fig. 3 presents the population difference $\langle\sigma_{z}(t)\rangle$ between state $|1\rangle$ and $|2\rangle$ in Model II with three different electron-phonon coupling strength $\gamma/\Delta$. Here, we compare the dynamics obtained from NRPMD (red) with multi-trajectory Ehrenfest dynamics (MTEF) (green), PLDM (blue), and numerically exact simulations (black). Panel (a) presents the results in a weak electron-phonon coupling regime with $\gamma/\Delta$=0.1. The temperature $\beta=16$ is low enough such that the initial quantum distribution is significantly different compared to the classical distribution. Thus, quantum mechanical treatment of the nuclear DOF is required for trajectory-based approaches, through either a Wigner initial distribution (for MTEF and PLDM) or the ring polymer quantization (for NRPMD). One can see that all three approximate methods behave accurately compared to the numerically exact results and reproduce correct oscillations and damping patterns up to $t$ =15 a.u. After that, both the PLDM and MTEF approach failSato, Kelly, and Rubio (2018) to accurately describe the longer time recurrence of oscillation in $\langle\sigma_{z}(t)\rangle$. This less accurate longer time dynamics might be caused by the zero point energy (ZPE) leakage problem,Liu and Miller (2011); Habershon and Manolopoulos (2009) which is typical for linearized path-integral approaches based on the classical Wigner dynamics.Poulsen, Nyman, and Rossky (2003); Shi and Geva (2004); W.H.Miller (2009) This ZPE leakage originates from the fact that classical dynamics does not preserve the ZPE incorporated in the initial Wigner distribution,Liu and Miller (2011); Habershon and Manolopoulos (2009) causing an incorrect energy flow from the nuclear DOF to the electronic DOF,Muller and Stock (1999) equalizing the longer time populations and giving $\langle\sigma_{z}(t)\rangle$=0. Compared to the classical Wigner dynamics of the nuclear DOF,W.H.Miller (2009); Shi and Geva (2004); Poulsen, Nyman, and Rossky (2003) quantizing the nuclear DOF with ring polymer can effectively incorporate nuclear ZPE and describe tunneling effectsShakib and Huo (2017); Shushkov, Li, and Tully (2012); Tao, Shushkov, and Miller (2018, 2019) and alleviate ZPE leaking problem, thus reliably provide the longer time recurrence of the oscillating population.

Fig. 3b-c presents the population dynamics for stronger electron-phonon couplings. We can see that NRPMD method reproduces the exact result fairly well up to $t$=5 a.u., especially for the model calculation presented in panel (c). At a longer time, however, NRPMD becomes less accurate compared to the exact results, missing the recurrence of the oscillations. These deviations at a longer time in panel (b)-(c) might be due to the intrinsic quantum coherence of nuclear dynamics, which is missed by NRPMD but can be well captured by methods that employ coupled trajectories.Sato, Kelly, and Rubio (2018) Nevertheless, as an independent trajectory-based approach, NRPMD still outperforms both MTEF and PLDM for all model calculations presented here.

We emphasize that the success of any RPMD-based approach relies on the separation of the time-scale between the high-frequency vibrations of the ring polymer and the dynamics of physical interest.Habershon et al. (2013) The high-frequency ring polymer oscillations could contaminate the real-time dynamics of the nuclei, which in turn influence the electronic quantum dynamics. This issue could potentially impact the accuracy of our population dynamics, for example, slightly shifting the electronic Rabi oscillation of $\langle\sigma_{z}(t)\rangle$ in Fig. 3. To address this issue, we incorporate the recently developed thermostatting techniqueM.Ceriotti et al. (2010); M.Rossi, M.Ceriotti, and Manolopoulos (2014); Hele (2016) to NRPMD, with details provided in Appendix B. Thermostatting has been shown to successfully remove these contaminations from spurious high-frequency oscillations,Hele (2016) and the model calculations provided in Appendix B indeed confirms that correct electronic Rabi frequency is recovered when applying a Langevin thermostat on the nuclear ring polymer normal modes. Thus, thermostating offers a valuable and practical approach to remove potential contaminations from the high-frequency ring polymer vibrations, enables the promising NRPMD approach to provide more accurate excited states non-adiabatic dynamics.

IV Conclusion

In this paper, we present a rigorous derivation of the NRPMD Hamiltonian which was originally proposed in the NRPMD approach.Richardson and Thoss (2013) Our derivation uses the MMST mapping representation for the electronic DOF, and the ring polymer path-integral quantization for the nuclear DOF. The NRPMD Hamiltonian, together with the previously derived CS-RPMD Hamiltonian,S.N.Chowdhury and Huo (2017) can be viewed as a unified theory for classical electronic states (through the MMST mapping description) and classical nuclei (through the ring polymer quantization).

We further propose to compute excited state non-adiabatic quantum dynamicsJ.R.Duke and Ananth (2015); R.Welsch et al. (2016) with NRPMD Hamiltonian. Numerical results with the coupled Morse potential and one-dimensional spin-boson model suggest that NRPMD can provide accurate short-time branching dynamics and a reliable longer time dynamics. The NRPMD Hamiltonian does not contain inter-bead coupling terms associated with the electronic mapping variables; thus it is capable to accurately capture the electronic quantum dynamics.Richardson and Thoss (2013); Richardson et al. (2017) In this particular model calculation, NRPMD outperforms the recently developed MV-RPMD methodsJ.R.Duke and Ananth (2015); Ananth (2013) which does include these inter-bead coupling terms that might contaminate the electronic dynamics. Compared to the linearized semi-classical methods based on Wigner quantization of nuclei,W.H.Miller (2009); Shi and Geva (2004); Huo and Coker (2011) quantizing nuclear DOF with ring polymer can effectively incorporate nuclear zero-point energy (ZPE) and alleviates ZPE leakage problem.

This work opens up new possibilities of using state-dependent RPMD approachesRichardson and Thoss (2013); Richardson et al. (2017); S.N.Chowdhury and Huo (2017) to accurately simulate electronic non-adiabatic dynamics. These approaches are potentially well-suited theoretical methods to investigate photochemical reactions, especially when nuclear quantum effects play an important role, such as in the photoinduced proton-coupled electron transfer reactions.Venkataraman, Soudackov, and Hammes-Schiffer (2010); Hazra, Soudackov, and Hammes-Schiffer (2010, 2011); Goyal and Hammes-Schiffer (2015); Goyal et al. (2016); Goyal and Hammes-Schiffer (2017); Mandal, Shakib, and Huo (2018) We note that the original RPMD method is limited to one electron non-adiabatic process,Habershon et al. (2013); Menzeleev, Ananth, and Miller (2010, 2011); Kretchmer and Miller (2013, 2016); Kenion and Ananth (2016) whereas the vibronic quantization approach (which uses an explicit quantum-state description of proton) can be numerically expensive for a three-dimensional quantum treatment of many protons.Hammes-Schiffer (2015); Hazra, Soudackov, and Hammes-Schiffer (2010, 2011); Goyal et al. (2016); Goyal and Hammes-Schiffer (2017); Mandal, Shakib, and Huo (2018); Mandal et al. (2019) State-dependent RPMD approaches, such as NRPMDRichardson and Thoss (2013) and CS-RPMD,S.N.Chowdhury and Huo (2017) provide accurate electronic non-adiabatic dynamics and nuclear quantum effects without any limitations on the number of electrons and protons that can be explicitly described.

We note that most of the state-dependent RPMD approaches, including the NRPMD method discussed here, are formulated in the diabatic representation. To perform on-the-fly simulation with adiabatic electronic structure calculations, these approaches are usually reformulated back to the adiabatic representation. This process requires nontrivial theoretical efforts, and the adiabatic equation of motion are computationally inconvenient due to the presence of derivative couplings. These non-trivial tasks, however, can be avoided by using the recently developed quasi-diabatic propagation scheme.Mandal, Yamijala, and Huo (2018) The QD scheme uses the adiabatic states associated with a reference geometry as the local diabatic states during a short-time propagation step, and dynamically update the definition of the diabatic states along the time-dependent nuclear trajectory. This scheme thus allows a seamless interface between diabatic dynamics approaches (such as NRPMD) with adiabatic electronic structure calculations, providing new frameworks to accurately and efficiently perform non-adiabatic on-the-fly simulations.

Future investigations will also focus on analytic derivation of the NRPMD dynamics based on rigorous theoretical frameworks,Shi and Geva (2003); Jang, Sinitskiy, and Voth (2014); Hele et al. (2015); T. J. H. Hele and Althorpe (2015) such as the Matsubara dynamicsHele et al. (2015); T. J. H. Hele and Althorpe (2015) and the exact mapping Liouvillian.Hele and Ananth (2016) These formal theoretical derivations will help to assess the validity and the accuracy of the state-dependent RPMD approaches for simulating excited states non-adiabatic dynamics.

V Acknowledgement

This work was supported by the National Science Foundation CAREER Award under Grant No. CHE-1845747. Computing resources were provided by the Center for Integrated Research Computing (CIRC) at the University of Rochester.

VI Appendix A: Details of the derivation for NRPMD partition function

In this appendix, we provide details of derivations for (i) the expression of $[{\hat{H}}_{e}({\bf R}_{\alpha})]_{\mathrm{w}_{\alpha}}$ in Eqn. 12, (ii) the expression of the electronic partition function $\mathrm{Tr_{e}}[\prod_{\alpha=1}^{N}e^{-\beta_{N}\hat{H}_{e}({\bf R}_{\alpha})}]$ in Eqn. 13, and (iii) the expression of $\boldsymbol{\Gamma}$ in Eqn. 15.

First, we provide the derivation of the $[{\hat{H}}_{\mathrm{e}}({\bf R}_{\alpha})]_{\mathrm{w}_{\alpha}}$ expression. The Wigner transformed mapping Hamiltonian operator (Eqn. 6) can be expressed as follows

[TABLE]

Between the first and the second equality, we have used the fact $\langle q-\frac{\Delta}{2}|\hat{p}|q+\frac{\Delta}{2}\rangle=\int dpp\langle q-\frac{\Delta}{2}|p\rangle\langle p|q+\frac{\Delta}{2}\rangle=(\frac{1}{2\pi\hbar})\int dppe^{-\frac{i}{\hbar}p\Delta}=(\frac{1}{2\pi\hbar})i\hbar\frac{\partial}{\partial\Delta}\int dpe^{-\frac{i}{\hbar}p\Delta}=i\hbar\frac{\partial}{\partial\Delta}\int dp\langle q-\frac{\Delta}{2}|p\rangle\langle p|q+\frac{\Delta}{2}\rangle=i\hbar\frac{\partial}{\partial\Delta}\langle q-\frac{\Delta}{2}|q+\frac{\Delta}{2}\rangle$. Between the second and the third equality, we use integration by parts for $\int d{\bf\Delta}_{\alpha}$. Analytically perform the rest integrals, we arrived at the final expression of $[{\hat{H}}_{e}({\bf R}_{\alpha})]_{\mathrm{w}_{\alpha}}$ in Eqn. 12. A similar derivation has been provided in the previous work by Kapral and co-workersBonella and Kapral (2010) for developing the Wigner mapping mixed quantum-classical Liouville (MQCL) equation.

Second, we provide the details of the derivation for the electronic trace expression in Eqn. 13. We start by inserting a resolution of identity $\int d{\bf q}_{2}\mathcal{P}|{\bf q}_{2}\rangle\langle{\bf q}_{2}|$ into $\left[\prod_{\alpha=2}^{N}e^{-\beta_{N}\hat{H}_{e}({\bf R}_{\alpha})}\mathcal{P}\right]_{\mathrm{w}_{1}}$, resulting in the second line of Eqn. VI. Next, we alter the order of terms in line 3-4 of Eqn. VI, leading to an integral over ${\bf q}_{2}$. Finally, replacing $\int d{\bf q}_{2}$ by the corresponding Wigner transform,J.R.Duke and Ananth (2015) we can use the property in Eqn. 9 again and repeat the same procedure outlined in Eqn. 10 to factorize the total Wigner function as a product of two Wigner functions inside the $\int d{\bf q}_{2}d{\bf p}_{2}$ integral, arriving at the last line of Eqn. VI.

[TABLE]

A similar procedure of the above derivation has been recently used to derive the population estimator in the MV-RPMD approach.J.R.Duke and Ananth (2015)

Third, we derive the expression of $\boldsymbol{\Gamma}$ (Eqn. 15). We start from the last line of Eqn. 13 and define

[TABLE]

where $\mathcal{P}=\sum_{n}|n\rangle\langle n|$ is the projection operator in the SEO basis. Recall that the SEO mapping wavefunction is the product of $(\mathcal{K}-1)$ ground state harmonic oscillator wavefunctions and one excited state harmonic oscillator wavefunction

[TABLE]

With the above relation, we can rewrite Eqn. 33 as

[TABLE]

Rearranging the prefactors of the above equation and grouping terms associated with $\boldsymbol{\Delta}_{\alpha}$, we have

[TABLE]

Analytically performing the integration over $\boldsymbol{\Delta}_{\alpha}$ (a Gaussian integral), we obtain the final expression in Eqn. 15. Similar derivations can also be found in the previous work of MV-RPMDAnanth (2013) as well as in the recently derived exact mapping variable Liouvillian.Hele and Ananth (2016)

VII Appendix B: Thermosttated NRPMD (T-NRPMD)

In this appendix, we investigate the effect of thermostatting nuclear ring polymer on the NRPMD non-adiabatic dynamics. It is known that ring polymer quantization often introduces spurious frequencies in RPMD dynamics due to the presence of the high-frequency normal mode vibrations,Habershon et al. (2013) and causes the “spurious resonance problem” for computing spectraM.Ceriotti et al. (2010); M.Rossi, M.Ceriotti, and Manolopoulos (2014); Hele (2016) and introducing incorrect frequency in time-correlation functions for nonlinear operators.Hele (2016) Thus, the success of any RPMD-type approach relies on the separation of the time-scale between the high-frequency normal mode vibrations of the ring polymer and the dynamics of physical interest.Habershon et al. (2013) Various thermostatted RPMD (TRPMD) approachesM.Ceriotti et al. (2010); M.Rossi, M.Ceriotti, and Manolopoulos (2014) are proposed to achieve this.M.Ceriotti et al. (2010); M.Rossi, M.Ceriotti, and Manolopoulos (2014)

Based on the Matsubara dynamics framework, it is recently shown that this frequency contamination arises due to discarding the imaginary term of the Matsubara LiouvillianHele et al. (2015); T. J. H. Hele and Althorpe (2015) when deriving the RPMD approach. This formal analysisHele (2016) shows that TRPMD can be justified by replacing the imaginary Matsubara Liouvillian with a friction term, such as the Fokker-Planck operator, Bussi and Parrinello (2007); Gardiner (2003); Risken (1989) instead of just discarding the imaginary part of the Matsubara Liouvillian as done in RPMD. Here, we apply a Langevin thermostat that couples to the nuclear normal mode in NRPMD. We briefly introduce the normal mode representation of the ring polymer, before we provide the equation of motion for thermostatting.

The free ring-polymer Hamiltonian (see Eqn. 3 and below) is defined as follows

[TABLE]

Often, the dynamical propagation of RPMD (and PIMD) can be simplified by transforming $H_{\mathrm{rp}}$ from the above bead representation (or so-called the primitive nuclear coordinate) to the normal mode representation, which is the eigenstate of the Hessian matrix of $H_{\mathrm{rp}}$. Diagonalizing the Hessian matrix of $H_{\mathrm{rp}}$ provides the eigenvalue, i.e., the following normal mode frequency

[TABLE]

where $\mu\in[0,N-1]$ represents the index of the normal mode. The same diagonalization process also gives the eigenvector $T_{\alpha\mu}$ of the Hessian matrix, which provides the relation between the primitive coordinate $\{{\bf R}_{\alpha}\}$ and the normal mode coordinate $\{\widetilde{\bf R}_{\mu}\}$, as well as the corresponding relation for momenta under two representations. These relations are expressed as follows

[TABLE]

The above transformation matrix elements have the following values

[TABLE]

Under the normal mode representation, the free ring polymer Hamiltonian $H_{\mathrm{rp}}$ in Eqn. 37 becomes

[TABLE]

where the normal mode frequency $\widetilde{\omega}_{\mu}$ is described in Eqn. 38. Note that the inter-bead coupling terms of the ring polymer become a set of simple quadratic terms with the normal mode frequencies. The nuclear equation of motion of NRPMD described in Eqn. 17 under the normal mode representation is expressed as follows

[TABLE]

where simple chain rule is used to establish the last equality, $-\nabla_{{\bf R}_{\alpha}}H_{N}(\{{\bf R}_{\alpha}\})=\dot{\bf P}_{\alpha}$ is the nuclear force in Eqn. 17, and ${T}_{\alpha\mu}={\partial{\bf R}_{\alpha}}/{\partial\widetilde{\bf R}_{\mu}}$ is the Jacobian matrix element of the transformation between the primitive nuclear variables $\{{\bf R}_{\alpha}\}$ and the normal mode coordinates $\{\widetilde{\bf R}_{\mu}\}$ described in Eqn. 40. Note that ${T}_{\alpha\mu}$ is same for all $\mathcal{F}$ nuclear DOF.

The normal mode NRPMD nuclear force $-\nabla_{\widetilde{\bf R}_{\mu}}H_{N}(\{{\bf R}_{\alpha}\})$ in Eqn. 43 contains three types of terms: (i) the force contribution from the free ring polymer, $-{\bf M}\omega_{\mu}^{2}\widetilde{\bf R}_{\mu}$, (ii) the state independent force, $-\sum_{\alpha}\nabla_{{\bf R}_{\alpha}}V_{0}({\bf R}_{\alpha}){T}_{\alpha\mu}$, and (iii) the state dependent force, -\frac{1}{2\hbar}\sum_{\alpha}\sum_{nm}\nabla_{{\bf R}_{\alpha}}V_{nm}({\bf R}_{\alpha})\big{(}[{\bf q}_{\alpha}]_{n}[{\bf q}_{\alpha}]_{m}+[{\bf p}_{\alpha}]_{n}[{\bf p}_{\alpha}]_{m}-\delta_{nm}\hbar\big{)}{T}_{\alpha\mu}.

Following the previous work of TRPMD,M.Ceriotti et al. (2010); Hele (2016) the nuclear ring polymer normal mode $\{\widetilde{\bf R}_{\mu}\}$ in $H_{N}$ is coupled to an Langevin thermostat, giving a method that we referred as the Thermostatted NRPMD (T-NRPMD). In T-NRPMD, the mapping equations of motion remain the same as described in Eqn 18, whereas the nuclear equation of motion in Eqn. 43 is replacedM.Ceriotti et al. (2010); Hele (2016) by the following one

[TABLE]

The first term is the force from the NRPMD Hamiltonian (Eqn. 43), the second term is the friction force acting on $\widetilde{\bf P}_{\mu}$ with ${\boldsymbol{\eta}}_{\mu}$ as the bead-specific normal mode friction matrix, and the last term is the random force, with ${\boldsymbol{\xi}}_{\mu}(t)$ representing an uncorrelated, Gaussian-distributed random forceM.Ceriotti et al. (2010) with unit variance $\langle{\boldsymbol{\xi}}_{\mu}(0){\boldsymbol{\xi}}_{\mu}(t)\rangle=\delta(t)$, and zero mean $\langle{\boldsymbol{\xi}}_{\mu}(t)\rangle=0$. Based on the recent analysis of TRPMD from the Matsubara dynamics framework,Hele (2016) we choose the same friction constant for all $\mathcal{F}$ nuclear DOF associated with the $\mu_{\mathrm{th}}$ normal mode, with the $\mu$-specific friction term, $\eta_{\mu}=2\lambda|\widetilde{\omega}_{\mu}|$, where $\lambda$ is viewed as a parameter.Hele (2016) The Langevin equation can be numerically propagated based on the algorithm in previous works,M.Ceriotti et al. (2010); Bussi and Parrinello (2007) whereas the mapping equation of motion is integrated with a symplectic integrator.Kelly et al. (2012); Church et al. (2018)

Fig. 4 presents the results of Model II obtained from T-NRPMD. Here, we investigate the effects of nuclear thermostatting on the excited state non-adiabatic dynamics and explore the impact of various friction constant $\lambda$ on the population dynamics. The results are obtained by using $\lambda=0$ (red), i.e., the NRPMD approach (same as results shown in Fig. 3), $\lambda=1.0$ (green), and $\lambda=5.0$ (blue), together with the numerically exact results (black dashed lines). Despite that NRPMD provides more accurate population dynamics compared to classical Wigner based methods as discussed in Fig. 3, the population tends to oscillate with a slightly shifted frequency compared to the electronic Rabi frequency. This can be clearly seen by comparing the results of NRPMD (red) and the exact ones (black) in panels (b),(d) and (f) which depict the magnified plots corresponding to the square regions in panels (a), (c), and (e). The higher frequency normal modes of the nuclear ring polymer might be the source of these spurious oscillations,Hele (2016) which has shown to contaminate the nuclear dynamics.Hele (2016)

Using $\lambda=1.0$ (green), T-NRPMD recovers the correct oscillation frequencies of the electronic population. This value of $\lambda$ is chosen based on a friction parameter that is derived from achieving the correct nuclear oscillation frequency in a harmonic potential.Hele (2016) By recovering the correct nuclear oscillations of $R$, the coupled electronic dynamics is also improved. We also observed that by applying a small friction parameter $0.5\leqslant\lambda\leqslant 1$, T-NRPMD already improves the dynamics and recover the correct oscillation period. The population dynamics continues to oscillate with the correct frequency when further increase the friction parameter to the overdamped regime with $\lambda=5.0$ (solid blue line). The correct electronic oscillation is likely stemmed from the correct nuclear oscillations when applying thermostat which has been demonstrated in the previous work.Hele (2016)

These investigations demonstrate that T-NRPMD is a valuable tool for providing accurate excited state dynamics and alleviate spurious frequency problem associated with the ring polymer quantization. Future studies include rigorous derivation of T-NRPMD approach through the Matsubara dynamics frameworkHele et al. (2015); T. J. H. Hele and Althorpe (2015) with the mapping LiouvillianHele and Ananth (2016) to treat electronic states explicitly.