Exact Quantum Dynamics in Structured Environments

TL;DR

This paper demonstrates that the time-evolving matrix product operator method can accurately simulate the complex dynamics of two quantum systems in structured environments, revealing how environmental tuning influences system behavior.

Contribution

It introduces the application of the time-evolving matrix product operator approach to structured environments with spatial correlations, enabling precise simulation of complex quantum dynamics.

Findings

Tuning system separation to environmental wavelength alters dynamics significantly.

Full bath dynamics can be inferred from system dynamics.

The method handles long-lived bath correlations effectively.

Abstract

The dynamics of a wide range of technologically important quantum systems are dominated by their interaction with just a few environmental modes. Such highly structured environments give rise to long-lived bath correlations that induce complex dynamics which are very difficult to simulate. These difficulties are further aggravated when spatial correlations between different parts of the system are important. By modeling the dynamics of a pair of two-level quantum systems in a common, structured, environment we show that a recently developed numerical approach, the time-evolving matrix product operator, is capable of accurate simulation under exactly these conditions. We find that tuning the separation to match the wavelength of the dominant environmental modes can drastically modify the system dynamics. To further explore this behavior, we show that the full dynamics of the bath can be calculated directly from those of the system, thus allowing us to develop intuition for the complex system dynamics observed.

Tables

Table 1. Table 1: Convergence parameters used for figures in the main text.

Result	Time-step, $\Delta t$	Precision, $\chi$
Fig. 2	0.2/$\Omega$	${10}^{-8}$
Fig. 4	0.1/${\omega }_0$	${10}^{-7}$

Equations

$$\begin{gathered} H=\sum _i{ϵ}_i\text{outerproduct}{X}_i{X}_i+\frac{\Omega }{2}\left(\text{outerproduct}{X}_1{X}_2+\text{outerproduct}{X}_2{X}_1\right) \\ +\sum _i{\sigma }_{z,i}\sum _k(g_{i,k}{a}_k+g_{i,k}^{\ast }{a}_k^{†})+\sum _k{\omega }_k{a}_k^{†}{a}_k. \end{gathered}$$

$$J(\omega )=\sum _k|g_k{|}^2\delta (\omega -{\omega }_k),$$

$$C(t)={\int }_0^{\infty }d\omega J(\omega )\left(\coth \left(\frac{\omega }{2T}\right)\cos (\omega t)-i\sin (\omega t)\right).$$

$$J_0(\omega )=\frac{\alpha \Gamma {\omega }_0^2\omega }{({\omega }_0^2-{\omega }^2{)}^2+{\Gamma }^2{\omega }^2}.$$

$$H=\frac{\Omega }{2}{\sigma }_x+\frac{ϵ}{2}{\sigma }_z+{\sigma }_z\sum _k({\tilde{g}}_k{a}_k+{\tilde{g}}_k^{\ast }{a}_k^{†})+\sum _k{\omega }_k{a}_k^{†}{a}_k,$$

$$J(\omega )=2J_0(\omega )(1-\cos (\omega R))$$

$$⟨a_k(t)⟩=-i{g}_k{e}^{-i{\omega }_{k}t}{\int }_0^{t}d{t}^{\prime }{e}^{i{\omega }_k{t}^{\prime }}⟨{\sigma }_z(t^{\prime })⟩.$$

$$\Phi (x,t)=\sum _k\frac{1}{\sqrt{{\omega }_k}}(⟨a_k(t)⟩e^{ikx}+⟨a_k^{†}(t)⟩e^{-ikx}).$$

$$\begin{gathered} \Phi (x,t)=8{\int }_0^{\infty }d\omega [\sqrt{\frac{J_0(\omega )}{\omega }}\sin \left(\frac{\omega R}{2}\right)\sin (\omega x) \\ \times {\int }_0^{t}d{t}^{\prime }\sin (\omega (t-t^{\prime }))⟨{\sigma }_z(t^{\prime })⟩] \end{gathered}$$

$$H_I={\sigma }_z\sum _k{g}_k(a_k^{†}+a_k)\to \lambda {\sigma }_z(c^{†}+c).$$

$$H_{RC}=H_S+H_I+H_B+H_C$$

$$H_S=\Omega {\sigma }_x+ϵ{\sigma }_z+\lambda {\sigma }_z(c^{†}+c)+{\Omega }_{RC}{c}^{†}c$$

