Quantum dynamics in a tiered non-Markovian environment

We introduce a new analytical method for studying the open quantum systems problem of a discrete system weakly coupled to an environment of harmonic oscillators. Our approach is based on a phase space representation of the density matrix for a system coupled to a two-tiered environment. The dynamics of the system and its immediate environment are resolved in a non-Markovian way, and the environmental modes of the inner environment can themselves be damped by a wider `universe'. Applying our approach to the canonical cases of the Rabi and spin-boson models we obtain new analytical expressions for an effective thermalisation temperature and corrections to the environmental response functions as direct consequences of considering such a tiered environment. A comparison with exact numerical simulations confirms that our approximate expressions are remarkably accurate, while their analytic nature offers the prospect of deeper understanding of the physics which they describe. A unique advantage of our method is that it permits the simultaneous inclusion of a continuous bath as well as discrete environmental modes, leading to wide and versatile applicability.


Introduction
The field of open quantum systems, originally devised for quantum optics problems, has recently gained significant traction in the study of condensed matter systems: this is due to the exquisite level of quantum control that is becoming available over increasingly mesoscopic solid state systems, as well as the tantalizing prospect that nature itself may be harnessing quantum effects under adverse 'warm and wet' conditions, e.g. in photosynthesis [1,2] and the avian compass [3,4]. In current literature there is a range of methods to evaluate the evolution of a general open quantum system, from the straightforward but approximate weak-coupling master equation approach [5] through to the fully-numerical path integral based on quasi-adiabatic propagator path integral (QUAPI) [6][7][8][9][10]. It is important to find ways of treating quantum systems embedded in environments that are realistically complex, both in terms of their structure and their non-Markovian nature (i.e. environments which have a 'memory'). When a new approach is analytic rather than numerical, there is the considerable benefit that one gains a route to intuitive insight as well as a simulation tool.
In this paper we introduce a method based on a sequence of three steps: first, we introduce the 'P matrix', which allows a phase space description of a multilevel system coupled to complex environment. Second, we perform a perturbative expansion of the resulting dynamical solution. Finally, we express the reduced dynamics in terms of an influence functional, a quantity which allows new insights into the behaviour of open systems. Our method is intuitive, highly accurate as long as the system environment coupling does not get too large, and works for general spectral densities. In contrast to many conventional open quantum system approaches, such as those mentioned above, we consider a hierarchical environment consisting of two tiers. The outer tier represents a zero-correlation-time heat bath that acts on an inner tier that is the immediate environment of the system. The inner tier may consist of a single harmonic oscillator, a continuous bath of oscillator modes, or any additive combination thereof.
Previous works such as [11][12][13] consider similarly tiered environments for a different conceptual reason: in those cases a single environmental tier is subdivided with the purpose of capturing more accurate, non-Markovian dynamics. In a similar manner, [14] considers a second tier which is constantly randomized for gaining a numerical advantage in simulating a singly tiered environment. By contrast, our approach here is not motivated by 'mathematical' convenience but is rather designed to capture a commonly occurring 'physical' reality. This latter motivation had already been applied to some specific models such as the damped Jaynes-Cummings model [15,16] and fictitious harmonic oscillators [12], and the idea has led to the theory of pseudomodes [17] (intrinsically restricted to zero temperature). A similar idea underlies the so-called 'reaction coordinate' method, where the inner tier is a single harmonic oscillator that is coupled to a wider environment [18][19][20], an approach that is often referred to as a 'structured environment' in the literature [21][22][23][24]. This method employs a mapping between the original environment and a spin-boson model with an effective spectral density [18].
The method we introduce here applies to a general choice of system and bosonic environment at finite temperature, and the two environmental tiers typically represent different environmental influences. There also exist methods for modelling a long or infinite chain of identical environmental tiers, for example, the problem of a quantum system coupled to the end of a linear chain of fermions [25] or bosons [26]. We remark that our method remains applicable when there is no natural division into separate tiers and only a single environment is considered (or when both tiers arise from the same environment). In this case we still obtain non-Markovian contributions to the dynamics, and when applied to canonical cases, we recover known results from the literature. However, our method is more distinctive when two different environmental influences are present.
Another active area of research on open quantum systems is that of hierarchical equations, which was pioneered by Tanimura [27][28][29] in the late 80ʼs. This includes hierarchical equations for both the density matrix [30][31][32][33] and wavefunction [34], generally relying on a specific form of the memory kernel of the bath. Non-Markovian state quantum state diffusion [35][36][37] also makes use of a hierarchy of abstract functionals and has recently been used to study energy transfer in molecular aggregates [38]. Note, however, that the technique presented in this paper is conceptually quite different from any of these hierarchical approaches, since our interest focusses on a doubly tiered physical environment instead of mathematical hierarchies of equations.
Our approach of using a two-tiered environment makes our technique particularly suitable for modelling several of today's most intensely studied experimental systems: this includes many examples of discrete quantum systems interacting with an optical or mechanical resonator, such as, e.g., NV − centres on diamond cantilevers [39,40], quantum dots on carbon nanotubes [41,42], nanomechanical resonators coupled to quantum dots [43] or superconducting qubits [44], and superconducting circuit QED [45,46]. Each of these systems features a high quality resonator, some with extremely high-though of course finite-Q factors, as well as a discrete system whose interaction with the environment will in general not be entirely restricted to the resonator.
Additionally, our technique can be applied to the study of nanoscale energy transfer. For example, the interplay of vibrational modes and the excitonic states in molecular structures are thought to be key to fully understanding photosynthesis [1]. Indeed, a dominant coupling of an energy transfer complex to a small number of discrete vibrational modes may be responsible for efficient energy transfer [47], and previous work has shown how a continuous spectrum of modes can be mapped onto a bath plus one or more coupled and discrete oscillator modes [48,49]. However, new theoretical developments, and further experiments, are needed to understand the functional role of discrete modes in energy transfer systems. The theoretical framework we describe here is ideal for studying this kind of system-discrete mode-bath system and is applicable across a wide range of parameter space. For example, it can accurately reproduce the energy transfer dynamics occurring in the FMO complex [50].
To illustrate our method, we show that it delivers a highly accurate description of the ubiquitous Rabi model, even when the oscillator is damped by a larger environment. As a second example, we take the spin-boson model, showing how our method reduces to the weak-coupling results in the appropriate limit, whilst in general giving better agreement with exact QUAPI calculations than traditional weak-coupling techniques. Moreover, since we do not restrict ourselves to the Markovian limit with a static environment, we are able to explore the case where the bath oscillators are themselves coupled to a larger environment, and we derive analytical expressions for the decoherence and dephasing rates for this case.
This paper is organized as follows: in section 2 we define our model and give a brief introduction to the coherent state representation, and introduce the influence functional. Section 3 introduces the perturbative solution to the case where the environment is a single damped vibrational mode. In section 4 we examine the case of a more complex environment which is defined via a general spectral density, and show that up to second order in perturbation, each mode contributes independently to the dynamics. Section 4.1 studies the spin-boson model, comparing our method to other approaches, and finally, in section 5, we summarize our results and discuss the validity of our technique.

Model
We start with the Hamiltonian where  S is the Hamiltonian of the governing the system of interest. We shall take the 'system Hamiltonian' to be defined on a discrete, finite-dimensional Hilbert space, on which measurements can be performed. No other assumptions are necessary, and in particular  S does not need to be time-independent. The term represents an environment consisting of harmonic oscillators, where a k † (a k ) is the creation (annihilation) operator for a mode with angular frequency ω k . The term = ∑ is the interaction coupling the system (via the system operator V) to the environment.
Equation (1) also includes terms that allow our environment to be coupled to the rest of the universe denoted by  U . When such a wider environment is present, we assume that it is well approximated by an infinite heat bath that is kept in a thermal state. The oscillator modes of the immediate environment are then dynamically driven towards a thermal state by virtue of the environment to universe coupling term  EU . However, unlike conventional Born-Markov weak coupling approaches which commonly keep the entire environment fixed in thermal equilibrium, the inner tier modes will in general deviate from the thermal state. We shall show this adds an exponential cut-off to the response kernel. Figure 1 gives an illustration of our model.
Instead of explicitly treating the coupling between the environment and the rest of the universe with a microscopic derivation, we make the simplifying assumption that  EU is small enough that each mode ω k of the environment simply experiences damping with rate γ k via standard Lindblad operators (for a derivation see, e.g., [5]). For this to be consistent, two conditions must be satisfied: firstly, the damping rate γ ω ≪ k k must be small for each mode, because this is the parameter regime assumed in the derivation of the damped harmonic oscillator master equation. Secondly, the system-environment coupling described by  I may not become too large either, otherwise the damping Lindblad operators acting on each mode are influenced by the presence of the system and our simple independent choice ceases to be a good approximation [51] (also see [15] for a discussion of this approximation in the context of the resonant damped Jaynes-Cummings model).
Finally, we assume that the initial density matrix can be factorized as ρ E , (where  is the appropriate normalization factor).

Coherent representation
To represent the density matrix of a single harmonic oscillator we use the coherent state or P representation [52], which has been extensively studied in quantum optics. The coherent state representation maps between the density matrix of a harmonic oscillator ρ and a function of two continuous variables α α P ( , *) via . The mapping yields the following operator correspondence [52]: † Figure 1. An illustration of the model under study. The system of interest is coupled to an immediate environment, which is in turn coupled to the wider 'universe'. The environment is modelled as a set of harmonic oscillators, whereas the 'universe' weakly dampens each of these oscillators to a thermal state.
For a system with states| 〉 i coupled to an oscillator, instead of a P function we now need a P matrix to represent the density matrix i j i j ,

,
Generalizing from a single mode to a set of modes is straightforward, with the corresponding set of variables A partial trace over the oscillator space is given by i j k For notational ease, from hereon we switch to a vectorized form of the density matrix and operators, mapping n × n matrices A i j , to vectors A i of dimension n 2 . Further, we use the generalized Gell-Mann matrices with the notation from [53]. For an n-site system, these consist of − , the generalized Gell-Mann matrices satisfy: . Any n × n matrix P can be written as a vector P i : (15) n 2 Using this vectorized form we can write the density matrix as , and we are interested in the partial trace over the environment

The influence functional
At this stage, we use the following form for writing down the full dynamics of the reduced system: where U(t) is the propagator (in the vectorized representation) of the system without the environment, and the influence of the rest of the world on the system is encoded in the influence functional Θ t ( ). The motivation for this comes from the Feynman-Vernon influence functional [54] of the same form. Further, we anticipate that 4 For n = 2 (a qubit) ν σ = i i are the Pauli matrices, and for n = 3 we get the Gell-Mann matrices ν λ = i i . this form will be a convenient one for recovering the known exponential decay in the weak-coupling limit. The main result of this paper is that it is possible to find an exact expansion of Θ t ( ) as a perturbation series with respect to the interaction  I , and expansion up to second order recovers the known dephasing and relaxation rates given by standard Born-Markov weak master-equation techniques, but with an added non-Markovian contribution.

A single mode
Let us first examine the case where the environment ω =  a a E † consists of only a single mode. When taking a two-level system (2LS) as the system (a limitation which is not required in the following), then this is just the well-known Rabi model.
In its vectorized form, the system-environment part of Hamiltonian (1) can be decomposed to Then the operator correspondence between ρ and ⃗ P , with the vector is simply the corresponding P representation Fokker-Plank operator [5], i.e. for a single damped oscillator the master equation would read In the vectorized representation, the terms − ×  P i S and gA P g take the place of g n n n , Note that ×  S is Hermitian, and the propagator U(t) satisfies The central strategy of this paper now is to solve equation (23) perturbatively with g being the small parameter, based on the form (18) of the full solution in order to estimate the influence functionalΘ t ( ).

Perturbation series
For the perturbation treatment, we use the expansion hence equation (23) translates to: The solution for the uncoupled system P 0 is simply given by In principle it is possible to solve this series term by term. However, we are interested in the state of the system and not the oscillator, which makes things much easier: we use the boundary condition where α α ⟶ α→∞ P ( ) 0 k n for all k n , . This is justified since the oscillator can be expected not to deviate by too much from a thermal, Gaussian state, and it certainly also should not occupy extreme high-energy states. Therefore performing the integration ∫ ∫ i n n i n n , , , is the matrix equivalent to the superoperator □ V , i.e. at time t = 0 the qubit and the mode are factorized, and the mode is in the thermal state, which gives 1 for all times. The first contribution in the expansion therefore comes from ∫ , which is 2 nd order in the coupling constant g. This is in analogy to the usual QME treatment, where the influence of the environment also enters at the 2 nd order in the coupling constant. In order to solve equation (40) we first need to evaluate 1 , which can be done by invoking the following mathematical procedure: (i) multiply equation (35) by α or α* from the left; (ii) perform the ∫ α integral; (iii) integrate by parts all terms possessing a derivative. The sequence of these steps yields the following two equations: which after a bit of algebra and ODE solving yield a solution for ∫ α P 1 . Substituting this solution into equation (40) then results in Here, the notation × • V Ṽ ,˜denotes operators in the Heisenberg picture, n n n 2 2 2 At this point we note that the influence functionalΘ t ( ) up to second-order in g is then given by equation (47) and We proceed by showing that this provides a highly accurate solution for the single mode case in the weakcoupling limit. We shall then generalize the technique to an environment consisting of a (quasi)continuous bath of oscillators. In appendix B we sketch the derivation of higher-order terms in the perturbation series.

Example: the (damped) Rabi model
The Rabi model, consisting of a coupled 2LS to a harmonic oscillator, represents perhaps the most basic and ubiquitous compound quantum system. Focussing only on the dynamics of the 2LS and tracing over the oscillator then results in arguably the conceptually most simple and yet a highly non-trivial open systems problem. Let us consider the Rabi Hamiltonian where σ i are the usual Pauli matrices referring to the 2LS. In this case, we immediately find that the matrices ). In this basis, the real terms on the diagonal of Θ t ( ) that are proportional to t and correspond to the finite eigenvalues, are both equal to the dephasing rate. The one corresponding to the vanishing eigenvalue is the relaxation rate. These rates are given by where Ω ϵ ω = + 2 2 is the Rabi frequency. Note that in the limit γ → 0, i.e. no damping on the oscillator from the wider environment or universe, we recover the standard Born-Markov ME result for relaxation and dephasing, given in equations (C11)-(C12). The imaginary parts on the diagonal of Θ t ( ) correspond to the Lamb shift Hamiltonian, given by where σ z is given by writing the system Hamiltonian, i.e. the first two terms in equation (54) in its diagonal basis z S Again, in the limit γ → 0 we recover the 'standard' Lamb shift given in equation (C7). Furthermore, we can extract the steady state of the system at long times: at times much larger than the relaxation time, the system tends to the state This is indeed only the expected thermal system state when γ → 0 and ω Ω → , i.e. no damping and when oscillator and system are resonant. However, one should take this limit with caution, because for vanishing damping, γ → 0 the relaxation time Γ − relax 1 tends to infinity and the system will thus never actually reach this state. In figure 2 we plot the effective temperature, that is, the temperatureT eff given by equating − k T exp [˜] s b eff with equation (62). On the same figure we plot the relaxation rate for the same parameters, showing a Lorentzian peak in efficiency near resonance.
We note that in general the effective temperature differs from the temperature of the universe. In order to explain this apparent discrepancy, we examine equation (62): the universe is only directly coupled to the oscillator which has energy levels spacing of ω, this accounts for the term is maximized when on resonance (ω Ω = ). Detuning suggests that in order to extract energy from the qubit, the universe exchanges energy with the oscillator to match the detuning. This adds uncertainty to the system effectively increasing the temperature. The system-environment coupling γ adds additional uncertainty.
We also note that in this scheme we do not keep track of the environment, only trace over it. The thermal state of system + environment is proportional to S E I , i.e. the system and environment are entangled, and defining a temperature of just one subsystem is questionable.
The example we discuss in this section is formally equivalent to the reaction coordinate [18][19][20] or structured environment [21][22][23][24] model in the weak coupling and weak damping regime. Here, the reaction coordinate model employs an effective spectral density with a Lorenzian peak, yielding the same rates as equations (58)-(60) except for the 'counter rotating' terms Ω ω ∼ + − ( ) n (which are typically small). Interestingly however, this nice agreement only extends to the real part of the response function, D(t), which determines the damping rates. By contrast, the modified spectral density of the reaction coordinate method does not account for corrections to the imaginary part D t ( ) 1 , which yields the long time asymptotic behaviour of the system. To ensure that our approach does indeed deliver the correct steady state, we have made a comparison with an exact numerical simulation of the dynamics given by Hamiltonian (54) (with +   EU U replaced by a Lindblad dissipator). We obtain perfect agreement between equation (62) and a purely numerical simulation in the weak coupling regime.
In figure 3 we plot a comparison between equation (18) withΘ t ( ) approximated by equation (53), and exact numerical simulation, showing that for the weak-coupling regime there is a very good agreement between the two.

Extending the analysis to a multimode environment
In the previous section the 'environment' consisted of only one single harmonic oscillator. However, adding multiple oscillators is straightforward, and in the weak coupling limit, where environmental influence is assumed to be small, each environmental mode contributes to the influence functional Θ t ( ) independently. The difference is that now the environment Hamiltonian  E has a set of modes, and in our vectorized form the equivalent of equations (19)-(22) becomes The derivation for this case is very similar to the single mode case and is given in full detail in appendix A. Once more, the influence of the bath on the system's dynamics is given by equation (18), where now  Here Ã 1,2 are given by equation (48), and we adapt our notation to match that common in the literature on phonon baths, introducing the (damped) phonon response function defined as are the (damped) dissipation and response kernels, respectively. In terms of the spectral density function . We note that for the case of γ ω = ( ) 0, i.e. when there is no external universe, the result (68) is exactly coincides with the well-studied time-convolutionless projection operator technique (TCL) from the literature when the TCL generator is expanded to second order in the system-environment coupling, see [11].
It is interesting to note that the thermalization of the immediate environment by the wider universe is fully captured by switching to the above generalized form of the response kernel (69) (within a perturbative treatment to second order, higher orders give additional corrections, see appendix B). At T = 0 our expression is in full agreement with the previously derived zero temperature response function of the damped spin-boson model given in [12]. We suggest that the same kernel redefinition might also be applicable to other methods of studying open quantum systems, giving a simple recipe to adding a wider universe on top of a standard open system.

Example: the spin-boson model
To apply our generalized multimode technique to a particular example, we look at the well studied case of the (biased) spin-boson model with the following Hamiltonian: In this case, just like for the Rabi model, the system is two-dimensional and its P vector has four components (σ σ σ  , , , x y z ), and × × •  V V , , S are again given by equations (55)-(57). Since we have already calculated the relaxation and dephasing rates for the single mode case, showing that the different modes contribute independently for Θ t ( ) in the weak-coupling regime, we can immediately write down the following expressions for the relaxation rates: we only need to add a summation ∑ k over the different modes to equations (58)-(59): We note that, as discussed at the end of section 4, in the limit of γ → 0 k , we recover the known weak-coupling rates, see [55] or appendix C. The second part of equation (74) is known as the pure dephasing constant.
Below we study the no-bias case, setting ϵ = 0: the system Hamiltonian ( ×  S in our language) is static, hence the propagator U is given by To calculate Θ t ( ), we can make a change of variables in the to get the expression: In the above expression, Θ relax induces the relaxation and decoherence, Θ LS induces the Lamb-shift, and Θ th steers the system towards the thermal state. Θ RW is usually ignored under the rotating wave approximation. If one is interested in times τ ≫ t b much longer than the memory of the bath , it is justified to let the upper limit of the integrals go to infinity. For this case it is most insightful to examine this result in light of the standard quantum-optical master equation approach: in the standard approach, remarkably one gets exactly the same expressions as the above equation (75) [without equation (79)], but with an interesting change: The terms which are not proportional to t capture non-Markovian contributions, giving information about the bath's reorganization time. Interestingly, each of the environmental effects possesses its own timescale, and these are estimated by It is noteworthy that the reorganization times can be negative. This could happen when, for example, initially for τ ≲ t b the dephasing process, which includes a non-Markovian component, is more aggressive than at later times when it assumes a stable value. Then, as the aggressive decay stops, the population of the system has fallen by a greater amount than it would have done under the stable, long lived decay process. Thus the system appears as if it has been evolving under the stable dephasing rate for a longer time than it actually has, and hence the negative reorganization time. We note that the terms (81)-(83) in the limit γ → 0 are known in the literature as those leading to the slippage of initial conditions, and are important for preserving the positivity of the reduced density matrix [56,57].
The steady-state of the system is given by x relax 1 0 1 0 A comparison between the standard Markovian Master equation, the current method and exact numerical simulation for the case of a super-Ohmic environment is shown in figure 4. The QUAPI technique [6][7][8] is used as an exact numerical benchmark curve: our calculation uses nine kernel time steps, covering a total kernel memory time of 2 ps and is fully converged. The standard Born-Markov weak-coupling approach is given in appendix C. Clearly, our method's non-Markovian nature and lack of Born approximation results in an impressive improvement over the standard Born-Markov weak coupling ME approach. For this particular comparison, since there is no wider universe involved, γ ω = ( ) 0, the current method is equivalent to the second-order TCL approach, which also does not employ any approximations beyond a perturbation in the system-environment coupling. However, a key strength of the current formulation is that it is trivial to include a wider universe, which simply enters in the form of an exponential cut-off to the response function.
We note that this method allows us to easily study the case where the spectral density has several discrete sharp peaks as well as a smooth background, which is believed to be the case in many (if not all) systems studied in quantum biology [10,58]. In this case the response function vanishes very slowly, which makes an exact numerical treatment extremely demanding, as a long history of the system needs to be tracked. In some papers, such as [58] this issue is resolved by approximating a delta-function peak in the spectral density as a Lorentzian with a finite width. We note that if one allows this single peak to be damped, then in light of equation (62), this mode drives the system to an effective temperature different from the initial temperature of the environment T. Hence replacing discrete modes with Lorentzian distributions added to a continuous spectral density may in some parameters regimes become a questionable approximation. By contrast, the additive property of modes to the influence functionalΘ t ( ) here allows us to combine a discrete set of modes with a smooth background by taking smooth discrete As an example for this, let us study the spin-boson model with a smooth background of oscillators plus a more strongly coupled discrete peak of frequency ω s in the environment. We single out this peak and label it henceforth with a subscript s, writing the system-environment Hamiltonian as a a a a g a a g a a 1 2 In figure 5 we start with the system in its ground state and plot the excited state population ρ xx as a function of time, for the cases where the system is only coupled to a smooth environment ( → g 0 s ), only coupled to a single mode ( → g { } 0 k ), and for the combined case. Due to the non-Markovian nature of this method, we are able to capture the revival effect [59] for the Rabi model. These revivals can be damped via a combination of two mechanisms: either the mode itself is coupled to a wider environment damping it, or there might be an additional continuous bath directly damping the system. In figure 6 we plot the first case, where the environment consists of a single damped mode. The damping of the mode induces relaxation rate given by Γ= 1 equation (58). We also plot the decay envelope for this case, as well as the decay envelope produced by coupling of the system to a continuous bath and no  damping on the mode, choosing parameters such that the relaxation rate induced by the bath equation (74)  We note that the second case yields an exponential envelope to the dynamics for times ≫ t t R relax , while for a single damped mode with damping rate γ, the envelope only becomes exponential for times γ ≫ t 1 , which could be much longer. We note that the Lamb-shift given by equation (77) also differs between the two cases, albeit in the plotted parameter regime this difference is very subtle and not shown.

Discussion and conclusion
We have introduced a novel method for studying a ubiquitous open quantum systems problem. Our approach differentiates between the immediate environment of the system of interest and a wider universe which effectively serves as a heat bath for this environment; this hierarchy of environments corresponds to many practical situations and is-remarkably-accomplished by a simple redefinition of the response kernel. The expressions resulting from our method are easy to evaluate numerically, and scale favourably with increasing system size. Moreover, the method still leads to soluble equations when the system of interest possesses a general time dependent Hamiltonian.
Whilst our method is limited to the weak coupling regime, it performs favourably when compared with traditional Born-Markov weak coupling master equations. Its approximate analytical expressions scale well with increasing system size and permit valuable physical insight, in contrast to some numerically exact approaches. Like many recent developments in the field of open systems, see e.g. [38,49,60,61] (and with the notable exception of [26]), we do not presently have stringent criteria demarcating its precise regime of validity, which must thus be established by comparison with exact numerics. As a general guideline, however, our technique can be expected to perform well whenever other weak coupling approaches such as the timeconvolutionless or the Nakajima-Zwanzig projection operator expansions [5] are valid for the system-toimmediate-environment coupling. As an additional criterion, our treatment of the wider universe (if present) assumes that γ ω ≪ k k , i.e. each mode is weakly coupled to its heat bath. We have benchmarked our technique against the well-studied spin-boson model and the Rabi model, finding it leads to expressions that are indeed highly accurate when compared with numerically converged solutions. This remains true even for coupling strengths where a conventional standard second order Born-Markov master equation begins to performs poorly, and exactly recovers the time-convolutionless solution when no wider universe is present. For cases when the system-environment coupling is not sufficiently weak for the second order expansion of the interaction, we provide an explicit recipe to calculate higher orders in the perturbation series. Perhaps a unique advantage of this approach is that these two models, i.e. the Rabi and the spin-boson models can easily be combined even for long-time dynamics. This makes our method eminently suitable for studying the exciton energy transfer in photosynthetic or artificial molecular systems, since the coupling of the excitonic degree of freedom to both the vibrational quasi-continuum of the wider protein scaffolding as well as to specific localized vibronic modes is believed to be of crucial functional importance. Figure 6. Long time population of a TLS (blue), illustrating the revivals which occur when a discrete system is coupled to a single (damped) oscillator mode. The corresponding relaxation envelope (purple) and that of an undamped mode but where the system is coupled to a bath (yellow) are also shown. Here, we have chosen a bath coupling strength to obtain the same average relaxation rate for both cases (see inset), even though this does not become apparent during the first two revivals. The parameters in this figure are chosen to show revivals, so that the mode is almost resonant with the TLS and the damping is weak, Δ π = and γ γ ω = ( ) k k is the damping rate of mode k. The matrices A g (k) are given by  The solution for the uncoupled system P 0 is then equivalent to the single mode case, and is given by (assuming a factorized initial state): We assume that α α ⟶ α→∞ P ( ) 0 k l n for all k n l , , for the same reasons given in the main text. Performing the integration ∫ α on equations (A11)-(A13) yields where just as before, × V is the equivalent of □ V [ , ]and is given by equation (42), and the initial condition is and a corresponding equation for α k * . Crucially, there is no sum over k here, which means each k gives rise to exactly two equations of the type of equations (45) and (46), which we have already solved. The first nonvanishing term is hence given by which is just equation (47) with an added sum over all modes, and where × • Ṽ , are given by equation (48). From here we continue to equation (68).

Appendix B. Higher orders calculation
In this Appendix, we show how to calculate higher orders of the influence functional Θ t ( ) defined in equation (18), where the main text only gives the 2 nd order expression. We also show that, in analogy to the known result of the non-hierarchichal case [5], all of the odd orders vanish when the initial state factorises ρ ρ ρ = ⨂ (0) (0) s E th . We start by giving a formal expression of the quantity  In the above expression we usedV t ( ) which is defined by equation (48), and the sloppy notation