Simulating quantum vibronic dynamics at finite temperatures with many body wave functions at 0K

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.


I. 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, lowtemperature phenomena [1].Indeed, the through-space tunneling of electrons, protons and their coupled dynamics critically determine how either ambient energy is transduced, or stored energy is utilised 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 [2,3].
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 [4,5].A powerful way of modelling the microscopic physics at work during these irreversible dynamics is the concept of an 'open' quantum system [6,7].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 behaviour 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 behaviour of the system's observable variables.The effective behaviour 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 [3,6,7].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 B T , 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 [8,9].However, computationally demanding methods based on the propagation of a large wave function from a definite ini-tial state will typically struggle when dealing with finitetemperature 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 realtime 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 [10].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 zerotemperature [11][12][13][14][15][16].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. zerotemperature wave function [17].
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 II we will summarise Tamascelli et al.'s T-TEDOPA mapping.In section III 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 IV we will present results for a model system inspired by electron transfer in a multidimensional 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 nonequilibrium dynamics in a wide variety of vibronic materials and molecular reactions.g k < l a t e x i t s h a 1 _ b a s e 6 4 = " E 0 G 8 K F T R G + G s F P n c q K k H o q + l J q g = " > A A A B 6 n i c b V B N S 8 N A E J 3 U r 1 q / q h 6 9 L B b B U 0 l E 0 G P R i 8 e K 9 g P a U D b b T b p 0 s w m 7 E 6 G E / g Q v H h T x 6 i / y 5 r 9 x 2 + a g r Q 8 G H u / N M D M v S K U w 6 L r f T m l t f W N z q 7 x d 2 d n d 2 z + o H h 6 1 T Z J p x l s s k Y n u B t R w K R R v o U D J u 6 n m N A 4 k 7 w T j 2 5 n f e e L a i E Q 9 4 i T l f k w j J U L B K F r p I R q M B 9 W a W 3 f n I K v E K 0 g N C j Q H 1 a / + M G F Z z B U y S Y 3 p e W 6 K f k 4 1 C i b 5 t N L P D E 8 p G 9 O I 9 y x V N O b G z + e n T s m Z V Y Y k T L Q t h W S u / p 7 I a W z M J A 5 s Z 0 x x Z J a 9 m f i f 1 8 s w v P Z z o d I M u W K L R W E m C S Z k 9 j c Z C s 0 Z y o k l l G l h b y V s R D V l a N O p 2 B C 8 5 Z d X S f u i 7 r l 1 7 / 6 y 1 r g p 4 i j D C Z z C O X h w B Q 2 4 g y a 0 g E E E z / A K b 4 5 0 X p x 3 5 2 P R W n K K m W P 4 A + f z B 0 m y j c k = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " E 0 G 8 K F T R G + G s F P n c q K k H o q + l J q g = " > A A A B 6 n i c b V B N S 8 N A E J 3 U r 1 q / q h 6 9 L B b B U 0 l E 0 G P R i 8 e K 9 g P a U D b b T b p 0 s w m 7 E 6 G E / g Q v H h T x 6 i / y 5 r 9 x 2 + a g r Q 8 G H u / N M D M v S K U w 6 L r f T m l t f W N z q 7 x d 2 d n d 2 z + o H h 6 1 T Z J p x l s s k Y n u B t R w K R R v o U D J u 6 n m N A 4 k 7 w T j 2 5 n f e e L a i E Q 9 4 i T l f k w j J U L B K F r p I R q M B 9 W a W 3 f n I K v E K 0 g N C j Q H 1 a / + M G F Z z B U y S Y 3 p e W 6 K f k 4 1 C i b 5 t N L P D E 8 p G 9 O I 9 y x V N O b G z + e n T s m Z V Y Y k T L Q t h W S u / p 7 I a W z M J A 5 s Z 0 x x Z J a 9 m f i f 1 8 s w v P Z z o d I M u W K L R W E m C S Z k 9 j c Z C s 0 Z y o k l l G l h b y V s R D V l a N O p 2 B C 8 5 Z d X S f u i 7 r l 1 7 / 6 y 1 r g p 4 i j D C Z z C O X h w B Q 2 4 g y a 0 g E E E z / A K b 4 5 0 X p x 3 5 2 P R W n K K m W P 4 A + f z B 0 m y j c k = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " E 0 G 8 K F T R G + G s F P n c q K k H o q + l J q g = " > A A A B 6 n i c b V B N S 8 N A E J 3 U r 1 q / q h 6 9 L B b B U 0 l E 0 G P R i 8 e K 9 g P a U D b b T b p 0 s w m 7 E 6 G E / g Q v H h T x 6 i / y 5 r 9 x 2 + a g r Q 8 G H u / N M D M v S K U w 6 L r f T m l t f W N z q 7 x d 2 d n d 2 z + o H h 6 1 T Z J p x l s s k Y n u B t R w K R R v o U D J u 6 n m N A 4 k 7 w T j 2 5 n f e e L a i E Q 9 m 6 e 9 F + / d + 5 i 3 r n j 5 z B H 8 g f f 5 A 3 k q j w w = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " d e H L w Y 3 G 4 O 3 P P d I F 9 3 w s u H 1 9 y D 4 = " > A m 6 e 9 F + / d + 5 i 3 r n j 5 z B H 8 g f f 5 A 3 k q j w w = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " d e H L w Y 3 G 4 O 3 P P d I F 9 3 w s u H 1 9 y D 4 = " > A m 6 e 9 F + / d + 5 i 3 r n j 5 z B H 8 g f f 5 A 3 k q j w w = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " d e H L w Y 3 G 4 O 3 P P d I F 9 3 w s u H 1 9 y D 4 = " > A m 6 e 9 F + / d + 5 i 3 r n j 5 z B H 8 g f f 5 A 3 k q j w w = < / l a t e x i t > +! < l a t e x i t s h a 1 _ b a s e 6 4 = " F X u p 3 s e i t e D l M 8 f w B 9 7 n D / d N j 1 A = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " F X u p 3 s e i t e D l M 8 f w B 9 7 n D / d N j 1 A = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " F X u p 3 s e i t e D l M 8 f w B 9 7 n D / d N j 1 A = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " F X u p