$$H_I=(c^{†}+c)\sum _k{d}_k(b_k^{†}+b_k)$$

$$H_B=\sum _k{\nu }_k{b}_k^{†}{b}_k$$

$$H_C=(c^{†}+c{)}^2\sum _k\frac{c_k^2}{{\omega }_k}$$

$$\begin{gathered} \frac{\partial \rho }{\partial t}=-i[H_S,\rho (t)]-\gamma {\int }_0^{\infty }d\tau {\int }_0^{\infty }d\omega \omega \cos (\omega \tau )\coth \left(\frac{\omega }{2T}\right)[\hat{C},[\hat{C}(-\tau ),\rho (t)]] \\ -\gamma {\int }_0^{\infty }d\tau {\int }_0^{\infty }d\omega \omega \cos (\omega \tau )[\hat{C},\{[\hat{C}(-\tau ),H_S],\rho (t)\}] \end{gathered}$$

$$\hat{C}(t)=\sum _{j,k=1}^{2n}{C}_{jk}{e}^{i{\omega }_{jk}t}\text{outerproduct}{\varphi }_j{\varphi }_k$$

$$\frac{\partial \rho }{\partial t}=-i[H_S,\rho (t)]-[\hat{C},[\hat{\chi },\rho (t)]]+[\hat{C},\{\hat{\Xi },\rho (t)\}]$$

$$\hat{\chi }$$

$$\hat{\Xi }$$

$$A^{j_N\dots j_1}(t_N)=\prod _{n=2}^N{B}_{i_{n-1},...,i_1}^{j_n,j_{n-1},...,j_1}{\rho }_{i_1}(t_1).$$

$$A^{j_N\dots j_1}(t_N)=\sum _{{\alpha }_1\dots {\alpha }_N}[a^{j_N}{]}_{{\alpha }_N}[a^{j_{N-1}}{]}_{{\alpha }_N,{\alpha }_{N-1}}\dots [a^{j_1}{]}_{{\alpha }_1}.$$

$$H_I(t)+H_B=\hat{A}(t)\sum _k(g_k{a}_k^{†}+g_k^{\ast }{a}_k)+\sum _k{\omega }_k{a}_k^{†}{a}_k$$

$$\frac{d}{dt}{a}_k(t)$$

$$=-i{g}_k\hat{A}(t)-i{\omega }_k{a}_k(t)$$

$$e^{i{\omega }_{k}t}\frac{d}{dt}{a}_k(t)+i{\omega }_k{e}^{i{\omega }_{k}t}{a}_k(t)$$

$$\frac{d}{dt}\left[e^{i{\omega }_{k}t}{a}_k(t)\right]$$

$$⟨a_k(t)⟩$$

$$⟨a_k^{†}(t)⟩$$

$$\Phi (x,t)=\sum _k\frac{1}{\sqrt{{\omega }_k}}\left[⟨a_k(t)⟩+⟨a_{-k}^{†}(t)⟩\right]e^{ikx}.$$

$$\Phi (x,t)$$

$$=-\sum _k\frac{2{\tilde{g}}_k}{\sqrt{{\omega }_k}}{e}^{ikx}{\int }_0^{t}d{t}^{\prime }\sin [{\omega }_k(t-t^{\prime })]⟨{\sigma }_z(t^{\prime })⟩$$

$$=-4i{\int }_{-\infty }^{\infty }dk\frac{g_k}{\sqrt{{\omega }_k}}\sin \left(\frac{kR}{2}\right)\left[\cos (kx)+i\sin (kx)\right]{\int }_0^{t}d{t}^{\prime }\sin [{\omega }_k(t-t^{\prime })]⟨{\sigma }_z(t^{\prime })⟩$$

$$=8{\int }_0^{\infty }d\omega \sqrt{\frac{J_0(\omega )}{\omega }}\sin \left(\frac{\omega R}{2}\right)\sin (\omega x){\int }_0^{t}d{t}^{\prime }\sin [\omega (t-t^{\prime })]⟨{\sigma }_z(t^{\prime })⟩$$

Full text

