Efficient non-Markovian quantum dynamics using time-evolving matrix product operators

In order to model realistic quantum devices it is necessary to simulate quantum systems strongly coupled to their environment. To date, most understanding of open quantum systems is restricted either to weak system–bath couplings or to special cases where specific numerical techniques become effective. Here we present a general and yet exact numerical approach that efficiently describes the time evolution of a quantum system coupled to a non-Markovian harmonic environment. Our method relies on expressing the system state and its propagator as a matrix product state and operator, respectively, and using a singular value decomposition to compress the description of the state as time evolves. We demonstrate the power and flexibility of our approach by numerically identifying the localisation transition of the Ohmic spin-boson model, and considering a model with widely separated environmental timescales arising for a pair of spins embedded in a common environment.

T he theory of open quantum systems describes the influence of an environment on the dynamics of a quantum system 1 . It was first developed for quantum optical systems 2 , where the coupling between system and environment is weak and unstructured. In such situations, one can almost always assume that the environment is memoryless and uncorrelated with the system-that is, the Markov and Born approximations holdallowing a time-local equation of motion to be derived for the open system. The resulting Born-Markov master equation works because the environment-induced changes to the system dynamics are slow relative to the typical correlation time of the environment.
There are now a growing number of quantum systems where a structureless environment description is not justified, and memory effects 3 play a significant role. These include micromechanical resonators 4 , quantum dots 5,6 and superconducting qubits 7 , and can underpin emerging quantum technologies such as the singlephoton sources needed for quantum communication 8 . In addition, structured environments are ubiquitous in problems involving the strong interplay of vibrational and electronic states. For example, those involving the photophysics of natural photosynthetic systems 9,10 , complex organic molecules used for light emission or solar cells 11 , or semiconductor quantum dots [12][13][14][15] . Similar problems arise when considering non-equilibrium energy transport in molecular systems 16 or non-adiabatic processes in physical chemistry 17 . Non-Markovian effects can even be a resource for quantum information 18,19 .
Various approaches exist for dealing with non-Markovian dynamics 1,3 . Some particular problems have exact solutions 20 . For others, unitary transformations can uncover effective weak coupling theories, and perturbative expansions beyond the Born-Markov approximations 12,21 ; these techniques typically yield time-local equations and are limited to certain parameter regimes. Diagrammatic formulations of such perturbative expansions can also form the basis for numerically exact approaches, for example, the real-time diagrammatic Monte Carlo as implemented in the Inchworm algorithm 22,23 . Finally, there are non-perturbative methods that enlarge the state space of the system. This can be through hierarchical equations of motion 24 , through capturing part of the environment within the system Hilbert space [25][26][27] or by using augmented density tensors (ADTs) to capture the system's history 28,29 . These can be very powerful but require either specific assumptions about the environments 24,27 or resources that scale poorly with bath memory time.
In this Article, we describe a computationally efficient, general and yet numerically exact approach to modelling non-Markovian dynamics for an open quantum system coupled to an harmonic bath. Our method, which we call the time-evolving matrix product operator (TEMPO), exploits the ADT 28,29 to represent a system's history over a finite bath memory time τ c . If the bath is well behaved, then using a singular value decomposition (SVD) to compress the ADT on the fly is expected to enable accurate calculations with computational resources scaling only polynomially with τ c . We demonstrate the power of TEMPO by exploring two contrasting problems: the localisation transition in the spin-boson model (SBM) 30 and spin dynamics with an environment that has both fast and slow correlation timescales-a problem for which other methods are not available. For both these problems we observe polynomial scaling with memory time.