II. T-TEDOPA
In this section we shall summarise the essential features of the T-TEDOPA approach, closely following the original notation and presentation of Tamascelli et al [17].Our starting point is the generic Hamiltonian for a system coupled to a bosonic environment consisting of a continuum of harmonic oscillators where 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) [3].A S is the system operator which couples to the bath.For the bath operators we take the displacements 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 [18][19][20].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 ρ       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 neighbour interactions.The non-perturbative wave function dynamics for such a manybody 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.
where ρ S (0) is an arbitrary density matrix for the system and ρ E (0) = exp(−H E β)/Z, with the environment partition function given by Z = Tr{exp(−H E β)}.Such a product state is commonly realised 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 [3,21].Indeed, this can also occur following any sufficiently rapid non-adiabatic event, just as ultra-fast charge separation at a donor-acceptor interface [22,23].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 Fig. 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 [24,25].
The initial thermal condition of the environmental oscillators is also a Gaussian state, for which is it further known that the influence functional [7] -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 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 [7], 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. [26].
The time dependence in eq. 4 refers to the interaction picture so that the bath operators evolve under the free bath Hamiltonian: O ω (t) = e iH E t O ω (0)e −iH E t .Using eq. 3 and a † ω a ω = n β (ω) we have Making use of the relation we can write eq. 5 as an integral over all positive and negative ω ))e −iωt .( But eq. 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 Thus, we find that our open system problem is completely equivalent to the one governed by the Hamiltonian in which the system couples to an extended environment, where 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 Fig. 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!Analysing the effective spectral density of Eq. 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.
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 a † ω (t)a ω (t) β .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 [27][28][29].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 Vegas and Banulus [30].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 manybody system+environment wave function.This is to apply a unitary transformation to the bath modes which converts the star -like geometry of H ext I into a chainlike geometry, thus allowing the use of MPS methods [12,13,31].We thus define new modes c ω , known as chain modes, via the unitary transformation U n (ω) = J β (ω)p n (ω) 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 neighbours [31].Our interaction and bath Hamiltonians thus become The chain coefficients appearing in eq. 12 are related to the three-term recurrence parameters of the orthonormal polynomials and can be computed using standard numerical techniques [31].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 manybody environmental state is initialised 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 renormalisation group and ML-MCTDH [8,9].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.

III. NUMERICAL TESTS AND COMPUTATIONAL EFFICIENCY
All numerical results in the following sections are obtained by representing the many body systemenvironment 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 Refs.[10,[32][33][34][35].
A. 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-neighbour nature of eq. 12 (see Fig. 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 N th 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.
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 [36].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 [37].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.

B. 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) [6,38].This is a model of pure dephasing, defined by H S = ω0 2 σ 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 = 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 bonddimension D max = 4.We find that the results obtained using the T-TEDOPA method agree very well with the exact solution (see fig. 4) and correctly reproduce the transition from under-damped to over-damped decay as the temperature is increased [6,38].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 [7].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 Refs.[7,39] 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 Fig. 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 a † ω a ω 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 eq. 9 rather than the original bosonic bath and thus the mode energies run from −ω c to ω c .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 = 0.2ω c ).This agrees well with the renormalized gap ω r 0 = ω 0 (ω 0 /ω c ) α 1−α predicted by the non-perturbative variational polaron theory of Silby & Harris [40], which for the parameters used here gives ω r 0 = 0.167ω 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 nω +1 n−ω = e β 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.
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 [21].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 Optimised Boson Basis (OBB) which was introduced into 1TDVP for open systems by Schroeder et al. [37,41].

IV. 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 [2,3,42].The problem of modeling electron transfer between molecules comes down to accurately treating  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.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 [16,20,43].In many molecular systems -and particularly biological systems where the transfer between electronic states is affected by coupling to chromophore and protein modes -the systembath physics is highly non-perturbative and J(ω) has very sharp frequency-dependence [3,13,44,45].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 Fig. 6.In this model the electron transfer process is modelled 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 = 2 σ z + λ R 1+σx 2 , and the coupling operator as A S = 1+σ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 6(a) 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 J(ω)/ω.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.
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 semiclassical (Born-Oppeheimer) picture, we see that the electron must overcome a kind of effective energy barrier E b that scales with the total reorganisation energy of the entire environment λ R in order for the reaction to progress.We thus might well expect to see thermally activated (exponential) behaviour whereby the tunneling rate ∝ exp(−βE b ).However, at low temperatures this behaviour should be dramatically quenched and dissipative quantum tunneling should become dominant and strongly dependent on the spectral function of the environment [7].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.sec III).The inset shows an enlarged view of the initial fast dynamics which appear to be broadly independent of temperature.

