Using the Environment to Understand non-Markovian Open Quantum Systems

Tracing out the environmental degrees of freedom is a necessary procedure when simulating open quantum systems. While being an essential step in deriving a tractable master equation it represents a loss of information. In situations where there is strong interplay between the system and environmental degrees of freedom this loss makes understanding the dynamics challenging. These dynamics, when viewed in isolation, have no time-local description: they are non-Markovian and memory effects induce complex features that are difficult to interpret. To address this problem, we here show how to use system correlations, calculated by any method, to infer any correlation function of a Gaussian environment, so long as the coupling between system and environment is linear. This not only allows reconstruction of the full dynamics of both system and environment, but also opens avenues into studying the effect of a system on its environment. In order to obtain accurate bath dynamics, we exploit a numerically exact approach to simulating the system dynamics, which is based on the construction and contraction of a tensor network that represents the process tensor of this open quantum system. Using this we are able to find any system correlation function exactly. To demonstrate the applicability of our method we show how heat moves between different modes of a bosonic bath when coupled to a two-level system that is subject to an off-resonant drive.


Introduction
When modelling a realistic quantum system the effect of its environment must be captured in some form [1]. The environment consists of many degrees of freedom and is often approximated as an infinite collection of harmonic oscillators. Thus, an explicit description of the environment is usually impossible. Most methods for describing open quantum systems rely on tracing out the environment to obtain an effective description of the system involving only a tractable number of degrees of freedom. When the systemenvironment coupling is strong, and/or when the spectral density of the environment is structured, the full dynamics involves a significant interplay between the system and environment degrees of freedom [2][3][4][5][6][7][8][9][10]. When the environment is traced out, the dynamics of the system cannot be described by a time-local master equation and in this sense are non-Markovian [11,12], whose simulation generally requires sophisticated numerical methods. In addition, such system dynamics are often complex to interpret.
Analyzing the behaviour of the environment is useful for providing insights into processes where the system-environment interaction is non-trivial or for any application where the goal is to manipulate the environment via the system. An important example of this class of problem occurs in modelling quantum thermal machines [13,14] where distinct thermodynamic effects beyond weak coupling can be observed [15,16]. Access to the state of the environment could also be used to explicitly track the formation of bound states such as polarons [17][18][19][20] or polaritons [17,21]. There is also potential for simulating phase transitions in the environment, such as to lasing or superradiant states [22] in multimode cavity QED [23]. As well as these, the technique we introduce here allows for calculation of out-of-timeorder correlations of the environment; these are known to characterise the degree of information scrambling in many-body quantum systems [24][25][26].
Some existing techniques allow varying levels of insight into the behaviour of the environment. These generally rely on a direct calculation of the dynamics of the environment itself. For example, the reaction coordinate mapping explicitly calculates the time-dependence of a collective coordinate that gives a sense of the behaviour of the environment [27,28]. It has also been shown that the auxiliary density operators calculated in the hierarchy equations of motion method can be used to calculate the same collective coordinate dynamics [29] as well as higher order moments of the total bath coupling operator [30][31][32]. Counting field techniques can also be used to calculate moments of the total bath coupling operator [33][34][35]. Techniques such as TEDOPA [36] map the star topology of an environment to a chain of modes with nearest neighbour couplings. Dynamics within this chain can also be used to characterise the behaviour of the environment [37] and these dynamics can be passed through an inverse chain-mapping to extract the full environmental dynamics [38,39].
In contrast here we give a prescription that allows calculation of any multi-time bath correlation function using only system moments of equal order: no auxiliary modes or additional degrees of freedom are required. Our method can therefore be used in conjunction with any technique that allows for the accurate calculation of system correlation functions, in order to yield bath dynamics and correlations. This is even possible using those methods where the entire environment is traced out.
Accurately calculating the required correlation functions for a system undergoing non-Markovian evolution is numerically challenging for all but the simplest cases. It is necessary to use a technique which is capable of recording the history of the system dynamics while allowing for the application of arbitrary sequences of system operators. Here we employ a version of the timeevolving matrix product operator (TEMPO) algorithm [40]. This has been shown to efficiently capture the required time non-local influences by representing the discretised trajectories of the system as a matrix product state (MPS). This is very similar to the use of MPS in simulating 1D quantum many-body problems which can be compressed by truncating the sizes of the matrices, keeping only the most relevant contributions. An imaginary time version of TEMPO can be used to determine the thermal state of a system with a strongly-coupled environment [41]. However, even with the TEMPO algorithm, calculating all of the necessary correlation functions would require many separate numerical simulations. This difficulty can be overcome by building on the approach developed in [42] where a modification of the network contracted in the TEMPO algorithm allowed efficient generation of an object called the process tensor. From a single process tensor it is possible to access not just the system dynamics but also all system correlation functions at arbitrary times [43]. This significantly reduces the computational overhead of calculating the behaviour of the environment.
In a previous work [44] we showed how the displacement of a bosonic bath that is linearly coupled to a spin system can be calculated from a simple transformation of the system dynamics. We used this to gain insight into the origin of some of the non-Markovian features in the dynamics. The approach used there is limited to the calculation of first moments. In Sec. 2 of this paper we reformulate the problem in terms of a moment generating functional, to show how we can calculate all higher order bath correlation functions from system correlations alone. An overview of the process tensor formulation [43,45] of open quantum systems is given in Sec. 3 where we demonstrate how it can be used to efficiently calculate the system moments necessary for calculation of higher order bath correlation functions. Finally, in Sec. 4 we use this approach to study how the biased spin-boson model [46] can be tuned to redistribute energy in its environment.

