Simulating Quantum Vibronic Dynamics at Finite Temperatures With Many Body Wave Functions at 0 K

For complex molecules, nuclear degrees of freedom can act as an environment for the electronic “system” variables, allowing the theory and concepts of open quantum systems to be applied. However, when molecular system-environment interactions are non-perturbative and non-Markovian, numerical simulations of the complete system-environment wave function become necessary. These many body dynamics can be very expensive to simulate, and extracting finite-temperature results—which require running and averaging over many such simulations—becomes especially challenging. Here, we present numerical simulations that exploit a recent theoretical result that allows dissipative environmental effects at finite temperature to be extracted efficiently from a single, zero-temperature wave function simulation. Using numerically exact time-dependent variational matrix product states, we verify that this approach can be applied to vibronic tunneling systems and provide insight into the practical problems lurking behind the elegance of the theory, such as the rapidly growing numerical demands that can appear for high temperatures over the length of computations.

1. Introduction

The dissipative quantum dynamics of electronic processes play a crucial role in the physics and chemistry of materials and biological life, particularly in the ultra-fast and non-equilibrium conditions typical of photophysics, nanoscale charge transfer and glassy, low-temperature phenomena (Miller et al., 1983). Indeed, the through-space tunneling of electrons, protons and their coupled dynamics critically determine how either ambient energy is transduced, or stored energy is utilized in supramolecular “devices,” and real-time dynamics are especially important when the desired processes occur against thermodynamical driving forces, or at the single-to-few particle level (Devault, 1980; May and Kühn, 2008).

In many physio-chemical systems, a reaction, energy transfer, or similar event proceeds in the direction of a free energy gradient, necessitating the dissipation of energy and the generation of entropy (Dubi and Dia Ventra, 2011; Benenti et al., 2017). A powerful way of modeling the microscopic physics at work during these irreversible dynamics is the concept of an “open” quantum system (Breuer and Petruccione, 2002; Weiss, 2012). Here a few essential and quantized degrees of freedom constituting the “system” are identified and explicitly coupled to a much larger number of “environmental” degrees of freedom. Equations of motion for the coupled system and environment variables are then derived and solved, with the goal of obtaining the behavior of the “system” degrees of freedom once the unmeasureable environmental variables are averaged over their uncertain initial and final states. It is in this “tracing out” of the environment that the originally conservative, reversible dynamics of the global system gives way to apparently irreversible dynamics in the behavior of the system's observable variables. The effective behavior of the system “opened” to the environment is entirely contained within its so-called reduced density matrix, which we shall later define. Important examples of the emergent phenomenology of reduced density matrices include the ubiquitous processes of thermalization, dephasing, and decoherence.

In the solid state, a typical electronic excitation will interact weakly with the lattice vibrations of the material, particularly the long-wavelength, low frequency modes. Under such conditions it is often possible to treat the environment with low-order perturbation theory and—given that the lattice “environment” relaxes back to equilibrium very rapidly—it is possible to derive a Markovian master equation for the reduced density matrix, such as the commonly used Bloch-Redfield theory (Breuer and Petruccione, 2002; May and Kühn, 2008; Weiss, 2012). However, in sufficiently complex molecular systems, such as organic bio-molecules, the primary environmental degrees of freedom acting on electronic states are typically the stochastic vibrational motions of the atomic nuclear coordinates. Unlike the solid state, these vibrations can: (1) couple non-perturbatively to electronic states, (2) relax back to equilibrium on timescales that are longer than the dynamics they induce in the system, and (3) have frequencies ω such that ℏω ≫ K_BT, where T is the environmental temperature, and so must be treated quantum mechanically (zero-point energy and nuclear quantum effects). In this regime, the theory and numerical simulation of open quantum systems becomes especially challenging, as the detailed dynamics of the interacting system and environmental quantum states need to be obtained, essentially requiring the solution of a correlated (entangled) many body problem.

One well-known and powerful approach to this problem in theoretical chemistry is the Multi-layer Multiconfigurational Time-dependent Hartree (ML-MCTDH) technique, which enables vibronic wave functions to be efficiently represented and propagated without the a priori limitations due to the “curse of dimensionality” associated with many body quantum systems (Lubich, 2015; Wang and Shao, 2019). However, computationally demanding methods based on the propagation of a large wave function from a definite initial state will typically struggle when dealing with finite-temperature environments (vide infra), as the probability distribution of initial states requires extensive sampling. For this reason, the majority of ML-MCTDH studies have been effectively on zero-temperature systems.

In this article we will explore a recent and intriguing development in an alternative approach to real-time dynamics and chemical rate prediction. This approach is based on the highly efficient representation and manipulation of large, weakly entangled wave functions with DMRG, Matrix-Product, and Tensor-Network-State methods (Orus, 2014). These methods, widely used in condensed matter, quantum information and cold atom physics, have recently been applied to a range of open system models, including chemical systems, but—as wave function methods—are typically used at zero-temperature (Prior et al., 2010, 2013; Chin et al., 2013; Alvertis et al., 2019; Schröder et al., 2019; Xie et al., 2019). However, a remarkable new result due to Tamascelli et al. shows that it is indeed possible to obtain the finite-temperature reduced dynamics of a system based on a simulation of a “pure,” i.e., zero-temperature wave function (Tamascelli et al., 2019).

In principle, this opens the way for many existing wave function methods to be extended into finite temperature regimes, although the present formulation of Tamascelli et al.'s T-TEDOPA mapping is most easily implemented with matrix product states (MPS). In this article, we shall investigate this extension to finite temperature in the regime of relevance for molecular quantum dynamics, that is, non-perturbative vibrational environments, and present numerical data that verifies the elegance and utility of the method, as well as some of the potential issues arising in implementation.

The structure of the article is as follows. In section 2, we will summarize Tamascelli et al.'s T-TEDOPA mapping. In section 3, we verify the theory by comparing numerical simulations against an exactly solvable open system model, and also employ further numerical investigations to provide some insight into the manner in which finite temperatures are handled within this method. By looking at the observables of the environment, we find that the number of excitations in the simulations grows continuously over time, which may place high demands on computational resources in some problems. In section 4, we will present results for a model system inspired by electron transfer in a multi-dimensional vibrational environment, and show how the temperature-driven transition from quantum tunneling to classical barrier transfer are successfully captured by this new approach. This opens a potentially fruitful new phase for the application of tensor network and related many body approaches for the simulation of non-equilibrium dynamics in a wide variety of vibronic materials and molecular reactions.

2. T-TEDOPA

In this section we shall summarize the essential features of the T-TEDOPA approach, closely following the original notation and presentation of Tamascelli et al. (2019). Our starting point is the generic Hamiltonian for a system coupled to a bosonic environment consisting of a continuum of harmonic oscillators

$\begin{array}{ll}{H}_{SE}=H_S+H_E+H_I,& (1)\end{array}$

where

$\begin{array}{ll}{H}_I=A_S\otimes {\displaystyle {\int }_0^{\infty }}d\omega {Ô}_{\omega },H_E={\displaystyle {\int }_0^{\infty }}d\omega \omega a_{\omega }^{†}{a}_{\omega }.& (2)\end{array}$

The Hamiltonian H_S is the free system Hamiltonian, which for chemical systems, molecular photophysics and related problems will often be a description of a few of the most relevant diabatic states at some reference geometry of the environment(s) (May and Kühn, 2008). A_S is the system operator which couples to the bath. For the bath operators we take the displacements

$\begin{array}{ll}{O}_{\omega }=\sqrt{J(\omega )}(a_{\omega }+a_{\omega }^{†}),& (3)\end{array}$

