Photo-induced dynamics with continuous and discrete quantum baths

.


I. INTRODUCTION
The capability to simulate the ultrafast quantum dynamics of photophysical processes in complex molecules is essential to understanding phenomena such as exciton transfer and energy conversion, whose practical control is the foundation of the design of functional compounds with fascinating properties.Their possible applications range from tailored chemical catalysts to synthetic light-harvesting materials and solar cells with efficiencies beyond the Shockley-Queisser limit [1][2][3][4][5] .However, the corresponding computational problems involve Markovian embedding methods FIG. 1.Current methods for simulating the quantum dynamics of excited states coupled to vibrational degrees of freedom target different types of environments.Non-Markovian methods are well-suited to describe continuous spectral densities (top-left), while unitary methods work best for discrete spectral densities (top-right), which capture long-memory effects in the environment.We introduce a pure-state unraveled hybrid-bath scheme using a Markovian-embedding (bottom) to efficiently treat both types of environments.
the dynamics of excitonic states coupled to vibrational modes, each of which is a challenging problem on its own.Exploiting various simplifications, the computational task can be reduced to faithfully describe excitonic degrees of freedom (DOFs) representing electronic excited states, which are (non-locally) coupled to a manifold of vibrational DOFs.Here, especially the vibrational system of complex molecules poses various numerical challenges and there has been a great effort to deal with the exponential scaling of the computational complexity.In general, two conceptually different approaches were established.(i) A discretization of the vibrational DOFs is combined with selection rules to choose the most relevant modes, such that the unitary dynamics of the resulting, closed quantum system can be described using efficient state representations and integration schemes [6][7][8][9][10] .(ii) The excitonic DOFs are treated exactly while the vibrational DOFs are traced out and encoded in an effective environment.The complexity of the resulting, open quantum-system dynamics then mainly depends on the choice of the environment 11,12 .The dynamics in the presence of memoryless, Markovian environments can be treated efficiently using Lindblad master equations [13][14][15][16][17] .Structured, non-Markovian environments introduce time-non-local correlations between system and environment, which increase the computational cost [18][19][20][21][22][23][24][25][26][27][28][29] but are crucial to simulate the inherently mesoscopic vibrational systems of molecules, exhibiting long-time memory effects.
To efficiently simulate the quantum dynamics of both, coupled electronic and vibrational, as well as environmental DOFs, the decomposition of the many-body wave function into local tensors has proven to be an extremely fruitful approach [30][31][32][33][34][35][36][37][38] .The unitary evolution of excitonic states coupled to O(100) discrete vibrational DOFs is a remarkable milestone facilitating, for instance, the unbiased investigation of ultrafast photo-induced dynam-ics of large molecules [39][40][41] .Recently, wave functionbased tensor network (TN) methods for open quantumsystems with both Markovian and non-Markovian environments have been developed [42][43][44] to describe continuous and smooth spectral densities.
However, these approaches leave a methodological gap, if the system's spectral density contains both aspects, discrete peaks immersed into a smooth background.This situation, generated for instance by non-vanishing memory effects, is inevitable in realistic materials, such as organic semiconductors, which involve both dispersive modes with continuous spectral densities and local vibrational modes with sharp peaks.While in these cases the unitary approach would require a daunting amount of discrete vibrational modes, generic open-quantum system methods become unstable in the presence of discrete peaks 45 .To fill this gap, pseudomode approaches have been introduced to approximate the spectral density by introducing damped auxiliary modes 22,23,[46][47][48] and solving the equations of motion for the density operator 49,50 .Here, we propose a TN-based pure-state unraveled hybrid-bath (PS-HyB) scheme computing quantum trajectories in terms of matrix-product state (MPS) of an enlarged system with Markovian embedding, where both continuous and discrete contributions to the spectral density can be handled (see Fig. 1).Simulating two paradigmatic examples, i.e., singlet fission in a threestate model with Debye spectral density, as well as finitetemperature exciton-dynamics in a Fenna-Matthews-Olson (FMO) complex with an additional sharp peak in the spectral density, we demonstrate that PS-HyB-MPS can efficiently capture the regime of long-time memory effects at significantly smaller computational cost than unitary methods, while it circumvents instabilities from a direct evolution of the density operator.This is achieved by exploiting recent advances for TN representations of systems with large local Hilbert spaces, rendering a purestate unraveling practically feasible, also in the presence of strong exciton-phonon couplings 40,[51][52][53][54] .