Results
Time-evolving matrix product operators. In this section, we outline how the TEMPO algorithm works; further details are provided in the Methods section. We start by introducing the ADT. To define the notation and our graphical representation of it, we first consider the evolution of a Markovian system, which can be described by a density operator that contains d 2 numbers for a d-dimensional Hilbert space. Usually, the density operator is written as a d × d matrix, but we instead use a length d 2 vector with elements ρ i (t). To evolve by a timestep Δ, we write where L is the Liouvillian 1 . The graphical representation of this is shown in Fig. 1a. The red circle represents the density operator, with the protruding 'leg' indicating this is a tensor of rank one, that is, a vector. This leg is indexed by an integer i = 1,…,d 2 . The blue square with two legs represents the propagator e ΔL , written as a d 2 × d 2 superoperator 1 . The matrix-vector multiplication in Eq. (1) is shown by joining a leg of the propagator to the density operator, indicating tensor contraction. This contraction generates the density operator at time t + Δ.
In order to capture non-Markovian dynamics, we extend our representation of the state at time t from a vector to an ADT, representing the history of the system. This is motivated by the path integral of a system interacting linearly with a bosonic environment. After integrating out the environment, the influence of the environment on the system can be captured by an 'influence functional' of the system paths alone 1 . The influence functional couples the current evolution to the history, and captures the non- Input state Grow P r o p a g a t e Fig. 1 Schematic description of the TEMPO algorithm. a Pictorial representation of matrix-vector multiplication. In b we show how the ADT can be decomposed into an MPS. c The full tensor network starting from an initial standard density operator which is grown to an ADT with K legs, as shown in d, where we have contracted the contents of the green box. To propagate forward one step, we contract the ADT with the next row of the propagator, as in e. A schematic representation of the spin-boson model is shown in f Markovian dynamics. Makri and Makarov 28,29 showed that by considering discrete timesteps, and writing the sum over system states in a discrete basis, the path integral could be reformulated as a propagator for the ADT, written as a discrete sum over paths. The influence functional becomes a series of influence functions I k (j, j′) that connect the evolution of the amplitude of state j to the amplitudes of states j′ an integer number, k, of timesteps ago. This approach is known as the quasi-adiabatic path integral (QUAPI).
As described so far, the ADT grows at each timestep, to record the lengthening system history. However, the influence functions have no effect once kΔ exceeds the bath correlation time τ c . One can therefore propagate an ADT containing only the previous K = τ c /Δ steps: this is the finite memory approximation. This means we consider an ADT of rank K, written as A i 1 ;i 2 ; ;i K ðtÞ, where each index runs over i k = 1,…,d 2 . The explicit construction of this tensor is described in the Methods section. In general A i 1 ;i 2 ; ;i K ðtÞ contains d 2K numbers, which scales exponentially with the correlation time τ c . If the full tensor is kept, one quickly encounters memory problems, and typical simulations are restricted to K < 20 31,32 . Improved QUAPI algorithms 33,34 show that (for some models) typical evolution does not explore this entire space, leading us to seek a minimal representation of the ADT.
Matrix product states (MPS) 35,36 are natural tools to represent high-rank tensors efficiently where correlations are constrained in some way. Examples include the ground state of one-dimensional (1D) quantum systems with local interactions 37 , steady state transport in 1D classical systems 38 or time-evolving 1D quantum states 39 . Inspired by these results, we show how an ADT can be efficiently represented and propagated using standard MPS methods. One may decompose high-rank tensors into products of low-rank tensors using SVDs and truncation. By combining indices, the tensor A can be written as 36 : Here, U, V are unitary matrices, and λ α denotes a singular value of the matrix A. Truncation corresponds to throwing away singular values λ α smaller than some cutoff λ c , consequently reducing the size of the matrices U, V. This procedure can be iterated by sweeping k across the whole tensor. The result of this is shown graphically in Fig. 1b, and can be written as: This provides an efficient representation of the state, with a precision controlled by λ c . A i 1 ;i 2 ; ;i K ðtÞ can be time locally propagated using a tensor B j i ; ;j K i 1 ; ;i K . Crucially, this propagation can be performed directly on the matrix product representation of A. Moreover, the tensor product description of B j i ; ;j K i 1 ; ;i K , shown as the connected blue squares in Fig. 1c, has a small dimension, d 2 , for the internal legs. Similarly to the time evolution shown in Fig. 1a, the state A(t + Δ) is generated by contracting the legs of A(t) with the input legs of B. Contracting a tensor network with a MPS, and truncating the resulting object by SVDs is a standard operation 36 . In all the applications we discuss below, we find that as time propagates we are able to maintain an efficient representation of A i 1 ;i 2 ; ;i K ðtÞ with precision determined by λ c .
The structure of the propagator depends on the influence functions I k (j, j′) as shown in Fig. 1c (see also Methods section). We use darker colours to represent influence functions corresponding to more recent time points, which are expected to generate stronger correlations in the ADT. The input and output legs of the propagator are offset in the figure, so time can be viewed as propagating from left to right. In effect, at each step the register is shifted so that the right-most output index corresponds to the new state: events that occurred more than τ c ago are dropped, as illustrated by the white semicircles in Fig. 1, since they do not influence the future evolution. Evolution over a series of timesteps is depicted in Fig. 1c-e. In Fig. 1c we show the full tensor network. Assuming the initial state of the system is uncorrelated with its environment means it can be drawn as a regular density operator. In the 'grow' phase, a series of asymmetric B propagators are applied, which allow the relevant system correlations to extend in time. Once the system has grown to an object with K legs, we enter the regular propagation phase, shown in Fig. 1d, e.
Spin-boson phase transition. To demonstrate the utility of the TEMPO algorithm, we apply it to two problems of a quantum system coupled to a non-Markovian environment. We first consider the unbiased SBM 30 , which has long served as the proving ground for open system methods. The generic Hamiltonian of this model is where the S i are the usual spin operators, a y i (a i ) and ω i are, respectively, the creation (annihilation) operators and frequencies of the ith bath mode, which couples to the system with strength g i . The behaviour of the bath is characterised by the spectral density function This model is known to show a rich variety of physics depending on the particular form of spectral density and system parameters chosen. When the spectral density is Ohmic, JðωÞ ¼ 2αω expðÀω=ω c Þ, the model is known to exhibit a quantum phase transition in the BKT universality class 40 , at a critical value of the system-environment coupling α = α c 30,41 . The transition takes the system from a delocalised phase below α c , where any spin excitation decays (〈S z 〉 = 0 in the steady state), to a localised phase above α c (〈S z 〉 ≠ 0 in the steady state). Most analytic results are restricted to the regime where the cutoff frequency ω c ) Ω. For example, when S describes a spin-1/2 particle, the phase transition occurs at α c ¼ 1 þ OðΩ=ω c Þ 30,40,42 .
We are able to explore the dynamics around this phase transition using TEMPO. In Fig. 2a we show the polarisation dynamics of the spin-1/2 SBM for a range of α at K = 200. This memory length is an order of magnitude larger than standard ADT implementations 30 and is required to reach the asymptotic limit of the dynamics in the vicinity of the phase transition. We achieve convergence by varying the timestep Δ and SVD cutoff λ c . We take an initial condition 〈S z 〉 = + 1/2 with no excitations in the environment, and find 〈S z (t)〉.
Before reaching the localisation transition at α = α c , one first reaches a crossover at α ' 0:5 from coherent decaying oscillations to incoherent decay 29 . For α > 0.5, we find 〈S z 〉 always decays to zero asymptotically as hS z ðtÞi / expðÀγtÞ to a very good approximation; fits to this function are shown as dashed lines in Fig. 2a. Decay to zero for all α > 0.5 conflicts with the existence of a localised phase at large α, where 〈S z 〉 should asymptotically approach a non-zero value. The origin of this discrepancy is the finite memory approximation, which produced a time-local equation in the enlarged space of K timesteps. Timelocal dynamics of a finite system typically generates a gapped spectrum of the effective Liouvillian 43 . In the localised phase, α > α c , the spectral gap should vanish asymptotically as we increase the memory cutoff τ c → KΔ. We should thus examine how the extracted decay rate, γ, depends on the memory cutoff. For α < α c , γ should remain finite as τ c → ∞, while for α > α c it should vanish. In Fig. 2b we plot γ as a function of 1/K =Δ/τ c for different values of α around the phase transition. At small α, γ does appear to remain finite as K → ∞, while at large α the behaviour appears consistent with localisation.
We may estimate the location of the phase transition by extrapolating 1/K → 0 for each α, and finding the smallest value of α consistent with γ → 0. To do this, we use cubic fits in Fig. 2b (solid lines), and extract the constant part, with the restriction that the extracted γ cannot be negative. In order to find the phase transition as accurately as possible, we must perform simulations up to very large values of K: we here perform simulations up to K = 200, something that would be simply impossible without the tensor compression we exploit. Errors in our fits are assessed by monitoring the sensitivity of the best-fit result to truncation precision λ c . These errors are all <10 −4 and so are smaller than the points in Fig. 2. This allows us to find an error in the extracted K → ∞ limit. The extracted values for γ are displayed in Fig. 2d where we show our estimate for its 68 and 95% confidence intervals. These suggest that α c ' 1:25, consistent with the known analytic results 30,40,42 . We note that identifying α c precisely from the time dependence of 〈S z 〉 is particularly challenging: since the localisation transition is in the BKT class 40 , the order parameter approaches zero continuously.
The efficiency of TEMPO enables consideration of models with a larger local Hilbert space. To demonstrate this, we examine the localisation transition in the spin-1 SBM. Physically this could either arise from a spin-1 impurity or from a pair of spin-1/2 particles interacting with a common environment 44 . On switching to this problem, the local dimension of each leg of our state tensor increases from d 2 = 4 to d 2 = 9, reducing the values of K we can reach. However, we also find convergence occurs for larger timesteps, allowing access to similar values of τ c .
In Fig. 3a we show the dynamics of this model, after initialising to 〈S z 〉 = 1. In this case, on both sides of the localisation transition, the dynamics shows complex oscillatory behaviour before settling down to an exponential decay. This introduces more uncertainty to our exponential fits. However, as shown in Fig. 3b the extracted decay rate vanishes at α c ' 0:28, indicative of the phase transition and agreeing with numerical renormalisation group results 44,45 , but in contrast to the results found using a variational ansatz 46 .
Two spins in a common environment. We next demonstrate the flexibility of TEMPO by applying it to a dynamical problem for which other methods are not available. We consider a pair of identical spins-1/2, at positions r a and r b , which couple directly to each other through an isotropic Heisenberg coupling Ω, and which both couple to a common environment, see Fig. 4a. The Hamiltonian reads: The system-bath coupling constants have a position-dependent phase, g i;ν ¼ g i e Àik i Ár ν , where k i is the wavevector of the ith bosonic mode. We assume linear dispersion ω i = c|k i | and c = 1. This model exhibits complex dissipative dynamics on two different timescales. The faster timescale describes dissipative dynamics of the spins due to interactions with their nearby environment, typically set by the ω c defined earlier. The other timescale is set by the spin separation R = |r a − r b | over which there is an environment-mediated spin-spin interaction. By The dependence of the decay rate of the exponential fit on 1/K. This allows us to analyse the behaviour as K → ∞. c The change in the decay rate as we go through the transition by varying α for the values of K indicated. In d, we give 68% (blue) and 95% (red) confidence intervals for the extrapolated decay rate which crosses zero at around α c ' 1:25. The bath cutoff frequency is ω c = 5 and everything is measured in units of the Hamiltonian driving term Ω changing R we can control the ratio of these timescales. The dimension, D, of the bath also has an effect: the intensity of environmental excitations propagating from one spin to the other will be stronger for lower D.
When the spins are close together, R<ω À1 c , it is difficult to distinguish local dissipative effects from the environmentmediated interaction and both master equation techniques 13 and the standard ADT method 47 generate accurate dynamics. Instead, we consider large separation R>ω À1 c , about which little is known. The ADT then requires both a small timestep Δ ( ω À1 c to capture the fast local dissipative dynamics and a large cutoff time τ c = KΔ > R to capture environment-induced interactions; hence, a very large K is needed. Using TEMPO we are able to investigate these dynamics without even having to go beyond the tensor growth stage shown in Fig. 1c, and thus avoid any error caused by a finite memory cutoff K. We project onto the S z,a + S z,b = 0 subspace of the system, consisting of the two anti-aligned spin states, since this is the only sector with non-trivial dynamics. The effective Hamiltonian for this 2d subspace can then be mapped onto the spin-1/2 SBM, Eq. (4), albeit with a modified spectral density that depends on R. Details of this procedure are given in the Methods section.
In Fig. 4b, c we show dynamics for different R for environments with D = 1 and D = 3. Insets show the effective spectral densities, J(ω), and real part of the bath autocorrelation functions, C(t), which we define in the Methods section. We initialise the spins in a product state with 〈S z,a 〉 = 1/2, 〈S z,b 〉 = −1/2 and calculate the probability, P(t), of finding the system in this state at time t. The bath is initialised in thermal equilibrium at temperature T. For D = 1, after initial oscillations decay away over a timescale $ ω À1 c , there are revivals at t = R. This is due to the strongly oscillating spectral density which results in a large peak at C(t = R). As expected for a one-dimensional environment, the profile of these secondary oscillations is independent of R when R ) ω À1 c . Additionally for R = 20 more small amplitude oscillations appear at t ≈ 40, due to the effective interaction of the spins at t ≈ 20 sending more propagating excitations into the environment. For D = 3 the spectral density still has an oscillatory component, though it is much less prominent. The resulting peaks at C(t = R) are thus much smaller than the t = 0 peak and have only a small effect on the dynamics. Small amplitude oscillations can be seen at t ≈ R when R = 8, but with R = 16 it is difficult to see any significant features in the dynamics.