Bath Dynamics
Here we provide a method to calculate any correlation function of bath operators in terms of those of system operators, for a system linearly coupled to a Gaussian environment. We do this by expressing the relevant expectation values of environment operators as derivatives of a generating functional. By adding source terms to the propagator that remain easy to integrate out for the Gaussian bath, the derivatives required to calculate the relevant moments can then be carried out explicitly.
The system-environment Hamiltonian we consider has the general form where a q destroys an excitation in bosonic environment mode q with frequency ω q . Here, s is a general operator acting exclusively on the system Hilbert space and g q is the coupling strength of mode q to the system. Formally, the full dynamics, in the interaction picture, can be written as [47]: The full density operator of the system and environment, is ρ(t) and To proceed we assume that the initial state is factorizable and can be written as ρ(0) = ρ S (0) ⊗ ρ B (0). We can then trace out the bath and arrive at an expression for the evolution of the system: where · B denotes an expectation value with respect to the initial state of the bath.
If we further assume that the initial state of the bath is Gaussian we can apply the relation e X = e X 2 /2+ X which is valid for expectations taken over any X that has a Gaussian distribution. This condition holds for L I , which is linear in the a ( †) q , when the expectation is taken with respect to ρ B . We have defined the Hamiltonians such that L I = 0 and so the environment part of Eq. (3) can be simplified to: The right-hand side of the above equation is exactly the Feynman-Vernon influence functional [48] in superoperator form. This relation holds for any L I (t) that is linear in a q and a † q , and so we can add extra such terms that will allow us to use a moment generating function approach to calculate expectation values of bath operators, as we will now explain. We include the extra terms by making the replacement L I (t) → L I (Λ, Λ * , t) with: where Λ α( * ) q (t) are scalar valued functions associated with superoperators A α( †) q (t). These superoperators can be expressed in terms of a ( †) acting either to the left or right of ρ: Using Eq. (5) we can define a moment generating functional where the angular brackets now denote the expectation be taken with respect to the total initial density matrix ρ. We have included a scalar valued function φ in the definition of G; the exact form of φ is dependent on the order of the ladder operators in the quantity of interest. This arises due to Baker-Campbell-Hausdorff factors that appear in the construction of the generating functional, see App. B for more details. As an example any normal-ordered expression results in The expectation of any function of left-acting ladder operators f and right-acting ladder operators g can then be expressed as functional derivatives of Eq. (7): In Eq. (9) we have implicitly performed the derivative before taking either of the traces over the system or bath. These operations commute so we can either: trace over both system and environment and then take a numerical derivative (this is the counting field approach [33,35,49]); or take an analytic derivative after tracing out the bath to get an expression purely in terms of system correlation functions. We use the latter approach as it allows us to calculate multiple different bath observables from a single influence functional.
As an example, we can use the technique introduced here to derive the following expression for the occupation of a specific mode (see App. B for details): where n q (t) = a † q (t)a q (t) and T is the temperature of the bath. Expressions such as Eq. (10) can be obtained for any bath correlation function, using the general technique outlined above. We give a further example of how this approach can be used to calculate multi-time correlation functions in App. C. These expressions will involve integrals over system correlation functions that are of the same order or lower. Any technique can be used to calculate the required system correlation functions, though naturally any approximations made in obtaining these will be reflected in the inferred bath dynamics. We include details of how n q (t) can be calculated for a simple model including just two modes in the environment in App. D. Further, in App. E we show how these results can be generalised to find higher order statistics of environment modes such as g (2) .
To calculate a single correlation function Tr[s 1 (t )s 2 (t )ρ] in the integrand of Eq. (10) would typically require evolving ρ to time t , acting with s 2 , evolving to t and finally taking the expectation of s 1 . To calculate the time evolution for each required t and t needed to evaluate the integral in Eq. (10) is a computationally expensive procedure. To address this problem we use the process tensor approach as this requires a single (albeit individually more expensive) calculation from which any system observable can be determined with little further cost. Below, we introduce the idea of a process tensor conceptually and give an overview of how it can be calculated efficiently by contracting a tensor network.
To summarise, in this section we have described a method for finding bath correlation functions from system correlations alone. In the next section we will introduce process tensors, which are ideal for finding system correlation functions efficiently.

