Direct probing of the simulation complexity of open quantum many-body dynamics

TL;DR

This paper investigates the computational complexity of simulating open quantum many-body systems, revealing how dissipation influences correlation length, mixing time, and the disparity between quantum and classical simulation resources.

Contribution

It provides a detailed analysis of how dissipation affects the complexity of open quantum system simulations, highlighting differences between quantum and classical methods.

Findings

Dissipation impacts correlation length and mixing time differently over time.

Classical tensor network complexity does not decrease with stronger dissipation.

Quantum and classical resource scaling diverge in the presence of dissipation.

Abstract

Simulating open quantum systems is key to understanding non-equilibrium processes, as persistent influence from the environment induces dissipation and can give rise to steady-state phase transitions. A common strategy is to embed the system-environment into a larger unitary framework, but this obscures the intrinsic complexity of the reduced system dynamics. Here, we investigate the computational complexity of simulating open quantum systems, focusing on two physically relevant parameters -- correlation length and mixing time -- and explore whether it can be comparable (or even lower) to that of simulating their closed counterparts. In particular, we study the role of dissipation in simulating open-system dynamics using both quantum and classical methods, where the classical complexity is characterised by the bond dimension and operator entanglement entropy. Our results show that dissipation affects correlation length and mixing time in distinct ways at intermediate and long timescales. Moreover, we observe numerically that in classical tensor network simulations, classical complexity does not decrease with stronger dissipation, revealing a separation between quantum and classical resource scaling.

Figures

Equations

$$\mathcal{L}[\rho ]=-i[H,\rho ]+\gamma \mathcal{N}[\rho ],$$

$$\mathcal{N}[\rho ]=$$

$$+(1-{\eta }_k)L_k\rho L_k^{†}+{\eta }_k{⟨L_k{L}_k^{†}⟩}_{\rho }\rho ),$$

$$K_0=|0⟩⟨0{|}_R\otimes \left(\mathbb{I}-iH\delta t-\frac{1}{2}\sum _k\gamma {\alpha }_k{L}_k^{†}{L}_k\delta t\right),$$

$$K_{k,\mathrm{jp}}=|k,\mathrm{jp}⟩⟨0{|}_R\otimes \sqrt{\gamma {\alpha }_k(1-{\eta }_k)\delta t}{L}_k,$$

$$K_{k,\mathrm{disc}}=|k,\mathrm{disc}⟩⟨0{|}_R\otimes \sqrt{\gamma {\alpha }_k{\eta }_k\delta t}{L}_k.$$

$${\rho }^{\mathrm{un}}(t+\delta t)$$

$$+\gamma \sum _k(1-{\eta }_k){\mathcal{D}}_k[\rho ])+\mathcal{O}(\delta t^2),$$

$$C^{\mathrm{env}}(t)=\sum _{k=1}^K{d}_k{e}^{-z_{k}t},$$

$${\hat{a}}_k^{†}|\mathbf{m}⟩$$

$${\hat{b}}_k^{†}|\mathbf{n}⟩$$

$${\hat{\mathcal{L}}}_S[{\rho }_{\mathcal{S}}]=-i\left[{\hat{H}}_{\mathrm{pm}},{\rho }_{\mathcal{S}}\right]+\gamma {\mathcal{D}}_{\mathrm{pm}}[{\rho }_{\mathcal{S}}],$$

$${\hat{H}}_{\mathrm{pm}}=H+\sum _k\left[\mathrm{Im}(z_k){\hat{a}}_k^{†}{\hat{a}}_k-\sqrt{\mathrm{Re}(d_k)}{q}_S\left({\hat{a}}_k^{†}+{\hat{a}}_k\right)\right],$$

$$Q(t+\tau )=\underset{\delta t\to 0}{\lim }\prod _j^{\tau /\delta t}(1-P_{\text{jp}}(t+j\delta t))=e^{-{\int }_t^{t+\tau }{P}_{\text{jp}}({\tau }^{\prime })d{\tau }^{\prime }}.$$

$$H_{\text{eff}}=H+\frac{i\gamma }{2}\sum _k\left(L_k^{†}{L}_k-⟨L_k^{†}{L}_k⟩\right),$$

$$\mathrm{Cor}(\hat{A},\hat{B})\le \Vert \hat{A}\Vert \Vert \hat{B}\Vert e^{-\frac{d}{\xi }}$$

$$\mathrm{Cor}(\hat{A},\hat{B}):=⟨{\Psi }_S|\hat{A}\otimes \hat{B}|{\Psi }_S⟩-⟨{\Psi }_S|\hat{A}\otimes \mathbb{I}|{\Psi }_S⟩⟨{\Psi }_S|\mathbb{I}\otimes \hat{B}|{\Psi }_S⟩,$$

$$C_{A,B}(t):=\inf \left\{C>0:\left|\mathrm{Cor}(\hat{A},\hat{B})\right|\le \Vert \hat{A}\Vert \Vert \hat{B}\Vert e^{-\frac{d}{C}}\right\},$$

$$\overline{C(d,t)}:=\frac{1}{|{\mathcal{P}}_d|}\sum _{(A,B)\in {\mathcal{P}}_d}{C}_{A,B}(t).$$

$$\overline{\mathrm{Cor}(d,t)}\le \overline{e^{-d/C_{A_i,B_i}}}{|}_{A_i-B_i=d}\le e^{-d/\bar{C}}$$

$$H_I=J\sum _i{X}_i{X}_{i+1}+h_x\sum _i{X}_i+h_z\sum _i{Z}_i,$$

$$L_1$$

$$L_2$$

$$L_3$$

$$L_4$$

$$c_{ij}=\mathrm{Tr}({\sigma }_i^z{\sigma }_j^z\rho )-\mathrm{Tr}({\sigma }_i^z\rho )\mathrm{Tr}({\sigma }_j^z\rho ),$$

$$\epsilon ={\left(\frac{{\sum }_{ij}{\left(c_{ij}-{\tilde{c}}_{ij}\right)}^2}{{\sum }_{ij}{\tilde{c}}_{ij}^2}\right)}^{1/2},$$

$$\hat{\rho }=\sum _i{\lambda }_i{\hat{o}}_i^{[L]}{\hat{o}}_i^{[R]},$$

$$|\rho ⟩⟩=\sum _{i,j}{\rho }_{ij}|i⟩\otimes |j⟩.$$

$$\frac{d}{dt}|\rho ⟩⟩=\hat{\mathcal{L}}|\rho ⟩⟩,$$

$$\hat{\mathcal{L}}=$$

$$+\sum _k\left(L_k\otimes L_k-\frac{1}{2}\mathbb{I}\otimes L_k^{†}{L}_k-\frac{1}{2}{L}_k^{†}{L}_k\otimes \mathbb{I}\right).$$

$$\hat{\mathcal{L}}=-i({\lambda }_1{H}_1-i{\lambda }_2{H}_2),$$

$$|\rho (t)⟩⟩=\exp [-i({\lambda }_1{H}_{1}t-i{\lambda }_2{H}_{2}t)]|\rho (0)⟩⟩.$$

$$|\rho (t)⟩⟩=$$

Full text

I Introduction

In a way, quantum systems are in constant contact with a noisy environment [1, 2]. Such openness gives rise to processes absent in closed systems, like decoherence and dissipation, which fundamentally shape the system dynamics and often drive irreversibility. Open quantum dynamics has therefore gained considerable attention due to its relevance in understanding non-equilibrium processes that are pervasive in real-world scenarios across biology [3, 4], materials science [5], thermodynamics [6, 7], and control [8, 9]. Simulating such open quantum systems is generally hard [10], as the system-environment interactions lead to nonunitary dynamics, driven by dissipative effects like relaxation and thermalisation [11, 12, 13, 14]. A common way to deal with nonunitaries is by embedding into a unitary evolution (on an extended space), such as in Stinespring dilation [15, 16, 17], purification [18], unitary decomposition [19, 20], or in a more general context, by adding ancilla qubits [21, 22, 23, 24, 25, 26, 27]. In this way, the nonunitary dynamics of the reduced system is recovered by tracing out the environmental degrees of freedom. This idea has been instrumental in the algorithmic design of open quantum systems, enabling, for instance, the modelling of noisy quantum channels [28] and realisations based on a linear combination of unitaries or channels [23, 24, 25, 26, 29, 30, 31, 32].

Nevertheless, the embedding inherent to dilation often obscures the intrinsic complexity of the reduced system evolution [33], making it hard to attribute observed complexity to the structure of the open quantum system itself, rather than to artifacts of the dilation procedure. Meanwhile, the role of dissipation in the simulation complexity of open quantum dynamics has attracted interest on its own due to its nuanced nature. Dissipation can in fact lead to an increase in complexity in various scenarios, which is a bit counter-intuitive, as it is typically linked to the suppression of coherence and entanglement [34, 35, 36, 37]. As a concrete example, Prosen and Pizorn [38] showed that for an open Heisenberg XY spin chain, quantum phase transitions occur with small magnetic field, which exhibit long-range correlations in the non-equilibrium steady state and linear growth of the operator space entanglement entropy. Similarly, dephasing may lead to surprising behaviour at different time scales; specifically, the operator entanglement, while expected to be suppressed by dephasing, may rise, fall and rise again under time evolution [39]. From a complexity point of view, Aloisio et al. [33] studied the complexity of sampling an open quantum system and demonstrated the existence of a stochastic quantum process that makes the multitime sampling as complex as sampling from a many-body state, which is classically hard. Therefore, a question arising from this is: Can we directly probe the complexity of an open quantum system and the contribution due to dissipation, via physically relevant parameters, without resorting to dilation.