Discussion
We have presented a highly efficient method for modelling the non-Markovian dynamics of open quantum systems. Our method is applicable to a wide variety of situations. In well-established ADT methods, non-Markovianity is accounted for by encoding the system's history in a high-rank tensor; we have overcome the restrictive memory requirements of storing this tensor by representing it as an MPS. We can then efficiently calculate open system dynamics by propagating this MPS via iterative application of an MPO. To test our technique we used it to find the localisation transition in the SBM, for both spin-1/2 and spin-1, and found estimates for the critical couplings, consistent with other techniques. We then applied our method to a pair of interacting spins embedded within a common environment, in a regime where a large separation of timescales prevents the use of other methods. Precisely locating the phase transition is a rigorous test of any numerical method: as we found, very large memory times, up to K = 200 were required to precisely locate this point. Other improved numerical methods 22,23,33 have demonstrated a degree of enhanced efficiency when considering conditions away from the critical coupling. As yet, other such general methods have not been used to precisely locate the transition.
The key to our technique is that tensor networks provide an efficient representation of high-dimensional tensors encoding restricted correlations. As well as the widespread application of such methods in low-dimensional quantum systems [35][36][37][38][39] , they have also been applied to sampling problems in classical statistical physics 48 , and analogous techniques (under the name 'Tensor trains') have been developed in computer science 49 . Moreover, there has been a recent synthesis showing how techniques developed in one context can be extended to others, such as machine learning 50 , or Monte Carlo sampling of quantum states 51 . Our work defines a further application for these methods, and future work may yet yield even more efficient approaches.
The methods described in this article are already very powerful in their ability to model general non-Markovian environments. They also enable easy extension to study larger quantum systems, by adapting other methods from tensor networks such as the optimal boson basis 52 -these will be the subject of future work. They may also be combined with approaches such as the tensor transfer method described in Ref. 53 . This method allows efficient long time propagation of dynamics, so long as an exact map is known up to the bath memory time: TEMPO enables efficient calculation of the required exact map. With such tools available, the study of the dynamics of quantum systems in non-Markovian environments 3 can now move from studying isolated examples to elucidating general physical principles, and modelling real systems.