Process Tensors
To be able to show the full utility of the approach derived above we need to make use of a concrete method which is able to accurately provide us with the required correlation functions of system observables. To that end we employ the processes tensor formalism which we outline below. We emphasise at this point that any method capable of calculating these correlation functions can be used.
In this section we review the process tensor formalism, which is a general operational approach to non-Markovian open quantum systems [45]. It supposes a finite set of control operations act on the system at a sequence of time points, such as the preparation of a particular state or a projective measurement. The process tensor has two indices for each time at which a control operation can be applied. It therefore allows the computation of any multi-time correlation function of system operators by inserting control operations at the corresponding times.
For Gaussian environments it has been shown that the process tensor can be efficiently calculated by contracting a tensor network [42,50] (see also Fig. 1 which we will discuss in detail shortly). This is done by an alternative contraction order of the TEMPO network originally derived in Ref.
[40] as a representation of the quasi-adiabatic path integral (QUAPI) method [51]. The process tensor fully captures the non-Markovian nature of the interaction of the system with its environment and is therefore challenging to compute. Once obtained, however, it can be used repeatedly to compute any multitime correlation function at very little added numerical cost. With this, we have a tool at hand to evaluate the occupation, Eq. (10), with moderate computational effort. In the following, and with reference to Fig. 1, we outline the form of the TEMPO tensor network and then explain how it can be adapted to compute the process tensor. We refer the reader to [40,52] for more details on the derivation of the network and how it can be efficiently contracted.
The TEMPO network can be understood in terms of superoperators beginning from a discretisation of Eq. (3) with Eq. (4) substituted: We move back into the Schrödinger picture and perform a symmetrised second-order Trotter splitting between the system and interaction parts: (12) where we have discretised the evolution into N timesteps of length ∆t and t k = k∆t. We have also defined the system density matrix in the Schrödinger picture as ρ (S) S (t). Here we assume H S is time independent such that U = e L S ∆t/2 , with the half-timesteps arising due to the Trotter splitting we use. For each timestep t k the terms k k =0 I k−k (t k , t k ) account for the influence of the past system states on t k via the environment: We can represent this propagation as a tensor network. For example, the leftmost network in Fig. 1(a) corresponds to propagating a state forward in time by 4 steps. In the superoperator description an initial density operator in a Ddimensional Hilbert space is represented as a D 2 component vector (purple circle in Fig 1). This can be regarded as a one-site matrix product state (MPS). Within this description the system propagators U are rank-2 tensors (matrices) and the influence functions I are rank-4 tensors.
At each timestep the current state is contracted with a matrix product operator (MPO), corresponding to a row in the leftmost network Here the I tensors account for the non-Markovian influence of the environment and the U tensors the time-local system Hamiltonian evolution. Moving through the steps gives how the tensor network representation of the process tensor relates to that used for propagation. First the legs between each timestep are cut and left open and then the network is contracted column-by-column with standard tensor network routines. (b) Calculation of the correlation function B(t 3 )A(t 1 ) using the process tensor shown in the final step of (a). For a product Hilbert space spanned by {|e i ⊗ |e j } the trace cap is defined as i |e i ⊗ |e i , the vectorised Hilbert space identity matrix and a null eigenvector of any Liouvillian.
in Fig. 1(a), that grows the MPS by one site. In other literature the growing MPS is referred to as the augmented density tensor (ADT) and contains information about both the current system state and its correlations with the environment. The state can be extracted at each timestep by summing over all but the rightmost legs of the ADT. For example, the fully contracted network at the left in Fig. 1(a) would be the state after four timesteps. We show how this network can be converted into a process tensor following the steps in Fig. 1(a): The process tensor can be identified as the full propagation network with the U propagators between each timestep disconnected and the legs left open. This network can then be contracted from left to right and the process tensor stored as an MPO.
Alternatively, one can view the system propagators U as control operations and exclude them from the process tensor. Such a division can be exploited to efficiently explore the dynamics for different system Hamiltonians [43]. In this work, however, we include the system propagators in the process tensor, and we are then able to insert arbitrary operations between the open legs (for example the operations in Fig.1(b) shown as orange and blue circles). This process allows us to efficiently calculate multitime correlation functions. Alternatively, we can retrieve the system evolution by inserting identity operations.
Carrying out the contraction step in Fig. 1(a), column-by-column, would result in an object that grows exponentially with each column contracted. Many of the degrees of freedom of this object, however, contribute negligibly and can be discarded without significantly sacrificing accuracy [53,54]. These irrelevant degrees of freedom can be systematically identified by performing a singular value decomposition (SVD) on each tensor in the chain and dispensing with any components that have singular values below a threshold value. It is difficult to make any concrete statements on what the discarded components correspond to. Regardless of this we can be confident that we have retained enough singular values by ensuring the results are converged with increasing the SVD threshold value. These SVD sweeps are carried out once in each direction after contracting a column. For all of the results in this paper we discarded components with singular values of relative magnitude less than 10 −8 . The other main approximation we make is a finite memory cutoff where the influence of the bath is only stored for K timesteps into the past [51,55]. The value of K depends on ∆t: for convergence we require that K be increased until K∆t exceeds the correlation time of the bath. However ∆t must in turn be reduced to min-imise the error from the Trotter splitting. For all the results presented below in Sec. 4 the parameters Ω∆t = 0.05 and K = 50 were sufficient for convergence. These parameter values lead to more intensive calculations than would typically required for this form of coupling. This is because the calculation of bath dynamics is highly sensitive to small errors in system expectation values. Since these are integrated over, minor errors which propagate in time can lead to large errors in the resulting bath observables.
The final process tensor is of the form shown in the final step of Fig. 1(a) which for N time points is N − 1 sites long; the bond dimension is related to the number of relevant bath degrees of freedom [56].
From this a correlation such as B(t 3 )A(t 1 ) , where A is an arbitrary system superoperator, can be calculated by contracting the open legs of the process tensor with A and B inserted at the corresponding timesteps as shown in Fig. 1(b). This process is straightforward to generalize to higher order correlation functions, through the contraction of further system superoperators.
To summarize this section: we have shown how the process tensor is ideally suited to finding system correlation functions of an open system that is strongly coupled to its environment. In particular, the process tensor needs only to be calculated once, before contracting with relevant system operators enable the calculation of any desired correlation function (up to the final time chosen for the process tensor). Multitime bath correlation functions can then be found from these system correlation functions by using the method described in Sec. 2. In the next section we will give an example of this procedure for a well known model.

