Room Temperature Quantum Coherence vs. Electron Transfer in a Rhodanine Derivative Chromophore
Duvalier Madrid-Úsuga Cristian E. Susa John H. Reina




duvalier.madrid@correounivalle.edu.co
Centre for Bioinformatics and Photonics (CIBioFi), Calle 13 No. 100-00, Edificio E20 No. 1069, Universidad del Valle, 760032 Cali, Colombia
Departamento de Física, Universidad del Valle, 760032 Cali, Colombia

cristiansusa@correo.unicordoba.edu.co
Centre for Bioinformatics and Photonics (CIBioFi), Calle 13 No. 100-00, Edificio E20 No. 1069, Universidad del Valle, 760032 Cali, Colombia
Departamento de Física y Electrónica, Universidad de Córdoba, 
230002 Montería, Colombia

john.reina@correounivalle.edu.co
Centre for Bioinformatics and Photonics (CIBioFi), Calle 13 No. 100-00, Edificio E20 No. 1069, Universidad del Valle, 760032 Cali, Colombia
Departamento de Física, Universidad del Valle, 760032 Cali, Colombia


  

Understanding electron transfer in organic molecules is of great interest in quantum materials for light harvesting, energy conversion and integration of molecules into solar cells. This, however, poses  the challenge of designing specific optimal molecular structure for which the processes of ultrafast quantum coherence and electron transport are not so well understood. In this work, we investigate subpicosecond time scale quantum dynamics and electron transfer in an efficient electron acceptor Rhodanine chromophoric complex. We consider an open quantum system approach to model the complex-solvent interaction,  and compute the crossover from weak to strong dissipation on the reduced system dynamics for both a polar (Methanol) and a non polar solvent (Toluene). We show that the electron transfer rates are enhanced in the strong chromophore-solvent coupling regime, being the highest transfer rates those found at room temperature. Even though the computed dynamics are highly non-Markovian, and they
 may exhibit a quantum character up to hundreds of femtoseconds, we  show that quantum coherence does not necessarily optimise the electron transfer in the chromophore.


03.65.Yz, 
87.15.ht,       
82.20.Xr.       
87.15.ag,      
04.25.-g       




Room Temperature Quantum Coherence vs. Electron Transfer in a Rhodanine Derivative Chromophore
    John H. Reina
    
==============================================================================================





§ INTRODUCTION


Electron transfer (ET) mechanisms and dynamics are ubiquitous to fundamental processing of energy in biological and chemical systems that are essential to life. They are at the core of light harvesting for energy production in a variety of systems; e.g., biomolecular photosynthetic complexes, biosensors, and solar cells <cit.> to name but a few. From a perspective of novel materials and technologies for solar energy conversion, organic photovoltaic (OPV) sytems have attracted significant attention <cit.> due to the  enhancement of their power conversion efficiency (PCE) <cit.>, flexibility,manufacturing reduction costs, and high-temperature production of cells <cit.>.
On a different facet, but complementary to this development, quantum coherence has, over the past decade, been experimentally linked to early-stage energy transfer in light harvesting photosynthetic complexes <cit.>;  in particular, time-resolved ultrafast optical and electronic spectroscopies have  evidenced oscillations of exciton state populations lasting up to a few hundred femtoseconds, which have been attributed to quantum coherence emerging as a result of initially prepared laser pulses <cit.>. This two-way breakthrough poses an intriguing question about the possible functional role of quantum coherent effects as a precursor of efficient ET in materials and biological systems under structural and energetic disorder at physiological conditions.

The above physical insight requires that ET reactions in such chemical or biological light harvesters must consider the effects due to electronic coupling, interactions with the environment, and intramolecular vibrational modes.In particular, photoinduced ultrafast energy and ET processes in OPV  devices require a quantum-mechanical treatment of such fundamental microscopic processes. Here, electron transfer occurs between electron-donor and electron-acceptor molecules immersed in a solvent, and such process can be modelled by considering a quantum system in interaction with a local environment (solvent), whose dynamics, control, and correlations are amenable to a description within  the framework of open quantum systems <cit.>. Donor-acceptor chromophoric complexes convey electronic coupling between sites on a scale that relies on the solvent reorganisation energy and defines different regimes for the complex coherent/incoherent dynamics. Due to the interplay between energy fluctuation and transport, here we account for the crossover from weak to strong dissipation and consider memory effects for the chromophore-solvent quantum dynamics <cit.>; these can be described by means of intensive computational techniques such as path-integral based formalisms <cit.>, Monte Carlo algorithms <cit.>, Hierarchical Equations of Motion (HEOM) <cit.>, and reaction-coordinate methods <cit.>.

