Exploiting the causal tensor network structure of quantum processes to efficiently simulate non-Markovian path integrals

TL;DR

This paper introduces a tensor network-based algorithm leveraging the structure of the influence functional and process tensor to efficiently simulate non-Markovian quantum dynamics, enabling accurate computation of phonon emission spectra in strongly coupled systems.

Contribution

It establishes a novel connection between influence functionals and process tensors, leading to an efficient tensor network algorithm for non-Markovian quantum process simulation.

Findings

Achieved orders-of-magnitude speedup in simulations.

Accurately computed phonon emission spectra in the spin-boson model.

Demonstrated significant deviations from memoryless assumptions.

Abstract

In the path integral formulation of the evolution of an open quantum system coupled to a Gaussian, non-interacting environment, the dynamical contribution of the latter is encoded in an object called the influence functional. Here, we relate the influence functional to the process tensor -- a more general representation of a quantum stochastic process -- describing the evolution. We then use this connection to motivate a tensor network algorithm for the simulation of multi-time correlations in open systems, building on recent work where the influence functional is represented in terms of time evolving matrix product operators. By exploiting the symmetries of the influence functional, we are able to use our algorithm to achieve orders-of-magnitude improvement in the efficiency of the resulting numerical simulation. Our improved algorithm is then applied to compute exact phonon emission spectra for the spin-boson model with strong coupling, demonstrating a significant divergence from spectra derived under commonly used assumptions of memorylessness.

Equations

$${\rho }_k(\{{\mathcal{A}}_j\})={\text{tr}}_E\left\{{\mathcal{U}}_{\delta t}{\mathcal{A}}_{k-1}...{\mathcal{U}}_{\delta t}{\mathcal{A}}_0\left[{\chi }_0\right]\right\},$$

$${{\rm Y}}_{k:0}=$$

$$\times \text{outerproduct}{s}_k^{\prime }{s}_{k-1}\dots s_1^{\prime }{s}_0{s}_0^{\prime }{r}_k^{\prime }{r}_{k-1}\dots r_1^{\prime }{r}_0{r}_0^{\prime },$$

$${\widetilde{{\rm Y}}}_{k:0}={\left({\mathcal{V}}_{\delta t}^{1/2}\otimes {{\mathcal{V}}_{\delta t}^{\ast }}^{1/2}\right)}^{\otimes k}\left[{\mathcal{F}}_{k:0}\right]\otimes {\rho }_0,$$

$${\mathcal{F}}_{k:0}=$$

$$\times \text{outerproduct}{s}_k{s}_k\dots s_1{s}_1{r}_k{r}_k\dots r_1{r}_1.$$

$${\mathcal{F}}_{k:0}^{{\alpha }_k...{\alpha }_1}=\prod _{i=1}^k\prod _{j=1}^i{\left[b_{(i-j)}\right]}^{{\alpha }_i{\alpha }_j},$$

$${\mathcal{F}}_{k:0}^{{\alpha }_k...{\alpha }_1}$$

$$\displaystyle=\prod_{i=1}^{k}\ \ \begin{minipage}[h]{0.0pt} \vspace{0pt} \includegraphics[width=113.81102pt]{tempo_propagator.png} \end{minipage}\qquad\qquad\qquad\qquad\qquad\quad,$$

$${\mathcal{F}}_{k:0}^{{\alpha }_k...{\alpha }_1}={\mathcal{G}}_{{\beta }_{k-1}...{\beta }_1}^{{\alpha }_k{\alpha }_{k-1}...{\alpha }_1}{\mathcal{F}}_{k-1:0}^{{\beta }_{k-1}...{\beta }_1},$$

$${\mathcal{F}}_{k:0}^{{\alpha }_k...{\alpha }_1}$$

$$\displaystyle=\prod_{i=1}^{k}\ \ \begin{minipage}[h]{0.0pt} \vspace{0pt} \includegraphics[width=113.81102pt]{process_propagator.png} \end{minipage}\qquad\qquad\qquad\qquad\qquad\quad,$$

$${\tilde{\mathcal{F}}}_{(k:i+1)}^{{\alpha }_k...{\alpha }_1}:={\mathcal{C}}_{{\beta }_k...{\beta }_{i+1}}^{{\alpha }_k...{\alpha }_{i+1}}{\tilde{\mathcal{F}}}_{(k:i)}^{{\beta }_k...{\beta }_{i+1}{\alpha }_i...{\alpha }_1},$$

$${\rho }_k(\{{\mathcal{A}}_j{\}}_{j=0}^{k-1})={\text{tr}}_E\left\{{\mathcal{U}}_{k:k-1}{\mathcal{A}}_{k-1}...{\mathcal{U}}_{1:0}{\mathcal{A}}_0\left[{\chi }_0\right]\right\}$$

$$U_{j:j-1}=T_{\leftarrow }\exp \left\{-i{\int }_{t_{j-1}}^{t_j}dsH(s)\right\}.$$

$$\mathsf{A}=\mathcal{A}\otimes \mathcal{I}[\Psi ]=\sum _{sr{s}^{\prime }{r}^{\prime }}{\mathcal{A}}^{(s,r,s^{\prime },r^{\prime })}\text{outerproduct}s{s}^{\prime }r{r}^{\prime },$$

$${\rho }_k(\{{\mathcal{A}}_j{\}}_{j=0}^{k-1})={\text{tr}}_{k-1:0}\left\{{{\rm Y}}_{k:0}(\text{mathbbm}{1}_k\otimes {\mathbf{A}}_{k-1:0}^T)\right\},$$

$${{\rm Y}}_{k:0}=$$

$$\times \text{outerproduct}{s}_k^{\prime }{s}_{k-1}\dots s_1^{\prime }{s}_0{s}_0^{\prime }{r}_k^{\prime }{r}_{k-1}\dots r_1^{\prime }{r}_0{r}_0^{\prime },$$

$${\mathcal{U}}_{j:j-1}^{(s^{\prime },r^{\prime },s,r)}\simeq$$

$$=$$

$$\times ⟨s|{\mathcal{V}}_{j:j-1}^{\ast 1/2}\left[\text{outerproduct}tu\right]|r⟩{\mathcal{W}}_{j:j-1}^{(t^{\prime },u^{\prime },t,u)},$$