When a quantum system interacts with a structured environment in which a narrow band of modes dominate, the resulting dynamics can be very complex and difficult to simulate accurately. Such environments can retain a memory of their interactions with the system on a timescale comparable to that on which the state of the system changes, and in the face of such memory effects standard open system techniques can fail. However, reliable simulation of such environments is vital in order to understand the behavior of an ever-increasing number of experimental platforms and quantum devices. For example, in photosynthetic systems the observation of quantum coherence comes from strong electronic and vibrational coupling in biomolecules Turner et al. (2011), where the interplay of these degrees of freedom has been shown to enhance energy transport Thorwart et al. (2009); Chin et al. (2013); del Rey et al. (2013); Mohseni et al. (2014). Similar physics can also be found in e.g. trapped ions Maier et al. (2019) and molecular junctions Segal and Agarwalla (2016). Further, bath memory plays a key role in the function of micromechanical resonators, where adjustment of the environmental noise spectrum is possible Gröblacher et al. (2015) and photonic crystals, where band gaps in the spectral density give rise to localized modes and dissipationless oscillations Zhang et al. (2012); Goban et al. (2014); Sundaresan et al. (2019). The effect of structured environments has also been explored with superconducting qubits subject to different kinds of noise Potočnik et al. (2018).

Memory effects can be even more pronounced when the system coupled to the structured environment described above is extended. For example if it consists of several spatially separated sites. The dynamics then become sensitive to groups of environmental modes with a narrower bandwidth, and so a longer correlation time. This sensitivity was recently observed in an experiment on a ‘giant’ superconducting atom Andersson et al. (2018). Spatial correlation effects have also been studied in the context of quantum transport Nalbach et al. (2010); Nazir (2009); McCutcheon and Nazir (2011); Abramavicius and Mukamel (2011); Sarovar et al. (2011), and spatially correlated noise in quantum registers has been shown to affect quantum error propagation Aharonov et al. (2006).

In order to model environmental memory effects, we must go beyond the standard Born-Markov approximations Breuer and Petruccione (2002) which lead to a simple time-local master equation for the system density matrix. The development of techniques to simulate these non-Markovian dynamics has therefore been the subject of much relatively recent theoretical effort de Vega and Alonso (2017). These techniques broadly fall into two categories: First, there are approximate methods which change the boundary between the system and environment such that the Born-Markov approximations are valid for the new system Silbey and Harris (1984); Harris and Silbey (1985); Jang et al. (2008); Iles-Smith et al. (2014); Iles-Smith (2015); Iles-Smith et al. (2016); second, there are numerically exact methods which utilize some particular structure of the bath Hamiltonian to provide exact dynamics Tanimura and Kubo (1989); Makri and Makarov (1995a, b); Prior et al. (2010); Strathearn et al. (2018); Mascherpa et al. (2019). These approaches usually work best for models where the system Hilbert space is quite small, or has a specific form of spectral density. However, it can often be difficult to know a priori the range of validity of some of these approaches, and so accurate and efficient benchmarking procedures are essential.

In this Letter, we demonstrate that a recently developed exact approach, the Time-Evolving Matrix Product Operator (TEMPO) Strathearn et al. (2018) can be used to efficiently find non-local system dynamics for highly structured environments. To do this, we study a model consisting of a pair of spatially separated two-level subsystems interacting with a common environment, in which a narrow band of modes dominate the interaction. The complex environmental structure and spatial correlations give rise to highly non-Markovian dynamics, which we are able to explore using TEMPO. In addition, we find that by tuning the separation between the subsystems, it is possible to control their interaction with the environment.

An advantage of techniques that include some environmental modes in the system description is that they enable the extraction of some of the bath’s dynamics, which can in turn lead to insight into the system’s behavior through understanding system-environment correlations Iles-Smith et al. (2014); Chin et al. (2013). However, we will show here that exact environmental dynamics can be easily and efficiently extracted from TEMPO, leading us to a full picture of the complex quantum dynamics of this model.

The Hamiltonian for our model is given by:

[TABLE]

The two-level system (TLS) at site $i$ and position $r_{i}$ has energy splitting $\epsilon_{i}$ and an excited (ground) state denoted $\ket{X_{i}}(\ket{0_{i}})$. The pair of TLSs form a dimer with a coherent coupling of strength $\Omega$, and the environment consists of bosons confined to one dimension: $a_{k}^{\dagger}$ creates an excitation in the mode with wave vector $k$. We also introduce the Pauli operators $\sigma_{z,i}=\outerproduct{X_{i}}{X_{i}}-\outerproduct{0_{i}}{0_{i}}$. The coupling $g_{i,k}=g_{k}\exp(-ikr_{i})$ consists of a position independent part $g_{k}$ and a phase that depends on the site $i$. This kind of model could underpin a wide range of physical systems, for example biological or molecular systems undergoing energy transport and interacting with vibrational modes del Rey et al. (2013), superconducting qubits in microwave resonators Zheng and Baranger (2013), or quantum dots interacting with a micromechanical resonator Yeo et al. (2013).