In this work, we consider an OPV light harvester (an efficient Rhodanine derivative chromophore <cit.>),  and focus on theoretical aspects regarding its ultrafast carrier dynamics and the interplay between the solvent, electronic coupling, and temperature conditions on both quantum coherence and electronic transfer. We consider Methanol, the polar solvent used in <cit.>,  and, for completeness, we also quantify the effects due to a non polar solvent, Toluene. In both cases, the chromophore-solvent coupling is found to be in the strong regime and the corresponding non trivial non-Markovian dynamics is treated  within a non-perturbative open quantum system approach, with the numerically exact  HEOM method <cit.>, on a spin-boson model of the complex-solvent <cit.>.

This paper is organised as follows: the chromophore and solvent  properties are developed in Sec. <ref>; the chemical structure and main parameters of the organic molecule are presented in Sec. <ref>, and the mathematical model and chromophore+solvent physical parameters are given in Sec. <ref>. The main results and discussions are presented in Sec. <ref>: solvent and temperature effects on the quantum dynamics (Sec. <ref>), and electronic transfer (Sec. <ref>). Finally, conclusions are given in Sec. <ref>.




§ CHROMOPHORE+SOLVENT SYSTEM: CHEMICAL AND PHYSICAL PROPERTIES


The considered push-pull chromophore derived from Rhodanine, the Dicyanorhodanine+2-formyltetratiafulvalene (D2F) complex <cit.>, has a  chemical structure as shown in Fig. <ref>(a): it consists of a 2-formyltetratiafulvalene as the electron-donor (D) and a Dicyanorhodanine as the electron-acceptor (A) linked by Propylene. The compound exhibits appealing photo- and electron-chemical properties; a D-A intramolecular charge transfer in the absorption band of the donor at λ_max=498 nm (visible region). This was confirmed by solvatochromics studies using Acetone, Ethanol, and Methanol as solvents <cit.>: the “push-pull”  compound displays a clear electrochemically amphoteric behaviour and  reveals a significant D-A electronic communication through the π-conjugated core, a feature that can be useful for ET integration into molecular systems of greater complexity.





 §.§ D2F Chemical Structure
 

The geometry of the D2F molecule was optimised by means of Gaussian 09 <cit.>. The hybrid functional B3LYP and the 6-31G^++ basis are used to estimate the energy values of D and A electronic states, as well as the coupling between them. We consider a scenario where the system is firstly photoexcited from the ground electronic state to a π-π^* state localised on D. This is followed by a nonradiative process that corresponds to an electronic transition from the π-π^* state to a charge-transfer (CT) state that involves a significant D-A charge transfer. The molecular orbital analysis of both ground and excited states of the optimised geometry indicates that the HOMO (LUMO) is mainly located on D (A). The two aforementioned states are illustrated in Fig. <ref>(b); the first transition that gives the excited state (EX) corresponds to the photoexcitation D–A → D^*–A at an excitation energy of 2.85 eV. This exhibits an oscillator strength f_osc=0.2348 and a percentage contribution of 83% (HOMO → LUMO+5 transition in Fig. 1(b)). The second one that generates the CT state (D^-–A^+) has an energy of 2.73 eV, and this state has f_osc=0.0002, with a percentage contribution of 75% for the HOMO → LUMO+4 transition (Fig. 1(b)). These values were computed for Methanol via the time dependent density functional theory, by using the B3LYP functional. This scenario corresponds to an ET process from the 2-formyltetratiafulvalene compound to the Dicyanorhodanine  acceptor.

The energies associated to the two main states involved in the ET mechanism are sketched in Fig. <ref>. According to the B3LYP/6-31G^++ optimisation, the estimated site energies for Methanol are E_D=-4323 cm^-1 and E_A=-5291 cm^-1.  On the other hand, for Toluene E_D=-5201 cm^-1 and E_A=-5838 cm^-1.  Other physical parameters of interest are given in Fig. <ref> and Table <ref>.
















 §.§ Hamiltonian and  Chromophore Quantum Dissipative Dynamics
 

