# Quantum trajectory entanglement in various unravelings of Markovian dynamics

Tatiana Vovk<EMAIL_ADDRESS>Institute for Theoretical Physics,
University of Innsbruck, 6020 Austria  Institute for Quantum Optics and
Quantum Information, Austrian Academy of Sciences, 6020 Austria  Hannes
Pichler  Institute for Theoretical Physics, University of Innsbruck, 6020
Austria  Institute for Quantum Optics and Quantum Information, Austrian
Academy of Sciences, 6020 Austria

###### Abstract

The cost of classical simulations of quantum many-body dynamics is often
determined by the amount of entanglement in the system. In this paper, we
study entanglement in stochastic quantum trajectory approaches that solve
master equations describing open quantum system dynamics. First, we introduce
and compare adaptive trajectory unravelings of master equations. Specifically,
building on Ref. [Phys. Rev. Lett. 128, 243601 (2022)], we study several
greedy algorithms that generate trajectories with a low average entanglement
entropy. Second, we consider various conventional unravelings of a one-
dimensional open random Brownian circuit and locate the transition points from
area- to volume-law-entangled trajectories. Third, we compare various
trajectory unravelings using matrix product states with a direct integration
of the master equation using matrix product operators. We provide concrete
examples of dynamics, for which the simulation cost of stochastic trajectories
is exponentially smaller than the one of matrix product operators.

##  I Introduction

Entanglement plays a central role in modern quantum science as a resource for
several quantum information processing tasks Wootters (1998a); Amico  _et al._
(2008); Horodecki  _et al._ (2009); Calabrese  _et al._ (2009); Eisert  _et
al._ (2010). At the same time it poses a substantial obstacle to the efficient
simulations of quantum many-body systems on classical computers Vidal (2003);
Eisert and Osborne (2006); Prosen and Znidaric (2007); Schuch  _et al._
(2008). For instance, simulating quench dynamics of closed many-body systems
can lead to the accumulation of extensive entanglement, resulting in
prohibitively large cost of classical simulations Calabrese and Cardy (2007);
Prosen and Znidaric (2007); Schuch  _et al._ (2008); Polkovnikov  _et al._
(2011); Eisert  _et al._ (2015); Nahum  _et al._ (2017). The hardness of
classically simulating such dynamics has been rigorously established via
arguments grounded in complexity theory Aaronson (2005); Aaronson and Chen
(2016); Bouland  _et al._ (2019); Napp  _et al._ (2022) and has spurred recent
interest in developing and implementing protocols aimed at achieving quantum
advantage Bremner  _et al._ (2016); Harrow and Montanaro (2017); Boixo  _et
al._ (2018); Arute  _et al._ (2019); Daley  _et al._ (2022).

Experimental implementations of quantum many-body dynamics on contemporary
quantum hardware are however inherently affected by noise, and thus they are
best described in terms of open quantum systems De las Cuevas  _et al._
(2013); Yung and Gao (2017); Preskill (2018); Cross  _et al._ (2019); Weimer
_et al._ (2021). When the quantum system is open, the complexity of its
classical simulation requires reevaluation, since simulating noisy quantum
dynamics on a classical computer can potentially be more efficient than its
noiseless counterpart. Several classical algorithms for open system dynamics
leverage this fact, albeit mostly indirectly. Among them are algorithms based
on matrix product or neural network density operators that represent and
evolve the full system density operator Werner  _et al._ (2016); White  _et
al._ (2018); Yoshioka and Hamazaki (2019); Vicentini  _et al._ (2019); Carisch
and Zilberberg (2023). Another category encompasses stochastic methods, where
the presence of noise is incorporated explicitly as a stochastic element in
the algorithm. In particular, in quantum trajectory techniques the system
density operator and its dynamics are obtained by statistical averaging over
pure-state wave functions Zoller  _et al._ (1987); Dum  _et al._ (1992);
Dalibard  _et al._ (1992); Gardiner  _et al._ (1992); Gisin and Percival
(1992); Tian and Carmichael (1992); Castin  _et al._ (1993); Van Dorsselaer
and Nienhuis (1998); Plenio and Knight (1998). Tensor network methods such as
matrix product states (MPSs) can sometimes provide an efficient means to
represent each of these many-body pure-state trajectories Schollwöck (2011);
Orús (2014); Daley (2014); Orús (2019).

In Ref. Vovk and Pichler (2022) we have introduced a novel way to directly
leverage noise in trajectory-based stochastic methods. The central premise
relies on the fact that the same system dynamics can be obtained by different
unravelings, which give different ensembles of pure-state trajectories.
Specifically, Ref. Vovk and Pichler (2022) puts forward the idea to adaptively
optimize the unraveling choice to minimize the average entanglement, which
acts as a proxy of the cost of classically representing trajectories. The
physical mechanism underlying this idea is reminiscent of the phenomenon of
measurement-induced phase transitions Li  _et al._ (2018); Skinner  _et al._
(2019); Chan  _et al._ (2019); Li  _et al._ (2019); Choi  _et al._ (2020);
Zabalo  _et al._ (2020); Jian  _et al._ (2020); Turkeshi  _et al._ (2020); Bao
_et al._ (2020); Van Regemortel  _et al._ (2022); Sharma  _et al._ (2022);
Chen  _et al._ (2023); Cheng and Ippoliti (2023); Hauser  _et al._ (2023).
This association stems from the interpretation of the quantum trajectory
methods as simulations of the dynamics of (continuously) monitored quantum
systems Holland  _et al._ (1996); Gardiner and Zoller (2015); Fuji and Ashida
(2020).

In this paper, we expand upon the concept introduced in Ref. Vovk and Pichler
(2022), analyzing entanglement in various quantum trajectory schemes. We
complement our discussion with explicit examples of one-dimensional open
quantum dynamics, demonstrating that some trajectory-based methods employing
MPSs can yield an exponential reduction in classical computational cost
compared to other MPS trajectory-based methods or compared to conventional
matrix product operator (MPO) techniques. Finally, we note that our findings
are interesting not only from a computational point of view, but also from a
fundamental quantum-information-theoretic perspective. This follows the fact
that our analysis gives rise to heuristic algorithms for finding upper bounds
on mixed-state entanglement measures, such as the entanglement of formation
(EoF), a task that holds an independent and intrinsic interest Audenaert  _et
al._ (2001); Guevara and Viviescas (2014); Arceci  _et al._ (2022); Holevo
(2022, 2023); Carisch  _et al._ (2023).

This paper is structured as follows. In Section II we review basic concepts
and introduce the formalism this work is built upon. This includes a brief
review of the notion of entanglement in pure and mixed states, as well as the
discussion on stochastic propagation and quantum trajectories within the
formalism of quantum channels. In Section III we show the results obtained
using the concepts presented in Section II. These include a comparison between
various optimization strategies to simulate noisy quantum systems, as well as
the results of the numerical experiments on relevant many-body models, with
which we explore the practicality of developed ideas.

##  II Formalism and Methods

In this section we review important basic concepts, introduce the notation and
develop the formalism that we employ for the rest of the paper. In Subsection
II.1 we briefly review standard entanglement measures for pure states and
efficient tensor network representations thereof and consider the notions of
pure-state ensembles and ensemble-averaged functions of mixed states. In
Subsection II.2 we introduce the concept of a quantum channel as well as
variations of the quantum trajectory approach as a means to simulate quantum
channels. Then in Subsection II.3 we consider a special case of continuous
Markovian quantum channel, for which in Subsection II.4 we develop the concept
of an adaptive entanglement optimization.

###  II.1 Bipartite entanglement and pure-state decompositions

####  II.1.1 Pure states

Consider a pure state $|\psi\rangle$ of a one-dimensional111Throughout this
text we consider one-dimensional systems for simplicity. However, many
considerations can be generalized to higher dimensions. quantum many-body
system, consisting of $L$ subsystems with local Hilbert space dimension $d$.
Throughout this work we are interested in the entanglement between subsystem
$A$ and its complement $B$. Important objects for quantifying this bipartite
entanglement are the rank $\chi$ and the eigenvalues
$\\{\lambda_{1},\dots,\lambda_{\chi}\\}$ of the reduced state of either of the
partitions. For instance, the von Neumann entanglement entropy is given by
Nielsen and Chuang (2010):

$\displaystyle
E\left(|\psi\rangle\right)=-\sum_{i=1}^{\chi}\lambda_{i}\log_{2}\lambda_{i}.$
(1)

The rank $\chi$, also known as the bond dimension, is of central importance
for tensor-network-based methods, such as MPS techniques, as it governs the
computational cost of representing and processing such states on classical
computers Vidal (2003). The von Neumann entropy of the state (1) is often used
as a proxy for this cost. For a comprehensive review on MPSs we refer the
reader to Refs. Schollwöck (2011); Orús (2014).

####  II.1.2 Mixed states and pure-state decompositions

Imagine now that our system is in a mixed state given by the density matrix
$\rho$, a positive semi-definite Hermitian operator with a unit trace. This
operator can always be decomposed into an ensemble of pure states:

$\displaystyle\rho=\sum_{\alpha=1}^{k}p_{\alpha}|{\psi}_{\alpha}\rangle\langle{\psi}_{\alpha}|.$
(2)

The decomposition (2) represents $\rho$ as a statistical mixture of pure,
normalized states $|\psi_{\alpha}\rangle$ with probabilities $p_{\alpha}$. It
is also useful to introduce the unnormalized states
$|\tilde{\psi}_{\alpha}\rangle=\sqrt{p_{\alpha}}|\psi_{\alpha}\rangle$, such
that
$\rho=\sum_{\alpha=1}^{k}|\tilde{\psi}_{\alpha}\rangle\langle\tilde{\psi}_{\alpha}|$,
and the notation
$\tilde{\Psi}\equiv(|\tilde{\psi}_{1}\rangle,|\tilde{\psi}_{2}\rangle,\dots,|\tilde{\psi}_{k}\rangle)$
for a compact representation of the pure-state ensemble. Note that the
ensemble size $k$ is bounded from below by the rank of the density matrix,
$k\geq\mathrm{rank}\left(\rho\right)$.

It is important to notice that the decomposition (2) is not unique Hughston
_et al._ (1993). Two different ensembles of pure states give the same density
matrix

$\displaystyle\rho=\sum_{\alpha=1}^{k}|\tilde{\psi}_{\alpha}\rangle\langle\tilde{\psi}_{\alpha}|=\sum_{\alpha=1}^{r}|\tilde{\phi}_{\alpha}\rangle\langle\tilde{\phi}_{\alpha}|,$
(3)

if these two pure-state ensembles are isometrically related:

$\displaystyle|\tilde{\phi}_{\alpha}\rangle=\sum_{\beta=1}^{k}T_{\alpha,\beta}|\tilde{\psi}_{\beta}\rangle,$
(4)

where $T\in\mathcal{T}_{k}^{r}$ is an instance of an $r\times k$ isometry
(right unitary matrices) Audenaert  _et al._ (2001). We use the short-hand
notation $\tilde{\Phi}=T\tilde{\Psi}$ to write Eq. (4) more compactly.

Since the pure-state ensembles $\tilde{\Psi}$ and $\tilde{\Phi}$ are
equivalent in terms of their density matrix $\rho$, all ensemble-averaged
linear functions of the pure-state projectors, e.g., expectation values of
operators, coincide. However, the ensemble averages of nonlinear functions are
in general different from one ensemble to another. An important example is an
ensemble-averaged entanglement entropy (EAEE):

$\displaystyle\bar{E}[\tilde{\Psi}]=\sum_{\alpha=1}^{k}p_{\alpha}E\left(|{\psi}_{\alpha}\rangle\right),$
(5)

where $E(|{\psi}_{\alpha}\rangle)$ is the von Neumann entanglement entropy
(1). The dependence of the EAEE (5) on the pure-state ensemble $\tilde{\Psi}$
that decomposes $\rho$ in Eq. (2) is the key concept behind the mixed-state
entanglement measure known as the entanglement of formation (EoF) Bennett  _et
al._ (1996), which is the result of the minimization of the EAEE over all
possible ensembles of all possible sizes $r$:

$\displaystyle\min_{r}\inf_{T\in\mathcal{T}_{k}^{r}}\bar{E}[T\tilde{\Psi}]\equiv
E_{f}\left(\rho\right).$ (6)

According to Carathéodory’s theorem Rockafellar (1997); Uhlmann (1998, 2010),
the optimum (6) is attained by ensembles of size
$r\leq\mathrm{rank}^{2}\left(\rho\right)$. Computing the EoF for a generic
many-body quantum state is known to be NP-hard Gurvits (2003); Gharibian
(2008). Nevertheless, finding a pure-state decomposition with low EAEE is
practically useful, even if the global minimum can not be found. In
particular, in the case when the pure states from (2) are expressed as MPSs, a
lower EAEE (5) should allow for a more economic state representation.

###  II.2 Quantum channels and trajectories

We describe the dynamics of quantum systems using a quantum channel
$\mathcal{E}$ Nielsen and Chuang (2010), which is a completely positive trace-
preserving map that maps an input state of the system $\rho_{\mathrm{in}}$ to
an output state $\rho_{\mathrm{out}}$:

$\displaystyle\rho_{\mathrm{out}}=\mathcal{E}(\rho_{\mathrm{in}})=\sum_{\alpha=1}^{k}K_{\alpha}\rho_{\mathrm{in}}K_{\alpha}^{\dagger},$
(7)

where $K_{\alpha}$ are the Kraus operators that satisfy
$\sum_{\alpha=1}^{k}K^{\dagger}_{\alpha}K_{\alpha}=\mathbb{1}$. It is also
useful to introduce the notation
$\mathcal{K}=\left(K_{1},K_{2},\dots,K_{k}\right)$ as a compact way to denote
the Kraus representation of a quantum channel.

It is important to notice that the Kraus representation (7) is not unique. Two
different sets of Kraus operators represent the same quantum channel
$\mathcal{E}$:

$\displaystyle\rho_{\mathrm{out}}=\sum_{\alpha=1}^{k}K_{\alpha}\rho_{\mathrm{in}}K_{\alpha}^{\dagger}=\sum_{\alpha=1}^{r}R_{\alpha}\rho_{\mathrm{in}}R_{\alpha}^{\dagger},$
(8)

if they are isometrically related:

$\displaystyle R_{\alpha}=\sum_{\beta=1}^{k}T_{\alpha,\beta}K_{\beta},$ (9)

with $T\in\mathcal{T}_{k}^{r}$ being an isometric (right unitary)
transformation. We use the short-hand notation $\mathcal{R}=T\mathcal{K}$ to
write Eq. (9) more compactly. One can notice the analogy between the
transformation of the Kraus representations (9) and the transformation of the
ensemble states (4). As we will explore below, this transformation between
different Kraus representations is directly related to the possibility to
construct different stochastic propagation operators for a given quantum
channel.

Throughout this work we are interested in sequential quantum channels that are
composed of many channels applied in sequence. That is, if the system is
initially in a state $\rho_{0}$, its output state $\rho_{n}$ after a sequence
of $n$ quantum channels can be written as:

$\displaystyle\rho_{n}=\mathcal{E}^{(n)}(\rho_{n-1})=\mathcal{E}^{(n)}(\dots\mathcal{E}^{(i)}(\dots\mathcal{E}^{(1)}(\rho_{0}))),$
(10)

where $\mathcal{E}^{(i)}$ is the $i^{\mathrm{th}}$ quantum channel given by
the Kraus representation
$\mathcal{K}^{(i)}=(K_{1}^{(i)},\dots,K_{k_{i}}^{(i)})$. This model includes
two important examples. First, it can be used to describe discrete quantum
circuits, where a single quantum channel $\mathcal{E}^{(i)}$ represents the
the $i^{\mathrm{th}}$ layer of the circuit. Importantly, the quantum channel
description allows to include the case of imperfect quantum gates with
uncorrelated noise, an instance of which is considered in Subsubsection
III.3.1. Second, this model also includes continuous dynamics governed by
Markovian master equations, which is described in detail in Subsection II.3
and is used in the most parts of Section III. In the remainder of this
subsection we consider the connection between sequential quantum channels and
quantum trajectories.

####  II.2.1 Stochastic propagation and quantum trajectories

A sequential quantum channel implies a quantum dynamics that can be solved by
the quantum trajectory method Gardiner and Zoller (2015). In this method,
instead of keeping track of the full system state $\rho_{n}$ from one step to
the next, one sequentially samples pure states generated at each step and
recovers the solution by statistical averaging. To illustrate this, let us
consider a system initialized in a pure state
$\rho_{0}=|\psi\rangle\langle\psi|$. The state $\rho_{n}$ generated by the
sequence of quantum channels (10) can be expressed as

$\displaystyle\rho_{n}$
$\displaystyle=\sum_{\alpha_{n},\dots\alpha_{1}}K^{\left(n\right)}_{\alpha_{n}}\dots
K^{\left(1\right)}_{\alpha_{1}}|\psi\rangle\langle\psi|K^{\left(1\right){\dagger}}_{\alpha_{1}}\dots
K^{\left(n\right){\dagger}}_{\alpha_{n}}$ (11)
$\displaystyle=\sum_{\alpha_{n},\dots,\alpha_{1}}p_{\alpha_{n},\dots,\alpha_{1}}|\psi_{\alpha_{n},\dots,\alpha_{1}}\rangle\langle\psi_{\alpha_{n},\dots,\alpha_{1}}|,$
(12)

where $\alpha_{i}\in\left\\{1,...,k_{i}\right\\}$ for each
$i\in\left\\{1,...,n\right\\}$. In the spirit of decomposition (2), the state
(12) can be interpreted as a mixture of normalized pure states

$\displaystyle|\psi_{\alpha_{n},\dots,\alpha_{1}}\rangle\propto
K^{\left(n\right)}_{\alpha_{n}}\dots
K^{\left(1\right)}_{\alpha_{1}}|\psi\rangle$ (13)

with probabilities

$\displaystyle
p_{\alpha_{n},\dots\alpha_{1}}=\langle\psi|K^{\left(1\right){\dagger}}_{\alpha_{1}}\dots
K^{\left(n\right){\dagger}}_{\alpha_{n}}K^{\left(n\right)}_{\alpha_{n}}\dots
K^{\left(1\right)}_{\alpha_{1}}|\psi\rangle.$ (14)

The quantum trajectory approach is a method to sample the states (13) referred
to as quantum trajectories according to the probability distributions (14).
This is achieved by sequentially sampling pure states at each step of the
sequence (10). For instance, after application of the first quantum channel
from the sequence, $\mathcal{E}^{(1)}$, the full state of the system is
described by

$\displaystyle\rho_{1}=\mathcal{E}^{(1)}(\rho_{0})=\sum_{\alpha_{1}=1}^{k_{1}}K_{\alpha_{1}}^{\left(1\right)}|\psi\rangle\langle\psi|K_{\alpha_{1}}^{{\left(1\right)}\dagger},$
(15)

which again, as in (2), can be viewed as a pure-state decomposition:

$\displaystyle\rho_{1}=\sum_{\alpha_{1}=1}^{k_{1}}|\tilde{\psi}_{\alpha_{1}}\rangle\langle\tilde{\psi}_{\alpha_{1}}|=\sum_{\alpha_{1}=1}^{k_{1}}p_{\alpha_{1}}|{\psi}_{\alpha_{1}}\rangle\langle{\psi}_{\alpha_{1}}|,$
(16)

where
$|\tilde{\psi}_{\alpha_{1}}\rangle=K_{\alpha_{1}}^{\left(1\right)}|\psi\rangle$
are pure states with corresponding probability weights
$p_{\alpha_{1}}=\langle\tilde{\psi}_{\alpha_{1}}|\tilde{\psi}_{\alpha_{1}}\rangle$.
Using these weights, one can stochastically sample a single pure state from
the ensemble
$\tilde{\Psi}_{1}=\left(|\tilde{\psi}_{1}\rangle,...,|\tilde{\psi}_{k_{1}}\rangle\right)$,
thus stochastically propagating the state. That is, with probability
$p_{\alpha_{1}}$ one sets the quantum trajectory after the first step to be
$|{\psi}_{\alpha_{1}}\rangle$.

To continue the stochastic propagation of this trajectory, one executes the
same procedure for the next quantum channel,
$\mathcal{E}^{(2)}\left(|{\psi}_{\alpha_{1}}\rangle\langle{\psi}_{\alpha_{1}}|\right)$.
Iterating this stochastic algorithm for the whole sequence of channels
generates a single quantum trajectory after $n$ steps,
$|\psi_{\alpha_{n},\dots,\alpha_{1}}\rangle$, as given in Eq. (13). It is
important to observe that the probability to obtain this trajectory is given
by (14).

Stochastically generating $M$ independent trajectories according to the
algorithm described above, one recovers the state of the system $\rho_{n}$
from the following unbiased estimator:

$\displaystyle\rho_{n}=\lim_{M\rightarrow\infty}\frac{1}{M}\sum_{\ell=1}^{M}|\psi_{\ell}\rangle\langle\psi_{\ell}|,$
(17)

where $\ell=\left(\alpha_{n},\dots,\alpha_{1}\right)$ labels the sampled
trajectories. In the case when the initial state is mixed, one can simply
consider each pure state from the mixture separately and then take the
average.

####  II.2.2 Adaptive stochastic propagation

In the standard trajectory sampling method described above the Kraus
representations are fixed for all of the quantum channels and trajectories. In
a more general setting, one is allowed to change Kraus representations from
one trajectory to another. This simply follows from the fact that mixtures of
different pure-state decompositions of the same density operator also form
valid decompositions. We can further generalize this by choosing Kraus
representations in a way that depends on the state of current trajectory, that
is, in an adaptive way Van Dorsselaer and Nienhuis (1998). The corresponding
adaptive stochastic propagators generate valid quantum trajectories that
recover the density matrix via (17), as we discuss in the following.

To demonstrate the validity of adaptive trajectory sampling, let us consider
again the first step in the sequence (10). Following (9) we choose a
transformed Kraus representation
$R_{\alpha_{1}}^{(1)}=\sum_{\beta_{1}=1}^{k_{1}}T_{\alpha_{1},\beta_{1}}K_{\beta_{1}}^{(1)}$
to obtain:

$\displaystyle\rho_{1}=\sum_{\alpha_{1}=1}^{r}R^{(1)}_{\alpha_{1}}|\psi\rangle\langle\psi|R^{(1){\dagger}}_{\alpha_{1}}.$
(18)

Since the isometry $T$ is arbitrary, we can choose it in such a way that the
corresponding pure-state ensemble
$\tilde{\Phi}_{1}=(R^{(1)}_{1}|\psi\rangle,\dots,R^{(1)}_{r}|\psi\rangle)$ has
favourable properties. For example, as mentioned in Subsection II.1, the
states in an ensemble with a lower EAEE (5) can be represented more
efficiently with MPSs. Importantly, the transformation $T$ can depend on the
initial state $|\psi\rangle$. After the choice of $T$ is made, one samples a
state $|\tilde{\phi}_{\alpha_{1}}\rangle=R^{(1)}_{\alpha_{1}}|\psi\rangle$
with probability
$q_{\alpha_{1}}=\langle\tilde{\phi}_{\alpha_{1}}|\tilde{\phi}_{\alpha_{1}}\rangle$
and identifies it with the state of the trajectory after the first step.

To continue the propagation of the trajectory $|{\phi}_{\alpha_{1}}\rangle$,
one considers the next quantum channel,
$\mathcal{E}^{(2)}\left(|{\phi}_{\alpha_{1}}\rangle\langle{\phi}_{\alpha_{1}}|\right)$.
In the adaptive propagation scheme one chooses the Kraus representation of
$\mathcal{E}^{(2)}$ and the corresponding stochastic propagator based on the
state of the current trajectory $|{\phi}_{\alpha_{1}}\rangle$. To see that
this is valid and that the resulting trajectories faithfully recover the
density matrix (17), let us consider different Kraus representations of
$\mathcal{E}^{(2)}$ for each of the possible trajectories after the first
step. That is, we consider Kraus representations that explicitly depend on
$\alpha_{1}$:

$\displaystyle
R_{\alpha_{2}|\alpha_{1}}^{(2)}=\sum_{\beta_{2}}T_{\alpha_{2},\beta_{2}}^{\alpha_{1}}K_{\beta_{2}}^{(2)}.$
(19)

Here $T^{\alpha_{1}}$ is an isometric transformation that depends on the index
$\alpha_{1}$. It is easy to see that the state after the first two time steps
can be equivalently written as:

$\displaystyle\rho_{2}$
$\displaystyle=\sum_{\alpha_{2}\alpha_{1}}K^{(2)}_{\alpha_{2}}K^{(1)}_{\alpha_{1}}|\psi\rangle\langle\psi|K^{(1){\dagger}}_{\alpha_{1}}K^{(2){\dagger}}_{\alpha_{2}}$
(20)
$\displaystyle=\sum_{\alpha_{2}\alpha_{1}}R^{(2)}_{\alpha_{2}|\alpha_{1}}R^{(1)}_{\alpha_{1}}|\psi\rangle\langle\psi|R^{(1){\dagger}}_{\alpha_{1}}R^{(2){\dagger}}_{\alpha_{2}|\alpha_{1}}.$
(21)

Eq. (21) indeed shows that the adaptive stochastic propagation scheme
described above generates trajectories that faithfully recover the state of
the system from (17). This argument can be easily generalized to prove that an
$n$-step adaptive stochastic propagation scheme generates trajectories from a
valid decomposition of $\rho_{n}$ with the proper probability distributions
Vovk and Pichler (2022).

###  II.3 Continuous dynamics

In this work we are particularly interested in the continuous time dynamics of
Markovian quantum processes. If $t$ is the process duration, then the
corresponding Markovian quantum channel reads:

$\displaystyle\rho_{t}=\mathcal{E}_{t}\left(\rho_{0}\right)=e^{\mathcal{L}t}\rho_{0},$
(22)

where $\mathcal{L}$ is the generator of the system dynamics known as the
Lindbladian Lindblad (1976), which can be written as222Here we assume a time-
independent Lindbladian for notational simplicity. Generalizations to time-
dependent master equations are straightforward. In Section III the time-
dependence is indicated explicitly when necessary.

$\displaystyle\mathcal{L}~{}\boldsymbol{\cdot}=-i\left[H,\boldsymbol{\cdot}\right]+\sum_{j=1}^{m}\gamma_{j}\left(c_{j}\boldsymbol{\cdot}c_{j}^{\dagger}-\frac{1}{2}\left\\{c_{j}^{\dagger}c_{j},\boldsymbol{\cdot}\right\\}\right),$
(23)

where $H$ is a Hamiltonian, $c_{j}$ are arbitrary operators known as jump (or
Lindblad) operators, $\gamma_{j}\geq 0$ are the associated decoherence rates
and $m$ is a positive integer. In order to integrate the dynamical map (22),
it is often convenient to use the semigroup property and discretize the
quantum channel into small time steps $dt$ Gorini  _et al._ (1976):

$\displaystyle\rho_{t}=\mathcal{E}_{t}\left(\rho_{0}\right)=\underbrace{\mathcal{E}_{dt}(\mathcal{E}_{dt}...(\mathcal{E}_{dt}(}_{\text{$n$
times}}\rho_{0}))),$ (24)

with $t=ndt$ and $\mathcal{E}_{dt}=\exp\left(\mathcal{L}dt\right)$. Thus, a
Markovian quantum channel can be written as a sequential quantum channel (10)
composed of $n$ identical quantum channels
$\mathcal{E}^{(i)}=\mathcal{E}_{dt}~{}\forall i$. For sufficiently small $dt$
one can express $\mathcal{E}_{dt}$ with a representation consisting of $m+1$
Kraus operators:

$\displaystyle\left.\begin{array}[]{ll}K_{0}&=\mathbb{1}-i\left(H-\frac{i}{2}\sum_{j=1}^{m}\gamma_{j}c_{j}^{\dagger}c_{j}\right)dt,{}\\\
K_{j}&=\sqrt{\gamma_{j}dt}c_{j}\quad~{}\mathrm{for}~{}j\in\left[1,m\right].\end{array}\right.$
(27)