In this work, we study this question by exploring the complexity of simulating open quantum systems, focusing on physically meaningful parameters—correlation length and mixing time—in both short- and long-time regimes [40, 41, 42, 43, 44], as suggested in the context of quantum imaginary-time evolution (QITE) algorithms [45]. Specifically, we employ two complementary approaches: (i) classical tensor network (TN) simulations, and (ii) quantum simulation algorithms based on unitary approximation of nonunitary operators [45] (see subsequent development in [46, 47, 48, 49, 50, 51]). We are particularly interested in how these spatial and temporal resources scale with respect to dissipation strength, time, and initial conditions of the system. Here, we consider a general open quantum system setup. We start with Markovian dynamics, including both time-dependent [52] and nonlinear Lindbladian master equations [53] which can interpolate to non-Hermitian dynamics [54]. Then, we discuss the simulation of non-Markovian dynamics [55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67], in which we are focused on unravelling the influence functional of the environment with complex-valued paths proposed in [67] and simulation in an extended bosonic Hilbert space. To simulate the (non-)Markovian dynamics, we develop a stochastic simulation method based on quantum jump formulation [68]. We analyse how dissipation affects the simulation complexity, where our technical contribution is linking it to the average correlation length, which is region- and time-dependent, and strictly smaller than the correlation length that is commonly used in literature (e.g., in [45]).

We show both analytically and numerically that increasing dissipation reduces the correlation length over time for quantum spin systems. This indicates that the dynamics of open quantum spins may be easier to simulate at long times than their closed counterpart. In the quantum case, we simulate the evolution of local spin observables in a dissipative spin chain with transversal and longitudinal fields up to 25 spins. The classical simulation complexity is characterised by the bond dimension and operator entanglement entropy, which quantify the growth of correlations and the cost of representing the system state efficiently. In contrast to the quantum case, we observe numerically that the complexity of classical TN simulations remains non-decreased with increasing dissipation strength, thereby showing a separation between quantum and classical scaling. This suggests the existence of regions, where quantum methods may be easier to find advantages over classical techniques in dissipative chaotic spin systems. Finally, we consider an extension to ground-state preparation under time-dependent Lindbladians, where the jump operators are applied randomly at different sites over time in contrast to the dissipative methods in [69, 70].

II Results

II.1 Formulation of open quantum dynamics

Many quantum processes can be captured by quantum Markov dynamics, where the evolution is recast into a quantum dynamical semigroup that generalises the Schrödinger equation. We consider a quantum system $S$, described by Hamiltonian $H$ and density operator $\rho$, which are Hermitian in the operator space $\mathcal{A}(\mathcal{H})$ with Hilbert space $\mathcal{H}_{S}$; and suppose the evolution of $S$ is monitored by a detector with efficiency $\eta_{k}\in[0,1]$, see [53].

In particular, a spin system coupled linearly to a bosonic environment. The total Hamiltonian is defined on the composite Hilbert space $\mathcal{H}_{\mathrm{tot}}=\mathcal{H}_{S}\otimes\mathcal{H}_{B}$, where $\mathcal{H}_{B}$ is the Hilbert space of the bosonic environment (or called bath). We consider both Markovian and non-Markovian regimes, where formally (in the weak coupling limit) the dynamics of $S$ is governed by a Lindblad master equation $\partial_{t}\rho=\mathcal{L}[\rho]$, with Lindbladian given by

[TABLE]

which consists of both unitary and dissipative components, denoted by $-i[H,\rho]$ and $\gamma\mathcal{N}[\rho]$, respectively, and $\gamma\geq 0$ is the interaction strength which dictates the rate of dissipation. The detector efficiency implies a nonlinear term in the Lindbladian. Explicitly

[TABLE]

where $\|\mathcal{N}\|=1$. When $\eta_{k}=0$, it reduces to the standard Lindbladian master equation, while for $\eta_{k}=1$, it describes evolution under an effective non-Hermitian Hamiltonian [54]. The set of jump operators $\left\{L_{k}\right\}\subset\mathcal{A}(\mathcal{H})$ describe the dissipation in $S$, i.e., the decay of the off-diagonal terms of the positive semidefinite matrix $\rho$, while $\left\{\alpha_{k}\right\}\in\mathbb{R}$ are the corresponding decay coefficients; in the Markovian regime $\alpha_{k}>0$. Here, each term in Eq. 1 is time-dependent and the coefficients are omitted for notational simplicity. Unlike the linear case, this Lindbladian can capture non-Hermitian dynamics. From a computational perspective, its formulation does not allow for straightforward vectorisation.

This nonlinear process can be described in several ways. One way is to formalise it within a quantum instrument framework as $\rho(t+dt)=\mathcal{E}(\rho(t))$, where $\mathcal{E}=\sum_{k}\mathcal{E}_{k}$ denotes a family of dynamical maps [53]. Here, we consider an extended Hilbert space $S\otimes R$, where $R$ stores the classical record of the detector with orthogonal labels $\{\ket{0},|k,\textrm{jp}\rangle,|k,\textrm{disc}\rangle\}_{k}$, denoting the three possible outcomes at step $k$: no jump, a jump takes place but remains undetected, and a jump is detected but is then discarded, respectively. Each of the three outcomes is associated with a corresponding operator acting on $S\otimes R$. Specifically

[TABLE]

It is easy to check that the reduced map on the system space $\Lambda(\rho)=\operatorname{Tr}_{R}(\sum_{k}K_{k}(\rho\otimes|0\rangle\langle 0|)K_{k}^{\dagger})$ is completely positive trace-preserving (CPTP) in the first order of $\delta t$. Given $S$ is in state $\rho$ at time $t$, the unnormalised post-selected state (labelled by $r\in\{0\}\cup\{(k,\textrm{jp})\}$) at $t+\delta t$ is

[TABLE]

where $\mathcal{D}_{k}[\rho]=\alpha_{k}(L_{k}\rho L_{k}^{\dagger}-\frac{1}{2}\{L_{k}^{\dagger}L_{k},\rho\})$ is the dissipator. The corresponding success probability is then $p_{\text{keep }}=1-\sum_{k}\gamma\alpha_{k}\eta_{k}\delta t\langle L_{k}^{\dagger}L_{k}\rangle+\mathcal{O}(\delta t^{2})$. Normalising $\rho^{\text{un}}$ yields the evolution given by Eq. 2, where the last term is exactly the nonlinear drift that appears solely due to the renormalisation after discarding observed-click runs.

Another way to describe Eq. 2 makes use of the fact that it admits a natural interpretation in terms of quantum jumps [71, 72, 73]. We again use a detector with efficiency $\eta_{k}$ to monitor the occurrence of quantum jumps, where we suppose that a fraction of the data (given by $\eta_{k}$) is discarded as post-processing. The role of $\eta_{k}$ (and the post-selection) here is to interpolate the evolution between the normal dissipative and non-Hermitian dynamics. Considering the normalisation, a jump happens with probability $\frac{(1-\eta_{k})\delta p_{k}}{1-\eta_{k}\delta p_{k}}\approx(1-\eta_{k})\delta p_{k}$, where $\delta p_{k}(t):=\delta t\alpha_{k}\operatorname{Tr}(L_{k}^{\dagger}L_{k}\rho)$, whereas a “no jump” event occurs with probability $\frac{1-\delta p_{k}}{1-\eta_{k}\delta p_{k}}\approx 1-(1-\eta_{k})\delta p_{k}$. From here, it is clear that this can be described by a quantum jump process. A system $S$ in pure state $\ket{\psi(t)}$ evolves for each small timestep $\delta t$ according to the stochastic dynamics: with probability $(1-\eta_{k})\delta p_{k}(t)$ the state jumps to $L_{k}\ket{\psi(t)}/{\|L_{k}\ket{\psi(t)}\|}$, while with probability $1-(1-\eta_{k})\delta p_{k}(t)$ it evolves to $e^{-\frac{1}{2}L_{k}^{\dagger}L_{k}\delta t}{\ket{\psi(t)}}/{\|e^{-\frac{1}{2}L_{k}^{\dagger}L_{k}\delta t}\ket{\psi(t)}\|}$, which leads to Eq. 2 [54].

The above discussion uses a master equation to describe the dynamics of $S$. In general, the system-bath coupling induces finite environmental memory, and the dynamics cannot be captured by time-local Lindbladian operators alone. Let us consider the system-bath coupling $H_{\rm int}={q}_{S}\otimes{q}_{B}$, with bath operator ${q}_{B}\in\mathcal{H}_{B}$ and a free Gaussian bath. There are various approaches to unravel the influence of the bath, such as using pseudomodes [55, 74, 56, 75, 59], hierarchy of the equations of motion (HEOM) [57, 66, 76], dissipaton equation of motion [65, 77], dissipaton-embedded quantum master equation (DQME) [64], and quantum dissipation with minimal state space (QD-MESS) [67]. The simulation of non-Markovian dynamics on a quantum computer has also been discussed [63, 78, 62]. A common approach is to characterise the environment by the bath correlation function $C^{\rm env}(t)=\braket{{q}_{B}(t){q}_{B}(0)}$, and then decompose it into a finite number of complex exponential terms

[TABLE]

where we choose to constrain the amplitude to be real $d_{k}\in\mathbb{R}^{1}$, while $z_{k}\in\mathbb{C}^{1}$ is a complex number. In this way, one can represent the dynamics of states in an extended space $\Gamma=\mathbb{L}_{S}\otimes\mathbb{F}_{2K}$, where $\mathbb{L}_{S}$ is the Liouville space of $S$, and $\mathbb{F}_{2K}$ is the Fock space of $2K$-mode harmonic system. In Supplementary Section S1.2, we discuss the unravelling of the influence functional based on the extended state space introduced in [67] and how to map to a quasi-Lindblad equation.

Formally, one can introduce a generator $\hat{\Gamma}_{k}=-z_{k}\hat{a}_{k}^{\dagger}\hat{a}_{k}-z_{k}^{*}\hat{b}_{k}^{\dagger}\hat{b}_{k}$ to describe the dynamics of the unperturbed auxiliary bosons in $\mathbb{F}_{2K}$, with $\hat{a}_{k}$ and $\hat{b}_{k}$ defined, respectively, as

[TABLE]