We consider a spin-boson model to treat the chromophore-solvent ultrafast dynamics in a quantum fashion <cit.>. The total Hamiltonian can be written as H = H_S + H_B + H_SB, where the subscripts stand for system (the chromophore), bath (the solvent), and for the system-bath interaction, respectively. The solvent is modelled as a bath of harmonic oscillators with H_B = ∑_jħω_j( a_j^†a_j + 1/2), where a_j^† (a_j) denotes the creation (annihilation) operator of the j-th bosonic mode with frequency ω_j. 

The N-partite donor-bridge-acceptor chromophore Hamiltonian is described in terms of the local (diabatic) electronic states |φ_n⟩ at site n:


    H_S = ∑_n=1^NE_n|φ_n⟩⟨φ_n| + ∑_n≠ m =1^NV_nm|φ_n⟩⟨φ_m|,


where E_n corresponds to the energy of the n-th site and V_nm is the coupling energy between sites n and m. Here, we consider a reduced two-site model (N=2) between the EX and CT states with energies as shown in Fig. <ref>, following the prominent D-A charge transfer due to the π-π^* state described in Sec. <ref>. Thus, an effective two-level (D-A) system Hamiltonian can be identified for H_S:


    H_TLS=ϵ/2σ_z+V_DAσ_x,


where V_DA=ϵμ_12/[(μ_1-μ_2)^2+4μ_12^2]^1/2
is the electronic coupling, here computed by means of the generalised Mulliken-Hush model <cit.>, ϵ≡|E_D-E_A| is the energy difference between donor and acceptor orbitals,  and σ_x,z are Pauli's operators. For Methanol, μ_12=6.6990 D is the transition dipole moment connecting the two adiabatic states in the charge transition, and Δμ_12= μ_1-μ_2=18.8447 D is the difference in adiabatic state dipole moments.  A similar procedure gives the corresponding values for Toluene.

The solvent is modelled by considering that each mode interacts independently with the electronic sites. The explicit interaction in the simplified two-site scenario, where g_j accounts for the coupling strength to the j-th mode, is written as:


    H_SB = ħσ_z⊗∑_jg_j(a_j^† + a_j).


 The information about the solvent and its interaction with the chromophore is encoded in the spectral density J(ω)=∑_jg_j^2δ(ω-ω_j) and the bath correlation functions (see Appendix <ref> for details). A precise description of this information through the spectral density is usually not a straightforward task. Here, we consider Onsager's reaction model for describing the D2F-solvent coupling. Geometrically, the D2F chromophore is assumed to occupy a region of radius a_0 (∼ the van de Waals size) surrounded by the solvent <cit.>. Within this model, a spectral density with Lorentzian shape can be derived <cit.>:


    J(ω)=(Δμ)^2/2πε_0a_0^36(ε_s-ε_∞)/(2ε_s+1)(2ε_∞+1)ωτ_s/ω^2τ_s^2+1,


where Δμ=μ_e-μ_g is the relative dipole moment between the excited and ground state dipole moments of the D-A molecule. ε_0, ε_s and ε_∞ stand for vacuum, solvent static and high frequency dielectric constants, respectively. τ_s is the solvent relaxation time (cutoff frequency ω_c≡1/τ_s). The reorganisation energy λ:=2/π∫J(ω)/ωdω is simply given by λ=2ħαω_c, where α = (Δμ)^2/4πϵ_0a_0^36(ϵ_s-ϵ_∞)/(2ϵ_s+1)(2ϵ_∞+1)ω_c is the frequency-independent factor of Eq. (<ref>).

The quantum dissipative dynamics of the chromophore can be described within the density operator formalism <cit.>. If ρ_T(t) is the time-dependent density operator associated to the total (chromophore+solvent) system (considered to be isolated), the reduced density operator for the chromophore is given by ρ_S(t)=Tr_B(ρ_T(t)). Here, we employ the HEOM technique <cit.> to numerically compute the temporal evolution of ρ_S(t), as explained in Appendix <ref>. 

The coherent nature of the ET depends on several factors such as temperature and chromophore-solvent coupling. Here, we consider two different solvents; Methanol (polar) and Toluene (non polar) and present numerical results for their effects on the D2F dynamics and the electronic transfer. Important physical parameters that characterise the D2F system are given in Table <ref>.