II. METHODOLOGY
We aim to describe the dynamics of L excitonic states, each interacting with a set of vibrational DOFs, acting as independent environments.For this system-environment partitioning, the spectral density J I (ω) (I = 1, 2, . . .L) is the key quantity encoding the relevant information about the coupling between excitonic and the vibrational DOFs.Our goal is to capture the general situation of a continuous background and discrete, sharp peaks describing long-time memory effects, as pictured in Fig. 1.For that purpose, we adopt a method called mesoscopic leads or extended reservoirs [22][23][24][25][26][27][28] .In this framework, the general idea is to approximate each bath spectral density ) by fitting it with N I Lorentzians: where the couplings between excitonic and vibrational DOFs are given by g I m = 2J I (ω I m )γ I m /π, .Each Lorentzian contributing to the spectral density can be represented by a harmonic oscillator with frequency ω I m coupled to a memoryless environment, where the spectral broadening is given by the potentially non-uniform discretization γ I m .This way, we can represent a continuous spectral density with sharp features by a set of I N I bosonic modes coupled to Markovian environments.

A. Pure-state unraveled hybrid-baths
Let us now work with the explicit Hamiltonian for the combined excitonic and vibrational system Ĥ = Ĥexc + Ĥvib + Ĥex−vib . ( Here, Ĥexc is an arbitrary excitonic Hamiltonian.In the following, we set k B , ℏ = 1.The vibrational modes are described by harmonic oscillators with frequencies ω The sum on the right-hand side models the interaction of the system with a Markovian environment.In our case, these terms describe the thermalization of each vibrational subsystem I via the Davies map 55,56 where the thermalization rate is proportional to the spacing and f I m (β) = 1/(exp βω I m − 1) is the Bose-Einstein distribution function at inverse temperature β.Here, we introduced the dissipator specifying the couplings to the environment via the jump operators L1 = bI † m and L2 = bI m .For a large number of bath modes N I , the dissipative evolution of the excitonic subsystem Eq. (3) becomes equivalent to the unitary evolution with discretized bath spectral densities 57 , described for instance via a chain mapping 28 .We will show the mesoscopic leads approach to be numerically much more efficient when solving the dynamics using TN methods.This is based on the observation that one obtains the unitary case from taking the zero-broadening limit of Eq. (3).In fact, when J(ω) is continuous and smooth ( Fig. 1(top left)) a large number of discrete modes would be required to accurately fit J(ω).Moreover, unitary dynamics suffer from severe finite-size errors caused by the scattering of excitations at the systems' boundaries, which limits the reachable simulation time 28 .
Nonetheless, compared to unitary time evolution, Lindblad dynamics pose challenges to numerical treatments not only because of the dimension but also because the truncation of a mixed state requires particular care 58,59 .Thus, it is convenient to resort to so-called pure-state unravelings, which consist of reconstructing the Lindblad dynamics from an ensemble of stochastic pure-state time evolutions.Different pure-state unraveling schemes, such as quantum state diffusion (QSD) 60 and the quantum jumps (QJ) 16,17 approaches, have been developed.Here, we apply the QJ scheme, which is the more widely used variant due to its good stability and convergence properties.The starting point for the QJ scheme is to rewrite Eq. (3) as where Ĥeff = Ĥ − i/2 l L † l Ll is an effective, non-Hermitian Hamiltonian.If the initial state is pure (ρ(0) = |Ψ(0)⟩ ⟨Ψ(0)|) the first term of Eq. (C6) describes a simple (non-unitary) Schrödinger evolution.The action of the second term can be mapped to the stochastic application of Ll to a pure state.We can achieve this by introducing a collection of random numbers defining trajectories q(t) (see for instance 16) such that at any time t the ensemble average E over N T pure-state evolutions yields an approximation of the Lindblad-evolved density matrix E N T [|Ψ(t) q ⟩ ⟨Ψ(t) q |] = ρN T (t) (which formally becomes exact for an infinite number of trajectories, i.e. ρN T (t) ). Usually, one avoids computing the density matrix explicitly and evaluates the expectation value of an observable of interest directly as The quick convergence of local observables in the number of trajectories is shown in Fig. 2(top) as the example of a minimal model for singlet fission, which is discussed in more detail in Sec.III.Using the independence of different trajectories, one can estimate the statistical error ϵ O N T for the expectation value of an observable Ô as 16 We will adopt Eq. ( 7) as error bars for all QJ results.
We note that dissipation-assisted matrix-product factorization (DAMPF), an approach similar to the one discussed in this manuscript, was introduced in Ref. 49. DAMPF is based on fitting the bath correlation function (which is connected to J(ω) by a Fourier transform) in the time domain and on the direct propagation of the density matrix.While the direct propagation of the density matrix using DAMPF avoids the sampling error present in the QJ method, it can suffer from positivity issues 61 and the evolution of matrix-density operators is in general numerically less efficient than advanced time-evolution schemes for MPS, which we use to compute the QJ trajectories.Essentially, the advantage of PS-HyB-MPS is based on the computationally beneficial properties of MPS, which are made accessible by the pure state unraveling and allow for a taylored MPS representation of the vibrational degrees of freedom and to invoke high-performance time-evolution schemes for MPS, which we discuss in Sec.II B.

B. Projected-purified matrix-product statess
We represent the excitonic-vibronic system as a MPS 6,8 .In general, this allows to employ efficient truncation schemes for |Ψ⟩ that are controlled by the bond dimension and/or the truncated weight.However, representing the vibrational modes as an MPS is an arduous task because one has to truncate the formally infinite bosonic Hilbert space.Moreover, the interaction Hamiltonian breaks the U (1) symmetry for the bosonic excitations.We can address both challenges by working with projected-purified matrix-product states (PP-MPS) 51,52 , which were successfully applied to simulate unitary photo-induced dynamics in different setups 40,53,54 .
The idea of PP-MPS is to pair every physical bosonic degree of freedom |n m ⟩ with a fictitious bath mode |n B;m ⟩, where n m , n B;m ∈ {0, . . ., d − 1} and d denotes the cutoff dimension of the bosonic Hilbert spaces.The exchange of excitations between physical and bath degrees is constrained by imposing local gauge constraints n m + n B;m = d − 1, which can be satisfied easily when preparing the initial state of the time evolution.In order to conserve the local gauge constraints throughout the dynamics, bosonic excitations of the physical degrees of freedom have to be moved from or to the bath sites.This is easily achieved by introducing balancing operators β( †) B;m that counteract the change of excitations in the physical, bosonic system.The local gauge constraints have various advantageous properties.They introduce an artificial, global U (1) symmetry because for each pair of physical and bath degree of freedom, the total number of excitations is conserved, which allows the employment of efficient, block-diagonal representations of site tensors.As a consequence, the representation of a local bosonic operator is transformed from one d-dimensional block to d one-dimensional blocks and one can exploit efficient (Top) The convergence of PS-HyB-MPS with the number of trajectories at the example of the CT-state occupation for the minimal model of singlet fission studied in Sec.III.(Middle and bottom) We compare the convergence of occupation of site 1 in the FMO complex at finite temperatures using either a naive pure-state unraveling of the thermal density operator (middle) or the thermofield and thermal vacuum representation described in Sec.II C. Note that when using the naive pure state unraveling, we observe a very slow convergence, requiring extremely small time steps δt ≤ 0.5fs.Instead, the thermofield and thermal vacuum representation exhibits a significantly increased stability, allowing for computationally beneficial, larger time steps δt ≥ 5fs.
parallelization schemes over different quantum number sectors.Even more importantly, it can be shown that under the introduction of the gauge constraints, the diagonal elements of the bosonic single-site reduced density-matrix (1RDM) ρ n m ,n ′ m are directly related to the Schmidt values when introducing a bipartition between site m and m + 1.This allows to directly truncate at the level of the bosonic 1RDM by using only standard density matrix renormalization group (DMRG) truncation schemes.In fact, if the probabilities ρ nm,nm to have n m excitations at vibrational mode m are below a certain truncation threshold δ > ρ nm,nm , discarding these matrix elements implies an actual reduction of the local Hilbert space dimension d → d eff 51 .In practice, this allows us to choose d rather large since the required local Hilbert space dimension d eff is controlled by tuning the truncation threshold δ of excitation probabilities for the physical, bosonic degrees of freedom.The resulting truncation of the local Hilbert space dimension is similar to the effect of a local basis optimization 62,63 .However, in PP-MPS, there is no need for an additional basis transformation on the local bosonic Hilbert spaces, which alleviates typical problems of dynamical basis adaptions and truncations, such as reintroducing previously discarded basis states 52 .
To calculate the time-evolution of a MPS, one of the most efficient methods is a TN variant of the time-dependent variational principle (TDVP) 9,10,64 .It is based on the Dirac-Frenkel variational principle 65 and yields an approximation to the solution of the Schrödinger equation projected to the tangent space of the variational manifold generated by the state parametrization (here the MPS tensor elements).The numerical efficiency of TDVP in its TN formulation stems from the fact that the projector to the tangent space of the MPS manifold can be split into local projectors and thereby induces a local integration scheme that can be solved efficiently 9 .This comes at the cost of projection errors typical for TN algorithms, which originates from the restriction of the local state space and can introduce uncontrolled errors.Consequently, convergence with respect tp the MPS bond dimension is crucial, and a careful analysis for the systems simulated in this article can be found in 66 .
It should be noted that the so-called projectorsplitting constitutes the major difference between classical multi-layer multi-configurational time-dependent Hartree (ML-MCTDH) implementations and the TN formulation of TDVP.Moreover, a projector-splitting for ML-MCTDH 67 was introduced recently avoiding the inversion of potentially ill-conditioned transformation matrices 68,69 .The remaining conceptual difficulty is how to deal with projection errors, which are generated by lowest-rank approximations of the wavefunction.For bosonic degrees of freedom, we extended the single-site variant of TDVP with a local subspace expansion 70,71 , which exhibits the computationally most favorable scaling to resolve this issue.It combines single-site TDVP updates with local basis expansions, increasing the states' bond dimension and thereby reducing projection errors such that they are not the dominant source of computational errors.The resulting expansion scheme captures the spread of correlations during the time evolution, while the dominating contractions scale only quadratically in the cutoff dimension of the local bosonic Hilbert spaces.

C. Thermofield doubling and thermal vacuum representation
To describe vibrational environments at finite temperatures, we perform a thermofield doubling of the bosonic system and map the bosonic thermal state to the vacuum via a Bogoljubov transformation, as detailed in Ref. 72.This allows for an efficient representation of thermal fluctuations, if the vibrational modes are described by independent harmonic oscillators and the coupling to the ex-citonic system is linear.Here, the idea is to introduce a new set of bosonic operators âk , linked to the original operators bk by a unitary transformation, which annihilate the thermal state ρβ vib .By transforming the Hamiltonian to the new basis, the thermal state can be represented as a vacuum state.This not only avoids thermal sampling of initial states, but also most importantly, reduces the number of quantum jumps from the environment to the vibrational modes.Thereby, the convergence of the QJ method is strongly improved, which is shown exemplarily in Fig. 2 (middle and bottom), where we also compare the convergence against a sampling of the thermal state using an additional pure-state unraveling of the density operator 73 .

III. EXAMPLE 1: SINGLET FISSION
To benchmark the capability of PS-HyB-MPS to reduce the necessary complexity of a faithful description of vibrational subsystems, we first studied a paradigmatic, minimal model for singlet fission at zero temperature 31,34,53 .The excitonic part of the Hamiltonian consists of L = 3 excitonic states usually referred to as the (i) locally excited (LE), (ii) charge transfer (CT), and (iii) two-triplet (TT) states.Each state is coupled to its own bath of harmonic oscillators characterized by a Debye spectral density We set λ = 807cm −1 and ω c = 1452cm −1 and draw the explicit model parameter of Ĥexc from 34 .The details of the numerical implementation together with a convergence analysis can be found in App. A. Following Sec.II, for each excitonic state we fit the corresponding spectral density J D (ω) and we found N I = 10 vibrational modes to be sufficient to faithfully reproduce J D (ω).At t = 0, we prepared the excitonic system in the LE state and the vibrational environment in the T = 0K vacuum state.Choosing a timestep δt = 5 fs and a cutoff dimension d = 16 for the bosonic sites, we evaluated the time-dependent excitonic populations, averaging over 1000 independent trajectories.In Fig. 3a, we show the resulting dynamics (solid lines, error bars indicated every 50fs) compared to reference data (crosses) extracted from Ref. 34.We find an excellent agreement also at very long simulation times using only N I = 10 bosonic modes per excitonic state, which should be compared to N I = 61 modes per excitonic state used in the unitary time-evolution scheme in Ref. 34.Small deviations can be attributed to the number of sample trajectories, which is discussed in more detail in 66 .In the insets, we display the fitting of the spectral density for the example of FMO.We find excellent agreement with the reference data using a rather small number of Lorentzians N I = 10 for (a) and N I = 19 for (b) and (c).

IV. EXAMPLE 2: LIGHT-HARVESTING COMPLEX
Simulating complex excitonic dynamics at finite temperatures is an even more challenging problem, particularly when using wave function-based methods.As a second example, we therefore studied the photo-induced dynamics of the FMO-complex at finite temperature.It is a pigment-protein complex found in green sulfur bacteria, which is involved in the early stages of photosynthesis 74 .
Following the literature, we used L = 7 excitonic states for the FMO complex.Each state is coupled to an independent vibrational environment 12 , which we initialized at a finite temperature 72 (see App. D).The excitonic Hamiltonian is given by Ĥexc where we adopted the couplings t IJ from Ref. 12.The standard approach to describe the vibrational system is by assuming a Debye spectral density Eq. ( 8), which we fit with N I = 19 Lorentzians per vibrational environment.We extend this approach by adding a sharp peak at a frequency ω 0 = 29cm −1 , which is generated from a vibrational mode of the individual pigments.We emphasize that this modification introduces long-time memory effects in the environment that can significantly alter the photo-induced dynamics.Moreover, previous studies that used a non-Markovian description of the environments, only considered finite broadenings 12 of this mode, drastically reducing memory-effects.We performed simulations of the photo-induced dynamics at the practically relevant temperatures T = 77 and 300K measuring the time-dependent exciton population, starting from an initial exciton configuration n 1 = 1, n I>1 = 0, and using 1000 independent trajectories.In Fig. 3b, we show the resulting dynamics for the case of a pure Debye spectral density, i.e., without the additional peak at ω 0 , at T = 300K.The solid lines represent the time-dependent excitonic state occupations and the error bars are now indicated every 100fs.We again find excellent agreement with the reference data (crosses), reaching a very long simulation time of 1ps.In the inset, we show the quality of the fit of the Debye spectral density Eq. ( 8), using N I = 19 Lorentzians.We checked that the small oscillations around the exact shape of J D (ω) have no impact on the obtained dynamics.
Taking into account the additional δ-peak in J D (ω), we performed computations at T = 77K.The resulting, time-dependent excitonic occupations are depicted in Fig. 3c.Note the additional δ-peak in the spectral density indicated by the grey dashed line of the inset.The effect of the long-term memory induced by the dissipationless additional mode immediately translates into the excitonic dynamics.Compared to the Debye-only case, we find long-lived oscillations in the excitonic occu-  Convergence of the LE population ⟨LE|ρ|LE⟩ for the SF-example obtained from unitary evolution using N I = 30 (red dashed) and N I = 61(black dotted) modes, representing the most relevant vibrational modes 41 compared to the dynamics obtained from PS-HyB-MPS.For the case of N I = 30 modes, the unitary dynamics exhibit stronger oscillations caused by finite-size effects.Note that using PS-HyB-MPS, the unitary dynamics with N I = 61 modes is wellreproduced but with a significantly smaller number of bosonic modes N I = 10.In the inset, we show a runtime comparison of the unitary dynamics (orange bar) with PS-HyB-MPS (blue bars).The computational speed up is nearly an order of magnitude, due to the drastically reduced number of bosonic modes required.
pations, with an increased coherence time compared to Ref. 12, due to the exact treatment of the δ-peak.

V. BENCHMARK
Besides ensuring the consistency of PS-HyB-MPS, we also performed benchmark computations to compare the efficiency of PS-HyB-MPS to that of standard, unitary dynamics, here at the example of SF( Sec.III) 75 .In Fig. 4, we show the dynamics of the LE state (solid lines), computed using different numbers of N I = 5, 10, 15, 20 Lorentzians to approximate Eq. ( 8).The unitary dynamics was calculated taking into account N I = 30 (dashed line) and N I = 61 (dotted line) relevant vibrational modes where it has been shown that N I = 61 is required to faithfully describe the excitonic dynamics 41 .Due to the discrete nature of the vibrational environment in the unitary approach, finite-size effects are observed that translate to stronger oscillations in the initial dynamics.Instead, PS-HyB-MPS reproduces both, the initial and long-time dynamics but with a much smaller number of bosonic modes.This translates into a significant reduction of the runtime, which is shown in the inset 76 .Using only N I = 10 modes, the observed computational speedup is a factor of ∼ 7.5.We also note that the finite lifetime of the bosonic excitations generated by the dissipative evolution in PS-HyB continuously decreases both, the bond and effective bosonic Hilbert space dimensions, in the scope of the dynamics 77 .

SUMMARY AND OUTLOOK
In this work, we introduced PS-HyB-MPS, a method that allows us to efficiently describe the influence of both continuous and discrete vibrational quantum baths on exciton dynamics.To this end, we combine a Markovian-embedding scheme with a pure-state unraveling and state-of-the-art tensor network techniques to represent and evolve bosonic baths with very large Hilbert space dimensions.PS-HyB-MPS is capable of describing both short-and long-time memory effects in exciton dynamics at finite temperatures while requiring a much smaller amount of bosonic modes compared to commonly used unitary descriptions.The resulting reduction in computational complexity translates into significantly lower numerical costs; for instance, we obtain a speedup of nearly an order of magnitude for a paradigmatic model of singlet fission.We believe that PS-HyB-MPS will be a valuable tool to tackle the description of intricated excitonic dynamics in vibrational environments, opening the possibility to describe mesoscopic environments efficiently by using wave function-based time-evolution schemes.Moreover, the prospect of overcoming the current limitation of independent vibronic baths using non-Markovian descriptions of the environment (see App. E) opens the possibility of increasing the computational efficiency of more complex problems, such as photo-induced excitonic dynamics in the presence of coupled vibrational modes.modes at every time evolution step.For all results presented in the main text, we allowed for a maximum local dimension d max of = 16.For the time evolution, we used the TDVP algorithm 9,64 .Specifically, we employed the local subspace expansion time-dependent variational principle (LSE-TDVP) 70,71 , which is very efficient, especially for systems with large local dimensions, like the vibrational modes in our case.This consists of single-site TDVP steps with some maximal bond dimension m alternated by expansions that increase the state's bond dimension by m e .Unless otherwise stated, we set m = 50 and m e = 10.For the TN calculations, we used the SyTen toolkit 78,79 .We unraveled the Lindblad equation using the QJ method (see also App.C) with a time step of dt = 5 fs and either N trajs = 500 or N trajs = 1000 trajectories (see the figure captions).The QJ method was implemented in evos, a Python package for open quantum system dynamics that uses SyTen as the backend.All computations were performed on an Intel Xeon Gold 6130 CPU (2.10GHz) with 16 threads.For all singlet fission calculations, we fitted the spectral density using an equally spaced discretization.Instead, for the FMO dynamics we found it more convenient to use a non-equally-spaced discretization, as shown in table I.In the Markovian embedding scheme, the mesoscopic leads method [22][23][24][25][26][27][28] decomposes the spectral densities into Lorentzians (with center ωm and width γ).In analogy to Eq. (C1), the coupling between an excitonic mode I and a vibrational mode with frequency ωm is given by where the Lorentzian widths are chosen to be γ = ωm+1 − ωm = ∆ω.The vibrational system is embedded in Markovian environments characterized by a dissipation strength γ.Therefore, if the same discrete frequencies are chosen, these two schemes evolve the same excitonic-vibrational system either unitarily as a closed system where the rest of the continuous bath is completely ignored (namely, the broadening is completely absorbed into the discrete peak) or as an open system with the Markovian embedding.
In the following, we compare the unitary dynamics to the dissipative Markovian-embedding dynamics unraveled with QJ.The evolution of a closed system with initial pure state |Ψ(0)⟩ obeys the Schrödinger equation where Ĥ is the Hermitian Hamiltonian of the enlarged system discussed above (exctions and vibrational modes).The induced time-dependence is given by and in practical numerical computations, one typically splits the unitary time evolution operator into several time steps In the Markovian embedding scheme, the evolution of the enlarged system obeys the Lindblad master equation where Ĥeff = Ĥ − i/2 l L † l Ll ≡ Ĥ − i Ĥ′ .The QJ method unravels the evolution of the density matrix ρ(t) into an ensemble of quantum trajectories 16 .Each trajectory follows a non-unitary evolution, accompanied by quantum jumps happening at random times.At every timestep, with a probability of 1 − δp, a trajectory either evolves under a non-Hermitian Hamiltonian The second method consists in unraveling the thermal state ρβ vib = k p(k) |ψ k ⟩ vib ⟨ψ k | vib , by sampling pure states |ψ k ⟩ vib according to the probability distribution p(k).This additional sampling of the initial state, which adds to the sampling of the QJ method, increases the total number of trajectories required to converge the simulations.Since the vibrational degrees of freedom are described by independent harmonic oscillators, the probability p I m (n) of finding n vibrational excitations for the m-th mode relative to the I-th exciton is In Fig. 9 we compare the performance of the thermofield doubling against the unraveling of the thermal state.We simulate the occupation of the first excitonic site of the FMO at 77K using 50 trajectories for different timesteps dt.It can be seen that while the thermofield method yields good results for dt = 5fs the unraveling method is not converged even for a ten times smaller timestep.This comes from the fact that at 77K the probability of quantum jumps happening between the thermal environment and the vibrational modes is much higher than at zero temperature, and the increased number of jumps strongly slows down the real-time evolution.Therefore, the mapping of the thermal state to a vacuum state 72 is crucial.consider the interaction Hamiltonian studied in the main text (although more general interactions are possible 21 ): where g is the coupling strength, n is the excitonic number operator, and b † ( b) creates (annihilates) one vibrational mode.Differently from the mesoscopic leads approach, one describes the system-bath interaction by working in the time domain and considering the bath correlation function Here J I (ω) is the frequency-dependent bath spectral density, the superscript I labels the different baths, and β is the inverse bath temperature.Instead of of fitting J I (ω) with Lorentzians, we approximate α I (t − t ′ ) with a sum of complex exponentials, for instance via the Laplace-Pade method 81 , as: where γ and ω are the damping rate and the vibrational frequency, respectively.
In the following, we will denote a pure state on the combined exciton-phonon system by |Ψ⟩ and for the excitons only with |ψ⟩.In analogy to the Markovian case 60 , it is possible to formulate a stochastic equation of motion for pure-state trajectories q (see main text), so that the ensemble average E[|ψ q (t)⟩ ⟨ψ q (t)|] = ρex (t) (E5) yields the correct time-evolved mixed state ρex for the excitons.
Here, z(t) is a time-dependent complex colored noise with mean zero and correlations E[z(t)z * (t ′ )] = α(t − t ′ ) and E[z(t)z(t ′ )] = 0 82  An efficient strategy for solving Eq. (E6) was put forward by Suess, Eisfeld, and Strunz, who introduced the HOPS method 21 .The main idea of HOPS consists of approximating the time non-local Eq. (E6) with a hierarchy of timelocal equations, which are simpler to solve.In 43,44 , the authors showed that these hierarchies can be mapped to bosonic modes, facilitating their representation as MPS.For every trajectory q, the time evolution for the state on the combined excitonic-auxiliary moded |Ψ⟩ reduces to where Ĥq eff (t) is a non-Hermitian, time-dependent effective Hamiltonian 42,43 .
In 17 , a systematic comparison of a Markovian-embedding QJ description with HOPS was carried out for the Hubbard-Holstein model, which describes the interaction of electrons with locally-coupled Einstein phonons.Interestingly, HOPS was found to perform best in the presence of strong dissipation (i.e.strong damping γ), where the non-Markovian effects and thus the hierarchy depth are strongly suppressed.Instead, QJ was the method of choice in the case of weak dissipation, which is because for γ → 0 it reduces to simple unitary dynamics for the enlarged system.Based on this, we believe that in the context of optically-induced exciton dynamics, HOPS can prove to be beneficial compared to QJ when large damping rates γ are considered, i.e. either when few vibrational modes are used to discretize the bath or for large bath temperatures.A comparison between the mesoscopic leads and the HOPS methods for these setups, as well as their extension to correlated baths, will be the topic of future investigations.

2 .
annihilation (creation) operators.The interaction between excitons and vibrational modes takes the form Ĥex−vib = m,I g I m nI ( bI m + bI † m )/ √ Upon introducing a sufficient number of oscillators N I , the dynamics of the combined system obeys the Lindblad equation 13,14 dρ dt = −i Ĥ, ρ + I D I (ρ) .
FIG. 2.(Top) The convergence of PS-HyB-MPS with the number of trajectories at the example of the CT-state occupation for the minimal model of singlet fission studied in Sec.III.(Middle and bottom) We compare the convergence of occupation of site 1 in the FMO complex at finite temperatures using either a naive pure-state unraveling of the thermal density operator (middle) or the thermofield and thermal vacuum representation described in Sec.II C. Note that when using the naive pure state unraveling, we observe a very slow convergence, requiring extremely small time steps δt ≤ 0.5fs.Instead, the thermofield and thermal vacuum representation exhibits a significantly increased stability, allowing for computationally beneficial, larger time steps δt ≥ 5fs.

FIG. 3 .
FIG.3.Dynamics of the excitonic populations ρII = ⟨I|ρ|I⟩ for a three-state model of SF (a) at T = 0K and a photo-excited FMO complex at T = 300K (b) and T = 77K (c), using a Debye spectral density, which in (c) is extended by a vibrational mode of the FMO pigments.The solid lines are obtained using PS-HyB-MPS with a QJ pure state unraveling using N = 1000 independent trajectories; error bars are shown only every 50fs (a) and 100fs (b), (c).Reference values are denoted by crosses, extracted from Ref. 34 for (a) and Ref. 12 for (b).In the insets, we display the fitting of the spectral density for the example of FMO.We find excellent agreement with the reference data using a rather small number of Lorentzians N I = 10 for (a) and N I = 19 for (b) and (c).
FIG. 4.Convergence of the LE population ⟨LE|ρ|LE⟩ for the SF-example obtained from unitary evolution using N I = 30 (red dashed) and N I = 61(black dotted) modes, representing the most relevant vibrational modes 41 compared to the dynamics obtained from PS-HyB-MPS.For the case of N I = 30 modes, the unitary dynamics exhibit stronger oscillations caused by finite-size effects.Note that using PS-HyB-MPS, the unitary dynamics with N I = 61 modes is wellreproduced but with a significantly smaller number of bosonic modes N I = 10.In the inset, we show a runtime comparison of the unitary dynamics (orange bar) with PS-HyB-MPS (blue bars).The computational speed up is nearly an order of magnitude, due to the drastically reduced number of bosonic modes required.

FIG. 5 .
FIG. 5. Population of TT state in Markovian embedding dissipative dynamics(equally spaced N I = 10) with different numbers of trajectories.

FIG. 6 .FIG. 7 .
FIG.6.Total phononic occupation in the bath of each site(N I ph = k ⟨ bI † k b I k ⟩) using equally spaced N I = 30 unitary evolution (upper panel), equally spaced N I = 6 unitary evolution (middle panel) and equally spaced N I = 6 Markovian embedding scheme (lower panel).The phononic occupations of the unitary evolutions are larger than the ones of the Markovian embedding scheme by one order of magnitude.Here, solid lines indicate positive-frequency modes and dashed lines for negative-frequency modes in the thermofield doubling scheme 72 (see App. D).

FIG. 8 .
FIG.8.Local dimension d convergence of the unitary evolution (upper) and Markovian embedding (lower) using equally spaced N I = 30.We show the populations of site 1 (black, left panels) and site 3 (blue, right panels).The bond dimensions of dissipative dynamics are fixed to m = 10.

FIG. 10 .
FIG.10.Two equivalent descriptions of exciton dynamics in a vibrational environment.Non-Markovian methods like HOPS consider only the exciton subsystem explicitly, while Markovian-embedding techniques as the mesoscopic leads method include also the vibrational modes.

α
I (t − t ′ ) ≈ N m g I m e γ I m |t−t ′ |−iω I m (t−t ′ ) ,(E4) . The time evolution for each trajectory is described by the so-called non-Markovian quantum state diffusion (NMQSD)∂ t |ψ(t)⟩ = −i Ĥex |ψ(t)⟩ + m,I g I m z I * m (t)n I |ψ(t)⟩ − δz I * m indicatesa functional derivative.In practice, solving this equation directly is unfeasible because the third term is time non-local.