Investigating Photoinduced Proton Coupled Electron Transfer Reaction using Quasi Diabatic Dynamics Propagation
Arkajit Mandal Farnaz A. Shakib Pengfei Huo






pengfei.huo@rochester.edu
Department of Chemistry, University of Rochester, 120 Trustee Road, Rochester, New York 14627, United States


We investigate photoinduced proton-coupled electron transfer (PI-PCET) reaction through a recently developed quasi-diabatic (QD) quantum dynamics propagation scheme. This scheme enables interfacing accurate diabatic-based quantum dynamics approaches with adiabatic electronic structure calculations for on-the-fly simulations. Here, we use the QD scheme to directly propagate PI-PCET quantum dynamics with the diabatic Partial Linearized Density Matrix (PLDM) path-integral approach with the instantaneous adiabatic electron-proton vibronic states. Our numerical results demonstrate the importance of treating proton quantum mechanically in order to obtain accurate PI-PCET dynamics, as well as the role of solvent fluctuation and vibrational relaxation on proton tunneling in various reaction regimes that exhibit different kinetic isotope effects. This work opens the possibility to study the challenging PI-PCET reactions through accurate diabatic quantum dynamics approaches combined with efficient adiabatic electronic structure calculations. 


Investigating Photoinduced Proton Coupled Electron Transfer Reaction using Quasi Diabatic Dynamics Propagation
    Pengfei Huo
    
==============================================================================================================



§ INTRODUCTION

Photoinduced proton-coupled electron transfer (PI-PCET) reactions<cit.> involve the coupled transfer of both electron and proton upon photoexcitation. Thus, PI-PCET is fundamentally different from the extensively studied photoinduced proton transfer (PT) or electron transfer (ET) reactions. Several recent experimental and theoretical studies have revealed PI-PCET reactions in a wide range of systems,<cit.> such as simple hydrogen-bonded organic complexes,<cit.> organometallic complexes,<cit.> enol-keto tautomerization,<cit.> as well as photocatalytic water oxidation on small nanoparticles.<cit.>

Initiated through photoexcitation process,  PI-PCET reactions play a critical role in solar energy conversion processes.<cit.> At the same time, they are promising for providing new and unique reactivities which are not directly accessible in regular thermally activated PCET reactions.<cit.> Thus, understanding the fundamental mechanistic principles of PI-PCET will allow tuning and controlling this reaction and provide design principles for more efficient solar energy conversion devices. 

Accurately and efficiently simulating PI-PCET reactions, however, remains a challenging theoretical task as it requires an explicit quantum mechanical description of both electronic nonadiabatic transitions and nuclear quantum effects for proton (such as tunneling and zero-point energy). Further, the charge distribution in the excited state can significantly deviate from the ground state,<cit.> leading to a highly non-equilibrium configuration of the solvent upon initial photoexcitaiton. Relaxation of the solvent initial configuration can then drastically affect the reaction course of the PI-PCET.<cit.> Thus, compared to the well-explored thermal PCET reactions,<cit.> the proper description of the non-equilibrium PI-PCET process is beyond the equilibrium rate constant expressions<cit.> and requires either a detailed time-dependent dynamics<cit.> or non-equilibrium Fermi's golden rule.<cit.>

Recent theoretical studies based on fewest-switches surface-hopping (FSSH)<cit.> simulations with electron-proton vibronic basis have made significant contributions to elucidating PI-PCET dynamics.<cit.> As a mixed quantum-classical (MQC) method, however, FSSH treats quantum and classical degrees of freedom (DOF) on different footings. This can generate artificial coherence that leads to incorrect ET dynamics<cit.> or the breakdown of the detailed balance.<cit.> Several modified FSSH methods,<cit.> especially decoherence corrected FSSH,<cit.> can resolve these issues and give accurate charge transfer dynamics upon carefully chosen schemes.<cit.> Other recent simulations of PI-PCET based on numerical exact methods such as hierarchical equations of motion (HEOM),<cit.> matching-pursuit/split-operator Fourier transform (MP/SOFT),<cit.> or ab initio multiple spawning (AIMS)<cit.> can provide accurate time-dependent reaction dynamics. However, the numerical costs of these calculations will limit their scope of applications in simulating complex systems with many electronic states and nuclear DOF. It is thus ideal to use accurate yet efficient trajectory-based quantum dynamics approaches, which are essentially different from traditional MQC methods, to directly simulate PI-CPET reaction dynamics.

In this study, we apply a recently developed quasi-diabatic (QD)<cit.> scheme to directly propagate PI-PCET quantum dynamics with diabatic Partial Linearized Density Matrix (PLDM) path-integral method.<cit.> As an example of recently developed accurate diabatic trajectory-based quantum dynamics approaches, PLDM uses a consistent dynamical footing for describing all DOFs. Further, we treat both transferring electron and proton quantum mechanically by describing them in their adiabatic vibronic states, and directly use these adiabatic states as the quasi-diabatic states during the QD-PLDM propagation. Our numerical results demonstrate the importance of treating proton quantum mechanically for obtaining accurate PI-CPET dynamics, as well as the crucial role of solvent fluctuations and vibrational relaxations that dictate the reaction pathways. We also demonstrate how various reaction regimes can exhibit different kinetic isotope effects (KIE). The QD scheme outlined in this paper does not require any tedious efforts for building strict diabatic system-bath models,<cit.> and can be directly generalized to perform ab-initio on-the-fly simulations.<cit.> This work opens the possibility for studying PI-PCET reaction in realistic chemical systems through combining accurate diabatic quantum dynamics approaches with efficient adiabatic electronic structure calculations. 



§ THEORY AND METHOD



 §.§ Partial Linearized Density Matrix Method

We provide a brief outline of the Partial Linearized Density Matrix (PLDM) path-integral approach.<cit.> We begin with expressing the total Hamiltonian as follows

    Ĥ = T̂ + V̂_el(r̂,R̂),

where r̂ and R̂ represent the electronic and nuclear coordinate operators, T̂ is the nuclear kinetic energy operator, and V̂_el represents the “electronic Hamiltonian". Under a set of strict diabatic states {|i⟩, |j ⟩} which do not explicitly depend on the nuclear configuration, the total Hamiltonian can be expressed as follows

    Ĥ = T̂ + ∑_ij^N V_ij(R̂)|i ⟩⟨ j|,