The bath can be completely characterized by its spectral density

[TABLE]

from which we can calculate its autocorrelation function in thermal equilibrium at temperature $T$:

[TABLE]

If this function decays slowly compared to the system timescales then the bath memory is important.

We now consider a bath in which a narrow band of modes dominate the interaction with the system, giving rise to the spectral density Mukamel (1995)

[TABLE]

Here $\alpha$ gives the coupling strength, $\omega_{0}$ is the frequency characterizing the dominant mode and $\Gamma$ provides a measure of the width of the dominant band.

The Hamiltonian, Eq. (1), does not couple states with a different total number of system excitations and so we restrict ourselves to the single excitation subspace, the only one where non-trivial dynamics are observed. Within this subspace it is possible to map the problem onto that of a single spin coupled to a bosonic environment Stace et al. (2005). The mapped Hamiltonian is then a spin-boson Hamiltonian

[TABLE]

where $\epsilon=\epsilon_{2}-\epsilon_{1}$ and $\tilde{g}_{k}=2ig_{k}\sin(k(r_{1}-r_{2})/2)$. We also define new Pauli operators: $\sigma_{z}=(\sigma_{z,2}-\sigma_{z,1})/2$ and $\sigma_{x}=\outerproduct{X_{2}}{X_{1}}+\outerproduct{X_{1}}{X_{2}}$. This modification to the coupling constants results in a renormalization of the spectral density

[TABLE]

where $R=|r_{1}-r_{2}|$ is the dimer separation and we have assumed a linear dispersion $\omega(k)=c|k|$ with $c=1$. The cosine term arises from the phase factors in the original coupling terms.

We begin by studying a simplified case of a single spin described by Eq. (5) interacting with a bath with the unmodified spectral density, Eq. (4). It is straightforward to simulate this case with the reaction-coordinate (RC) mapping, which has been rigorously benchmarked for similar problems Iles-Smith et al. (2014); Iles-Smith (2015); Iles-Smith et al. (2016), allowing us to verify that TEMPO is able to accurately simulate dynamics in structured environments. This mapping takes a single collective mode of the environment (the RC) into the system definition; then the rest of the bath, which is assumed to be Markovian, couples to the now-augmented system. More details can be found in Ref. 111Supplementary material which cites Refs. Garg et al. (1985); Leggett (1984) and contains details of the RC mapping, numerical methods used and derivation of the expression for the bath displacement.. The RC mapping can be performed analytically for the spectral density in Eq. (4) Iles-Smith et al. (2014), and is expected to give accurate results in the cases we present here. To simulate more complex environments would require more RCs in the augmented system Hilbert space, whose dimension grows exponentially with the number of included RCs.

An alternative, numerically exact way of simulating the system dynamics arising from Eq. (5) is to use a Feynman sum-over-histories approach, using an influence functional to capture environmental memory effects. By writing this path integral representation as a sum over discrete time steps, we can calculate the system quantum dynamics by contracting a tensor network: this is the TEMPO approach Strathearn et al. (2018). The memory effects are then encoded in a matrix product state which contains information about the history of the system. This is propagated forward in time by successive contraction with a matrix product operator. By performing singular value decompositions and truncation after each timestep Orús (2014); Schollwöck (2011) we can capture bath memory times orders-of-magnitude beyond those possible with previous path integral approaches Makri and Makarov (1995a, b). The similarities between TEMPO and the process tensor Pollock et al. (2018) have recently been explored Jørgensen and Pollock (2019), and TEMPO has also been used to model dynamics in optomechanical systems Minoguchi et al. (2019). Further details on TEMPO can be found in Ref. Note (1).