Biased Spin Boson Model
In this section we study the simple example system of a driven spin-1/2 particle coupled to a continuum of bosonic modes, the biased spin-boson model [46,57]. Studies of the spin-boson model are ubiquitous in the field of open quantum systems for two reasons. First, the Hamiltonian is relatively simple while still allowing non-trivial behaviour, making the spin-boson model inviting for those looking to test new methods. Second, it provides the paradigmatic model for dissipative two-level systems meaning it can be applied to a variety of physical processes for example quantum dot emission in the presence of phonons [52,58], energy transfer in biological or molecular systems [5,28,44,59], the interaction of superconducting qubits with microwave waveguides [60] quantum dots interacting with micromechanical resonators [61] and energy transport in solid state systems [62]. The Hamiltonian of the system is: where we have introduced the spin-1/2 operators s z = (|1 1| − |0 0|)/2 and s x = (|1 0| + |0 1|)/2 which act on the states {|0 , |1 }. The transition between these states is driven classically with strength Ω and bias . The system is in turn coupled with strength g q to a bath of bosonic modes of energy ω q . The bosonic degrees of freedom are described by ladder operators a q and a † q which satisfy the canonical commutation relations: [a q , a † q ] = δ qq and [a q , a q ] = [a † q , a † q ] = 0. The environment in this model is characterised by its spectral density which captures both the coupling to each mode g q and the density of states of the environment. It is defined as: This can be related to the auto-correlation function of the bath coupling operator B = q g q (a q + a † q ), where we have introduced and T is the temperature of the bath. We will also find it useful to define the derivative of this quantityC(ω, t) = ∂C(ω, t)/∂t.
The spectral density is typically given a functional form informed by experimental measurements. Here we use an Ohmic spectral density with an exponential cutoff [46]: where ω c is the cut-off frequency and α a dimensionless coupling constant.