One can check that applying these Kraus operators in (7) with
$\alpha\in[0,~{}m]$ indeed recovers the action of the channel
$\mathcal{E}_{dt}$ up to $\mathcal{O}(dt^{2})$.

The discretized propagation explained above can be simplified even further by
Trotterization Trotter (1959) of each single-step quantum channel
$\mathcal{E}_{dt}$:

$\displaystyle\mathcal{E}_{dt}=\mathcal{E}_{dt}^{\mathrm{coh}}\mathcal{E}_{dt}^{\mathrm{incoh}},$
(28)

where the coherent channel $\mathcal{E}_{dt}^{\mathrm{coh}}$ is unitary and
thus has only one Kraus operator $K^{\mathrm{coh}}=\exp{\left(-iHdt\right)}$,
while the incoherent channel can be further split into a sequence of $m$
individual channels
$\mathcal{E}_{dt}^{\mathrm{incoh}}=\mathcal{E}_{dt}^{(1)}...\mathcal{E}_{dt}^{(j)}...\mathcal{E}_{dt}^{(m)}$
Daley (2014), where every channel $\mathcal{E}_{dt}^{(j)}$ can be represented
with two Kraus operators:

$\displaystyle\left.\begin{array}[]{ll}K^{(j)}_{0}&=\mathbb{1}-\left({\gamma_{j}dt}/{2}\right)c_{j}^{\dagger}c_{j},\\\
K^{(j)}_{1}&=\sqrt{\gamma_{j}dt}c_{j},\end{array}\right.$ (31)

that form the Kraus representation $\mathcal{K}^{(j)}_{dt}$. With the
decompositions performed above, the Markovian channel (22) is now expressed as
a sequential application of $n(m+1)$ simple channels, each of them given
either by a unitary channel $\mathcal{E}_{dt}^{\mathrm{coh}}$ of rank $k=1$ or
by a individual channel $\mathcal{E}^{(j)}_{dt}$ of rank $k=2$. Thus, the
application of any of these channels to a pure state results in a mixed state
whose rank is at most two. This simplicity allows for an optimization of the
adaptive stochastic propagation introduced in II.2.2, which we explore in the
next subsection.

###  II.4 Adaptive entanglement optimization

In this section we discuss strategies to adaptively choose propagators to find
trajectories with a low value of the EAEE (5). As outlined in Subsection II.1,
such pure-state ensembles are more economic in terms of their MPS
representations. Since finding the global EAEE optimum (6) is NP-hard Gurvits
(2003); Gharibian (2008), we devise simpler, heuristic strategies that find
ensembles with low (albeit not necessarily minimal) EAEE. Specifically, we
exploit the fact that we consider sequential quantum channel to design greedy
algorithms that break down the global optimization problem into a sequence of
simpler problems.

The central idea underlying our greedy optimization algorithms is to propagate
pure-state trajectories adaptively by choosing stochastic propagators for each
individual channel in such a way that the EAEE is minimised. Let us specify
this idea on the example of a continuous Markovian dynamics (22), which can be
broken down into a sequence of simple channels, as explained in Subsection
II.3. For each of these simple channels the pure state of the trajectory,
which we denote as $|\psi\rangle$, is stochastically propagated by Kraus
operators that are specified by an isometry $T$ (see Subsection II.2.2). We
thus choose $T$ such that the EAEE of the trajectory after this propagation
step is as small as possible. Therefore, in analogy with Eq. (6), we define
the following greedy entanglement optimization (GEO) problem333In what follows
we refrain from writing the individual decoherence channel index $j$, whenever
there is no potential for confusion.:

$\displaystyle\inf_{T\in\mathcal{T}_{2}^{4}}\bar{E}[T\tilde{\Psi}_{dt}],$ (32)

where the reference ensemble
$\tilde{\Psi}_{dt}=\mathcal{K}_{dt}|\psi\rangle=(K_{0}|\psi\rangle,K_{1}|\psi\rangle)$
is generated by the stochastic propagator defined in Eq. (31) and optimization
is performed over the isometries of the class $\mathcal{T}_{2}^{4}$, which,
according to Carathéodory’s theorem Rockafellar (1997); Uhlmann (1998, 2010),
is sufficient to attain the minimum.

Since we work with a Trotterization of a Markovian quantum channel, it is
natural to consider the limit of small time steps $dt$. This motivates the
definition of the following quantity:

$\displaystyle\lim_{dt\rightarrow
0}\inf_{T\in\mathcal{T}_{2}^{4}}\frac{\bar{E}[T\tilde{\Psi}_{dt}]-E\left(|\psi\rangle\right)}{dt}\equiv\dot{\bar{E}}_{\mathrm{GEO}}.$
(33)

Solving the minimization problem (32) at each step of the propagation sequence
thus defines an adaptive trajectory propagation algorithm. Note that each
state propagation is optimized independently, with the goal to minimize the
entanglement immediately after the propagation step. In this sense the
algorithm optimizes locally and in a greedy way. Below we consider several
simplifications that are obtained by replacing the cost function (32) with
proxies that are numerically more tractable.

A way to simplify the optimization (32) is to consider a limited class of
isometries to optimize over. The minimal non-trivial class of isometries is
$\mathcal{T}^{2}_{2}\subset\mathcal{T}^{4}_{2}$, the instances of which are
$2\times 2$ matrices $T(\theta,\varphi)$ that can be written as Hedemann
(2013):

$\displaystyle
T(\theta,\varphi)=\begin{pmatrix}\cos(\theta)&&\sin(\theta)e^{-i\varphi}\\\
\sin(\theta)&&-\cos(\theta)e^{-i\varphi}\\\ \end{pmatrix}.$ (34)

With this transformation the optimization reduces to a simple two-parameter
optimization over $\theta$ and $\varphi$. Clearly, the result of this
optimization gives an upper bound on (33). In the following we refer to this
optimization simply as 2-GEO. The Supplementary Materials (SM) I and II
contain the details regarding the numerical implementations of the
optimization method.

An even simpler algorithm can be obtained by considering the time derivative
of the EAEE for a given isometry $T$:

$\displaystyle\dot{\bar{E}}\left(T\right)\equiv\lim_{dt\rightarrow
0}\frac{\bar{E}[T\tilde{\Psi}_{dt}]-E\left(|\psi\rangle\right)}{dt}.$ (35)

Minimizing this time derivative with respect to $T$ defines an alternative
greedy algorithm to construct trajectories with small average entanglement,
which was explored in Ref. Vovk and Pichler (2022). Note that the
corresponding greedy entanglement derivative optimization (GEDO) problem

$\displaystyle\inf_{T\in\mathcal{T}_{2}^{4}}\dot{\bar{E}}\left(T\right)\equiv\dot{\bar{E}}_{\mathrm{GEDO}},$
(36)

is obtained formally by exchanging the limit and minimization in Eq. (33),
which implies that the quantity (36) is an upper bound to the minimum (33):

$\displaystyle\dot{\bar{E}}_{\mathrm{GEDO}}\geq\dot{\bar{E}}_{\mathrm{GEO}}.$
(37)

Even though GEO (33) can give trajectories with lower average entanglement,
GEDO (36) has a key feature: $\dot{\bar{E}}_{\mathrm{GEDO}}$ can be calculated
efficiently when the state $|\psi\rangle$ is represented by an MPS. The reason
for this is twofold. First, in the GEDO the minimum is attainable with
isometries from $\mathcal{T}^{2}_{2}\subset\mathcal{T}^{4}_{2}$. Second, the
minimization landscape in this case is trivial and can be explored
analytically. Specifically, the minimization reduces to choosing between the
transform parameters $\theta=0$ and $\theta=\pi/4$ and evaluating the optimal
phase parameter $\varphi$ (see Vovk and Pichler (2022) for details). The
evaluation of the analytic expressions of GEDO has a cost of
$\mathcal{O}(\chi^{3}d)$, where $\chi$ is the bond dimension and $d$ is the
local Hilbert space dimension. In contrast, calculating
$\dot{\bar{E}}_{\mathrm{GEO}}$ in general is very costly, since it requires
gradient-descent-type numerical optimization, for which we need a
normalization sweep at every optimizer iteration. This operation has a cost of
$\mathcal{O}(\chi^{3}d^{3}L\kappa)$, where $L$ is the system length and
$\kappa$ is the number of gradient descent iterations (see SM I for technical
details).

##  III Results

In this section we analyse three distinct aspects of the quantum trajectory
algorithms introduced above. First, in Subsection III.1, we apply the adaptive
unraveling schemes to simple two-qubit systems. This serves to illustrate
differences between optimization schemes as well as their greedy nature.
Second, in Subsection III.2, we show that different unraveling schemes can
lead to qualitatively different behaviour in the entanglement of many-body
systems. For this we consider noisy random Brownian circuits and determine the
parameter regimes that correspond to different entanglement scaling laws. We
locate the transition between area- and volume-law phases and show that the
transition points depend on the unraveling scheme. Third, in Subsection III.3,
we compare trajectory-based MPS calculations with standard MPO-based approach
and show that the former can be more efficient. For this we construct an
explicit example of a noisy many-body models, that can be unraveled using MPSs
with a finite bond dimension, while the MPO solution of the same model
requires exponential bond dimensions.

###  III.1 Comparison of optimization strategies

####  III.1.1 Minimal example: Bell state dephasing

Figure 1: Comparison of optimization strategies. (a) The entanglement
probability density distributions over $M=1000$ trajectories for a two-qubit
dephasing process with $c_{j}=\sigma_{j}^{z}/2$. The colour plot corresponds
to the distribution of trajectories generated by the GEDO method (36) with the
average indicated by the purple line. The solid black line corresponds to the
distribution of the trajectories generated by the 2-GEO method (33). The
entanglement entropy obtained by the latter method is equal to the EoF (6) for
all the trajectories, therefore its distribution looks like a line. (b) The
EAEE change rates obtained by the 2-GEO vs. GEDO method for a single-channel
dephasing of two qubits with $c=\sigma^{z}/2$ and $\gamma=1$. Each dot
represents an instance of a Haar random state (total of 20000 instances), the
black rims around the dots indicate instances from Eq. (38) (total of 200
instances). The inset shows the EoF vs. 2-GEO EAEE change rates obtained for
the same conditions. The color of the dots corresponds to the initial state
entanglement, indicated by the color bars on the right. (c) The EAEE change
rates of the 2-GEO (solid colored lines) and GEDO (dashed black line) methods
vs. isometry parameter $\theta$ from Eq. (34). In contrast to the 2-GEO
method, the optimization landscape is constant in $\theta$ except the point
$\theta=0$. The inset shows the dependence of the optimal $\theta$ of the
2-GEO method on the value of the time step $dt$, which is in agreement with
the analytic formula (39). (d) The two-qubit EoF vs. the 2-GEO EAEE change
rates in the case of two incoherent channels with $c_{j}=\sigma_{j}^{z}/2$ and
$\gamma_{j}=1$ ($j=1,2$). Each dot represents an instance of a Haar random
state (total of 20000 instances). The color of the dots corresponds to the
initial state entanglement, indicated by a color bar on the right.

Let us start with a minimal example of two qubits initially in a Bell state,
$|\psi\rangle=|\Phi^{+}\rangle=(|00\rangle+|11\rangle)/\sqrt{2}$. We consider
a situation where both qubits dephase with a rate $\gamma$. That is, their
dynamics is described by a Lindblad master equation (23) with $H=0$ and
$c_{j}=\sigma_{j}^{z}/2$ for $j=1,2$. The corresponding Kraus operators, as
defined in Eq. (31), generate trajectories that jump between two Bell states
$|\Phi^{+}\rangle$ and $|\Phi^{-}\rangle=(|00\rangle-|11\rangle)/\sqrt{2}$,
which are both maximally entangled. Hence, the ensemble generated by such an
unraveling has the maximal EAEE at all times. In contrast, if we transform the
Kraus representation according to Eq. (9) with
$T(\theta,\varphi)\in\mathcal{T}_{2}^{2}$, then it is easy to see that the
states in all trajectories are of the form

$\displaystyle|\psi(x,\nu)\rangle=\cos(x)|11\rangle+\sin(x)e^{i\nu}|00\rangle,$
(38)

where $x\in\left[0,\pi/2\right]$ and $\nu\in\left[0,2\pi\right)$. For these
states the optimization over $T(\theta,\varphi)$, as defined in Eq. (32), can
be explicitly performed. The optimal $\theta$ depends on the state of the
trajectory (38) and is given by

$\displaystyle\theta(x,dt)=\arcsin\left[\sec\left(2x\right)\sqrt{\gamma
dt}/2\right],$ (39)

and the optimal $\varphi=0$. It is interesting to observe that all the
trajectories obtained by this adaptively optimized unraveling have the same
entanglement, which matches the EoF during the whole evolution time Wootters
(1998b), as indicated in Fig. 1(a) by the black line. We note that this is
specific to this minimal example and is not the case in general. In Fig. 1(a)
we also show the entanglement probability density distribution $p(E)$ for
trajectories obtained if we instead minimize the EAEE derivative, given in Eq.
(36). In this case the entanglement distribution initially spreads and drifts
to smaller values before accumulating around zero for longer times. As
expected from Eq. (37), in this case the mean entanglement, indicated in Fig.
1(a) by the dark purple line, is larger than the EoF.

####  III.1.2 Numerical analysis of single-step optimization strategies