§ RESULTS AND DISCUSSION
 

The spectral density (<ref>) can be written in a simpler way as J(ω)=2αω/(ω^2τ_s^2+1), where the  α values are determined from the parameters in Table <ref>, for each solvent. In fact, this chromophore-solvent coupling gives reorganisation energies of 86.18 cm^-1 and 2.20 cm^-1 for Methanol and Toluene, respectively. In order to directly compare the effects of both solvents on the D2F dynamics, we introduce a dimensionless system-bath coupling:


    η := α/ħ=λ/2ħω_c.


A system-bath coupling falls within the weak regime if η≪1; in the strong regime η∼1 <cit.>. We obtained, according to the solvent properties given in Table <ref>,  η=0.25 and 0.33 for Methanol and Toluene, respectively: these indicate that the strong coupling regime dictates the chromophore dynamics and ET in both cases. 

 




 §.§ Solvent and Temperature Effects
 

Having determined the main properties of the considered chromophore-solvent system, we investigate the solvent effects on the D2F dynamics in both weak and strong regimes. For doing so, we vary the value of η and compute the population inversion ⟨σ_z(t)⟩=ρ_D(t)-ρ_A(t), as well as the real-part coherence (ρ_DA(t)) associated to the chromophore, where ρ_D, ρ_A (populations) and ρ_DA=ρ_AD^* (coherence) are the matrix elements of the density operator ρ_S(t). 
In all the figures, populations are considered to be localised on the donor at t=0, i.e., ⟨σ_z(0)⟩=1, and a logarithmic time scale is used in Figs. <ref>, <ref> and <ref>.

We begin by considering room temperature effects on the D2F dynamics, since possible applications of this chromophore as an efficient electron acceptor are expected to work at such condition. As shown in Fig. <ref> (weak regime) and Fig. <ref> (strong regime), these are computed for Methanol and Toluene  in terms of the chromophore-solvent coupling η. As expected, the oscillatory behaviour of population inversion and coherence persists for longer times in the weak regime. However, it is worth noting that those oscillations are more favoured in Toluene (than in Methanol), regardless the strength of the coupling regime. 

In the strong regime, say η=0.25 (see violet line in Fig. <ref>) the population inversion and coherence oscillations disappear around 100 fs in Methanol, while they continue to occur up to ∼ 220 fs in Toluene. The process of electronic energy transfer is more coherent in Toluene as the interplay between the reorganisation energy and the inter-chromophore coupling lies in the coherent transfer regime λ≪ V_DA<cit.>. Indeed,  λ/V_DA∼0.7×10^-2 for Toluene while λ/V_DA∼0.3 for Methanol (values of λ are taken from Table <ref>).

In spite of the fact that the dynamics in Methanol becomes incoherent faster than in Toluene, the final population inversion is more efficient in the former solvent; i.e., the stronger the coupling η the higher the D-to-A transfer. In particular, the stationary population inversion in Toluene is ∼25%, and independent of η (panel (c) in Figs. <ref> and <ref>). For Methanol, although the standing population inversion starts at a ∼15% in the weak regime, it increases with η and reaches about a ∼40% in the strong regime (panel (a) of Figs. <ref> and <ref>). This implies that a higher ET can be obtained in Methanol at room temperature. This is in agreement with the ET rate computed below in Section <ref> in terms of the density operator. We stress that this phenomenon takes place at room temperature and that the ET does not increase at lower temperatures, as shown below in Fig. <ref>.

The studied D2F complex exhibits a long-term coherence that is comparable to the ones experimentally found  in some light-harvesting and organic materials <cit.>. As shown in Fig. <ref>(b) for Methanol, coherence persists above ∼100 fs in the strong coupling regime η_M=0.25,  where the system is found according to the computed reorganisation energy. Furthermore, the coherence lifetime is longer for Toluene (≳ 220 fs, see Fig. <ref>(d)). As expected, the decoherence process is slower for both solvents in the weak coupling regime (Figs. <ref>(b) and <ref>(d)). 
 
 In addition to the time scale introduced by the solvent relaxation (τ_s ≈ 100 fs; ħ_c=334 cm^-1 for Methanol, and τ_s≈ 220 fs; ħ_c=167 cm^-1 for Toluene), other relevant physical scales of interest are set by thermal fluctuations, the reorganisation energy,  and the intrinsic electronic coupling. In particular, thermal fluctuations set by k_BT/ħ introduce a time scale τ_T= _T^-1≡ħ/k_BT that affects the chromophore dynamics for t > τ_T. On the other hand, quantum vacuum fluctuations contribute to the dephasing process in the regime τ_s <t <τ_T. Solvent polarity also plays  a crucial role as it affects the electronic energy levels of the states involved in the transfer process. In our particular case, Methanol exhibits a separation between the considered electronic states higher than Toluene (ϵ_M>ϵ_T). A full numerical comparison of the different time scales is given in Table <ref> by means of the associated energy  ratios.
 
 


 