Steady State Properties
To analyze the behaviour of the system by looking at both the system and environment degrees of freedom we wish to look at how energy is redistributed during the dynamics induced by the spin boson model. To this end we consider the change in energy of each bath mode, defined as: This quantity gives the difference in energy of a particular mode q compared to its initial, thermal equilibrium, value. However, this is ill-defined in the continuum limit. We cannot refer to the energy change of a specific mode of the environment as the coupling to any single mode is infinitesimal. Instead we consider the change over a narrow band of modes: The frequency interval is of width δ centered at ω. We can then use the machinery described above to compute this quantity in the steady state of the spin boson model. For the results in this section we choose a relatively weak system-bath coupling. This is not due to any inherent limitations of the methods used but rather because the behaviour in this limit is more interesting. At stronger coupling contributions from the re-organisation energy of the environment dominate and the behavior reduces to that of the exactly solvable independent boson model. This effect has been seen in previous calculations involving the total heat of the bath [33].
In the main panel of Fig. 2 we show the periodaveraged steady state value for ∆Q as a function of the bias with the two-level system initialised in the low energy state, |0 . This initial state is uncorrelated with the environment and so as the system moves towards thermal equilibrium energy needs to flow between the system and environment: it is this that we track. The periodaveraging is performed over a single oscillation of each mode, such that we calculate where t SS is chosen to be long enough that the system density operator is stationarẏ ρ S (t ≥ t SS ) 0. We do this because even after the system degrees of freedom reach a stationary state there are still oscillations in the mode occupations which need to be averaged over to find the appropriate quantity. This change in energy is then able to probe the way in which energy is removed from or added to particular environment modes to populate the spin-1/2 in its stationary configuration. We see that for small bias the initial system energy is higher than that in thermal equilibrium and hence the environment gains energy from the system, while at large bias the opposite is true. This might seem at odds with our later analysis which shows that the total bath energy increases at all values of . This is a result of the energy gained from the interaction being switched on at t = 0, known as the re-organisation energy. The contribution from the re-organisation energy appears in the main panel of Fig. 2 as a pale red background heating. It is only through the mode resolution of our approach that we see that it is possible for the direction of energy flow into/out of the bath to be frequency-dependent. In the following we will expand on the intuition for these observed effects.
To confirm the accuracy of this approach we also checked that the total energy was conserved. This required calculation of the total change in H S , H I and H B which we will denote ∆Q S , ∆Q I and ∆Q B respectively. The change in system energy, ∆Q S is trivial to calculate from the dynamics. The total energy exchanged with the bath, ∆Q B , is given by