A. 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 6(b) 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 nonexponential behavior is observed for lower temperatures.In figure 7(b) 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 τ ≈ ω −1 c , 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 [46], 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.
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 reorganisation 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 [7].
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 nonmonotonic with a transition from power law growth (quantum, 2α − 1 > 0) to power-law decay (classical, ∝ √ β) as the temperature increases from T = 0. We note that for the parameters we present here, the intermediate regime where thermally activated behaviour is predicted βω c ∼ 1 is not observed for the Ohmic en- vironment, and one essentially switches from tunneling limited by the effect of friction on the attempt frequency to the low-temperature polaronic tunneling of Eqs. 13.

V. 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 [47], 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 [48].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 [16].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 [29], 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 [49], or even the creation of new excitations, such as polaritons [50][51][52], 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.

VI. APPENDIX
The chain mapping used in section II is based on the theory of orthogonal polynomials.A polynomial of degree n is defined by The space of polynomials of degree n is denoted P n and is a subset of the space of all polynomials P n ⊂ P. 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 p, q µ = b a dµ(x)p(x)q(x).
This inner product gives rise to a unique set of orthonormal polynomials {p n ∈ P n , n = 0, 1, 2, ...} which all satisfy pn , pm = δ n,m .
This set forms a complete basis for P, and more specifically the set {p n ∈ P n , n = 0, 1, 2, ...m} is a complete basis for m r=1 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 pn (x) whose leading coefficient is 1 (a n = 1) The key property of orthogonal polynomials for the construction of the chain mapping is that they satisfy a three term recurrence relation where it can be easily shown that Now that we have defined the orthogonal polynomials we may use them to construct the unitary transformation that will convert the star Hamiltonian of Eq. 9 with into the chain Hamiltonian of Eq. 12.The transformation is given by where and the polynomials pn (ω) 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 For the zeroth order monic polynomial we have π 0 = 1 and so we may insert this into the above expression Recognising the inner product in the above expression and making use of the orthogonality of the polynomials we have and thus, in the new basis, only one mode now couples to the system.Now for the environment part of the Hamiltonian we have Substituting for ωπ n (ω) from the three term recurrence relation of Eq. 18 yields Again, evaluating the inner products we have where in the second line we have used the fact that We thus arrive at the nearest-neighbour coupling Hamiltonian of Eq. 12 and are able to identify the chain coefficients as Note that in Eq. 12 the chain sites are labeled starting from n = 1 and not n = 0 as in Eq. 28.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 interatively using Eqs 18 and 19 and numerically evaluating the inner product integrals using a quadrature rule.
a N O p 2 B C 8 5 Z d X S f u i 7 r l 1 7 / 6 y 1 r g p 4 i j D C Z z C O X h w B Q 2 4 g y a 0 g E E E z / A K b 4 5 0 X p x 3 5 2 P R W n K K m W P 4 A + f z B 0 m y j c k = < / l a t e x i t > !k < l a t e x i t s h a 1 _ b a s e 6 4 = " 8 9 h w 5 Z I A g F M A + Q k d j n j O p L I D 3 I 4 = " > A A A B 7 3 i c b V D L S g N B E O y N r x h f U Y 9 e B o P g K e y K o M e g F 4 8 R z A O S J c x O O s m Q m Z 1 1 Z l Y I S 3 7 C i w d F v P o 7 3 v w b J 8 k e N L G g o a j q p r s r S g Q 3 1 v e / v c L a + s b s h a 1 _ b a s e 6 4 = " 8 9 h w 5 Z I A g F M A + Q k d j n j O p L I D 3 I 4 = " > A A A B 7 3 i c b V D L S g N B E O y N r x h f U Y 9 e B o P g K e y K o M e g F 4 8 R z A O S J c x O O s m Q m Z 1 1 Z l Y I S 3 7 C i w d F v P o 7 3 v w b J 8 k e N L G g o a j q p r s r S g Q 3 1 v e / v c L a + s b s h a 1 _ b a s e 6 4 = " 8 9 h w 5 Z I A g F M A + Q k d j n j O p L I D 3 I 4 = " > A A A B 7 3 i c b V D L S g N B E O y N r x h f U Y 9 e B o P g K e y K o M e g F 4 8 R z A O S J c x O O s m Q m Z 1 1 Z l Y I S 3 7 C i w d F v P o 7 3 v w b J 8 k e N L G g o a j q p r s r S g Q 3 1 v e / v c L a + s b s h a 1 _ b a s e 6 4 = " 8 9 h w 5 Z I A g F M A + Q k d j n j O p L I D 3 I 4 = " > A A A B 7 3 i c b V D L S g N B E O y N r x h f U Y 9 e B o P g K e y K o M e g F 4 8 R z A O S J c x O O s m Q m Z 1 1 Z l Y I S 3 7 C i w d F v P o 7 3 v w b J 8 k e N L G g o a j q p r s r S g Q 3 1 v e / v c L a + s b t e x i t s h a 1 _ b a s e 6 4 = " d e H L w Y 3 G 4 O 3 P P d I F 9 3 w s u H 1 9 y D 3 s e i t e D l M 8 f w B 9 7 n D / d N j 1 A = < / l a t e x i t > !< l a t e x i t s h a 1 _ b a s e 6 4 = " K x / h Y x V 7 C K d 3 C I T r H i H D T u r W g P o = " > A A A B 7 n i c b V D L S g N B E O y N r x h f U Y 9 e B o P g x b A r g h 6 D X j x G M A 9 I l j A 7 6 S R D Z m a X m V k h L P k I L x 4 U 8 e r 3 e P N v n C R 7 0 M S C h q K q m + 6 u K B H c W N / / 9 g p r 6 x u b W 8 X t 0 s 7 u 3 v 5 B + f C o a e J U M 2 y w W s h a 1 _ b a s e 6 4 = " 8 9 h w 5 Z I A g F M A + Q k d j n j O p L I D 3 I 4 = " > A A A B 7 3 i c b V D L S g N B E O y N r x h f U Y 9 e B o P g K e y K o M e g F 4 8 R z A O S J c x O O s m Q m Z 1 1 Z l Y I S 3 7 C i w d F v P o 7 3 v w b J 8 k e N L G g o a j q p r s r S g Q 3 1 v e / v c L a + s b s h a 1 _ b a s e 6 4 = " 8 9 h w 5 Z I A g F M A + Q k d j n j O p L I D 3 I 4 = " > A A A B 7 3 i c b V D L S g N B E O y N r x h f U Y 9 e B o P g K e y K o M e g F 4 8 R z A O S J c x O O s m Q m Z 1 1 Z l Y I S 3 7 C i w d F v P o 7 3 v w b J 8 k e N L G g o a j q p r s r S g Q 3 1 v e / v c L a + s b

