Polynomial complexity of open quantum system problems

TL;DR

This paper proves that the complexity of modeling open quantum systems grows at most polynomially with time and bath size, enabling efficient algorithms for simulating their dynamics using tensor networks.

Contribution

It establishes a polynomial complexity theorem for open quantum systems and introduces tensor network algorithms with linear bond dimension growth.

Findings

Number of equations grows linearly with time

Tensor network bond dimension increases linearly with time

Algorithms successfully applied to spin-boson and Gaudin models

Abstract

Open quantum systems (OQS's) are ubiquitous in non-equilibrium quantum dynamics and in quantum science and technology. Solving the dynamics of an OQS in a quantum many-body bath has been considered a computationally hard problem because of the dimensionality curse. Here, considering that full knowledge of the bath dynamics is unnecessary for describing the reduced dynamics of an OQS, we prove a polynomial complexity theorem, that is, the number of independent equations required to fully describe the dynamics of an OQS increases at most linearly with the evolution time and polynomially with the bath size. Therefore, efficient computational algorithms exist for solving the dynamics of a small-sized OQS (such as a qubit or an atom). We further prove that, when the dynamics of an OQS and the bath is represented by a tensor network, a tensor contraction procedure can be specified such that the bond dimension (i.e., the range of tensor indices contracted in each step) increases only linearly (rather than exponentially) with the evolution time, providing explicitly efficient algorithms for a wide range of OQS's. We demonstrate the theorems and the tensor-network algorithm by solving two widely encountered OQS problems, namely, a spin in a Gaussian bath (the spin-boson model) and a central spin coupled to many environmental spins (the Gaudin model). This work provides approaches to understanding dynamics of OQS's, learning the environments via quantum sensors, and optimizing quantum information processing in noisy environments.

Equations

$$H(t)=\sum _{\alpha =1}^{d^2-1}{S}_{\alpha }(t)\otimes B_{\alpha }(t),$$

$$\chi (T)=e^{ϵ{\mathcal{H}}^{-}(T)}\cdots e^{ϵ{\mathcal{H}}^{-}(2ϵ)}{e}^{ϵ{\mathcal{H}}^{-}(ϵ)}\chi (0).$$

$$e^{ϵ{\mathcal{H}}^{-}(t)}\approx {\mathcal{I}}_{\mathrm{S}}\otimes {\mathcal{I}}_{\mathrm{B}}+ϵ\sum _{\sigma =\pm }\sum _{\alpha =1}^{d^2-1}{\mathcal{S}}_{\alpha }^{\sigma }(t)\otimes {\mathcal{B}}_{\alpha }^{\bar{\sigma }}(t)\equiv {\mathcal{S}}^{\varphi }(t)\otimes {\mathcal{B}}_{\varphi }(t),$$

$$\rho (T)\equiv {\text{Tr}}_{\mathrm{B}}[\chi (T)]={\mathbb{C}}_{{\mathbf{\Phi }}^{(T)}}{\mathbb{S}}^{{\mathbf{\Phi }}^{(T)}}\rho (0),$$

$${\mathbb{C}}_{{\mathbf{\Phi }}^{(T)}}\equiv {\mathrm{Tr}}_{\mathrm{B}}\left[{\mathcal{B}}_{{\varphi }_T}(T)\cdots {\mathcal{B}}_{{\varphi }_2}(2){\mathcal{B}}_{{\varphi }_1}(1)\Omega \right].$$

$${\mathbb{C}}_{{\mathbf{\Phi }}^{(t)}}\equiv \sum _{\mathbf{\Theta }\subseteq {\mathbf{\Phi }}^{(t)}}{\delta }_{\bar{\mathbf{\Theta }}}^0\sum _{\left\{{\mathbf{\Theta }}_i\right\}}{\tilde{C}}_{{\mathbf{\Theta }}_1}{\tilde{C}}_{{\mathbf{\Theta }}_2}\cdots {\tilde{C}}_{{\mathbf{\Theta }}_n},$$

$${\rho }^{{\mathbf{\Phi }}^{(t)}}\equiv {\mathbb{C}}_{{\mathbf{\Phi }}^{\prime (t)}}\left[\prod _{\tau =1}^t{P}_{{\varphi }_{\tau }^{\prime \prime }}^{{\varphi }_{\tau },{\varphi }_{\tau }^{\prime }}\right]{\mathbb{S}}^{{\mathbf{\Phi }}^{\prime \prime (t)}}\rho (0)\equiv {\mathbb{C}}_{{\mathbf{\Phi }}^{\prime (t)}}^{{\mathbf{\Phi }}^{(t)}}{\mathbb{S}}^{{\mathbf{\Phi }}^{\prime (t)}}\rho (0),$$

$${\rho }^{{\mathbf{\Phi }}^{{}^{\prime }(t)}}={\mathbb{K}}_{{\varphi }_t;{\mathbf{\Phi }}^{(t-1)}}^{{\varphi }_t^{\prime };{\mathbf{\Phi }}^{{}^{\prime }(t-1)}}{\mathcal{S}}^{{\varphi }_t}{\rho }^{{\mathbf{\Phi }}^{(t-1)}}.$$

$${\rho }^{{\mathbf{\Phi }}^{\prime (t)}}$$

$$={\text{Tr}}_{\mathrm{B}}\left[{\mathcal{S}}^{{\varphi }_t^{\prime }}{e}^{ϵ{\mathcal{H}}^{-}(t)}\cdots {\mathcal{S}}^{{\varphi }_1^{\prime }}{e}^{ϵ{\mathcal{H}}^{-}(ϵ)}\Omega (0)\right]\rho (0).$$

$$N_T\sim d^{2}T/ϵ,$$

$${\mathbf{R}}^a={\left[r\right]}_{{\mathbf{\Phi }}^{(T)}}^a{\mathbf{\rho }}^{{\mathbf{\Phi }}^{(T)}},$$

$$\rho (t)={\left[A_{{\varphi }_t}\right]}_{a_{t-1}}{\left[A_{{\varphi }_{t-1}}\right]}_{a_{t-1},a_{t-2}}\cdots {\left[A_{{\varphi }_2}\right]}_{a_2,a_1}{\left[A_{{\varphi }_1}\right]}_{a_1}{\mathbb{S}}^{{\mathbf{\Phi }}^{(t)}}\rho (0),$$

$${\rho }^{{\mathbf{\Phi }}^{(t)}}={\left[A_{{\varphi }_t^{\prime }}^{{\varphi }_t}\right]}_{a_{t-1}}{\left[A_{{\varphi }_{t-1}^{\prime }}^{{\varphi }_{t-1}}\right]}_{a_{t-1},a_{t-2}}\cdots {\left[A_{{\varphi }_2^{\prime }}^{{\varphi }_2}\right]}_{a_2,a_1}{\left[A_{{\varphi }_1^{\prime }}^{{\varphi }_1}\right]}_{a_1}{\mathbb{S}}^{{\mathbf{\Phi }}^{{}^{\prime }(t)}}\rho (0),$$

$$D_{\tau }\le N_{\tau }\sim d^2\tau /ϵ.$$

$$R_t^a={\left[r\right]}_{{\mathbf{\Phi }}^{(t)}}^a{\rho }^{{\mathbf{\Phi }}^{(t)}},\mathrm{and}{\rho }^{{\mathbf{\Phi }}^{(t)}}={\left[\bar{r}\right]}_a^{{\mathbf{\Phi }}^{(t)}}{R}_t^a.$$

$${\mathbb{C}}_{{\mathbf{\Phi }}^{{}^{\prime }(t)}}\sim {\left[A_{{\varphi }_t^{\prime }}\right]}_{a_{t-1}}{\left[A_{{\varphi }_{t-1}^{\prime }}\right]}_{a_{t-1},a_{t-2}}\cdots {\left[A_{{\varphi }_2^{\prime }}\right]}_{a_2,a_1}{\left[A_{{\varphi }_1^{\prime }}\right]}_{a_1},$$

$${\mathbb{C}}_{{\mathbf{\Phi }}^{{}^{\prime }(t)}}^{{\mathbf{\Phi }}^{(t)}}\sim {\left[A_{{\varphi }_t^{\prime }}^{{\varphi }_t}\right]}_{a_{t-1}}{\left[A_{{\varphi }_{t-1}^{\prime }}^{{\varphi }_{t-1}}\right]}_{a_{t-1},a_{t-2}}\cdots {\left[A_{{\varphi }_2^{\prime }}^{{\varphi }_2}\right]}_{a_2,a_1}{\left[A_{{\varphi }_1^{\prime }}^{{\varphi }_1}\right]}_{a_1}.$$

$$D_{\tau }\le \min (N_{\tau },N_{T_m}).$$

$$H(t)=S_z(t)\sum _k{g}_k\left(a_k^{†}{e}^{i{\omega }_{k}t}+h.c.\right)\equiv S_z(t)B(t),$$

$${\tilde{C}}_{\mathbf{\Theta }}=\sum _{k=1}^K{\tilde{C}}_{k;\mathbf{\Theta }},$$

$${\rho }_{L+1}^{{\mathbf{\Phi }}^{(T)}}={\mathbb{C}}_{L+1;{\mathbf{\Phi }}^{{}^{\prime }(T)}}\left[\prod _{\tau =1}^T{P}_{{\varphi }_{\tau }^{\prime \prime }}^{{\varphi }_{\tau },{\varphi }_{\tau }^{\prime }}\right]{\rho }_L^{{\mathbf{\Phi }}^{{}^{\prime \prime }(T)}},$$

$$H(t)=\sum _{k=1}^K{g}_k\mathbf{S}\cdot {\mathbf{J}}_k,$$

$${\mathbb{C}}_{{\mathbf{\Phi }}^{(t)}}^{{\mathbf{\Phi }}^{{}^{\prime }(t)}}$$

$${\tilde{K}}_{{\mathbf{\Phi }}^{(t)}}^{{\mathbf{\Phi }}^{{}^{\prime }(t)}}\left(\mathbf{\Theta }\right)$$

$$P_{\phi }^{{\varphi }^{\prime },{\varphi }^{\prime \prime \prime }}{P}_{\varphi }^{\phi ,{\varphi }^{\prime \prime }}=P_{\varphi }^{{\varphi }^{\prime \prime \prime },\phi }{P}_{\phi }^{{\varphi }^{\prime \prime },{\varphi }^{\prime }},$$

$${\mathbb{K}}_{{\varphi }_{t+1}^{\prime \prime };{\mathbf{\Phi }}^{{}^{\prime \prime }(t)}}^{{\varphi }_{t+1}^{\prime };{\mathbf{\Phi }}^{{}^{\prime }(t)}}\equiv {\left[\prod _{\mathbf{\Theta }\subseteq {\mathbf{\Phi }}^{(t)}}\tilde{K}\left({\varphi }_{t+1};\mathbf{\Theta }\right)\right]}_{{\mathbf{\Phi }}^{{}^{\prime \prime }(t+1)}}^{{\mathbf{\Phi }}^{{}^{\prime }(t+1)}},$$

$${\mathbb{C}}_{{\mathbf{\Phi }}^{(t+1)}}^{{\mathbf{\Phi }}^{{}^{\prime }(t+1)}}={\mathbb{K}}_{{\varphi }_{t+1};{\mathbf{\Phi }}^{{}^{\prime \prime }(t)}}^{{\varphi }_{t+1}^{\prime };{\mathbf{\Phi }}^{{}^{\prime }(t)}}{\mathbb{C}}_{{\mathbf{\Phi }}^{(t)}}^{{\mathbf{\Phi }}^{{}^{\prime \prime }(t)}},$$

$${\mathbb{P}}_{{\mathbf{\Phi }}^{(t)}}^{{\mathbf{\Phi }}^{{}^{\prime }(t)},{\mathbf{\Phi }}^{{}^{\prime \prime }(t)}}\sum _{{\mathbf{\Theta }}^{\prime \prime }\subseteq {\mathbf{\Phi }}^{{}^{\prime \prime }(t)}}{\delta }_{{\bar{\mathbf{\Theta }}}^{\prime \prime }}^0\sum _{\{{\mathbf{\Theta }}_i^{\prime \prime }\}}{\tilde{C}}_{{\mathbf{\Theta }}_1^{\prime \prime }}{\tilde{C}}_{{\mathbf{\Theta }}_2^{\prime \prime }}\cdots {\tilde{C}}_{{\mathbf{\Theta }}_n^{\prime \prime }}={\left[\prod _{\mathbf{\Theta }\subseteq {\mathbf{\Phi }}^{(t)}}\tilde{K}\left(\mathbf{\Theta }\right)\right]}_{{\mathbf{\Phi }}^{(t)}}^{{\mathbf{\Phi }}^{{}^{\prime }(t)}},$$

$${\left[\tilde{K}\left({\mathbf{\Theta }}_u\right)\tilde{K}\left({\mathbf{\Theta }}_d\right)\right]}_{{\mathbf{\Phi }}^{(t)}}^{{\mathbf{\Phi }}^{{}^{\prime }(t)}}=$$

Full text

I Introduction

Quantum systems inevitably couple with their environments or baths, becoming open quantum systems (OQS’s) [1, 2]. The OQS problem is of fundamental interest for understanding quantum decoherence [3, 4]. It is ubiquitous in non-equilibrium quantum dynamics, including quantum transport [5, 6, 7] and quantum optical processes [8]. It has become a critical problem for the development of modern quantum technologies [9, 10, 11, 12] such as quantum metrology, quantum communication, and quantum computing, where the coupling of quantum systems to their environments is a main obstacle to maintaining quantum superpositions and scaling up the systems. In quantum sensing, on the other hand, an OQS can act as a quantum probe to its environment [13, 14, 15, 16]. To design quantum systems and to optimize quantum controls of such systems [17, 18, 19, 20] for mitigating the decoherence effects and utilizing the quantum advantages it requires good understanding and often the computation of the dynamics of OQS’s.

The solution of the dynamics of an OQS in general is believed to be a computationally hard problem (time consumption and/or memory demand increasing at least exponentially with the evolution time and the bath size), since it is assumed to be equivalent to solving the dynamics of the quantum many-body system including the bath. In cases where environments can be taken as white noise, their effects are captured by the Lindblad master equation, known as Markovian dynamics [21]. However, in general, memory (i.e., non-Markovian) effects are important [22, 23], especially when the quantum system is strongly coupled to its environment or the environment itself has a structure. Notable examples include a localized electron spin coupled with nuclear spins [24, 25], an emitter in a photonic crystal [26], and a multi-qubit system with long-range correlated noises [27]. The memory effects mean that extra many-body correlations can build up in the baths via coupling to the central system. For example, even a group of non-interacting spins can accumulate many-body correlations when they are commonly coupled to a central spin (such as in the Gaudin model [28, 29, 30, 31, 32]). Therefore, the non-Markovian dynamics of OQS’s is highly non-trivial.

Previous studies indicate that numerically efficient methods may exist for certain OQS’s. In particular, when system evolution and environment dynamics have very different timescales, the problem can be efficiently solved using perturbative theories such as the time-convolutionless projection operator technique [1] and the cluster expansion methods [33]. Beyond the perturbation theories, several exact methods have been developed to study OQS’s where the environments exhibit simple structures. The pseudo-mode method [34, 35, 36, 37, 38, 39], the stochastic Schrödinger equation [40, 41, 42, 43], the hierarchical equations of motion [44, 45], and the hierarchy of pure states [46, 47] are applicable when the environment is Gaussian and has certain spectral characteristics. Path-integral influence functionals [48, 49, 50, 51, 52, 53] and time-evolving matrix product operators [54, 55, 56] can deal with Gaussian environments that have special system-environment coupling. In addition, the process tensor networks [57, 58] provide frameworks for solving general OQS problems but the sampling of process tensors can be classically hard [59]. On the other hand, advanced numerical algorithms developed for quantum many-body systems, such as the numerical renormalization group method [60], the density matrix renormalization group method [61, 30], and the tensor network method [62, 63], have been applied to investigate OQS’s in strongly correlated environments [64, 65, 66, 67], which have significantly mitigated the difficulties of non-Markovian dynamics [61, 58, 50, 51, 52, 53, 68]. Furthermore, it has been established (though only formally) that the reduced dynamics of a central quantum system can be fully captured by the correlations in the bath before it is coupled to the central system [69, 14]; in particular, the dynamics of a spin coupled to a number of harmonic oscillators is fully determined by the spectral density of the harmonic oscillators and is independent of the details of the couplings (for instance, a bath of one harmonic oscillator and that of $10^{3}$ harmonic oscillators of the same frequency make no difference in the central spin dynamics as long as the summation of the coupling strengths squared is equal in the two cases, but the dynamics of the whole spin-bath system would be much more complicated in the latter case). All these results suggest that the large number of degrees of freedom of a bath may not necessarily prevent an efficient solution of the non-Markovian dynamics of an OQS [70, 71, 72]. However, what the ultimate complexity of the non-Markovian dynamics is and how efficient algorithms, if they exist, can be constructed systematically are still open questions.

Here, we establish a polynomial complexity theorem that the non-Markovian dynamics of a general OQS can be fully described by a set of independent equations of extended density matrices (EDMs), and the size of the set grows at most linearly with the evolution time and polynomially with the bath size. This means that the complexity of the non-Markovian dynamics increases at most as a polynomial function of the evolution time and the bath size. Moreover, we propose an explicit tensor network algorithm to solve these EDMs when the correlations of environmental noise are known (from computation or measurement). We demonstrate the theorem and the efficiency of the tensor network algorithm by solving two widely encountered classes of OQS problems, namely a spin in a Gaussian bath (the spin-boson model) and a central spin coupled to many environmental spins (the Gaudin model). The numerical results show that the bond dimension of the EDM tensor increases at most linearly with time evolution, regardless of the coupling strength (in the spin-boson model) and the bath size (in the Gaudin model), verifying the efficiency of the tensor network algorithm.

II OQS dynamics in terms of bath correlations

For a general OQS, the system-bath coupling can be written in the interaction picture with respect to the bath Hamiltonian ${H}_{\rm B}$ and the system Hamiltonian $H_{\rm S}$ as

[TABLE]

where $d$ is the dimension of the system’s Hilbert space, $\left\{S_{\alpha}(t)\equiv e^{iH_{\rm S}t}S_{\alpha}e^{-iH_{\rm S}t}\right\}$ is a basis of the system operators (with the trivial identity operator $S_{0}\equiv{I}_{S}$), and $B_{\alpha}(t)\equiv e^{i{H}_{\rm B}t}B_{\alpha}e^{-i{H}_{\rm B}t}$ is the bath operator coupled to ${S}_{\alpha}$. In principle, the reduced dynamics of the central system can be obtained by first solving the coupled system-bath dynamics (given by the density operator $\chi(t)$) and then carrying out the partial trace of the bath’s degrees of freedom [1]. By discretizing the time as $0,\epsilon,2\epsilon,\ldots,T$ with a small time step $\epsilon\equiv T/N$, the evolution of the total system reads

[TABLE]

Here, we have introduced the superoperators $\mathcal{A}^{-}\chi\equiv-i(A\chi-\chi A)$ (essentially a commutator) and $\mathcal{A}^{+}\chi\equiv\frac{1}{2}(A\chi+\chi A)$ (essentially an anti-commutator) for an operator $A$. A small step of evolution can be expanded to the first order of $\epsilon$ as

[TABLE]

where $\bar{\sigma}\equiv-\sigma$ and a summation over $\phi=0,1,\ldots,2(d^{2}-1)$ is implied in the tensor contraction convention with ${\mathcal{S}}^{2{\alpha}-1}\equiv\epsilon{\mathcal{S}}_{\alpha}^{+}$, ${\mathcal{S}}^{2{\alpha}}\equiv\epsilon{\mathcal{S}}_{\alpha}^{-}$, ${\mathcal{B}}_{2{\alpha}-1}\equiv{\mathcal{B}}_{\alpha}^{-}$, ${\mathcal{B}}_{2{\alpha}}\equiv{\mathcal{B}}_{\alpha}^{+}$, and the null evolutions ${\mathcal{S}}^{0}\equiv{\mathcal{I}}_{\rm S}$ for the system and ${\mathcal{B}}_{0}\equiv{\mathcal{I}}_{\rm B}$ for the bath. The physical meanings of the superoperators are clear by checking their effects on the density operators [73]. For example, ${\mathcal{S}_{\alpha}^{-}}\rho=-i\left[S_{\alpha},\rho\right]$ causes the evolution of the system and ${\rm Tr}[{\mathcal{B}}_{\alpha}^{+}\Omega]=\left\langle B_{\alpha}\right\rangle$, so the term ${\mathcal{S}}_{\alpha}^{-}\otimes{\mathcal{B}}_{\alpha}^{+}\chi$ corresponds to the evolution of the system driven by the $\alpha$-th component of the force from the bath. Similarly, ${\mathcal{S}}_{\alpha}^{+}\otimes{\mathcal{B}}_{\alpha}^{-}\chi$ means the evolution of the bath under the $\alpha$-th component of the force from the system. A small step of evolution can be presented by a tensor network diagram, as illustrated in Fig. 1a.

For the sake of simplicity, we consider an initial product state $\chi(0)=\rho(0)\otimes\Omega$ with $\Omega$ and $\rho(0)$ being the initial bath and system states, respectively. The state of the system can be written as [69, 14]

[TABLE]

with ${\bm{\Phi}}^{(T)}\equiv\left\{\phi_{T},\ldots,\phi_{2},\phi_{1}\right\}$, where the time-ordered tensor ${\mathbb{S}}^{{\bm{\Phi}^{(T)}}}\equiv{\mathcal{S}}^{\phi_{t}}(T)\cdots{\mathcal{S}}^{\phi_{2}}(2){\mathcal{S}}^{\phi_{1}}(1)$ denotes the evolution of the system driven by the corresponding correlation in the bath denoted by the correlation tensor

[TABLE]

Hereafter we set the time step $\epsilon$ as unity so that $t$ stands for $\epsilon t$ and $T=N$. The corresponding tensor-network representation is shown in Fig. 1b. Note that the result in Eq. (4) can be generalized to a correlated initial state $\chi(0)=\sum_{\alpha}\rho_{\alpha}(0)\otimes\Omega_{\alpha}$ with $\Omega_{\alpha}$ being the $\alpha$-th state of the bath and $\rho_{\alpha}$ a system operator [74]. An important feature of Eq. (4) is that the reduced dynamics of the OQS is fully determined by the correlation functions of the bath without including the system-bath interaction as can be seen from Eq. (5) [69, 14]. Thus, given the knowledge of the correlation functions of the bath, the dynamics of the system can be calculated with the complexity independent of the size of the bath. However, solving ${\mathbb{C}}_{\bm{\Phi}^{(T)}}$ is still a many-body problem as it captures the dynamics of the bath [75, 59]. In fact, the number of correlation functions, $(2d^{2}-1)^{T}$, increases exponentially with the evolution time. Therefore, the direct evaluation of the tensor network in Fig. 1b is computationally hard.

Usually, the correlation tensor defined in Eq. (5) has a high degree of redundancy. For example, very often, the correlation functions can be reduced to a few irreducible correlation functions, i.e., cumulants and the cumulants have finite memory time. The cumulants $\tilde{C}_{{\bm{\Theta}}_{i}}$ are recursively defined as

[TABLE]

where $\left\{{\bm{\Theta}}_{i}\right\}$ are non-overlapping non-empty subsets of indices selected from ${\bm{\Phi}}^{(t)}$ with $\bigcup_{i=1}^{n}{\bm{\Theta}}_{i}={\bm{\Theta}}\subseteq{\bm{\Phi}}^{(t)}$, $\bar{\bm{\Theta}}\equiv{\bm{\Phi}}^{(t)}\setminus{\bm{\Theta}}$ is complement of the subset ${\bm{\Theta}}$, and the summation runs over all subsets ${\bm{\Theta}}$ and all possible partitions of ${\bm{\Theta}}$. For example, the cumulnant expansion up to the third-order correlation tensors reads ${\mathbb{C}}_{\phi_{1}}\equiv\tilde{C}_{\phi_{1}}+\delta^{0}_{\phi_{1}}$, ${\mathbb{C}}_{\phi_{2},\phi_{1}}\equiv\tilde{C}_{\phi_{2},\phi_{1}}+\tilde{C}_{\phi_{2}}\tilde{C}_{\phi_{1}}+\delta^{0}_{\phi_{2}}\tilde{C}_{\phi_{1}}+\delta^{0}_{\phi_{1}}\tilde{C}_{\phi_{2}}+\delta^{0}_{\phi_{2}}\delta^{0}_{\phi_{1}}$, and ${\mathbb{C}}_{\phi_{3},\phi_{2},\phi_{1}}\equiv\tilde{C}_{\phi_{3},\phi_{2},\phi_{1}}+\tilde{C}_{\phi_{3},\phi_{2}}\tilde{C}_{\phi_{1}}+\cdots+\delta^{0}_{\phi_{3}}\delta^{0}_{\phi_{2}}\delta^{0}_{\phi_{1}}$. Note that the cumulant $\tilde{C}_{\bm{\Theta}}=0$ when any index $\phi_{\tau}\in{\bm{\Theta}}$ takes zero value ($\phi_{\tau}=0$) since if $\phi_{\tau}\in{\bm{\Theta}}$ and $\phi_{\tau}=0$ then the correlation $C_{\bm{\Theta}}=C_{{\bm{\Theta}}\setminus\{\phi_{\tau}\}}$ (reduced to a lower-order correlation). The noise from the bath at a certain time $\mathcal{B}_{\phi_{t}}(t)$ can affect the non-Markovian dynamics of the system via different types of cumulants, as a time-local mean field (via the first-order irreducible correlation function) or in correlations with noises at earlier and/or later times (via second- or higher-order irreducible correlations), as illustrated in Fig. 1c. We introduce the picking tensor $P^{\phi^{\prime},\phi^{\prime\prime}}_{\phi}=\delta^{\phi^{\prime}}_{\phi}\delta^{\phi^{\prime\prime}}_{0}+\delta^{\phi^{\prime\prime}}_{\phi}\delta^{\phi^{\prime}}_{0}-\delta^{0}_{\phi}\delta^{\phi^{\prime}}_{0}\delta^{\phi^{\prime\prime}}_{0}$ to select how a noise at a certain time participates in non-Markovian dynamics (diagram shown in the circle in Fig. 1c). When $\phi=0$, the picking tensor selects the null evolution with $\phi^{\prime\prime}=\phi^{\prime}=0$; otherwise, it selects one cumulant by setting either $\phi^{\prime}$ or $\phi^{\prime\prime}$ equal to $\phi$ and the other upper index equal to 0. The cumulants in Eq. (6) are interconnected to form a tensor network, as illustrated in Fig. 1d by an example of three-step evolution. Note that to transfer the addition of cumulants in Eq. (6) into a tensor product form, we have replaced the cumulants with the irreducible kernel tensors $\tilde{K}_{{\bm{\Theta}}}\equiv\delta^{{\mathbf{0}}}_{{\bm{\Theta}}}+\tilde{C}_{{\bm{\Theta}}}$. One can demonstrate that the contraction of these kernel tensors interconnected by the picking tensors correctly yields the correlation tensor [see Appendix A for the proof].

An important feature of the non-Markovian dynamics of an OQS is that the noises from the bath at present and earlier times may take delayed effects in correlations with noises at later times. Thus, the reduced density operator $\rho(t)$ is insufficient to keep track of the dynamics of the OQS. Instead, one needs to bookmark different possibilities of how the noise at a certain time would be correlated with future noises, that is, in the tensor network diagram as illustrated in Fig. 1d, to keep the upward open arms stemming from the picking tensors. We denote a tensor network with such upward open arms as an extended density matrix (EDM) at a given time $t$

[TABLE]

where the upper indices ${{\bm{\Phi}}^{(t)}}$ keep track of the history (to account for the memory effect), and $\mathbb{C}^{{\bm{\Phi}}^{(t)}}_{\bm{\Phi}^{\prime(t)}}\equiv\mathbb{C}_{\bm{\Phi}^{\prime\prime(t)}}\prod^{t}_{\tau=1}{P}^{\phi_{\tau},\phi^{\prime}_{\tau}}_{\phi^{\prime\prime}_{\tau}}$ (referred to also as correlation tensor with a little abuse of notation). For example, the tensor in Fig. 1d represents the EDM $\rho^{\phi^{\prime}_{3},\phi^{\prime}_{2},\phi^{\prime}_{1}}$. We can further group all kernel tensors ended at time $t$ and the associated picking tensors, keeping track of the upward open arms, into a combined kernel tensor $\mathbb{K}^{\phi^{\prime}_{t};{\bm{\Phi}}^{{}^{\prime}(t-1)}}_{\phi_{t};{\bm{\Phi}}^{(t-1)}}$ [shadowed diagrams in Figs. 1d and  1e for illustration of ${\mathbb{K}}_{\phi_{3};\phi_{2},\phi_{1}}^{\phi^{\prime}_{3};\phi^{\prime}_{2},\phi^{\prime}_{1}}$].

The evolution of EDMs is driven by the combined kernel tensors as

[TABLE]

The tensor network representation of this evolution is shown in Fig. 1e and Fig. 2a. Equation (8) serves the role of differential equations for the EDMs, capturing fully the memory effects of the bath [22, 23]. The reduced dynamics of the OQS is obtained from the EDMs through $\rho(t)=\delta^{\mathbf{0}}_{{\bm{\Phi}}^{(t)}}\rho^{{\bm{\Phi}}^{(t)}}$ (by terminating the upward open arms with the closing tensor $\delta^{0}_{\phi}$, as shown in Fig. 2a). Note that the number of EDMs and hence the corresponding number of time-local equations still increases exponentially with the evolution time.

The EDMs have a physical meaning as the state of an OQS evolving under a sequence of weak interventions. Using the formula $\mathcal{S}^{\phi_{t}}\mathcal{S}^{\phi^{\prime}_{t}}=P^{\phi_{t},\phi^{\prime}_{t}}_{\phi^{\prime\prime}_{t}}\mathcal{S}^{\phi^{\prime\prime}_{t}}+O\left(\epsilon^{2}\right)$, we can demonstrate that the EDM tensor can be approximated as [see Fig. 2b for illustration by tensor-network diagrams]

[TABLE]

Therefore, the trace of the EDMs corresponds to different types and orders of the time-ordered correlation functions of the OQS defined by the sequence of super-operators ${\mathcal{S}}^{\phi^{\prime}_{1}},{\mathcal{S}}^{\phi^{\prime}_{2}},\ldots,{\mathcal{S}}^{\phi^{\prime}_{t}}$.

III Complexity of non-Markovian dynamics

III.1 Polynomial complexity theorem

The EDMs offer a general formulation of OQS dynamics using cumulants of baths, but solving the dynamics is still challenging since the number of EDMs increases exponentially with the evolution time. Here, we prove a polynomial complexity (PC) theorem, which establishes that the number of independent equations needed to fully describe the dynamics of an OQS increases at most linearly with the evolution time and polynomially with the bath size.

Theorem 1** (Polynomial complexity theorem).**

For a central quantum system (S) coupled to a generic quantum bath (B), the EDMs that describe the non-Markovian dynamics of the system S in a finite time $T$ can be expanded by a set of linearly independent vectors, and the size of the set increases at most linearly with the evolution time, i.e.,

[TABLE]

where $d$ is the dimension of the Hilbert space of the system S and $\epsilon$ represents the minimal time in which the bath B induces a discernible change in the reduced dynamics of the system S.

Proof.

We write the evolution path of an EDM as an operator-valued vector ${\bm{\rho}}^{{\bm{\Phi}}^{(T)}}\equiv\left[{\rho}^{{\bm{\Phi}}^{(T)}},\ldots,{\rho}^{{\bm{\Phi}}^{(2)}},{\rho}^{{\bm{\Phi}}^{(1)}}\right]$. Since each element ${\rho}^{{\bm{\Phi}}^{(t)}}$ in the vector is an Hermitian matrix and can be expanded by $d^{2}$ Hermitian system operators with $d^{2}$ independent real parameters (including the trace), each operator-valued EDM vector ${\bm{\rho}}^{{\bm{\Phi}}^{(T)}}$ can be expressed as a vector in a linear space of real dimension $T\times d^{2}$. Consequently, at most $N_{T}=Td^{2}$ independent vectors can be formed using all the EDMs, given by

[TABLE]

where ${\mathbf{R}}^{{a}}=\left[R_{T}^{{a}},\ldots,R^{{a}}_{2},R^{{a}}_{1}\right]$ with each component $R^{{a}}_{t}$ being a $d$-dimensional Hermitian matrix at time $t$, $\left[r\right]^{{a}}_{{{\bm{\Phi}}^{(T)}}}$ are real coefficients for ${a}\in\left\{0,1,\ldots,N_{T}-1\right\}$ (see Appendix C for detailed definitions), and the summation over the indices ${{\bm{\Phi}}^{(T)}}$ is implied. These $N_{T}$ independent EDM vectors $\left\{{\mathbf{R}}^{{a}}\right\}$ are sufficient to fully describe the non-Markovian dynamics in the finite time range $[0,T]$. ∎

Theorem 1 means that the ultimate complexity of OQS problems increases at most as a polynomial function of the evolution time. Note that the number of independent EDMs $N_{T}$ does not explicitly depend on the size of the bath. With increasing bath size, the energy range of the whole system (or, mathematically, the norm of the total Hamiltonian) may grow and therefore the time step $\epsilon$, which is inversely proportional to the norm of the Hamiltonian, needs to be reduced. This means an implicit dependence of $N_{T}$ on the bath size. However, such a dependence would be only polynomial (instead of exponential), since the total energy increases with the bath size only polynomially (actually, at most linearly for local interactions). Therefore, in principle, efficient solutions exist for general OQS problems.

III.2 Efficient MPS representation of OQS dynamics

The PC theorem upper-bounds the ultimate complexity of the OQS problems, but finding a set of independent EDM vectors $\left\{{\mathbf{R}}^{{a}}\right\}$ (and hence an explicit polynomial algorithm) is still challenging. Here, we propose an efficient tensor contraction algorithm to solve the EDM tensor. This algorithm is based on the matrix product state (MPS) representation of the EDMs $\rho^{{\bm{\Phi}}^{(t)}}$, in which the correlation tensors ${\mathbb{C}}^{{\bm{\Phi}}^{(t)}}_{{\bm{\Phi}}^{{}^{\prime}(t)}}$ (and hence the combined kernel tensors ${\mathbb{K}}^{\phi_{t};{\bm{\Phi}}^{(t-1)}}_{\phi^{\prime}_{t};{\bm{\Phi}}^{{}^{\prime}(t-1)}}$) are effectively represented by matrix product operator (MPO) tensor networks [64].

Theorem 2** (linearly increasing bond dimension theorem).**

The density matrix of an OQS has a time-local MPS representation as (see Fig. 3a)

[TABLE]

and in turn (with ${A}^{\phi_{t}}_{\phi^{\prime}_{t}}\equiv{A}_{\phi^{\prime\prime}_{t}}P^{\phi_{t},\phi^{\prime\prime}_{t}}_{\phi^{\prime}_{t}}$) the EDM

[TABLE]

and the MPS representation is efficient in the sense that the dimension $D_{\tau}$ (defined as the range of the lateral index $a_{\tau}$) of the bonds between the local tensors ${A}$’s increases at most linearly in time, that is

[TABLE]

Proof.

Figure 3b and 3c illustrate the proof. The EDMs can be obtained from the independent EDM vectors by ${\bm{\rho}}^{{\bm{\Phi}}^{(T)}}=\left[\bar{r}\right]^{{\bm{\Phi}}^{(T)}}_{a}{\mathbf{R}}^{{a}}$, where the inverse transform $\bar{r}$ is defined by $\left[r\right]_{{\bm{\Phi}}^{(T)}}^{a}\left[\bar{r}\right]^{{\bm{\Phi}}^{(T)}}_{{b}}=\delta^{{a}}_{{b}}$. Equation (11) holds for each component of the vectors. Therefore,

[TABLE]

This means the transforms $r$ and $\bar{r}$ can be viewed as isometry tensors as shown in Fig. 3b. Consider the density matrix $\rho(t)$, which can be obtained by applying the kernel tensor ${\mathbb{K}}_{{\bm{\Phi}}^{(t)}}\equiv{\mathbb{K}}^{0;{\mathbf{0}}}_{\phi_{t};{\bm{\Phi}}^{(t-1)}}$ to the EDM $\rho^{{\bm{\Phi}}^{(t-1)}}$, as shown in diagram 1 in Fig. 3c. By Theorem 1, the EDM can be reduced to $R_{t-1}^{a}$ using the isometry $r$ and then be recovered by $\bar{r}$ (diagram 2 in Fig. 3c), with the bond dimension between the two isometries $r$ and $\bar{r}$ (the range of $a$) being $D_{t-1}\leq N_{t-1}$. Then one can perform the singular value decomposition (SVD) of the combined tensor ${\mathbb{K}}\bar{r}$ as indicated by the vertical blue dashed line in the diagram, with the result ${\mathbb{K}}\bar{r}=U\Lambda V^{{\dagger}}$, in which $U$ and $V^{{\dagger}}$ are semi-unitary and the positive-definite diagonal matrix $\Lambda$ contains all non-zero singular values (see diagram 3) [65]. The number of non-zero singular values and hence the dimension of the bond between $\Lambda$ and $V^{{\dagger}}$ are at most $D_{t-1}$. The local tensor can be obtained as ${A}\equiv U\Lambda$ (diagram 4). The combined tensor $V^{{\dagger}}r{\mathbb{K}}$ applied to the EDM at $t-2$ (diagram 5) can be separated using the same procedure as that transforming diagram 2 to 3 to obtain another local tensor ${A}$ as shown in diagram 6. This procedure is continued until all local tensors ${A}$’s are obtained. Therefore, Eq. (12) is obtained (tensor diagram shown in upper row of Fig. 3a) and the upper limit of bond dimension in Eq. (14) is guaranteed. By applying the picking tensors (illustrated in lower row of Fig. 3a), Eq. (13) is also proved. ∎

Remarks:

As seen from the proof, the effects of a correlation tensor on the dynamics of the OQS are equivalent to those of a time-local MPO, that is

[TABLE]

and in turn

[TABLE]

The MPO representation dramatically reduces the complexity of the correlation tensors. Although the correlations of fluctuations in a many-body environment can be very complex (requiring $2^{{\mathcal{O}}(L)}$ free parameters to describe all $L$-th order correlations), they need to go through a very much limited number of channels to affect the dynamics of an OQS since the central system has a finite number of degrees of freedom, and therefore they can be organized into a small number of “modes” represented by an MPO. It should be noted that the equivalence in Eqs. (16) and (17) is only effective, that is, it does not mean that in general the correlation tensors have efficient MPO representation. 2. 2.

The combined kernel tensor ${\mathbb{K}}^{\phi_{t};{\bm{\Phi}}^{(t-1)}}_{\phi^{\prime}_{t};{\bm{\Phi}}^{{}^{\prime}(t-1)}}$ can be regarded as a special case of a correlation tensor (${\mathbb{K}}^{\phi_{t};{\bm{\Phi}}^{(t-1)}}_{\phi^{\prime}_{t};{\bm{\Phi}}^{{}^{\prime}(t-1)}}={\mathbb{C}}^{{\bm{\Phi}}^{(t)}}_{{\bm{\Phi}}^{{}^{\prime}(t)}}$ when ${\mathbb{C}}^{{\bm{\Phi}}^{(t-1)}}_{{\bm{\Phi}}^{{}^{\prime}(t-1)}}$ is an identity tensor), so it also has an efficient MPO representation when its effects on the OQS are considered. 3. 3.

The time-local MPO structure of the correlation tensors does not mean that the existence of a time-local equation of motion for the OQS dynamics, since the local tensors ${A}$’s need to be updated for each step of time evolution. Indeed, the evolution equation for EDMs in Eq. (8) is non-local in time.

Theorem 2 does not assume any constraint on the memory time in the non-Markovian process. When the bath has a finite memory time, the bond dimension in the MPS representation is reduced according to the following corollary.

Corollary 2.1**.**

When the bath has a finite memory time $T_{m}$, i.e., if all the culmulants $\tilde{C}_{\phi_{t_{l}},\ldots,\phi_{t_{2}},\phi_{t_{1}}}=0$ for $t_{l}-t_{1}>T_{m}$, the bond dimension between the local tensors ${A}$’s in Eqs. (12) and (13) is upper-bounded by

[TABLE]

This is because the memory time limits the number of independent EDM vectors involved in the evolution [see Appendix D for a proof].

Theorem 2 and Corollary 2.1 indicate that, once the cumulants of the noises from the bath are known (by computation or measurement) or if the OQS evolution has a tensor network representation, which are the case in a broad range of models (with examples discussed in Sec. IV), one can recursively solve the EDM tensor with the computation hour and memory expenditure increasing polynomially with the evolution time and the bath size. In the proof of Theorem 2, the MPS representation has been obtained for each time $t$ from the end to the beginning of the evolution. To save the computational time consumption, the local tensors ${A}$ can actually be recursively calculated from the beginning to the end of evolution as illustrated in Fig. 3d. The procedure is as follows. First, the calculation of ${A}$ tensor for $t=1$ is trivial. Then, if the ${A}$ tensors have been obtained for $t$, the EDM for $t+1$ is reached by applying the combined kernel function ${\mathbb{K}}^{\phi^{\prime}_{t+1};{\bm{\Phi}}^{{}^{\prime}(t)}}_{\phi_{t+1};{\bm{\Phi}}^{(t)}}{\mathcal{S}}^{\phi_{t+1}}$ (see diagram 1 in Fig. 3d). According to Theorem 2, the ${\mathbb{K}}$ tensor has an efficient MPO representation as shown in diagram 2 in Fig. 3d. The time-local ${A}$ tensors at the same time (enclosed by dashed hexagons in diagram 2 in Fig. 3d) can be combined into a single ${A}$ tensor as shown in diagram 3 in Fig. 3d. The ${A}$ tensors in diagram 3 in Fig. 3d can be reduced by transforming the tensor network into the right canonical form [65] and Corollary 2.1 guarantees that the bonds between them have dimensions bounded by Eq. (18). This procedure can be repeated until the EDM at a desired time $T$ is obtained.

Note that other numerical tensor-network methods have been proposed to deal with non-Markovian dynamics, such as the time-evolving MPO (TEMPO) [54] and the process tensor in the MPO (PT-MPO) form [57, 58, 68]. The TEMPO method is designed to deal with OQS’s in Gaussian bosonic environments. The PT-MPO method is formulated for general OQS’s by solving the process tensor, which in general, as noted in Ref. [59], is classically hard due to the potentially divergent bond dimension. In contrast, the method presented in this paper solves the EDM tensor, whose bond dimension is guaranteed to increase only linearly with evolution time.

IV Application of EDM tensor networks

IV.1 Spin-boson model

Here we use the spin-boson model [76] to demonstrate the efficiency of the EDM method. The Hamiltonian in the interaction picture reads

[TABLE]

where ${S}_{z}(t)=\cos\left(\mu t\right){S}_{z}+\sin\left(\mu t\right){S}_{y}$ describes the $z$-direction spin-1/2 operator under the tunneling term $\mu{S}_{x}$, ${a}^{\dagger}_{k}$ (${a}_{k}$) is the creation (annihilation) operator of the $k$-th bosonic mode with frequency $\omega_{k}$, and $g_{k}$ is the coupling strength between the spin and the $k$-th bosonic mode. The distribution of the coupling strengths is described by the spectral density $J(\omega)\equiv\sum_{k}g^{2}_{k}\delta(\omega-\omega_{k})$. In the numerical calculation, we set the initial state of the system to be $\rho(0)={S}_{z}+1/2$ (fully polarized along the $z$-axis) and the environment in the zero-temperature state (vacuum), i.e. $\Omega=\left|{\rm vacuum}\right\rangle\left\langle{\rm vacuum}\right|$ and take the Ohmic spectrum $J(\omega)=2{J_{0}}\omega e^{-\omega/\omega_{c}}$ with a dimensionless coupling constant ${J_{0}}$ and a cutoff frequency $\omega_{c}$. The quantum noise $B(t)$ has only second-order cumulants. From $J(\omega)$ we analytically obtain the combined kernel tensor and then numerically evaluate the EDM tensor using Eq. (8) and the tensor-network contraction procedure shown in Fig. 3d [see Appendix E for details].

The spin polarization $\left\langle S_{z}(t)\right\rangle=\text{Tr}\left[{S}_{z}(t)\rho(t)\right]$ is illustrated in Fig. 4a for various coupling strengths $J_{0}$. In the calculation, we set the time step as ${\epsilon}=0.01\mu^{-1}$ and introduce a self-adaptive truncation of the bond dimension by discarding the singular values $\{\lambda_{\alpha}\}$ that satisfy $\lambda_{\alpha}/\lambda_{d^{2}+1}\leq\xi$ with the truncation precision $\xi$ being a small number (e.g., $10^{-5}$). The spin relaxation dynamics exhibits a crossover from oscillation to overdamped decay at $J_{0}\approx 0.5$, indicating the competition between the coherent evolution and dissipation [76, 54]. Figure 4b shows the maximum bond dimension $D_{\rm max}\equiv\max\left(D_{\tau}\right)$ increasing with the evolution time $t$. In short times, $D_{\rm max}$ increases linearly for all coupling strengths, which means that the number of EDMs involved in the relaxation dynamics increases linearly. However, $D_{\rm max}$ saturates at later times, suggesting that a finite number of EDMs is sufficient for describing the non-Markovian dynamics. The maximum bond dimension increases faster and saturates at higher values for stronger coupling, which implies a dependence on the bath size since $J_{0}$ is proportional to the density of bosonic modes. However, the dependence of the bond dimension on $J_{0}$ (or the mode density) is sublinear (with a saturating behavior for strong coupling). These results validate the upper bound of the bond dimension given by Eq. (18) (dashed line in Fig. 4b). In addition, we verify that the linearly increasing bond dimension theorem holds regardless of the truncation precision [see Appendix E for details].

IV.2 Separable baths

IV.2.1 General formalism

There are cases where the quantum baths are separable, that is, can be decomposed into (approximately) independent sub-baths, as illustrated in Fig. 5a. Examples include a bath composed of independent bosonic modes [76], distinct spins in a spin bath [77, 78], and multi-bath OQS’s [79, 80]. In such cases, the noise fields are sums of those from different sub-baths, $B_{\alpha}(t)=\sum^{K}_{k=1}B_{k;\alpha}(t)$, with $K$ being the number of sub-baths. The bath being separable means that the noises from different sub-baths have no correlations, or mathematically, $\left[B_{k;\alpha}(t),B_{k^{\prime};\alpha^{\prime}}(t^{\prime})\right]=0$ for $k\neq k^{\prime}$ and $\Omega=\otimes_{k=1}^{K}\Omega_{k}$ with $\Omega_{k}$ being the initial density matrix of the $k$-th sub-bath. For such separable baths, the cumulants can be decomposed as

[TABLE]

where $\tilde{C}_{k;{\bm{\Theta}}}$ denotes the cumulant of the $k$-th sub-bath. This decomposition leads to a recursive formula of the EDMs [tensor diagrams illustrated in Fig. 5(b, c)]

[TABLE]

where $\rho^{{\bm{\Phi}}^{{}^{\prime\prime}(T)}}_{L}$ denotes the EDM tensor of the system coupled to the first $L$ sub-baths with $\rho^{{\bm{\Phi}}^{(T)}}_{0}={\mathcal{S}}^{{\bm{\Phi}}^{(T)}}\rho(0)$. Using the cumulants or the combined kernel functions of the $(L+1)$-th sub-bath ${\mathbb{K}}_{L+1;{\bm{\Phi}}^{(t)}}^{{\bm{\Phi}}^{{}^{\prime}(t)}}$ and the EDM $\rho^{{\bm{\Phi}}^{(t)}}_{L}$, one can calculate recursively the EDMs in a series of sub-baths following the same procedure as in Fig. 3d. Note that the linearly increasing bond dimension theorem holds for every EDM tensor $\rho^{{\bm{\Phi}}^{(T)}}_{L}$.

IV.2.2 Gaudin model

As a typical and important example of separable baths, we consider the Gaudin model [28], in which a central spin interacts with $K$ independent bath spins. The Hamiltonian reads

[TABLE]

where $g_{k}$ is the strength of the coupling between the central spin-1/2 ${{\mathbf{S}}}\equiv({S}_{x},{S}_{y},{S}_{z})$ and the $k$-th bath spin-1/2 ${{\mathbf{J}}}_{k}\equiv({J}_{k;x},{J}_{k;y},{J}_{k;z})$. The solution of the Gaudin model is in general hard. Solutions are available only for some special cases, such as reduction of the bath to single large spins when the coupling is uniform ($g_{k}={\rm constant}$), exact diagonalization for small $K$ [29], Bethe ansatz for highly polarized spin bath [81, 32], and density matrix renormalization group method for short time dynamics [30, 31].

An important feature of the Gaudin model is that the bath spins are time-independent, which leads to an infinite memory time. Furthermore, each bath spin constitutes a highly non-Gaussian quantum noise with an infinite number of high-order cumulants, which makes the cumulant expansion infeasible. Nevertheless, the decomposition procedures shown in Eq. (21) offer an efficient approach to solving the EDM tensor. The procedure is as follows (as illustrated in Fig. 5c). First, we compute the MPS EDM of the central spin $\rho^{{\bm{\Phi}}^{(T)}}_{L=1}$ with only the first bath spin considered. Then, using $\rho^{{\bm{\Phi}}^{(T)}}_{L=1}$ from the previous step, we compute the MPS EDM of the central spin including the effects of the first 2 bath spins, $\rho^{{\bm{\Phi}}^{(T)}}_{L=2}$. This is repeated until all bath spins are included and the MPS EDM $\rho^{{\bm{\Phi}}^{(T)}}=\rho^{{\bm{\Phi}}^{(T)}}_{L=K}$ is obtained. Each step is essentially an exactly solvable single-spin bath problem [see Appendix F for details].

In the numerical calculation, we assumed that the central spin is initially polarized along the $z$ axis, and that the bath spins are at an infinitely high temperature, that is, unpolarized. The couplings are set to have a linearly decreasing distribution $g_{k}=g\sqrt{\frac{6K}{2K^{2}+3K+1}}\frac{K+1-k}{K}$ [30]. The numerical results are shown in Fig. 6. In the calculation, the bond dimension $D_{t}$ is truncated by discarding the singular values $\{\lambda_{\alpha}\}$ that satisfy $\lambda_{\alpha}/\lambda_{d^{2}+1}\leq\xi$ with the truncation precision $\xi=10^{-6}$. In addition, a hard cutoff $D_{c}=400$ of the bond dimension is introduced to manage the computational time. The error introduced by this cutoff is analyzed in Appendix F. As shown in Fig. 6b, the bond dimension $D_{t}$ in the EDM tensor $\rho^{{\bm{\Phi}}^{(T)}}_{L}$ in the bath including the strongest $L$ spins increases linearly with time $t$ until it reaches the cutoff, which demonstrates the linearly increasing bond dimension theorem. It is interesting to note that the bond dimension increases nearly independently of the bath size as a function of the scaled time $\bar{g}_{L}t$, where the effective coupling strength $\bar{g}_{L}\equiv\sqrt{\sum^{L}_{k=1}g^{2}_{k}}$ for the $L$-spin bath (see Fig. 6b). Appendix F.4 presents more data to demonstrate the universal behavior of bond dimension increasing with the scaled time. This universal behavior means that for a given real time $t$ the number of independent equations needed to describe the central spin dynamics increases at most linearly with the bath size.

IV.2.3 Possible generalization to inseparable baths

The formalism presented here for computing EDMs in separable baths where the sub-baths have no interactions or initial correlations can be generalized to non-separable baths as long as the sub-baths can be arranged in a one-dimensional chain or tree-like structure and the interactions and initial correlations, if exist, are only short-ranged. The possibility of such generalization results from the chain-structure feature of Eq. (21) as shown in Fig. 5b. To generalize the EDM tensor network, one needs to modify the picking tensors connecting the sub-baths to be MPOs taking into account the short-range interactions or initial correlations. Remarkably, even when solving the chain-like quantum bath itself is computationally hard, the PC theorem guarantees that the OQS problem is solvable.

V Discussions and conclusions

In conclusion, the problem of non-Markovian dynamics of OQS’s in general is not computationally hard according to the polynomial complexity theorem. Furthermore, the problem can be formulated in terms of extended density matrix tensor networks. Guaranteed by the linearly increasing bond dimension theorem, the computational complexity of evaluating the EDM tensors increases only as a polynomial function of the evolution time and the bath size. We demonstrated the theorems and the method using two typical models, namely, the spin-boson model and the Gaudin model. The theory and numerical examples show that the OQS problems are solvable as long as the cumulants of noises in baths are known (by measurement or computation). Several numerical methods are available for computing the noise correlations and cumulants, such as the link-cluster expansion [82] and cluster-correlation expansion [83, 84, 85]. It is also possible to experimentally measure cumulants for many complex quantum baths through, e.g., quantum nonlinear spectroscopy [14, 15]. Another scenario of existence of explicit efficient algorithms for OQS problems is the cases where a two-dimensional tensor network can be constructed for the EDMs (including the temporal dimension for evolution), such as separable baths and one-dimensional baths with short-range interactions.

The significance of the PC theorem goes beyond the EDM tensor-network formalism, because the theorem guarantees the existence of an efficient solution of the OQS problem regardless of the complexity of the quantum many-body bath. The question is how to find an explicit algorithm. In the EDM formalism, the noise cumulants are needed, but the computation of the noise correlations in quantum many-body baths is in general a difficult problem [59]. However, in many cases, the OQS problem can be formulated without explicitly involving the noise correlations in the baths. For example, bath states can be formulated as tensor networks for a broad range of many-body systems [86, 66] and the reduced dynamics of the OQS can be obtained by tensor contraction, which, if executed in a proper procedure, is guaranteed by the PC theorem to have no exponentially growing bond dimension.

The PC theorem and the EDM formalism also offer new approaches to studying many-body systems. The evaluation of local observables in quantum many-body systems can be formulated as equivalent to an OQS problem. A direct consequence is that the complexity of many-body problems can be largely reduced according to the PC theorem when only the local observables are of interest.

Acknowledgements.

This work was supported by the Innovation Program for Quantum Science and Technology Project No. 2023ZD0300600, the Hong Kong Research Grants Council - General Research Fund (Project No. 14302121), the National Natural Science Foundation of China/Hong Kong Research Council Collaborative Research Scheme Project CRS-CUHK401/22, the Hong Kong Research Grants Council Senior Research Fellow Scheme Project SRFS2223-4S01, and the New Cornerstone Science Foundation.

Appendix A Tensor network representation of correlations

The correlation in Eq. (6) is in the form of summation of cumulant decomposition. Here we transform the correlation tensor into a tensor product of kernel functions $\tilde{K}_{{\bm{\Theta}}}\equiv\delta^{{\mathbf{0}}}_{{\bm{\Theta}}}+\tilde{C}_{{\bm{\Theta}}}$ (an example tensor network shown in Fig. 1d). Using the picking tensors $P^{\phi^{\prime},\phi^{\prime\prime}}_{\phi}$ and ${\mathbb{P}}^{{\bm{\Theta}}^{\prime},{\bm{\Theta}}^{\prime\prime}}_{{\bm{\Theta}}}\equiv\prod_{\phi_{\alpha}\in{\bm{\Theta}}}P^{\phi^{\prime}_{\alpha},\phi^{\prime\prime}_{\alpha}}_{\phi_{\alpha}}$, we define the open-armed correlation tensors and kernel tensors as

[TABLE]

where ${\bm{\Theta}}\subseteq{\bm{\Phi}}^{(t)}$ is an index subset and $\bar{\bm{\Theta}}\equiv{\bm{\Phi}}^{(t)}\setminus{\bm{\Theta}}$ is the complement. Note that the tensor $\tilde{K}^{{\bm{\Phi}}^{{}^{\prime}(t)}}_{{\bm{\Phi}}^{(t)}}\left({\bm{\Theta}}\right)$ simply passes the indices not involved in the cumulant $\tilde{C}_{\bm{\Theta}}$ to the next level. By the meaning of the picking tensor or direct calculation, one can show that the picking tensor is symmetric, i.e., $P_{\phi}^{\phi^{\prime},\phi^{\prime\prime}}=P_{\phi}^{\phi^{\prime\prime},\phi^{\prime}}$ and the product of two picking tensors is exchangeable, i.e.,

[TABLE]

which is illustrated in Fig. 7a. Therefore, the tensor product of two kernel tensors is exchangeable, i.e., $\left[\tilde{K}\left({\bm{\Theta}}_{1}\right)\tilde{K}\left({\bm{\Theta}}_{2}\right)\right]^{{\bm{\Phi}}^{{}^{\prime}(t)}}_{{\bm{\Phi}}^{(t)}}=\left[\tilde{K}\left({\bm{\Theta}}_{2}\right)\tilde{K}\left({\bm{\Theta}}_{1}\right)\right]^{{\bm{\Phi}}^{{}^{\prime}(t)}}_{{\bm{\Phi}}^{(t)}}$. The combined kernel tensor that contains all kernel functions with the first index being $\phi_{t+1}$ is defined in the tensor product form as

[TABLE]

which is illustrated as shadowed parts in Fig. 1d and e for $t=2$. To replace the summation form in Eq. (6) with the tensor product form, suffices it to demonstrate

[TABLE]

or equivalently [with Eq. (23a) and Eq. (6) applied to the l.h.s. and Eq. (25) and Eq. (24) recursively applied to the r.h.s.],

[TABLE]

where $\{{\bm{\Theta}}^{\prime\prime}_{1},{\bm{\Theta}}^{\prime\prime}_{2},\ldots,{\bm{\Theta}}^{\prime\prime}_{n}\}$ are all possible partitions of the subset of indices ${\bm{\Theta}}^{\prime\prime}$ and $\bar{\bm{\Theta}}^{\prime\prime}\equiv{\bm{\Phi}}^{{}^{\prime\prime}(t)}\setminus{\bm{\Theta}}^{\prime\prime}$ is the complement.

Let us consider the tensor product of two terms in r.h.s. of Eq. (26). Using Eq. (23b), we get

[TABLE]

When the two cumulants $\tilde{C}_{{\bm{\Theta}}_{u}}$ and $\tilde{C}_{{\bm{\Theta}}_{d}}$ have at least one overlapping index, e.g., $\phi_{\tau}$ (see Fig. 7b), the last term in Eq. (27) must vanish since $\delta^{\phi^{\prime\prime}_{\tau}}_{0}\tilde{C}_{{\bm{\Theta}}^{\prime\prime}_{u}}=0$. The physical meaning of this result is that noise at any time can affect the OQS only once through one cumulant (see Fig. 1c). When the two cumulants have no overlapping indices, i.e., ${\bm{\Theta}}_{u}\cap{\bm{\Theta}}_{d}=\varnothing$, the tensor product of the two kernel functions can be written as

[TABLE]

where ${\bm{\Theta}}_{u+d}\equiv{\bm{\Theta}}_{u}\bigcup{\bm{\Theta}}_{d}$ and $\bar{\bm{\Theta}}_{u+d}\equiv{\bm{\Phi}}^{(t)}\setminus{\bm{\Theta}}_{u+d}$. The tensor product with more kernel functions can be carried out similarly, and in the expansion, all products of cumulants with overlapping indices must vanish. Therefore, the r.h.s. of Eq. (26) can be written as

[TABLE]

where $\left\{{\bm{\Theta}}_{1},{\bm{\Theta}}_{2},\ldots,{\bm{\Theta}}_{k}\right\}$ is a non-overlapping partition of the indices ${\bm{\Theta}}$. The cumulant $\tilde{C}_{{\bm{\Theta}}_{i}}=0$ if any index $\phi\in{\bm{\Theta}}_{i}$ takes zero value, so ${\delta}^{{\bm{\Theta}}^{\prime}_{i}}_{\mathbf{0}}\tilde{C}_{{\bm{\Theta}}_{i}}={\mathbb{P}}^{{\bm{\Theta}}^{\prime}_{i},{\bm{\Theta}}^{\prime\prime}_{i}}_{{\bm{\Theta}}_{i}}\tilde{C}_{{\bm{\Theta}}^{\prime\prime}_{i}}.$ Inserting this formula and ${\delta}^{\bar{\bm{\Theta}}^{\prime}}_{\bar{\bm{\Theta}}}={\mathbb{P}}^{\bar{\bm{\Theta}}^{\prime},\bar{\bm{\Theta}}^{\prime\prime}}_{\bar{\bm{\Theta}}}\delta^{\mathbf{0}}_{\bar{\bm{\Theta}}^{\prime\prime}}$ into expression (29), we get the l.h.s. of Eq. (26). Thus, the correlation tensor is transformed into a tensor product of kernel functions as in the r.h.s. of Eq. (26).

Appendix B MPO representation of kernel tensors

Here we derive the MPO form of combined kernel tensors $\mathbb{K}^{\phi^{\prime}_{t+1};{\bm{\Phi}}^{{}^{\prime}(t)}}_{\phi_{t+1};{\bm{\Phi}}^{(t)}}$. Each kernel tensor can be written as a matrix product state (MPS) through a series of QR decomposition [65]

[TABLE]

where ${\bm{\Theta}}\equiv\left\{\phi_{\tau_{m}},\cdots,\phi_{\tau_{1}}\right\}$ and $\left\{{a}_{i}\right\}$ are left/right indices of the state tensors $Q$ and $R$.

We define the operator tensors $Q^{\phi^{\prime}}_{\phi}\equiv P^{\phi^{\prime},\phi^{\prime\prime}}_{\phi}Q_{\phi^{\prime\prime}}$, $R^{\phi^{\prime}}_{\phi}\equiv P^{\phi^{\prime},\phi^{\prime\prime}}_{\phi}R_{\phi^{\prime\prime}}$, and $\left[I^{\phi^{\prime}}_{\phi}\right]_{{a},{b}}\equiv\delta^{\phi^{\prime}}_{\phi}\delta_{{a},{b}}$, and obtain the MPO representation of kernel tensors as

[TABLE]

where the contraction of lateral indices $\{{a}_{i}\}$ between neighboring MPOs is implied. The combined kernel tensor $\mathbb{K}^{\phi^{\prime}_{t+1};{\bm{\Phi}}^{{}^{\prime}(t)}}_{\phi_{t+1};{\bm{\Phi}}^{(t)}}$ is then obtained by the contraction of all kernel tensors $\tilde{K}\left(\phi_{t+1};{\bm{\Theta}}\right)$ using Eq. (24). The MPO representation is preserved during the contraction of different kernel tensors through the following techniques

[TABLE]

[TABLE]

where the compact singular value decomposition $K_{{b}^{\prime},{b}}=\sum_{{a=1}}^{D}U_{{b}^{\prime},{a}}\lambda_{a}V^{\dagger}_{{a},{b}}\equiv U_{{b}^{\prime}}\Lambda\ V^{\dagger}_{{b}}$ is employed, the open arms in colors are indices kept intact (abbreviated as ${b}$ and ${b}^{\prime}$) in the decomposition, $U$ and $V$ are semi-unitary, i.e., $\left(U^{{\dagger}}U\right)_{{a}^{\prime},{a}}=\delta_{{a}^{\prime},{a}}$ and $\left(V^{{\dagger}}V\right)_{{a},{a}^{\prime}}=\delta_{{a},{a}^{\prime}}$, and $\lambda_{1}\geq\lambda_{2}\cdots\geq\lambda_{D}>0$ are all of the non-zero singular values with the bond dimension $D\leq\min\left(||{b}^{\prime}||,||{b}||\right)$ being the rank of the tensor $K$. In the calculation, the singular values that are smaller than a given truncation precision $\xi$ are neglected to reduce the bond dimension (bond dimension $D$ set by $\lambda_{D}/\lambda_{1}\geq\xi$). The contraction of these kernel tensors is executed from left to right, after which the result is restored to the QR representation as shown in Eq. (31) by QR decomposition of $\left[\tilde{K}(\bm{\Theta}_{u})\tilde{K}(\bm{\Theta}_{d})\right]^{\bm{0}}_{\bm{\Phi}^{{}^{\prime\prime}(t)}}$ and $\left[\tilde{K}(\bm{\Theta}_{u})\tilde{K}(\bm{\Theta}_{d})\right]^{\bm{\Phi}^{{}^{\prime}(t)}}_{\bm{\Phi}^{(t)}}\equiv{\mathbb{P}}^{{\bm{\Phi}}^{{}^{\prime}(t)},{\bm{\Phi}}^{{}^{\prime\prime}(t)}}_{{\bm{\Phi}}^{(t)}}\left[\tilde{K}(\bm{\Theta}_{u})\tilde{K}(\bm{\Theta}_{d})\right]^{\bm{0}}_{\bm{\Phi}^{{}^{\prime\prime}(t)}}$. Repeating the above contraction processes, we finally get

[TABLE]

Note that in the QR decomposition in Eq. (30) and hence in the MPO representation of kernel tensors in Eq. (31), the bond dimension of ${a}_{\tau_{i}}$ increases exponentially with $\min(i,t-i)$, making the QR decomposition of high-order tensors infeasible. Using Eq. (24), we can consider the contributions from different kernel tensors one by one and carry out the QR decomposition from right to left (the direction of time increasing) to construct the MPS representation of the EDMs. As guaranteed by the linearly increasing bond dimension theorem (Theorem 2), the bond dimension of ${a}_{i}$ in such MPS construction of EDMs will increase at most linearly with $i$. But the number of cumulants and hence that of the kernel tensors could still increase exponentially with the evolution time. However, for many types of quantum baths, only a limited number of low-order kernel tensors are involved. For example, Gaussian baths have only first- and second-order cumulants, $\tilde{C}_{\phi_{t+1}}$ and $\tilde{C}_{\phi_{t+1},\phi_{\tau}}$. Moreover, very often quantum baths have finite memory times, which upper bounds the number of cumulants needed to construct the tensor networks of EDMs.

Appendix C Progressive construction of independent EDM vectors

Formally, one can construct the independent EDM vectors step by step with time evolution. First, for $t=1$, we choose $\left[r\right]^{{a}}_{\phi_{1}}$ such that

[TABLE]

with ${a}=0,1,\ldots,N_{t=1}-1$. For $t=2$, we set $\left[r\right]^{{a}}_{\phi_{2},\phi_{1}}=\delta^{0}_{\phi_{2}}\left[r\right]^{{a}}_{\phi_{1}}$ for ${a}<N_{1}$ such that

[TABLE]

and construct $N_{2}-N_{1}$ new vectors

[TABLE]

Continuing this process, for the EDMs at time $t+1$ we can set $\left[r\right]^{{a}<N_{t}}_{\phi_{t+1},{\bm{\Phi}}^{(t)}}=\delta^{0}_{\phi_{t+1}}\left[r\right]^{{a}}_{{\bm{\Phi}}^{(t)}}$ such that

[TABLE]

and choose $N_{t+1}-N_{t}$ new independent vectors

[TABLE]

Thus we can formally construct the independent vectors step by step of the time evolution through the procedure of choosing coefficients

[TABLE]

Without loss of generality, we choose

[TABLE]

The EDMs can be obtained from the independent vectors through

[TABLE]

for any $t\leq T$, where the inverse coefficient matrices satisfy

[TABLE]

Appendix D Bond dimension for baths with finite memory times

Here we prove Corollary 2.1. The memory time $T_{m}$ of a bath is defined by the condition that the cumulants $\tilde{C}_{\phi_{t_{l}},\cdots,\phi_{t_{1}}}$ of the bath vanishes if $t_{l}-t_{1}>T_{m}$. From Eq. (24), one can get that the combined kernel tensors satisfy

[TABLE]

for $t>T_{m}$, where ${\bm{\Phi}}^{(t:t-T_{m}+1)}\equiv\left\{\phi_{t},\phi_{t-1},\ldots,\phi_{t-T_{m}+1}\right\}$. Therefore, the set of EMDs $\left\{\rho^{{\bm{\Phi}}^{(t:t-T_{m}+1)},{\mathbf{0}}}\equiv\rho^{{\bm{\Phi}}^{(t)}}\delta^{\mathbf{0}}_{{\bm{\Phi}}^{(t-T_{m})}}\right\}$ are sufficient for solving $\rho(t^{\prime})$ for $t^{\prime}>t$. All the EDM vectors ${\bm{\rho}}^{{\bm{\Phi}}^{(t:t-T_{m}+1)},{\mathbf{0}}}$ form a linear space with dimension $\leq N_{T_{m}}$ and hence can be expanded by at most $N_{T_{m}}$ independent vectors ${\mathbf{R}}^{{a}}={\left[{r}\right]}_{{\bm{\Phi}}^{(t:t-T_{m}+1)}}^{{a}}{\bm{\rho}}^{{\bm{\Phi}}^{(t:t-T_{m}+1)},{\mathbf{0}}}$ with ${a}=0,1,\ldots,N_{T_{m}}-1$, that is,

[TABLE]

The procedure for obtaining the MPS representation of the EDM tensor is illustrated in Fig. 8 by taking $T_{m}=2$ as an example. It follows the same SVD procedure as in Fig. 3c, except that the isometries $\bar{r}^{\bm{\Phi}^{(t)}}_{a}$ and $r^{a}_{\bm{\Phi}^{(t)}}$ stand for $\bar{r}^{\bm{\Phi}^{(t:t-T_{m}+1)},\bm{0}}_{a}$ and $r^{a}_{\bm{\Phi}^{(t:t-T_{m}+1)},\bm{0}}$, respectively, due to the finite memory time. The decomposition process shown in Fig. 8 is repeated until $t-T_{m}=0$, after which the same procedure as in Fig. 3c is followed and the bond dimension is upper bounded by $N_{t}$. Consequently, the bond dimension between local tensors defined in Eq. (12) and Eq. (13) is upper-bounded by $D_{\tau}\leq\min(N_{\tau},N_{T_{m}})$.

Appendix E Spin-boson model

The bath is Gaussian. At zero temperature, the noise ${B}(t)=\sum_{k}g_{k}\left({a}^{\dagger}_{k}e^{i\omega_{k}t}+{\rm h.c.}\right)$ has only the second-order cumulants being non-zero, which are

[TABLE]

where $f(\tau)\equiv\int d\omega J(\omega)e^{-i\omega\tau}$, with $\mathcal{B}_{\phi_{t}=1}\equiv\mathcal{B}^{-}(t)$ and $\mathcal{B}_{\phi_{t}=2}\equiv\mathcal{B}^{+}(t)$. The kernel tensor has QR decomposition $\tilde{K}_{\phi_{t+1},\phi_{\tau+1}}=\left[Q_{\phi_{t+1}}\right]_{{a}}\left[R_{\phi_{\tau+1}}\right]_{a}$ with lateral bond dimension $D_{a}=2$ and the combined kernel tensor

[TABLE]

where $Q^{\phi^{\prime}}_{\phi}=P^{\phi^{\prime},\phi^{\prime\prime}}_{\phi}Q_{\phi^{\prime\prime}}$ and $R^{\phi^{\prime}}_{\phi}=P^{\phi^{\prime},\phi^{\prime\prime}}_{\phi}R_{\phi^{\prime\prime}}$ with

[TABLE]

[TABLE]

and

[TABLE]

for $0<\tau<t$.

In programming, we introduce super-operator-valued local tensors ${\mathcal{L}}^{\phi_{\tau}}=A^{\phi_{\tau}}_{\phi^{\prime}_{\tau}}{\mathcal{S}}^{\phi^{\prime}_{\tau}}$ for MPS representation of the EDMs shown in Eq. (13), which results in

[TABLE]

The diagram illustration is shown in Fig. 9a. The standard MPO-MPS algorithms are used to contract the combined kernel tensor and the EDM tensor to evaluate the EDMs [65]. Such a process is shown in Fig. 9b. With the contraction carried out from the left boundary to the right boundary, for example, the following SVD technique is used

[TABLE]

where in the approximation the singular values that are too small ($\lambda_{a}<\xi\lambda_{d^{2}+1}$) are dropped with $\xi$ being the truncation precision.

In the numerical simulation, to reduce the error due to the discretization of time, we expand a small step of evolution to the second order of time step $\epsilon$,

[TABLE]

which can be written as the tensor product form

[TABLE]

where $\mathcal{S}_{1}^{\phi_{t}=0}\equiv\mathcal{S}_{2}^{\phi_{t}=0}\equiv\mathcal{I}_{\rm S}$, $\mathcal{S}^{\phi_{t}\neq 0}_{1}\equiv\frac{1-i}{2}\mathcal{S}^{\phi_{t}}$, and $\mathcal{S}^{\phi_{t}\neq 0}_{2}\equiv\frac{1+i}{2}\mathcal{S}^{\phi_{t}}$. Using this formula, we generalize the tensor network algorithm developed in the main text to the second order of $\epsilon$, with each term $\mathcal{B}_{\phi_{t}}\mathcal{S}^{\phi_{t}}$ replaced with $\mathcal{B}_{\phi_{t}}\mathcal{B}_{\psi_{t}}\mathcal{S}^{\phi_{t}}_{2}\mathcal{S}^{\psi_{t}}_{1}$. It is straightforward to generalize the second-order expansion of $\epsilon$ to higher orders.

Figure 10 compares the calculated spin dynamics and the maximal bond dimension $D_{\rm\max}$ needed to evaluate $\rho^{\bm{\Phi}^{(t)}}$ for truncation precision $\xi=10^{-5}$ (the same as in Fig. 4 in the main text) and $\xi=10^{-6}$. The calculation of the spin dynamics for truncation precision $\xi=10^{-5}$ is sufficiently accurate. As expected, smaller truncation precision induces a faster increase of the bond dimension. However, for all coupling strengths and both truncation precisions, the increase is not faster than linear.

Appendix F Gaudin model

F.1 MPS representation of EDMs

As a feature of non-interacting spin-1/2’s, the correlation tensor for the $k$-th bath spin at infinite temperature can be decomposed into a product of local tensors as

[TABLE]

where $\left[{A}_{k;\phi_{t}}\right]_{a_{t},a_{t-1}}\equiv\frac{1}{2}\text{Tr}\left[\sigma_{a_{t}}\mathcal{B}_{k;\phi_{t}}\sigma_{a_{t-1}}\right]$ with $B_{k;1/2/3}\equiv g_{k}{J}_{k;x/y/z}$, $\sigma_{0}\equiv I$ and $\sigma_{1/2/3}\equiv 2J_{k;x/y/z}$ with the lateral bond dimension $D_{a}=4$ for $a_{\tau}$. The EDM of the central spin coupled to the first spin is

[TABLE]

whose MPS can be obtained with the SVD described in Appendix E; then the MPS representation of the EDM coupled to the first $L+1$ bath spins can be recursively obtained using Eq. (21).

F.2 Central spin dynamics

The dynamics of the central spin at any given time $t$ can be obtained from

[TABLE]

since

[TABLE]

where ${\bm{\Phi}}^{(T)}/\{\phi_{t}\}$ is the complement of $\{\phi_{t}\}$ and Eq. (9) has been used in the last step.

In numerical calculation, we expand each step of evolution to the second order using Eq. (E). For this second-order expansion, the EDM tensor is generalized to

[TABLE]

where $\mathbb{C}_{\bm{\Phi}^{(T)},\bm{\Psi}^{(T)}}\equiv\text{Tr}_{B}\left[\mathcal{B}_{\phi_{T}}(T)\mathcal{B}_{\psi_{T}}(T)\cdots\mathcal{B}_{\phi_{1}}(1)\mathcal{B}_{\psi_{1}}(1)\Omega\right]$ and $\mathbb{S}^{\bm{\Phi}^{(T)}}_{i}\equiv\mathcal{S}^{\phi_{T}}_{i}(T)\cdots\mathcal{S}^{\phi_{1}}_{i}(1)$ for $i=1$ or $2$. Then, the evaluation of spin dynamics in Eq. (44) is modified to be

[TABLE]

To obtain this expression, we expand small steps of evolution $e^{\epsilon{\mathcal{H}}^{-}(\tau)}$ to the second order using Eq. (42) and get

[TABLE]

where the last equation above is obtained by inserting tensor products like ${\mathcal{B}}_{\phi}{\mathcal{S}}_{i}^{\phi}$ in the beginning of the expression within the trace operation since trace of commutators always vanishes. Using the substitution $\mathcal{S}^{2{\alpha}-1}_{1}=P^{\psi_{t+1},2{\alpha}-1}_{\psi^{\prime}_{t+1}}{\mathcal{B}}_{\psi_{t+1}}\mathcal{S}^{\psi^{\prime}_{t+1}}_{1}$ for $2\alpha-1>0$, we transfer Eq. (F.2) to

[TABLE]

which leads to Eq. (45). Note that in the model considered here, we have $S_{\alpha}(t)=S_{\alpha}(t+1)$. In general, there could be evolution of the system operators due to dynamical control. In such cases, one just need to modify the formula by writing $S_{\alpha}(t)$ as a linear combination of $S_{\beta}(t+1)$.

F.3 Error induced by bond dimension cutoff

In the calculation, we introduce a cutoff of the bond dimension $D_{c}$ to manage the computation cost, which would induce error in addition to the truncation precision $\xi$ in the SVD and the expansion error $\mathcal{O}(\epsilon^{3})$ in Eq. (E). Figure 11 compares the results obtained with various bond dimension cutoff, with the exact numerical results also shown for small numbers of bath spins ($L\leq 13$). One can find that when $L=5$, different cutoffs yield almost the same decoherence dynamics, all in good agreement with the exact solution, which indicates that the non-Markovian dynamics in this case can be described by a small number of EDMs. However, for $L=13$ bath spins considered, the deviation from the exact solution is noticeable when the bond dimension cutoff is small $D_{c}\leq 300$. The deviation decreases with increasing the cutoff $D_{c}$. For larger spin baths ($L\geq 25$), no exact solution is available. However, the convergent dynamics indicates that the EDM calculation with bond dimension cutoff $D_{c}=400$ (adopted in the main text) is reliable up to errors $\sim 10^{-3}$.

F.4 Universal behavior of bond dimension increase

In the main text, we have seen that the bond dimension $D_{t}$ in different EDM tensors $\rho^{\bm{\Phi}^{(T)}}_{L}$ as functions of the rescaled evolution time $\bar{g}_{L}t$ are nearly identical for different $L$. Figure 12 shows more data. The increase of the bond dimension collapses to a universal curve until it meets the cutoff (set as $D_{c}=400$) or saturates for small baths ($L\leq 7$). This phenomenon demonstrates that with a given coupling strength between the system and the bath, the growth of bond dimension is independent of the bath size. As the number of bath spins increases, the degrees of freedom associated with the bath spins increase exponentially. However, the growth of the bond dimension in EDMs is nearly linear in time, which confirms the linearly increasing bond dimension theorem.