×C(t − t ) dt dt . (22)
The bath will generally have a finite memory time t K such that the auto-correlation function satisfies C(t > t K ) 0 and hence the derivative also vanishes. Therefore if the total energy exchanged is the relevant quantity then only correlation functions that span this memory time are required.
A similar expression can be derived for ∆Q I again using the generating functional formalism. This is given by: All of the bath observables calculated here can be derived from the same set of correlation functions. The only variation is in the form of the kernel used for the integral transformation. Generally the dynamics of an n th order bath correlation function up to time t depends on the set of n th order system correlation functions up to that point.
In the upper panel of Fig. 2 we compare the numerical calculation of the heat gained by the bath, ∆Q B with the remaining contribution −(∆Q I +∆Q S ). We see that, up to small numerical inaccuracies, these quantities are effectively equal as would be expected from energy conservation.
To gain a better understanding of the results seen here we proceed to develop semi-analytic expressions for the changes in energy observed. In Ref. [33] it was shown that for weak coupling ∆Q B can be well approximated by where E r is the reorganisation energy of the bath given by and ∆Q G S is approximate change in system energy defined as where ρ SS (T ) is the reduced Gibbs state of the system Hamiltonian, Since the coupling to the environment is weak we expect there to be no significant correlations between the system and environment such that the system simply reaches a state in approximate thermal equilibrium with respect to H S at the temperature of the environment. From this we can then find an analytic expression for ∆Q G S : whereΩ = √ Ω 2 + 2 is the generalised Rabi frequency and we have used ρ S (0) = |0 0|. This in turn allows us to derive an (approximate) analytic expression for ∆Q B by substituting (28) and (25) into (24) giving In the upper panel of Fig. 2 we compare the analytic approximation given by Eq. (29) with the numerical results and find good agreement as would be expected at weak system environment coupling. From the form of Eq (28) we can identify a crossover temperature for the environment, T * , above which the dynamics induced by the coupling to the environment increases the system energy and vice-versa. This is given by the solution of For initial states that are incoherent mixtures of the two basis states ρ S (0) = p |0 0| + (1 − p) |1 1| we can solve this to find In Fig. 3 we show the steady state heat change in the environment ∆Q(ω) as a function of temperature for three values of p. In each case we see that for temperatures below T * the bath absorbs energy from the system and at temperatures above T * the opposite occurs. The main transfer occurs at a frequency close to the Rabi frequency of the system. However, as temperature increases the heat transfer occurs across a broader distribution of modes. This is linked to the broadening of the Bose-Einstein distribution meaning that the system can absorb energy from a wider range of frequencies. From this we can see that the analytic expression is a good predictor of whether the system has absorbed or emitted energy into the bath at a specific frequency.

Dynamics
We now move on to show how the methods described above can also be used to access the complete dynamics of the bath. To do this we implement a time-dependent driving sequence with the goal of demonstrating the absorption of energy from the bath over a range of frequencies followed by emission into the bath over a different frequency range. This shows how, by controlling the dynamics of the system, we can control the state of the environment. We initialize the system in the |0 state at t = 0 and allow the bias to vary over time with the system Hamiltonian We use a linear ramp for the bias and hence set For times t < t 1 this is identical to the problem studied above. Hence, we choose t 1 to be large enough to allow the system to first reach thermal equilibrium with the environment and 1 is chosen such that this equilibrium state is higher in energy than the initial state and thus absorbs energy from the bath. The bias is then linearly increased, increasing the system energy along with it. The system is then in a state such that energy must be emitted into the bath for the system to be in equilibrium with respect to the new system Hamiltonian, which now has a much higher value for . This behaviour is shown in Fig. 4. As can be seen in the upper panel of the Figure during the linear ramp the energy of the system decreases, emitting energy into the bath. This can be seen in the heat distribution of the bath as a band of modes which increase in energy. The modes targeted correspond to frequencies close to the (time-dependent) Rabi fre-quencyΩ(t) = (t) 2 + Ω 2 which changes over the course of the linear ramp. This shows that by using a time varying drive, along with the intuition we have built up about the behaviour of the environment, we are able to control the (frequency dependent) absorption and emission of energy into and out of the bath.
The spin-boson model has provided a concrete example of how we can calculate the dynamics of bath modes, by combining the moment generating functional method introduced in Sec. 2 and the process tensor approach of Sec. 3. We saw how the steady state behaviour can be calculated as well as how we can control the population of different modes by using an appropriate time-dependent drive.

Conclusions & Outlook
In this paper we have introduced a method for calculating any dynamical observable of a Gaussian bath purely in terms of correlation functions of a linearly coupled system. For systems coupled to infinite baths we demonstrated that the process tensor formalism provides a framework to calculate, cheaply, the large number of system correlation functions needed. In general it is numerically costly to calculate a process tensor. However, once it exists it can then be used to completely characterise any time dependent observable or any multi-point correlation function of the system, bath, or a combination of both.
Specifically, we were able to calculate how energy is moved between different environment modes and a spin-1/2 particle in the biased spin boson model. We showed that, for weak system environment coupling, the spin reaches thermal equilibrium with the bath removing or adding energy to modes which correspond to the system frequency. We also saw how the techniques described here can allow us to track the flow of energy into and away from the system when it is subject to a time dependent driving.
The framework we presented is not specific to any technique for calculating the required correlation functions and opens up many avenues to analyze a broad range of open quantum systems. This includes simulations of any Gaussian environment, including, for example, fermionic environments. For non-Gaussian environments, the moment generating functional approach is still valid, but Wick's theorem does not hold so the derivatives in Eq. 7 must to carried out numerically.
Recent advances have extended the applicability of numerically exact process tensor techniques to many-body systems [63] and multiple environments [64]. Those techniques could be combined with the moment generating functional method we have described here to expose the dynamics of environments in such situations.
When the environment induces complex non-Markovian dynamics in a system this can lead to advantageous effects, e.g. in increased quantum channel capacities [65], entanglement distribution [66] and quantum control [67]. Our method has the potential to reveal the origin of such complicated non-Markovian dynamical effects, and to aid the design of new systems that could control and harness these. During the preparation of this work we have implemented these methods into an open source package [50] in an effort to enable others to easily apply these results in their own fields.