In Fig. 2 we compare dynamics generated by these two approaches for the simplified model of the single spin interacting with the structured spectral density, for three different values of $\omega_{0}$. We find excellent agreement between the two algorithms for all sets of parameters. The RC approach is accurate for these parameters and this form of spectral density, and so able to capture the complex dimer dynamics. In this simple case RC requires significantly reduced numerical resources compared to TEMPO, though we find that both techniques take longer to converge at lower mode frequencies and higher temperatures. For the RC algorithm the reason for this is clear: the occupation of the mode included in the system increases in these regimes and so a larger system Hilbert space is needed for convergence. For TEMPO we find that the number of singular values that need to be retained for convergence increases – i.e. more system paths gain a significant amplitude, as expected when the environment becomes more occupied. These results show both that the RC approach is able to accurately simulate dynamics for this model Iles-Smith et al. (2016) and that TEMPO is able to deal with such structured spectral densities.

We now turn to analyzing the full model including the spatial modifications to the spectral density (i.e. using Eq. 6). Previous studies have shown that spatial correlations in the bath can have a significant effect on population transfer between localized systems Nalbach et al. (2010); Nazir (2009); McCutcheon and Nazir (2011); Strathearn et al. (2018). However, these studies focused on the case where the original spectral density $J_{0}(\omega)$ is much broader than that given in Eq. (2) – i.e. they do not focus on systems where just a narrow band of modes dominate the system interaction – and so the modification due to the spatial structure of the system is not as pronounced.

Since we are now using the general mapped spectral density in Eq. (6), treating the dynamics using the RC formalism becomes much more difficult: to accurately capture the dynamics it would be necessary to include more modes in the system than is computationally feasible. This means that for the remainder of this Letter all results are obtained using TEMPO.

In this model there are now three possible resonance conditions which can be met by matching two of: the bare system frequency $\Omega$, the characteristic frequency of the environment $\omega_{0}$, or the frequency corresponding to the separation between the TLSs $\omega_{R}=2\pi/R$. In Fig. 3 we show how choosing $\omega_{R}=\omega_{0}$ leads to a complete suppression of the main peak of the bare spectral density, and choosing $\omega_{R}=(1\pm 0.1)\omega_{0}$ leads to a very asymmetric lineshape. In the single line plots of Fig. 4 we show how these differences in spectral density manifest themselves in the system dynamics. We show results for the three TLS spacings $\omega_{R}=\omega_{0},(1\pm 0.1)\omega_{0}$ and for bare system frequencies $\Omega$ tuned to the same three possible values, resulting in nine sets of results.

The system dynamics are complex and difficult to interpret, with multiple oscillation frequencies appearing in the dynamics. These are not only due to the peaked structure of the original spectral density (as occurred for the simplified single spin model) but also because of the spatial correlations. However, we next show that TEMPO can be used not only to extract system dynamics, as discussed in Ref. Strathearn et al. (2018), but also to find evolution of bath degrees of freedom. This will enable us to gain more insight into how the interplay between the bath and system leads to the complex dynamics shown in Fig. 4. To do this we formally solve the Heisenberg equations of motion Gardiner and Collett (1985); Gardiner and Zoller (2004) for the expectation value of e.g. the annihilation operator for a particular bath mode, and obtain an expression entirely in terms of the expectation of the system operator which couples to the bath. For our model, Eq. (5), we have

[TABLE]

We thus have access to the dynamics of all the bath modes and can construct the real-space displacement of the bath by creating an appropriately weighted sum of the individual mode displacements. The quantity we consider is the displacement of a one-dimensional field Mahan (2013) given by:

[TABLE]

In the continuum limit this takes the following form:

[TABLE]

where we have used $\omega=|k|$. See Ref. Note (1) for a derivation of this expression.

We show the real-space displacement of the bath as a function of time and position as a color map beneath the corresponding system dynamics in Fig. 4. The black dotted lines indicate the positions of the two TLSs. In the fully resonant case (the centre plot) we see that most of the bath excitation gets trapped between the two TLSs which simply show oscillations which decay away at a very slow rate. We also see slow decay in the other two plots where $\Omega=\omega_{R}$; this is because $J(\Omega)=0$ in these cases and hence a simple Markovian treatment would predict no decay at all.

When at least one of the frequencies is detuned from this condition we see a more pronounced propagation of bath excitations away from the system as the trapping effect is reduced. In these cases the overall decay rate of the TLS is gradually enhanced as the detuning is increased, since now the value of $J(\Omega)$ is larger. The detuning can also introduce a beat frequency into the bath dynamics, most clearly seen in the bath population trapped between the two TLSs, and most pronounced in the most detuned cases, i.e. in the bottom right and top left panels of Fig. 4. This beating gives rise to distinctive revivals in the TLS dynamics.