thus defining the spectral density J(ω). This has been written here as a continuous function, but coupling to a discrete set of vibrational modes in, say, a molecular chromophore, can be included within this description by adding suitable structure to the spectral density, i.e., sets of Lorentzian peaks or Dirac functions (Wilhelm et al., 2004; Schulze and Kuhn, 2015; Mendive-Tapia et al., 2018). The state of the system+environment at time t is described by a mixed state described by a density matrix ρ_SE(t). The initial condition is assumed to be a product of system and environment states ρ_SE(0) = ρ_S(0) ⊗ ρ_E (0) where ρ_S(0) is an arbitrary density matrix for the system and ${\rho }_E(0)=exp(-H_E\beta )/{Z}$, with the environment partition function given by ${Z}=\text{Tr}\left\{exp(-H_E\beta )\right\}$. Such a product state is commonly realized in photophysics, where the reference geometry for the environment is the electronic ground state and the electronic system is excited according to the Franck-Condon principle into some manifold of electronic excited states without nuclear motion (Mukamel, 1995; May and Kühn, 2008). Indeed, this can also occur following any sufficiently rapid non-adiabatic event, just as ultra-fast charge separation at a donor-acceptor interface (Gélinas et al., 2019; Smith and Chin, 2015). The environment thus begins in a thermal equilibrium state with inverse temperature β, and the energy levels of each harmonic mode are statistically populated, as shown in Figure 1A. For a very large (continuum) of modes, the number of possible thermal configurations of the initial probability distribution grows extremely rapidly with temperature, essentially making a naive sampling of these configurations impossible for full wave function simulations. We note, however, that some significantly better sampling methods involving sparse grids and/or stochastic mean-field approaches have been proposed and demonstrated (Alvermann and Fehske, 2009; Binder and Burghardt, 2019).

FIGURE 1

Figure 1. (A) A generic open quantum system contains a few-level “system” (S) that interacts with a much larger thermal heat bath of bosonic oscillators (the environment, E). The continuum of oscillator modes are initially uncorrelated with the system and each is thermally occupied with characteristic temperature T = β⁻¹. Coupling and stochastic fluctuations of the environment lead to the effective thermalization of the system, once the environmental states have been traced over. (B) In the T-TEDOPA approach, the harmonic environment is extended to include modes of negative frequency, and all modes (positive and negative frequency) are initially in their ground states. It can be formally demonstrated that the thermalization of S in (A) can always be obtained from the pure zero-temperature state in (B), provided the spectral density of the original environment is known.

The initial thermal condition of the environmental oscillators is also a Gaussian state, for which is it further known that the influence functional (Weiss, 2012)—which is a full description of the influence of the bath on the system—will depend only on the two-time correlation function of the bath operators

$\begin{array}{ll}S(t)={\displaystyle {\int }_0^{\infty }}d\omega \langle O_{\omega }(t)O_{\omega }(0)\rangle .& (4)\end{array}$

Any two environments with the same S(t) will have the same influence functional and thus give rise to the same reduced system dynamics, i.e., the same ρ_S(t) = Tr{ρ_SE(t)}. That the reduced systems dynamics are completed specified by the spectral density and temperature of a Gaussian environment has been known for a long time (Weiss, 2012), but the key idea of the equivalence—and thus the possibility of the interchange—of environments with the same correlation functions has only recently been demonstrated by Tamascelli et al. (2018).

The time dependence in Equation (4) refers to the interaction picture so that the bath operators evolve under the free bath Hamiltonian: $O_{\omega }(t)=e^{i{H}_{E}t}{O}_{\omega }(0)e^{-i{H}_{E}t}$. Using Equation (3) and $\langle a_{\omega }^{†}{a}_{\omega }\rangle =n_{\beta }(\omega )$ we have

$\begin{array}{ll}S(t)={\displaystyle {\int }_0^{\infty }}J(\omega )\left[e^{-i\omega t}(1+n_{\beta }(\omega ))+e^{i\omega t}{n}_{\beta }(\omega )\right].& (5)\end{array}$

Making use of the relation