A Superoperator formalism
To calculate the evolution of the total system-bath density matrix we start with the von Neumann equation (in the interaction picture) and use the superoperator formalism: where L I (t) is the superoperator associated with taking the commutator with −iH I (t) i.e.
. We can then formally integrate this and arrive at the expression: To find the dynamics of reduced system we trace over the bath to leave where we have made the assumption of a separable initial state such that |ρ(0) = |ρ S (0) ⊗ |ρ B (0) and · B represents an expectation taken with respect to |ρ B (0) . Throughout this paper we represent operators in Hilbert space as lowercase letters in normal italics (e.g. s and x) and superoperators in Liouville space superoperators as uppercase calligraphic letters (e.g. L). If a superoperator corresponds to simple extension of an operator then it will be labelled as such, for example: In a similar vein we can define the left-acting superoperator aρ → A L |ρ , the right-acting superoperator ρa → A R |ρ and the anti-commutator superoperator {a, ρ} → A + |ρ . The only exception is that when referring to superoperators which correspond to Hamiltonians in Hilbert space we will use the typical L symbol. Although the time-ordering is carried out on superoperators we reserve the calligraphic T for the trace operation. In Hilbert space the trace of an operator returns a scalar. In Liouville space the density matrix is vectorised and thus to return a scalar the trace must correspond to an inner product. To see how we define this we give a brief overview of the formalism of superoperators.
Consider a Hilbert space spanned by the vectors {|e i } a density matrix can be written as a vector in the corresponding Liouville space by using the mapping An equivalent mapping can be carried out on any operator in Hilbert space. We can then define the trace as a dual vector in the Liouville space which acts as follows From this we can see that the correct dual vector is Now consider a bipartite system H A ⊗ H B spanned by {|f i ⊗ |h k }. The corresponding Liouville space is in turn spanned by {|f i ⊗ |f j ⊗ |h k ⊗ |h l }. The partial trace over H B is then We now define the partial trace over a series of superoperators as: We will reserve the use of angular brackets to refer to expectations of superoperators. If there is no subscript after the angular brackets then the trace is over the entire space. If we consider the bath to initially be in a Gaussian state, e.g. thermal equilibrium, then the trace over the bath can be simplified using Wick's theorem to give: The key result of this paper relies on an equivalent relation holding when source terms linear in the ladder operators are added to the exponent.

B Number operator derivation
In this section we derive Eq. (10) of the main text. For simplicity we consider a bath that consists of just a single mode such that H I = gs(a † + a) and H B = ωa † a, where s is an arbitrary system operator. The generalisation to a multimode bath is straightforward. We wish to calculate n(t) = Tr[a † (t)a(t)ρ] and write this quantity in terms of expectation values of superoperators as follows, Now we can write We can combine the exponentials in Eq. (46) without using the BCH formula by virtue of the identity which holds for arbitrary B, C, and integral limits. This is because the portion of the integral that contains C(t), which in general does not commute with B(t) under time ordering, is infinitesimal, i.e. has measure zero. Thus, Eq. (46) becomes with To evaluate G(Λ, Λ * ) we now split the trace over the full system in Eq. (48) into a partial trace over the bath, followed by a partial trace over the system where we have used idempotency of time ordering, Since the operational part of L I (Λ, Λ * , t ) is linear in bath operators we can use standard Gaussian integral results to evaluate the trace over the bath analytically. This gives is the usual Feynman-Vernon influence phase (i.e. the exponent in Eq. (43)) and n(0) = Tr[a † (0)a(0)ρ]. Since this is the only part of G(Λ, Λ * ) that carries dependence on Λ and Λ * we can now take the derivatives and set Λ = Λ * = 0 to obtain To evaluate Eq. (52) we now use which results in where we have used the fact that for any α and α . We have also defined to indicate expectation values taken with respect to the full system density matrix under full system evolution. To be clear, the subscript H means superoperators within the brackets are in the full system Heisenberg picture and the trace is taken with respect to the full system density matrix, the subscript S means superoperators are in the interaction picture and the trace is taken with respect to the reduced system density matrix, and the brackets without subscript mean superoperators are in the interaction picture and the trace is taken with respect to the full system density matrix. For instance, if t > t then Now we evaluate the bath expectation values in Eq. (54). For example, in the first line of Eq. (54) these evaluate to where we have used cyclicity of the trace, , and made the interaction picture time evolution of bath superoperators explicit, A α (t) = A α e −iωt . The overall result is To help us enforce time-ordering upon system superoperators we split the rectangular integral domain into two triangular domains, where the first domain has t > t and the second had t > t . With this, and again using cyclicity of the trace, e.g. S L (t )S R (t ) = S L (t )S L (t ) , we find that terms propotional to n(0) 2 cancel out and we are left with (62) as used in the main text.