where $\mathbf{m}=\{m_{1},...,m_{k},...,m_{K}\}$ and $\mathbf{n}=\{n_{1},...,n_{k},...n_{K}\}$ label the HEOM $\rho_{\mathbf{m},\mathbf{n}}$. They are the quantum occupation numbers that are aligned with the multi-index in free-pole HEOM to label the auxiliary boson states, where $\ket{\mathbf{n}}$ denotes the pseudo-Fock states. In the main text, the hat notation is used to denote the terms in the extended space. One may then apply a similarity transformation $\mathcal{S}$ with $\det(\mathcal{S})=1$ in Fock space that mixes $\hat{a}_{k}$ and $\hat{b}_{k}$ to map the time evolution of the density matrix to a Lindblad form when the amplitude is constrained to the real domain. Now to make it compatible for implementation on a quantum computer, one may further define the state by only using half of the auxiliary bosons that reside in a dual Fock space $\{\ket{\mathbf{m}}\bra{\mathbf{n}}\}$. In doing so, it was shown in [67] that the dynamics is mapped to a Lindblad form $\partial_{t}{\rho}_{\mathcal{S}}=\hat{\mathcal{L}}_{S}[\rho_{\mathcal{S}}]$, where

[TABLE]

and the pseudomode Hamiltonian takes the form

[TABLE]

with $\mathcal{D}_{\rm pm}=\sum_{k=1}^{K}\alpha_{k}\left[\hat{a}_{k}\rho_{\mathcal{S}}\hat{a}_{k}^{\dagger}-\frac{1}{2}\left\{\hat{a}_{k}^{\dagger}\hat{a}_{k},\rho_{\mathcal{S}}\right\}\right]$, $\gamma=2\sum_{k=1}^{K}|\operatorname{Re}(z_{k})|$, and $\alpha_{k}=2\operatorname{Re}(z_{k})/\gamma$. The dynamics of the state of $S$ is given by tracing out the degrees of freedom of the extended Hilbert space, $\rho(t)=\operatorname{Tr}_{\backslash S}\rho_{\mathcal{S}}(t)$. This is similar to the pseudomode formulation of non-Markovian dynamics [58, 79, 80, 74], where the only difference compared to the Markovian case is that the dynamics of $S$ is defined on an extended space $\Gamma$ with auxiliary bosons characterised by $\hat{a}_{k}$. In Supplementary Section S1.2, we show how to get the Lindblad dynamics from the auxiliary paths defined on the product space. Compared to the formulation by DQME [63], the dissipative operators and density operators are defined upon the Hilbert space $\mathcal{H}_{S}\otimes\mathcal{H}_{F}$. If the amplitude is not constrained to be real, then the reduced system dynamics is described by a quasi-Lindblad master equation defined on the extended space. In this case, the vectorised form of the generator $\hat{\mathcal{L}}_{S}$ can still be written as ${\mathcal{L}}_{S}=-i(H_{1}-iH_{2}),$ where $H_{1}$ and $H_{2}$ denote, respectively, the Hermitian and anti-Hermitian matrices, see Eq. S15 in Appendix.

II.2 Simulating the open-system dynamics

From the above discussion, simulating open-system dynamics can be mapped onto the simulation of a Lindbladian (or quasi-Lindbladian) evolution. Both Markovian and non-Markovian dynamics take a form similar to Eq. 2 with the latter requiring an enlarged (dilated) space. Such dynamics can be simulated either by directly evolving the dilated system or through stochastic approaches based on pure-state trajectories. Here, we mainly focus on the latter. To that end, we consider the stochastic Schrödinger equation’s approach to simulate the Lindblad equation. To do so, one may employ either a discretised or a continuous method. The circuit complexity in the discretised approach inevitably depends on the number of steps. To avoid this, we present a stochastic Monte Carlo sampling method. The algorithm is described below.

The jump time is computed classically based on the total jump probability (related to $\delta p_{k}$ and $\eta_{k}$) using the sampling method in [68], where the key task is to compute the probability that there will be no jump in a given time period. The total jump probability in the short duration from $t$ to $t+\delta{t}$ is $P_{\textrm{jp}}(t)={\gamma}\sum_{k}(1-\eta_{k})\delta p_{k}(t)$. Given the state undergoes a quantum jump at time $t$, the probability that there will be no jump until $t+\tau$ is

[TABLE]

From this, we can compute the sequence of jump times ${t_{i}}$, which in turn determines the ensemble of trajectories; specifically, one can generate a random number $q\in[0,1]$. Solve $\tau=\mathrm{argmax}_{\tau}(Q(t+\tau)\leq q)$, which is equivalent to finding the largest $\tau$ such that $\int_{t}^{t+\tau}P_{\textrm{jp}}(\tau^{\prime})d{\tau}^{\prime}\leq-\ln(q)$.

In case there is no jump in a given time period, the state continues to evolve under the effective operator

[TABLE]

where the subtraction of $\braket{L_{k}^{\dagger}L_{k}}$ is due to the normalisation after applying $e^{-\frac{1}{2}L_{k}^{\dagger}L_{k}\delta t}$. In the case of Markovian dynamics, $H=H_{S}$ is the system Hamiltonian and $L_{k}$ is introduced in Eq. 2. For non-Markovian dynamics, $H=\hat{H}_{\mathrm{pm}}$ and $L_{k}=\alpha_{k}\hat{a}_{k}$ in Eq. 7. For each trajectory, we start from an initial state $\ket{\psi(0)}$ which we evolve under Eq. 9 by using the QITE algorithm [45] until the jump time. This process is repeated until the final simulation time $T$ is reached, generating a single trajectory. Observable quantities are then estimated by averaging over many such trajectories.

The complexity of the simulation is mainly determined by the cost of simulating the non-Hermitian operator $H_{\text{eff}}$ and quantum jumps on a digital quantum computer. Quantum jumps interrupt the dynamics only through interleaving with the jump operator $L_{k}$, which has a singular value decomposition $L_{k}=U_{k}D_{k}V_{k}$, with unitaries $U_{k}$ and $V_{k}$, and a diagonal operator $D_{k}=\sum_{j}a_{k,j}\ket{j}\bra{j}$. As shown in [68], $D_{k}$ can be realised by a nonunitary exponentiated operator $D_{k}\approx\exp(-H^{D_{k}}T^{D_{k}})$, with $-H^{D_{k}}T^{D_{k}}=\sum_{a_{k,j}\neq 0}\log(a_{k,j})\ket{j}\bra{j}-b_{D_{k}}\sum_{a_{k,j}=0}\ket{j}\bra{j}$. Here, $b_{k}$ is an appropriately chosen constant, determined by $\braket{D_{k}^{2}}$, that bounds the approximation error, see Supplementary Section S2 for details and examples. As $L_{k}$ usually acts locally, $H^{D_{k}}$ is also local, making it efficient to approximate $D_{k}$ by the QITE algorithm.

We remark that the average number of jump events is $\mathcal{O}(\gamma T\sum_{k}\alpha_{k})$, which is linear in $T$. This realisation of the jump integrates naturally with the QITE framework. Therefore, in what follows, we primarily consider the cost associated with simulating $H_{\text{eff}}$ and consider a continuous evolution up to time $T$.

Sketch of the quantum simulation procedure. From the above discussion, it is easy to see that the Markovian 1 and non-Markovian dynamics 7 can be simulated using the stochastic method. The procedure is as follows. Evolve the state under $H_{\rm eff}$ until the next jump happens. The probability that the next jump happens at time $t+\tau$ is $Q(t+\tau)$, based on which the jump time $\tau$ can be determined. Then, we update the quantum state by applying the quantum jump $L_{k}$ which is selected according to the weight $\alpha_{k}$.

For non-Markovian dynamics with pseudomode decomposition, given the bath correlation function, truncate pseudomodes to $K$ with which to construct the effective Hamiltonians $\hat{H}_{\mathrm{eff}}(t)\in\mathcal{H}_{S}\otimes\mathcal{H}_{F}$. If we choose to constrain $d_{k}$ to be real, then we can use the stochastic method with pure states only to simulate the dynamics $\hat{L}_{S}$. In practice, we truncate the energy level of the pseudo-Fock space to $N_{B}$ levels. Then, the total number of qubits is $N+K\log(N_{B})$. Encode the initial state on a quantum computer. In the general case where $\operatorname{Im}(d_{k})\neq 0$, the dynamics is given by Eq. S15. In this situation, one can perform vectorisation to simulate the Hermitian and anti-Hermitian components in $\hat{L}_{S}$, described in Eq. S16. The total number of qubits is $2(N+K\log(N_{B}))$. The rest of the procedure is thus similar to that for simulating the linear Lindblad master equation.

Remarks on the role of dissipation. QITE-based simulations of dissipation typically enlarge the effective support of the evolution operator, suggesting a higher complexity compared to purely unitary dynamics. Still, this increase is not universal. In the case of strong dissipation—where the dissipative term acts effectively like a projector—the complexity is determined mainly by the projector’s support and becomes independent of the total evolution time. Consequently, for a fixed $t$, increasing $\gamma$ can, in fact, reduce the required support, and thereby lower the QITE simulation complexity, Fig. S1. Moreover, even for weaker dissipation, its effects accumulate over time and can eventually dominate the dynamics, again leading to reduced complexity. This supports our conjecture that there exists a timescale beyond which the complexity of simulating open quantum systems is below that of Hamiltonian systems. Yet, at very short times, the situation is more subtle: the use of a unitary approximation to represent the nonunitary evolution can temporarily increase the complexity.