where V_ij(R̂)=⟨ i|V̂_el(r̂,R̂)|j⟩ is the state-dependent potential for the electronic Hamiltonian operator, and N is the total number of electronic states.

Using the Meyer-Miller-Stock-Thoss<cit.> (MMST) mapping representation to transform the discrete electronic states into continuous variables, we have |i ⟩⟨ j | →â_i^†â_j, where â^†_i =(q̂_i - ip̂_i)/√(2). With this transformation, the non-adiabatic transitions between electronic states are exactly mapped onto the classical motion of fictitious harmonic oscillators.<cit.> Thus, MMST mapping Hamiltonian provides a consistent classical footing for both electronic and nuclear DOFs.

Expressing the full density matrix operator with the real-time path-integral expression, then applying a partial linearization approximation<cit.> selectively to the nuclear DOF and keeping the explicit propagation of the electronic mapping DOF for both forward and backward paths, we arrive at the PLDM expression for computing reduced density matrix<cit.>

    ρ_ij(t)=Tr_ R[ρ̂(0)e^iĤt/ħ|i⟩⟨ j|e^-iĤt/ħ]
       ≈∑_kl∫ d Rd P/ 2πħd q d p d q' d p' G_0G'_0 [ρ̂(0)^W_kl]  T_ki(t) T'_jl(t),

where T_ki(t)=1/2(q_i(t)+ip_i(t))(q_k(0)-ip_k(0)) and T'_jl(t)=1/2(q_l(0)+ip_l(0))(q_j(t)-ip_j(t)) are the electronic transition amplitudes, and [ρ̂(0)^W_kl] is the partial Wigner transform (with respect to the nuclear DOF) of the kl^th matrix element of the initial total density operators ρ̂(0).  The initial distribution of electronic DOF is G_0( q, p)= e^-1/2∑_m (q_m^2 + p_m^2 ) and  G'_0( q', p') =  e^-1/2∑_n (q'_n^2 + p'_n^2 ).

Classical trajectories are used to evaluate the approximate time-dependent reduced density matrix in Eqn. <ref>. The forward mapping variables are evolved based on the Hamilton's equations of motion<cit.>

    q̇_i = ∂ h / ∂ p_i;    ṗ_i = -∂ h / ∂ q_i,

where h is the classical mapping Hamiltonian<cit.> with the following expression

    h( p, q, R)=12∑_ijV_ij(R)(p_ip_j+q_iq_j).

The backward mapping variables are propagated with the similar equations of motion governed by h( p', q', R). The nuclei are propagated with the PLDM force<cit.> with the following expression

    F_PL=-1 4∑_ij∇ V_ij(R)[p_ip_j+q_iq_j+p'_ip'_j+q'_iq'_j].


PLDM uses consistent dynamical footing for both electronic and nuclear DOFs and thus, accurately describes their coupled motion. In contrast,  widely used mixed quantum-classical methods such as Ehrenfest or FSSH<cit.> treat quantum and classical DOFs on different footings, which causes the breakdown of detailed balance<cit.> or creating the artificial electronic coherence.<cit.> In addition, compared to the closely related methods that fully linearize both mapping and nuclear DOFs,<cit.> PLDM retains full dynamical propagation along both forward and backward paths for the mapping DOF, thus achieving a more accurate description of the electronic dynamics.<cit.> PLDM has already been successfully applied to simulate a broad range of non-adiabatic processes, including excitation energy transfer dynamics,<cit.> electron transfer reactions,<cit.> singlet fission quantum dynamics,<cit.> and nonlinear optical spectroscopy calculations.<cit.>

It is worth mentioning that the numerical cost of PLDM scales as N^2, where N is the total number of electronic states. This is similar to recently developed methods, including symmetrical windowing quasi-classical (SQC)<cit.> and forward-backward trajectory solution to the quantum-classical Liouville equation (QCLE).<cit.> In contrast, the numerical cost of Redfield theory<cit.> or generalized quantum master equation (GQME)<cit.> scales as N^4, though more accurate results can be obtained under several particular parameter regimes.



 §.§ Quasi-Diabatic Propagation Scheme

Most of the routinely available electronic structure methods are formulated in the adiabatic representation. However, a large number of recently developed non-adiabatic dynamics methods,<cit.> including PLDM,<cit.> are formulated in the diabatic representation. Thus, the typical strategy for applying these new methods to “real” molecular systems is to reformulate them in the adiabatic representation,<cit.> which usually requires tedious theoretical efforts. Moreover, the adiabatic version of these methods are computationally inconvenient due to the presence of the first and second order derivative couplings,<cit.> which could potentially lead to numerical instabilities during dynamical propagations. 

To address this discrepancy between the accurate diabatic quantum dynamics approaches and routinely available electronic structure calculations in the adiabatic states, we have developed Quasi-Diabatic (QD) propagation scheme,<cit.> a general approach which allows interfacing adiabatic electronic structure calculations with diabatic trajectory-based quantum dynamics methods. Here, we provide a brief summary of this scheme, whereas the details of the algorithm can be found in Ref. . 

Consider a short-time propagation of the nuclear DOF during t∈[t_1, t_2], where the nuclear positions evolve from R(t_1) to R(t_2), with the corresponding adiabatic states {|Φ_α( R(t_1))⟩} and {|Φ_μ( R(t_2))⟩}. These adiabatic states are defined as the eigenstates of electronic part of the Hamiltonian

    V̂_el(r̂,R̂)|Φ_α( R)⟩=E_α ( R)|Φ_α( R)⟩,

where |Φ_α( R)⟩ are the adiabatic states with the corresponding eigenenergies E_α ( R), both of which explicitly depend on the nuclear coordinates.

The central idea of the QD propagation scheme<cit.>  is to use the nuclear geometry at time t_1 as the reference geometry, R_0≡ R(t_1), and the adiabatic basis {|Φ_α( R(t_1))⟩} as the quasi-diabatic basis during this short-time quantum dynamics propagation, such that

    |Φ_α( R_0)⟩≡|Φ_α( R(t_1))⟩,  for t∈[t_1,t_2].