In this subsection we consider performance differences and features of the
optimization strategies presented above. For this we take several typical
states, let them evolve under a decoherence channel for a single time step and
compare the resulting EAEE change rates for different channel unravelings.
Specifically, we take Haar-random two-qubit states Hayden  _et al._ (2006) and
compare the average entanglement after evolution with 2-GEO and GEDO
unravelings of a single decoherence channel. In Fig. 1(b) we make such a
comparison for the jump operator $c=\sigma^{z}/2$. We can see that the 2-GEO
method that directly optimizes the entanglement gives lower EAEE values than
the GEDO approach that optimizes the entanglement derivative, which is
consistent with the inequality (37). We note that the EAEE change rate is
larger for the states with larger initial entanglement. Interestingly, the
states taken from (38) highlighted by black rims in Fig. 1(b) give the largest
difference for the EAEE change rate. In addition, as shown in the inset of
Fig. 1(b), we also observe that, for the case of two qubits with a single
decoherence channel, the 2-GEO method (32) is sufficient to reach the EoF
Wootters (1998b), which implies that in this case going from GEO to 2-GEO is
not a limitation. In Section III of the SM we provide the same analysis for
states of more qubits.

To gain more insights into the difference between the two optimization
methods, it is useful to consider the optimization landscapes in both cases.
For this it is important to note that the optimization landscape of the 2-GEO
method (33) depends on the time step $dt$. In contrast, the GEDO approach (36)
does not depend on $dt$ due to the fact that the limit of $dt\rightarrow 0$ is
taken before the minimization. To illustrate this, in Fig. 1(c) we take a
particular instance of the state (38) with $\cos(x)=0.98$ as an initial state
and consider the EAEE optimization landscapes444We consider specific cuts of
optimization landscapes with a fixed phase parameter at $\varphi=0$, since in
the considered case this contains the global minimum. for different values of
the time step $dt$. Notably, the minimal EAEE values differ for the 2-GEO and
GEDO methods. In particular, the minimal value of
$\dot{\bar{E}}_{\text{2-GEO}}$ is below $\dot{\bar{E}}_{\mathrm{GEDO}}$. This
highlights the fact that the the optimization and the limit $dt\rightarrow 0$
do not commute. Additionally, we observe that the GEDO optimization landscape
is trivial, as the optimization landscape is constant for every value of
$\theta$ except when $\theta=0$, and the minimum in this case is obtained by
comparing the cases of $\theta=0$ and $\theta=\pi/4$ (see also Ref. Vovk and
Pichler (2022)).

####  III.1.3 Two-channel propagation analysis

An important aspect of the adaptive propagation strategies discussed in
Subsection II.4 is their greedy nature, i.e., the fact that the optimization
is done separately for each decoherence channel (31) in the (Trotter)
sequence. This becomes crucial in the case of multiple decoherence channels,
where optimizing each of them separately does not necessarily produce the
ensemble with minimal EAEE. We illustrate this again with an example of two
independently dephasing qubits. For this setup we compare the entanglement
obtained from sequentially optimized unraveling to the global optimum, which
is given by the EoF Wootters (1998b). Specifically, we evolve the two qubits
for a single time step $dt$, which consists of a sequential application of a
dephasing channels to each of the qubits. That is, we first apply a dephasing
channel to the first qubit and then propagate the resulting state with a
dephasing channel that acts on the second qubit. Each propagation is optimized
independently according to the 2-GEO method. From the resulting states we
extract the change rate of the EAEE, which we plot against the entanglement of
formation of the mixed quantum state in Fig. 1(d). As expected, the sequential
optimization of the two channels does not give the minimal configuration.

###  III.2 Entanglement in noisy random Brownian circuits

So far we only considered two-qubit systems and evolution for a single time
step. For larger system sizes and longer evolution times more interesting
phenomena can occur. The competition between unitary dynamics and incoherent
processes can give rise to distinct phases, characterized by different
entanglement properties of the trajectories. If the unitary dynamics dominates
in a given unraveling scheme, the entanglement of the trajectories can grow
extensively with the system size (volume-law phase). In contrast, when the
stochastic processes dominate, the entanglement can be bounded by the size of
the subsystem boundary (area-law phase) Eisert and Osborne (2006); Bao  _et
al._ (2020). This transition, typically referred to as measurement-induced
phase transition, has been observed in several models in recent studies Li
_et al._ (2018); Skinner  _et al._ (2019); Chan  _et al._ (2019); Li  _et al._
(2019); Choi  _et al._ (2020); Zabalo  _et al._ (2020); Jian  _et al._ (2020);
Turkeshi  _et al._ (2020); Bao  _et al._ (2020); Van Regemortel  _et al._
(2022); Sharma  _et al._ (2022); Chen  _et al._ (2023); Cheng and Ippoliti
(2023). Importantly, as shown in Ref. Vovk and Pichler (2022), the transition
point can change for different unraveling schemes. Due to the connection
between the entanglement and computational cost for classical simulations with
tensor networks, different trajectory propagation approaches can thus differ
significantly in terms of their computational cost for solving the same master
equation.

Figure 2: Entanglement distributions in various unravelings of the continuous
RBC model with dephasing. (a) The initial state $|\psi\rangle=|11\dots
1\rangle$ is propagated by an open continuous RBC consisting of alternating
Trotter layers corresponding to coherent and incoherent part of the master
equation. For each incoherence channel $\mathcal{E}_{dt}$ (the channel index
$j$ is omitted) we consider the “quantum jump” unraveling that uses the
standard Kraus representation $\mathcal{K}_{dt}$, “quantum state diffusion”
unraveling that uses the transformed Kraus representation $\mathcal{R}_{dt}$
(see text), the GEDO method (36) and the 2-GEO method (33). Panels (b) – (e)
show the corresponding entanglement probability density distributions of the
entanglement exponent $p\left(e^{E}\right)$ in the trajectories modelling the
RBC dephasing process. For panels (c), (d) and (e) we additionally plot the
average of the entanglement exponent on top of the distributions (black
lines). The inset in (d) illustrates the average proportion of choices between
cases $\theta=0$ and $\theta=\pi/4$ done by the GEDO algorithm (see text). The
number of qubits in the system is $L=16$, the jump operators are
$c_{j}=\sigma_{j}^{z}/2$ and the decoherence rate relative to the variance of
the Gaussian variable in the RBC (40) is $\gamma=10$ ($\gamma_{j}=\gamma$).
The number of trajectories with various RBC realizations is $M=200$ for each
panel, the bond dimension is $\chi=256$, which is exact for $L=16$. We
emphasize a different vertical scale in the panel (b) compared to the panels
(c), (d) and (e).

To see this, let us consider the model studied in Ref. Vovk and Pichler
(2022), i.e., the noisy Random Brownian circuit (RBC) that is schematically
depicted in Fig. 2(a). It describes a one-dimensional chain of spin-$1/2$
systems, whose coherent evolution is given by the Hamiltonian:

$\displaystyle
H(t)=\sum_{j=1}^{L-1}\sum_{o,p=0}^{3}g^{o,p}_{j}(t)~{}\sigma_{j}^{o}\otimes\sigma_{j+1}^{p},$
(40)

where
$\sigma_{j}^{o}\in\left\\{\mathbb{1}_{j},\sigma_{j}^{x},\sigma_{j}^{y},\sigma_{j}^{z}\right\\}$
are the Pauli matrices acting on qubit $j$. The parameters $g^{o,p}_{j}$ are
Gaussian stochastic variables with $\left<\left<g^{o,p}_{j}\right>\right>=0$
and
$\left<\left<g^{o,p}_{j}g^{o^{\prime},p^{\prime}}_{j^{\prime}}\right>\right>=\alpha\delta_{j,j^{\prime}}\delta_{o,o^{\prime}}\delta_{p,p^{\prime}}$,
where $\left<\left<...\right>\right>$ denotes the average over Hamiltonian
realizations and $\alpha$ is the variance of the Gaussian variable555Without
the loss of generality, in what follows we set $\alpha=1$.. The incoherent
part of the evolution is described by jump operators $c_{j}=\sigma_{j}^{z}/2$
with $\gamma_{j}=\gamma$.

In panels (b – e) of Fig. 2 we show the histograms of the entanglement in the
various trajectories that solve the dynamics of the noisy RBC. Specifically,
we plot the exponential $e^{E}$ of the entanglement entropy (1) of every
trajectory, since this quantity serves as a proxy for the bond dimension
$\chi$ required to represent the trajectory in the MPS form Vidal (2003). In
panel (b) of Fig. 2 we depict the histogram of trajectories generated by the
Kraus representation $\mathcal{K}_{dt}$ defined in Eq. (31), which corresponds
to the standard quantum trajectory method known as the “quantum jump” approach
Zoller  _et al._ (1987); Dum  _et al._ (1992); Dalibard  _et al._ (1992);
Gardiner  _et al._ (1992); Castin  _et al._ (1993); Plenio and Knight (1998).
In this case, both of the propagation operators from Eq. (31) are local
unitary operations, and thus the incoherent part of the evolution does not
remove any entanglement, despite the large decoherence rate ($\gamma=10$). The
dynamics of the entanglement is thus determined by the coherent part, which
leads generically to linear entanglement growth. This growth continues until a
saturation value proportional to the system size is reached. In Fig. 2(c) we
plot the histogram generated by the Kraus representation
$\mathcal{R}_{dt}=T(\theta,\varphi)\mathcal{K}_{dt}$, where
$T(\theta,\varphi)$ defined in Eq. (34) has fixed parameters $\theta=\pi/4$
and $\varphi=0$. This unraveling is known as the “quantum state diffusion”
Gisin and Percival (1992); Tian and Carmichael (1992); Van Dorsselaer and
Nienhuis (1998). In contrast to the “quantum jump” method, the “quantum state
diffusion” approach results in entanglement distributions that saturate at
much smaller values. Importantly, these values here are not determined by the
system size, but by the decoherence rate (see also Ref. Vovk and Pichler
(2022)). In Fig. 2(d) we depict the histogram obtained by the GEDO method,
which minimizes the EAEE change rate (36) by switching between transformed
Kraus representations with $\theta=0$ and $\theta=\pi/4$ (with free phase
parameter $\varphi$) and selecting the best one at every propagation step.
Comparing panels (c) and (d) we observe a similarity in the trajectory
distributions. This is due to the fact that the GEDO approach in this case
almost exclusively selects $\theta=\pi/4$ with $\varphi=0$, which corresponds
to the “quantum state diffusion” unraveling. This can be also seen from the
inset in Fig. 2(d), where we show the average proportion of $\theta$ parameter
choices made by the GEDO algorithm666As in the example in Fig. 1(c), in this
case the GEDO minimum is found always at $\varphi=0$, therefore this parameter
is not shown.. At last, in Fig. 2(e) we plot the histogram obtained by the
2-GEO method (32), which minimizes the EAEE by adaptively choosing the
transform $T(\theta,\varphi)$. In this case the parameters $\theta$ and
$\varphi$ are determined numerically and the optimization cannot be reduced to
switching between two cases as it happens in the GEDO method. Among the four
unraveling schemes shown in Fig. 2, the 2-GEO one results in the trajectories
with the smallest average entanglement and also with the smallest weight on
trajectories with higher entanglement.

Figure 3: Determination of the critical decoherence rates $\gamma^{*}$ in the
continuous RBC model with dephasing for various unravelings. (a) The
dependence of the distributions of the steady-state EAEE
$\langle\langle\bar{E}_{\mathrm{SS}}^{L_{A}}\rangle\rangle$ over the bipartite
cut location $L_{A}$ on the decoherence rate $\gamma$ for various system sizes
$L$ and decoherence channel representations $\tilde{\mathcal{K}}_{dt}$,
$\mathcal{R}_{dt}$ and $\tilde{\mathcal{R}}_{dt}$ (see text for more details).
The arrows indicate the increasing values of $\gamma$ different for each of
the Kraus representations: $\gamma=\left\\{9,~{}10,~{}11,~{}12\right\\}$ for
$\tilde{\mathcal{K}}_{dt}$, $\gamma=\left\\{7,~{}8,~{}9,~{}10\right\\}$ for
$\mathcal{R}_{dt}$ and $\gamma=\left\\{14,~{}16,~{}18,~{}20\right\\}$ for
$\tilde{\mathcal{R}}_{dt}$. (b) The dependence of the middle-cut steady-state
EAEE divided by the system size
$\langle\langle\bar{E}_{\mathrm{SS}}^{L/2}\rangle\rangle/L$ on the inverse
system size $1/L$. According to the ansatz (46), the slope of the plot
corresponds to the area-law coefficient $a(\gamma)$. The arrows indicate the
increasing values of $\gamma$. (c) The dependence of the inverse values of the
area-law coefficients $1/a(\gamma)$ on the decoherence rate $\gamma$. The data
is compatible with linear behaviour, therefore we use the linear extrapolation
to determine the critical decoherence rates $\gamma^{*}$ written in
corresponding colors in the plot. The grey shadings indicate the standard
errors of the numerical lines, while the dashed lines illustrate the
extrapolations. The system sizes used to generate the data:
$L=\left\\{44,~{}48,~{}52,~{}56,~{}60,~{}64,~{}68,~{}72\right\\}$, the number
of trajectories with various RBC realizations is $M=200$ and the bond
dimension is $\chi=300$. The bond dimension convergence is checked against
instances with $\chi=512$.

Importantly, the entanglement distributions do not only differ quantitatively
across unraveling schemes, but they also change qualitatively at long times
Vovk and Pichler (2022). In particular, for the unraveling in Fig. 2(b) the
entanglement distribution concentrates around a saturation value that is
proportional to the system size. In contrast for unravelings in panels (c – e)
of Fig. 2, the saturation value is independent on the system size, indicating
that the system is in the area-law phase. This area-law phase is expected to
persist for sufficiently large decoherence rate $\gamma$. At the critical
value $\gamma=\gamma^{\ast}$, one expects a so-called measurement-induced
phase transition from area- to volume-law entanglement. Hence, interesting
quantities to determine are the critical decoherence rates $\gamma^{*}$ for
various unraveling schemes.

In the following we focus on determining the critical points in unravelings
congenial to some of the ones presented in Fig. 2. The “quantum jump”
unraveling shown in Fig. 2(b) demonstrates a volume-law scaling and hence does
not possess a critical point. We therefore consider a variation of “quantum
jump” unraveling with the following Kraus representation:

$\displaystyle\left.\begin{array}[]{ll}\tilde{K}^{(j)}_{0}&=\mathbb{1}-\left({\gamma
dt}/{2}\right)|1\rangle_{j}\langle 1|,\\\ \tilde{K}^{(j)}_{1}&=\sqrt{\gamma
dt}|1\rangle_{j}\langle 1|,\end{array}\right.$ (43)