$${\widetilde{{\rm Y}}}_{k:0}=\underset{j=1}{\overset{k}{⨂}}\left({\mathcal{V}}_{j:j-1}^{1/2}\otimes {\mathcal{V}}_{j:j-1}^{\ast 1/2}\right)\left[{\mathcal{F}}_{k:0}\right]\otimes {\rho }_0,$$

$${\mathcal{F}}_{k:0}=$$

$$\times \text{outerproduct}{s}_k^{\prime }{s}_{k-1}\dots s_1{s}_1^{\prime }{s}_0{r}_k^{\prime }{r}_{k-1}\dots r_1{r}_1^{\prime }{r}_0.$$

$${\mathcal{W}}_{j:j-1}^{(s,r)}[{\rho }^E]=W_{j:j-1}^{(s)}{\rho }^E{W}_{j:j-1}^{(r)†}$$

$$W_{j:j-1}^{(s)}=T_{\leftarrow }\exp [-i{\int }_{t_{j-1}}^{t_j}dx{B}_s(x)].$$

$${\mathcal{F}}_{k:0}^{{\alpha }_k...{\alpha }_1}=$$

$$=$$

$$-{\zeta }_{i,j}^{(r_i,s_j)}-{\zeta }_{i,j}^{(s_i,r_j)\ast })],$$

$${\zeta }_{i,j}^{(u,v)}={\int }_{t_{i-1}}^{t_i}dx{\int }_{t_{j-1}}^{t_j}dy\text{tr}\left\{{\hat{B}}_u(x){\hat{B}}_v(y)\tau \right\}$$

$${\zeta }_{i,i}^{(u,v)}={\int }_{t_{i-1}}^{t_i}dx{\int }_{t_{i-1}}^{x}dy\text{tr}\left\{{\hat{B}}_u(x){\hat{B}}_v(y)\tau \right\}.$$