With the above QD basis, the derivative couplings vanish in a trivial way during this short-time propagation, and V̂_el(r̂; R) has off-diagonal elements. We emphasize that there is always a non-removable part of the derivative coupling over the entire configurational space for polyatomic systems.<cit.>  This is a well-known result in literature.<cit.> Here, the QD scheme circumvents this challenge by requiring locally-defined diabatic states, such that the derivative couplings vanish in this configurational subspace during a given short-time propagation.

Because of the diabatic nature of the QD basis during this short-time propagation, one can use any diabatic based approach to propagate the quantum dynamics. These approaches usually require diabatic energies, electronic couplings, and nuclear gradients, which can be conveniently computed<cit.> under the QD basis, {| Φ_α( R_0)⟩}. For example, one can easily evaluate the matrix elements V_αβ( R(t))  = ⟨Φ_α ( R_0)| V̂_el (r̂;  R(t))|  Φ_β( R_0)⟩ at both R(t_1) and R(t_2) as follows 

    V_αβ( R(t_1))    =   ⟨Φ_α ( R_0)| V̂_el (r̂;  R(t_1))| Φ_β( R_0)⟩
    
    V_αβ( R(t_2))    =   ∑_μν  b_αμ⟨Φ_μ( R(t_2))| V̂_el (r̂; R(t_2))|Φ_ν( R(t_2)) ⟩ b^†_βν,

where the first expression is simply equal to E_α( R(t_1))δ_αβ, and the second expression contains the basis transformation matrix elements b_αμ= ⟨Φ_α( R_0)|Φ_μ( R(t_2))⟩ and b^†_βν = ⟨Φ_ν( R(t_2))|Φ_β( R_0)⟩.  The time-dependent matrix elements V_αβ( R(t)) can then be obtained by a linear interpolation between V_αβ( R(t_1)) and V_αβ( R(t_2)) as follows<cit.>

    V_αβ( R(t))= V_αβ(𝐑(t_1))+(t - t_1)/(t_2 - t_1)[V_αβ(𝐑(t_2)) - V_αβ(𝐑(t_1))].
 
Similarly, the nuclear gradients on electronic Hamiltonian matrix elements ∇ V_αβ( R(t_2))≡∂ V_αβ( R(t_2))/∂ R are evaluated as<cit.>

    ∇ V_αβ( R(t_2)) = ∇⟨Φ_α( R_0)| V̂_el(r̂; R(t_2))|  Φ_β ( R_0)⟩
          =⟨Φ_α ( R_0)| ∇V̂_el (r̂; R(t_2))| Φ_β ( R_0)⟩
          =∑_kl  b_αμ⟨Φ_μ( R(t_2))|∇V̂_el (r̂;𝐑(t_2))|Φ_ν( R(t_2))⟩ b^†_βν,
  
where ⟨Φ_μ( R(t_2))|∇V̂_el (r̂;𝐑(t_2))|Φ_ν( R(t_2))⟩ is the nuclear gradient evaluated under the adiabatic basis {|Φ_μ( R(t_2))⟩}, which is readily available from most electronic structure methods.