To enable a fair comparison between open quantum system and unitary dynamics in later section, we introduce the operator $\tilde{H}_{\text{eff}}=H_{\text{eff}}/\|H_{\text{eff}}\|$, i.e., a normalised version of Eq. 9. Explicitly, we divide Eq. 9 into the Hermitian and non-Hermitian part, $\tilde{H}_{\text{eff}}=\lambda_{1}H_{1}-i\lambda_{2}H_{2}$, where $\lambda_{1}^{2}+\lambda_{2}^{2}=1$, and $\|H_{1}\|=\|H_{2}\|=1$. By normalising the effective operator, the final time can be regarded as stretched and we denote this rescaled final time by $T$. Overall, the procedure involves decomposing $H_{\text{eff}}$ into $H_{1}$ and $H_{2}$, allowing the representation of the time-evolution operator as a product of exponentials $\ket{\psi(t)}=\exp[-i(\lambda_{1}H_{1}t-i\lambda_{2}H_{2}t)]\ket{\psi(0)}$. Using Trotterisation, the nonunitary dynamics is approximated by an alternating sequence of unitary operations, generated by $H_{1}$, and QITE-based steps that incorporate the effects of $H_{2}$ as $\ket{\psi(t)}=[\exp(-i\lambda_{1}H_{1}\delta{t})\exp(-\lambda_{2}H_{2}\delta{t})]^{T/\delta{t}}\ket{\psi(0)}+\mathcal{O}(T\delta{t})$. Here a first-order Trotter decomposition is applied as an example, while $K$th-order Trotter decomposition is used in deriving the analytical bound in Supplementary Section S2. The first part $H_{1}$ can be implemented on a quantum computer by applying standard simulation techniques for Trotterisation [81, 82]. In contrast, the nonunitary component $H_{2}$ cannot be directly realised on quantum hardware and is instead simulated using QITE proposed in [45].

Before we continue, we ought to compare our approach with others in the literature. There are two key aspects that distinguish the above continuous stochastic method from that of the discretised approach [53]. First, within the QITE framework, the described quantum jump process enables both nonlinear and linear Lindbladians in Eq. 9 to be simulated deterministically, in contrast to the probabilistic realisation of the nonlinear Lindbladian simulation with two ancillary qubits per timestep [53]. Second, this approach does not explicitly involve Trotterisation but rather takes a stochastic, random-sampling approach that recovers the ideal trajectory on average. It is thereby compatible with digital-analogue types of quantum simulation, which may not necessarily incur the additional cost associated with Trotterisation. Moreover, the above method directly applies to time-dependent Lindbladians, as we show in Fig. S2 of the Supplementary Section S4.

Finally, we comment on the possible strategies for simulating the Lindbladian master equation by density matrix vectorisation, where the density matrix is mapped to a pure state in a Liouville space [30]. The dynamics is then simulated by applying the Lindbladian superoperator to this vectorised state [83]. As with the generalised evolution operator, the Lindbladian can also be decomposed into normalised Hermitian and non-Hermitian parts. Here, one can find that the vectorised Lindbladian in [47] takes a similar form to Eq. 9, which can be divided into Hermitian and non-Hermitian parts, but acting on two copies of the state. Our stochastic approach is different from the methods for open quantum systems proposed in [46] and [47], where the vectorised form of the system’s density matrix $\rho$ is encoded in a larger state vector. A key consequence is that they can only preserve $\operatorname{Tr}(\rho^{2})$, rather than $\operatorname{Tr}(\rho)$. As a result, the observable expectation value is calculated by $\braket{O}/\operatorname{Tr}(\rho)$, with $\operatorname{Tr}(\rho)$ being measured at each timestep. Note that the complexity analysis of the full evolution of the vectorised density matrix remains the same for the pure-state evolution, as it adopts a similar form with unitary and nonunitary components, see Supplementary Section S1.1.

For non-Markovian dynamics, the formulation given by Eq. 7 is norm-conserving. In contrast, in DQME-SQ the state dynamics are non-Hermitian and not trace-preserving. The norm of the state may be decreased with both time $t$ and the number of dissipatons in the extended space, due to the non-Hermitian term of $\hat{a}_{k}^{\dagger}$ alone (instead of $\hat{a}_{k}^{\dagger}+\hat{a}_{k}$). Therefore, the resulting states carry a failure probability after post-selection, which will introduce additional sample complexity when implementing on a quantum computer.

II.3 Physical measures of quantum algorithmic complexity

We now wish to analyse the quantum algorithmic complexity in terms of physically interpretable quantities. It is well known that in the short-time regime, open quantum systems exhibit growing spatial correlations due to information propagation, which influences the computational resources required to simulate them [84, 85].111In closed many‑body systems, this growth is constrained by the Lieb–Robinson bounds [86]. For a quantum system $S$ in state $\ket{\Psi_{S}}$, usually the following definition is adopted (see [87, 88]): we say that $\ket{\Psi_{S}}$ exhibits $\xi$-exponential decay of correlations if there exists a correlation length $\xi>0$ such that for all pairs of disjoint regions $(A,B)$ with separation $d(A,B)\geq l_{0}$, where $l_{0}\in\mathbb{R}^{+}$ is a fixed distance threshold, and for all operators $\hat{A}\in\mathcal{A}(\mathcal{H}_{A})$, $\hat{B}\in\mathcal{A}(\mathcal{H}_{B})$, the bound

[TABLE]

holds. Here, $\xi$ corresponds to the slowest decaying correlation in $S$ across all region pairs and all times [89], and the correlation function reads

[TABLE]

where for a pair $(A,B)$ at distance $d$, one can always define a correlation length $C_{A,B}:=-\frac{d(A,B)}{\log\left(\frac{|\mathrm{Cor}(A,B)|}{\|{A}\|\|{B}\|}\right)}\leq\xi$ that satisfies the definition of correlation length for a given choice of regions. This bound applies for a fixed state $\ket{\Psi_{S}}$.

In this work, we wish to establish a tighter bound that can track correlation decay during the time evolution of $S$. Particularly, we consider a pair-specific, distance- and time-dependent correlation length as a replacement of $\xi$. More formally, we take $C_{A,B}(t)$ to be the infimum over all such $\mathrm{Cor}(\hat{A},\hat{B})$

[TABLE]

where $C$ controls how fast the correlations between $A$ and $B$ decay with $d$. From all pair-specific correlation lengths, we can thereby recover $\xi$ as the worst-case across all pairs of regions $\xi(t):=\sup_{A,B}C_{A,B}(t)$.

To define the average of Eq. 11, let $\mathcal{P}_{d}$ denote the set of all pairs of regions $(A,B)$ at distance $d$, then

[TABLE]

From the concavity of the function $f(x)=e^{-d/x}$ when $x\geq{1}/{2}$, and using Jensen’s inequality, we find that

[TABLE]

when $\min_{i}C_{i}\geq 1/2.$ The average correlation length is introduced to capture the general case when $C_{A,B}(t)$ varies across different pairs of regions. The correlation length $\xi(t)$ captures the worst-case behaviour, whereas the average behaviour is captured by $C_{A,B}(t)\leq\xi(t)$, which gives a tighter bound for that pair as it is defined by the smallest value of $C$ that satisfies the inequality 11. These quantities will later be used to relate the error bounds and algorithmic complexity. As we shall see later, the algorithmic complexity is mainly determined by the average correlation length $\overline{C}$. In Supplementary Section S2, we describe its scaling, including Trotter errors and the dependence of the total simulation error on $\overline{C}$.

The above discussion focuses on methods that operate within a fixed Hilbert space. Before we move to the numerical simulation, we would like to comment on how the complexity scales when we allow for systematic expansion via additional ancilla qubits. Due to locality and finite correlation length, each of the $N$ sites in $S$ is influenced primarily by bath modes within a spatial radius $\xi$. Therefore, for each site per each timestep, we only need to keep track of the part of the environment memory within distance $\xi$. Putting all together, the total complexity (i.e., ancilla $\times$ circuit depth) is $Nm\xi N_{\text{steps}}=\mathcal{O}(N\xi(\lambda T)^{1+o(1)}/\varepsilon^{o(1)})$, which is polynomial in correlation length and simulation time, where $N_{\text{steps}}$ is the number of Trotter steps. This may support the conjecture that under reasonable physical assumptions, the cost of simulating an open quantum system with a finite correlation length is efficiently bounded. We refer to Supplementary Section S3 for details. On the other hand, in the asymptotic long-time limit, $S$ approaches a unique steady state $\rho_{\infty}$. Where the rate of convergence is controlled by the spectral gap $\Delta$ of the Liouvillian $\mathcal{L}$, and the mixing time satisfies: $t_{\mathrm{mixing}}\sim\Delta^{-1}.$ After $t_{\mathrm{mixing}}$, the quantum simulation no longer requires evolving time-dependent correlations. The required ancilla size saturates to $C_{\mathrm{mixing}}$, and the number of steps becomes independent of the initial state. Thus, the quantum complexity is saturated up to a certain timescale.

II.4 Investigation on quantum and classical runtime

Quantum runtime analysis. In this section, we study the influence of dissipation on two physically relevant parameters in both the short- and long-time regimes, which will inform us about the quantum algorithmic complexity of simulating the system of interest $S$. In particular, we will focus on the mixing time, which characterises how quickly $S$ reaches its steady state, and the correlation length, which quantifies the spatial extent of quantum correlations within $S$. Both provide physical insight into the dissipative dynamics and serve as key indicators for assessing the simulability of $S$. We numerically analyse the dependence of simulation complexity on the dissipation strength $\gamma$, evolution time $t$, and system size $N$.

We begin by describing our setup. We consider a quantum spin chain system $S$ with $N$ sites, described by a one-dimensional Ising Hamiltonian with transverse and longitudinal fields

[TABLE]

where $J$ is the nearest-neighbour coupling strength, and $h_{x}$ and $h_{z}$ denote the magnetic fields in the $X$ and $Z$ directions, respectively. We set $J=1$ and choose $h_{x}=h_{z}=-2$, as in this regime the expectation values of the local operators typically exhibit oscillatory behaviour, making the model suitable for testing our conjecture about the Lindbladian evolution. Here, $S$ is in contact with an external environment, as described by the time-independent version of Eq. 1, see Fig. 1. We couple $S$ to two baths (modelled as local Markovian reservoirs) at its two edges, $i=0$ and $j=N-1$, and make a simple but common choice

[TABLE]

where $\gamma$ is the system-bath coupling strength which determines the dissipation rate, $\mu$ is the driving parameter, and $\sigma^{\pm}=\sigma^{x}\pm i\sigma^{y}$. In Supplementary Section S5 we present additional results with $h_{x}=0.5$ and $h_{z}=-1.05$, corresponding to a chaotic regime, see Fig. S3. The correlation length there exhibits the same trend, becoming increasingly damped with growing dissipation strength and total evolution time, but more noticeably than in the case studied in the main text.