$${\eta }_{i-j}=\{\begin{array}{ll}{\int }_{t_{i-1}}^{t_i}{\int }_{t_{j-1}}^{t_j}d{t}^{\prime }d{t}^{\prime \prime }C(t^{\prime }-t^{\prime \prime }),& i\ne j\\ {\int }_{t_{i-1}}^{t_i}{\int }_{t_{i-1}}^{t^{\prime }}d{t}^{\prime }d{t}^{\prime \prime }C(t^{\prime }-t^{\prime \prime }),& i=j\end{array},$$

$$C(t)=\frac{1}{\pi }{\int }_0^{\infty }d\omega J(\omega )\frac{\cosh \left[\omega \left(\beta /2-it\right)\right]}{\sinh \left[\beta \omega /2\right]}.$$

$${\left[b_{(i-j)}\right]}^{{\alpha }_i{\alpha }_j}=e^{-({\lambda }_{s_i}-{\lambda }_{r_i})({\eta }_{i-j}{\lambda }_{s_j}-{\eta }_{i-j}^{\ast }{\lambda }_{r_j})}.$$

Full text

Introduction. – All nanoscale quantum systems are open, meaning they inevitably interact with their environments, exchanging energy and generating correlations. If the system and its environment remain approximately uncorrelated, then the reduced system dynamics is well described by a Markovian model Breuer and Petruccione (2002); Weiss (2012); Carmichael (2003). However, in physical systems such as photosynthetic complexes, nanoscale lasers and quantum thermal machines Mujica-Martinez et al. (2013); McCutcheon (2016); Newman et al. (2017), the need to go beyond a Markovian description has long been recognized, and techniques accounting for non-Markovian physics have been developed, with greater or lesser breadth of applicability. Analytical methods involving time-local equations of motion exist, but tend to be highly restricted to specific parameter regimes Vacchini and Breuer (2010); Lee et al. (2012); Fruchtman et al. (2016). Exact simulation often requires numerical methods, e.g. discrete path integrals Makri and Makarov (1995a, b); Nalbach et al. (2011); Dattani (2012); Strathearn et al. (2017), time non-local memory kernels Shi and Geva (2003); Cohen and Rabani (2011); Cerrillo and Cao (2014); Rosenbach et al. (2016); Buser et al. (2017); Gelzinis et al. (2017); Pollock and Modi (2018), hierarchical equations of motion Tanimura (2006); Strümpfer and Schulten (2012) and others Bulla et al. (2008); Chen et al. (2017a, b); Cerrillo and Cao (2014). Overall these methods tend to scale unfavourably with both the simulation time and the system size de Vega and Alonso (2017), making them inapplicable to important processes involving large complexes or when long time dynamics is important.

Recently, tensor network methods have been applied to the simulation Chin et al. (2010); Prior et al. (2010); Schröder and Chin (2016); Wall et al. (2016) and characterization Pollock et al. (2018a); Luchnikov et al. (2019) of open quantum dynamics. Physically, these methods incorporate the fact that typical open quantum systems are only finitely correlated with their environments, massively reducing their description Orús (2014). In particular, Strathearn et al. Strathearn et al. (2018) reformulated the discrete path integral for open systems with Gaussian environments in terms of matrix product operators; the resulting time-evolving matrix product operator (TEMPO) algorithm is numerically exact and has an efficiency comparable to other state of the art methods. By effectively only considering the most important non-Markovian contributions to the dynamics, the algorithm circumvents the exponential memory scaling of the bare path integral representation, in a similar spirit to earlier path-filtering techniques Sim and Makri (1997); Sim (2001); Dattani (2013). Motivated by this success, it is natural to ask if tensor network methods can be efficiently generalized from the simulation of reduced system density operators, to the simulation of general non-Markovian processes and multi-time correlations, which typically require many realizations of the dynamics to characterize.

In this Letter, we propose such a generalization, by making a formal connection between the path integral structure and the recently developed process tensor framework for characterizing general non-Markovian quantum processes Pollock et al. (2018a). We then use this to argue for an alternative formulation of the TEMPO algorithm, where we exploit the symmetry of the underlying tensor network to better account for the causal structure inherent in the dynamics. This not only allows us to efficiently compute multi-time correlation functions – the simulation need only run once to extract all multi-time observable properties – but also opens the door to the simulation of more general models. Our alternative formulation is demonstrated to significantly improve the efficiency of the method, which we use to straightforwardly compute non-Markovian emission spectra for the spin-boson model, beyond the point where the commonly used quantum regression theorem breaks down Li et al. (2018).

Process tensor framework. – We consider stationary unitary dynamics of the system $S$ we are interested in along with its environment $E$, and suppose the system is transformed by superoperators $\mathcal{A}_{j}$ at the discrete times $\left\{t_{k-1},...,t_{0}\right\}$, which we take to be at evenly spaced intervals of $\delta t=t_{j}-t_{j-1}$. In an experiment, these superoperators could correspond to actual interventions on the system as it evolves, i.e. unitary rotations, measurements etc., in which case they are completely-positive and, if the interventions are not conditional on a particular measurement outcome, trace preserving. Otherwise, the set $\{\mathcal{A}_{j}\}$ could represent more abstract transformations useful in the computation of physical quantities such as operator expectation values or emission spectra. The reduced, and potentially subnormalized, state of the system at time $t_{k}$ is given by

[TABLE]

where $\mathcal{U}_{\delta t}$ is a superoperator representation of the unitary evolution of duration $\delta t$, i.e. $\mathcal{U}_{\delta t}[\rho]=U_{\delta t}\rho U_{\delta t}^{\dagger}$, with $U_{\delta t}$ a unitary matrix, and $\chi_{0}$ is the initial system-environment state. The inclusion of intermediate transformations makes it possible for us to consider a much broader class of physical properties than free evolution of the density operator (corresponding to $\mathcal{A}_{j}=\mathcal{I}$ the identity superoperator $\forall j$) would allow.

Since the state at time $t_{k}$ in Eq. (1) is linearly related to each of the set of superoperators $\{\mathcal{A}_{j}\}$ it can be written as a linear function of the tensor product of their Choi state representations $\textbf{A}_{k-1:0}=\mathsf{A}_{k-1}\otimes\dots\otimes\mathsf{A}_{1}\otimes\mathsf{A}_{0}$, with $\mathsf{A}_{j}:=\sum_{sr}\mathcal{A}_{j}[\outerproduct{s}{r}]\otimes\outerproduct{s}{r}$ obtained via the Choi-Jamiołkowski isomorphism Milz et al. (2017a); Nielsen and Chuang (2011); Wilde (2017); here, $\{\ket{s}\}$ forms an orthonormal basis for $S$. Specifically, $\rho_{k}(\{\mathcal{A}_{j}\})={\rm tr}_{k-1:0}\{\Upsilon_{k:0}(\mathbbm{1}_{k}\otimes\textbf{A}_{k-1:0}^{T})\}$, with the trace over all subsystems on which $\textbf{A}_{k-1:0}$ acts. As we detail explicitly in Appendix A,

[TABLE]

with environment superoperators $\mathcal{U}_{\delta t}^{(s^{\prime},r^{\prime},s,r)}[\rho^{E}]=\bra{s^{\prime}}U_{\delta t}(\outerproduct{s}{r}\otimes\rho^{E})U^{\dagger}_{\delta t}\ket{r^{\prime}}$ and operators $\chi^{(r^{\prime},s^{\prime})}_{0}=\bra{r^{\prime}}\chi_{0}\ket{s^{\prime}}$, is the Choi representation of the process tensor Pollock et al. (2018a), a many-body operator (on $2k+1$ copies of $S$) containing all information about the system’s evolution that is independent of the transformations $\{\mathcal{A}_{j}\}$. Correlations between subsystems of $\Upsilon_{k:0}$ correspond to temporal correlations between observables, and a representation in terms of process tensors has been shown to consistently generalize stochastic processes, and related notions such as the Markov property and Markov order, to the quantum case Milz et al. (2017b); Pollock et al. (2018b); Taranto et al. (2019a, b). The process tensor is illustrated graphically in Fig. 1, and can be thought of simply as a sequence of correlated maps on the system Sakuldee et al. (2018). Unlike in a conventional open quantum systems picture, where density operators are mapped to density operators, this operational formulation stresses that the proper input to a quantum process is the set of interventions $\textbf{A}_{k-1:0}$, and that the intermediate dynamics, and the initial state, are features of the process itself.

Gaussian influence functional. – Here, we consider the specific structure of the process tensor for systems interacting with Gaussian environments, where the system-environment Hamiltonian and initial state depend at most quadratically on environment creation and annihilation operators. For concreteness, we focus on spin-boson type models, but our results would extend to fermionic baths as well Atland and Simon (2010). Working in natural units ($\hbar=k_{B}=1$), we consider a spin system, with Hilbert space dimension $d$, interacting linearly with a bath of harmonic oscillators described by the Hamiltonian $H=H_{0}+H_{B}$. Here, $H_{0}$ describes the free spin system and the full bath influence is collected in $H_{B}=\hat{s}\sum_{n}\left(g_{n}\hat{a}_{n}+g_{n}^{*}\hat{a}_{n}^{\dagger}\right)+\sum_{n}\omega_{n}\hat{a}_{n}^{\dagger}\hat{a}_{n}$. A bath mode $n$ has energy $\omega_{n}$, and is created (annihilated) by the bosonic operator $\hat{a}_{n}^{\dagger}$ ($\hat{a}_{n}$). The system operator $\hat{s}$ interacts with the bath with coupling constants $g_{n}$. Additional linear interaction terms to different system operators could be included, as long as all these system operators commute. For simplicity, we take the initial state to be product, such that $\chi_{0}=\rho_{0}\otimes\tau_{\beta}$, with the environment initially described by a thermal state $\tau_{\beta}=\exp\left[-\beta\sum\omega_{n}a^{\dagger}_{n}a_{n}\right]/\mathcal{Z}$ at inverse temperature $\beta$, where $\mathcal{Z}=\tr\{\exp\left[-\beta\sum\omega_{n}a^{\dagger}_{n}a_{n}\right]\}$. This choice is not essential, and other, possibly correlated, Gaussian initial states of the environment could be considered.

In the limit that the time difference $\delta t$ is small, the generated unitary dynamics can be approximately separated into contributions arising from $H_{0}$ and $H_{B}$ as $\mathcal{U}_{\delta t}\simeq\mathcal{V}^{1/2}_{\delta t}\mathcal{W}_{\delta t}\mathcal{V}^{1/2}_{\delta t}$, where $\mathcal{V}_{\delta t}$ describes the free dynamics of $S$ and $\mathcal{W}_{\delta t}$ describes the environment influence. The discrepancy between the approximate unitary maps and the actual ones vanishes as $\mathcal{O}(\delta t^{3})$ for this symmetric decomposition Trotter (1959). Note that the approximate model does not break unitarity, and so corresponds to a valid physical process independently of step size. Moreover, since the Hamiltonian only contains a single interaction term, the interaction unitary preserves the eigenbasis of the corresponding system operator $\hat{s}=\sum_{s}\lambda_{s}\outerproduct{s}{s}$: $\bra{s^{\prime}}\mathcal{W}_{\delta t}[\outerproduct{s}{r}]\ket{r^{\prime}}=\delta_{ss^{\prime}}\delta_{rr^{\prime}}\mathcal{W}_{\delta t}^{(s,r)}$. Together with the decomposition of unitary maps, this allows us to write the approximate process tensor Choi state as

[TABLE]

where $\mathcal{F}_{k:0}$ is an operator representation of the discretized Feynman-Vernon influence functional Feynman and Vernon (1963) encoding environment induced correlations

[TABLE]

For Gaussian environments, the bath degrees of freedom can be traced over analytically using standard path integral techniques Makri and Makarov (1995a, b); Diósi and Ferialdi (2014); Strathearn et al. (2017). In this case, introducing the $d^{2}$ compound indices $\alpha=(s,r)$, an element of the influence functional $\mathcal{F}_{k:0}^{\alpha_{k}...\alpha_{1}}:=\bra{s_{k}s_{k}\dots s_{1}s_{1}}\mathcal{F}_{k:0}\ket{r_{k}r_{k}\dots r_{1}r_{1}}$ can be decomposed as

[TABLE]

where $b_{(i-j)}$ is called an influence tensor; the exact form, which can often be approximated by an analytic function Dattani et al. (2012), is given in Appendices B and C along with a full derivation of Eqs. (4) and (5). The influence tensors connect the dynamics around time step $i$ with that around step $j$, quantifying the temporal correlations mediated by the environment between those two points; that is, they describe memory effects. Since the Hamiltonian is time-independent, the individual tensors $\left[b_{l}\right]$ depend only on the temporal separation $l\delta t$, simplifying the potential complexity considerably. However, since the influence functional is a $k$ index tensor, it is still potentially exponentially complex; we now show how viewing Eq. (5) as a tensor network can make its calculation more tractable.

Tensor network simulation. – In many cases, the environment interaction produces only finite length correlations in $\mathcal{F}_{k:0}$, a fact used by the authors of Ref. Strathearn et al. (2018) to circumvent the exponential complexity growth by representing it efficiently in terms of matrix product operators (MPOs) Schollwöck (2011). To introduce this representation we first extend our two-index influence tensors into three-index tensors as $\left[b_{(i-j)}\right]^{\gamma\alpha_{i}}{}_{\alpha_{j}}:=\delta^{\gamma}_{\alpha_{j}}\left[b_{(i-j)}\right]^{\alpha_{i}}{}_{\alpha_{j}}$, where by convention an upper and a lower repeated index in a product of tensors is summed over (otherwise, tensor elements differing only through raising or lowering are treated as equal). In terms of these, we then define the non-local time-evolving MPOs

[TABLE]

where the outgoing (ingoing) arrows in the graphical representation indicate upper (lower) indices, and lines running through a given row or column are fixed to have the same index through Kronecker deltas. At the right boundary of the tensors shown in the second line of Eq. (6), we end up with a redundant lower index which we can trace over, while at the left boundary we impose that the two upper indices be equal. The full influence functional can then be constructed by iteratively multiplying such MPOs. If we label the individual MPOs in the product by $G^{\alpha_{i}...\alpha_{1}}$, then we can express the iterative multiplication as

[TABLE]

with $\mathcal{G}^{\alpha_{k}\alpha_{k-1}...\alpha_{1}}_{\ \ \ \ \beta_{k-1}...\beta_{1}}:=G^{\alpha_{k}...\alpha_{1}}\delta^{\alpha_{k-1}}_{\beta_{k-1}}...\delta^{\alpha_{1}}_{\beta_{1}}$; this is represented graphically by the two-dimensional tensor network shown in Fig. 2a. Conceptually the use of time-evolving MPOs, allows the state of the system to be propagated by updating indices to encode memory effects from the past process. This type of propagation is analogous to a description in terms of a time non-local memory kernel, since the open legs are connected to tensors describing the influence of the state at various points in its history 111This correspondence is not precise, however, and our usage of ‘(non-)local’ should not be confused with that in the context of memory kernel convolution..

A key insight of this paper is that the decomposition of the influence functional into MPOs is not unique. Kronecker deltas implicit in Eqs. (6) and (7) mean that the open leg in a given row or column in Fig. 2 could be shifted to any tensor in that same row or column. In particular, the causal structure of the process tensor motivates an alternative definition in terms of local time-evolving MPOs

[TABLE]

where now we end up with a redundant index at the left boundary, which we trace over, and at the right boundary we impose the condition that the lower index must equal $\alpha_{i}$. As with the non-local propagators, the full influence functional is constructed iteratively by locally contracting MPOs. Labelling the individual MPOs in the product by $C^{\alpha_{k}...\alpha_{i}}$, the iterative multiplication can be expressed as (for $i\geq 1$)

[TABLE]

with $\mathcal{C}^{\alpha_{k}...\alpha_{i}}_{\gamma_{k}...\gamma_{i}}:=C^{\alpha_{k}...\alpha_{i}}\delta^{\alpha_{k}}_{\gamma_{k}}...\delta^{\alpha_{i}}_{\gamma_{i}}$ and $\tilde{\mathcal{F}}^{\alpha_{k}...\alpha_{1}}_{(k:1)}:=C^{\alpha_{k}...\alpha_{1}}$. The resulting network representation for the influence functional $\mathcal{F}_{k:0}=\tilde{\mathcal{F}}_{(k:k)}$ is shown in Fig. 2c. Conceptually, the local time-evolving MPOs propagate the state by updating a set of effective memory space indices. These indices describe how the environment is conditioned by the process at a given time, and this information on the conditioning can be propagated locally. Since the process tensor, and hence the influence functional, has a well-defined causal structure, this conditioning only occurs from the past to the future. This means that, for a fixed evolution time, the size of the tensor to be updated decreases with each iteration.

In contracting the network, efficiency is achieved by incorporating a tensor compression procedure of the obtained boundary in each iteration. In this work we make use of the singular value compression procedure (see Appendix D) Schollwöck (2011); Strathearn et al. (2018). Roughly speaking, the local tensors of the boundary are subjected to a singular value decomposition. The compression consists of discarding the eigenvalues below a specified singular value cutoff $\lambda_{c}$, quantifying the hardness of the compression. For the non-local algorithm (Fig. 2a,b), the tensors contracted in each iteration encode information about the influence of multiple time-steps on each other. When correlations become smaller at longer time scales, as is typically the case, not all this information is relevant for describing the process as a whole. The local algorithm (Fig. 2c,d) incorporates this insight, and separates out the most important contribution by including only the future influence of the environment at each timestep. Generally the most local contributions have the largest singular values, and therefore the separation means that the part of the boundary being propagated in the local case is less correlated, which translates into a more efficient algorithm.

Network complexity for a two level system. – We now turn to the specific simulation of the dynamics of a two-level system, and compare the performance of the non-local and local algorithms. Consider the free Hamiltonian $H_{0}=\Omega\sigma_{x}/2$ and $\hat{s}=\sigma_{z}/2$, where $\sigma_{x}$ and $\sigma_{z}$ are the usual Pauli operators. The environment is fully characterized by its spectral density defined as $J(\omega)=\sum_{n}|g_{n}|^{2}\delta\left(\omega-\omega_{n}\right)$ Breuer and Petruccione (2002). Here we consider a continuum bath model with the spectral density $J(\omega)=(\alpha\omega_{c}/2)(\omega/\omega_{c})^{\nu}\exp\left(-\omega/\omega_{c}\right)$ with coupling strength $\alpha$, cutoff frequency $\omega_{c}$ and Ohmicity $\nu$, where for an Ohmic spectral density $\nu=1$.

The computational complexity is quantified by the computation time required to contract a network of a certain size with a fixed singular value cutoff (see Appendix D for details on implementation). In general, this will depend on the overall magnitude of the influence functional, as well as the characteristic memory time quantifying how the elements of the influence tensors $b_{(i-j)}$ decrease in magnitude at large $|i-j|$. In Appendix E, we show that, for a fixed evolution time, the memory time goes as $\alpha/(\beta\omega_{c})$ when $\omega_{c}$ is large, and that the overall coupling goes as $\alpha\omega_{c}t_{\rm max}^{2}/\beta$ when $\omega_{c}$ is small. In Fig. 3a, we plot the computation time for the local and non-local algorithms as a function of coupling strength. We find that the local representation outperforms the non-local one by one-to-two orders of magnitude, even at weak coupling, and that the improvement increases at larger coupling strengths (in Appendix F we show that there is an advantage everywhere across a wide range of parameters). Furthermore, in Fig. 3b, we illustrate that this computational efficiency is maintained as the characteristic timescale of the bath (the inverse of the environment cutoff frequency) is varied, consistent with our predictions.

It should be kept in mind that, unlike in most other quantum simulation methods, including the original TEMPO algorithm, the object we are computing is the full process tensor, which encodes all multi-time correlations, and from which a host of properties can be extracted efficiently. In particular, we can compute the steady state emission spectrum $S(\Delta\omega)=\text{Re}\left[\int_{0}^{\infty}d\tau(g^{(1)}(\tau)-g^{(1)}(\infty))e^{-i\Delta\omega\tau}\right]$, defined in terms of the two-point correlation function $g^{(1)}(\tau)=\lim_{t\rightarrow\infty}\left\langle\sigma^{\dagger}(t+\tau)\sigma(t)\right\rangle$. The two-point correlation function is defined in terms of the raising and lowering operators on the spin system. To compute it, we take all superoperators which the process tensor acts on to be the identity superoperator $\mathcal{I}$ (with action $\mathcal{I}[\rho]=\rho$) except for two, which append a raising or lowering operator respectively. In Fig. 3c, we study the physical effects of system-environment correlations by looking at the phonon emission spectrum. We compare this with the spectrum computed using the quantum regression theorem, which approximates intermediate dynamics with that from an initial product state and is valid in the weak-coupling limit McCutcheon (2016). The regression theorem correlations are obtained by breaking all correlations in the full process tensor Choi state across time steps at which the raising and lowering operators are evaluated. In addition, we compare the exact spectrum with a fully Markovian process in which correlations are broken after each time step. Fig. 3c shows that non-Markovian effects produce a phonon sideband in the spectrum at positive frequencies, this is contrasted with the regression theorem result which gives a symmetric sideband structure. Furthermore the fully Markovian spectrum contains no phonon sidebands, but rather a resonant emission peak. These significant differences illustrates the importance of accounting for non-Markovian physics.

Conclusion. – In this Letter, we have established a direct connection between recent frameworks for characterizing general non-Markovian quantum processes and the path integral formulation of open quantum dynamics. By relating the influence functional to a process tensor on an infinitesimal time grid, which has an explicit causal structure, we were able to build on recent progress in the classical simulation of Gaussian open quantum systems in terms of tensor networks. Specifically, we showed that the speed of the TEMPO algorithm, when computing the multi-time properties encapsulated in the process tensor, can be improved by orders-of-magnitude by shifting the boundary of the corresponding tensor network from a temporally non-local to a local one. Our contribution is immediately applicable to the efficient simulation of realistic complex open systems, and additionally illustrates the utility of thinking about a time-local propagation of a conditioned environment space, rather than a description resembling the use of a non-local memory kernel.

The utility of the non-local TEMPO algorithm has been amply illustrated by computing the Ohmic localization transition, and the dynamics of complex problems with multiple separated timescales Strathearn et al. (2018). The improved algorithm presented here is capable of exploring the same physics more efficiently, and, in addition, easily extends to the computation of multi-time observables, of the sort crucial to describing, for example, multi-dimensional spectroscopy experiments Yuen-Zhou et al. (2014). Moreover, relating path integral techniques to the more general process tensor formalism indicates how they might be generalized to more complex system-environment interactions, or even beyond the Gaussian regime. Even within the spin-boson model, the freedom of boundary choice in Fig. 2 that we have identified could be further exploited in other contexts. While the local choice appears optimal here, it may be that for other, structured spectral densities, where there are recurrent correlations, different boundary choices are more efficient, a point whose exploration we leave for future work.

Appendix A Process tensor formalism.

As described in the main text, we consider the scenario characterized by Eq. (1), where an open quantum system $S$ is periodically interrogated as it evolves. The slightly more general case involves an arbitrary time separation between interventions and a possibly time-dependent $SE$ Hamiltonian $H(t)$, such that the state at the $k^{\rm th}$ time is given by

[TABLE]

where $\mathcal{U}_{j:j-1}$ is the time-evolution superoperator from time $t_{j-1}$ to $t_{j}$ with action $\mathcal{U}_{j:j-1}[\rho]=U_{j:j-1}\rho U_{j:j-1}^{\dagger}$ in terms of the time-ordered exponential

[TABLE]

As in the main text, $\chi_{0}$ is the potentially correlated initial $SE$ state and each $\mathcal{A}_{j}$ can be any superoperator, though only those that are completely positive correspond to physically realizable transformations Milz et al. (2017a).

Any superoperator can be represented in terms of an operator sum as $\mathcal{A}[\rho]=\sum_{n}X_{n}\rho Y_{n}^{\dagger}$ (completely positive maps are characterized by $X_{n}=Y_{n}$, in which case the latter are called Kraus operators). As such, by taking the trace over the final state in Eq. (10), any multi-time correlation function $\langle B_{0}(t_{0})\dots B_{k-1}(t_{k-1})A_{k-1}(t_{k-1})\dots A_{0}(t_{0})\rangle_{\chi_{0}}$, where $\{A_{j}(t_{j})=U_{j:0}^{\dagger}A_{j}U_{j:0}\}$ and $\{B_{j}=U_{j:0}^{\dagger}B_{j}U_{j:0}\}$ are Heisenberg picture operators on $S$, can be obtained by choosing $\mathcal{A}_{j}[\rho]=A_{j}\rho B_{j}$. In this way, with a sufficiently dense set of times $\{t_{j}\}$, any dynamically observable property of $S$ can be represented in the form of Eq. (10) (correlation functions involving fewer observables can be obtained by choosing some of the $\{A_{j}\}$ and $\{B_{j}\}$ to be the identity operator). In addition, the freely evolved final state $\rho_{k}$ can be obtained by choosing all $\mathcal{A}_{j}=\mathcal{I}$, where the latter is the identity superoperator with action $\mathcal{I}[\rho]=\rho$.

As we will now show, all this information can be encoded in a single object, the process tensor. While it is often introduced in terms of an abstract multi-linear map, we will here express it solely in terms of a concrete matrix representation via a version of the Choi-Jamiołkowski isomorphism. The Choi state (or Choi matrix) $\mathsf{A}$ of a superoperator $\mathcal{A}$ is defined in terms of its action on one half of the (unnormalized) maximally entangled state $\Psi=\sum_{sr}\outerproduct{ss}{rr}$, where $\{\ket{s}\}$ forms an orthonormal basis for the $d$-dimensional system, as

[TABLE]

with $\mathcal{A}^{(s,r,s^{\prime},r^{\prime})}=\bra{s}\mathcal{A}[\outerproduct{s^{\prime}}{r^{\prime}}]\ket{r}$. Labelling the first and second copies of the original system’s Hilbert space $\mathtt{o}$ (for output) and $\mathtt{i}$ (for input) respectively, the action of the superoperator on an initial state $\rho$ can be written in terms of this representation as $\mathcal{A}[\rho]=\tr_{\mathtt{i}}[\mathsf{A}\,\mathbbm{1}_{\mathsf{o}}\otimes\rho^{T}]$, where the trace is over the input subsystem.

By expanding out the action of the superoperators in Eq. (10) and inserting a resolution of the identity on $S$ to the left and right of every unitary matrix, one arrives at the equivalent expression

[TABLE]

with $\mathbf{A}_{k-1:0}=\mathsf{A}_{k-1}\otimes\cdots\otimes\mathsf{A}_{1}\otimes\mathsf{A}_{0}$ and the trace over all (input and output) Hilbert spaces on which $\mathbf{A}_{k-1:0}$ acts. The positive operator $\Upsilon_{k:0}$ is (the Choi state of) the process tensor, and it can be expressed in terms of the underlying $SE$ dynamics as

[TABLE]

with $\mathcal{U}_{j:j-1}^{(s^{\prime},r^{\prime},s,r)}[\rho^{E}]=\bra{s^{\prime}}U_{j:j-1}(\outerproduct{s}{r}\otimes\rho^{E})U^{\dagger}_{j:j-1}\ket{r^{\prime}}$ and $\chi^{(r^{\prime},s^{\prime})}_{0}=\bra{r^{\prime}}\chi_{0}\ket{s^{\prime}}$. This object can be directly constructed by swapping the system with one half of a maximally entangled state at each point where a superoperator $\mathcal{A}$ is to be applied Pollock et al. (2018a). When the $SE$ Hamiltonian is time independent and $t_{j}-t_{j-1}=\delta t$ for all $j$, Eq. (A) is equivalent to Eq. (Exploiting the Causal Tensor Network Structure of Quantum Processes to Efficiently Simulate Non-Markovian Path Integrals) of the main text.

A representation in terms of the process tensor separates the process and external interventions, as illustrated in Fig. 1 of the main text. The process tensor’s properties reflect the necessary features of any physical open dynamics: $\Upsilon_{k:0}$ is positive if and only if the process is completely positive, and causality is encoded in the hierarchy of trace conditions $\tr_{j}\Upsilon_{j:0}=\Upsilon_{j-1:0}\otimes\mathbbm{1}_{\mathtt{o}_{j-1}}$. While we have expressed it in terms of $SE$ quantities, it is an operator only on copies of the Hilbert space of $S$ and it has a natural matrix product form, allowing for an efficient representation in many cases, a fact we exploit in this paper.

Appendix B Connection with influence functional

In the case where the time spacing is small (and the Hamiltonian varies relatively slowly), the Trotter formula can be used to approximate the time evolution superoperators as $\mathcal{U}_{j:j-1}\simeq\mathcal{V}^{1/2}_{j:j-1}\mathcal{W}_{j:j-1}\mathcal{V}^{1/2}_{j:j-1}$, with $\mathcal{V}_{j:j-1}$ generated by the $S$ part of the Hamiltonian and $\mathcal{W}_{j:j-1}$ by the remainder. Expanding out the superoperators appearing inside the trace in Eq. (A) and introducing further resolutions of the identity, one finds

[TABLE]

where $\mathcal{W}^{(s^{\prime},r^{\prime},s,r)}_{j:j-1}:=\bra{s^{\prime}}\mathcal{W}_{j:j-1}\left[\outerproduct{s}{r}\right]\ket{r^{\prime}}$, and we have used that $\bra{t}\mathcal{V}^{1/2}_{j:j-1}\left[\outerproduct{s}{r}\right]\ket{u}=\bra{s}{\mathcal{V}^{*\,1/2}_{j:j-1}}\left[\outerproduct{t}{u}\right]\ket{r}$. Assuming a factorizing initial condition $\chi_{0}=\rho_{0}\otimes\tau$ and substituting Eq. (B) into Eq. (A) leads to the following slightly more general version of Eq. (5) for the approximate process tensor $\tilde{\Upsilon}_{k:0}$:

[TABLE]

with

[TABLE]

In the special case that the bath coupling part Hamiltonian can be written in the form $H_{B}(t)=\sum_{s}\outerproduct{s}{s}\otimes B_{s}(t)$ where $\sum_{s}\outerproduct{s}{s}=\mathbbm{1}_{S}$ (the Hamiltonian of the main text takes this form in the interaction picture with respect to the bath), $\mathcal{W}_{j:j-1}^{(s^{\prime},r^{\prime},s,r)}=\delta_{ss^{\prime}}\delta_{rr^{\prime}}\mathcal{W}_{j:j-1}^{(s,r)}$ with

[TABLE]

and

[TABLE]

For a time-independent Hamiltonian with even time spacing $\delta t$, Eq. (B) then reduces to the operator representation of the discretized Feynman-Vernon influence functional in Eq. (4) of the main text.

Appendix C Explicit form of influence tensors in the spin-boson model

Further decomposing the influence functional into a product of the form of Eq. (5) requires that the bath be composed of field modes coupled linearly to the system, and that the Hamiltonian and initial state are quadratic in the corresponding creation and annihilation operators $\{\hat{a}^{\dagger}_{n}\}$ and $\{\hat{a}_{n}\}$. In other words, the environment must be Gaussian and the operators coupling to the system must take the form $\hat{B}_{s}(t)=\sum_{n}(g_{s,n}\hat{a}_{n}e^{-i\omega_{n}t}+g_{s,n}^{*}\hat{a}_{n}^{\dagger}e^{i\omega_{n}t})$ in the interaction picture. In this case, Wick’s theorem can be applied to express Eq. (B) as a product of exponentiated two point correlation functions. Specifically, we use the fact that for any linear functional of bath operators $\hat{X}$, $\tr\{T\exp[\hat{X}]\rho\}=\exp[\frac{1}{2}\tr\{T\hat{X}^{2}\rho\}]$, with $T$ any time ordering operator, when $\rho$ is Gaussian Diósi and Ferialdi (2014). Treating the left and right appended operators in Eq. (18) as a single contour ordered exponential under the trace, this results in the following expression:

[TABLE]

where

[TABLE]

for $i\neq j$, else for $i=j$:

[TABLE]

Restricting to the spin-boson type Hamiltonian considered in the main text, with a single interaction term $\hat{s}\sum_{n}(g_{n}\hat{a}_{n}+g_{n}^{*}\hat{a}_{n}^{\dagger})$ and with a thermal (and hence Gaussian) initial bath state $\tau_{\beta}$, we can compute the influence tensors explicitly. Here, the interaction picture bath operators appearing in Eqs. (21) and (22) take the simple form $\hat{B}_{u}(t)=\lambda_{u}\sum_{n}(g_{n}\hat{a}_{n}e^{-i\omega_{n}t}+g_{n}^{*}\hat{a}_{n}^{\dagger}e^{i\omega_{n}t})$, written in terms of the eigenvalues of the system operator $\hat{s}=\sum_{u}\lambda_{u}\outerproduct{u}{u}$.

In this case, $\zeta^{(u,v)}_{i,j}=2\lambda_{u}\lambda_{v}\eta_{i-j}$, where the memory kernel elements

[TABLE]

depend only on the difference between time steps and are expressed in terms of the environment auto-correlation function Strathearn et al. (2018, 2017)

[TABLE]

Here $J(\omega)=\sum_{n}|g_{n}|^{2}\delta(\omega-\omega_{n})$ is the spectral density, as defined in the main text. In terms of these quantities, the elements of the influence tensors $b_{(i-j)}$ that correspond to each of the terms in the exponentiated sum in Eq. (20) (such that $\mathcal{F}_{k:0}^{\alpha_{k}...\alpha_{1}}=\prod_{i\geq j}[b_{(i-j)}]^{\alpha_{i}\alpha_{j}}$), can be written

[TABLE]

These are the values that enter directly into our algorithm.

Appendix D Tensor network compression

In contracting the network, efficiency is achieved by finding a minimal approximate representation for the boundary matrix product operator in each iteration Schollwöck (2011). This is obtained by replacing high rank tensors with small singular values by lower rank approximations. These are found by performing a singular value decomposition (SVD) of the local tensors in the matrix product state, and discarding the singular values below a cutoff $\lambda_{c}$. Indicating an index partition by raised and lowered indices, the SVD decomposed tensor takes the form

[TABLE]

where the diagonal matrix $\Lambda^{\gamma}{}_{\delta}$ contains the singular values, and $U$, $V$ are rectangular isometric matrices satisfying $U^{\dagger}U=\mathbbm{1}$ and $V^{\dagger}V=\mathbbm{1}$ (see Fig. 4a). Truncating the singular values reduces the sizes of $U$ and $V$ (we will refer to the truncated versions as $\bar{U}$ and $\bar{V}$), and introduces a corresponding truncation error, whose magnitude is determined by the cutoff. We truncate the singular values such that

[TABLE]

where $\tilde{\Lambda}$ denotes the truncated diagonal matrix. Having truncated the singular values, we contract the $\bar{U}$ and $\tilde{\Lambda}$ to give a new matrix $Q$ with which to express the newly compressed local tensor:

[TABLE]

The singular value compression proceeds from one end of the boundary matrix product operator. Say we begin the compression at the right boundary (see Fig. 2), then initially the first (farthest to the right) local tensor is singular value decomposed and compressed. The matrix $\bar{U}$ in the above decomposition is then contracted with the diagonal matrix $\bar{\Lambda}$ to give $Q$, which is subsequently contracted with the second local tensor to the left. The compressed $\bar{V}$ tensor is stored as the new first local tensor (see Fig.4b) and this procedure is repeated for the second local tensor, and so on until the left boundary is reached. This constitutes a left sweep of the SVD compression procedure. After this left sweep, an equivalent right sweep is performed, where the left-most tensor is singular value decomposed and compressed, followed by the next left-most and so on. To produce the figures in this paper, we implemented one left sweep and one right sweep in each stage of the algorithm. Including more sweeps back and forth would in principle improve the quality of the compression; however, as we discuss below our implementation is sufficient to demonstrate an improvement of the local over the non-local algorithm.

Appendix E Scaling of memory effects with physical parameters

We now proceed to estimate the complexity of contracting the network. This depends crucially on how quickly the memory decays, and hence the effective depth of the tensor network in Fig. 2c (as we will see, the overall size of memory effects is also important). From Eq. (25), it is clear that non-trivial contributions of the $b_{(i-j)}$ tensors to the influence functional depend on the magnitude of the memory kernel elements $\eta_{i-j}$. Specifically, the effect of truncating the network in Fig. 2c at a depth $m$, with resulting influence functional elements $[\mathcal{F}^{\alpha_{k}\dots\alpha_{1}}_{k:0}]_{m}$, is to introduce a relative error:

[TABLE]

where $\|\hat{s}\|_{\rm op}:=\max{\{|\lambda_{r}|\}}$ (we will henceforth take $\|\hat{s}\|_{\rm op}=1$, effectively absorbing it into the coupling strength); in the second line we have assumed $m$ is sufficiently large that the error is small. For a fixed error $\varepsilon$, the memory time $t_{m}=m\delta t$, and hence the complexity of our algorithm (we expect the error due to SVD compression to scale similarly), will therefore depend on how quickly $|\eta_{l}|$ decays with $l$. If it decays exponentially with rate $c$, then in the limit of large $k$, it is relatively straightforward to show that the memory time scales as $t_{m}\sim\delta t(\log k+\log\varepsilon^{-1})/c$. However, as we will now see, for the spectral density we have chosen, the memory kernel decays as a power law.

When $\delta t=t_{j}-t_{j-1}$ is sufficiently small, we have

[TABLE]

For the spectral density introduced in the main text $J(\omega)=(\alpha\omega_{c}/2)(\omega/\omega_{c})^{\nu}\exp\left(-\omega/\omega_{c}\right)$ with $\nu=1$ (i.e. the Ohmic case), the integral in Eq. (24) can be evaluated explicitly, giving

[TABLE]

with $\psi^{(1)}[z]:=\int_{0}^{\infty}\,dx\,xe^{-zx}/(1-e^{-x})$ the order-1 polygamma function. For large $|z|$, the latter goes as $\psi^{(1)}[z]\sim 1/z+1/(2z^{2})$ Abramowitz and Stegun (1964). Hence, in the limit that $t\gg\omega_{c}^{-1}$ and $t\gg\beta$, we can expand out Eq. (E) and combine with Eq. (30) to arrive at

[TABLE]

to leading order in $|i-j|^{-1}$. Therefore, for sufficiently large $m$ and $k=t_{\rm max}/\delta t$, we can perform the sum in Eq. (E), finding $\varepsilon_{m}\lesssim\alpha k\psi^{(1)}[m]/(\pi\beta\omega_{c})\simeq\alpha k/(\pi\beta\omega_{c}m)$. For fixed error $\varepsilon$, we therefore have that the bound on the memory time, and hence the complexity of the algorithm scales as

[TABLE]

This explains the behaviour of Fig. 3a and the large $\omega_{c}$ behaviour of Fig. 3b in the main text. However, the onset of this limit depends on the parameter combinations $\omega_{c}t$ and $\sqrt{(t^{2}+\omega_{c}^{-2})/\beta^{2}}$ both being large (these elicit expansions for the first and second terms of Eq. (E) respectively). When $\omega_{c}\lesssim t_{\rm max}^{-1}$, the former limit is never reached, and only the term involving the polygamma function contributes significantly to the error. Specifically, the magnitude of the memory kernel is approximately constant, going as

[TABLE]

Therefore, for small $\delta t$, Eq. (E) leads to $\varepsilon_{m}\lesssim(\alpha\omega_{c}/\pi\beta)(t_{\rm max}-t_{m})^{2}$ and one can see that, as long as $t_{\rm max}$ is fixed, the error is bounded by a number that goes to zero as $\omega_{c}$ does, even for very small memory times. Hence, in this limit, the effective coupling to the bath is weak overall and even the local influence tensors do not contribute significantly to the dynamics, explaining the behaviour at small $\omega_{c}$ in Fig. 3b.

Appendix F Scaling comparison

In Fig. 4c we compare the computation time per iteration of the non-local and the local network representations for a fixed network size and SVD cutoff. We see that the time per iteration of the non-local TEMPO algorithm increases approximately linearly. For a finite network, the local TEMPO algorithm has a time per iteration which rapidly goes to a non-increasing value. The growth of complexity in the non-local case is mainly due to the build up of irrelevant information, rather than a genuine build-up of temporal correlations. The observed decrease in the computation time per iteration for the local case, is a consequence of working with a finite network. In the original TEMPO proposal Strathearn et al. (2018), it was argued that one could implement a truncation of the number of tensors in the propagators and obtain a constant scaling at long times. The same method could be applied with the local TEMPO algorithm, only with a significantly reduced time per iteration. More rigorously, we could combine the tools developed here with the transfer tensor approach Cerrillo and Cao (2014); Pollock and Modi (2018), which infers long-time correlations from a short-time simulation. The problem would then become efficiently contracting the full network up to a sufficiently long-time.

The advantage in contracting the network persists across a wide range of parameters, as depicted in Fig. 4d. Even in the easier regime of weak coupling and small cutoff frequency, the non-local network takes longer to contract than the local one.

All the simulations presented in this work were carried out on a 2011 MacBook Pro with a 2,4 GHz Intel Core i5 processor and a 4GB 1333 MHz DDR3 memory. The algorithm has been implemented using the programming language Python, no specialized packages, beyond NumPy, were used. There is extensive scope for optimization, and utilizing packages for tensor network manipulations would supposedly make the implementation quite simple.