$\begin{array}{ll}\frac{1}{2}(1+\coth (\omega \beta /2))\equiv \{\begin{array}{c}{n}_{\omega }(\beta ),\omega \ge 0\hfill \\ -(n_{\left|\omega \right|}(\beta )+1),\omega <0\hfill \end{array}& (6)\end{array}$

we can write Equation (5) as an integral over all positive and negative ω

$\begin{array}{ll}S(t)={\displaystyle {\int }_{-\infty }^{\infty }}d\omega \text{sign}(\omega )\frac{J(|\omega |)}{2}(1+coth(\frac{\omega \beta }{2}))e^{-i\omega t}.& (7)\end{array}$

But Equation (7) is exactly the two-time correlation function one would get if the system was coupled to a bath, now containing positive and negative frequencies, at zero temperature, with a temperature weighted spectral density given by

$\begin{array}{ll}{J}_{\beta }(\omega )=\text{sign}(\omega )\frac{J(|\omega |)}{2}(1+coth(\frac{\omega \beta }{2})).& (8)\end{array}$

Thus, we find that our open system problem is completely equivalent to the one governed by the Hamiltonian

$\begin{array}{ll}H=H_S+H_E^{\text{ext}}+H_I^{\text{ext}},& (9)\end{array}$

in which the system couples to an extended environment, where

$\begin{array}{ll}\begin{array}{c}{H}_I^{\text{ext}}=A_S\otimes {\displaystyle {\int }_{-\infty }^{\infty }}d\omega \sqrt{J_{\beta }(\omega )}(a_{\omega }+a_{\omega }^{†}),\\ H_E^{\text{ext}}={\displaystyle {\int }_{-\infty }^{\infty }}d\omega \omega a_{\omega }^{†}{a}_{\omega },\hfill \end{array}& (10)\end{array}$

and which has the initial condition ρ_SE(0) = ρ_S(0) ⊗ |0〉_E 〈0|. The system now couples to a bath consisting of harmonic oscillators of positive and negative frequencies which are initially in their ground states, as shown in Figure 1B. This transformed initial condition is now far more amenable to simulation as the environment is now described by a pure, single-configuration wave function, rather than a statistical mixed state, and so no statistical sampling is required to capture the effects of temperature on the reduced dynamics!

Analyzing the effective spectral density of Equation (8), it can be seen that the new extended environment has thermal detailed balance between absorption and emission processes encoded in the ratio of the coupling strengths to the positive and negative modes in the extended Hamiltonian, as opposed to the operator statistics of a thermally occupied state of the original, physical mode, i.e.,

$\begin{array}{ll}\frac{J_{\beta }(\omega )}{J_{\beta }(-\omega )}=\frac{{\langle a_{\omega }{a}_{\omega }^{†}\rangle }_{\beta }}{{\langle a_{\omega }^{†}{a}_{\omega }\rangle }_{\beta }}=e^{\beta \omega }& (11)\end{array}$

Indeed, from the system's point of view, there is no difference between the absorption from an occupied, positive energy, bath mode and the emission into an unoccupied, negative energy, bath mode.

In fact, the equivalence between these two environments goes beyond the reduced system dynamics as there exists a unitary transformation which links the extended environment to the original thermal environment. This means that one is able to reverse the transformation and calculate thermal expectations for the actual bosonic bath such as ${\langle a_{\omega }^{†}(t)a_{\omega }(t)\rangle }_{\beta }$. This is particularly useful for molecular systems in which environmental (vibrational) dynamics are also important observables that report on the mechanisms and pathways of physio-chemical transformations (Musser et al., 2015; Schnedermann et al., 2016, 2019). This is a major advantage of many body wave function approaches, as the full information about the environment is available, c.f. effective master equation descriptions which are obtained after averaging over the environmental state. We note that the idea of introducing a second environment of negative frequency oscillators to provide finite temperature effects in pure wave functions was previously proposed in the thermofield approach of de Vega and Bañuls (2015). This approach explicitly uses the properties of two-mode squeezed states to generate thermal reduced dynamics, but the original thermofield approach, unlike the T-TEDOPA mapping, considered the positive and negative frequency environments as two separate baths.

Following this transformation a further step is required to facilitate efficient simulation of the many-body system+environment wave function. This is to apply a unitary transformation to the bath modes which converts the star-like geometry of $H_I^{\mathrm{\text{ext}}}$ into a chain-like geometry, thus allowing the use of MPS methods (Chin et al., 2010, 2013; Prior et al., 2013). We thus define new modes $c_n^{(†)}={\int }_{-\infty }^{\infty }{U}_n(\omega )a_{\omega }^{(†)}$, known as chain modes, via the unitary transformation $U_n(\omega )=\sqrt{J_{\beta }(\omega )}{p}_n(\omega )$ where p_n(ω) are orthonormal polynomials with respect to the measure dωJ_β(ω). Thanks to the three term recurrence relations associated with all orthonormal polynomials p_n(ω), only one of these new modes, n = 1, will be coupled to the system, while all other chain modes will be coupled only to their nearest neighbors (Chin et al., 2010). Our interaction and bath Hamiltonians thus become

$\begin{array}{ll}\begin{array}{c}{H}_I^{\text{chain}}=\kappa A_S(c_1+c_1^{†}),\hfill \\ H_E^{\text{chain}}={\displaystyle \sum _{n=1}^{\infty }}{\omega }_n{c}_n^{†}{c}_n+{\displaystyle \sum _{n=1}^{\infty }}(t_n{c}_n^{†}{c}_{n+1}+h.c).\end{array}& (12)\end{array}$

The chain coefficients appearing in Equation (12) are related to the three-term recurrence parameters of the orthonormal polynomials and can be computed using standard numerical techniques (Chin et al., 2010). The full derivation of the above Hamiltonian is given in the Appendix. Since the initial state of the bath was the vacuum state, it is unaffected by the chain transformation.

We have thus arrived at a formulation of the problem of finite-temperature open systems in which the many-body environmental state is initialized as a pure product of trivial ground states, whilst the effects of thermal fluctuations and populations are encoded in the Hamiltonian chain parameters and system-chain coupling. These parameters must be determined once for each temperature but—in principle—the actual simulation of the many body dynamics is now no more complex than a zero-temperature simulations. This thus opens up the use of powerful T = 0K wave function methods for open systems, such as those based on MPS, numerical renormalization group and ML-MCTDH (Lubich, 2015; Wang and Shao, 2019). However, while this seems remarkable—and we believe this mapping to be a major advance—there must be a price to be paid elsewhere. We shall now demonstrate with numerical examples where some of the computational costs for including finite-T effects may appear and discuss how they might effect the feasibility and precision of simulations. We also propose a number of ways to mitigate these potential problems within the framework of tensor network approaches.

3. Numerical Tests and Computational Efficiency

All numerical results in the following sections are obtained by representing the many body system-environment wave function as a MPS and evolving it using time-dependent variational methods. All results have been converged w.r.t. the parameters of MPS wave functions (bond dimensions, local Hilbert space dimensions, integrator time steps), meaning that the results and discussion should—unless explicitly stated—pertain to the essential properties of the T-TEDOPA mapping itself. Extensive computational details and background theory can be founds in Orus (2014), Schollwock (2011), Lubich et al. (2015), Paeckel et al. (2019), and Haegeman et al. (2016).

3.1. Chain Dynamics and Chain-Length Truncation

Before looking at the influence of thermal bath effects on a quantum system, we first investigate the effects of the changing chain parameters that appear due to the inclusion of temperature in the effective spectral density J_β(ω). As a consequence of the nearest-neighbor nature of Equation (12) (see Figure 2), the chain mapping establishes a kind of causality among the bath modes which is extremely convenient for simulation. Starting from t = 0 the system will interact first with the chain mode n = 1 which, as well as acting back on the system, will in turn excite the next mode along the chain and so on. The dynamics thus have a well-defined light-cone structure in which a perturbation travels outwards from the system along the chain to infinity. This means that we may truncate the chain at any distant mode n = N without causing an error in the system or bath observables up to a certain time T_LC(N) which is the time it takes for the edge of the light-cone to reach the Nth chain mode. Beyond T_LC(N) there will be reflections off the end of the chain leading to error in the bath observables, however these reflections will not cause error in the system observables until the time t ≈ 2T_LC(N). Figure 3 shows a snapshot of the chain mode occupations for the Ohmic spin-boson model considered in the next section. One can see that the velocity of the wave-front that travels outward from the system depends on temperature, with hotter baths leading to faster propagation and thus requiring somewhat longer chains.

FIGURE 2

Figure 2. (A) The extended proxy environment of Figure 1A is described by an effective, temperature-dependent spectral density J_β(ω). Once the effective J_β(ω) has been specified, new oscillator modes can be found that provide a unitary transformation to a linear chain representation of the environment with nearest neighbor interactions. The non-perturbative wave function dynamics for such a many-body 1D system can be very efficiently simulated with MPS methods. (B) J_β(ω) for a physical Ohmic environment at three representative temperatures. At very low temperature (ω_c β ≫1) there is essentially no coupling to the negative frequency modes, as excitation of these modes leads to an effective absorption of heat from the environment. At higher temperatures, J_β(ω) becomes increasingly symmetric for the positive and negative modes.

FIGURE 3

Figure 3. Chain mode occupations $\langle c_n^{†}{c}_n\rangle$ at time ω_ct = 45 for baths of several temperatures. The system, which in this case is the Ohmic SBM, with ω₀ = 0.2ω_c and α = 0.1, is attached at site n = 1 of the chain.

To enable simulation we are also required to truncate the infinite Fock space dimension of each chain mode to a finite dimension d, introducing an error for which there exist rigorously derived bounds (Woods et al., 2015). The initial state |Ψ(0)〉_SE = |ψ(0)〉_S ⊗ |0〉_E (here we specialize to the case where the system is initially in a pure state) can then be encoded in an MPS and evolved under one of the many time-evolution methods for MPS. We choose to use the one-site Time-Dependent-Variational-Principle (1TDVP) as it has been shown to be a efficient method for tracking long-time thermalization dynamics and has previously been shown to give numerically exact results for the zero-temperature spin-boson model in the highly challenging regime of quantum criticality (Schröder and Chin, 2016). In our implementation of 1TDVP the edge of the light-cone is automatically estimated throughout the simulation by calculating the overlap of the wave-function |Ψ(t)〉_SE with its initial value |Ψ(0)〉_SE at each chain site. This allows us expand the MPS dynamically to track to expanding light-cone, providing roughly a 2-fold speed-up compared to using a fixed length MPS.

3.2. Two-Level System Dynamics: Dephasing and Divergence of Chain Occupations Due to Energy Exchange

To confirm the accuracy of this approach in terms of reduced system dynamics we now explore the effects of a dissipative environment on a quantum two-level system. First, we compare the numerical results against the analytically solvable Independent-Boson-Model (IBM) (Mahan, 2000; Breuer and Petruccione, 2002). This is a model of pure dephasing, defined by $H_S=\frac{{\omega }_0}{2}{\sigma }_z$ and A_S = σ_z, where {σ_x, σ_y, σ_z} are the standard Pauli matrices. We take an Ohmic spectral density with a hard cut-off J(ω) = 2αωΘ(ω−ω_c) and choose a coupling strength of α = 0.1 and a gap of ω₀ = 0.2ω_c for the two level system (TLS). The initial state of the system is a positive superposition of the spin-up and spin-down states, and we monitor the decay of the TLS coherence, which is quantified by 〈σ_x(t)〉. All results were converged using a Fock space dimension of d = 6 for the chain modes and maximum MPS bond-dimension D_max = 4. We find that the results obtained using the T-TEDOPA method agree very well with the exact solution (see Figure 4) and correctly reproduce the transition from under-damped to over-damped decay as the temperature is increased (Mahan, 2000; Breuer and Petruccione, 2002).

FIGURE 4

Figure 4. Comparison of T-TEDOPA (Black crosses) with the exact solution for the Independent-Boson-Model at β = 100 (Red), β = 10 (Blue), and β = 1 (Green). $H_S=\frac{{\omega }_0}{2}{\sigma }_z$, A_S = σ_z, $J(\omega )=2\alpha {\omega }_c{(\frac{\omega }{{\omega }_c})}^s\Theta (\omega -{\omega }_c)),$, α = 0.1, s = 1, ω₀ = 0.2ω_c.