Figure 1 .
Figure 1.(a)A generic open quantum system contains a fewlevel '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 = β −1 .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.
c < l a t e x i t s h a 1 _ b a s e 6 4 = " 2 e Q r U C h G m 5 c B z 7 D 9 u c l H I S + N 0 X M = " > A A A B + H i c b V D L S g M x F L 1 T X 7 U + O u r S T b A I r u q M C L o s u n F Z w T 6 g H Y Z M m m l D k 8 y Q Z I Q 6 9 E v c u F D E r Z / i z r 8 x b W e h r Q c u 9 3 D O v e T m R C l n 2 n j e t 1 N a W 9 / Y 3 C p v V 3 Z 2 9 / a r 7 s F h W y e Z I r R F E p 6 o b o Q 1 3 7 y 9 r j Z s i j j I c w w m c g Q 9 X 0 I A 7 a E I L C G T w D K / w 5 j w 5 L 8 6 7 8 7 E Y L T n F z h H 8 g f P 5 A 4 V 7 k v w = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " 2 e Q r U C h G m 5 c B z 7 D 9 u c l H I S + N 0 X M = " > A A A B + H i c b V D L S g M x F L 1 T X 7 U + O u r S T b A I r u q M C L o s u n F Z w T 6 g H Y Z M m m l D k 8 y Q Z I Q 6 9 E v c u F D E r Z / i z r 8 x b W e h r Q c u 9 3 D O v e T m R C l n 2 n j e t 1 N a W 9 / Y 3 C p v V 3 Z 2 9 / a r 7 s F h W y e Z I r R F E p 6 o b o Q 1 3 7 y 9 r j Z s i j j I c w w m c g Q 9 X 0 I A 7 a E I L C G T w D K / w 5 j w 5 L 8 6 7 8 7 E Y L T n F z h H 8 g f P 5 A 4 V 7 k v w = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " 2 e Q r U C h G m 5 c B z 7 D 9 u c l H I S + N 0 X M = " > A A A B + H i c b V D L S g M x F L 1 T X 7 U + O u r S T b A I r u q M C L o s u n F Z w T 6 g H Y Z M m m l D k 8 y Q Z I Q 6 9 E v c u F D E r Z / i z r 8 x b W e h r Q c u 9 3 D O v e T m R C l n 2 n j e t 1 N a W 9 / Y 3 C p v V 3 Z 2 9 / a r 7 s F h W y e Z I r R F E p 6 o b o Q 1 3 7 y 9 r j Z s i j j I c w w m c g Q 9 X 0 I A 7 a E I L C G T w D K / w 5 j w 5 L 8 6 7 8 7 E Y L T n F z h H 8 g f P 5 A 4 V 7 k v w = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " 2 e Q r U C h G m 5 c B z 7 D 9 u c l H I S + N 0 X M = " > A A A B + H i c b V D L S g M x F L 1 T X 7 U + O u r S T b A I r u q M C L o s u n F Z w T 6 g H Y Z M m m l D k 8 y Q Z I Q 6 9 E v c u F D E r Z / i z r 8 x b W e h r Q c u 9 3 D O v e T m R C l n 2 n j e t 1 N a W 9 / Y 3 C p v V 3 Z 2 9 / a r 7 s F h W y e Z I r R F E p 6 o b o Q 1 3 7 y 9 r j Z s i j j I c w w m c g Q 9 X 0 I A 7 a E I L C G T w D K / w 5 j w 5 L 8 6 7 8 7 E Y L T n F z h H 8 g f P 5 A 4 V 7 k v w = < / l a t e x i t > J (!) < l a t e x i t s h a 1 _ b a s e 6 4 = " A H e b z y 6 f 4 3 C p 9 H E t l K z F 4 r v 6 e y I n Q e i I i 2 y m I G e p l b y b + 5 3 U z E 9 + E O U v S z E B C F 4 v i j G M j 8 S w G 3 G c K q O E T S w h V z N 6 K 6 Z A o Q o 0 N q 2 x D 8 J d f X i W t y 5 r v 1 f z H q 0 r 9 t o i j h E 7 R G a o i H 1 2 j O r p H D d R E F I 3 R M 3 p F b 0 7 u v D j v z s e i d c 0 p Z k 7 Q H z i f P 0 a j k s I = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " Y c E g B 9 h Z H c z f d G W J f u X g v Q N e f G k = " > A A A B + X i c b V B N S 8 N A E N 3 4 W e t X 1 K O X x S L U S 0 l E 0 G P R i 3 i q Y D + g C W W z n b R L d 7 N h d 1 M o o f / E i w d F v P p P v P l v 3 L Y 5 a O u D g c d 7 M 8 z M i 1 L O t P G 8 b 2 d t f W N z a 7 u 0 U 9 7 d 2 z 8 4 d I + O W 1 pm i k K T S i 5 V J y I a O E u g a Z j h 0 E k V E B F x a E e j u 5 n f H o P S T C Z P Z p J C K M g g Y T G j x F i p 5 7 o P v S A C Q 6 q B F D A g F 7 j n V r y a N w d e J X 5 B K q h A o + d + B X 1 J M w G J o Z x o 3 f W 9 1 I Q 5 U Y Z R D t N y k G l I C R 2 R A X Q t T Y g AH e b z y 6 f 4 3 C p 9 H E t l K z F 4 r v 6 e y I n Q e i I i 2 y m I G e p l b y b + 5 3U z E 9 + E O U v S z E B C F 4 v i j G M j 8 S w G 3 G c K q O E T S w h V z N 6 K 6 Z A o Q o 0 N q 2 x D 8 J d f X i W t y 5 r v 1 f z H q 0 r 9 t o i j h E 7 R G a o i H 1 2 j O r p H D d R E F I 3 R M 3 p F b 0 7 u v D j v z s e i d c 0 p Z k 7 Q H z i f P 0 a j k s I = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " Y c E g B 9 h Z H c z f d G W J f u X g v Q N e f G k = " > A A A B + X i c b V B N S 8 N A E N 3 4 W e t X 1 K O X x S L U S 0 l E 0 G P R i 3 i q Y D + g C W W z n b R L d 7 N h d 1 M o o f / E i w d F v P p P v P l v 3 L Y 5 a O u D g c d 7 M 8 z M i 1 L O t P G 8 b 2 d t f W N z a 7 u 0 U 9 7 d 2 z 8 4 d I + O W 1 p m i k K T S i 5 V J y I a O E u g a Z j h 0 E k V E B F x a E e j u 5 n f H o P S T C Z P Z p J C K M g g Y T G j x F i p 5 7 o P v S A C Q 6 q B F D A g F 7 j n V r y a N w d e J X 5 B K q h A o + d + B X 1 J M w G J o Z x o 3 f W 9 1 I Q 5 U Y Z R D t N y k G l I C R 2 R A X Q t T Y gA H e b z y 6 f 4 3 C p 9 H E t l K z F 4 r v 6 e y I n Q e i I i 2 y m I G e p l b y b + 5 3U z E 9 + E O U v S z E B C F 4 v i j G M j 8 S w G 3 G c K q O E T S w h V z N 6 K 6 Z A o Q o 0 N q 2 x D 8 J d f X i W t y 5 r v 1 f z H q 0 r 9 t o i j h E 7 R G a o i H 1 2 j O r p H D d R E F I 3 R M 3 p F b 0 7 u v D j v z s e i d c 0 p Z k 7 Q H z i f P 0 a j k s I = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " Y c E g B 9 h Z H c z f d G W J f u X g v Q N e f G k = " > A A A B + X i c b V B N S 8 N A E N 3 4 W e t X 1 K O X x S L U S 0 l E 0 G P R i 3 i q Y D + g C W W z n b R L d 7 N h d 1 M o o f / E i w d F v P p P v P l v 3 L Y 5 a O u D g c d 7 M 8 z M i 1 L O t P G 8 b 2 d t f W N z a 7 u 0 U 9 7 d 2 z 8 4 d I + O W 1 p m i k K T S i 5 V J y I a O E u g a Z j h 0 E k V E B F x a E e j u 5 n f H o P S T C Z P Z p J C K M g g Y T G j x F i p 5 7 o P v S A C Q 6 q B F D A g F 7 j n V r y a N w d e J X 5 B K q h A o + d + B X 1 J M w G J o Z x o 3 f W 9 1 I Q 5 U Y Z R D t N y k G l I C R 2 R A X Q t T Y g AH e b z y 6 f 4 3 C p 9 H E t l K z F 4 r v 6 e y I n Q e i I i 2 y m I G e p l b y b + 5 3U z E 9 + E O U v S z E B C F 4 v i j G M j 8 S w G 3 G c K q O E T S w h V z N 6 K 6 Z A o Q o 0 N q 2 x D 8 J d f X i W t y 5 r v 1 f z H q 0 r 9 t o i j h E 7 R G a o i H 1 2 j O r p H D d R E F I 3 R M 3 p F b 0 7 u v D j vz s e i d c 0 p Z k 7 Q H z i f P 0 a j k s I = < / l a t e x i t > = 100!c < l a t e x i t s h a 1 _ b a s e 6 4 = " N n i m g Y 8 6 E p j M u o u T G g L y l G w 4 + x s = " > A A A B + n i c b V B N S 8 N A E N 3 4 W e t X q k c v w S J 4 K h s R 9 C I U v X i s Y D + g D W G z n b R L d 7 N h d 6 O U 2 J / i x Y M i X v 0 l 3 v w 3 b t s c t P X B w O O 9 G W b m R S l n 2 m D 8 7 a y s r q 1 v b J a 2 y t s 7 u 3 v 7 b u W g p W W m K D S p 5 F J 1 I qK B s w S a h h k O n V Q B E R G H d j S 6 m f r t B 1 C a y e T e j F M I B B k k L G a U G C u F b q U X g S F X P s Y 9 K W B A Q h q 6 V V z D M 3 j L x C 9 I F R V o h O 5 X r y 9 p J i A x l B O t u z 5 O T Z A T Z R j l M C n 3 M g 0 p o S M y g K 6 l C R G g g 3 x 2 + s Q 7 s U r f i 6 W y l R h v p v 6 e y I n Q e i w i 2 y m I G e p F b y r + 5 3 U z E 1 8 G O U v S z E B C 5 4 v i j H t G e t M c v D 5 T Q A 0 f W 0 K o Y v Z W j w 6 J I t T Y t M o 2 B H / x 5 W X S O q v 5 u O b f n V f r 1 0 U c J X S E j t E p 8 t E F q q N b 1 E B N R N E je k a v 6 M 1 5 c l 6 c d + d j 3 r r i F D O H 6 A + c z x 8 h u p N C < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " N n i m g Y 8 6 E p j M u o u T G g L y l G w 4 + x s = " >A A A B + n i c b V B N S 8 N A E N 3 4 W e t X q k c v w S J 4 K h s R 9 C I U v X i s Y D + g D W G z n b R L d 7 N h d 6 O U 2 J / i x Y M i X v 0 l 3 v w 3 b t s c t P X B w O O 9 G W b m R S l n 2 m D 8 7 a y s r q 1 v b J a 2 y t s 7 u 3 v 7 b u W g p W W m K D S p 5 F J 1 I q K B s w S a h h k O n V Q B E R G H d j S 6 m f r t B 1 C a y e T e j F M I B B k k L G a U G C u F b q U X g S F X P s Y 9 K W B A Q h q 6 V V z D M 3 j L x C 9 I F R V o h O 5 X r y 9 p J i A x l B O t u z 5 O T Z A T Z R j l M C n 3 M g 0 p o S M y g K 6 l C R G g g 3 x 2 + s Q 7 s U r f i 6 W y l R h v p v 6 e y I n Q e i w i 2 y m I G e p F b y r + 5 3 U z E 1 8 G O U v S z E B C 5 4 v i j H t G e t M c v D 5 T Q A 0 f W 0 K o Y v Z W j w 6 J I t T Y t M o 2 B H / x 5 W X S O q v 5 u O b f n V f r 1 0 U c J X S E j t E p 8 t E F q q N b 1 E B N R N E j e k a v6 M 1 5 c l 6 c d + d j 3 r r i F D O H 6 A + c z x 8 h u p N C < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " N n i m g Y 8 6 E p j M u o u T G g L y l G w 4 + x s = " >A A A B + n i c b V B N S 8 N A E N 3 4 W e t X q k c v w S J 4 K h s R 9 C I U v X i s Y D + g D W G z n b R L d 7 N h d 6 O U 2 J / i x Y M i X v 0 l 3 v w 3 b t s c t P X B w O O 9 G W b m R S l n 2 m D 8 7 a y s r q 1 v b J a 2 y t s 7 u 3 v 7 b u W g p W W m K D S p 5 F J 1 I q K B s w S a h h k O n V Q B E R G H d j S 6 m f r t B 1 C a y e T e j F M I B B k k L G a U G C u F b q U X g S F X P s Y 9 K W B A Q h q 6 V V z D M 3 j L x C 9 I F R V o h O 5 X r y 9 p J i A x l B O t u z 5 O T Z A T Z R j l M C n 3 M g 0 p o S M y g K 6 l C R G g g 3 x 2 + s Q 7 s U r f i 6 W y l R h v p v6 e y I n Q e i w i 2 y m I G e p F b y r + 5 3 U z E 1 8 G O U v S z E B C 5 4 v i j H t G e t M c v D 5 T Q A 0 f W 0 K o Y v Z W j w 6 J I t T Y t M o 2 B H / x 5 W X S O q v 5 u O b f n V f r 1 0 U c J X S E j t E p 8 t E F q q N b 1 E B N R N E j e k a v 6 M 1 5 c l 6 c d + d j 3 r r i F D O H 6 A + c z x 8 h u p N C < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " N n i m g Y 8 6 E p j M u o u T G g L y l G w 4 + x s = " > A A A B + n i c b V B N S 8 N A E N 3 4 W e t X q k c v w S J 4 K h s R 9 C I U v X i s Y D + g D W G z n b R L d 7 N h d 6 O U 2 J / i x Y M i X v 0 l 3 v w 3 b t s c t P X B w O O 9 G W b m R S l n 2 m D 8 7 a y s r q 1 v b J a 2 y t s 7 u 3 v 7 b u W g p W W m K D S p 5 F J 1 I q K B s w S a h h k O n V Q B E R G H d j S 6 m f r t B 1 C a y e T e j F M I B B k kL G a U G C u F b q U X g S F X P s Y 9 K W B A Q h q 6 V V z D M 3 j L x C 9 I F R V o h O 5 X r y 9 p J i A x l B O t u z 5 O T Z A T Z R j l M C n 3 M g 0 p o S M y g K 6 l C R G g g 3 x 2 + s Q 7 s U r f i 6 W y l R h v p v 6 e y I n Q e i w i 2 y m I G e p F b y r + 5 3 U z E 1 8 G O U v S z E B C 5 4 v i j H t G e t M c v D 5 T Q A 0 f W 0 K o Y v Z W j w 6 J I t T Y t M o 2 B H / x 5 W X S O q v 5 u O b f n V f r 1 0 U c J X S E j t E p 8 t E F q q N b 1 E B N R N E j e k a v6 M 1 5 c l 6 c d + d j 3 r r i F D O H 6 A + c z x 8 h u p N C < / l a t e x i t > = 4!c < l a t e x i t s h a 1 _ b a s e 6 4 = " 5 8 z E h J / C z 2 R 9 J i Q A W q E 2 T M P T B R c = " > A A A B + H i c b V B N S 8 N A E N 3 U r 1 o / W v X o Z b E I n k o i g l 6 E o h e P F e w H N C F s t p N 2 6 e 4 m 7 G 6 E G v p L v H h Q x K s / x Z v / x m 2 b g 7 Y + G H i 8 N 8 P M v C j l T B v X / X Z K a + s b m 1 v l 7 c r O 7 t 5 + t X Z w 2 N F J p i i 0 a c I T 1 Y u I B s 4 k t A 0 z H H q p A i I i D t 1 o f D v z u 4 + g N E v k g 5 m k E A g y l C x m l B g r h b W q H 4 E h 1 x d + I m B I Q h r W 6 m 7 D n Q O v E q 8 g d V S g F d a + / E F C M w H S U E 6 0 7 n t u a o K c K M M o h 2 n F z z S k h I 7 J E P q W S i J A B / n 8 8 C k + t c o A x 4 m y J Q 2 e q 7 8 n c i K 0 n o j I d g p i R n r Z m 4 n / e f 3 M x F d B z m S a G Z B 0 s S j O O D Y J n q W

Figure 2 .
Figure 2. (a) The extended proxy environment of Fig.1 (a) 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 neighbour interactions.The non-perturbative wave function dynamics for such a manybody 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 c † n cn at time ωct = 45 for baths of several temperatures.The system, which in this case is the Ohmic SBM, with ω0 = 0.2ωc and α = 0.1, is attached at site n = 1 of the chain.

Figure 5 .
Figure 5. Bath mode occupations nω = a † ω aω for the extended environment after the TLS has decayed.The TLS is governed by a Hamiltonian HS = ω 0 2 σz where ω0 = 0.2ωc and is coupled to an Ohmic bath with a hard cut-off via AS = σx.The coupling strength is α = 0.1.Left inset: total mode occupation as a function of time n tot = ∞ −∞ dω nω .Right inset shows nω p +1 nω n plotted on a log scale against the inverse temperature, demonstrating the detailed balance of the absorption and emission rates.

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.

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.sec III).The inset shows an enlarged view of the initial fast dynamics which appear to be broadly independent of temperature.