to which we refer to as $\tilde{\mathcal{K}}_{dt}$ for brevity777This
representation can be obtained by applying the isometric transformation (34)
with parameters $\theta\left(dt\right)=\arcsin\left(\sqrt{\gamma dt}/2\right)$
and $\varphi=0$ to the original representation. Note that in this case the
transformation is time-step-dependent.. We also study “quantum-state-
diffusion” unravelings:

$\displaystyle\left.\begin{array}[]{ll}R^{(j)}_{\pm}&=\frac{1}{\sqrt{2}}\left[\left(1-{\gamma
dt}/{8}\right)\mathbb{1}_{j}\pm e^{-i\varphi}\sqrt{\gamma
dt}\sigma_{j}^{z}/2\right],\end{array}\right.$ (45)

of which we take two instances with $\varphi=0$ and $\varphi=\pi/4$ and
consider how changing the phase parameter affects the critical point. We refer
to these representations as $\mathcal{R}_{dt}$ and $\tilde{\mathcal{R}}_{dt}$,
respectively.

In what follows we do not consider adaptive unravelings shown in Fig. 2(d) and
(e). The reason for this is twofold. First, as already mentioned, in the case
of the open RBC the GEDO approach chooses mostly the “quantum state diffusion”
unraveling. We therefore expect this method to have similar critical
decoherence rate. Second, as mentioned in the end of Subsection II.4, our
implementation of the 2-GEO unraveling is too computationally demanding for
large system sizes that are required to establish the critical point (see SM I
for details).

Figure 4: Discrete noisy quantum circuit and entanglement dynamics therein for
different simulation methods. (a) Discrete quantum SWAP circuit with dephasing
applied to infinite chain of Bell pairs (47). Each discrete circuit step $n$
consists of nearest-neighbour SWAP gates applied alternatively to odd and even
neighbouring qubit pairs and individual decoherence channels
$\mathcal{E}^{(j)}$ defined in Eq. (48). The system’s bipartition $AB$ is
shown by the grey dashed line. (b) The (operator) entanglement dynamics during
the discrete noisy circuit evolution with various decoherence rates $p$. The
lines’ colors that change from bright orange to dark purple correspond to
different values of the decoherence rate,
$p=0,~{}0.005,~{}0.01,~{}0.02,~{}0.05$, and $0.1$. Colored solid lines
indicate the operator entanglement (49) and the black dotted line shows the
EAEE dynamics for the standard Kraus representation $\mathcal{K}^{(j)}$ from
Eq. (53) that is independent on $p$. Both of these cases exhibit an unbounded
entanglement growth. Colored solid lines with circles indicate the EAEE (57)
of the trajectories generated with the transformed Kraus representation
$\mathcal{R}^{(j)}$, while colored solid lines with crosses correspond to EAEE
of the 2-GEO method (total of $M=50$ trajectories). The colored dashed lines
illustrate the EoF minimum (6).

To determine the critical decoherence rates $\gamma^{*}$ for different
unravelings, we use a finite-size scaling argument. We consider the
entanglement entropy averaged over trajectories and RBC realizations,
$\langle\langle\bar{E}^{L_{A}}_{\mathrm{SS}}\rangle\rangle$, where $L_{A}$ is
the size of the subsystem. The subscript SS indicates that the value is
evaluated at late times, when the entanglement distribution becomes
stationary. When $L_{A}=L/2$, the above quantity is constant in the area-law
phase and grows linearly with the system size $L$ in the case of the volume-
law scaling. Thus one can make the following ansatz:

$\displaystyle\langle\langle\bar{E}^{L/2}_{\mathrm{SS}}\rangle\rangle=a\left(\gamma\right)+b\left(\gamma\right)L,$
(46)

where $a\left(\gamma\right)$ is the EAEE saturation value in the area-law
phase and $b\left(\gamma\right)$ is the slope of the volume-law scaling. The
area-law phase is thus identified when $b\left(\gamma\right)=0$ and
$a\left(\gamma\right)$ is finite. When $\gamma$ approaches the critical point
$\gamma^{*}$ from above (i.e., from the area-law phase), we expect
$a\left(\gamma\right)$ to diverge, while $b\left(\gamma\right)$ should remain
zero. In Fig. 3 we use this to locate the critical point $\gamma^{*}$.
Specifically, we start with the system initially in the state
$|\psi\rangle=|11\dots 1\rangle$, which we evolve with the random circuit (40)
with dephasing channels as depicted in Fig. 2(a). The dephasing channels in
the circuit are expressed by the mentioned representations
$\tilde{\mathcal{K}}_{dt}$, $\mathcal{R}_{dt}$ and $\tilde{\mathcal{R}}_{dt}$,
which result in different trajectory unravelings.

In Fig. 3(a) we plot the EAEE in the system after sufficiently long evolution
time, $\langle\langle\bar{E}^{L_{A}}_{\mathrm{SS}}\rangle\rangle$, for various
system sizes, unraveling schemes and decoherence rates. Using this data, in
Fig. 3(b) we plot the middle-cut steady-state EAEE against the inverse of the
system size, from which, using the ansatz (46), we can extract the area-law
coefficients $a\left(\gamma\right)$ and confirm that the volume-law
coefficient $b\left(\gamma\right)=0$ within the fitting error. In Fig. 3(c) we
plot the inverse of area-law coefficients against the decoherence rate and
extrapolate this data to get the critical rates $\gamma^{*}$.

From the results obtained in Fig. 3(c) we observe that the critical points
indeed depend on the particular unraveling of the decoherence channels.
Specifically, unravelings with Kraus representations
$\tilde{\mathcal{K}}_{dt}$ and $\mathcal{R}_{dt}$ give lower critical values
compared to the unraveling with $\tilde{\mathcal{R}}_{dt}$. Another
observation stems from comparing the “quantum state diffusion” unravelings
$\mathcal{R}_{dt}$ and $\tilde{\mathcal{R}}_{dt}$, where different values of
the phase parameter $\varphi$ significantly affect the location of the
critical point. Specifically, within the standard error the obtained critical
decoherence rates are related to each other as
$\gamma^{*}|_{\mathcal{R}_{dt}}=\cos^{2}\left(\pi/4\right)\gamma^{*}|_{\tilde{\mathcal{R}}_{dt}}$,
as predicted in Ref. Vovk and Pichler (2022).

###  III.3 Comparison with other approaches to solve the many-body dynamics

We study entanglement properties in various quantum trajectory unravelings,
because this is crucial for determining the computational cost of solving the
master equation (23) using quantum trajectories represented with MPSs (which
we refer to as QT MPS in the following). Given this motivation, it is
important to compare this approach with other existing methods to solve master
equations. Specifically, in this subsection we compare the QT MPS method with
an MPO approach Bonnes and Läuchli (2014); Wellnitz  _et al._ (2022); Preisser
_et al._ (2023). The latter consists in vectorizing the density matrix and
corresponding equations of motion (23) and representing them as tensors. In
general, it is unclear which of the two methods, QT MPS or MPO, is
computationally less expensive. The computational cost of a simulation based
on MPO scales as $\mathcal{O}\left(\chi^{3}d^{6}Ln\right)$, where $\chi$ is
the bond dimension, $d$ is the local Hilbert space dimension, $L$ is the
system length and $n$ is the total number of (continuous) time steps. In turn,
for the QT MPS method (with a fixed statistical accuracy) the computational
cost scales as $\mathcal{O}\left(\chi^{3}d^{3}Ln\right)$. Note that the bond
dimension $\chi$ has different physical meaning for the MPS and MPO. As
discussed above, in the MPS approach the bond dimension $\chi$ is related to
the entanglement in individual trajectories. On the other hand, for the MPO
the bond dimension represents both quantum and classical correlations in the
density matrix. To illustrate this difference we consider a simple model of a
system that features unbounded growth of classical correlations during its
evolution, while the amount of entanglement in individual trajectories can
remain finite.

####  III.3.1 Discrete case

Let us consider an infinite one-dimensional chain of qubits being initially a
product of nearest-neighbour Bell pairs:

$\displaystyle|\psi\rangle=\dots\otimes|\Phi^{+}\rangle\otimes|\Phi^{+}\rangle\otimes|\Phi^{+}\rangle\otimes|\Phi^{+}\rangle\otimes\dots.$
(47)

We consider time evolution of this initial state under the discrete circuit
depicted in Fig. 4(a), which consists of alternating layers of nearest-
neighbor SWAP gates Nielsen and Chuang (2010) and individual local dephasing
channels of the form:

$\displaystyle\mathcal{E}^{(j)}\left(\rho\right)=\left(1-p\right)\rho+p\sigma_{j}^{z}\rho\sigma_{j}^{z},$
(48)

where $p$ is the decoherence strength.

To understand the dynamics of this model, first note that the coherent part of
the evolution (induced by the SWAP gates) simply leads to a ballistic
separation of the the initially correlated nearest-neighbor pairs across the
chain, with the distance between each partner of a pair growing as $4n$. The
incoherent part does not affect this ballistic dynamics, but instead leads to
a dephasing of each pair, which transforms entanglement in each pair into
classical correlations. The dynamics thus simply reduces to that of
independent pairs, and measures of entanglement or correlations across a given
bipartition (and the associated computational costs of MPS or MPS
representations) can thus be evaluated analytically.

As a proxy of the computational cost of the MPO method we consider its
operator entanglement entropy Wellnitz  _et al._ (2022):

$\displaystyle
E_{\mathrm{MPO}}=-\sum_{i}\mu^{2}_{i}\log_{2}\mu^{2}_{i},~{}\text{with}~{}\mu_{i}=\frac{\tilde{\mu}_{i}}{\sqrt{\sum_{i^{\prime}}\tilde{\mu}_{i^{\prime}}^{2}}},$
(49)

where $\tilde{\mu}_{i}$ ($\mu_{i}$) are the (normalized) singular values of
the MPO’s middle cut. In the above setting the corresponding operator
entanglement grows with the circuit depth $n$, and has two characteristic
regimes. Initially, when $pn\ll 1$, the dephasing processes do not destroy the
quantum entanglement in Bell pairs, hence the operator entanglement initially
scales as $4n$. For later evolution times, when $pn\gg 1$, only the classical
correlations remain, which is reflected in the operator entanglement scaling
as $2n$. More generally, the operator entanglement is given by:

$\displaystyle E_{\mathrm{MPO}}(p,n)$
$\displaystyle=2n+\frac{2n}{1+\left(1-p\right)^{8n}}\bigg{\\{}\log_{2}\left[1+\left(1-p\right)^{8n}\right]$
$\displaystyle+\left(1-p\right)^{8n}\log_{2}\left[1+\left(1-p\right)^{-8n}\right]\bigg{\\}}.$
(50)

To analyse the entanglement dynamics in trajectories, we consider two types of
unravelings. A naive choice corresponds to a Kraus representation of the
channel (48):

$\displaystyle\left.\begin{array}[]{ll}&K_{0}^{(j)}=\sqrt{1-p}~{}\mathbb{1}_{j},\\\
&K_{1}^{(j)}=\sqrt{p}\sigma_{j}^{z},\end{array}\right.$ (53)

to which we refer as $\mathcal{K}^{(j)}$ for brevity. In this case it is easy
to see that the middle-cut EAEE (5) is linearly growing with $n$ as
$\bar{E}(n)=2n$, since the Kraus operators (53) are both unitary and do not
remove any entanglement from the system. A better choice for unraveling this
dynamics is given by an alternative Kraus representation
$\tilde{\mathcal{K}}^{(j)}$ of the channel (48):

$\displaystyle\left.\begin{array}[]{ll}&\tilde{K}_{0}^{(j)}=|0\rangle_{j}\langle
0|+\left(1-2p\right)|1\rangle_{j}\langle 1|,\\\
&\tilde{K}_{1}^{(j)}=2\sqrt{p\left(1-p\right)}~{}|1\rangle_{j}\langle
1|,\end{array}\right.$ (56)

where one of the Kraus operators is a projector onto the excited state,
$\tilde{K}_{1}^{(j)}\propto|1\rangle_{j}\langle 1|$. Applying such an operator
on one qubit removes its entanglement with other qubits. This means that the
resulting EAEE is reduced (compared to the EAEE in the naive choice case) and
is given by:

$\displaystyle\bar{E}_{\tilde{\mathcal{K}}}(p,n)=2n\left[{q}_{p,n}\log_{2}{q}_{p,n}-\tilde{q}_{p,n}\log_{2}\tilde{q}_{p,n}+\frac{1}{2}\right],$
(57)
$\displaystyle\text{with}~{}q_{p,n}=\frac{1+(1-2p)^{4n}}{2}~{}\text{and}~{}\tilde{q}_{p,n}=\frac{(1-2p)^{4n}}{2}.$

In Fig. 4(b) we plot the evolutions of the entanglement measures for various
values of $p$ and for different simulation methods as well as the EoF Wootters
(1998b). From this plot we can draw several conclusions. First, as discussed
above, the MPO method as well as the naive trajectory unraveling result in the
unbounded growth of the (operator) entanglement, which implies an exponential
growth of the bond dimension $\chi$. Second, the alternative unraveling, as
well as the 2-GEO adaptive approach, generate trajectories that obey area-law
scaling in EAEE. Finally, the 2-GEO unraveling results in the EAEE that is
smallest among the considered unraveling schemes and closest to the global
minimum $E_{f}$.

####  III.3.2 Continuous case

The essential features of the above model can also be seen in its continuous-
time generalization given by a spin-$1/2$ XX model Sachdev (1999):

$\displaystyle H_{\mathrm{XX}}=\Omega\sum_{j=1}^{L-1}$
$\displaystyle\left(\sigma_{j}^{+}\sigma_{j+1}^{-}+\mathrm{h.c.}\right)$ (58)