Temperature effects on the population inversion and coherence for a fixed D2F-solvent strong coupling are shown in Fig. <ref>. We consider the system-bath coupling that naturally identifies each solvent, i.e., η_M=0.25 for Methanol and η_T=0.33 for Toluene. In these Figures, the temperature is varied from 10 K to 300 K. ⟨σ_z(t)⟩ and (ρ_DA) exhibit oscillations showing the quantum character of the D2F dynamics even at room temperature along tens of femtoseconds. Despite the shortness of this quantum effect, experimental techniques involving femtosecond pump-probe spectroscopy and quantum coherence control are in place to detect and control the reported ultrafast quantum dynamics <cit.>.



 §.§ Electron Transfer in the D2F Complex
 

Electron transfer in the chromophore can be analysed throughout the chemical reaction rate, which is defined in terms of the flux-flux correlation function between the population of reactant and product states <cit.>. Thus, it is possible to relate the ET speed with the time-dependent density matrix elements previously computed in Section <ref>. In general, for an N-2-bridge system, the charge transfer from site m to site n can be described, within a charge hopping model <cit.>, by a jump rate k_mn≡ k_m → n, and the occupation probability (population) for each site can be written in terms of the ET rates by means of the following set of coupled differential equations:


    dρ_D/dt  = -k_D2ρ_D + k_2Dρ_2    
    
          ⋮
          dρ_N-1/dt = k_N-2N-1ρ_N-2  + k_AN-1ρ_A 
                       - (k_N-1N-2+k_N-1A)ρ_N-1
          dρ_A/dt = k_N-1Aρ_N-1 - k_AN-1ρ_A .


In our N=2-site model, this set reduces to: 


    dρ_D/dt  = -k_DAρ_D + k_ADρ_A,


for the donor, and the corresponding equation for the acceptor is given by the constrain ρ_D+ρ_A=1. The behaviour of the ET rate is quantified by considering the deviation from equilibrium of the density operator, which for the donor reads Δρ_D(t)=ρ_D(t)-ρ_D(t→∞). It is known that a time-dependent ET rate can be obtained for the expected exponential decay of the deviation approaching to zero. By taking the time derivative to Δρ_D(t), applying Eq. (<ref>), and assuming k_AD=0 for the sake of simplicity, according to the D2F chromophore site energies, the ET rate can be written as:


    k_DA=-1/ρ_D(t)dρ_D(t)/dt.


Although the value for the ET rate comes from taking the limit to infinity in Eq. (<ref>), this equation allows us to interpret the ET speed in terms of the quantum-mechanical dynamics of the chromophoric system. In so doing, we plot  the time dependence of k_DA as a function of the chromophore-solvent coupling strength in Fig. <ref>. An oscillatory behaviour is observed for k_DA before the ET rate reaches its stationary value. We point out that beyond the hopping model here described, the asymptotic k_DA relies on the physical properties of the system: the D-A coupling, the reorganization energy, etc. Indeed, we have compared our findings for the charge hopping model with the results found by means of the flickering resonance model, and corroborated that these converge to the same result for the ET rate (for a discussion on both methods see e.g. <cit.>). Thus, the matching with the computed values in the insets of Fig. <ref> is confirmed (numerics not shown).
As anticipated in Sec. <ref>, the ET rate is favoured by a strong D2F-solvent coupling, for both solvents. In particular, for Methanol, the ET rate increases up to one order of magnitude as the strength of the D2F-solvent coupling increases, as shown in the insets of Figs. <ref>(a) and (b). On the other hand, such a variation is less appreciable for Toluene, as can be seen in Figs. <ref>(c) and (d), which is in agreement with the quantum dynamics reported in Figs. <ref> and <ref>. 