We fix the driving part of the Lindbladian and are primarily interested in how $\gamma$ impacts the dynamics of $S$. The jump operators $\left\{L_{k}\right\}$ are assumed to act locally at the boundaries of the spin chain, leaving the coherent Hamiltonian-driven evolution in the bulk unaffected. Such boundary-driven setups are conceptually appealing as they are analogous to classical nonequilibrium lattice models. While a rigorous microscopic derivation of these local operators from weak system-bath coupling typically yields nonlocal terms, this approach remains physically insightful. It captures key features of non-equilibrium steady states, quantum transport, and dissipation-induced dynamics. This makes the framework well-motivated for exploring questions in quantum thermodynamics and boundary-driven open quantum systems [90, 91]. In the simplest case of single-site boundary driving, one typically introduces two Lindblad operators that flip a spin up or down with different probabilities, thereby inducing a net magnetisation at the boundary site.

To simulate the open quantum system dynamics, the vectorised density matrix is encoded as a matrix product state (MPS) with physical dimension $4$ for each of the $N$ sites, corresponding to a state in the Liouville space of dimension $4^{N}$, see Supplementary Section S1.1. The Liouvillian superoperator ${\mathcal{L}}$ is represented as a matrix product operator (MPO) with the same enlarged dimension. Time evolution is carried out using the time-evolving block decimation (TEBD) algorithm, where the exponential $e^{{\mathcal{L}}\delta t}$ is approximated via a second-order Trotter-Suzuki decomposition into sequentially applied local exponentials. The timestep $\delta t$ is dynamically adjusted to ensure consistent numerical accuracy throughout the simulation. The classical dependence is related to the bond dimension $\chi$. To ensure the simulation accuracy, we keep track of $\chi$ and check if the error converges with increasing $\chi$.

Before analysing global quantities, it is useful to first examine how dissipation manifests locally in $S$. Local magnetisation observables offer a direct window into the interplay between coherent dynamics and environmental noise: they reveal how quickly local excitations lose their quantum character and relax towards equilibrium. By studying their time evolution for different values of $\gamma$, we can build physical intuition for the mechanisms that ultimately determine the simulability of $S$. Here, for each site $r$ and spin component $a\in\{X,Y,Z\}$, the expectation value is given by $\langle\sigma^{a}_{r}(t)\rangle=\operatorname{Tr}[\sigma^{a}_{r}\rho(t)],$ where $\rho(t)$ is the density matrix evolved under Eq. 1 with Eq. 13 as the Hamiltonian. Fig. 2(a)-(c) show the time evolution of the expectation values of these local Pauli operators for various $\gamma$. As expected, increasing $\gamma$ results in stronger damping of the oscillations, reflecting how stronger environmental coupling suppresses coherent dynamics and drives $S$ towards a steady state at long times. In the closed case, oscillations persist over long times, reflecting coherent spin exchange in the Heisenberg model. As $\gamma$ increases, these oscillations decay more rapidly and their amplitude is strongly suppressed, signalling the progressive loss of phase coherence and the dominance of relaxation processes. Physically, this means that information about the initial state is erased faster, $S$ approaches a featureless steady state sooner, and the dynamical complexity is reduced. In the broader context of simulability, such damping hints at shorter correlation lengths and mixing times, which can make long-time dynamics more accessible to classical methods and reduce the potential advantage of quantum algorithms in this regime.

The damping of local magnetisation oscillations observed in Fig. 2(a)–(c) already suggests that dissipation accelerates the relaxation of $S$ towards equilibrium. To quantify this effect in a more global and model-independent way, we turn to the mixing time—the timescale over which $S$ effectively reaches a steady state. This provides a direct indicator of how long coherent or correlated dynamics persist, and therefore how challenging the long-time regime may be to simulate. We quantify this by computing the trace distance $D(t)=\|\rho(t)-\rho(t-\delta t)\|$ between density matrices at successive timesteps for various values of $\gamma$. Defining a small threshold $\varepsilon$, $S$ is considered to have effectively reached its steady state when $D<\varepsilon$. The results in Fig. 2(d) show that stronger dissipation leads to faster convergence, i.e., shorter mixing time, consistent with the relaxation of local observables in Fig. 2(a)-(c).

We then proceed to the relation between the correlation length and simulation complexity. Spin correlations serve as key indicators of the required computational resources as the presence of long-range correlations across distant sites generally means increased complexity [38]. Crucially, the correlation length bounds the number of measurements. The spin correlation function encodes these spatial dependencies

[TABLE]

where $i$ and $j$ denote different sites. Fig. 3(a) shows the correlation length $C^{z}_{d}(t)$ at distances $d\in{1,\dots,N-1}$ for a fixed long evolution time. The results demonstrate a systematic decrease in $C^{z}_{d}(t)$ as $\gamma$ increases. This behaviour indicates suppression of long-range correlations, and consequently, a reduction in the complexity of simulating strongly dissipative open quantum systems. Besides, to provide a broader perspective, Fig. 3 shows the average spatial correlations, which more clearly highlight the reduction of correlation length—and hence simulation complexity—with increasing dissipation at fixed time. Here we present results at relatively long times, where this behaviour is most pronounced, although a similar (albeit less drastic) trend is already visible at shorter times. Since, according to section II.3, the running time for simulating open quantum system dynamics on a quantum computer is bounded by the maximum correlation length across all distances at each timestep, we also compute this quantity throughout the evolution of $S$ for each $\gamma$, see Fig. 3(b). We note that correlations are not necessarily suppressed by long-time evolution alone: due to the oscillatory behaviour of the Hamiltonian, they can persist or even revive at later times. However, dissipation introduces a damping effect that becomes increasingly dominant with time, leading to a systematic reduction of correlations. This trend is consistent with the average correlation results in Fig. 3(a), where dissipation clearly shortens the correlation length. Together with the analytical scaling in section II.3, which links time complexity to the correlation length, these results indicate that dissipation effectively reduces the simulation complexity—even if correlations themselves do not always decay monotonically with total time evolution.

We also investigate time-dependent Lindbladian dynamics in Supplementary Section S4, where either the dissipation strength, the Pauli operators, or the sites on which the jump operators act vary during evolution. Overall, we observe qualitatively similar behaviour to the results discussed here, with no significant deviations. At the same time, we find that randomised application of jump operators can accelerate relaxation processes, for instance by enhancing the decay of correlations, thereby modifying the system’s dynamical pathways. This further supports the generality of our conclusions and highlights the potential of such protocols for applications such as dissipative quantum state preparation, where time-dependent Lindbladians are often required to engineer the desired steady state [69]. Finally, we wish to remark that while the discussion in this section has focused on linear Lindbladian dynamics, nonlinear master equations could be even more challenging to simulate classically, as their explicit matrix representation is already difficult to construct. Importantly, the quantum algorithm we introduce in this work can efficiently simulate both linear and nonlinear evolutions, making our conclusions directly relevant to the broader setting.

Classical runtime analysis. We now analyse the classical resource requirements—runtime and bond dimension—and provide simulation results to illustrate how simulation errors depend on time $t$, system size $N$, and bond dimension $\chi$. Generally, we expect simulation accuracy to improve for larger $\chi$ at the cost of increased computational runtime. The objective here is twofold: (i) to determine whether the growth trend differs between quantum and classical algorithms, and (ii) to assess whether classical algorithms exhibit rapid growth in complexity with increasing $N$ and $t$.

Figure 4 illustrates the classical simulability of open quantum systems through that accuracy-complexity trade-off. For smaller systems, we quantify simulation errors by comparing results at finite $\chi$ against the ground-truth solution computed using a sufficiently large bond dimension, chosen here as $\chi=250$ to ensure convergence of the expectation values of the observable (see Fig. S4 in Supplementary Section S5). The accuracy of TEBD simulations at given $\chi$ with respect to the ground truth is measured via spin-spin correlation functions, with deviations quantified as the correlation error

[TABLE]

where $\tilde{c}_{ij}$ are the ground-truth correlations. In Fig. 4(a), we show how the correlation error $\varepsilon$ evolves over $t$ for $N=25$ and $\chi\in[50,225]$. As expected, larger $\chi$ consistently yield better accuracy, especially at longer times. Since $S$ eventually reaches a steady state, $\varepsilon$ saturates to a constant value, which is why we think the error remains bounded over long simulation times. In a similar manner, in Fig. 4(b) we fix $t=3$ and examine $\varepsilon$ with respect to $N$ for $N\in[10,25]$. The results confirm the linear growth observed earlier. The main contribution to the complexity of the algorithm, however, arises from its scaling with $\chi$, as it determines the fidelity between our MPS state and the evolved quantum state. In Fig. 4(c), $\varepsilon$ is plotted on a logarithmic scale as a function of $\chi$ for $N\in[10,25]$ and fixed $t=3$. Evidently, for a given $\chi$, smaller systems exhibit smaller $\varepsilon$, making them easier to simulate. For practical purposes, in Fig. 4(d) we set an error threshold $\epsilon_{\chi_{\min}}\leq 0.05$ and identify the minimal bond dimension $\chi_{\text{min}}$ necessary to reach that threshold for every $N$. Naturally, $\chi_{\text{min}}$ depends on $t$, as errors consistently increase with longer simulation times, illustrated clearly in Fig. 4(d). Overall, as expected, we see that larger $\chi$ yield more precise results, progressively approaching exact solutions. This confirms that dynamics can be approximated with arbitrary accuracy given sufficiently large $\chi$. In our case, we find that $\chi=200$ is sufficient to reach the error threshold for the system sizes considered, so working with samples up to $\chi=250$ is well justified. This is consistent with the results in Supplementary Section S5, where we observe clear convergence of the observables with increasing bond dimension, further supporting that $\chi=250$ provides a reliable estimate.