At very short times all of the dynamics are very similar. On timescales $t<\omega_{R}^{-1}$ each TLS only senses its own local environment and hence behaves as if it were interacting with bosons described by the unmapped spectral density. However as soon as the influence of the other TLS is felt the dynamics become radically different to those predicted by an independent environment model. These dynamics are highly sensitive to the relative detunings of the three frequencies described above, becoming significantly different even for the very small detunings studied here.

In conclusion, we have shown that TEMPO can provide accurate simulations of quantum systems in structured environments. We first validated the technique by simulating the quantum dynamics of a simple model that can be solved accurately using the RC approach, and comparing the results. We then presented TEMPO simulations of the more complex dynamics that result from the interplay between a highly structured spectral function and spatial correlations between different parts of the system. TEMPO is not limited to specific forms for the spectral density and so provides a versatile method for simulating systems coupled to an environment with arbitrary structure. In addition, we have shown that it is possible to find the full dynamics of the environment directly from the system dynamics. This enabled us to explore how it is possible to tune the separation between the spins such that one can effectively remove the coupling to the dominant environmental mode and how the resultant dynamics are highly sensitive to slight variations in parameters around this point. Such exact environmental tracking can help to explain the behavior of open quantum systems in general, and so aid in the design of future quantum devices. This greater insight into the dynamics of open quantum systems could also lead to the development of new approximation schemes and to the self-consistent verification of others.

I Reaction Coordinate Mapping

In this section we outline how the Reaction Coordinate (RC) mapping is implemented to obtain a time-local master equation for the two-level system-RC (TLS-RC) reduced density matrix. The RC is the coordinate corresponding to the direction the bath equilibrium position is displaced when the system transitions between $\sigma_{z}$ eigenstates. It corresponds to a superposition of normal modes of the bath, and is defined by the following mapping:

[TABLE]

Here, $c^{\dagger}$ ($c$) create (destroy) an excitation in the reaction coordinate mode. The form of the transformation results in coupling between the RC and the residual modes, $b_{k}$. The spin-boson Hamiltonian is then mapped to the following form:

[TABLE]

where the last term exists to counter a renormalization of the RC’s potential that arises from the interaction term. Here the TLS-RC coupling is given by $\lambda^{2}=\sum_{k}g_{k}^{2}$ and this form ensures that the canonical bosonic commutation relation is satisfied. This coupling leads to an RC frequency of $\Omega_{RC}=\lambda^{-2}\sum_{k}\omega_{k}g_{k}^{2}$. The residual environment is characterised by a new spectral density $J_{RC}(\omega)=\sum_{k}|d_{k}|^{2}\delta(\omega-\omega_{k})$ which can be found by replacing the TLS in both the mapped and un-mapped Hamiltonians with a classical coordinate Garg et al. (1985). The spectral density does not contain any information about the system itself and therefore can be found, in both cases, through the classical equations of motion of this coordinate Leggett (1984); Iles-Smith et al. (2014). The under-damped spectral density considered in the main text leads to $\Omega_{RC}=\omega_{0}$, $\lambda=\sqrt{\pi\alpha_{\text{UD}}\omega_{0}/2}$ and $J_{RC}=\gamma\omega$ where $\gamma=\Gamma/2\pi\omega_{0}$.

Once this mapping has been performed we follow the procedure outlined in Ref. Iles-Smith et al. (2014) applying both Born and Markov approximations to arrive at the Schrödinger picture master equation for the reduced TLS-RC density operator, $\rho(t)$:

[TABLE]

where we have defined $\hat{C}=c^{\dagger}+c$ and $\hat{C}(t)$ is the same operator in the interaction picture at time $t$. Truncating the Hilbert space of the RC down to $n$ basis states, i.e. permitting a maximum of only $n$ excitations, allows us to numerically diagonalize $H_{S}$. This gives us a set of basis states $\ket{\phi_{j}}$ which satisfy $H_{S}\ket{\phi_{j}}=\phi_{j}\ket{\phi_{j}}$ allowing us to express the interaction picture operators as:

[TABLE]

where $C_{jk}=\bra{\phi_{j}}\hat{C}\ket{\phi_{k}}$ and $\omega_{jk}=\phi_{j}-\phi_{k}$. We can then re-write Eq. (S7) as:

[TABLE]

with rate operators:

[TABLE]