whose coupling to Markovian environments is described by the jump operators
$c_{j}=\sigma_{j}^{z}/2$ and decoherence rates $\gamma_{j}=\gamma$ and where
$\sigma_{j}^{\pm}=\sigma_{j}^{x}\pm i\sigma_{j}^{y}$. This model can be solved
efficiently numerically after a Jordan-Wigner transformation to a free fermion
model Cao  _et al._ (2019); Mbeng  _et al._ (2020).

Figure 5: Correlation and entanglement dynamics in the XX model (58) for
different numerical methods. The system consists of $L=20$ spin-$1/2$
particles and we use jump operators $c_{j}=\sigma_{j}^{z}/2$ and
$\gamma=0.4\Omega$. (a), (b) and (c): $\sigma^{z}$-correlations of the first
qubit with the fifth, tenth and fifteenth qubits in the chain simulated with
the maximal bond dimension $\chi=256$ by (a) the MPO method, (b) the “quantum
jump” QT MPS (total of $M=2000$ trajectories) and (c) the “quantum state
diffusion” QT MPS (total of $M=15000$ trajectories). The dotted black lines
correspond to the analytical solution. The statistical error bars in
trajectory simulations are denoted by gray filling. Panel (d) compares the
dynamics of the (operator) entanglement of the MPO and QT MPS simulations for
bond dimensions $\chi=256$ (solid lines), $\chi=128$ (dashed lines) and
$\chi=64$ (dotted lines). The EAEE lines of “quantum state diffusion”
trajectories overlap, which indicates their convergence even for smaller
truncation values of the bond dimension.

In Fig. 5 we solve this model with the MPO and QT MPS methods and compare it
to the exact solution888Due to the dispersion processes in this model the
steady-state correlation values are all equal to $1/L$ and vanish in the
thermodynamic limit.. Specifically, we consider the correlation and
entanglement dynamics of $L=20$ qubits initially in the state (47). In panels
(a), (b) and (c) we fix the truncation bond dimension at $\chi=256$ and
compare the $\sigma^{z}$-correlation function obtained by the MPO method and
the QT MPS unravelings with the “quantum jump” Kraus representation
$\mathcal{K}_{dt}$ from Eq. (31) and with the “quantum state diffusion”
representation $\mathcal{R}_{dt}=T(\theta,\varphi)\mathcal{K}_{dt}$, where
$\theta=\pi/4$ and $\varphi=0$. We can see that for both the MPO and “quantum
jump” QT MPS approaches the bond dimension is not sufficient to correctly
capture the correlation dynamics999We note that the MPO simulation
nevertheless correctly captures the steady state values.. In contrast, the
trajectories obtained by the “quantum state diffusion” match the exact results
at this value of the bond dimension (within the statistical accuracy). The
origin of this difference is visible in panel (d), where we show the dynamics
of the middle-cut (operator) entanglement for various bond dimensions. We
observe that the MPO method exhibits an entanglement barrier Reid and Bertini
(2021); Wellnitz  _et al._ (2022), a well-known computational issue occurring
in quantum many-body system simulations. A similar behaviour can be also seen
for the “quantum jump” QT MPS method, for which the average entanglement
rapidly grows during the evolution and saturates at a value set by the bond
dimension. In both cases the lack of convergence in the bond dimension is
apparent. In contrast, the average entanglement of “quantum state diffusion”
trajectories is small and converges with bond dimension. In the SM IV we
consider a more general case of the XX model with fields and give a more
detailed discussion on the convergence of observables with the bond dimension.

## Conclusions

In this work we studied strategies to simulate open quantum many-body system
dynamics by combining quantum trajectory methods with tensor-network
techniques. The key concept underlying our work is the fact that the same
system dynamics can be simulated by infinitely many quantum trajectory
unravelings. These unravelings however can significantly differ in the
ensemble-averaged entanglement entropy (EAEE) that acts as a proxy of the
average computational cost associated with tensor-network representation of
individual trajectories. To illustrate this crucial difference, we considered
several examples, to which we applied various trajectory unravelings as well
as another widely used approach to simulate open quantum many-body dynamics
based on matrix product operators (MPOs). These included models that exhibit
qualitative differences across different unraveling schemes, as well as models
that feature exponential advantage of the trajectory unravelings over MPO
simulations. These considerations motivated the development of adaptive
trajectory unravelings that minimize the EAEE. We proposed several strategies
that utilize this concept by directly minimizing the EAEE (greedy entanglement
optimization, GEO) or by minimizing the instantaneous EAEE change rate (greedy
entanglement derivative optimization, GEDO, proposed in Ref. Vovk and Pichler
(2022)). We performed a comprehensive comparison of the GEO and GEDO
approaches and concluded that the GEDO method allows to efficiently generate
trajectories with low EAEE and that the GEO approach (specifically, 2-GEO)
provides trajectory configurations with a lower EAEE value at a cost of less
efficient numerical performance. In summary, our work paves the way towards
better understanding of the classical simulatability of open quantum many-body
systems and defines new avenues for the development of algorithms for noisy
quantum many-body dynamics.

There are several directions left for future work. First, our methods are
limited to the optimization of individual decoherence channels. It would be
interesting to generalize our approach to the collective optimization. Second,
it would be interesting to compare or combine the developed methods with other
approaches for open system dynamics Nagy and Savona (2018, 2019); Carisch and
Zilberberg (2023); Gravina and Savona (2023); Müller  _et al._ (2024). Lastly,
in the context of measurement-induced phase transitions it would be
interesting to determine whether adaptive unravelings of master equations can
change the value of the critical point, and, e.g., give area-law trajectories
in regimes where non-adaptive unravelings can not.

## Acknowledgements

We thank Lukas Sieberer for carefully reading the manuscript draft and giving
constructive feedback. We thank Peter Zoller, Barbara Kraus, Tom Westerhout,
Matteo Magoni, Luca Arceci, and Gabriele Calliari for fruitful discussions.
The computational results presented here have been achieved using the LinuX
Cluster of the Institute for Theoretical Physics of the University of
Innsbruck and LEO HPC infrastructure of the University of Innsbruck. This work
is supported by the European Union’s Horizon Europe research and innovation
program under Grant Agreement No. 101113690 (PASQuanS2.1), the ERC Starting
grant QARA (Grant No. 101041435), the EU-QUANTERA project TNiSQ (N-6001), and
by the Austrian Science Fund (FWF): COE 1 and quantA. For open access
purposes, the author has applied a CC BY public copyright license to any
author accepted manuscript version arising from this submission.

## References

  * Wootters (1998a) W. K. Wootters, Philos. Trans. Royal Soc. A 356, 1717 (1998a).
  * Amico  _et al._ (2008) L. Amico, R. Fazio, A. Osterloh, and V. Vedral, Rev. Mod. Phys. 80, 517 (2008).
  * Horodecki  _et al._ (2009) R. Horodecki, P. Horodecki, M. Horodecki, and K. Horodecki, Rev. Mod. Phys. 81, 865 (2009).
  * Calabrese  _et al._ (2009) P. Calabrese, J. Cardy, and B. Doyon, J. Phys. A 42, 500301 (2009).
  * Eisert  _et al._ (2010) J. Eisert, M. Cramer, and M. B. Plenio, Rev. Mod. Phys. 82, 277 (2010).
  * Vidal (2003) G. Vidal, Phys. Rev. Lett. 91, 147902 (2003).
  * Eisert and Osborne (2006) J. Eisert and T. J. Osborne, Phys. Rev. Lett. 97, 150404 (2006).
  * Prosen and Znidaric (2007) T. Prosen and M. Znidaric, Phys. Rev. E 75, 015202 (2007).
  * Schuch  _et al._ (2008) N. Schuch, M. M. Wolf, F. Verstraete, and J. I. Cirac, Phys. Rev. Lett. 100, 030504 (2008).
  * Calabrese and Cardy (2007) P. Calabrese and J. Cardy, J. Stat. Mech. Theory Exp. 2007, P06008 (2007).
  * Polkovnikov  _et al._ (2011) A. Polkovnikov, K. Sengupta, A. Silva, and M. Vengalattore, Rev. Mod. Phys. 83, 863 (2011).
  * Eisert  _et al._ (2015) J. Eisert, M. Friesdorf, and C. Gogolin, Nat. Phys. 11, 124 (2015).
  * Nahum  _et al._ (2017) A. Nahum, J. Ruhman, S. Vijay, and J. Haah, Phys. Rev. X 7, 031016 (2017).
  * Aaronson (2005) S. Aaronson, Proc. R. Soc. A: Math. Phys. Eng. Sci. 461, 3473 (2005).
  * Aaronson and Chen (2016) S. Aaronson and L. Chen, arXiv preprint arXiv:1612.05903 (2016), 10.48550/arXiv.1612.05903.
  * Bouland  _et al._ (2019) A. Bouland, B. Fefferman, C. Nirkhe, and U. Vazirani, Nat. Phys. 15, 159 (2019).
  * Napp  _et al._ (2022) J. C. Napp, R. L. La Placa, A. M. Dalzell, F. G. Brandao, and A. W. Harrow, Phys. Rev. X 12, 021021 (2022).
  * Bremner  _et al._ (2016) M. J. Bremner, A. Montanaro, and D. J. Shepherd, Phys. Rev. Lett. 117, 080501 (2016).
  * Harrow and Montanaro (2017) A. W. Harrow and A. Montanaro, Nature 549, 203 (2017).
  * Boixo  _et al._ (2018) S. Boixo, S. V. Isakov, V. N. Smelyanskiy, R. Babbush, N. Ding, Z. Jiang, M. J. Bremner, J. M. Martinis, and H. Neven, Nat. Phys. 14, 595 (2018).
  * Arute  _et al._ (2019) F. Arute, K. Arya, R. Babbush, D. Bacon, J. C. Bardin, R. Barends, R. Biswas, S. Boixo, F. G. Brandao, D. A. Buell, _et al._ , Nature 574, 505 (2019).
  * Daley  _et al._ (2022) A. J. Daley, I. Bloch, C. Kokail, S. Flannigan, N. Pearson, M. Troyer, and P. Zoller, Nature 607, 667 (2022).
  * De las Cuevas  _et al._ (2013) G. De las Cuevas, N. Schuch, D. Pérez-García, and J. I. Cirac, New J. Phys. 15, 123021 (2013).
  * Yung and Gao (2017) M.-H. Yung and X. Gao, arXiv preprint arXiv:1706.08913 (2017), 10.48550/arXiv.1706.08913.
  * Preskill (2018) J. Preskill, Quantum 2, 79 (2018).
  * Cross  _et al._ (2019) A. W. Cross, L. S. Bishop, S. Sheldon, P. D. Nation, and J. M. Gambetta, Phys. Rev. A 100, 032328 (2019).
  * Weimer  _et al._ (2021) H. Weimer, A. Kshetrimayum, and R. Orús, Rev. Mod. Phys. 93, 015008 (2021).
  * Werner  _et al._ (2016) A. H. Werner, D. Jaschke, P. Silvi, M. Kliesch, T. Calarco, J. Eisert, and S. Montangero, Phys. Rev. Lett. 116, 237201 (2016).
  * White  _et al._ (2018) C. D. White, M. Zaletel, R. S. K. Mong, and G. Refael, Phys. Rev. B 97, 035127 (2018).
  * Yoshioka and Hamazaki (2019) N. Yoshioka and R. Hamazaki, Phys. Rev. B 99, 214306 (2019).
  * Vicentini  _et al._ (2019) F. Vicentini, A. Biella, N. Regnault, and C. Ciuti, Phys. Rev. Lett. 122, 250503 (2019).
  * Carisch and Zilberberg (2023) C. Carisch and O. Zilberberg, Quantum 7, 954 (2023).
  * Zoller  _et al._ (1987) P. Zoller, M. Marte, and D. Walls, Phys. Rev. A 35, 198 (1987).
  * Dum  _et al._ (1992) R. Dum, P. Zoller, and H. Ritsch, Phys. Rev. A 45, 4879 (1992).
  * Dalibard  _et al._ (1992) J. Dalibard, Y. Castin, and K. Mølmer, Phys. Rev. Lett. 68, 580 (1992).
  * Gardiner  _et al._ (1992) C. W. Gardiner, A. S. Parkins, and P. Zoller, Phys. Rev. A 46, 4363 (1992).
  * Gisin and Percival (1992) N. Gisin and I. C. Percival, J. Phys. A Math. Gen. 25, 5677 (1992).
  * Tian and Carmichael (1992) L. Tian and H. Carmichael, Phys. Rev. A 46, R6801 (1992).
  * Castin  _et al._ (1993) Y. Castin, J. Dalibard, and K. Mølmer, in  _AIP Conf. Proc._, Vol. 275 (1993) pp. 143–156.
  * Van Dorsselaer and Nienhuis (1998) F. Van Dorsselaer and G. Nienhuis, Eur. Phys. J. D 2, 175 (1998).
  * Plenio and Knight (1998) M. B. Plenio and P. L. Knight, Rev. Mod. Phys. 70, 101 (1998).
  * Schollwöck (2011) U. Schollwöck, Ann. Phys. 326, 96 (2011).
  * Orús (2014) R. Orús, Ann. Phys. 349, 117 (2014).
  * Daley (2014) A. J. Daley, Adv. Phys. 63, 77 (2014).
  * Orús (2019) R. Orús, Nat. Rev. Phys. 1, 538 (2019).
  * Vovk and Pichler (2022) T. Vovk and H. Pichler, Phys. Rev. Lett. 128, 243601 (2022).
  * Li  _et al._ (2018) Y. Li, X. Chen, and M. P. Fisher, Phys. Rev. B 98, 205136 (2018).
  * Skinner  _et al._ (2019) B. Skinner, J. Ruhman, and A. Nahum, Phys. Rev. X 9, 031009 (2019).
  * Chan  _et al._ (2019) A. Chan, R. M. Nandkishore, M. Pretko, and G. Smith, Phys. Rev. B 99, 224307 (2019).
  * Li  _et al._ (2019) Y. Li, X. Chen, and M. P. Fisher, Phys. Rev. B 100, 134306 (2019).
  * Choi  _et al._ (2020) S. Choi, Y. Bao, X.-L. Qi, and E. Altman, Phys. Rev. Lett. 125, 030505 (2020).
  * Zabalo  _et al._ (2020) A. Zabalo, M. J. Gullans, J. H. Wilson, S. Gopalakrishnan, D. A. Huse, and J. Pixley, Phys. Rev. B 101, 060301 (2020).
  * Jian  _et al._ (2020) C.-M. Jian, Y.-Z. You, R. Vasseur, and A. W. Ludwig, Phys. Rev. B 101, 104302 (2020).
  * Turkeshi  _et al._ (2020) X. Turkeshi, R. Fazio, and M. Dalmonte, Phys. Rev. B 102, 014315 (2020).
  * Bao  _et al._ (2020) Y. Bao, S. Choi, and E. Altman, Phys. Rev. B 101, 104301 (2020).
  * Van Regemortel  _et al._ (2022) M. Van Regemortel, O. Shtanko, L. P. García-Pintos, A. Deshpande, H. Dehghani, A. V. Gorshkov, and M. Hafezi, Phys. Rev. Res 4, L032021 (2022).
  * Sharma  _et al._ (2022) S. Sharma, X. Turkeshi, R. Fazio, and M. Dalmonte, SciPost Phys. Core 5, 023 (2022).
  * Chen  _et al._ (2023) Z. Chen, Y. Bao, and S. Choi, arXiv preprint arXiv:2306.17161 (2023), 10.48550/arXiv.2306.17161.
  * Cheng and Ippoliti (2023) Z. Cheng and M. Ippoliti, PRX Quantum 4, 040326 (2023).
  * Hauser  _et al._ (2023) J. Hauser, Y. Li, S. Vijay, and M. Fisher, arXiv preprint arXiv:2304.13198 (2023), 10.48550/arXiv.2304.13198.
  * Holland  _et al._ (1996) M. Holland, S. Marksteiner, P. Marte, and P. Zoller, Phys. Rev. Lett. 76, 3683 (1996).
  * Gardiner and Zoller (2015) C. Gardiner and P. Zoller,  _The Quantum World of Ultra-Cold Atoms and Light Book II: The Physics of Quantum-Optical Devices_ (World Scientific Publishing Company, 2015).
  * Fuji and Ashida (2020) Y. Fuji and Y. Ashida, Phys. Rev. B 102, 054302 (2020).
  * Audenaert  _et al._ (2001) K. Audenaert, F. Verstraete, and B. De Moor, Phys. Rev. A 64, 052304 (2001).
  * Guevara and Viviescas (2014) I. Guevara and C. Viviescas, Phys. Rev. A 90, 012338 (2014).
  * Arceci  _et al._ (2022) L. Arceci, P. Silvi, and S. Montangero, Phys. Rev. Lett. 128, 040501 (2022).
  * Holevo (2022) A. Holevo, Lobachevskii J. Math. 43, 1646 (2022).
  * Holevo (2023) A. Holevo, Lobachevskii J. Math. 44, 2033 (2023).
  * Carisch  _et al._ (2023) C. Carisch, A. Romito, and O. Zilberberg, Phys. Rev. Res. 5, L042031 (2023).
  * Nielsen and Chuang (2010) M. A. Nielsen and I. L. Chuang,  _Quantum computation and quantum information_ (Cambridge university press, 2010).
  * Hughston  _et al._ (1993) L. P. Hughston, R. Jozsa, and W. K. Wootters, Phys. Lett. A 183, 14 (1993).
  * Bennett  _et al._ (1996) C. H. Bennett, D. P. DiVincenzo, J. A. Smolin, and W. K. Wootters, Phys. Rev. A 54, 3824 (1996).
  * Rockafellar (1997) R. T. Rockafellar,  _Convex analysis_ , Vol. 11 (Princeton university press, 1997).
  * Uhlmann (1998) A. Uhlmann, Open Syst. Inf. Dyn. 5, 209 (1998).
  * Uhlmann (2010) A. Uhlmann, Entropy 12, 1799 (2010).
  * Gurvits (2003) L. Gurvits, in  _Proceedings of the thirty-fifth annual ACM symposium on Theory of computing_ (2003) pp. 10–19.
  * Gharibian (2008) S. Gharibian, arXiv preprint arXiv:0810.4507 (2008), 10.48550/arXiv.0810.4507.
  * Lindblad (1976) G. Lindblad, Comm. Math. Phys. 48, 119 (1976).
  * Gorini  _et al._ (1976) V. Gorini, A. Kossakowski, and E. C. G. Sudarshan, J. Math. Phys. 17, 821 (1976).
  * Trotter (1959) H. F. Trotter, Proc. Am. Math. Soc. 10, 545 (1959).
  * Hedemann (2013) S. R. Hedemann, arXiv preprint arXiv:1303.5904 (2013), 10.48550/arXiv.1303.5904.
  * Wootters (1998b) W. K. Wootters, Phys. Rev. Lett. 80, 2245 (1998b).
  * Hayden  _et al._ (2006) P. Hayden, D. W. Leung, and A. Winter, Commun. Math. Phys. 265, 95 (2006).
  * Bonnes and Läuchli (2014) L. Bonnes and A. M. Läuchli, arXiv:1411.4831 (2014).
  * Wellnitz  _et al._ (2022) D. Wellnitz, G. Preisser, V. Alba, J. Dubail, and J. Schachenmayer, Phys. Rev. Lett. 129, 170401 (2022).
  * Preisser  _et al._ (2023) G. Preisser, D. Wellnitz, T. Botzung, and J. Schachenmayer, Phys. Rev. A 108, 012616 (2023).
  * Sachdev (1999) S. Sachdev, Phys. World 12, 33 (1999).
  * Cao  _et al._ (2019) X. Cao, A. Tilloy, and A. De Luca, SciPost Phys. 7, 024 (2019).
  * Mbeng  _et al._ (2020) G. B. Mbeng, A. Russomanno, and G. E. Santoro, arXiv preprint arXiv:2009.09208 114 (2020), 10.48550/arXiv.2009.09208.
  * Reid and Bertini (2021) I. Reid and B. Bertini, Phys. Rev. B 104, 014301 (2021).
  * Nagy and Savona (2018) A. Nagy and V. Savona, Phys. Rev. A 97, 052129 (2018).
  * Nagy and Savona (2019) A. Nagy and V. Savona, Phys. Rev. Lett. 122, 250501 (2019).
  * Gravina and Savona (2023) L. Gravina and V. Savona, arXiv preprint arXiv:2312.13676 (2023), 10.48550/arXiv.2312.13676.
  * Müller  _et al._ (2024) A. Müller, T. Ayral, and C. Bertrand, arXiv preprint arXiv:2403.00152 (2024), 10.48550/arXiv.2403.00152.
  * Broyden (1970) C. G. Broyden, IMA J. Appl. Math. 6, 76 (1970).
  * Fletcher (1970) R. Fletcher, Comput. J. 13, 317 (1970).
  * Goldfarb (1970) D. Goldfarb, Math. Comput. 24, 23 (1970).
  * Shanno (1970) D. F. Shanno, Math. Comput. 24, 647 (1970).
  * van Nieuwenburg  _et al._ (2017) E. P. van Nieuwenburg, J. Y. Malo, A. J. Daley, and M. H. Fischer, Quantum Sci. Technol. 3, 01LT02 (2017).