The classical simulability of TNs is known to be fundamentally constrained by the entanglement entropy of the quantum state. To assess this limitation, we analyse the growth of entanglement entropy as an indicator of computational complexity. In Fig. 5, we study the operator space entanglement entropy of the time-evolved state. We find that although the entanglement entropy initially increases, it then exhibits qualitatively different behaviour depending on $\gamma$. For $\gamma$ above a certain threshold, the entanglement entropy decreases over time, consistent with our expectations and the results presented by [39]. However, when considering simulation complexity, a more relevant metric is the bond dimension $\chi$ required to represent the state. Surprisingly, we find that $\chi$ remains largely constant across different values of $\gamma$. This result is somewhat counter-intuitive compared to the behaviour of the entanglement entropy. The key difference lies in the fact that the entanglement entropy characterizes the final state, whereas simulating the full dynamics requires tracking the entire evolution, which can be more demanding. As a result, $\chi$ does not decrease with $t$ even for different values of $\gamma$.

To give a brief definition of the operator entanglement, we define our density matrix in MPO form as a Schmidt decomposition between the left $L$ and right $R$ block

[TABLE]

where $\operatorname{Tr}(\hat{o}^{[L/R]}_{i}\hat{o}^{[L/R]}_{j})=\delta_{ij}$ and $\lambda_{i}$ denotes the Schmidt coefficients. Then, define the operator entanglement as $S_{\text{OP}}=-\sum_{i}\lambda_{i}^{2}\log(\lambda_{i}^{2}).$ $S_{\text{OP}}$ is not strictly connected with quantum entanglement between blocks but it is an important quantity nonetheless that gives insights into the efficiency of the MPO representation. Recall that, for the purpose of our simulation, we have encoded the density matrix MPO into a MPS with double the quantum dimension. In this case, the operator entropy of $\rho$ is equivalent to calculating the entanglement entropy of the evolved MPS, carefully handling the extra physical dimension.

A noteworthy observations is that the complexity of TNs does not decrease with $t$ and $\gamma$ in certain regimes, in stark contrast to the complexity associated with quantum algorithms. This suggests there may exist a range of dissipation strength, where the separation between quantum and classical computational complexity is maximised.

III Discussion

In this work, we examined the algorithmic complexity of simulating open quantum systems under the conjecture that, beyond a threshold timescale, Lindbladian dynamics may be easier to simulate than Hamiltonian dynamics. Instead of relying on dilation methods that obscure reduced-system complexity [15, 16, 18], we focused on correlation length and mixing time, employing a quantum jump approach—via QITE for the quantum simulations and TNs for the classical simulations. We introduced two post-selection-free algorithms that simulate both linear and nonlinear Lindbladians without ancilla overhead and normalisation issue [46, 47] of vectorisation by averaging the trajectories. The stochastic, continuous simulation method enables stochastic simulations of open dynamics, which may serve as a promising protocol for digital-analogue quantum simulations. Our framework extends naturally to non-Markovian dynamics defined on an extended Hilbert space.

A key feature of open quantum systems is the competition between dissipative and unitary dynamics, where dissipation can either increase or suppress correlations or entanglement, thereby directly influencing the system’s complexity [92]. The empirical results found so far suggest this behaviour is not universal, but in fact can highly depend on the type of state, the underlying dynamics, and the specific form of dissipation, as numerically found in [38, 39]. To quantify this, we consider a normalised Lindbladian and track the correlation decay during the evolution. We analyse the algorithmic complexity in relation to the averaged region- and time-dependent correlation length in Eq. 12, which is dynamic, and strictly smaller than the standard correlation length (e.g., [87, 88]). From a computational perspective, our analytical result does not preclude the possibility that an intermediate time scale may exist, where dissipation increases the complexity beyond that of its closed counterpart. In particular, in the main part of this work we showed numerically that stronger dissipation shortens mixing time and reduces correlation length. This result is similar to findings in the SYK model, where Krylov complexity was used to quantify the complexity growth, and dissipation was shown to suppress it compared with the closed counterpart [93].

An noteworthy observation is that we observed a contrast between classical and quantum resource scaling: while TN simulations exhibit super-polynomial growth, quantum resources scale much favourably. Interestingly, dissipation suppresses correlations but does not reduce the TN bond dimension, showing there is a separation between quantum and classical computational complexity. This finding reveals the competition between dynamical correlations and transient entanglement and the distinction between them.

Our work provides insight into how dissipation shapes both classical and quantum algorithmic complexity via physically meaningful measures. These results may inform us in designing more efficient algorithms for simulating open quantum systems, and identifying regimes where quantum algorithms can outperform classical ones. Some interesting directions for future investigation emerge. Specifically, in preparing many-body ground states, which is known to be challenging even for quantum computers [94], a key question for dissipative quantum algorithms is how to reduce the complexity of open dynamics and the mixing time [69, 70]. In Supplementary Section S4, we highlighted the role of randomness in the steady-state properties of local observables. Specifically, we demonstrated that applying jump operators randomly, rather than at fixed boundary sites, accelerates dissipative state preparation, suggesting that randomised dissipation protocols can be more efficient. This result indicates that one may exploit such dynamics to algorithmically design randomised dissipation protocols that systematically reduce mixing times, and thus enable efficient preparation of many-body ground states. On a related note, it is also worthwhile to determine whether the observed polynomial scaling of quantum algorithms persists in realistic settings like those with nonlocal jump operators or correlated finite-temperature environments, which is an interesting research direction for ground-state preparation.

Acknowledgements.

We would like to thank Yuguo Liu, Yukai Guo and Weitang Li for the useful discussion related to open quantum dynamics simulation. This research receives support from Schmidt Sciences, LLC. L.V.E. was supported by the Clarendon Fund. We would also like to acknowledge funding from the UK EPSRC through EP/Z53318X/1.

Author contribution

J.S. conceived the idea. L.V.E., A.Y., and J.S. developed the theoretical aspects of the project. L.V.E. developed the source code and performed numerical simulations. J.S. carried out the theoretical analysis to support the study. L.V.E., A.Y. and J.S. wrote the manuscript.

Data availability

All data needed to evaluate the conclusions in the paper are present in the paper and the Supplementary Information. The source data is available by request from the authors.

Code availability

The code is available by request from the authors.

Competing interests

The authors have no competing interests.

Supplementary Information

Contents

1. I Introduction

2. II Results

3. II.1 Formulation of open quantum dynamics

4. II.2 Simulating the open-system dynamics

5. II.3 Physical measures of quantum algorithmic complexity

6. II.4 Investigation on quantum and classical runtime

7. III Discussion

8. S1 Formulation of open quantum dynamics

9. S1.1 Vectorisation of the reduced density matrix under linear Lindblad evolution

10. S1.2 Pseudomode formulation of non-Markovian dynamics and beyond

11. S2 Error bound of the evolution and the algorithmic complexity

12. S3 Space-time complexity bounds for simulating Lindbladian dynamics

13. S4 Time-dependent Lindbladian

14. S5 Tensor network simulation results for chaotic systems

Appendix S1 Formulation of open quantum dynamics

In Supplementary Information, the hat notation is used to denote the terms in the extended space and the operators as well.

S1.1 Vectorisation of the reduced density matrix under linear Lindblad evolution

To simulate the quantum dynamics governed by a linear Lindbladian, we encode the density matrix as a pure state through vectorisation. This transformation allows us to represent mixed-state evolution in a form suitable for quantum circuits. The same technique is also applied for classical TN simulations of open quantum systems, as described in the main text The mapping sends $\rho$ to a vector in the Liouville space

[TABLE]

Under this transformation, the Lindblad master equation takes a Schrödinger-like form

[TABLE]

where the Liouvillian superoperator $\hat{\mathcal{L}}$ acts linearly as

[TABLE]

To compare open-system ($\gamma>0$) with unitary ($\gamma=0$) dynamics on equal footing, we consider a normalised Liouvillian

[TABLE]

where $||H_{1}||=||H_{2}||=1$ and $\lambda_{1}^{2}+\lambda_{2}^{2}=1$, with $H_{1}$ and $H_{2}$ is the matrix representation of the superoperator. In analogy with the operator $H_{\text{eff}}$ in Eq. 9, the time evolution of the vectorised density matrix can be written as

[TABLE]

The derivation then proceeds by applying a first-order Trotter decomposition, as an example

[TABLE]

which adopts the same form as in the quantum jump situation. In Supplementary Section S2, $K$th-order Trotter decomposition is applied. The dynamics of the vectorised state is then approximated as alternating applications of unitary operations and operations derived from the QITE approach.

S1.2 Pseudomode formulation of non-Markovian dynamics and beyond

In open quantum dynamics, when the environment has a non-negligible memory, we need to go beyond the Lindblad master equation to describe the reduced system dynamics. A numerically exact way to deal with the non-Markovian dynamics is to use the HEOM [57], which contains a hierarchy that connects to different auxiliary density operators. Similarly, Yan et al. [77, 65] proposed the dissipaton equation of motion, where quasi-bosonic particles (i.e., bosonic dissipatons) are introduced as auxiliary density operators to encode the non-negligible system-environment interactions. Later, this was extended to the DQME formalism, where instead of a hierarchical structure, the reduced system dynamics is incorporated into a single dynamic equation, and then extended to a second-quantised form of DQME, where the reduced density tensor of the dissipation embedded system $\tilde{\rho}\in\mathcal{H}_{S}\otimes\mathcal{H}_{S}^{*}\otimes\mathcal{H}_{F}$ (with the Fock space $\mathcal{H}_{F}$) can be implemented by using linear combination of unitaries algorithms.

As we discussed in the main text, one way of unravelling the influence functional of the system-bath interaction is to decompose the bath correlations into a finite set of modes, cf. Eq. 5. This also aligns with the idea that we can rewrite the memory kernel of the reduced dynamics into a set of time-local interactions between the system and a finite number of pseudo-modes, and use the pseudo-mode Lindbladian [55] as an approximation to the non-Markovian dynamics when $d_{k}\in\mathbb{R}^{1}$.

Along those lines, a few approaches have been proposed for simulating non-Markovian dynamics based on this extended space, see [63]. From now on, we mainly follow the method in [67] which unravels influence functions using complex-valued auxiliary paths. We will thus provide a self-contained procedure for representing non-Markovian dynamics using the Lindblad-type equation in an extended space.

The reduced density operator ${\rho}_{S}$ is expressed as a functional integral over paths supported by a Keldysh contour