The resulting master equation can then be solved using standard open systems approaches for Markovian systems Breuer and Petruccione (2002).

II Time-Evolving Matrix Product Operators

In this section we outline the implementation of the Time-Evolving Matrix Product Operators (TEMPO) algorithm. For full details see Strathearn et al. (2018).

The quasi-adiabatic path integral (QUAPI) method Makri and Makarov (1995) is used to calculate the non-Markovian evolution of the system up to time $t_{N}=N\Delta t$ by summing over all possible paths the system has taken to that point from initial time $t_{0}$. In practice this can be done through an iterative tensor propagation routine where the object propagated is the augmented density tensor (ADT). The ADT is grown up to time $t_{N}$ from an initial physical density matrix $\rho_{i_{1}}(t_{1})$ through iterative application of tensors:

[TABLE]

The reduced system density matrix at time $t_{N}$ is then given by $\rho_{j_{N}}(t_{N})=\sum_{j_{1}\dots j_{N-1}}A^{j_{N}\dots j_{1}}(t_{N})$. For a harmonic bath linearly coupled to the system, the $B$ tensor is composed of influences across all pairs of time points from $t_{1}$ to $t_{N}$, as well as time-local components due to the coherent Hamiltonian evolution. The discretisation time-step $\Delta t$ needs to both resolve the features of the correlation function and ensure the Trotter errors are negligible. Without further considerations the exponential growth of the ADT only allows memory lengths of order $\sim 20\Delta t$ to be simulated. This restricts QUAPI to baths with relatively short-lived and/or slowly varying auto-correlation functions.

In TEMPO the ADT is built in the form of a matrix product state (MPS) Orús (2014); Schollwöck (2011) and a singular value decomposition (SVD) sweep is carried out at each time-step. The $B$ tensors are composed of a product of time-local operators and keeping these components separate leads to the natural formation of an MPS representation for the $A$ tensor given by:

[TABLE]

where the tensors at the two ends of the chain are rank two, and those in between are rank three. In this notation the superscripts are the ‘physical’ index which connects to the propagation tensor while the subscripts link adjoining $a$ tensors. At each step of the propagation an SVD is performed on each $a$ and singular values below a fixed precision $\chi$ (relative to the largest) are discarded. This leads to a reduced dimension of the internal indices $\alpha$ and a significant improvement on the exponential scaling present in QUAPI. The truncation acts to remove the least important internal degrees of freedom that are not needed for achieving converged dynamics. This proceedure has been shown to reduce the scaling of the required computational resources to polynomial with memory length for smooth spectral densities Strathearn et al. (2018).

In practice there are two convergence parameters associated with TEMPO which must be adjusted to produce the exact results in the main text:

• •

The size of the discretization time-step $\Delta t$, decreased until convergence.

• •

The singular value precision $\chi$, increased until convergence.

Further to these, with TEMPO a memory cut-off is also possible such that the growth of the ADT is stopped once it covers enough of the system’s history to capture all relevant memory effects. For the results presented in the main text, however, the long correlation times associated with the structured environments make this kind of memory cut-off impossible and as such none were used for generating any of the results. Table 1 gives the convergence parameters used for each of the results.

III Bath Dynamics

In this section we show that it is possible to calculate bath operator expectation values in terms of system dynamics and specifically present the derivation of Eq. (9) in the main text.

Consider general interaction and bath Hamiltonians given by:

[TABLE]

for a general system operator $\hat{A(t)}$ in the system interaction picture. The time evolution of the ladder operator $a_{k}(t)$ can be found from the Heisenberg equation of motion:

[TABLE]

Formally integrating this and taking the expectation with respect to the total density matrix $\rho=\rho_{S}\otimes\rho_{B}$ gives:

[TABLE]

Thus the dynamics of any bath mode can be found from knowing the exact system dynamics up to the time $t$. Calculations of higher order moments are much more involved requiring multiple numerical integrals to be performed and so for now we only consider quantities that can be constructed from the expressions above. The quantity analysed in the main text is given by:

[TABLE]

Now for the $\hat{A}=\sigma_{z}$ coupling from Eq. (5) in the main text we have $\tilde{g}_{-k}^{*}=\tilde{g}_{k}$ and $\omega_{k}=|k|=\omega_{-k}$ such that:

[TABLE]

where between the last two lines we have discarded the odd component of the integrand and since everything remaining is even we can then halve the integration domain.