C Two-time correlation functions
In this section we show it is possible to extend the above approach to not just calculate equal time correlation functions but also calculate general two-point correlation functions of bath operators. For concreteness we'll consider where t > t 1 . The expression in terms of system correlation functions is derived in a similar fashion to the number operator shown in the previous section. It can also be expressed as derivatives of a generating functional as in Eq. (48) but with: The absence of the term ∝ |Λ| 2 that arose in Eq. Taking derivatives with respect to Λ and Λ * and setting Λ = Λ * = 0 then allows us to write The derivation now proceeds by substituting the same form of L I as defined in Eq. (53) and then timeordering. The time-ordering is done by again splitting the integrals into parts with a fixed ordering i.e. From this we can see that the integrand decays as t − t and t 1 − t increase. If t 1 is large enough (and t by the requirement of t > t 1 ) then the system correlation functions, e.g. Tr[s(t )s(t )ρ], will be stationary for all t and t where the integrand is of significant weight. If we denote the stationary system correlation functions as and make the shifts t → t − t and t → t 1 − t we can express this term as The integrand now decays with increasing t and t such that in the stationary limit we can take the integral limits to infinity. Focusing now on the last term of (73) we can make the shifts t → t − t 1 and t → t − t 1 to yield Putting all of this together we can express the stationary bath correlation function as

D Toy model
In this section we benchmark our technique by calculating the energy change in a particular mode, ∆Q q , for a two-level system coupled to an environment consisting of only two bosonic modes. With so few degrees of freedom we can simulate the dynamics of the modes explicitly and compare with the prediction given by Eq. (10). The model we use is the two-mode Rabi model, the Hamiltonian of which is given by and describes a two-level system driven with strength Ω and bias . The environment consists of bosons of frequency ω q where a † q creates an excitation in the mode indexed by q ∈ {0, 1}. The system is coupled to each of these modes with strength g q .
In Fig. 5 we plot the dynamics of each mode, calculated first by numerically integrating the Schrödinger equation for the full system-environment dynamics, and then from Eq. (61). We initialised the two-level system in the spin-up state and each bosonic mode in a thermal state at temperature T = 0.1Ω. It is clear that our new method correctly predicts the environment dynamics for both modes.
Let us now introduce some dissipation into this model such that the density matrix evolves according to d dt |ρ(t) = L R |ρ(t) + γ ↓ D(a) + γ ↑ D(a † ), where the dissipative terms D(·) are as defined in Eq. (71) and we set γ ↓ = γ[n(0) + 1] and γ ↑ = γn(0). In Fig. 6 we plot the two-time correlation function C q (τ ) (of each mode q) as defined in Eq. (78). The same initial state and temperature were used as the unitary non-dissipative case. We can again see good agreement between our result and the exact dynamics. The correlation functions calculated here can readily be Fourier transformed to calculate emission spectra for these degrees of freedom therefore opening the possibility of calculating an experimentally observable quantity.

E Higher order correlation functions
In calculating higher order correlation functions it is beneficial to express the generating functional in a more concise way. This can be done in the following way: Here we have defined the vectors the matrix and φ is the same scalar function as previously which depends on the order of the ladder operators in the correlation function.
It can often be easier to simplify things in terms of these matrix elements and then substitute the explicit values. For example, let's consider the second order coherence at propagation time t with zero time delay, g (2) (t, 0) defined as g (2) (t, 0) = Tr a † (t)a † (t)a(t)a(t)ρ (84)