The results obtained in Fig. <ref> are a strong indication that non-Markovianity assists the electronic transfer in the considered organic molecules, and we stress that such a transfer is optimised  at room temperature. On the other hand, our calculations also show that cooling down the system (see Figs. <ref>(b) and (d)), although enhances the chromophore quantum coherence, translates into a slight reduction of the coupling-assisted ET rate, as can directly be compared  between left (T=300 K) and right (T=50 K) panels, for both Methanol and Toluene. Thus, we conclude that ET rate enhancement in the chromophore is a direct consequence of the natural strong coupling between the complex and the solvent, and that this is not necessarily due to its quantum coherence properties (cfr. Figs. <ref> and <ref>).
 
Although a large number of molecules is required for direct applications in functional OPV devices, our results contribute to the understanding of fundamental ultrafast processes and reveal the importance of considering the quantum nature of the involved early-stage  electron transfer, whose dynamics and efficiency determine the function of the OPV material.   
Further investigations of solvent effects can take place by contemplating other relevant aspects, e.g.,  different spectral densities. In the case of many-body scenarios, techniques like machine learning <cit.> could also be applied to further  extend the scope of our findings. 



§ CONCLUSIONS


We have investigated the ultrafast quantum dissipative dynamics and electronic transfer processes in an electron acceptor push-pull Rhodanine derivative chromophore (D2F complex). In so doing, we have computed the effects due to the chromophore environment by considering polar (Methanol) and non polar (Toluene) solvents, under different temperature conditions. The D2F chemical structure (and the identification of suitable excited  and charge transfer states) allowed us to evaluate, within the Onsager model of solvation, the complex-solvent coupling strengths η, which clearly set the dynamics within the strong coupling regime.  The corresponding non-Markovian features of the solvent were treated  within a non-perturbative open quantum system approach, with the numerically exact  HEOM method <cit.>.

Our findings demonstrate the quantum character (coherent oscillations of the time-dependent density operator) of the D2F dynamics, even at room temperature, for both solvents. Particularly, that coherence is favoured in the case of Toluene thanks to the interplay between the reorganisation energy and the chromophore-chromophore coupling, in both complex-solvent coupling regimes. In contrast, Methanol exhibits a less coherent but higher ET rates, as it turns out that the stationary population inversion increases with coupling strength η. On the other hand, Toluene exhibits a steady population inversion that is invariant with respect to η.

To give a simple description of the ET process, we computed the ET rate k_DA as the deviation from equilibrium of the density operator probabilities. As mentioned before, this  ET rate is significantly  enhanced by Methanol. Our simulations  show that, at room temperature, the steady value of k_DA increases with the coupling strength η, and reaches values almost two orders of magnitude higher in Methanol than in Toluene. 
Furthermore, although prolonged coherence oscillations are observed at lower temperatures, the ET rate optimises at room temperature for both solvents as it reaches higher steady values (compare right and left panels in Fig. <ref>).

Finally, we conclude that the non polar solvent enhances the quantum properties of the organic chromophore dynamics when compared to the polar solvent; yet, and perhaps counterintuitive, it is the latter that exhibits the higher ET rates. Even though the question as to 
 whether quantum coherence is crucial/helps electron transfer mechanisms in light harvesting chromophoric systems remains open to debate, our findings reveal not to be the case in the D2F-solvent complex  here considered.







§ HEOM FOR THE REDUCED SYSTEM DYNAMICS



We briefly describe the Hierarchical Equations of Motion employed to obtain the dissipative dynamics of the reduced chromophoric system.
Considering the total system to be formed by the D2F molecule (reduced system, S) plus the solvent (bath, B), the time evolution of the associated density operator can be computed as <cit.>: iħρ̇_T(t)=ℒ_T(t)ρ_T(t)≡[H_T, ρ_T(t) ], where the total Hamiltonian H_T≡H (Sec. <ref>). Assuming an uncorrelated initial condition  ρ_S(0)⊗ρ_B, with the bath in the thermal state ρ_B=exp(-βH_B)/Tr(exp(-βH_B)), β=1/k_BT (k_B is the Boltzmann constant), the reduced system dynamics is obtained by tracing out over the bath degrees of freedom:


    ρ_S(t)= Tr_B(U(t)ρ_S(0)⊗ρ_B),