As a second numerical example we take the Spin-Boson-Model (SBM), identical to the IBM considered above except that now the TLS couples to the bath via A_S = σ_x. Unlike the previous case, the bath can now drive transitions within the TLS, so that energy is now dynamically exchanged between the TLS and its environment. Indeed, as A_S no longer commutes with H_S, no exact solution for this model is known (Weiss, 2012). It has thus become an important testing ground for numerical approaches to non-perturbative simulations of open systems and has been widely applied to the physics of decoherence, energy relaxation and thermalization in diverse physical, chemical and biological systems—see Weiss, 2012; De Vega and Alonso, 2017 for extensive references. In our example, we prepare the spin in the upper spin state (〈σ_z〉 = +1) and allow the bath to thermalize by environmental energy exchange (see Figure 1A). Here, instead of presenting the spin dynamics for this model we will here interest ourselves in the observables of the bath as these will provided insight into the manner in which a finite temperature bath is being mimicked by an initially empty tight-binding chain. In Figure 5, we plot the bath mode occupations $\langle a_{\omega }^{†}{a}_{\omega }\rangle$ for several temperatures. Each observation was taken after the spin had decayed into its thermal steady state and thus provides a kind of absorption spectrum for the system. We note that these data refer to the modes of the extended environment of Equation (9) rather than the original bosonic bath and thus the mode energies run from −ω_c to ω_c.

FIGURE 5

Figure 5. Bath mode occupations $\langle n_{\omega }\rangle =\langle a_{\omega }^{†}{a}_{\omega }\rangle$ for the extended environment after the TLS has decayed. The TLS is governed by a Hamiltonian $H_S=\frac{{\omega }_0}{2}{\sigma }_z$ where ω₀ = 0.2ω_c and is coupled to an Ohmic bath with a hard cut-off via A_S = σ_x. The coupling strength is α = 0.1. Left inset: total mode occupation as a function of time ${\langle n\rangle }_{\mathrm{\text{tot}}}={\int }_{-\infty }^{\infty }d\omega \langle n_{\omega }\rangle$. Right inset shows $\frac{\langle n_{{\omega }_p}\rangle +1}{\langle n_{{\omega }_n}\rangle }$ plotted on a log scale against the inverse temperature, demonstrating the detailed balance of the absorption and emission rates.

We find that for zero temperature (β = ∞) the bath absorption spectrum contains a single peak at a frequency around ω_p = 0.17ω_c, suggesting that the spin emits into the bath at a re-normalized frequency that is lower than the bare gap of the TLS (ω₀ = 0.2ω_c). This agrees well with the renormalized gap ${\omega }_0^r={\omega }_0{({\omega }_0/{\omega }_c)}^{\frac{\alpha }{1-\alpha }}$ predicted by the non-perturbative variational polaron theory of Silbey and Harris (1984), which for the parameters used here gives ${\omega }_0^r=0.167{\omega }_c$.

Moving to non-zero temperature we see that a peak begins to form at a corresponding negative frequency, which we interpret as being due the spin absorbing thermal energy from the bath by the emission (creation) of negative energy quanta. In accordance with detailed balance, the ratio between the positive and negative frequency peaks approaches unity as temperature is increased and by βω_c = 2 the two peaks have merged to form a single, almost symmetric, distribution, reflecting the dominance of thermal absorption and emission over spontaneous emission at high temperature. Indeed, as shown in the right inset of Figure 5 the ratio of the peak heights we extract obeys $\frac{\langle n_{\omega }\rangle +1}{\langle n_{-\omega }\rangle }=e^{ϵ\beta }$ with ϵ = 0.118. Thus we see that the chain is composed of two independent vacuum reservoirs of positive and negative energy which the system emits into at rates which effectively reproduce the emission and absorption dynamics that would be induced by a thermal bath.

However, the introduction of positive and negative modes has an interesting and important consequence for the computational resources required for simulation. Shown in the left inset of Figure 5 is the total mode occupation as a function of time for some of the different temperatures simulated. One sees that for β = ∞ (zero temperature) the total occupation of the bath modes increases initially and then plateaus at a certain steady state value corresponding to the total number of excitations created in the bath by the TLS during its decay. In contrast, for finite temperature, the total mode occupation increases indefinitely at a rate which grows with temperature. This is despite the fact that for the finite temperature baths the total excitation number will also reach a steady state once the TLS has decayed. The reason for this is clear. The thermal occupation of the physical bath mode with frequency ω is obtained by subtracting its negative, from its positive energy counterpart in the extended mode basis, i.e., 〈_nω〉_β = 〈_nω〉_|0〉E − 〈_n−ω〉_|0〉E. While 〈_nω〉β will reach a steady state, the components 〈_nω〉_|0〉E and 〈_n−ω〉_|0〉E will be forever increasing, reflecting the fact that the TLS reaches a dynamic equilibrium with the bath in which energy is continuously being absorbed from and emitted into the bath at equal rates, thus filling up the positive and negative reservoirs. Since, it is the modes of the extended environment that appear in the numerical simulation, one will always encounter potentially large errors once the filling of the modes exceeds their capacity set by the truncation to d Fock states per oscillator. The rate at which this filling occurs increases with temperature and is linear in time. However, as the relaxation time of the system is also broadly proportional to temperature for βω_c ≪ 1, this may not be a problem, if one is only interested in the short-time transient dynamics. Where this may pose problems is for the extraction of converged properties of relaxed, i.e., locally thermalized excited states, such as their (resonance) fluorescence spectra, or multidimensional optical spectra (Mukamel, 1995). While these ever-growing computational resources must—as argued above—be present in any simulation approach, we note that one possible way to combat the growth of local dimensions could be to use the dynamical version of Guo's Optimized Boson Basis (OBB) which was introduced into 1TDVP for open systems by Schroeder et al. (Guo et al., 2012; Schröder and Chin, 2016).

4. Electron Transfer

Having established that the T-TEDOPA mapping allows efficient computational access to finite temperature open dynamics, we now study the chemically relevant problem of tunneling electron transfer. Electron transfer is a fundamental problem in chemical dynamics and plays an essential role in a vast variety of crucial processes including the ultra-fast primary electron transfer step in photosynthetic reaction centers and the electron transport that powers biological respiration (Devault, 1980; Marcus, 1993; May and Kühn, 2008). The problem of modeling electron transfer between molecules comes down to accurately treating the coupling between the electronic states and environmental vibrational modes, and often involves the use of first principle techniques to parameterize the total spectral functions of the vibrational and outer solvent, or protein environment (Mendive-Tapia et al., 2018; Schröder et al., 2019; Zuehlsdorff et al., 2019). In many molecular systems—and particularly biological systems where the transfer between electronic states is affected by coupling to chromophore and protein modes—the system-bath physics is highly non-perturbative and J(ω) has very sharp frequency-dependence (May and Kühn, 2008; Womick et al., 2011; Kolli et al., 2012; Chin et al., 2013). Until recently, and even at zero temperature, a fully quantum mechanical description of the coupling to a continuum of environmental vibrations was challenging due to the exponential scaling of the vibronic wave functions. However, with advances in numerical approaches driven by developments in Tensor-Networks and ML-MCTDH, the exact quantum simulation of continuum environment models can now be explored very precisely at zero temperature. Given this, we now explore how the T-TEDOPA mapping can extend this capability to finite temperature quantum tunneling.

Here, we will again adapt the spin-boson model to analyse a typical donor-acceptor electron transfer system, as shown in Figure 6. In this model the electron transfer process is modeled using two states representing the reactant and product states which we take to be the eigenstates of σ_x with |↓〉 representing the reactant and |↑〉 the product. We take our system Hamiltonian to be $H_S=\frac{ϵ}{2}{\sigma }_z+{\lambda }_R\frac{1+{\sigma }_x}{2}$, and the coupling operator as $A_S=\frac{1+{\sigma }_x}{2}$, where λ_R is the reorganization energy which for an Ohmic bath is λ_R = 2αω_c. The electron tunnels from the environmentally relaxed reactant state to the product state by moving through a multi-dimensional potential energy landscape along a collective reaction coordinate which is composed of the displacements of the ensemble of bath modes (this is effectively the coordinate associated with the mode that is directly coupled to the system in the chain representation of the environment). Figure 6A shows two potential energy surfaces—Marcus parabolas—of the electronic system for ϵ = 0. Although in the actual model we simulate the reaction coordinate is composed of the displacements of an infinite number of modes, in Figure 6 we present a simplified picture in which the electron moves along a single reaction coordinate, x. The potential minimum of the reactant state corresponds to the bath in its undisplaced, vacuum state, whereas at the potential minimum of the product state each bath mode is displaced by an amount depending on its frequency and the strength of its coupling to the TLS $\sqrt{J(\omega )}/\omega$. The presence of the reorganization energy in H_S ensures that these two minima are degenerate in energy and thus detailed balance will ensure an equal forward and backward rate.