[TABLE]

and the bare action factor $e^{iS\left[q_{S}^{+},q_{S}^{-}\right]}$ captures the quantum dynamics in the absence of a bath. Here, $S[q^{+}_{S}]$ and $S[q^{-}_{S}]$ are the corresponding actions associated with, respectively, forward and backward paths $q^{\pm}_{S}(t)$. Formally

[TABLE]

In Sec. II in the main text, we presented a common approach in this regime, where we represented the environment through its bath correlation function Eq. 5, and expanded it as a finite sum of complex exponentials with real amplitudes $d_{k}$ and complex frequencies $z_{k}$. This decomposition enabled the mapping of the system’s dynamics onto an extended space $\Gamma$. Building on that, we have

[TABLE]

As derived in [67] the actions are

[TABLE]

and $H_{k}\left(\phi_{k},\phi_{k}^{*}\right)=-iz_{k}\phi_{k}^{*}\phi_{k}-i\sqrt{d_{k}}\left(q_{S}^{+}-q_{S}^{-}\right)\phi_{k}+i\sqrt{d_{k}}q_{S}^{+}\phi_{k}^{*},$ and the intermediate interaction terms are introduced $\bar{H}_{k}\left(\psi_{k},\psi_{k}^{*}\right)=-iz_{k}^{*}\psi_{k}^{*}\psi_{k}-i\sqrt{d_{k}^{*}}\left(q_{S}^{-}-q_{S}^{+}\right)\psi_{k}+i\sqrt{d_{k}^{*}}q_{S}^{-}\psi_{k}^{*}.$

The idea behind dissipation-embedded quantum master equation [63] or QD-MESS [67] is to introduce an extended space $\Gamma$ and derive the dynamics of states in the extended space. Below, we provide a few details in arriving at the evolution of states given by Eq. 7.

For a bosonic mode with annihilation and creation operators $\hat{a},\,\hat{a}^{\dagger}$ and a normal-ordered Hamiltonian $\widehat{H}:=H_{N}(\hat{a}^{\dagger},\hat{a})$, the coherent-state matrix element of the short-time propagator is

[TABLE]

Since $\hat{a}|\phi_{n}\rangle=\phi_{n}|\phi_{n}\rangle$ and $\langle\phi_{n+1}|\hat{a}^{\dagger}=\langle\phi_{n+1}|\phi_{n+1}^{*}$, expanding to first order in $\Delta t$ gives

[TABLE]

and we have the standard short-time kernel identity

[TABLE]

up to second order. Now let’s define

[TABLE]

which is the propagator acting on the extended density-like object $W(0)$, but compared to the Schrödinger equation, it includes auxiliary-mode degrees of freedom. Differentiating with respect to time yields

[TABLE]

which is the operator form of the evolution of the extended object $W(t)$ in Liouville-Fock space. Consider now one auxiliary mode with complex frequency $z_{k}$. The quadratic (free) part of the coherent-state action for this mode is

[TABLE]

with the normal symbol $H_{N,k}^{\mathrm{free}}=-iz_{k}\phi_{k}^{*}\phi_{k}.$ Hence, the Lagrangian corresponding to this free part is

[TABLE]

where the term $i\phi_{k}^{*}\dot{\phi}_{k}$ arises from the coherent-state overlap between neighbouring time slices (see Eq. S10), while the term $iz_{k}\phi_{k}^{*}\phi_{k}$ comes from $H_{N,k}^{\mathrm{free}}=-iz_{k}\phi_{k}^{*}\phi_{k}$. Replacing the c-number fields by their corresponding ladder operators, $\phi_{k}\mapsto\hat{a}_{k},\ \phi_{k}^{*}\mapsto\hat{a}_{k}^{\dagger}$, gives the normal-ordered Hamiltonian operator

[TABLE]

According to Eq. S11, the contribution of this mode to the time derivative of $W$ is

[TABLE]

Similarly, for the backward branch, one obtains $-z_{k}^{*}\,\hat{b}_{k}^{\dagger}\hat{b}_{k}W$. Therefore, we have the generator $\hat{\Gamma}_{k}=-z_{k}\hat{a}_{k}^{\dagger}\hat{a}_{k}-z_{k}^{*}\hat{b}_{k}^{\dagger}\hat{b}_{k}$ that is introduced in the main text. Now the forward coherent fields correspond to operators acting from the left on the extended object. Specifically, the mapping to operators acting on $W$ is $\phi_{k}\mapsto\hat{a}_{k},\phi_{k}^{*}\mapsto\hat{a}_{k}^{\dagger}\quad(\text{left action}),$ $\psi_{k}\mapsto\hat{b}_{k},\psi_{k}^{*}\mapsto\hat{b}_{k}^{\dagger}\quad(\text{right action}).$ Forward system coordinate $q_{S}^{+}$ maps to left multiplication by $\hat{q}_{S}$, while backward coordinate $q_{S}^{-}$ maps to right multiplication by $\hat{q}_{S}$, which is

[TABLE]

It is easy to check that $\left(q_{S}^{+}-q_{S}^{-}\right)\phi_{k}=[\hat{q}_{S},\hat{a}_{k}{W}]$. The dynamics of an extend state ${W}$ then reads

[TABLE]

When considering the decomposition

[TABLE]

we will arrive at the HEOM of ${\rho}_{\mathbf{m},\mathbf{n}}(t)$, which could be shown to be equivalent to that in [66]. One may then apply the similarity transformation $\mathcal{S}$ discussed in the main text that mixes $\hat{a}_{k}$ and $\hat{b}_{k}$, allowing the dynamics to be expressed in Lindblad form. Besides, recall that to enable implementation on a quantum computer, the state was further defined using only half of the auxiliary bosonic modes in a dual Fock space ${\rho}=\sum_{\mathbf{m},\mathbf{n}}\ket{\mathbf{m}}\rho_{\mathbf{m},\mathbf{n}}\bra{\mathbf{n}}$. In doing so, the dynamics is mapped to a quasi-Lindblad

[TABLE]

One can find that the density matrix is not positive but is Hermitian and trace-preserving $\operatorname{Tr}({\rho}_{\mathcal{S}})=1$. It is easy to check that the reduced density is also trace-preserved $\operatorname{Tr}({\rho}_{S})=1$. When $d_{k}$ is constrained to be real, then we recover the Lindbladian 7 in the main text. It is moreover interesting to compare this with the one obtained from pseudomode theory in [55, 74, 56], where $d_{k}$ is real. The dynamics of the state of $S$ is given by

[TABLE]

where it directly encodes the real amplitude $d_{k}$ and the complex frequency $z_{k}$ into the Lindbladian. In the more general case with $\operatorname{Im}(d_{k})\neq 0$, we can still simulate the dynamics with a vectorised form. Specifically, c.f. Eq. (37) in [67]

[TABLE]

where the individual superoperators are given as

[TABLE]

and

[TABLE]

In this form, we can still express $\hat{\mathcal{L}}_{S}$ as a normalised Liouvillian

[TABLE]

where $H_{1}=\mathbb{I}\otimes\hat{H}_{\mathrm{eff}}-\hat{H}_{\mathrm{eff}}^{\mathsf{T}}\otimes\mathbb{I}$ and $H_{2}=\sum_{k}\mathcal{L}_{C,k}+\mathcal{L}_{D,k}$ is the matrix representation of the superoperators. Note that this is different from the DQME-SQ in [63], where the generator is not Hermitian but trace preserving. It may be worth noting that one may adopt other unravelling methods to get the effective Hamiltonian, e.g., [76].

Appendix S2 Error bound of the evolution and the algorithmic complexity

Below we analyse the simulation complexity of open quantum system dynamics. The simulation is composed of the dynamics governed by $H_{\rm eff}$ and the quanutm jump. The dynamics under $H_{\rm eff}$ is composed of both real-time evolution and imaginary-time evolution. As we shall see, the realisation of the quantum jump also uses imaginary-time evolution.

We first discuss the complexity in simulating the dynamics under $H_{\rm eff}$. We divide the total time evolution $T$ into $N_{2}$ segments of duration $\delta{t}=T/N_{2}$. For the non-Hermitian component $H_{2}$, non-commuting terms introduce a Trotter error in QITE. To control Trotter error, we further divide each $\delta{t}$ into $N_{1}$ sub-steps. Here, we consider a $K$th-order Trottersation at all levels (with even $K$). The error sources include Trotter error in dividing the evolution into Hermitian and non-Hermitian parts, Trotter error in real-time evolution, and for unitary decomposition of the imaginary-time evolution operator, each with its own error set by $\epsilon/3$. The number of sub-steps $N_{1}$ is set by the condition

[TABLE]

which gives

[TABLE]

The simulation error is governed by the average correlation length $\overline{C}$. For the effective non-Hermitian Hamiltonian $H_{\rm eff}$, we write

[TABLE]

Following the argument in [45], let $\Psi_{p}$ denote the pure state at the $p$th step

[TABLE]

with density matrix $\rho_{p}=\ket{\Psi_{p}}\bra{\Psi_{p}}$. Consider a region $R_{v}$ consisting of all sites within distance $v$ of the support of $h_{k}$. The distance between successive reduced states satisfies

[TABLE]

where $\rho_{p,i}:=c_{p,i}^{-1}\prod_{k=i}^{M_{k}}e^{-2\delta{t}h_{2,k}}\rho_{p-1}$, and $c_{p,i}$ denotes the normalisation factor. In the last line, we have used the result, Lemma 9 of [87], which is used as a key lemma for the result in [45]. The last inequality relates the region-dependent $C_{i}$ to the average correlation length. To be specific, we define $C_{i}$ as the minimal length such that Eq. 11 holds given all the possible pairs $(A_{i},B_{i})$, where $A_{i}$ is the region on which $h_{2,k}$ acts and $B_{i}$ is the region that has distance $v$ to $A_{i}$. Therefore, the relevant parameter here is the average correlation length $\overline{C}$, rather than the maximal one. Choosing the domain size to be $v=\mathcal{O}\left(\bar{C}\ln\left({M_{k}}N_{2}\varepsilon^{-1}\right)\right),$ we can bound the error from approximating the Lindbladian by a sequence of unitaries. Note the reason why we pair the terms in the second line is that we hope to have a large domain $v$ such that each individual term can be bounded by $e^{-v/C_{i}}$.