where <cit.> U(t)=exp(-i/ħ∫_0^tdτℒ_T(τ)). Hence, bath effects are encoded in the correlation functions <cit.>


    C_a(t)=⟨ u_a(t)u_a(0)⟩_B=1/π∫_0^∞ dω J_a(ω)e^-iω t/1-e^βħω,


where u_a(t)=∑_ξb_aξχ_ξ(t) are the time-dependent bath coupling operators. They arise from the system-bath Hamiltonian (<ref>) written in terms of the bath position coordinates χ_ξ. In general, this can be written as H_SB=∑_a,ξb_aξ|φ_a⟩⟨φ_a|χ_ξ. For the two-site model used in the main text, the coupling strength of Eq. (<ref>) simply reads g_ξ=b_ξ/√(2ħ m_ξω_ξ).

For a  Drude type 
spectral density (Eq.(<ref>)), the correlation functions can be expanded as:


    C_a^∞(t)=∑_k=0^∞c_ake^-ν_akt,


where ν_a0 = ω_c and ν_ak= 2π k/βħ (k≥ 1) are the Matsubara frequencies <cit.>, while the correlation coefficients <cit.> 
c_a0 = λ_aω_c/ħ[ (βħω_c/2) - i ] and 
c_ak=4λ_aω_c/βħ^2( ν_ak/ν_ak^2-ω_c^2), for k≥ 1. The sum in Eq.  (<ref>) must be truncated at some finite level K; ν_aKe^-ν_aK t≈δ(t).
Each exponential term in Eq. (<ref>) introduces a set of auxiliary density operators ρ_n that  are governed by the hierarchical structure of coupled equations of motion <cit.>:



    dρ_n/dt   =    -i/ħ[H_S,ρ_n]-∑_a=1^N∑_k=0^Kn_akν_akρ_n-i∑_a=1^N∑_k=0^K[|φ_n⟩⟨φ_n|,ρ_n_ak^+] 
          
    -∑_a=1^M(2λ_a/βħ^2ω_c -∑_k=0^Kc_ak/ν_ak) [ |φ_a⟩⟨φ_a|,[ |φ_a⟩⟨φ_a|,ρ_n] ] 
          
    - ∑_a=1^N∑_k=0^Kn_ak(c_a,k|φ_a⟩⟨φ_a|ρ_n_ak^- + ρ_n_ak^-|φ_a⟩⟨φ_a| c_a,k^*),



where each auxiliary density operator with nonnegative subscript n=(n_10, …, n_1K, …, n_N0, …, n_NK) is coupled to the auxiliary density operators with subscripts n_ak^±=(n_10, …, n_ak±1, …, n_NK). Each density operator in the hierarchy is assigned to a hierarchy level L=∑_a=1^N∑_k=0^K n_ak. The density operator ρ_S of the reduced system of interest corresponds to that one with subscript n=(0, …, 0).

The amount of auxiliary density operators (which rapidly increases with level L and in principle is infinite) accounts for the non-Markovian character of the system dynamics. Hence, a truncation such that L_max≫ω_max/min(ν_ak), with ω_max being a characteristic frequency of the system of interest, must be taken. Convergence of the HEOM is finally achieved after specifying the way to cutoff the terms beyond the truncation L_max. For doing so, the so-called time-local truncation can be used such that the Markovian approximation is employed for auxiliary density operators at level L_max. This implies that, for all auxiliary density operators with L=L_max-1,


    ∑_k=0^Kρ_n_ak^+≈    -i( Q_a^K(t)ρ_n-ρ_nQ_a^K(t)^†),


where Q_a^K(t)= ∫_0^tC_a^K(τ)e^( -i/ħH_Sτ)|φ_a⟩⟨φ_a|e^(i/ħH_Sτ)dτ.




§ ACKNOWLEDGEMENTS


We acknowledge J. Arce and A. Ortiz for helpful discussions. This work was supported by the Colombian Science, Technology and Innovation Fund-General Royalties System (Fondo CTeI-Sistema General de Regalías) under contract BPIN 2013000100007. C.E.S. thanks Universidad de Córdoba for partial support (Grant No. CA-097).

ieeetr