During the next short-time propagation segment t∈[t_2,t_3], we adapt a new reference geometry R'_0≡ R(t_2) and quasi-diabatic basis |Φ_μ( R'_0)⟩≡|Φ_μ( R(t_2))⟩. With the nuclear geometry close to the reference geometry at every single step, the QD representation remains to be a convenient and compact basis in each short-time propagation segment. Between [t_1,t_2] propagation and [t_2,t_3] propagation segments, all of these quantities will be transformed from {|Φ_α( R_0)⟩} to {|Φ_μ( R'_0)⟩} basis. In particular, we use the following expressions to transform the mapping variables from the previous to the current QD basis between every two consecutive propagation steps

    q_μ   ←   ∑_α q_α⟨Φ_α( R(t_1))| Φ_μ( R(t_2))⟩
    
    p_μ   ←   ∑_α p_α⟨Φ_α( R(t_1))| Φ_μ( R(t_2))⟩.


Here, we use PLDM outlined in Sec. <ref> as the diabatic dynamics method in the QD propagation scheme and refer this approach as QD-PLDM.<cit.> This approach provides the same accuracy for non-adiabatic dynamics as obtained from straight diabatic PLDM, with the additional capability to use adiabatic states and nuclear gradients, obtained from electronic structure calculations, for dynamics propagation.

Previous theoretical work with diabatic quantum dynamics approaches<cit.> for simulating PCET reactions usually require parametrizing the original model system into a strict diabatic system-bath model.<cit.> This parametrization process requires tedious efforts  and remains a highly non-trivial task and significant challenge for atomistic simulations.<cit.> Here, the QD scheme allows directly propagating quantum dynamics by using adiabatic vibronic basis with diabatic dynamics approaches, thus explicitly avoid any additional efforts for building strict diabatic models.

Further, for the PI-PCET dynamics investigated in this study, the QD scheme provides an additional advantage of reducing the number of electronic states required for dynamical propagation. In the case where the diabatic electron-proton basis (see section III) are treated with MMST mapping variables, one can directly use diabatic PLDM for the dynamical propagation due to the strict diabatic nature of this basis. However, the computational cost is significantly increased compared to the QD scheme, due to the large number of the diabatic states (N∼ 100 in this work). Such technical challenges can be resolved by using a more “compact" adiabatic vibronic basis along a given trajectory, i.e., the QD basis, which only requires the few low-lying adiabatic states (N ∼ 10 in this work) which directly participate in the non-adiabatic transitions. 

Thus, the QD scheme provides a convenient propagation framework and a seamless interface for diabatic quantum dynamics approaches with adiabatic electronic structure calculations.<cit.>



§ DETAILS OF MODEL CALCULATIONS



 §.§ Model system
 
The PI-PCET model used in this study is expressed as 

    Ĥ=Ĥ_ep+Ĥ_sb,

where Ĥ_ep describes the electron-proton free-energy surfaces, and Ĥ_sb describes the solvent-bath interaction. In this work, we adapt two commonly used models for Ĥ_ep, with one explicitly contains a collective solvent coordinate (ET coordinate)<cit.> which we refer as Model I, and the other that does not contains a special collective solvent coordinate<cit.> which we refer as Model II. Note that these two models are related with each other through a simple coordinate transformation,<cit.> and are equivalent if the latter one has a brownian spectral density.<cit.>

Here, we provide detailed expressions of Model I, with the parameters for solvent-bath interactions provided in Appendix A. The Hamiltonian of Model II is provided in Appendix B. All of the results presented in this work is based on Model I, except those presented in Fig. <ref> for benchmarking purpose. 

The electron-proton Hamiltonian Ĥ_ep in Eqn. <ref> of Model I is expressed as 

    1.!Ĥ_ep=T̂_p+[
      U^D(r̂_p)+1/2M_sω^2_s R_s^2   V_DA
    V_DA   U^A(r̂_p)+1/2M_sω^2_s (R_s- R_s^0)^2-Δ
    ],

where the first term T̂_p represents the kinetic energy operator of the proton, and the second term represents the electron-proton interaction potential in the diabatic donor |D⟩ and acceptor |A⟩ electronic states, with V_DA=0.03 eV as the coupling between the two electronic states. Here, r̂_p is the proton coordinate, and U^D(r̂_p) and U^A(r̂_p) represent the proton free-energy profile associated with |D⟩ and |A⟩ states, with the following expressions

    U^D(r̂_p)=1/2m_pω_p^2 (r̂_p-r_p^D)^2;  U^A(r̂_p)=1/2m_pω_p^2 (r̂_p-r_p^A)^2.

In this work, we use r_p^D=0 and r_p^A=0.5 Å as the minima of proton free-energy profile associated with the electronic donor and acceptor states. m_p=1.0073 amu and ω_p=3000 cm^-1 are the mass and vibrational frequency of the proton. R_s represents the collective solvent coordinate for electron transfer where M_s and ω_s=√(f_0/M_s) are the mass and frequency of this coordinate with f_0 as the force constant. λ is the solvent reorganization energy, and Δ is the driving force (bias) of the reaction. In this model, the proton and the solvent DOF do not explicitly interact with each other; rather, r̂_p directly interacts with various electronic states, which in turn interact with the solvent. In addition, we choose three possible driving forces,<cit.> with Δ=0 (Model IA), Δ=1 eV (Model IB), and Δ=3.51 eV (Model IC). All the other parameters are provided in Appendix A.



Fig. <ref> illustrates the electron-proton potential Ĥ_ep, with the driving force Δ= 1 eV (Model IB). Fig. <ref>a presents the electronic adiabatic free energy surface of the first excited electronic state S_1, as a function of proton coordinate r̂_p and solvent coordinate R_s. Two dashed lines perpendicular to r̂_p coordinate indicate the minima for proton donor (U^D(r̂_p)) and acceptor (U^A(r_p)) free energy diabats. One dashed line perpendicular to R_s indicate the center of the initial solvent distribution. Fig. <ref>b presents the photoexcited electronic state |D⟩ (blue) and |A⟩ (red), which correspond to the diagonal elements of Ĥ_ep (Eqn. <ref>) evaluated at R^0_s =√(2λ/f_0) (dashed line perpendicular to R_s in panel (a)). Black arrow indicates the photoexcitation process of the proton from its harmonic ground state, U^0=1/2m_pω_p^2r̂_p^2. Fig. <ref>c depicts the adiabatic electron-proton vibronic free energy surfaces as functions of R_s, which are obtained by diagonalizing Ĥ_ep in the vibronic basis (with details described in section <ref>). In this panel, black arrow indicates the photoexcitation that promotes the initial solvent distribution centered at √(2λ/f_0), whereas gray arrows indicate the subsequent vibrational relaxation process.

The solvent-bath Hamiltonian Ĥ_sb in Eqn. <ref> is expressed as follows

    Ĥ_sb=P^2_s/2M_s+∑_k[P_k^2/2M_k + 1/2M_kω_k^2(R_k- c_kR_s/M_kω_k^2)^2].


In the above equation, R_k represents the k^th bath mode, with the corresponding coupling constant c_k and frequency ω_k sampled from the following spectral density 

    J(ω)=π2∑_kc^2_kM_kω_kδ(ω-ω_k)=f_0τ_Lω e^-ωω_c.

Here, τ_L is the solvent response time (see Appendix A), M_k is the mass of the k^th bath mode, and ω_c is the characteristic frequency of the bath that is much faster than the motion of R_s. Here, we choose ω_c=10ω_s and M_k=M_s for all k. 

One can thus perform QD-PLDM simulation with the above total Hamiltonian. Alternatively, we can perform the following equivalent Langevin dynamics<cit.> that treats the bath implicitly, with the equation of motion for the collective solvent coordinate R_s as follows

    M_sR̈_̈s̈= F_ep(R_s)-f_0 τ_LṘ_̇ṡ+ F_r(t).

In the above Langevin equation, F_ep(R_s) is the PLDM force (see Eqn. <ref>) associated with Ĥ_ep, the friction force is -f_0 τ_LṘ_̇ṡ with the friction constant f_0 τ_L, and F_r(t) is the random force bounded by the fluctuation-dissipation theorem through equation ⟨ F_r(t) F_r(0)⟩=2k_BTf_0 τ_Lδ(t).  Here, F_r(t) is modeled as a Gaussian random force with the distribution width<cit.> σ =√(2k_BT f_0 τ_L/dt), where k_B is the Boltzmann constant and dt is the nuclear time step. The details for generating τ_L for a given solvent is provided in Appendix A.

As a consistency check, we have verified that equivalent results (for time-dependent electron-proton reduced density matrix) are obtained with either explicit bath (dynamics with the full Hamiltonian Ĥ_el+Ĥ_sb) or implicit bath (Langevin dynamics in Eqn. <ref>) approach. The equivalency of both approaches have also been recently explored in the condensed-phase ET dynamics<cit.> as well as PI-PCET dynamics.<cit.> We should also note that when the bath DOF has high vibrational frequency such that ħω_s≫ k_BT, linearization approximation<cit.> or classical treatment for the bath can become less accurate.<cit.> Thus, performing implicit Langevin dynamics for the bath can provide more accurate results, especially for those approximate quantum dynamics approaches.<cit.> 



 §.§ Adiabatic Vibronic Surfaces

In this work, we treat both electron and proton quantum mechanically with the corresponding vibronic states. Thus, the “electronic part" of the Hamiltonian, i.e.,  V̂_el in Eqn. <ref> and Eqn. <ref>, is defined as Ĥ_ep in Eqn. <ref>, such that

    V̂_el≡Ĥ_ep(T̂_p,r̂_p, r̂_e,R_s),

Thus,  V̂_el includes proton kinetic energy, electronic potential, electron-proton and electron-solvent interactions. 

In order to obtain the adiabatic vibronic states |Φ_α(R_s)⟩ for the coupled electron-proton Hamiltonian Ĥ_ep, we express |Φ_α(R_s)⟩ with a set of two-particle basis functions as follows

    |Φ_α(R_s)⟩=∑_i,m c_im^α(R_s)|ϕ_e^i⟩|ϕ_p^m⟩,

where  |ϕ_e^i⟩∈{|D⟩,|A⟩} and |ϕ_p^m⟩ is chosen to be the m^th eigenfunction of a quantum harmonic oscillator, with the total Hamiltonian Ĥ=T̂_p+1 2m_pω_ p^2r̂_p^2. Thus, by using M harmonic basis functions for proton and two basis states for electron, the total number of vibronic basis is N=2M, and Ĥ_ep contains 2M× 2M Hamiltonian matrix ⟨ϕ_p^n|⟨ϕ_e^j|Ĥ_ep|ϕ_e^i⟩|ϕ_p^m⟩ under this representation. Because both U^D(r̂_p) and U^A(r̂_p) are just simple displaced harmonic oscillator potentials, the matrix elements of Ĥ_ep can be obtained analytically by recognizing the basic property of harmonic oscillator as follows 

    ⟨ϕ_p^n|T̂_p+1 2m_pω_p^2r̂_p^2|ϕ_p^m⟩=(n+12)ħω_pδ_nm
          ⟨ϕ_p^n|r̂_p|ϕ_p^m⟩ = √(ħ/m_pω_p)1/√(2)(√(m) δ_n,m-1+√(m+1) δ_n,m+1).


The eigenvalues and the eigenvectors (adiabatic vibronic basis) are then obtained through direct diagonalization of the Ĥ_ep matrix under the above two-particle basis. 



 §.§ Quantum dynamics propagation approaches
 
We perform PLDM quantum dynamics simulations with three different choices for treating the transferring electron and proton, which are summarized as follows 


  * QD-PLDM: describe the electron-proton adiabatic vibronic basis |Φ_α(R_s)⟩ as the quasi-diabatic (QD) states with MMST mapping variables; propagate the dynamics with QD-PLDM approach,<cit.>

  * vib-PLDM: describe the electron-proton diabatic vibronic (vib) basis |ϕ_e^i⟩|ϕ_p^m⟩ with MMST mapping variables; propagate the dynamics with straight diabatic PLDM approach,

  * el-PLDM: only describe the electronic (el) state  |ϕ_e^i⟩ with MMST mapping variables, whereas the proton is treated through linearization approximation (which gives rise to classical equation of motion and a Wigner initial distribution); propagate the dynamics with straight diabatic PLDM approach. 


Here, we briefly comment on the numerical cost of each approach. The convergence of results in vib-PLDM and QD-PLDM is obtained with a considerably less number of trajectories (∼ 10 times fewer) than in el-PLDM. This is because el-PLDM requires sampling of the Wigner distribution of both the solvent and the proton DOF (with a broad width of the distribution due to large ω_p), while in vib-PLDM and QD-PLDM the initial sampling is only required for the solvent DOF. 

In addition, the number of states required to be explicitly propagated are different in each one of these schemes. Here, el-PLDM only requires explicit propagation of two states, |D⟩ and |A⟩;  vib-PLDM requires explicit propagation of 60-120 diabatic vibronic basis {|ϕ_e^i⟩|ϕ_p^m⟩}; QD-PLDM only requires 5-20 adiabatic vibronic states {|Φ_α(R_s)⟩}, due to the compactness of these adiabatic states for describing the changing wavefunction. Recall that the numerical cost of PLDM propagation scales as N^2 (with N as the total number of states). Thus, the ultimate numerical cost for vib-PLDM is much larger than both QD-PLDM and el-PLDM. QD-PLDM, on the other hand, still requires explicitly diagonalizing Ĥ_ep matrix with a numerical cost of N^3, besides the cost for dynamical propagation.



 §.§ Simulation details

Here, we provide the simulation details for Model I, whereas the corresponding details for Model II are provided in Appendix B. 

The converged results for Model I are obtained with 2400 trajectories for QD-PLDM or vib-PLDM propagations, with a time step of dt=0.024 fs (1 a.u.). The total number of vibrational basis {|ϕ_p^m⟩} used in Model IA is 30, i.e., m=0,1,...29 for the Harmonic oscillator eigenstates. The total number of vibrational basis used in Model IB and IC are 40 and 60, respectively. We have carefully checked the convergence of our results with additional 10 vibrational basis, which generates numerically identical results. For QD-PLDM propagation, we only used the first 5, 10, and 20 time-dependent low-lying adiabatic vibronic states as the QD states for the electron-proton description. We have also carefully checked the convergence of the QD propagation scheme with additional 10 more QD basis, which also provides the identical results.

In all calculations with Model I, the system is initially prepared in the proton vibrational ground state |ϕ_p^0⟩ of the electronic ground state |S_0⟩. The system is then excited to the |D⟩ state through Franck-Condon process, which generate the initial state described by the following total density operator 

    ρ̂(0)=|Φ(0)⟩⟨Φ(0)|⊗ρ̂_s.

Here, the initial electron-proton quantum state is expressed as 

    |Φ(0)⟩=|D⟩|ϕ_p^0⟩,

and ρ̂_s is the density operator of the solvent. PLDM requires the partial Wigner transform of the total density operator ρ̂(0) which can be easily obtained (due to the its simple direct product form) as follows

    [ρ̂(0)^W]=|Φ(0)⟩⟨Φ(0)|⊗ρ_s^W,
 
with the following Wigner density for the solvent 

    ρ_s^W = ω_sΓ_s e^-Γ_s[P_s^22M_s + 12M_sω_s^2(R_s-R_s^0)^2].

Here,  Γ_s= (2/ω_s)tanh(ω_s/2k_BT) and ω_s = √(f_0/M_s). In this study, we choose R^0_s =√(2λ/f_0) that corresponds to the minimum of the proton acceptor free energy diabatic surface. Further, we use the focused initial conditions<cit.> to facilitate the convergence of the sampling for the mapping variables, which obey the distribution governed by G_0( p,q) and G'_0 ( p',q') (see Sec.<ref> for details). In vib-PLDM, this means that q_ξ=q'_ξ=δ_ξη and p_ξ=-p'_ξ=δ_ξη where |η⟩=|Φ(0)⟩=|D⟩|ϕ_p^0⟩ and |ξ⟩=|ϕ^i_e⟩|ϕ^m_p⟩. Whereas in QD-PLDM, the initial values of the corresponding mapping variables are obtained through the following expressions

    q_α=∑_im q_ξ c^α_im;  p_α=∑_im p_ξ c^α_im,

where c^α_im is the eigenfunction coefficients obtained from Eqn. <ref> and |ξ⟩=|ϕ^i_e⟩|ϕ^m_p⟩.



§ RESULTS AND DISCUSSIONS

Fig. <ref> presents the diabatic electronic  population of the |D⟩ state obtained from el-PLDM (dotted line), vib-PLDM (open circle), and QD-PLDM (solid line) for Model II. In addition, in order to assess the accuracy of these approaches, we present results obtained from reduced density matrix (RDM) formalism<cit.> which provides the exact results (black dashed line) under the weak system-bath coupling regime for Model II. It can be clearly seen from Fig. <ref>a that el-PLDM which propagates the motion of proton classically (from its initial Wigner distribution), does not provide accurate electronic dynamics and causes large deviation from the exact result at a longer time. Similar deviation has also been reported by using symmetrical quasi-classical (SQC) method<cit.> with the same classical treatment for proton. On the other hand, quantizing the proton with vibronic basis in both vib-PLDM and QD-PLDM approaches leads to accurate PI-PCET dynamics. 

We would like to mention that treating vibrational diabatic states with MMST mapping variables can certainly improve the accuracy of the dynamics, as been demonstrated with other recent theoretical works. In one example, vib-PLDM is used to compute 2-dimensional electronic spectra<cit.> in a Frenkel excitons model with a high-frequency vibrational mode. in another example, a new method, so called extended SQC approach<cit.> is developed based on similar strategy of vib-PLDM, which provides accurate non-adiabatic dynamics when a highly non-harmonic mode is explicitly quantized with its vibrational eigenbasis. In the third example, the lowest four electron-proton diabatic states are explicitly used in PCET quantum dynamics studies.<cit.> However, we want to emphasize that in general, vib-PLDM does require a large number of strict diabatic states to be explicitly propagated for investigating PI-PCET reaction. Further, strict diabatic states cannot be easily obtained for real systems, in addition to the computational disadvantage associated with propagating a large number of states. This feature will ultimately limit the scope and applicability of vib-PLDM for investigating PI-PCET. The QD scheme, on the other hand, only requires a small set of adiabatic vibronic states for time-dependent propagation, and thus provides an accurate and efficient theoretical framework for investigating PI-PCET dynamics.



Fig. <ref>b presents the kinetic isotope effect (KIE) for proton (red) and deuterium (blue) in terms of the donor population. For clarity, here we only present the results obtained from el-PLDM (dotted lines) and QD-PLDM (solid lines), whereas the results from the latter one are identical to vib-PLDM, and agree very well with exact results obtained from RDM<cit.> (not shown here). Because deuterium is essentially a classical particle, the el-PLDM and QD-PLDM approaches provide similar results (blue curves) despite some small deviations. This suggests that treating deuterium classically provide a reasonably accurate dynamics. whereas the quantum nature of the proton requires an explicit propagation with electron-proton vibronic states for accurate results.





Fig. <ref> presents the adiabatic vibronic population of the PI-PCET dynamics, with the adiabatic surfaces of models IA-IC provided in (a)-(c). The initial photoexcitation is illustrated with black arrows, and the subsequent vibrational relaxation pathways are indicated with gray arrows. The corresponding color-coded vibronic state populations are presented in panels (d)-(e), calculated using vib-PLDM (open circles) and QD-PLDM (solid lines) propagations, which are identical. Quantitatively similar results for these vibrational relaxation dynamics have also been obtained from FSSH<cit.> approach (results not shown), suggesting that the decoherence correction might not have large impact on the short-time vibronic dynamics (∼ 1 ps). 

In Fig. <ref>a-c, the initial photoexcitation leads to populating a set of high-lying vibronic excited states, followed by vibrational relaxation process that propagates the vibronic wavepacket into low-lying states through non-adiabatic transitions. The grey arrows indicate these vibrational relaxation dynamics during the first 1 ps timescale of the simulations. Compared to the symmetric case in Model IA (Δ=0), the non-zero energy bias in Models IB and IC leads to populating much higher vibronic states at the beginning of the reaction, as well as stabilizing the acceptor state over the donor state that impacts the longer time dynamics. For model IC, the solvent coordinates R_s directly relaxes to the electronic acceptor side during the first 1 ps (results not shown), whereas in Model IA and IB, R_s relaxes back to the electronic donor side, suggesting a much slower ET dynamics associated with IA and IB, as will be demonstrated in the next figure. 

Fig. <ref> presents the KIE with electronic population dynamics of Model IB and IC. The corresponding dynamics in Model IA is much slower compared to IB and IC, and thus is not shown here. Fig. <ref>a depicts the population decay of the donor state for proton (red) and deuterium (blue) in Model IB. It can be seen that there is a clear separation of time scale during the PI-PCET dynamics, with an initial fast vibrational relaxation process during the first 1 ps (that corresponds to results shown in Fig. <ref>d), followed by a second stage, much slower non-adiabatic dynamics that transfer electronic population from the donor to the acceptor state. In contrast, in the system with deuterium, the donor electronic population does not significantly transfer during the same time scales. These results indicate a large KIE that can be observed in Model IB.

Fig. <ref>b presents the same population dynamics for Model IC. One can also observe a similar two-stage dynamical process, with an ultrafast sub-picosecond relaxation process and a relatively slower (∼ 5 ps) charge population transfer dynamics. In this model system, KIE is negligible compared to the previous model, whereas the early stage relaxation process for deuterium is even faster than proton. Similar negligible or even slightly inverse KIE<cit.> has also been observed through recent theoretical investigations. The inverse KIE can be easily understood as follows. When tunneling effects are less important, the vibrational relaxation dictates the dynamics,<cit.> and with a larger nuclear mass, deuterium relaxes even faster than proton<cit.> because the vibrational states are closer in energy. Quantum mechanically, the vibrational gap of deuterium is much smaller compared to proton. With the same initial photoexcitation, more high-lying excited vibronic states can be populated for deuterium, and thus promotes the PCET process.<cit.> Further, as been previously discussed,<cit.> the lack of KIE in the initial stage of PI-PCET cannot exclude the possibility of the concerted transfer of both proton and electron. 

In order to understand the distinctly different KIE in the above two model systems, we compute the time-dependent probability density of the transferring proton/deuterium associated with the donor electronic states |Φ_D(r)|^2=∑_ηξρ_ηξ(t) ⟨ r|ϕ_p^m⟩⟨ϕ_p^n|r⟩, where ρ_ηξ(t) is the reduced density matrix in the electron-proton diabatic vibronic basis {|ϕ^i_e⟩|ϕ^m_p⟩,|ϕ^j_e⟩|ϕ^n_p⟩}, |η⟩=|D⟩|ϕ_p^m⟩, and |ξ⟩=|D⟩|ϕ_p^n⟩. Similar expression is used for computing |Φ_A(r)|^2, where |η⟩=|A⟩|ϕ_p^m⟩, and |ξ⟩=|A⟩|ϕ_p^n⟩.

Fig. <ref>c-e presents these time-dependent probability densities for the transferring proton/deuterium described above. In Fig. <ref>c, there is a significant transfer of the proton probability distribution from the donor to the acceptor states in Model IB, whereas there is no transfer for deuterium probability distribution presented in Fig. <ref>e. This different behavior suggests that tunneling of the proton between the donor and acceptor states dominate the PI-PCET dynamics at a longer time, after the initial short-time vibrational relaxation. In contrast, for Model IC, the probability densities for both proton (Fig. <ref>d) and deuterium (Fig. <ref>f) exhibit very similar time-dependent behavior, suggesting a predominant role of vibrational relaxation and less important role of tunneling in the PI-PCET dynamics.


To demonstrate the effect of proton tunneling on PI-PCET dynamics in Model IB, in Fig. <ref> we present the proton vibrational eigenfunctions Ψ_α(r_p;R_s)=⟨r_p|Φ_α(r̂_p;R_s)⟩ with different solvent configurations R_s. In Fig. <ref>a, the solvent configuration R_s is chosen to be at its initial value R_s^0=√(2λ/f_0) upon photoexcitation. The initial proton wavepacket is created as a linear combination of vibational states 5-7, which are blue, red and brown states in Fig. <ref>a. During the first 1 ps, solvent coordinate R_s has relaxed to the donor side (with R_s∼0), with the corresponding vibronic eigenfunctions presented in Fig. <ref>c. At the same time, vibrational relaxation process has induced the population transfer to vibronic states 3-4 (depicted as the pink and violet states in Fig. <ref>c). Note that there is a sizable barrier along the r̂_p coordinate for the relaxed solvent configuration (R_s∼0) in panel (c) compared to panel (a) (R_s^0=√(2λ/f_0)). Under such circumstances in (c), the populated proton vibrational states (pink and violet) are below the barrier, and the tunneling through the barrier becomes the predominant mechanism for the later stage PI-PCET dynamics. This explains the significant KIE observed in Model IB (Fig. <ref>a). In contrast, PI-CPET dynamics in Model IC is dictated by the early stage vibrational relaxation process, with the adiabatic vibronic states that are more similar to Fig. <ref>a, such that it is almost barrier-less for the reaction (akin to the activation less regime in Marcus electron transfer theory), resulting in a negligible KIE. 




§ CONCLUSIONS

We apply the recently developed QD propagation scheme<cit.> to investigate the non-adiabatic dynamics of photoinduced proton coupled electron transfer (PI-PCET) reaction. Using PLDM path-integral approach and the electron-proton adiabatic vibronic states as the time-dependent quasi-diabatic states, the outlined QD propagation scheme provides an accurate and efficient theoretical framework for simulating PI-PCET dynamics. Compared to approaches that only treat the transferring electron quantum mechanically but proton classically (such as classical Wigner models<cit.>), the QD scheme explicitly quantizes proton with vibronic adiabatic states, and thus provides accurate non-adiabatic dynamics and KIE. Compared to the approaches that directly propagate dynamics with the diabatic vibronic basis (such as extended SQC<cit.>), the QD scheme only requires a smaller set of adiabatic vibronic states that are directly involved during the non-adiabatic process, thus significantly reduces the computational costs associated with the number of states that need to be explicitly propagated. Previous theoretical studies with diabatic quantum dynamics approaches<cit.> for simulating PCET reactions usually require parametrizing the original model system into a strict diabatic system-bath model.<cit.> This process requires tedious efforts and remains a non-trivial task and a significant challenge for atomistic simulations.<cit.> Here, by using adiabatic vibronic basis that can be obtained with routinely available electronic structure calculations,<cit.> the QD scheme allows directly propagating quantum dynamics with a diabatic based approach. 

With this QD propagation scheme, we investigate the vibronic population transfer and the KIE of PI-PCET dynamics with various driving forces. For systems with a small driving force, while the vibrational relaxation process significantly impacts the early stage dynamics, solvent relaxation to the donor side of the free energy surface eventually creates a large barrier for PCET, such that proton tunnelling plays a predominant role. A significant KIE will be observed in such scenario. Whereas for a system with a much larger driving force, the vibrational relaxation completely dictates both early stage dynamics as well as longer-time PI-CPET reactions, resulting in a negligible KIE.

Further, we want to emphasize that the outlined QD scheme and the simulation protocol are general enough and not limited to PLDM; they can be directly applied to a wide range of diabatic trajectory-based quantum dynamics approaches. These approaches include but are not limited to symmetrical quasi-classical (SQC),<cit.> forward-backward quantum-classical Liouville equations (FB-QCLE),<cit.> generalized quantum master equation (GQME),<cit.> and quantum-classical path-integral (QCPI) dynamics.<cit.> The QD propagation scheme provides a transformative theoretical framework for studying challenging PI-PCET reactions through accurate diabatic quantum dynamics approaches with efficient adiabatic electronic structure calculations.

Finally, we want to outline three alternative approaches for quantizing proton besides the commonly used vibronic adiabatic states description<cit.> adapted in this study. The first one requires using time-dependent Gaussian basis (TDGB) function to explicitly expand the proton (as well as the other nuclear) wavefunction. This approach that has been utilized in AIMS<cit.> or MP/SOFT,<cit.> can be numerically expensive as it requires many TDGB functions associated with each nuclear DOF. The second one uses nuclear-electronic orbital (NEO) approach<cit.> in which wavefunctions are used for the transferring electron and proton with molecular orbital techniques, adding additional difficulty on top of the already complex electronic structure problems. The last one quantizes the proton with an imaginary-time path-integral framework<cit.> in the extended classical phase-space (i.e., so-called ring polymer). With an explicit description of electronic states, recently emerged state-dependent ring polymer molecular dynamics approaches<cit.> can potentially provide accurate electronic non-adiabatic dynamics with nuclear quantum effects, and thus are promising for investing PI-PCET dynamics once combined with the QD-propagation scheme.



§ ACKNOWLEDGEMENT

This work was supported by the University of Rochester startup funds. Computing resources were provided by the Center for Integrated Research Computing (CIRC) at the University of Rochester. FAS appreciates valuable discussions with Prof. Anirban Hazra.



§ APPENDIX A: SOLVENT PARAMETERS FOR MODEL I

We provide the details of the parameters used in Model I. The force constant for the collective solvent DOF (so-called “inverse Pekar factor") is f_0=4πϵ_0 ϵ_∞/(ϵ_0 - ϵ_∞), where ϵ_0 and ϵ_∞ are the inertial and optical dielectric constants characterizing the polarizability of the solvent. Here, we chose these parameters that correspond to water as the solvent.<cit.> 



Further, τ_L=ϵ_∞(τ_0+τ_D)/ϵ_0 is the longitudinal relaxation time accounting for the long-time solvent response function, where τ_D is the Debye relaxation time and τ_0 is the characteristic rotational time of the solvent molecules. All of the parameters used in this paper are tabulated in Table I and a full description of them could be found in Ref. .



§ APPENDIX B: HAMILTONIAN FOR MODEL II

Here we provide the details for the PI-PCET model that does not contain collective solvent coordinate,<cit.> which is referred to as Model II in this paper. This Hamiltonian is used to explore the accuracy of various recently developed non-adiabatic approaches as well as the role of proton quantization for PI-PCET reactions in this study. Note that only the results presented in Fig. <ref> are obtained with this model system. 

The total Hamiltonian is defined in Eqn. <ref>. The electron-proton Hamiltonian Ĥ_ep is expressed as following

    Ĥ_ep=T̂_p+
     [
      U^D(r_p)   V_DA
    V_DA   U^A(r_p)
     ]
 
Here r_p is the proton coordinate, U^D(r_p) and U^A(r_p) are the proton potential associated with electronic donor and acceptor states that have exactly the same expression as Model I (see Eqn. <ref>), with r_p^D=0 and r_p^A=-0.5 Å. Prior to photoexcitation, proton is on the vibrational ground state of the S_0 state, with the potential U^0=1/2m_pω_p^2 (r_p-r^0_p)^2 where r^0_p=-0.15 Å. The rest of parameters in the above Hamiltonian have the same values as used in Model IA, with Δ=0 eV, V_DA=0.03 eV, ω_p=3000 cm^-1, and m_p=1.0073 amu.

The bath Hamiltonian which describes the interaction between the electron-proton system and a condensed-phase solvent environment is modeled by coupling of the donor electronic state to a dissipative harmonic bath with the following expression

    H_sb= ∑_k=1^K[ P_k^2/2M_k + 1/2M_kω_k^2(R_k- c_k/M_kω_k^2|D⟩⟨ D|)^2],

where R_k and P_k represent the k^th bath coordinate and momentum, with M_k and ω_k as the corresponding mass and frequency. The bath is characterized by an Ohmic spectral density J(ω) =  1/2πξω e^-ω/ω_c, where  ξ is the unit-less Kondo parameter and ω_c is the cut-off frequency. Here, we use ξ=24 and ω_c=600 cm^-1. Discretizing this spectral density yields N harmonic oscillators with frequencies ω_k = -ω_c ln(1-kω_0/ω_c) and coupling constants, c_k= √(ξω_0M_k)ω_k. Here, ω_0 for a total of K bath modes is given by ω_0 = ω_c/K(1- e^-ω_m/ω_c),  and ω_m was chosen to be 3ω_c.

The initial condition for the PLDM simulation are provided as follows. The initial conditions for the bath modes are sampled from the Wigner distribution for harmonic oscillators' thermal density as follows

    ρ_b^W = Π_k=1^Kω_k Γ_k e^-Γ_k[P_k^22M_k + 12M_kω_k^2(R_k-R_k^0)^2],

where Γ_k = (2/ω_k)tanh(ω_k/2k_BT), and ω_k sampled from the spectral density, and R_k^0=c_k/(M_kω_k^2). Further, we use the focused initial conditions<cit.> to facilitate the convergence of the sampling for the mapping variable.  

For the el-PLDM calculation, we choose to treat proton classically, with the corresponding initial conditions sampled from the following function

    ρ_p^W = ω_pΓ_p e^-Γ_p[T_p + 12m_pω_p^2(r_p-r_p^0)^2],

where Γ_p = (2/ω_p)tanh(ω_p/2k_BT), T_p is the classical kinetic energy of the proton, r_p is the proton coordinate, ω_p is the proton vibrational frequency.

The converged results for Model II with el-PLDM method are obtained by propagating an ensemble of 10^4 trajectories, with a time step of dt=0.024 fs. For the same model, we use 2000 trajectories for vib-PLDM or QD-PLDM propagation. The total number of vibrational basis {|ϕ_p^m⟩} used in this model is 80. In the QD-PLDM propagation, we use the first 20 low-lying time-dependent vibronic states as the QD basis. 