Supplementary Material

##  I Optimization over $T\in\mathcal{T}_{2}^{2}$

In this section we provide the details of the numerical solution of the
optimization problem (32). Here we consider a limited class of isometries
$T(\theta,\varphi)\in\mathcal{T}_{2}^{2}\subset\mathcal{T}_{2}^{4}$ explicitly
given in (34). Applying such an isometry to the reference ensemble
$\tilde{\Psi}_{dt}=\mathcal{K}_{dt}|\psi\rangle=\left(K_{0}|\psi\rangle,K_{1}|\psi\rangle\right)$
results in a new ensemble of the form111For simplicity in what follows we
assume the Kraus operators and corresponding jump operators to be local.:

$\displaystyle\begin{pmatrix}|\tilde{\phi}_{0}\left(\theta,\varphi,dt\right)\rangle\\\
|\tilde{\phi}_{1}\left(\theta,\varphi,dt\right)\rangle\\\
\end{pmatrix}=T(\theta,\varphi)\begin{pmatrix}K_{0}|\psi\rangle\\\
K_{1}|\psi\rangle\\\ \end{pmatrix}.$ (S1)

Inserting it into the EAEE (5) allows to reformulate the optimization problem
(32) as

$\displaystyle\inf_{T\in\mathcal{T}_{2}^{4}}\bar{E}[T\tilde{\Psi}_{dt}]\xrightarrow{T=T(\theta,\varphi)}\min_{\theta,\varphi}\bar{E}\left(\theta,\varphi,dt\right),$
(S2)

the solution of which we aim to explain below. The EAEE
$\bar{E}\left(\theta,\varphi,dt\right)$ in Eq. (S2) is convex, which allows to
define unambiguously the minimum using the standard optimization methods. We
choose the Broyden-Fletcher-Goldfarb-Shanno (BFGS) algorithm – a quasi-Newton
gradient-based method, which efficiently solves nonlinear optimization
problems Broyden (1970); Fletcher (1970); Goldfarb (1970); Shanno (1970).
Below we provide an explicit form of the EAEE as a function of the
optimization parameters and contraction schemes of MPS networks needed for
BFGS optimization.

Let us write the EAEE as a function of the isometry parameters222From here on
we do not write $dt$ for notational clarity.:

$\displaystyle\bar{E}\left(\theta,\varphi\right)=\sum_{\alpha=0,1}p_{\alpha}\left(\theta,\varphi\right)E\left(|{\phi}_{\alpha}\left(\theta,\varphi\right)\rangle\right),$
(S3)

where
$p_{\alpha}\left(\theta,\varphi\right)=\langle\tilde{\phi}_{\alpha}\left(\theta,\varphi\right)|\tilde{\phi}_{\alpha}\left(\theta,\varphi\right)\rangle$
are the probability weights of the ensemble (S1) and
$E\left(|{\phi}_{\alpha}\left(\theta,\varphi\right)\rangle\right)$ is the von
Neumann entanglement entropy (1) of the normalized state
$|{\phi}_{\alpha}\left(\theta,\varphi\right)\rangle=|\tilde{\phi}_{\alpha}\left(\theta,\varphi\right)\rangle/\sqrt{p_{\alpha}\left(\theta,\varphi\right)}$.
If we define the partial trace
$\tilde{\phi}_{\alpha}\left(\theta,\varphi\right)=\mathrm{tr}_{B}\left(|\tilde{\phi}_{\alpha}\left(\theta,\varphi\right)\rangle\langle\tilde{\phi}_{\alpha}\left(\theta,\varphi\right)|\right)$,
we can write:

$\displaystyle\bar{E}\left(\theta,\varphi\right)=\sum_{\alpha=0,1}p_{\alpha}\left(\theta,\varphi\right)\log_{2}p_{\alpha}\left(\theta,\varphi\right)-\mathrm{tr}\left[\tilde{\phi}_{\alpha}\left(\theta,\varphi\right)\log_{2}\left(\tilde{\phi}_{\alpha}\left(\theta,\varphi\right)\right)\right].$
(S4)

The gradients of such a function have a simple form:

$\displaystyle\partial_{\theta,\varphi}\bar{E}\left(\theta,\varphi\right)=-\sum_{\alpha=0,1}\mathrm{tr}\left[\partial_{\theta,\varphi}\tilde{\phi}_{\alpha}\left(\theta,\varphi\right)\log_{2}\left({\phi}_{\alpha}\left(\theta,\varphi\right)\right)\right],$
(S5)

where
${\phi}_{\alpha}\left(\theta,\varphi\right)=\tilde{\phi}_{\alpha}\left(\theta,\varphi\right)/p_{\alpha}\left(\theta,\varphi\right)$.
One can find the explicit forms of the probability weights
$p_{\alpha}\left(\theta,\varphi\right)$ and the partial traces
$\tilde{\phi}_{\alpha}\left(\theta,\varphi\right)$ in terms of the
optimization parameters $\theta,\varphi$ and the initial ensemble
$\tilde{\Psi}_{dt}=(K_{0}|\psi\rangle,K_{1}|\psi\rangle)=(|\tilde{\psi}_{0}\rangle,|\tilde{\psi}_{1}\rangle)$:

$\displaystyle
p_{0}\left(\theta,\varphi\right)=\cos^{2}(\theta)q_{0}+\sin^{2}(\theta)q_{1}+2\cos(\theta)\sin(\theta)\cos(\varphi-\upsilon_{01})q_{01},$
(S6) $\displaystyle
p_{1}\left(\theta,\varphi\right)=\sin^{2}(\theta)q_{0}+\cos^{2}(\theta)q_{1}-2\cos(\theta)\sin(\theta)\cos(\varphi-\upsilon_{01})q_{01},$
(S7)
$\displaystyle\tilde{\phi}_{0}\left(\theta,\varphi\right)=\cos^{2}(\theta)\tilde{\psi}_{0}+\sin^{2}(\theta)\tilde{\psi}_{1}+\cos(\theta)\sin(\theta)\left(e^{-i\varphi}\tilde{\psi}_{10}+\mathrm{h.c.}\right),$
(S8)
$\displaystyle\tilde{\phi}_{1}\left(\theta,\varphi\right)=\sin^{2}(\theta)\tilde{\psi}_{0}+\cos^{2}(\theta)\tilde{\psi}_{1}-\cos(\theta)\sin(\theta)\left(e^{-i\varphi}\tilde{\psi}_{10}+\mathrm{h.c.}\right).$
(S9)

where $q_{\alpha}=\langle\tilde{\psi}_{\alpha}|\tilde{\psi}_{\alpha}\rangle$
are the probability weights of the initial ensemble,
$q_{01}=|\langle\tilde{\psi}_{0}|\tilde{\psi}_{1}\rangle|$ and
$\upsilon_{01}=\arg\left(\langle\tilde{\psi}_{0}|\tilde{\psi}_{1}\rangle\right)$
are amplitude and phase of the overlap of the initial ensemble states,
$\tilde{\psi}_{\alpha}=\mathrm{tr}_{B}\left(|\tilde{\psi}_{\alpha}\rangle\langle\tilde{\psi}_{\alpha}|\right)$
are the partial traces of the initial ensemble and
$\tilde{\psi}_{10}=\mathrm{tr}_{B}\left(|\tilde{\psi}_{1}\rangle\langle\tilde{\psi}_{0}|\right)$.
With this one can calculate the derivatives of the matrices:

$\displaystyle\partial_{\theta}\tilde{\phi}_{0}\left(\theta,\varphi\right)=2\cos(\theta)\sin(\theta)\left(\tilde{\psi}_{1}-\tilde{\psi}_{0}\right)+2\left[\cos^{2}(\theta)-\sin^{2}(\theta)\right]\left(e^{-i\varphi}\tilde{\psi}_{10}+\mathrm{h.c.}\right),$
(S10)
$\displaystyle\partial_{\theta}\tilde{\phi}_{1}\left(\theta,\varphi\right)=2\cos(\theta)\sin(\theta)\left(\tilde{\psi}_{0}-\tilde{\psi}_{1}\right)-2\left[\cos^{2}(\theta)-\sin^{2}(\theta)\right]\left(e^{-i\varphi}\tilde{\psi}_{10}+\mathrm{h.c.}\right),$
(S11)
$\displaystyle\partial_{\varphi}\tilde{\phi}_{0}\left(\theta,\varphi\right)=i\cos(\theta)\sin(\theta)\left(-e^{-i\varphi}\tilde{\psi}_{10}+e^{i\varphi}\tilde{\psi}_{01}\right),$
(S12)
$\displaystyle\partial_{\varphi}\tilde{\phi}_{1}\left(\theta,\varphi\right)=i\cos(\theta)\sin(\theta)\left(e^{-i\varphi}\tilde{\psi}_{10}-e^{i\varphi}\tilde{\psi}_{01}\right),$
(S13)

which can be inserted into (S5) to calculate the gradient.

In the case when the system state is represented in the MPS form, one can
calculate the EAEE (S4) and its gradients (S5) by contracting tensor networks
as shown in Fig. S1. There are two important aspects to notice. First, this
optimization is numerically costly when using MPS contractions. The reason for
this is that for every set of parameters $\theta,\varphi$ one has to calculate
singular values of the MPS for every new optimization iteration, since these
values depend implicitly on the optimizing parameters. This can be seen by
analyzing the second term of Eq. (S4). This term can be written using the
singular value decomposition as:

$\displaystyle\mathrm{tr}\left[\tilde{\phi}_{\alpha}\left(\theta,\varphi\right)\log_{2}\left(\tilde{\phi}_{\alpha}\left(\theta,\varphi\right)\right)\right]=\sum_{i}\tilde{\lambda}^{(i)}_{\alpha}\left(\theta,\varphi\right)\log_{2}\tilde{\lambda}^{(i)}_{\alpha}\left(\theta,\varphi\right),$
(S14)

where $\tilde{\lambda}^{(i)}_{\alpha}\left(\theta,\varphi\right)$ are the
eigenvalues of $\tilde{\phi}_{\alpha}\left(\theta,\varphi\right)$. Second, in
the case of the gradient calculation (S5), in addition to singular value
calculation one has to explicitly contract the corresponding tensor networks
instead of just using the singular values [see Fig. S1(g)].

Figure S1: The MPS tensor objects in Vidal’s $\Gamma\lambda$-form Vidal (2003)
and their contractions required to calculate the EAEE and it gradients. (a)
The initial state $|\psi\rangle$ with the bipartition $AB$. (b) The initial
ensemble state $|\tilde{\psi}_{\alpha}\rangle$ obtained by applying the
corresponding Kraus operator $K_{\alpha}$. (c) The transformed ensemble state
$|\tilde{\phi}_{\alpha}\left(\theta,\varphi\right)\rangle$ with
$R_{\alpha}(\theta,\varphi)=\sum_{\beta}T_{\alpha\beta}(\theta,\varphi)K_{\beta}$.
(d) The probability weight $p_{\alpha}\left(\theta,\varphi\right)$. (e) The
partial trace
$\tilde{\phi}_{\alpha}\left(\theta,\varphi\right)=\mathrm{tr}_{B}\left(|\tilde{\phi}^{(\alpha)}\left(\theta,\varphi\right)\rangle\langle\tilde{\phi}^{(\alpha)}\left(\theta,\varphi\right)|\right)$.
(f) The term
$\mathrm{tr}\left[\tilde{\phi}_{\alpha}\left(\theta,\varphi\right)\log_{2}\left(\tilde{\phi}_{\alpha}\left(\theta,\varphi\right)\right)\right]$
from (S14). (g) The term
$\mathrm{tr}\left[\tilde{\psi}_{0}\log_{2}\left({\phi}_{\alpha}\left(\theta,\varphi\right)\right)\right]$
as an example of a building block to calculate the EAEE gradient (S5).

##  II Comment on optimization over $T\in\mathcal{T}_{2}^{4}$

So far we covered the EAEE optimization over a limited class of isometries
$T\in\mathcal{T}_{2}^{2}$. However, in order to guaranteedly attain the
minimum (32), one has to consider the full class of isometries
$\mathcal{T}_{2}^{4}$ Hedemann (2013). This implies optimizations of $3\times
2$ right-unitary matrices with five parameters:

$\displaystyle
T(\theta_{1},\theta_{2},\theta_{3},\alpha,\beta)=\begin{pmatrix}\cos(\theta_{1})&&\sin(\theta_{1})\cos(\theta_{2})e^{i\alpha}\\\
\sin(\theta_{1})\cos(\theta_{3})&\begin{split}&-\cos(\theta_{1})\cos(\theta_{2})\cos(\theta_{3})e^{i\alpha}\\\
&-\sin(\theta_{2})\sin(\theta_{3})e^{-i\beta}\end{split}\\\
\sin(\theta_{1})\sin(\theta_{3})&\begin{split}&-\cos(\theta_{1})\cos(\theta_{2})\sin(\theta_{3})e^{i\alpha}\\\
&+\sin(\theta_{2})\cos(\theta_{3})e^{-i\beta}\end{split}\\\ \end{pmatrix},$
(S15)

and of $4\times 2$ right-unitary matrices with eight parameters:

$\displaystyle
T(\theta_{1},\theta_{2},\theta_{3},\theta_{4},\theta_{5},\alpha,\beta,\gamma)=$
$\displaystyle~{}~{}~{}~{}\begin{pmatrix}\cos(\theta_{1})&&\sin(\theta_{1})\cos(\theta_{2})e^{i\alpha}\\\
\sin(\theta_{1})\cos(\theta_{3})&&\sin(\theta_{2})\sin(\theta_{3})\cos(\theta_{5})e^{i\beta}-\cos(\theta_{1})\cos(\theta_{2})\cos(\theta_{3})e^{i\alpha}\\\
\sin(\theta_{1})\sin(\theta_{3})\cos(\theta_{4})&\begin{split}&\sin(\theta_{2})\sin(\theta_{4})\sin(\theta_{5})e^{i\gamma}-\cos(\theta_{1})\cos(\theta_{2})\sin(\theta_{3})\cos(\theta_{4})e^{i\alpha}\\\
&~{}~{}~{}~{}~{}~{}-\sin(\theta_{2})\cos(\theta_{3})\cos(\theta_{4})\cos(\theta_{5})e^{i\beta}\end{split}\\\
\sin(\theta_{1})\sin(\theta_{3})\sin(\theta_{4})&\begin{split}&~{}~{}~{}-\cos(\theta_{1})\cos(\theta_{2})\sin(\theta_{3})\sin(\theta_{4})e^{i\alpha}-\sin(\theta_{2})\cos(\theta_{3})\sin(\theta_{4})\cos(\theta_{5})e^{i\beta}\\\
&~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}-\sin(\theta_{2})\cos(\theta_{4})\sin(\theta_{5})e^{i\gamma}\end{split}\\\
\end{pmatrix}.$ (S16)

The EAEE and its gradients associated with the isometries (S15) and (S16) have
a more complicated form than the ones shown in Section I, but the calculations
are similar. The important aspect to notice in this case is that the high-
dimensional optimization landscapes for (S15) and (S16) make the optimization
process more complicated.

##  III Numerical analysis of single-step optimization: A case of ten qubits

Figure S2: Middle-cut EAEE change rates of the system of $L=10$ qubits
obtained by the 2-GEO vs. GEDO methods with a single jump operator
$c_{j}=\sigma_{j}^{z}/2$, $\gamma_{j}=1$ and $j=L/2$. Each dot represents a
random state (total of 1250 instances).

In Subsubsection III.1.2 of the main text we compare the performance of the
adaptive unraveling strategies, namely, the 2-GEO and GEDO methods (see
Subsection II.4 of the main text). Here we extend our analysis to larger
systems. For this we consider a system of $L=10$ qubits and perform the same
procedure as in Fig. 1(b). Namely, we take typical states of the 10-qubit
system, let them evolve under a decoherence channel for a single time step and
compare the resulting EAEE change rates for different adaptive unravelings.
Note that, in contrast to the two-qubit scenario discussed in the main text,
the 10-qubit system can be partitioned in multiple ways, and, hence, many
bipartite entanglement values may be optimized. Here we focus on the case of
the middle-cut bipartition and corresponding EAEE values.

As typical states we take a class of random states with various middle-cut
entanglement entropies. These states are generated by propagating a product
state of excited qubits $\bigotimes_{i=1}^{L}|1\rangle_{i}$ with a time-
dependent random Hamiltonian defined in Eq. (40). This fast-scrambling
Hamiltonian increases the entanglement in the qubit system linearly in time.
Therefore, one can explore both weakly and strongly entangled many-body states
by changing the propagation times. In what follows we take 25 different
propagation times starting from $t_{\mathrm{min}}=0.2/\gamma_{j}$ and up to
$t_{\mathrm{max}}=5/\gamma_{j}$, where $\gamma_{j}=1$. The propagation time
step is $dt=0.01/\gamma_{j}$. For each of the propagation durations we apply
50 different random Hamiltonian realizations, which in total gives us 1250
random many-body states with different entanglement entropies.

In Fig. S2 we compare the EAEE values of the 2-GEO and GEDO unravelings of the
10-qubit evolution under the decoherence channel with
$c_{j}=\sigma_{j}^{z}/2$, where $j=5$. From this plot we draw a conclusion
similar to the two-qubit case: the 2-GEO method finds configurations that give
lower EAEE values compared to the GEDO method. We note that, in contrast to
the two-qubit case, in the 10-qubit case the entanglement change rate is no
longer monotonically related to the value of initial-state entanglement
entropy, which can be seen in the distribution of colors in Fig. S2.

##  IV XX model with longitudinal and transverse fields

Complementing the study in Subsubsection III.3.2 of the main text, here we
consider the noisy XX model with fields:

$\displaystyle H_{\mathrm{XXF}}=\Omega\sum_{j=1}^{L-1}$
$\displaystyle\left(\sigma_{j}^{+}\sigma_{j+1}^{-}+\mathrm{h.c.}\right)+\Delta\sum_{j=1}^{L}\sigma_{j}^{z}+h\sum_{j=1}^{L}\sigma_{j}^{x},$
(S17)

which cannot be mapped to free fermions and therefore is not efficiently
numerically solvable. In Fig. S3 we solve this model with the MPO and QT MPS
methods and compare it to the exact diagonalization solution. Specifically, we
consider the correlation and entanglement dynamics in a system of $L=8$ qubits
initially in the state (47). In panels (a), (b) and (c) we plot the
$\sigma^{z}$-correlations obtained by the MPO method and the QT MPS “quantum
jump” and “quantum state diffusion” unravelings. In each panel we vary the
value of the maximal bond dimension $\chi$ and check the convergence of the
simulations towards the exact diagonalization solution. We set the maximal MPO
bond dimension to be the square of the analogous value for the MPSs for
comparison clarity van Nieuwenburg  _et al._ (2017).

We observe that for a fixed maximal bond dimension $\chi=16$ the QT MPS
simulations give exact results, whilst the MPO method significantly deviates
from the exact solution even for a larger maximal bond dimension $\chi=64$. In
turn, comparing the two trajectory approaches, we see that “quantum jump”
trajectories also significantly deviate from the correct results once the
maximal bond dimension is reduced. In contrast to that, the “quantum state
diffusion” trajectories are closer to the exact results for lower $\chi$
values. The origin of these differences is visible in panel (d), where we plot
the dynamics of the middle-cut (operator) entanglement for various bond
dimensions. We observe that for a fixed bond dimension $\chi=16$ the MPO
method is not exact in contrast to the trajectory methods. In addition,
reducing the maximal bond dimension to $\chi=8$ makes “quantum jump”
trajectories deviate from the exact line, which does not take place for
“quantum state diffusion” trajectories.

Figure S3: Correlation and entanglement dynamics in the noisy XX model with
fields (S17) for different numerical methods. The system consists of $L=8$
spin-$1/2$ particles and we use jump operators $c_{j}=\sigma_{j}^{z}/2$ and
$\gamma=0.4\Omega$, $\Delta=h=0.02\Omega$. (a), (b) and (c):
$\sigma^{z}$-correlations of the first qubit with the second qubit (solid
lines) and with the seventh qubit (dashed lines) simulated with various bond
dimensions $\chi$ by (a) the MPO method, (b) the “quantum jump” QT MPS (total
of $M=10000$ trajectories) and (c) the “quantum state diffusion” QT MPS (total
of $M=20000$ trajectories). The dotted black lines correspond to the exact
solution. In this case the error bars are not exceeding the lines’ width and
hence are not shown. Panel (d) compares the dynamics of the (operator)
entanglement of the MPO method, shown in bright orange color with $\chi=256$
(solid line), $64$ (dashed line), $16$ (dash-dotted line) and $4$ (dotted
line), and the QT MPS methods, shown in pink and dark purple colors with
$\chi=16$ (solid lines), $8$ (dashed lines), $4$ (dash-dotted lines) and $2$
(dotted lines).