FIGURE 6

Figure 6. (A) Potential energy surfaces (Marcus parabolas) for ϵ = 0 as a function of the reaction coordinate x. We consider only the case of zero bias, i.e., when the minima of the two wells are at the same energy. (B) Turning the electronic coupling ϵ leads to an avoided crossing and thus an energy barrier E_b for the reaction. Note that this is a simplified picture in which we treat the bath as being represented by a single mode of frequency ω and coupling strength g whereas in the actual model we simulate there is a similar surface for all bath modes.

Turning on the coupling ϵ between the two levels leads to an avoided crossing in the two energy surfaces in an adiabatic representation of the vibronic tunneling system, leading to two potential wells. In such a semi-classical (Born-Oppeheimer) picture, we see that the electron must overcome a kind of effective energy barrier E_b that scales with the total reorganization energy of the entire environment λ_R in order for the reaction to progress. We thus might well expect to see thermally activated (exponential) behavior whereby the tunneling rate ∝exp(−βE_b). However, at low temperatures this behavior should be dramatically quenched and dissipative quantum tunneling should become dominant and strongly dependent on the spectral function of the environment (Weiss, 2012).

4.1. Numerical Results

For our numerical investigation we take an Ohmic spectral density with α = 0.8 for which the dynamics are expected to be incoherent at all temperatures, i.e., the energy surfaces of Figure 6B are well-separated and friction is such that there will be no oscillatory tunneling dynamics between reactant and product. In Figure 7, we present results for this model at several temperatures using the T-TEDOPA mapping and 1TDVP. The expectation of σ_x can be taken to be a measure of the progress of the reaction, starting at the value of −1 when the system is entirely in the reactant state, and approaching 0 as the electron tunnels through the barrier and the populations thermalize. We find that as the temperature is increased the dynamics tend to an exponential decay to the steady state, whereas non-exponential behavior is observed for lower temperatures. In Figure 7B, we show the expectation of σ_y, which is the conjugate coordinate to the σ_x and which may thus be interpreted as a kind of momentum associated with the tunneling. We find that there is a sharp initial spike in 〈σ_y〉 which decays with oscillations which are increasingly damped at higher temperatures. As we might have predicted, these transient dynamics occur on a timescale of $\tau \approx {\omega }_c^{-1}$, which the fastest response time of an environment with an upper cut-off frequency of ω_c. This is approximately the timescale over which the environment will adjust to the sudden presence of the electron, and essentially sets the timescale for the formation of the adiabatic landscape (or, alternatively, for the formation of the dressed polaron states), after which the tunneling dynamics proceed. This period is related to the slippage of initial conditions that is sometimes used to fix issues of density matrix positivity in perturbative Redfield Theory (Gaspard and Nagaoka, 1999), although here the establishment of these conditions is described exactly and in real-time. We also see that the crossover to the tunneling regime happens faster as the temperature increases, meaning that the effective initial conditions—particularly 〈σ_y(t)〉—are temperature dependent.

FIGURE 7

Figure 7. (A) 〈σ_x(t)〉 for several temperatures, which represents the progress of the reaction. The decay to the steady state is exponential at high temperature. (B) 〈σ_y(t)〉, representing the momentum along the reaction coordinate. We encounter some noise beyond about ω_ct = 50 in the β = 2 data. This is as a result of the truncation of the local Hilbert spaces of the bath modes (cf. section 3). The inset shows an enlarged view of the initial fast dynamics which appear to be broadly independent of temperature.

We extract approximate reaction rates from the TLS dynamics by fitting each 〈σ_x(t)〉 to an exponential decay −e^−Γt on timescales t > τ. We thus obtain the rates Γ(ϵ, β) for the various values of β and ϵ simulated. The values of ϵ were chosen to be small compared to the characteristic vibrational frequency of the bath, ϵ ≪ ω_c and to the reorganization energy, ϵ ≪ λ_R and thus lie in the non-adiabatic regime which is the relevant regime for electron transfer. One may then perform a perturbative expansion in ϵ, otherwise known as the “Golden Rule” approach which, for an Ohmic bath, yields the following formulas for the high and low temperature limits corresponding respectively to the classical and quantum regimes (Weiss, 2012).