The sample complexity per timestep for the non-Hermitian evolution is

[TABLE]

At a higher level, $N_{2}$ is fixed by the Trotter error between $H_{1}$ and $H_{2}$

[TABLE]

so that

[TABLE]

Thus, the sampling cost per time step due to the non-Hermitian part is

[TABLE]

In many physically relevant cases, environmental interactions are described by local Lindbladians with commuting terms in $H_{2}$. Then, no $N_{1}$ subdivision is required. With some derivation, the total complexity is given by

[TABLE]

Finally, note that $\overline{C}\leq\xi$, where $\xi$ is the universal correlation length of the system. Physically, $\overline{C}$ measures correlations at a fixed distance $d$ with all possible operators, whereas $\xi$ captures correlations at all distances. Our formulation sharpens the dependence on correlation length, since $\overline{C}$ can be much smaller than $\max_{i}C_{i}$. We note that the above expression characterises the total sample complexity accumulated over all steps. It also reflects the computational complexity associated with the stepwise growth of the circuit (with each block of duration $\delta t$).

The realisation of the quantum jump $L_{k}=U_{k}D_{k}V_{k}$, in particular $D_{k}\approx\exp(-H^{D_{k}}T^{D_{k}})$, is similar. To approximate the diagonal operator $D_{k}$, we can choose to set the exponent $-H^{D_{k}}T^{D_{k}}=\sum_{a_{k,j}\neq 0}\log(a_{k,j})\ket{j}\bra{j}-b_{D_{k}}\sum_{a_{k,j}=0}\ket{j}\bra{j}$. Here, $b_{k}$ is an appropriately chosen constant. In case $\braket{D_{k}^{2}}$ is exponentially small, then one can still choose $b_{k}$ and output a zero vector (which is recorded for classical computation). Otherwise, when $\braket{D_{k}^{2}}=\braket{\Psi_{S}|D_{k}^{2}|\Psi_{S}}$ is lower bounded by $b_{k}$, then one can choose to set that constant as

[TABLE]

such that $\|D_{k}-\exp(-H^{D_{k}}T^{D_{k}})\|\leq\varepsilon$. Here, $\braket{D_{k}^{2}}$ can be directly measured at each timestep, and $H_{D_{k}}$ is local. For example, for $L_{k}=\sigma^{+}_{k}=\ket{1}\bra{0}$, it was shown in [68] that the effective Hamiltonian can be set as $H_{D_{k}}=(Z_{k}+\mathbb{I}_{k})/2$, while $U_{k}=\mathbb{I}$ and $V_{k}=Z$. With this, one can use the QITE algorithm to find a unitary that approximates $D_{k}$.

Appendix S3 Space-time complexity bounds for simulating Lindbladian dynamics

To simulate an open quantum system with correlation length $\xi$ over time $T$, we may simply add more ancilla to mediate the relevant memory effects. The space-time complexity—defined as the total number of ancilla qubits multiplied by circuit depth—is found to be polynomial in $\xi$ and $T$. Here, we provide more details on this.

Assume that the Lindbladian dynamics arises from a unitary evolution of a joint system-bath state followed by tracing out the bath. That is

[TABLE]

with system-bath interactions of the form $H_{SB}=\sum_{j}h_{j,B_{j}}$, where $h_{j,B_{j}}$ acts locally on site $j$ and its associated bath region $B_{j}$. Assume the bath is gapped or thermal, and that system-bath correlations decay exponentially in both space and time

[TABLE]

where $\xi$ is the correlation length and $\tau_{c}$ is the bath memory (correlation) time. Discretise the evolution into $T$ timesteps of duration $\delta t$, so that $t=T\delta t$. At each timestep, the system evolves under a quantum channel $\rho{(t+\delta t)}=\mathcal{E}_{\delta t}(\rho{(t)}).$ Due to the finite memory time $\tau_{c}$, the bath correlations are negligible beyond a temporal window of $\tau_{c}$, implying that each $\mathcal{E}_{\delta t}$ can be modelled with ancillas that only retain correlations within the last $\tau_{c}/\delta t$ steps. From the Stinespring dilation theorem, such a channel can be realised by a unitary interaction between the system and an ancilla of finite dimension, followed by a partial trace

[TABLE]

Finite bath memory time $\tau_{c}$ means the environment (bath) only retains information about the $S$’s past for up to time $\tau_{c}$. Any influence of $S$’s state before $t-\tau_{c}$ on its current evolution is exponentially suppressed. In other words, the bath acts approximately Markovian on timescales larger than $\tau_{c}$.

If the dynamics of $S$ were strictly Markovian, we would only need a single ancilla per timestep (and could even reuse it). However, in general cases, where $S$ exhibits non-Markovianity with a memory time $\tau_{c}$, we must retain information for up to $m=\tau_{c}/\delta t$ previous steps. If the bath memory lasts for a finite time $\tau_{c}=\mathcal{O}(1)$, the number of previous timesteps that affect the current step is bounded by

[TABLE]

where $\lambda$ denotes the interaction strength. This is the memory length for each timestep. The total ancilla usage across all $N_{\text{steps}}$ timesteps is

[TABLE]

Note that in simulating an entire system of $N$ spatial sites, each site evolves under open-system dynamics, meaning it may couple to its own environment (e.g., local bath modes). In the worst case, each site needs its own set of ancilla qubits (or ancilla modes) to mediate interaction with the environment. So the ancilla resource scales linearly with $N$ in the worst case. Due to locality and finite correlation length, each system site is influenced primarily by bath modes within a spatial region of radius $\xi$. Therefore, for each site at each timestep, we only need to keep track of a region of environment memory of size. Putting all this together, the total complexity (i.e., ancilla $\times$ circuit depth) is

[TABLE]

which is polynomial in correlation length and simulation duration. This validates the claim that, under assumptions of system-bath correlation decay, the cost of simulating open quantum systems with finite correlation length is efficiently bounded.

Appendix S4 Time-dependent Lindbladian

Here, we provide details and illustrative results for two alternative dissipative protocols, where the Lindbladian dynamics is time-dependent in some way, focusing on the effects on both local magnetisation and correlations.

In the first protocol, the system evolves under a time-dependent Lindblad equation, where the dissipation strength $\gamma(t)$ varies smoothly during the evolution according to

[TABLE]

Here, $\gamma_{\mathrm{max}}$ sets the maximum dissipation rate reached during the dynamics. The Lindblad operators are fixed as $\sigma^{\pm}$ acting at the two ends of the chain, as in the main text.

In the second protocol, we fix $\gamma$ to a constant value but randomise both the operator type (choosing randomly from $X$, $Y$, or $Z$ Pauli operators) and the sites at which the Lindblad operators act. For each simulation run, two sites are selected at random, and each is assigned a randomly chosen Pauli operator as the Lindblad operator. Figure S2 displays the time evolution of the local magnetisation $\langle\sigma^{z}_{i}(t)\rangle$ at the center of the chain (top row) and the connected correlation function $C^{z}_{0,N/2}(t)$ (bottom row), for both protocols and for two values $\gamma_{\max}=0.1$ and $\gamma_{\max}=0.5$. For the correlations, we focus on the distance $d$ between site $i=0$ and the center of the chain, but note that the observed behaviour is representative for other values of $d$ as well.

These results highlight how both the temporal profile and spatial structure of the dissipation channel can influence relaxation and steady-state properties of local observables in open quantum systems. Notably, we find that the time-dependence of $\gamma$ in the first protocol does not lead to significant differences in either local magnetisation or correlations when compared to the corresponding time-independent case with $\gamma=\gamma_{\mathrm{max}}$. This is expected, as the overall effect of $\gamma$ can be understood as accumulating over time, but never exceeding $\gamma_{\mathrm{max}}$, so that both the relaxation and steady-state behaviour are largely determined by $\gamma_{\mathrm{max}}$. The correlations show a similar insensitivity to the time dependence of the dissipation, being primarily controlled by $\gamma$ itself. This behaviour persists in different choices of $d$.

In contrast, the second protocol, where Lindblad operators act at random sites with randomly chosen Pauli types, produces more noticeable and consistent differences in the observable dynamics. Here, the randomness in the spatial and operator structure of the dissipation can substantially alter the system’s relaxation pathways, leading to a richer variety of dynamical behaviours in the observables, depending on the particular realisation of the dissipative channel. As discussed in the main text, in the context of ground-state preparation through dissipative protocols, this suggests that a time-dependent, randomised application of jump operators could accelerate convergence and enable more efficient state preparation.

Appendix S5 Tensor network simulation results for chaotic systems

Here, we present additional results from TN simulations. Although the main text focuses on the Heisenberg model with oscillatory behaviour in its Pauli observables, in Fig. S3 we show the case for $h_{x}=0.5$, $h_{z}=-1.05$, corresponding to a chaotic regime. Here, the coherent oscillations present at weak dissipation are rapidly suppressed as $\gamma$ increases, leading to faster relaxation towards steady state. This behaviour illustrates how dissipation tends to erase the signatures of chaos in local observables, effectively driving the system into featureless dynamics at long times. Evidently, the correlation length exhibits the same trend, becoming increasingly damped with $\gamma$—more noticeably than in the case studied in the main text—and with total evolution time, as shown in Fig. S1. This confirms that dissipation not only suppresses temporal oscillations but also washes out spatial correlations across the chain, reinforcing the picture that strong coupling to the environment reduces the effective complexity of the dynamics.

On a different note, in this work we set the bond dimension to $\chi=250$ as the ground truth. In Fig. S4 we show that this choice is well justified, with the observable dynamics clearly converging as $\chi$ increases. The convergence is particularly evident at intermediate times, where entanglement growth is more pronounced, yet the difference between $\chi=200$ and $\chi=250$ remains negligible. This validates our numerical accuracy and ensures that all qualitative conclusions drawn in the main text are robust with respect to truncation errors.