Methods
TEMPO algorithm. In this section, we will present the details of the TEMPO algorithm, paying particular attention to how the ADT and propagator are constructed in a matrix product form.
The generic Hamiltonian of the models we consider is where H 0 is the (arbitrary) free system Hamiltonian and H E contains both the bath Hamiltonian and the system-bath interaction. Here a y i (a i ) and ω i are the creation (annihilation) operators and frequencies of the ith environment mode. The system operator O couples to bath mode i with coupling constant g i . As outlined in the main text, we work in a representation where d × d density operators are given instead by vectors with d 2 elements. These vectors are then propagated using a Liouvillian as in Eq. (1) of the main text, L ¼ L 0 þ L E , where L 0 and L E generate coherent evolution caused by H 0 and H E , respectively. It has been shown recently that it is straightforward to include additional Markovian dynamics in the reduced system Liouvillian 54 in the ADT description.
If the total propagation over time t N is composed of N short time propagators e t N L ¼ ðe ΔL Þ N , we can use a Trotter splitting 55 e ΔL % e ΔL E e ΔL 0 þ OðΔ 2 Þ: ð9Þ We note that the following arguments can be easily adapted to use the higherorder, symmetrized, Trotter splitting 28,29,56 that reduces the error to Δ 3 . All the numerical results presented use this symmetrized splitting, but for ease of exposition we use the form of Eq. (9) here. We assume the initial density operator factorises into system and environment terms, with the environment initially in thermal equilibrium at temperature T. Time evolution can then be written as a path sum over system states, by inserting resolutions of identity between each e ΔL E e ΔL 0 and then tracing over environmental degrees of freedom. The result is the discretized Feynman-Vernon influence functional 28,29 , which yields the following form for the time evolved density matrix: The indexing here is in a basis where O is diagonal. Each j index runs from 1 to d 2 , and due to the order of the splitting in Eq. (9), the initial state of the system has been propagated forward a single timestep, ρ j 1 ðΔÞ ¼ e ΔL 0 Â Ã j 1 j 0 ρ j 0 ð0Þ. We have defined the influence functions I k ðj; j′Þ ¼ e ϕ k ðj;j′Þ ; k ≠ 1; Here O À j are the d 2 possible differences that can be taken between two eigenvalues of O and O þ j the corresponding sums. The coefficients, η k , quantify the non-Markovian correlations in the reduced system across k timesteps of evolution and are given by the integrals dt′′Cðt′ À t′′Þ; n ≠ n′; R t n t nÀ1 dt′ R t′ t nÀ1 dt′′Cðt′ À t′′Þ; n ¼ n′; where C(t) is the bath autocorrelation function with temperature measured in units of frequency and with the spectral density JðωÞ ¼ P i jg i j 2 δðω i À ωÞ. The summand of the discretised path integral in Eq. (10) can be interpreted as the components of an N-index tensor A j N ;j NÀ1 ; ;j 1 . This tensor is an ADT of the type originally proposed by Makri and Makarov 28,29 . We will show below that this N-index tensor can also be written as tensor network consisting of N(N + 1)/2 tensors with, at most, four legs each and that this network can be contracted using standard MPS-MPO contraction algorithms 35,36 , First we gather terms in the inner piece of the double product in Eq. (10) into a single object, which we write as components of an n-index tensor Next, we define the (2n − 1)-index tensors B j n ;j nÀ1 ; ;j 1 i nÀ1 ; ;i 1 for n > 1, and the 1-index initial ADT: We may now evolve this ADT in time iteratively by successive contraction of tensors. This process is shown graphically in Fig. 1c. The first contraction produces a 2-index ADT which describes the full state and history at the second time point: We next contract with B j 3 ;j 2 ;j 1 i 2 ;i 1 to produce a 3-index ADT and so on. The nth step of this process then looks like A j n ;j nÀ1 ; ;j 1 ¼ B j n ;j nÀ1 ; ;j 1 i nÀ1 ; ;i 1 A i nÀ1 ;i nÀ2 ; ;i 1 ; ð19Þ and the density operator for the open system at time t n = nΔ is recovered by summing over all but the j n leg, ρ j n ðt n Þ ¼ X j nÀ1 ; ;j 1 A j n ;j nÀ1 ; ;j 1 ; ð20Þ from which observables can be calculated. At each iteration the size of the ADT grows by one index, since up to now we have made no cutoff for the bath memory time: we are in the 'grow' phase depicted in Fig. 1c. To compress the state after each application of this B tensor, we sweep along the resulting ADT performing SVD's and truncating at each bond, throwing away the components corresponding to singular values smaller than our cutoff λ c . This gives an MPS representation of the ADT, as given in Eq. (3). As discussed in Ref. 57 , we must in fact sweep both left to right and then right to left to ensure the most efficient MPS representation is found. If no bath memory cutoff is made, this whole process is repeated until the final time point is reached at n = N. The (2n − 1)-index propagation tensor, B, can be represented as an MPO such that the above process of iteratively contracting tensors becomes amenable to standard MPS compression algorithms 35,36 . The form required is B j n ;j nÀ1 ; ;j 1 i nÀ1 ; ;i 1 where we define the rank-4 tensor and the rank-2 and rank-3 tensors appearing at the ends of the product are Upon substituting these forms, Eqs. (22)-(24), into Eq. (21) it is straightforward to verify that we recover the expression Eq. (16). The rank-(2n − 1) MPO, B j n ;j nÀ1 ; ;j 1 i nÀ1 ; ;i 1 , is represented by the tensor network diagram in Fig. 5.
We note it has recently been shown that if the spectrum of O has degeneracies, then part of the sum in Eq. (10) can be performed analytically, vastly reducing computational cost of the ADT method for systems where the environment only couples to a small subsystem 58 . Here we can further exploit the fact that, even when there is no degeneracy in the d eigenvalues of O, there is always degeneracy in the d 2 differences between its eigenvalues, O À j , that is, d of these differences are always zero. Using the same partial summing technique described in Ref. 57  The finite memory approximation can now be introduced by throwing away information in the ADT for times longer than τ c = KΔ into the system's history. To do this we write ½b k α; j Thus, when propagating A j n ; ;j 1 beyond the Kth timestep only indices j n to j n−K+1 have any relevance and we can sum over the rest. The way we do this in practice is to define the 2K-leg tensor MPO such that contraction with a rank-K MPS is equivalent to first growing the MPS by one leg and then summing over (i.e. removing) the leg which is earliest in time.
Repeating this contraction propagates an A-tensor MPS forward in time, but maintains its rank of K for all timesteps n > K. This is what we show in the 'propagate' phase of Fig. 1c. For some spectral densities, it is possible to improve the convergence with τ c by making a softer cutoff 59,60 , but since TEMPO can go to very large values of K this is not necessary here. For time-independent problems (as we study here), the 'propagate' phase involves repeated contraction with the same MPO, Eq. (26), which is independent of the timestep. To make this clear, it is convenient to change our index labelling (which, so far has referred to the absolute number of timesteps from t = 0). We will instead relabel the indices on the MPO and MPS as follows: B j Kþ1 ; ;j 2 i K ; ;i 1 ! B j 1 ; ;j K i 1 ; ;i K and A j n ; ;j nÀKþ1 ! A j 1 ; ;j K ðt n Þ. The indices now refer to the distance back in time from the current time point. To summarise, with the new labelling we first grow the initial state into a K-index MPS, A j 1 ; ;j K ðτ c Þ, and then propagate as: and the physical density operator is found via Having described the TEMPO algorithm we now briefly analyse the computational cost of applying it to the SBM of Eq. (4). In Fig. 6a we plot the total size, N tot , of the MPS and maximum bond dimension, λ max , used to obtain converged results in Fig. 2 against coupling strength with K = 200. We find the most computationally demanding regime to be around α = 0.5, the point of crossover from underdamped to overdamped oscillations of S z . We find the CPU time required is linear in the total memory requirement. For the largest memory required (at α = 0.5), the time to obtain 500 data points using TEMPO on the HPC Cirrus cluster was ≈20.5 h. In Fig. 2b we show how N tot grows with K for different values of α. For α = 0.1,0.5 we see quadratic growth with K, while for couplings near and above the phase transition, α = 1,1.5, the growth is only linear. Both cases ARTICLE thus represent polynomial scaling, a substantial improvement on the exponential scaling of the standard ADT method for which one has N tot = 4 K .
Mapping two spins in a common environment to a single spin model. We show here how to map Eq. (6) describing a pair of spin-1/2 particles in a common environment onto Eq. (4), a single spin-1/2 SBM. The Hamiltonian Eq. (6) has the property that the total z-component of the two-spin system is conserved, [S z,a + S z,b , H] = 0. Thus, the problem can be separated into three distinct subspaces: the two states with the spins anti-aligned (S z,a + S z,b = 0) form one subspace and the two aligned spin states (S z,a + S z,b = ±1) are the other two. The onedimensional subspaces with aligned spins cannot evolve in time; hence, all nontrivial dynamics in this model happen in the S z,a + S z,b = 0 subspace. We therefore focus on this subspace. By doing so, we may subtract a term proportional to S z,a + S z,b from the system-bath coupling in Eq. (6). The remaining system-bath interaction is given by The effective coupling here isg i j j ¼ g i;a À g i;b ¼ 2g i sin k i Á ðr a À r b Þ=2 ½ . These couplings lead to a modified effective spectral density 13,61 , where J p (ω) is the actual density of states of the bath. The function F D (ωR) arises from angular averaging in D-dimensional space, and so crucially depends on the dimensionality of the environment. Specifically we have: where J 0 (x) is a Bessel function. We note that F D (ωR) → 0 as R → ∞ for D > 1, due to the diminishing effect of the environment-induced coupling in higher dimensions.
(When considering R → ∞, we should note that in the original Hamiltonian we neglected any retardation in the Heisenberg interaction.) At small separations, R → 0, F D (ωR) → 1 and so J(ω) → 0 for all D due to the loss of relative phase shift between the couplings of the anti-aligned states to the environment. For the bare density of states J p (ω), we consider a simple model of e.g. a quantum dot in a phonon environment, for which the coupling constants appearing in the Hamiltonian, Eq. (6), have g i $ ffiffiffiffi ffi ω i p 19 . This means that in the continuum limit the spectral density for a D-dimensional environment is where ω c describes a high-frequency cutoff and α is the strength of the interaction with the environment.