$\begin{array}{ll}\Gamma (\beta )=\{\begin{array}{c}\frac{\sqrt{\pi }}{4\sqrt{\alpha }}{ϵ}^2{(\frac{\pi }{\beta {\omega }_c})}^{2\alpha -1},\beta {\omega }_c\gg 1\hfill \\ \frac{e^2}{4}\sqrt{\frac{\pi \beta {\omega }_c}{2\alpha }}\exp (-\frac{\alpha \beta {\omega }_c}{2}),\beta {\omega }_c\ll 1\end{array}.& (13)\end{array}$

The golden rule result is based on second-order perturbation in the tunneling coupling ϵ, but it is exact to all orders in the system-environment coupling α. Additionally, the Ohmic form of the spectral function generates a non-trivial power-law dependence of the tunneling rate on the temperature for βω_c ≫ 1 in which the rate may either decrease or increase as the temperature is lowered, depending on the value of α. We plot these formulas along with the numerically evaluated rates in Figure 8. There is a good agreement in the high and low temperature limits between the Golden Rule expressions and the T-TEDOPA results, and one clearly sees that the temperature dependence of the rate is non-monotonic with a transition from power law growth (quantum, 2α − 1 > 0) to power-law decay (classical, $\propto \sqrt{\beta }$) as the temperature increases from T = 0. We note that for the parameters we present here, the intermediate regime where thermally activated behavior is predicted βω_c ~ 1 is not observed for the Ohmic environment, and one essentially switches from tunneling limited by the effect of friction on the attempt frequency to the low-temperature polaronic tunneling of Equation (13).

FIGURE 8

Figure 8. Log plot of the rates, Γ, extracted from 〈σ_x(t)〉 for ϵ = 0.2 (Red), ϵ = 0.3 (Blue), and ϵ = 0.4 (Red) as a function of β. (Dashed lines) High temperature (T≫ω_c), classical, limit of Golden Rule formula. (Dotted lines) Low temperature (T ≪ ω_c), quantum, limit of Golden Rule formula.

5. Conclusion

In this article we have shown how the combination of the Tamasceli's remarkable T-TEDOPA mapping and non-perturbative variational Tensor-Network dynamics can be applied to chemical and photophysical systems under laboratory conditions. Through numerical experiments we have carefully investigated how the T-TEDOPA mapping allows the effects of finite temperatures to be obtained efficiently without any need for costly sampling of the thermal environment state, or the explicit use of density matrices. However, analysis of these environmental dynamics reveals how incorporating finite temperatures can lead to more expensive simulations, due to the filling-up of the chain modes and the longer chains that are needed to prevent recurrence dynamics. Yet, we believe that this method, and others like it, based on the exact quantum many-body treatment of vibrational modes (Somoza et al., 2019), could present an attractive complementary approach to the Multi-Layer Multi-Configurational Time-Dependent Hartree Method (MLMCTDH) commonly used in chemical dynamics. One possible direction for this would be to consider a problem in which a (discretized) potential surface for a reaction is contained within the system Hamiltonian, while the environment bath provides the nuclear thermal and quantum fluctuations that ultimately determine both real-time kinetics and thermodynamical yields for the process, as is currently captured in methods such as Ring Polymer Molecular Dynamics (Craig and Manolopoulos, 2004). Furthermore, the Tensor-Network structures are not limited to the simple chain geometries we consider here but can in fact adopt a tree structure, thus enabling the treatment of complex coupling to multiple independent baths (Schröder et al., 2019). Such trees tensor networks have recently been interfaced with ab initio methods to explore ultra-fast photophysics of real molecules and their pump-probe spectra (Schnedermann et al., 2019), but such efforts have so far been limited to zero temperature. Finally, the cooperative, antagonistic or sequential actions of different types of environments, i.e., light and vibrations (Wertnik et al., 2018), or even the creation of new excitations, such as polaritons (Memmi et al., 2017; Del Pino et al., 2018; Herrera and Owrutsky, 2020), could play a key role in sophisticated new materials for energy transduction, catalysis or regulation (feedback) of reactions, and T-TDEPODA-based tensor networks are currently being used to explore these developing areas.

Data Availability Statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: https://github.com/angusdunnett/MPSDynamics.

Author Contributions

AD implemented the T-TEDOPA mapping in a bespoke 1TDVP code and performed the numerical simulations. AD and AC wrote the manuscript. AC oversaw the project. All authors contributed to the article and approved the submitted version.

Funding

AD was supported by the Ecole Doctorale 564 Physique en Ile-de-France. AC was partly supported by ANR project No. 195608/ACCEPT.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

References

Alvermann, A., and Fehske, H. (2009). Sparse polynomial space approach to dissipative quantum systems: Application to the sub-ohmic spin-boson model. Phys. Rev. Lett. 102:150601. doi: 10.1103/PhysRevLett.102.150601

PubMed Abstract | CrossRef Full Text | Google Scholar

Alvertis, A. M., Schröder, F. A. Y. N., and Chin, A. W. (2019). Non-equilibrium relaxation of hot states in organic semiconductors: Impact of mode-selective excitation on charge transfer. J. Chem. Phys. 151:084104. doi: 10.1063/1.5115239

PubMed Abstract | CrossRef Full Text | Google Scholar

Benenti, G., Casati, G., Saito, K., and Whitney, R. (2017). Fundamental aspects of steady-state conversion of heat to work at the nanoscale. Phys. Rep. 694, 1–124. doi: 10.1016/j.physrep.2017.05.008

CrossRef Full Text | Google Scholar

Binder, R., and Burghardt, I. (2019). First-principles quantum simulations of exciton diffusion on a minimal oligothiophene chain at finite temperature. Faraday Discuss. 221, 406–427. doi: 10.1039/C9FD00066F

PubMed Abstract | CrossRef Full Text | Google Scholar

Breuer, H.-P., and Petruccione, F. (2002). The Theory of Open Quantum Systems. Oxford University Press.

Google Scholar

Chin, A., Prior, J., Rosenbach, R., Caycedo-Soler, F., Huelga, S. F., and Plenio, M. B. (2013). The role of non-equilibrium vibrational structures in electronic coherence and recoherence in pigment-protein complexes. Nat. Phys. 9, 113–118. doi: 10.1038/nphys2515

CrossRef Full Text | Google Scholar

Chin, A. W., Rivas, Á., Huelga, S. F., and Plenio, M. B. (2010). Exact mapping between system-reservoir quantum models and semi-infinite discrete chains using orthogonal polynomials. J. Math. Phys. 51:092109. doi: 10.1063/1.3490188

CrossRef Full Text | Google Scholar

Craig, I. R., and Manolopoulos, D. E. (2004). Quantum statistics and classical mechanics: real time correlation functions from ring polymer molecular dynamics. J. Chem. Phys. 121, 3368–3373. doi: 10.1063/1.1777575

PubMed Abstract | CrossRef Full Text | Google Scholar

De Vega, I., and Alonso, D. (2017). Dynamics of non-Markovian open quantum systems. Rev. Modern Phys. 89:015001. doi: 10.1103/RevModPhys.89.015001

CrossRef Full Text | Google Scholar

de Vega, I., and Bañuls, M.-C. (2015). Thermofield-based chain-mapping approach for open quantum systems. Phys. Rev. A 92:052116. doi: 10.1103/PhysRevA.92.052116

CrossRef Full Text | Google Scholar

Del Pino, J., Schröder, F. A., Chin, A. W., Feist, J., and Garcia-Vidal, F. J. (2018). Tensor network simulation of polaron-polaritons in organic microcavities. Phys. Rev. B 98:165416. doi: 10.1103/PhysRevB.98.165416

CrossRef Full Text | Google Scholar

Devault, D. (1980). Quantum mechanical tunnelling in biological systems. Q. Rev. Biophys. 13, 387–564. doi: 10.1017/S003358350000175X

PubMed Abstract | CrossRef Full Text | Google Scholar

Dubi, Y., and Dia Ventra, M. (2011). Colloquium: heat flow and thermoelectricity in atomic and molecular junctions. Rev. Modern Phys. 83, 131–155. doi: 10.1103/RevModPhys.83.131

CrossRef Full Text | Google Scholar

Gaspard, P., and Nagaoka, M. (1999). Slippage of initial conditions for the redfield master equation. J. Chem. Phys. 111, 5668–5675. doi: 10.1063/1.479867

CrossRef Full Text | Google Scholar

Gélinas, S., Rao, A., Kumar, A., Smith, S. L., Chin, A. W., Clark, J., et al. (2014). Ultrafast long-range charge separation in organic semiconductor photovoltaic diodes. Science 343, 512–516. doi: 10.1126/science.1246249

PubMed Abstract | CrossRef Full Text | Google Scholar

Guo, C., Weichselbaum, A., von Delft, J., and Vojta, M. (2012). Critical and strong-coupling phases in one-and two-bath spin-boson models. Phys. Rev. Lett. 108:160401. doi: 10.1103/PhysRevLett.108.160401

PubMed Abstract | CrossRef Full Text | Google Scholar

Haegeman, J., Lubich, C., Oseledets, I., Vandereycken, B., and Verstraete, F. (2016). Unifying time evolution and optimization with matrix product states. Phys. Rev. B 94:165116. doi: 10.1103/PhysRevB.94.165116

CrossRef Full Text | Google Scholar

Herrera, F., and Owrutsky, J. (2020). Molecular polaritons for controlling chemistry with quantum optics. J. Chem. Phys. 152:100902. doi: 10.1063/1.5136320

PubMed Abstract | CrossRef Full Text | Google Scholar

Kolli, A., O'Reilly, E. J., Scholes, G. D., and Olaya-Castro, A. (2012). The fundamental role of quantized vibrations in coherent light harvesting by cryptophyte algae. J. Chem. Phys. 137:174109. doi: 10.1063/1.4764100

PubMed Abstract | CrossRef Full Text | Google Scholar

Lubich, C. (2015). Time integration in the multiconfiguration time-dependent Hartree method of molecular quantum dynamics. Appl. Math. Res. Express 2015, 311–328. doi: 10.1093/amrx/abv006

CrossRef Full Text | Google Scholar

Lubich, C., Oseledets, I., and Vandereycken, B. (2015). Time integration of tensor trains. SIAM J. Num. Anal. 53, 917–941. doi: 10.1137/140976546

CrossRef Full Text | Google Scholar

Mahan, G. D. (2000). Many-Particle Physics. Boston, MA: Springer US. doi: 10.1007/978-1-4757-5714-9

CrossRef Full Text | Google Scholar

Marcus, R. A. (1993). Electron transfer reactions in chemistry. Theory and experiment. Rev. Modern Phys. 65:599. doi: 10.1103/RevModPhys.65.599

CrossRef Full Text | Google Scholar

May, V., and Kühn, O. (2008). Charge and Energy Transfer Dynamics in Molecular Systems. Weinheim: John Wiley & Sons.

Google Scholar

Memmi, H., Benson, O., Sadofev, S., and Kalusniak, S. (2017). Strong coupling between surface plasmon polaritons and molecular vibrations. Phys. Rev. Lett. 118:126802. doi: 10.1103/PhysRevLett.118.126802

PubMed Abstract | CrossRef Full Text | Google Scholar

Mendive-Tapia, D., Mangaud, E., Firmino, T., de la Lande, A., Desouter-Lecomte, M., Meyer, H.-D., et al. (2018). Multidimensional quantum mechanical modeling of electron transfer and electronic coherence in plant cryptochromes: the role of initial bath conditions. J. Phys. Chem. B 122, 126–136. doi: 10.1021/acs.jpcb.7b10412

PubMed Abstract | CrossRef Full Text | Google Scholar

Miller, W. H., Schwartz, S. D., and Tromp, J. W. (1983). Quantum mechanical rate constants for bimolecular reactions. J. Chem. Phys. 79, 4889–4898. doi: 10.1063/1.445581

CrossRef Full Text | Google Scholar

Mukamel, S. (1995). Principles of Nonlinear Optical Spectroscopy, Vol. 6. New York, NY: Oxford University Press.

Google Scholar

Musser, A. J., Liebel, M., Schnedermann, C., Wende, T., Kehoe, T. B., Rao, A., et al. (2015). Evidence for conical intersection dynamics mediating ultrafast singlet exciton fission. Nat. Phys. 11, 352–357. doi: 10.1038/nphys3241

CrossRef Full Text | Google Scholar

Orus, R. (2014). A practical introduction to tensor networks: matrix product states and projected entangled pair states. Ann. Phys. 349, 117–158. doi: 10.1016/j.aop.2014.06.013

CrossRef Full Text | Google Scholar

Paeckel, S., Kohler, T., Swoboda, A., Manmana, S. R., Schollwock, U., and Hubig, C. (2019). Time-evolution methods for matrix-product states. Ann. Phys. 411:167998. doi: 10.1016/j.aop.2019.167998

CrossRef Full Text | Google Scholar

Prior, J., Chin, A. W., Huelga, S. F., and Plenio, M. B. (2010). Efficient simulation of strong system-environment interactions. Phys. Rev. Lett. 105:050404. doi: 10.1103/PhysRevLett.105.050404

PubMed Abstract | CrossRef Full Text | Google Scholar

Prior, J., de Vega, I., Chin, A. W., Huelga, S. F., and Plenio, M. B. (2013). Quantum dynamics in photonic crystals. Phys. Rev. A 87:013428. doi: 10.1103/PhysRevA.87.013428

CrossRef Full Text | Google Scholar

Schnedermann, C., Alvertis, A. M., Wende, T., Lukman, S., Feng, J., Schröder, F. A., et al. (2019). A molecular movie of ultrafast singlet fission. Nat. Commun. 10, 1–11. doi: 10.1038/s41467-019-12220-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Schnedermann, C., Lim, J. M., Wende, T., Duarte, A. S., Ni, L., Gu, Q., et al. (2016). Sub-10 fs time-resolved vibronic optical microscopy. J. Phys. Chem. Lett. 7, 4854–4859. doi: 10.1021/acs.jpclett.6b02387

PubMed Abstract | CrossRef Full Text | Google Scholar

Schollwock, U. (2011). The density-matrix renormalization group in the age of matrix product states. Ann. Phys. 326, 96–192. doi: 10.1016/j.aop.2010.09.012

CrossRef Full Text | Google Scholar

Schröder, F. A., and Chin, A. W. (2016). Simulating open quantum dynamics with time-dependent variational matrix product states: towards microscopic correlation of environment dynamics and reduced system evolution. Phys. Rev. B 93:075105. doi: 10.1103/PhysRevB.93.075105

CrossRef Full Text | Google Scholar

Schröder, F. A. Y. N., Turban, D. H. P., Musser, A. J., Hine, N. D. M., and Chin, A. W. (2019). Tensor network simulation of multi-environmental open quantum dynamics via machine learning and entanglement renormalisation. Nat. Commun. 10:1062. doi: 10.1038/s41467-019-09039-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Schulze, J., and Kuhn, O. (2015). Explicit correlated exciton-vibrational dynamics of the FMO complex. J. Phys. Chem. B 119, 6211–6216. doi: 10.1021/acs.jpcb.5b03928

PubMed Abstract | CrossRef Full Text | Google Scholar

Silbey, R., and Harris, R. A. (1984). Variational calculation of the dynamics of a two level system interacting with a bath. J. Chem. Phys. 80, 2615–2617. doi: 10.1063/1.447055

CrossRef Full Text | Google Scholar

Smith, S. L., and Chin, A. W. (2015). Phonon-assisted ultrafast charge separation in the PCBM band structure. Phys. Rev. B 91:201302. doi: 10.1103/PhysRevB.91.201302

CrossRef Full Text | Google Scholar

Somoza, A. D., Marty, O., Lim, J., Huelga, S. F., and Plenio, M. B. (2019). Dissipation-assisted matrix product factorization. Phys. Rev. Lett. 123:100502. doi: 10.1103/PhysRevLett.123.100502

PubMed Abstract | CrossRef Full Text | Google Scholar

Tamascelli, D., Smirne, A., Huelga, S. F., and Plenio, M. B. (2018). Nonperturbative treatment of non-markovian dynamics of open quantum systems. Phys. Rev. Lett. 120:030402. doi: 10.1103/PhysRevLett.120.030402

PubMed Abstract | CrossRef Full Text | Google Scholar

Tamascelli, D., Smirne, A., Lim, J., Huelga, S. F., and Plenio, M. B. (2019). Efficient simulation of finite-temperature open quantum systems. Phys. Rev. Lett. 123:090402. doi: 10.1103/PhysRevLett.123.090402

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, H., and Shao, J. (2019). Quantum phase transition in the spin-boson model: a multilayer multiconfiguration time-dependent Hartree study. J. Phys. Chem. A 123, 1882–1893. doi: 10.1021/acs.jpca.8b11136

CrossRef Full Text | Google Scholar

Weiss, U. (2012). Quantum Dissipative Systems, 4th Edn. Singapore: World Scientific. doi: 10.1142/8334

CrossRef Full Text | Google Scholar

Wertnik, M., Chin, A., Nori, F., and Lambert, N. (2018). Optimizing co-operative multi-environment dynamics in a dark-state-enhanced photosynthetic heat engine. J. Chem. Phys. 149:084112. doi: 10.1063/1.5040898

PubMed Abstract | CrossRef Full Text | Google Scholar

Wilhelm, F., Kleff, S., and Von Delft, J. (2004). The spin-boson model with a structured environment: a comparison of approaches. Chem. Phys. 296, 345–353. doi: 10.1016/j.chemphys.2003.10.010

CrossRef Full Text | Google Scholar

Womick, J. M., Liu, H., and Moran, A. M. (2011). Exciton delocalization and energy transport mechanisms in r-phycoerythrin. J. Phys. Chem. A 115, 2471–2482. doi: 10.1021/jp111720a

PubMed Abstract | CrossRef Full Text | Google Scholar

Woods, M., Cramer, M., and Plenio, M. (2015). Simulating bosonic baths with error bars. Phys. Rev. Lett. 115:130401. doi: 10.1103/PhysRevLett.115.130401

PubMed Abstract | CrossRef Full Text | Google Scholar

Xie, X., Liu, Y., Yao, Y., Schollwock, U., Liu, C., and Ma, H. (2019). Time-dependent density matrix renormalization group quantum dynamics for realistic chemical systems. J. Chem. Phys. 151:224101. doi: 10.1063/1.5125945

PubMed Abstract | CrossRef Full Text | Google Scholar

Zuehlsdorff, T. J., Montoya-Castillo, A., Napoli, J. A., Markland, T. E., and Isborn, C. M. (2019). Optical spectra in the condensed phase: capturing anharmonic and vibronic features using dynamic and static approaches. J. Chem. Phys. 151:074111. doi: 10.1063/1.5114818

PubMed Abstract | CrossRef Full Text | Google Scholar

Appendix

The chain mapping used in section 2 is based on the theory of orthogonal polynomials. A polynomial of degree n is defined by

$\begin{array}{ll}{p}_n(x)={\displaystyle \sum _{m=0}^n}{a}_m{x}^m.& (\mathrm{A1})\end{array}$

The space of polynomials of degree n is denoted ℙ_n and is a subset of the space of all polynomials ℙ_n ⊂ ℙ. Given a measure dμ(x) which has finite moments of all orders on some interval [a, b], we may define the inner product of two polynomials

$\begin{array}{ll}{\langle p,q\rangle }_{\mu }={\displaystyle {\int }_a^b}d\mu (x)p(x)q(x).& (\mathrm{A2})\end{array}$

This inner product gives rise to a unique set of orthonormal polynomials $\left\{{\tilde{p}}_n\in {\mathbb{P}}_n,n=0,1,2,\dots \right\}$ which all satisfy

$\begin{array}{ll}\langle {\tilde{p}}_n,{\tilde{p}}_m\rangle ={\delta }_{n,m}.& (\mathrm{A3})\end{array}$

This set forms a complete basis for ℙ, and more specifically the set $\left\{{\tilde{p}}_n\in {\mathbb{P}}_n,n=0,1,2,\dots m\right\}$ is a complete basis for ${\bigcup }_{r=1}^m{\mathbb{P}}_r$.

It is often useful to express the orthonormal polynomials in terms of the orthogonal monic polynomials π_n(x) which are the unnormalized scalar multiples of ${\tilde{p}}_n(x)$ whose leading coefficient is 1 (a_n = 1)

$\begin{array}{ll}{\tilde{p}}_n(x)=\frac{{\pi }_n(x)}{\Vert {\pi }_n\Vert }.& (\mathrm{A4})\end{array}$

The key property of orthogonal polynomials for the construction of the chain mapping is that they satisfy a three term recurrence relation

$\begin{array}{ll}{\pi }_{k+1}(x)=(x-{\alpha }_k){\pi }_k(x)-{\beta }_k{\pi }_{k-1}(x),& (\mathrm{A5})\end{array}$

where it can be easily shown that

$\begin{array}{ll}{\alpha }_k=\frac{\langle x{\pi }_k,{\pi }_k\rangle }{\langle {\pi }_k,{\pi }_k\rangle },{\beta }_k=\frac{\langle {\pi }_k,{\pi }_k\rangle }{\langle {\pi }_{k-1},{\pi }_{k-1}\rangle }.& (\mathrm{A6})\end{array}$

Now that we have defined the orthogonal polynomials we may use them to construct the unitary transformation that will convert the star Hamiltonian of Equation (9) with

$\begin{array}{ll}{H}_I^{\text{ext}}=A_S\otimes {\displaystyle {\int }_{-\infty }^{\infty }}d\omega \sqrt{J_{\beta }(\omega )}(a_{\omega }+a_{\omega }^{†}),H_E^{\text{ext}}={\displaystyle {\int }_{-\infty }^{\infty }}d\omega \omega a_{\omega }^{†}{a}_{\omega },& (\mathrm{A7})\end{array}$

into the chain Hamiltonian of Equation (12). The transformation is given by

$\begin{array}{ll}{c}_n^{(†)}={\displaystyle {\int }_{-\infty }^{\infty }}{U}_n(\omega )a_{\omega }^{(†)},& (\mathrm{A8})\end{array}$

where

$\begin{array}{ll}{U}_n(\omega )=\sqrt{J_{\beta }(\omega )}{\tilde{p}}_n(\omega )=\sqrt{J_{\beta }(\omega )}\frac{{\pi }_n(\omega )}{\Vert {\pi }_n\Vert },& (\mathrm{A9})\end{array}$

and the polynomials ${\tilde{p}}_n(\omega )$ are orthonormal with respect to the measure dωJ_β(ω). The unitarity of U_n(ω) follows immediately from the orthonormality of the polynomials.

Applying the above transformation to the interaction Hamiltonian we have

$\begin{array}{ll}{H}_I^{\text{ext}}=A_S\otimes {\displaystyle \sum _{n=0}^{\infty }}{\displaystyle {\int }_{-\infty }^{\infty }}d\omega J_{\beta }(\omega )\frac{{\pi }_n(\omega )}{\Vert {\pi }_n\Vert }(c_n^{†}+c_n)& (\mathrm{A10})\end{array}$

For the zeroth order monic polynomial we have π₀ = 1 and so we may insert this into the above expression

$\begin{array}{ll}{H}_I^{\text{ext}}=A_S\otimes {\displaystyle \sum _{n=0}^{\infty }}{\displaystyle {\int }_{-\infty }^{\infty }}d\omega J_{\beta }(\omega )\frac{{\pi }_n(\omega ){\pi }_0}{\Vert {\pi }_n\Vert }(c_n^{†}+c_n).& (\mathrm{A11})\end{array}$

Recognizing the inner product in the above expression and making use of the orthogonality of the polynomials we have

$\begin{array}{ll}{H}_I^{\text{ext}}=A_S\otimes {\displaystyle \sum _{n=0}^{\infty }}\Vert {\pi }_n\Vert {\delta }_{n,0}(c_n^{†}+c_n)=A_S\otimes \Vert {\pi }_0\Vert (c_0^{†}+c_0),& (\mathrm{A12})\end{array}$

and thus, in the new basis, only one mode now couples to the system.

Now for the environment part of the Hamiltonian we have

$\begin{array}{ll}{H}_E^{\text{ext}}={\displaystyle \sum _{n,m=0}^{\infty }}{\displaystyle {\int }_{-\infty }^{\infty }}d\omega J_{\beta }(\omega )\omega \frac{{\pi }_n(\omega ){\pi }_m(\omega )}{\Vert {\pi }_n\Vert \Vert {\pi }_m\Vert }{c}_n^{†}{c}_m.& (\mathrm{A13})\end{array}$

Substituting for ωπ_n(ω) from the three term recurrence relation of Equation (A5) yields

$\begin{array}{l}{H}_E^{\text{ext}}={\displaystyle \sum _{n,m=0}^{\infty }}{\displaystyle {\int }_{-\infty }^{\infty }}d\omega \frac{J_{\beta }(\omega )}{\Vert {\pi }_n\Vert \Vert {\pi }_m\Vert }[{\pi }_{n+1}(\omega )+{\alpha }_n{\pi }_n(\omega )\\ +{\beta }_n{\pi }_{n-1}(\omega )]{\pi }_m(\omega )c_n^{†}{c}_m.& (\mathrm{A14})\end{array}$

Again, evaluating the inner products we have

$\begin{array}{l}{H}_E^{\text{ext}}={\displaystyle \sum _{n,m=0}^{\infty }}\frac{1}{\Vert {\pi }_n\Vert }[\Vert {\pi }_m\Vert {\delta }_{n+1,m}+{\alpha }_n\Vert {\pi }_m\Vert {\delta }_{n,m}\\ +{\beta }_n\Vert {\pi }_m\Vert {\delta }_{n-1,m}]c_n^{†}{c}_m\\ ={\displaystyle \sum _{n=0}^{\infty }}\sqrt{{\beta }_{n+1}}{c}_n^{†}{c}_{n+1}+{\alpha }_n{c}_n^{†}{c}_n+\sqrt{{\beta }_{n+1}}{c}_n^{†}{c}_{n-1},& (\mathrm{A15})\end{array}$

where in the second line we have used the fact that

$\begin{array}{ll}\frac{\Vert {\pi }_{n+1}\Vert }{\Vert {\pi }_n\Vert }=\sqrt{{\beta }_{n+1}}.& (\mathrm{A16})\end{array}$

We thus arrive at the nearest-neighbor coupling Hamiltonian of Equation (12) and are able to identify the chain coefficients as

$\begin{array}{l}\kappa =\Vert {\pi }_0\Vert ,\\ {\omega }_{n+1}={\alpha }_n,\\ t_n=\sqrt{{\beta }_n}.& (\mathrm{A17})\end{array}$

Note that in Equation (12) the chain sites are labeled starting from n = 1 and not n = 0 as in Equation (A15). All that remains now to calculate the chain coefficients for a particular spectral density J_β(ω) is to compute the recurrence coefficients, α_n and β_n, and this may done iteratively using Equations (A5) and (A6) and numerically evaluating the inner product integrals using a quadrature rule.