# Entanglement-Optimal Trajectories of Many-Body Quantum Markov Processes

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

###### Abstract

We develop a novel approach aimed at solving the equations of motion of open
quantum many-body systems. It is based on a combination of generalized wave
function trajectories and matrix product states. We introduce an adaptive
quantum stochastic propagator, which minimizes the expected entanglement in
the many-body quantum state, thus minimizing the computational cost of the
matrix product state representation of each trajectory. We illustrate this
approach on the example of a one-dimensional open Brownian circuit. We show
that this model displays an entanglement phase transition between area and
volume law when changing between different propagators and that our method
autonomously finds an efficiently representable area law unravelling.

Classical simulation of the evolution of quantum many-body systems is a
formidably hard task, in particular if the system is fully coherent Preskill
(2018). Most near-term intermediate scale quantum devices are however noisy,
which opens a possibility for the existence of efficient classical algorithms
for simulating the corresponding open system dynamics. Nonetheless it is often
unclear how to best exploit this potential.

Here we address this challenge by developing an algorithm that explicitly
harnesses the quantum noise inherent to an open quantum system to minimize the
computational cost of representing the many-body state. Our approach is based
on a combination of matrix product states (MPSs) White (1992); Vidal (2003);
Schollwöck (2011) and a generalization of the quantum trajectory (QT) method
Dum  _et al._ (1992); Dalibard  _et al._ (1992); Gardiner  _et al._ (1992);
Tian and Carmichael (1992); Castin  _et al._ (1993); Van Dorsselaer and
Nienhuis (1998). The latter identifies the dynamics of an open quantum system
as a stochastic evolution of pure quantum states, corresponding to a
continuous measurement of the environment Gardiner and Zoller (2015).
Importantly, the choice of the monitored environment observables results in
qualitatively different ensembles of QTs. Our method utilizes this flexibility
and continuously optimizes the monitored environment observables by predicting
and minimizing the expected entanglement in the trajectory wave function (Fig.
1), thus minimizing the computational cost of MPS representations.

We illustrate our approach by applying it to solve the Markovian master
equation (ME) of an open Brownian circuit, where the coherent part of the
evolution rapidly generates entanglement, while the dissipative part leads to
dephasing. We show that various types of QT methods lead to ensembles that
differ dramatically in their entanglement properties. This includes a phase
transition between area law and volume law entangled ensembles, depending on
the monitored environment observables. In addition to being an interesting
phenomenon per se, with connections to recently discussed measurement-induced
phase transitions Li  _et al._ (2018, 2019); Skinner  _et al._ (2019); Chan
_et al._ (2019); Gullans and Huse (2020); Choi  _et al._ (2020); Bao  _et al._
(2020); Zabalo  _et al._ (2020); Ippoliti  _et al._ (2021); Zabalo  _et al._
(2021); Botzung  _et al._ (2021); Müller  _et al._ (2021), it provides an
ideal test bed for our algorithm: We show that our entanglement predictor
allows to generate ensembles of QTs that keep the system in the area law phase
at all times.

Figure 1: (a) Schematic setting of a continuously monitored quantum many-body
system with adaptive measurements. (b) Propagation of a state from time $t$
into the future with two options for the stochastic propagator (color coded,
grey vs. green). The colored dashed lines are potential QTs for a given
unravelling, their stochastic averages are shown as solid lines. A QT (black)
is generated randomly at each time step, corresponding to the stochastic
propagator that gives the lowest expected entanglement (green).

We are interested in open quantum many-body systems with $n$ constituents and
$m$ dissipative channels described by a ME of a generic Lindblad form
($\hbar~{}=~{}1$) Lindblad (1976):

$\frac{d\rho}{dt}=-i\left[H_{\mathrm{sys}},\rho\right]+\sum_{j=1}^{m}\gamma_{j}\left(c_{j}\rho
c_{j}^{\dagger}-\frac{1}{2}\left\\{\rho,c_{j}^{\dagger}c_{j}\right\\}\right).$
(1)

Here, $\rho$ is the (many-body) density operator of the quantum system,
$H_{\rm sys}$ is the (so far unspecified, potentially time-dependent) system
Hamiltonian, and the $c_{j}$’s are the jump operators corresponding to the
$j$-th decay channel with associated decay rate $\gamma_{j}$. In the following
it is convenient to assume that $H_{\rm sys}$ is short-range interacting and
the jump operators are local. Such a ME describes the reduced dynamics of a
system coupled to an environment (or bath). The underlying joint evolution of
system and environment is unitary and is described by the Schrödinger equation
$i\partial_{t}|\Psi(t)\rangle=H(t)|\Psi(t)\rangle$, where $|\Psi(t)\rangle$
represents the joint state of the system and environment, and the Hamiltonian
is given by Gardiner and Zoller (2015)

$\displaystyle H(t)=H_{\rm
sys}+i\sum_{j=1}^{m}\sqrt{\gamma_{j}}\left[b_{j}(t)^{\dagger}c_{j}-c_{j}^{\dagger}b_{j}(t)\right].$
(2)

Here, $b_{j}(t)$’s are so-called quantum noise operators with bosonic
statistics
$[b_{j}(t),b_{j^{\prime}}(t^{\prime})^{\dagger}]=\delta_{j,j^{\prime}}\delta(t-t^{\prime})$.
Those operators act on the degrees of freedom of the environment, which are
initially in the vacuum state: $b_{j}(t)|\Psi(0)\rangle=0$. The ME (1) is
obtained from the full Schrödinger equation by tracing out these bosonic
degrees of freedom Gardiner and Zoller (2015).

The QT method solves the ME (1) by simulating a continuous monitoring of the
environment degrees of freedom in (2). It stochastically generates an ensemble
of pure states of the system, called quantum trajectories, from which the
density matrix dynamics can be retrieved as a statistical average:

$\displaystyle\rho(t)=\lim_{N\rightarrow\infty}\frac{1}{N}\sum_{k=1}^{N}|\phi^{(k)}(t)\rangle\langle\phi^{(k)}(t)|.$
(3)

Here the quantum trajectories, $|\phi^{(k)}(t)\rangle$, represent the
conditional state of the system for a particular (simulated) measurement
history of the monitored environment (enumerated by $k$). They satisfy a
stochastic equation of motion, i.e. a quantum stochastic Schrödinger equation
(QSSE) Gardiner and Zoller (2015). The form of the QSSE depends on the
environment monitoring schemes, in particular on the monitored observables.
For instance, the standard QT method, known as quantum jump approach, is
obtained by continuous monitoring of the photon numbers,
$b_{j}^{\dagger}(t)b_{j}(t)$ Dum  _et al._ (1992); Dalibard  _et al._ (1992);
Gardiner  _et al._ (1992); Castin  _et al._ (1993). Another approach is based
on the balanced homodyne measurement, where a phase-sensitive homodyne current
$b_{j}(t)e^{i\varphi_{j}}+b_{j}^{\dagger}(t)e^{-i\varphi_{j}}$ is continuously
measured Tian and Carmichael (1992). Here, the $\varphi_{j}$ specifies the
measured quadrature.

A QT is generated by integrating the QSSE in (small) time steps, $dt$. In
practice, one Trotterizes the corresponding propagator into a stochastic and a
deterministic part. The stochastic component can be further split into
independent components corresponding to $m$ decay channels Daley (2014).
Specifically, a single integration step
$|\phi^{(k)}(t)\rangle\rightarrow|\phi^{(k)}(t+dt)\rangle$ is achieved by a
sequential application of stochastic operators $K_{j}$ for each decay channel
$j$: Starting from the state
$|\phi^{(k)}_{1}(t)\rangle=|\phi^{(k)}(t)\rangle$, one generates the sequence
of normalized states

$\displaystyle|\phi^{(k)}_{j+1}(t)\rangle$ $\displaystyle\propto
K_{j}|\phi^{(k)}_{j}(t)\rangle,\quad(j=1,\dots,m)$ (4)

and then concludes the integration step with the unitary operation
$|\phi^{(k)}(t+dt)\rangle=e^{-iH_{\mathrm{sys}}dt}|\phi^{(k)}_{m+1}(t)\rangle.$
The form of the stochastic operator $K_{j}$ depends on the type of (simulated)
environment measurement. For example, the operators simulating a homodyne
measurement of the output field in channel $j$ is

$\displaystyle K_{j}^{\mathrm{hom}}=$ $\displaystyle
e^{-\gamma_{j}dtc_{j}^{\dagger}c_{j}/2}+\sqrt{\gamma_{j}}c_{j}e^{i\varphi_{j}}d\xi_{j}(t),$
(5)

where
$d\xi_{j}(t)=\sqrt{\gamma_{j}}\langle\phi^{(k)}_{j}(t)|c_{j}e^{i\varphi_{j}}\\!+c_{j}^{\dagger}e^{-i\varphi_{j}}|\phi^{(k)}_{j}(t)\rangle
dt+dW_{j}(t)$ and the $dW_{j}(t)$’s are independent, normally distributed
Gaussian variable with zero mean and variance $dt$, also known as Wiener
increment. For a simulated number measurement, the propagator is
probabilistically selected from the two options:

$\displaystyle
K_{j}^{\mathrm{num}}=\left[\begin{array}[]{ll}\sqrt{\gamma_{j}dt}c_{j}&\mathrm{~{}with~{}probability~{}}p_{j},\\\
e^{-\gamma_{j}dtc_{j}^{\dagger}c_{j}/2}&\mathrm{~{}with~{}probability~{}}1-p_{j},\end{array}\right.$
(8)

where
$p_{j}=\gamma_{j}dt\langle\phi^{(k)}_{j}(t)|c_{j}^{\dagger}c_{j}|\phi^{(k)}_{j}(t)\rangle$
corresponds to the probability of measuring a photon in channel $j$ between
$t$ and $t+dt$.

While these two schemes are by far the most common ones, we stress that a
multitude of stochastic propagators can be constructed by choosing different
measurements of the bath operators specified via the eigenbases of arbitrary
Hermitian functions of the bath operators,
$X_{j}(t)=f[b_{j}(t),b_{j}^{\dagger}(t)]$. This also includes measurement
strategies that change with time and depend on prior measurement outcomes
111Since the trajectories are by construction sampled according to Born’s
rule, all such strategies, including the adaptive time- and history-dependent
ones, generate trajectories that provide valid solutions of the ME [see
Section I of the Supplementary Materials (SM) for more detail].. Importantly,
the stochastic averages over linear functionals of the state projectors
$\phi^{(k)}=|\phi^{(k)}\rangle\langle\phi^{(k)}|$ do not depend on how the
environment is monitored. This includes in particular the density operator (3)
and the expectation values of linear observables. However the ensemble of QTs
itself, as well as non-linear functionals of the trajectories, do depend on
the unravelling. In particular this holds for the entanglement entropy (EE),
which is of central interest in our work. We consider the EE between two
partitions of the many-body system, $A\cup B$. In this case the EE of a single
QT $|\phi^{(k)}\rangle$ is defined as the von Neumann entropy of the reduced
state of one of the subsystems
{$\rho^{(k)}_{A}=\mathrm{tr}_{B}\left[{\phi^{(k)}}\right]$}:

$\displaystyle
E\left[\phi^{(k)}\right]=S\left[\rho^{(k)}_{A}\right]=-\mathrm{tr}\left[{\rho^{(k)}_{A}\log_{2}\rho^{(k)}_{A}}\right].$
(9)

The ensemble averaged EE (EAEE) is then defined as

$\displaystyle\overline{E}=\lim_{N\rightarrow\infty}\frac{1}{N}\sum_{k=1}^{N}E\left[\phi^{(k)}\right].$
(10)

It is bounded from below by the entanglement of formation Wootters (1998) of
the density operator, $E_{f}(\rho)$, and from above by the entropies of
$\rho_{A}=\mathrm{tr}_{B}\left(\rho\right)$ and
$\rho_{B}=\mathrm{tr}_{A}\left(\rho\right)$

$\displaystyle
E_{f}(\rho)\leq\overline{E}\leq\min\left[S\left(\rho_{A}\right),S\left(\rho_{B}\right)\right].$
(11)

The dependence of the EAEE on the unravelling becomes especially important
when the QT method is combined with MPS techniques to represent and propagate
the (stochastic) wave function of the quantum many-body system. MPS
representations are efficient as long as the bond dimension $\chi$, and thus
the entanglement, is small. Since the computational cost of an MPS simulation
grows exponentially with the entanglement, choosing the unravelling that
minimizes the EAEE is of paramount importance. A priori it is unclear which
stochastic propagation scheme leads to small values of the EAEE, without
constructing the trajectories Gharibian (2008). We can however predict how
various types of the stochastic propagators affect the change rate of the EAEE
at each instant of time and accordingly optimize it. This motivates a time-
local greedy algorithm that continuously adapts the stochastic propagators to
the conditional wave function as time progresses (Fig. 1a).

More specifically, given the pure state of the system $|\phi^{(k)}(t)\rangle$
of $k^{\rm th}$ trajectory at time $t$, we want to choose the stochastic
propagator that minimizes the expected instantaneous entanglement increase
rate, $\dot{\overline{E}}=d\overline{E}/dt$ (see Fig. 1b). For simplicity we
consider only propagators corresponding to independent measurements of each
bath channel $j$ and jump operators that do not couple the partitions $A$ and
$B$. As outlined above, the stochastic components of the propagator can then
be applied in sequences. This allows to optimize the corresponding change rate
of the EAEE, $\dot{\overline{E}}_{j}$, channel by channel. Remarkably, one can
analytically perform the minimization
$\min_{X_{j}(t)}\left(\dot{\overline{E}}_{j}\right)$ over all stochastic
propagators that correspond to measurements in eigenbases of arbitrary bath
observables $X_{j}(t)=f[b_{j}(t),b_{j}^{\dagger}(t)]$. Moreover, the minimum
is always obtained either for a number measurement,
$X_{j}(t)=b_{j}^{\dagger}(t)b_{j}(t)$, or for a homodyne measurement,
$X_{j}(t)=b_{j}(t)e^{i\varphi_{j}}+b_{j}^{\dagger}(t)e^{-i\varphi_{j}}$. For
the number measurement the EAEE change rate at a channel $j$ is

$\displaystyle\dot{\overline{E}}_{j}^{\textrm{num}}=\mathrm{tr}\left(c_{j}\phi
c_{j}^{\dagger}\right)\log_{2}\left[\mathrm{tr}\left(c_{j}\phi
c_{j}^{\dagger}\right)\right]+\mathrm{tr}\left(\mathrm{tr}_{B}\left(c_{j}\phi
c_{j}^{\dagger}\right)\left\\{\log_{2}\left[\mathrm{tr}_{B}\left(\phi\right)\right]-\log_{2}\left[\mathrm{tr}_{B}\left(c_{j}\phi
c_{j}^{\dagger}\right)\right]\right\\}\right).$ (12)

The measurement of the homodyne current gives instead:

$\displaystyle\dot{\overline{E}}_{j}^{\mathrm{hom}}=\frac{1}{2\ln
2}\Bigg{[}\bigg{|}e^{-i\varphi_{j}}\mathrm{tr}\left(c_{j}\phi\right)+e^{i\varphi_{j}}\mathrm{tr}\left(\phi
c_{j}^{\dagger}\right)\bigg{|}^{2}-\sum_{k,\ell}\frac{\ln\left(\xi_{k}\right)-\ln\left(\xi_{\ell}\right)}{\xi_{k}-\xi_{\ell}}\bigg{|}e^{-i\varphi_{j}}{\langle\xi_{k}|\mathrm{tr}_{B}\left(c_{j}\phi\right)|\xi_{\ell}\rangle}+e^{i\varphi_{j}}{\langle\xi_{k}|\mathrm{tr}_{B}\left(\phi
c_{j}^{\dagger}\right)|\xi_{\ell}\rangle}\bigg{|}^{2}\Bigg{]}.$ (13)

Here $\phi=|\phi_{j}^{(k)}(t)\rangle\langle\phi_{j}^{(k)}(t)|$ [see Eq. (4)]
and
$\mathrm{tr}_{B}\left(\phi\right)|\xi_{\ell}\rangle=\xi_{\ell}|\xi_{\ell}\rangle$.
We refer the reader to the Supplementary Materials (SM, Section III) for the
proof and derivation of the above statements.

This motivates the following entanglement-optimized quantum trajectory (EOQT)
algorithm (see Fig. 2a). For each discrete time step $dt$ propagate the state
$|\phi^{(k)}(t)\rangle\rightarrow|\phi^{(k)}(t+dt)\rangle$ as follows:

  1. 1.

Sequentially for each channel $j$ calculate the EAEE change rates
$\min_{\varphi_{j}}\dot{\overline{E}}_{j}^{\rm hom}$ and
$\dot{\overline{E}}_{j}^{\rm num}$ and update the state according to the
correspondingly optimal stochastic propagator;

  2. 2.

Complete the update of the resulting state by applying the coherent propagator
$e^{-iH_{\rm sys}dt}$.

This EOQT algorithm can be naturally combined with MPS methods, such as the
time evolving block decimation (TEBD) algorithm Daley  _et al._ (2004); White
and Feiguin (2004) (see Fig. 2a). Importantly, the computational cost of the
evaluation and optimization of (12) and (13) for an MPS with bond dimension
$\chi$ is $\mathcal{O}(\chi^{3}d)$, where $d$ is the local Hilbert space
dimension. This should be compared to the cost of the coherent propagation,
which for each of the time steps $dt$ is $\mathcal{O}(\chi^{3}d^{3}n)$. Thus
the EAEE optimization does not significantly add to the cost of the standard
propagation. On the other hand, the potential gain in the simulation
efficiency can be substantial, as we show in the remainder of this paper.

Figure 2: (a) Schematic illustration of the EOQT-MPS algorithm. To propagate a
state from $t$ to $t+dt$ we first apply a “stochastic layer”: We calculate the
unravelling that minimizes the expected EE in the next time step and apply the
corresponding stochastic propagator. This procedure is implemented
sequentially for all jump operators. Then we implement a “deterministic
layer”, propagating the system with $H_{\rm sys}$ via standard methods such as
TEBD. (b) Time dependence of the excess EAEE, $\overline{E}-E_{f}$, of a
continuously monitored Bell pair for the number measurement (blue), the
homodyne unravelling with $\varphi_{j}=0$ (red), and the EOQT (green, number
of trajectories $N=10^{4}$). The inset shows the statistics of measurement
choices of EOQT algorithm. The error bars are denoted by gray filling.

The simplest example illustrating the dependence of $\overline{E}$ on the
unravelling is obtained by considering two qubits initially in a Bell state,
$(|00\rangle+|11\rangle)/\sqrt{2}$, that are undriven ($H_{\rm sys}=0$) and
coupled to a bath with a jump operator $c_{j}=|1\rangle_{j}\langle 1|$ and
strength $\gamma_{j}=\gamma$ ($j=1,2$). In this case, $\overline{E}$ can be
analytically calculated for different unravellings. For the photon number
measurements, one obtains $\overline{E}^{\mathrm{num}}=\sigma(2\gamma t)$,
where we introduced the function:

$\displaystyle\sigma(s)=\frac{1}{2\ln
2}\left[\left(1+e^{-s}\right)\ln{\left(1+e^{-s}\right)}+se^{-s}\right].$ (14)

The expression for homodyne measurement is

$\displaystyle\overline{E}^{\mathrm{hom}}(\tau)$
$\displaystyle=\frac{1}{2\sqrt{2\pi\tau}}\int
ds\sigma(s)e^{-\frac{(s-2\tau)^{2}}{8\tau}}$ (15)

with $\tau=\gamma t\left(\cos^{2}\varphi_{1}+\cos^{2}\varphi_{2}\right)$. In
Fig. 2 we compare the excess entanglements, $\overline{E}-E_{f}$, for number,
homodyne, and optimal (EOQT) measurements. In this case $E_{f}$ evaluates to
$E_{f}(t)=-r_{+}\log_{2}(r_{+})-r_{-}\log_{2}(r_{-}),$ where
$r_{\pm}=\left(1\pm\sqrt{1-e^{-2\gamma t}}\right)/{2}$ Wootters (1998). We
note that for short times the homodyne unravelling with $\varphi_{j}=0$
($j=1,2$) saturates this bound, while for longer times the number unravelling
approaches $E_{f}$ faster (Fig. 2b). The optimized algorithm indeed finds
unravellings resulting in an EAEE that is always close to the theoretical
minimum, $E_{f}$. As visible in the inset of Fig. 2b, the optimizer chooses
mostly the homodyne propagator with $\varphi_{j}=0$ for early times and
switches over to predominantly choosing the number propagators around $\gamma
t\approx 1$.

To demonstrate the potential of our approach in a many-body setting, we
consider a one-dimensional open random Brownian circuit (RBC) for a chain of
spin-1/2 particles. We choose this example since the evolution under the RBC
leads to rapid EE growth. The coherent part of the evolution is given by the
time-dependent RBC Hamiltonian:

$\displaystyle H_{\rm
sys}(t)=\sum_{j=1}^{n-1}\sum_{k,\ell=0}^{3}g_{j}^{k,\ell}(t)\sigma_{j}^{k}\otimes\sigma_{j+1}^{\ell}.$
(16)

Here $n$ is the number of spins and
$\sigma_{j}^{k}\in\\{\mathbb{1}_{j},\sigma_{j}^{x},\sigma_{j}^{y},\sigma_{j}^{z}\\}$
are standard Pauli operators acting on spin $j$. The parameters
$g_{j}^{k,\ell}(t)$ are Gaussian stochastic variables with $\langle\langle
g_{j}^{k,\ell}(t)\rangle\rangle=0$ and $\langle\langle
g_{j}^{k,\ell}(t)g_{j^{\prime}}^{k^{\prime},\ell^{\prime}}(t^{\prime})\rangle\rangle=\alpha\delta_{j,j^{\prime}}\delta_{k,k^{\prime}}\delta_{\ell,\ell^{\prime}}\delta(t-t^{\prime})$,
where $\langle\langle\cdots\rangle\rangle$ denotes the average over
Hamiltonian realizations and $\alpha$ is the variance of RBC. The incohenert
part is given by the jump operators $c_{j}=\sigma_{j}^{z}$ with uniform decay
strength (measurement rate), $\gamma_{j}=\gamma$.

Figure 3: EAEE for the open RBC. (a) EAEE profile in the long time limit for
various unravellings. Solid lines denote homodyne unravellings with changing
phase from $\varphi_{j}=0$ (red) to $\varphi_{j}=\pi/2$ (grey) in increments
of $\pi/20$, dashed blue and green lines correspond to the number and EOQT
measurements. The insets show time evolution of EAEE for various unravellings
and measurement choices of the EOQT. Here $m=n=16$, $\chi=128$. (b) Half-chain
EAEE for larger systems (with $\chi=512$). Depending on the homodyne phase,
the ensembles display an area law at small $\varphi_{j}$ and a volume law at
$\varphi_{j}\approx\pi/2$. The EOQT method results in an area-law EAEE close
to the homodyne unravelling with $\varphi_{j}=0$, whereas the number
unravelling leads to an area law with larger EAEE. The phase diagram is shown
in the inset. For all data in this figure we used $\alpha/\gamma=0.1$
($\alpha=1$) and $N=200$. The statistical error bars lie within the lines’
thickness.

For sufficiently strong $\gamma$ this model exhibits a type of measurement-
induced phase transition that is solely driven by the type of unravelling (see
Fig. 3). Depending on the bath measurement operators, the character of EAEE
changes from area law to volume law. That is, even though all types of bath
measurements solve in the end the same ME (1), some unravellings lead to
ensembles with area law EAEE that can be efficiently computed on classical
computers, while other unravellings fail to do so.

This change from volume to area law can be shown explicitly by considering the
stochastic propagator for a homodyne unravelling, $K_{j}$, defined above [see
Eq. (S136)]. Since
${\sigma_{j}^{z}}^{\dagger}{\sigma_{j}^{z}}={\sigma_{j}^{z}}^{2}=\mathbb{1}$,
it takes the form 222See Section VI of the SM for more details.

$\displaystyle
K_{j}\propto\exp\left\\{e^{i\varphi_{j}}\sqrt{\gamma}\sigma^{z}_{j}d\xi_{j}(t)\right\\},$
(17)

where
$d\xi_{j}(t)=2\sqrt{\gamma}\left<\sigma^{z}_{j}\right>(t)\cos{\varphi_{j}}dt+dW_{j}(t)$.
It has a unitary part,
$\exp\left\\{i\sin(\varphi_{j})\sqrt{\gamma}\sigma^{z}_{j}d\xi_{j}(t)\right\\}$,
and a non-unitary part,
$\exp\left\\{\cos(\varphi_{j})\sqrt{\gamma}\sigma^{z}_{j}d\xi_{j}(t)\right\\}$.
The unitary component can be absorbed into the coherent part of the evolution,
as it leaves the ensemble of RBCs invariant. The non-unitary component can be
re-interpreted a stochastic propagator of a Markov process with an effective
decoherence rate $\gamma_{j}^{\rm eff}=\gamma\cos^{2}\varphi_{j}$. Changing
the unravelling at fixed $\gamma$ by changing $\varphi_{j}$ is thus equivalent
to changing $\gamma$ for a fixed unravelling. This equivalence is remarkable,
as these two interpretations describe two profoundly different physical
scenarios: The first refers to various representations of the solution of a
single ME, while the second refers to particular solutions of various MEs. The
latter has been studied extensively recently Li  _et al._ (2018, 2019);
Skinner  _et al._ (2019); Chan  _et al._ (2019); Gullans and Huse (2020); Choi
_et al._ (2020); Bao  _et al._ (2020); Ippoliti  _et al._ (2021); Zabalo  _et
al._ (2020, 2021); Botzung  _et al._ (2021); Müller  _et al._ (2021) and it is
known that the conditional states can undergo a so-called measurement-induced
phase transitions from area law to volume law entanglement, depending on the
coupling strength to the environment. In our case, this phase transition
occurs as a function of the unravelling, e.g. parametrized by $\varphi_{j}$.
Indeed, for $\varphi_{j}=\pi/2$ the QTs map to fully coherent RBC evolution,
which generates entanglement that obeys a volume law. In the other limiting
case, $\varphi_{j}=0$ (and sufficiently large $\gamma$), the measurement back-
action continuously leads to an effective projection onto product states,
resulting in entanglement within the system that satisfies an area law. These
two phases are separated by a critical point at
$\gamma\cos^{2}\varphi_{j}=\gamma_{\rm crit}$. Thus, when $\gamma>\gamma_{\rm
crit}$, the many-body ME can be efficiently solved using MPS via QTs if an
efficient unravelling is chosen. In Fig. 3 we show results of MPS simulations
of various unravellings of the open RBC demonstrating that the EOQT method
indeed finds efficiently computable unravelling.

In summary, we introduced a novel QT method to simulate the time evolution of
noisy quantum many-body system and showed that our method can enable an
efficient classical simulation, where standard QT methods may fail. While in
this work we focused on monitoring and optimizing each output channel
separately, more general schemes, based on collective monitoring and
optimization can be developed Holland  _et al._ (1996). Detailed performance
comparisons with other techniques for open quantum many-body systems Bonnes
and Läuchli (2014); Cui  _et al._ (2015); Werner  _et al._ (2016); Jin  _et
al._ (2016); White  _et al._ (2018); Biella  _et al._ (2018); Schwarz  _et
al._ (2018); Hartmann and Carleo (2019); Vicentini  _et al._ (2019); Yoshioka
and Hamazaki (2019); Nagy and Savona (2019); Weimer  _et al._ (2021) are left
for future work.

###### Acknowledgements.

We thank P. Zoller and B. Kraus for stimulating discussions. This work is
supported by an ESQ Discovery Grant and an ERC Starting Grant (No. 101041435).

## References

  * Preskill (2018) J. Preskill, Quantum 2, 79 (2018).
  * White (1992) S. R. White, Phys. Rev. Lett. 69, 2863 (1992).
  * Vidal (2003) G. Vidal, Phys. Rev. Lett. 91, 147902 (2003).
  * Schollwöck (2011) U. Schollwöck, Ann. Phys. 326, 96 (2011).
  * 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).
  * 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).
  * 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).
  * Li  _et al._ (2018) Y. Li, X. Chen, and M. P. A. Fisher, Phys. Rev. B 98, 205136 (2018).
  * Li  _et al._ (2019) Y. Li, X. Chen, and M. P. A. Fisher, Phys. Rev. B 100, 134306 (2019).
  * 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).
  * Gullans and Huse (2020) M. J. Gullans and D. A. Huse, Phys. Rev. Lett. 125, 070606 (2020).
  * Choi  _et al._ (2020) S. Choi, Y. Bao, X.-L. Qi, and E. Altman, Phys. Rev. Lett. 125, 030505 (2020).
  * Bao  _et al._ (2020) Y. Bao, S. Choi, and E. Altman, Phys. Rev. B 101, 104301 (2020).
  * Zabalo  _et al._ (2020) A. Zabalo, M. J. Gullans, J. H. Wilson, S. Gopalakrishnan, D. A. Huse, and J. H. Pixley, Phys. Rev. B 101, 060301 (2020).
  * Ippoliti  _et al._ (2021) M. Ippoliti, M. J. Gullans, S. Gopalakrishnan, D. A. Huse, and V. Khemani, Phys. Rev. X 11, 011030 (2021).
  * Zabalo  _et al._ (2021) A. Zabalo, M. J. Gullans, J. H. Wilson, R. Vasseur, A. W. W. Ludwig, S. Gopalakrishnan, D. A. Huse, and J. H. Pixley, arXiv:2107.03393 (2021).
  * Botzung  _et al._ (2021) T. Botzung, S. Diehl, and M. Müller, arXiv:2106.10092 (2021).
  * Müller  _et al._ (2021) T. Müller, S. Diehl, and M. Buchhold, arXiv:2105.08076 (2021).
  * Lindblad (1976) G. Lindblad, Comm. Math. Phys. 48, 119 (1976).
  * Daley (2014) A. J. Daley, Adv. Phys. 63, 77 (2014).
  * Note (1) Since the trajectories are by construction sampled according to Born’s rule, all such strategies, including the adaptive time- and history-dependent ones, generate trajectories that provide valid solutions of the ME [see Section I of the Supplementary Materials (SM) for more detail].
  * Wootters (1998) W. K. Wootters, Phys. Rev. Lett. 80, 2245 (1998).
  * Gharibian (2008) S. Gharibian, arXiv:0810.4507 (2008).
  * Daley  _et al._ (2004) A. J. Daley, C. Kollath, U. Schollwöck, and G. Vidal, J. Stat. Mech.: Theory Exp , P04005 (2004).
  * White and Feiguin (2004) S. R. White and A. E. Feiguin, Phys. Rev. Lett. 93, 076401 (2004).
  * Note (2) See Section VI of the SM for more details.
  * Holland  _et al._ (1996) M. Holland, S. Marksteiner, P. Marte, and P. Zoller, Phys. Rev. Lett. 76, 3683 (1996).
  * Bonnes and Läuchli (2014) L. Bonnes and A. M. Läuchli, arXiv:1411.4831 (2014).
  * Cui  _et al._ (2015) J. Cui, J. I. Cirac, and M. C. Bañuls, Phys. Rev. Lett. 114, 220601 (2015).
  * 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).
  * Jin  _et al._ (2016) J. Jin, A. Biella, O. Viyuela, L. Mazza, J. Keeling, R. Fazio, and D. Rossini, Phys. Rev. X 6, 031011 (2016).
  * White  _et al._ (2018) C. D. White, M. Zaletel, R. S. K. Mong, and G. Refael, Phys. Rev. B 97, 035127 (2018).
  * Biella  _et al._ (2018) A. Biella, J. Jin, O. Viyuela, C. Ciuti, R. Fazio, and D. Rossini, Phys. Rev. B 97, 035103 (2018).
  * Schwarz  _et al._ (2018) F. Schwarz, I. Weymann, J. von Delft, and A. Weichselbaum, Phys. Rev. Lett. 121, 137702 (2018).
  * Hartmann and Carleo (2019) M. J. Hartmann and G. Carleo, Phys. Rev. Lett. 122, 250502 (2019).
  * Vicentini  _et al._ (2019) F. Vicentini, A. Biella, N. Regnault, and C. Ciuti, Phys. Rev. Lett. 122, 250503 (2019).
  * Yoshioka and Hamazaki (2019) N. Yoshioka and R. Hamazaki, Phys. Rev. B 99, 214306 (2019).
  * Nagy and Savona (2019) A. Nagy and V. Savona, Phys. Rev. Lett. 122, 250501 (2019).
  * Weimer  _et al._ (2021) H. Weimer, A. Kshetrimayum, and R. Orús, Rev. Mod. Phys. 93, 015008 (2021).

Entanglement-Optimal Trajectories of Many-Body Quantum Markov Processes:
Supplementary material

This Supplementary Material is structured as follows:

  * •

In Section I we discuss adaptive quantum trajectory (QT) schemes on general
grounds to confirm that such schemes indeed allow for an unbiased solution of
master equations (MEs).

  * •

In Section II and in Section III we elaborate on technical details underlying
the main theoretical results of our paper: We derive the the formulae for the
ensemble-averaged entanglement entropy (EAEE) change rate for arbitrary
environment monitoring schemes, obtain the expressions for number and homodyne
measurement schemes, and prove the optimality result stated in the main text.

  * •

In Section IV we provide details about the matrix product state (MPS)
implementation of our entanglement-optimized quantum trajectory (EOQT)
algorithm.

  * •

In Section V we apply our algorithm to several physically relevant examples of
MEs that are chosen to illustrate, first, the convergence properties of the
different QT schemes, and, second, the superior performance of the EOQT method
compared to other schemes.

  * •

Section VI contains technical derivations underlying the discussion of the
open random Brownian circuit (RBC) in the main text.

##  I Adaptive quantum trajectory methods

In this section we discuss general aspects regarding the validity of adaptive
QT methods (such as the EOQT method introduced in the main text). While non-
adaptive QT methods have been extensively discussed in the literature,
adaptive schemes have received far less attention. For adaptive schemes it is
sometimes not immediately obvious that no bias is introduced. Our aim in this
section is therefore to clarify that also adaptive continuous monitoring
schemes give rise to valid trajectories, and thus that (properly constructed)
adaptive QT methods lead to valid, unbiased solutions of the ME.

We find it easiest to discuss this on very general grounds in the following.
Let us consider a generic state $|\Psi\rangle$ of a system $S$ and its
environment $E$. For the sake of this discussion, let us assume that the
environment consists of a collection of several degrees of freedom, such that
we can write $E=e_{1}\otimes e_{2}\otimes e_{3}\otimes\dots$. Note that this
includes the generic quantum-optical setting discussed in the main text, where
each of these degrees of freedom is represented by a bosonic mode with
creation and annihilation operators $b(t)$ and $b^{\dagger}(t)$ (for
notational simplicity we suppress the index $j$ here). The density operator
$\rho$, describing the reduced state of the system $S$, is obtained by tracing
over the degrees of freedom of $E$,
$\rho=\textrm{tr}_{E}\left(|\Psi\rangle\langle\Psi|\right)$. The notion of the
density operator implies an inherent statistical interpretation: It describes
the full measurement statistics an experimenter (called Alice) obtains from
repeated preparations of the same state followed by measurements of
observables in $S$, where any knowledge of $E$ is disregarded.

Now, instead of disregarding $E$ for each such repetition, let us assume that
another observer (called Bob) measures the degrees of freedom of the
environment. For each repetition $i$, Bob obtains a set of measurement
outcomes $\\{x_{1}^{(i)},x_{2}^{(i)},\dots\\}$, projecting the environment
degrees of freedom onto the state $|x_{1}^{(i)},x_{2}^{(i)},\dots\rangle$.
Here, $x_{\ell}^{(i)}$ denotes the measurement outcome of the $\ell$-th degree
of freedom of $E$ in the $i$-th run of the experiment. Importantly, the
probabilities for each outcome are determined by _Born’s rule_ , such that
$p^{(i)}\equiv p(x_{1}^{(i)},x_{2}^{(i)}\dots)=\textrm{tr}(\langle
x_{1}^{(i)},x_{2}^{(i)},\dots|\Psi\rangle\langle\Psi|{x_{1}^{(i)},x_{2}^{(i)},\dots}\rangle)$.
Thus, Bob obtains a different description of the state $S$ for each run,
namely, the state of $S$ conditional on the measurement outcomes of the
environment:

$\displaystyle|\psi^{(i)}\rangle\equiv|\psi(x_{1}^{(i)},x_{2}^{(i)}\dots)\rangle\equiv\frac{\langle
x_{1}^{(i)},x_{2}^{(i)},\dots|\Psi\rangle}{p(x_{1}^{(i)},x_{2}^{(i)}\dots)^{1/2}}.$
(S1)

At first glance, this is in conflict with the description of the state of the
system $S$ by Alice which assigns the density operator $\rho$ to every single
repetition. The conflict is of course resolved by noting that the measurement
of the environment precisely results in samples $|\psi^{(i)}\rangle$ that are
distributed according to the probability $p^{(i)}$, which form a decomposition
of $\rho=\sum_{i}p^{(i)}|\psi^{(i)}\rangle\langle\psi^{(i)}|$. Thus, upon
averaging over repetitions, Bob obtains the same description of the state as
Alice. Mathematically, this is guaranteed by Born’s rule:

$\rho=\mathrm{tr}_{E}\left(|\Psi\rangle\langle\Psi|\right)=\sum_{x_{1},x_{2},\dots}\langle
x_{1},x_{2},\dots|\Psi\rangle\langle\Psi|x_{1},x_{2},\dots\rangle=\sum_{x_{1},x_{2},\dots}p(x_{1},x_{2},\dots)|\psi(x_{1},x_{2},\dots)\rangle\langle\psi(x_{1},x_{2},\dots)|.$

This is the core idea behind QT methods, which numerically simulate the
measurements of the environment by randomly drawing the measurement outcomes
$\\{x_{1}^{(i)},x_{2}^{(i)}\dots\\}$ according to Born’s probability
distribution, thus obtaining samples $|\psi^{(i)}\rangle$ according to the
distribution $p^{(i)}$.

A subtle but crucial point in the above discussion is that the description of
the state of $S$ by Alice cannot (and does not) depend on the way that Bob
measures the environment degrees of freedom. In particular, Bob can apply
arbitrarily involved measurement strategies, including adaptive ones: While
Bob’s conditional states will depend on his measurement strategy, the ensemble
of states he obtains _due to Born’s rule_ is by construction always a valid
decomposition of the reduced density operator. In fact, if this was not the
case, Bob could send a signal to Alice by simply measuring the environment
degrees of freedom, which is not possible by the basic principles of quantum
mechanics.

Mathematically, we can show this as follows. Consider the following generic,
adaptive, strategy for Bob: Instead of measuring all the environment degrees
of freedom in a fixed basis, he chooses to first measure the first degree of
freedom, $e_{1}$, in a given basis, spanned by states $|x_{1}\rangle$ (where
$x_{1}=\\{1,2,\dots\dim(e_{1})\\}$ labels basis states and potential
measurement outcomes). The measurement outcome of this first measurement is
random and distributed according to Born’s rule,
$p(x_{1})=\mathrm{tr}\left(\langle
x_{1}|\Psi\rangle\langle\Psi|x_{1}\rangle\right)$. Bob then chooses to measure
the second degree of freedom, $e_{2}$, in a basis that depends on the
measurement outcome of the measurement, $e_{1}$. We label these basis states
by $|x_{2}(x_{1})\rangle$ (where $x_{2}=\\{1,\dots\dim(e_{2})\\}$). This
notation explicitly indicates that the measurement basis for $e_{2}$ is
determined by the measurement result of $e_{1}$, according to an _adaptive_
strategy specified by Bob. Of course, the outcome of $e_{2}$ is again random
and distributed according to the probability distribution that can be obtained
from Born’s rule:

$p(x_{2}|x_{1})=\mathrm{tr}\left(\langle
x_{1},x_{2}(x_{1})|\Psi\rangle\langle\Psi|x_{1},x_{2}(x_{1})\rangle\right)/p(x_{1})=p(x_{1},x_{2})/p(x_{1}).$

Bob then continues and measures sequentially each environment degree of
freedom, $e_{j}$, in a basis that he chooses adaptively, based on previous
measurement outcomes, i.e. a basis spanned by states
$|x_{j}(x_{j-1},\dots,x_{2},x_{1})\rangle$ (where
$x_{j}\in\\{1,2,\dots\dim(e_{j})\\}$, and we again highlight the dependence of
the basis on previous results). For each $e_{j}$, Born’s rule dictates the
probability distribution of the measurement outcomes. By measuring all $n$
environment degrees of freedom, Bob obtains the state of $S$ conditional on
the sequence of measurement outcomes:

$\displaystyle\frac{\langle x_{1},x_{2}(x_{1}),x_{3}(x_{2},x_{1}),\dots
x_{n}(x_{n-1},\dots,x_{2},x_{1})|\Psi\rangle}{\sqrt{p(x_{1},x_{2},x_{3},\dots,x_{n})}}$
(S2)

with a probability given by Born’s rule as

$\displaystyle p(x_{1},x_{2},x_{3},\dots,x_{n})=\mathrm{tr}\left[{\langle
x_{1},x_{2}(x_{1}),\dots
x_{n}(x_{n-1},\dots,x_{2},x_{1})|\Psi\rangle\langle\Psi|x_{1},\dots
x_{n}(x_{n-1},\dots,x_{2},x_{1})\rangle}\right].$ (S3)

We stress that the above expressions describe the adaptive measurement
strategy of Bob, which includes the one considered in our work. We can now
easily show that the ensemble of states that Bob obtains forms a valid
decomposition of $\rho$:

$\displaystyle\sum_{x_{1},\dots,x_{n}}p(x_{1},x_{2},x_{3},\dots,x_{n})\frac{\langle
x_{1},x_{2}(x_{1}),\dots
x_{n}(x_{n-1},\dots,x_{1})|\Psi\rangle\langle\Psi|x_{1},x_{2}(x_{1}),\dots
x_{n}(x_{n-1},\dots,x_{1})\rangle}{p(x_{1},x_{2},x_{3},\dots,x_{n})}$ (S4)
$\displaystyle=\sum_{x_{1},\dots,x_{n-1}}\sum_{x_{n}}\langle
x_{1},x_{2}(x_{1}),\dots
x_{n}(x_{n-1},\dots,x_{1})|\Psi\rangle\langle\Psi|x_{1},x_{2}(x_{1}),\dots
x_{n}(x_{n-1},\dots,x_{1})\rangle$ (S5)
$\displaystyle=\sum_{x_{1},\dots,x_{n-1}}\mathrm{tr}_{{e_{n}}}\left[{\langle
x_{1},x_{2}(x_{1}),\dots
x_{n-1}(x_{n-2},\dots,x_{1})|\Psi\rangle\langle\Psi|x_{1},x_{2}(x_{1}),\dots
x_{n-1}(x_{n-2},\dots,x_{1})\rangle}\right]$ (S6)
$\displaystyle=\sum_{x_{1},\dots,x_{n-2}}\mathrm{tr}_{{e_{n-1},e_{n}}}\left[{\langle
x_{1},x_{2}(x_{1}),\dots
x_{n-2}(x_{n-3},\dots,x_{1})|\Psi\rangle\langle\Psi|x_{1},x_{2}(x_{1}),\dots
x_{n-2}(x_{n-3},\dots,x_{1})\rangle}\right]$ (S7)
$\displaystyle=\dots=\mathrm{tr}_{e_{1},\dots,e_{n-1},e_{n}}\left(|\Psi\rangle\langle\Psi|\right)\equiv\mathrm{tr}_{E}\left(|\Psi\rangle\langle\Psi|\right)\equiv\rho.$
(S8)

Here we used the fact that the trace over the degrees of freedom can be taken
in an arbitrary basis: $\mathrm{tr}_{e_{j}}\left(O\right)=\sum_{x_{j}}\langle
x_{j}(x_{j-1},\dots x_{1})|O|x_{j}(x_{j-1},\dots x_{1})\rangle$, which holds
for any $\left(x_{j-1},\dots x_{1}\right)$, since the states $\langle
x_{j}(x_{j-1},\dots x_{1})|$ form a orthonormal basis for every
$(x_{j-1},\dots x_{1})$. In consequence, simulating such an adaptive
measurement scheme to obtain QTs is equivalent to sampling pure states from a
valid decomposition of the density operator (as long as the simulated
measurements respect Born’s rule).

###  I.1 Toy example

The above discussion is perhaps best illustrated by the following toy example.
Consider a quantum system $S$ that consists of two qubits initially in the
maximally entangled state
$|\psi_{S}\rangle=\frac{1}{2}\left(|00\rangle+|01\rangle+|10\rangle-|11\rangle\right)$
and the environment $E$ consisting of two qubits in a product state
$|00\rangle$ (see Fig. S1).

Figure S1: Illustration of a toy example.

Imagine that two consecutive CNOT-operations $U_{T}$ were applied to first and
second pair of system + environment qubits, as depicted in Fig. S1. These
operations can be considered as elementary joint evolution of $S+E$. One can
check that the state of $S$ after the evolution is a maximally mixed state:

$\displaystyle\rho=\mathrm{tr}_{E}\left(U_{T}|\psi_{S}\rangle\langle\psi_{S}|\otimes|00\rangle\langle
00|U_{T}^{\dagger}\right)=\mathbb{1}/4.$ (S9)

Now let us imagine that Bob measures the environment in some basis
$|x_{i}x_{j}\rangle$ according to one of the three following strategies.

####  I.1.1 First strategy

Let us first consider Bob uses the computational basis, i.e.
$\\{x_{i}x_{j}\\}=\left\\{00,01,10,11\right\\}$. This may be interpreted as a
minimal model for a number measurement in the QT language. Depending on the
measurement outcomes the conditional state of $S$ is

$\displaystyle|\psi_{00}\rangle=|00\rangle\textrm{ with probability
}p_{00}={1}/{4};$ (S10) $\displaystyle|\psi_{01}\rangle=|01\rangle\textrm{
with probability }p_{01}={1}/{4};$ (S11)
$\displaystyle|\psi_{10}\rangle=|10\rangle\textrm{ with probability
}p_{10}={1}/{4};$ (S12) $\displaystyle|\psi_{11}\rangle=|11\rangle\textrm{
with probability }p_{11}={1}/{4}.$ (S13)

In this case each of the measurement outcomes is a non-entangled product
state. Summing the pure-state projectors with corresponding weights
(probabilities) gives, as expected, the density matrix of the system after the
evolution [as in Eq. (S9)]:

$\displaystyle\frac{1}{4}\bigg{(}|00\rangle\langle 00|+|01\rangle\langle
01|+|10\rangle\langle 10|+|11\rangle\langle 11|\bigg{)}=\mathbb{1}/4.$ (S14)

####  I.1.2 Second strategy

Now let us consider instead Bob measuring the $E$ qubits in the states
$|\pm\rangle=\frac{1}{\sqrt{2}}(|0\rangle\pm|1\rangle)$, which gives the
following span of measurement outcomes:
$\\{x_{i}x_{j}\\}=\left\\{++,+-,-+,--\right\\}$. This may be considered as a
minimal model analogous to a homodyne measurement in the QT language. For each
of these measurement outcomes, one can also find the corresponding conditional
states

$\displaystyle|\psi_{++}\rangle=\frac{1}{4}\bigg{(}|00\rangle+|01\rangle+|10\rangle-|11\rangle\bigg{)}\textrm{
with probability }p_{++}={1}/{4};$ (S15)
$\displaystyle|\psi_{+-}\rangle=\frac{1}{4}\bigg{(}|00\rangle-|01\rangle+|10\rangle+|11\rangle\bigg{)}\textrm{
with probability }p_{+-}={1}/{4};$ (S16)
$\displaystyle|\psi_{-+}\rangle=\frac{1}{4}\bigg{(}|00\rangle+|01\rangle-|10\rangle+|11\rangle\bigg{)}\textrm{
with probability }p_{-+}={1}/{4};$ (S17)
$\displaystyle|\psi_{--}\rangle=\frac{1}{4}\bigg{(}|00\rangle-|01\rangle-|10\rangle-|11\rangle\bigg{)}\textrm{
with probability }p_{--}={1}/{4}.$ (S18)

As one can see, each of the measurement outcomes is a maximally-entangled
state. Again, the weighted sum gives the the expected maximally mixed state:

$\displaystyle\frac{1}{4}\bigg{(}|\psi_{++}\rangle\langle\psi_{++}|+|\psi_{+-}\rangle\langle\psi_{+-}|+|\psi_{-+}\rangle\langle\psi_{-+}|+|\psi_{--}\rangle\langle\psi_{--}|\bigg{)}=\mathbb{1}/4.$
(S19)

####  I.1.3 Third strategy

Now let us turn to a simple example of adaptive measurement. Imagine that Bob
first measures the first environment qubit and, depending on the outcome,
chooses the measurement basis for the second qubit. For example, if the first
qubit, measured in the computational basis, results in a outcome $0$, then the
second qubit is also measured in the computational basis. On the other hand,
if the first qubit gives $1$ as an outcome, the second qubit should be
measured in “$\pm$”-basis. This situation produces the following span of
outcomes: $\\{x_{i}x_{j}\\}=\left\\{00,01,1+,1-\right\\}$. As before, we can
find the corresponding trajectories:

$\displaystyle|\psi_{00}\rangle=\frac{1}{2}|00\rangle\textrm{ with probability
}p_{00}={1}/{4};$ (S20)
$\displaystyle|\psi_{01}\rangle=\frac{1}{2}|01\rangle\textrm{ with probability
}p_{01}={1}/{4};$ (S21)
$\displaystyle|\psi_{1+}\rangle=\frac{1}{2\sqrt{2}}\bigg{(}|10\rangle-|11\rangle\bigg{)}\textrm{
with probability }p_{1+}={1}/{4};$ (S22)
$\displaystyle|\psi_{1-}\rangle=\frac{1}{2\sqrt{2}}\bigg{(}|10\rangle+|11\rangle\bigg{)}\textrm{
with probability }p_{1-}={1}/{4}.$ (S23)

And, as before, this adaptive measurement gives the expected state of the
system, the maximally-mixed density matrix:

$\displaystyle\frac{1}{4}\bigg{(}|\psi_{00}\rangle\langle\psi_{00}|+|\psi_{01}\rangle\langle\psi_{01}|+|\psi_{1+}\rangle\langle\psi_{1+}|+|\psi_{1-}\rangle\langle\psi_{1-}|\bigg{)}=\mathbb{1}/4.$
(S24)

As one can see, independent of the measurement strategy, Bob always obtains
pure samples (trajectories) from various decompositions of the same density
operator. This example further illustrates that different measurement
strategies lead to trajectories with different entanglement properties.

##  II Entanglement of conditional ensembles

In this section we give a general expression of the average entanglement
entropy of an ensemble of conditional states, which will allow us to calculate
the EAEE change rate in QT in the next section. This forms the basis of one of
the main results of our work.

We consider a closed tripartite system $A\otimes B\otimes M$. The state of the
full system $A\otimes B\otimes M$ can be written as

$\displaystyle|\Psi\rangle=\sum_{\textbf{x}}|\tilde{\phi}(\textbf{x})\rangle|\textbf{x}\rangle,$
(S25)

where $|\tilde{\phi}(\textbf{x})\rangle$ is the unnormalized state of $A$ and
$B$, conditional on $M$ being in state $|\textbf{x}\rangle$. Here the states x
form an orthonormal basis of $M$,
$\langle\textbf{x}|\textbf{y}\rangle=\delta_{\textbf{x},\textbf{y}}$. That is,
if $M$ is measured in this basis, the system $A\cup B$ is projected onto the
state

$\displaystyle|\tilde{\phi}(\textbf{x})\rangle=\langle\textbf{x}|\Psi\rangle\textrm{
with probability
}p(\textbf{x})=\langle\tilde{\phi}(\textbf{x})|\tilde{\phi}(\textbf{x})\rangle.$
(S26)

It is convenient to introduce normalized states
$|\phi(\textbf{x})\rangle=|\tilde{\phi}(\textbf{x})\rangle/\sqrt{p(\textbf{x})}$,
and the notation $\Psi=|\Psi\rangle\langle\Psi|$,
$\tilde{\phi}(\textbf{x})=|\tilde{\phi}(\textbf{x})\rangle\langle\tilde{\phi}(\textbf{x})|$
and $\phi(\textbf{x})=|\phi(\textbf{x})\rangle\langle\phi(\textbf{x})|$.
Importantly, the ensemble of states (S26) depends on the choice of basis in
which $M$ is measured. Of course the density operator of $A\cup B$, i.e. the
ensemble average of the conditional states $\phi(\textbf{x})$, is manifestly
invariant under changes of the measurement basis:

$\displaystyle\rho=\sum_{\textbf{x}}p(\textbf{x})\phi(\textbf{x})=\sum_{\textbf{x}}\tilde{\phi}(\textbf{x})=\mathrm{tr}_{{M}}\left(\Psi\right),$
(S27)

However, the ensemble average of any non-linear functional of the conditional
states in general does depend on the basis in which $M$ is measured.

We are interested in the entanglement properties of the ensemble of states
(S26), specifically the EAEE of the bipartite splitting of $A$ and $B$:

$\displaystyle\overline{E}=\sum_{\textbf{x}}p(\textbf{x})E\left[\phi(\textbf{x})\right].$
(S28)

To proceed, it is convenient to define the operator

$\displaystyle\mu=\sum_{\textbf{x}}\tilde{\phi}(\textbf{x})\otimes|\textbf{x}\rangle\langle\textbf{x}|=\sum_{\textbf{x}}p(\textbf{x})\phi(\textbf{x})\otimes|\textbf{x}\rangle\langle\textbf{x}|.$
(S29)

Note that $\mu$ is dependent on the specific choice of basis in $M$. We define
$\mu_{M}=\mathrm{tr}_{A,B}\left(\mu\right)=\sum_{\textbf{x}}p(\textbf{x})|\textbf{x}\rangle\langle\textbf{x}|$,
$\mu_{A,M}=\mathrm{tr}_{B}\left(\mu\right)=\sum_{\textbf{x}}p(\textbf{x})\rho_{A}(\textbf{x})\otimes|\textbf{x}\rangle\langle\textbf{x}|$,
and $\rho_{A}(\textbf{x})=\mathrm{tr}_{B}\left(\phi(\textbf{x})\right)$. With
this we define the conditional entropy of the state $\mu$ as

$\displaystyle S(A|M)$ $\displaystyle=S(\mu_{AM})-S(\mu_{M}),$ (S30)

where $S(X)=-\mathrm{tr}\left[X\log_{2}\left(X\right)\right]$ is the von
Neumann entropy. The conditional entropy is an important quantity since it is
precisely equivalent to the EAEE of the ensemble (S26):

$\displaystyle S(A|M)$
$\displaystyle=-\mathrm{tr}\left[{\sum_{\textbf{x}}p(\textbf{x})\rho_{A}(\textbf{x})\otimes|\textbf{x}\rangle\langle\textbf{x}|\log_{2}\sum_{\textbf{y}}p(\textbf{y})\rho_{A}(\textbf{y})\otimes|\textbf{y}\rangle\langle\textbf{y}|}\right]+\mathrm{tr}\left[\sum_{\textbf{x}}p(\textbf{x})|\textbf{x}\rangle\langle\textbf{x}|\log_{2}\sum_{\textbf{y}}p(\textbf{y})|\textbf{y}\rangle\langle\textbf{y}|\right]$
(S31)
$\displaystyle=-\sum_{\textbf{x}}p(\textbf{x})\mathrm{tr}\left\\{{\rho_{A}(\textbf{x})\log_{2}\left[p(\textbf{x})\rho_{A}(\textbf{x})\right]}\right\\}+\sum_{\textbf{x}}p(\textbf{x})\log_{2}p(\textbf{x})$
(S32)
$\displaystyle=-\sum_{\textbf{x}}p(\textbf{x})\left\\{{\mathrm{tr}\left[{\rho_{A}(\textbf{x})\log_{2}\rho_{A}(\textbf{x})}\right]+\mathrm{tr}\left[{\rho_{A}(\textbf{x})\log_{2}p(\textbf{x})}\right]-\log_{2}p(\textbf{x})}\right\\}$
(S33)
$\displaystyle=-\sum_{\textbf{x}}p(\textbf{x})\mathrm{tr}\left[{\rho_{A}(\textbf{x})\log_{2}\rho_{A}(\textbf{x})}\right]=\sum_{\textbf{x}}p(\textbf{x})S\left[\rho_{A}(\textbf{x})\right]=\sum_{\textbf{x}}p(\textbf{x})E\left[\phi(\textbf{x})\right]\equiv\overline{E}.$
(S34)

Note that, since $S(\rho_{A})=S(\rho_{B})$ for a pure state on $A\cup B$, we
have $S(A|M)=S(B|M)$. In the following sections we apply this to derive
expressions for the change rate of the EAEE in a QT simulation.

##  III Entanglement entropy change rates in Monitored systems

###  III.1 Monitored systems

We are interested in an open quantum many-body system with continuously
measured output ports, described by the Hamiltonian

$\displaystyle H(t)=H_{\rm
sys}+i\sum_{j}\sqrt{\gamma_{j}}\left[b_{j}(t)^{\dagger}c_{j}-c_{j}^{\dagger}b_{j}(t)\right]=H_{\rm
sys}+i\sum_{j}H_{j}(t)$ (S35)

with
$[b_{j}(t),b_{j^{\prime}}(t^{\prime})^{\dagger}]=\delta_{j,j^{\prime}}\delta(t-t^{\prime})$
and $b_{j}(t)|\rm vac\rangle=0~{}\forall$j. Here the quantum many-body system
takes role of the $A\cup B$ in Section II, and the bosonic modes of the
environment play the role of $M$. We are specifically interested in the
situation, where we are given a pure state of $A\cup B$, $|\phi(t)\rangle$,
propagate the the state of system and environment
$|\Psi(t)\rangle=|\phi(t)\rangle|\rm vac\rangle$ by a small time step $dt$
using the Hamiltonian (S35), to obtain the state $|\Psi(t+dt)\rangle$, and
then measure the state of the modes of the environment (separately for each
channel), to obtain a state $|\phi(t+dt)\rangle$, of $A\cup B$, which depends
on the measurement basis and measurement outcome. As detailed in the main
text, this procedure can be decomposed in a coherent component and a sequence
of stochastic components. Specifically, one defines the sequence of states
$|\phi_{j}\rangle$ as follows: Staring form the state
$|\phi_{1}\rangle=|\phi(t)\rangle$, we sequentially (for $j=1,\dots m$) do the
following:

  1. 1.

Propagate the state $|\phi_{j}\rangle|\rm vac\rangle$ of the full system
$A\cup B\cup M$ with the unitary

$\displaystyle
U_{j}=\mathcal{T}\exp\left[\int_{t}^{t+dt}dt^{\prime}H_{j}(t^{\prime})\right]\cong\exp\left[\lambda(b_{j}^{\dagger}c_{j}-c_{j}^{\dagger}b_{j})\right],$
(S36)

where the last equation holds in the limit of small $dt$ and we introduced the
normalized quantum Ito increment
$b_{j}=\frac{1}{\sqrt{dt}}\int_{t}^{t+dt}dt^{\prime}b_{j}(t)$, (satisfying
$[b_{j},b_{j}^{\dagger}]=1$ and $b_{j}|\rm vac\rangle=0$). We also introduced
$\lambda=\sqrt{dt}$ for notation simplicity in what follows. This gives the
state $|\Psi_{j}\rangle=U_{j}|\phi_{j}\rangle|\rm vac\rangle$;

  2. 2.

Measure the state of the mode $b_{j}$ in an arbitrary orthonomal basis,
specified by basis vectors $|x\rangle$. Given the measurement outcome $x$, we
obtain the state of $A\cup B$ as $|\phi_{j+1}\rangle=\langle
x|\Psi_{j}\rangle$.

Finally, we obtain the state $|\phi(t+dt)\rangle=\exp(-iH_{\rm
sys}dt)|\phi_{m+1}\rangle$. The EE changes during each of the steps
$|\phi_{j}\rangle\rightarrow|\phi_{j+1}\rangle$, as well as in the last,
coherent step. While the coherent step is fixed, the steps
$|\phi_{j}\rangle\rightarrow|\phi_{j+1}\rangle$ depend on the choice of the
measurement basis. We are thus interested in the average change rate of the EE
in these steps $\dot{\overline{E}}_{j}=\lim_{dt\rightarrow
0}\frac{1}{dt}(\overline{E(\phi_{j+1})}-E(\phi_{j}))$, where the line
indicates the average over measurement outcomes.

###  III.2 Entanglement entropy change rate

To calculate $\dot{\overline{E}}_{j}$, we apply the formalism developed in
Section II to the state $|\Psi_{j}\rangle$. We identify the measured bosonic
mode (i.e. the mode with bath annihilation and creation operators $b_{j}$ and
$b_{j}^{\dagger}$) with $M$. To proceed we assume that the jump opearator
$c_{j}$ is fully supported in either $A$ or $B$ (without loss of generality we
consider the case when it’s support is in $B$). With this we identify
$\overline{E(\phi_{j+1})}=S(A|M)=S(B|M)$, which can be obtained via (S34),
using

$\displaystyle\mu_{AM}$ $\displaystyle=\sum_{x}\langle
x|\mathrm{tr}_{B}\left(U_{j}|\phi_{j}\rangle\langle\phi_{j}|\otimes|0\rangle\langle
0|U_{j}^{\dagger}\right)|x\rangle\otimes|x\rangle\langle x|$ (S37)
$\displaystyle=\sum_{x}\mathrm{tr}_{B}\left(\langle
x|U_{j}|0\rangle|\phi_{j}\rangle\langle\phi_{j}|\langle
0|U_{j}^{\dagger}|x\rangle\right)\otimes|x\rangle\langle x|$ (S38)
$\displaystyle=\sum_{x}\mathrm{tr}_{B}\left(M_{x}\phi
M_{x}^{\dagger}\right)\otimes|x\rangle\langle x|,$ (S39)

where we defined $M_{x}=\langle x|U_{j}|0\rangle$ and
$\phi=|\phi_{j}\rangle\langle\phi_{j}|$ (suppressing the index $j$ from here
on when clear from the context), and $|0\rangle$ is the vacuum state, i.e.
$b_{j}|0\rangle=0$. We also have:

$\displaystyle\mu_{M}$ $\displaystyle=\sum_{x}\mathrm{tr}\left(M_{x}\phi
M_{x}^{\dagger}\right)|x\rangle\langle x|.$ (S40)

For the entropies we can write:

$\displaystyle S(\mu_{M})$ $\displaystyle=-\sum_{x}\mathrm{tr}\left(M_{x}\phi
M_{x}^{\dagger}\right)\log_{2}\mathrm{tr}\left(M_{x}\phi
M_{x}^{\dagger}\right),$ (S41) $\displaystyle S(\mu_{AM})$
$\displaystyle=-\sum_{x}\mathrm{tr}\left[\mathrm{tr}_{B}\left(M_{x}\phi
M_{x}^{\dagger}\right)\log_{2}\mathrm{tr}_{B}\left(M_{x}\phi
M_{x}^{\dagger}\right)\right].$ (S42)

This suggest to introduce the notation $\phi_{x}=M_{x}\phi M_{x}^{\dagger}$.
We remark that $M_{x}$ depends on $\lambda$.

In the following we calculate $\dot{\overline{E}}_{j}$ from a Taylor expansion
of $\overline{E(\phi_{j+1})}$ up to second order in $\lambda$, i.e.
$\overline{E(\phi_{j+1})}=S(A|M)|_{\lambda=0}+\lambda\partial\lambda
S(A|M)|_{\lambda=0}+\frac{1}{2}\lambda^{2}\partial^{2}_{\lambda}S(A|M)|_{\lambda=0}$.
The explicit calculation of these terms below leads to
$\partial_{\lambda}S(A|M)|_{\lambda=0}=0$, and together with
$E(\phi_{j})=S(A|M)|_{\lambda=0}$ we find

$\displaystyle\dot{\overline{E}}_{j}=\frac{1}{2}\partial^{2}_{\lambda}S(A|M)|_{\lambda=0}.$
(S43)

To evaluate the terms in this Taylor expansion we use
$\partial_{\lambda}^{k}M_{x}=\langle
x|U(b^{\dagger}c-c^{\dagger}b)^{k}|0\rangle$. Since $U|_{\lambda=0}=1$, we
have

$\displaystyle M_{x}|_{\lambda=0}=\langle x|0\rangle,$ (S44)
$\displaystyle\partial_{\lambda}M_{x}|_{\lambda=0}=\langle x|1\rangle c,$
(S45) $\displaystyle\partial_{\lambda}^{2}M_{x}|_{\lambda=0}=-\langle
x|0\rangle c^{\dagger}c+\langle x|2\rangle\sqrt{2}c^{2},$ (S46)

where $|n\rangle=\frac{{b^{\dagger}}^{n}}{\sqrt{n!}}|0\rangle$. We then find
the derivatives of $\phi_{x}$:

$\displaystyle\phi_{x}|_{\lambda=0}=|\langle x|0\rangle|^{2}\phi,$ (S47)
$\displaystyle\partial_{\lambda}\phi_{x}|_{\lambda=0}=\langle
x|1\rangle\langle 0|x\rangle c\phi+\langle x|0\rangle\langle 1|x\rangle\phi
c^{\dagger},$ (S48)
$\displaystyle\partial_{\lambda}^{2}\phi_{x}|_{\lambda=0}=-|\langle
x|0\rangle|^{2}(c^{\dagger}c\phi+\phi c^{\dagger}c)+2|\langle
x|1\rangle|^{2}c\phi c^{\dagger}+\sqrt{2}\langle x|2\rangle\langle 0|x\rangle
c^{2}\phi+\sqrt{2}\langle x|0\rangle\langle 2|x\rangle\phi{c^{\dagger}}^{2}.$
(S49)

It is useful to introduce the following:

$\displaystyle\langle x|j\rangle\rightarrow\zeta_{j}(x),$ (S50)
$\displaystyle\langle
x|0\rangle\rightarrow\zeta_{0}(x)=|\zeta_{0}(x)|e^{i\mu_{0}(x)},$ (S51)
$\displaystyle\langle
x|1\rangle\rightarrow\zeta_{1}(x)=|\zeta_{1}(x)|e^{i\mu_{1}(x)},$ (S52)
$\displaystyle\textrm{and~{}}\mu_{01}(x)=\mu_{0}(x)-\mu_{1}(x).$ (S53)

####  III.2.1 First derivative: $\partial_{\lambda}S(A|M)$

With this we can calculate the derivatives of $S(\mu_{M})$ and $S(\mu_{AM})$,
which formally are obtained as

$\displaystyle\partial_{\lambda}S(\mu_{M})$
$\displaystyle=-\sum_{x}{\mathrm{tr}\left(\partial_{\lambda}\phi_{x}\right)\log_{2}\mathrm{tr}\left(\phi_{x}\right)+\partial_{\lambda}\mathrm{tr}\left(\phi_{x}\right)/\ln
2}=-\sum_{x}\mathrm{tr}\left(\partial_{\lambda}\phi_{x}\right)\log_{2}\mathrm{tr}\left(\phi_{x}\right)$
(S54) $\displaystyle\partial_{\lambda}S(\mu_{AM})$
$\displaystyle=-\sum_{x}{\mathrm{tr}\left[{\mathrm{tr}_{B}\left(\partial_{\lambda}\phi_{x}\right)\log_{2}\mathrm{tr}_{B}\left(\phi_{x}\right)}\right]+\partial_{\lambda}\mathrm{tr}\left(\phi_{x}\right)/\ln
2}=-\sum_{x}\mathrm{tr}\left[{\mathrm{tr}_{B}\left(\partial_{\lambda}\phi_{x}\right)\log_{2}\mathrm{tr}_{B}\left(\phi_{x}\right)}\right].$
(S55)

We can then evaluate those expressions at $\lambda=0$ and obtain:

$\displaystyle\partial_{\lambda}S(\mu_{M})|_{\lambda=0}$
$\displaystyle=-\sum_{x}\mathrm{tr}\left({\zeta_{1}(x)\zeta^{*}_{0}(x)c\phi+\zeta_{0}(x)\zeta^{*}_{1}(x)\phi
c^{\dagger}}\right)\log_{2}\mathrm{tr}\left(|\zeta_{0}(x)|^{2}\phi\right)$
(S56)
$\displaystyle=-\mathrm{tr}\left(c\phi\right)\sum_{x}\zeta_{1}(x)\zeta^{*}_{0}(x)\log_{2}|\zeta_{0}(x)|^{2}+\rm{c.c.}$
(S57) $\displaystyle\partial_{\lambda}S(\mu_{AM})|_{\lambda=0}$
$\displaystyle=-\sum_{x}\mathrm{tr}\left[{\mathrm{tr}_{B}\left(\zeta_{1}(x)\zeta^{*}_{0}(x)c\phi+\zeta_{0}(x)\zeta^{*}_{1}(x)\phi
c^{\dagger}\right)\log_{2}\mathrm{tr}_{B}\left(|\zeta_{0}(x)|^{2}\phi\right)}\right]$
(S58)
$\displaystyle=\partial_{\lambda}S(\mu_{M})|_{\lambda=0}-\sum_{x}\mathrm{tr}\left[{\mathrm{tr}_{B}\left(\zeta_{1}(x)\zeta^{*}_{0}(x)c\phi+\zeta_{0}(x)\zeta^{*}_{1}(x)\phi
c^{\dagger}\right)\log\mathrm{tr}_{B}\left(\phi\right)}\right]$ (S59)
$\displaystyle=\partial_{\lambda}S(\mu_{M})|_{\lambda=0}.$ (S60)

Here we used the property $\sum_{x}\langle x|j\rangle\langle
k|x\rangle=\langle k|j\rangle=\delta_{k,j}$. In summary we have

$\displaystyle\partial_{\lambda}S(A|M)|_{\lambda=0}=0,$ (S61)

as stated above.

####  III.2.2 Second derivative: $\partial^{2}_{\lambda}S(A|M)$

We now proceed to calculate

$\displaystyle\partial_{\lambda}^{2}S(\mu_{M})$
$\displaystyle=-\sum_{x}\mathrm{tr}\left(\partial_{\lambda}^{2}\phi_{x}\right)\log_{2}\mathrm{tr}\left(\phi_{x}\right)+\frac{\left[\mathrm{tr}\left(\partial_{\lambda}\phi_{x}\right)\right]^{2}}{\mathrm{tr}\left(\phi_{x}\right)\ln
2},$ (S62) $\displaystyle\partial_{\lambda}^{2}S(\mu_{AM})$
$\displaystyle=-\sum_{x}\mathrm{tr}\left[{\mathrm{tr}_{B}\left(\partial_{\lambda}^{2}\phi_{x}\right)\log_{2}\mathrm{tr}_{B}\left(\phi_{x}\right)}+\mathrm{tr}_{B}\left(\partial_{\lambda}\phi_{x}\right)\partial_{\lambda}\log_{2}\mathrm{tr}_{B}\left(\phi_{x}\right)\right].$
(S63)

Again, we are interested in evaluating this at $\lambda=0$. More precisely, we
are interested in the quantity

$\displaystyle\dot{\overline{E}}_{j}\equiv\frac{1}{2}\partial_{\lambda}^{2}S(A|M)|_{\lambda=0}=\frac{1}{2}\lim_{\lambda\rightarrow
0}$
$\displaystyle\left[\partial_{\lambda}^{2}S(\mu_{AM})-\partial_{\lambda}^{2}S(\mu_{M})\right]$
(S64) $\displaystyle=\frac{1}{2\ln 2}\lim_{\lambda\rightarrow 0}\bigg{\\{}$
$\displaystyle\sum_{x}\mathrm{tr}\left(\partial_{\lambda}^{2}\phi_{x}\right)\ln\mathrm{tr}\left(\phi_{x}\right)-\mathrm{tr}\left[{\mathrm{tr}_{B}\left(\partial_{\lambda}^{2}\phi_{x}\right)\ln\mathrm{tr}_{B}\left(\phi_{x}\right)}\right]$
(S65)
$\displaystyle+\frac{\left[\mathrm{tr}\left(\partial_{\lambda}\phi_{x}\right)\right]^{2}}{\mathrm{tr}\left(\phi_{x}\right)}-\mathrm{tr}\left[{\mathrm{tr}_{B}\left(\partial_{\lambda}\phi_{x}\right)R^{-1}_{\mathrm{tr}_{B}\left(\phi_{x}\right)}(\mathrm{tr}_{B}\left(\partial_{\lambda}\phi_{x}\right))}\right]\bigg{\\}}.$
(S66)

Here we introduced the logarithmic derivative, defining

$\displaystyle R^{-1}_{\rho}(X)=\sum_{j,k}\frac{\ln r_{j}-\ln
r_{j}}{r_{j}-r_{k}}\langle r_{j}|X|r_{k}\rangle|r_{j}\rangle\langle r_{k}|$
(S67)

with $\rho|r_{k}\rangle=r_{k}|r_{k}\rangle$. The above expression for the EAEE
time-derivative $\dot{\overline{E}}_{j}$ is valid for an arbitrary measurement
of the environement mode of the channel $j$. Note that evaluating the limit
$\lambda\rightarrow 0$ requires special attention if the measurement basis
contains a state $|x\rangle$ such that $\langle x|0\rangle=0$, which we refer
to as vacuum-resolving measurements, i.e. the measurements that exclude the
vacuum state of the bath. If there is no such state, we name the measurement
as non-vacuum-resolving. We consider these cases separately in the following.

###  III.3 Non-vacuum-resolving measurement

Let us consider the case where $\zeta_{0}(x)\neq 0$ for all $x$. And let us
look now at each term in the expression for $\dot{\overline{E}}_{j}$, i.e.
(S65) and (S66), individually. The first two terms in (S65) can be written as

$\displaystyle\sum_{x}$
$\displaystyle\mathrm{tr}\left(\partial_{\lambda}^{2}\phi_{x}\right)\ln\mathrm{tr}\left(\phi_{x}\right)-\mathrm{tr}\left(\mathrm{tr}_{B}\left(\partial_{\lambda}^{2}\phi_{x}\right)\ln\mathrm{tr}_{B}\left(\phi_{x}\right)\right)$
(S68) $\displaystyle=\sum_{x}$
$\displaystyle\left(|\zeta_{1}(x)|^{2}-|\zeta_{0}(x)|^{2}\right)\mathrm{tr}\left(2c\phi
c^{\dagger}\right)\ln{\left\\{|\zeta_{0}(x)|^{2}+\mathrm{tr}\left[{O(x,\lambda)}\right]\right\\}}$
(S69) $\displaystyle-\sum_{x}$
$\displaystyle\mathrm{tr}\left(\left[|\zeta_{1}(x)|^{2}\mathrm{tr}_{B}\left(2c\phi
c^{\dagger}\right)-|\zeta_{0}(x)|^{2}\mathrm{tr}_{B}\left(c^{\dagger}c\phi+\phi
c^{\dagger}c\right)\right]\ln{\left\\{|\zeta_{0}(x)|^{2}\mathrm{tr}_{B}\left(\phi\right)+\mathrm{tr}_{\mathrm{B}}\left[{O(x,\lambda)}\right]\right\\}}\right),$
(S70)

where
$O(x,\lambda)=\lambda\zeta_{1}(x)\zeta_{0}^{*}(x)c\phi+\textrm{h.c.}+\mathcal{O}(\lambda^{2})$,
which is not relevant in the limit $\lambda\rightarrow 0$ in the non-vacuum-
resolving case. In this limit we get

$\displaystyle\lim_{\lambda\rightarrow 0}$
$\displaystyle\sum_{x}\mathrm{tr}\left(\partial_{\lambda}^{2}\phi_{x}\right)\ln\mathrm{tr}\left(\phi_{x}\right)-\mathrm{tr}\left(\mathrm{tr}_{B}\left(\partial_{\lambda}^{2}\phi_{x}\right)\ln\mathrm{tr}_{B}\left(\phi_{x}\right)\right)=\sum_{x}\left(|\zeta_{1}(x)|^{2}-|\zeta_{0}(x)|^{2}\right)\mathrm{tr}\left(2c\phi
c^{\dagger}\right)\ln{\left[|\zeta_{0}(x)|^{2}\right]}$ (S71)
$\displaystyle-\sum_{x}\mathrm{tr}\left\\{\left[|\zeta_{1}(x)|^{2}\mathrm{tr}_{B}\left(2c\phi
c^{\dagger}\right)-|\zeta_{0}(x)|^{2}\mathrm{tr}_{B}\left(c^{\dagger}c\phi+\phi
c^{\dagger}c\right)\right]\ln{\left[|\zeta_{0}(x)|^{2}\mathrm{tr}_{B}\left(\phi\right)\right]}\right\\}$
(S72) $\displaystyle=$
$\displaystyle-\sum_{x}\mathrm{tr}\left\\{\left(|\zeta_{1}(x)|^{2}-|\zeta_{0}(x)|^{2}\right)\mathrm{tr}_{B}\left(2c\phi
c^{\dagger}\right)\ln{\left[\mathrm{tr}_{B}\left(\phi\right)\right]}\right\\}=0.$
(S73)

In the last line we first made use of the fact that $c$ acts non-trivially
only in $B$, which allows us to use the cyclic property of the partial trace,
and then we used $\sum_{x}\langle x|j\rangle\langle k|x\rangle=\langle
k|j\rangle=\delta_{k,j}$ [we remind the reader of the short-hand notations
from Eqs. (S51) and (S52)].

The first term in (S66) can be expressed as follows:

$\displaystyle\sum_{x}\frac{1}{2\ln 2}\lim_{\lambda\rightarrow
0}\frac{\left[\mathrm{tr}\left(\partial_{\lambda}\phi_{x}\right)\right]^{2}}{\mathrm{tr}\left(\phi_{x}\right)}=\sum_{x}\frac{|\zeta_{1}(x)|^{2}}{2\ln
2}\bigg{|}e^{-i\mu_{01}(x)}\mathrm{tr}\left(c\phi\right)+e^{+i\mu_{01}(x)}\mathrm{tr}\left(\phi
c^{\dagger}\right)\bigg{|}^{2}.$ (S74)

For the last term in (S66) we first note that

$\displaystyle\lim_{\lambda\rightarrow
0}R^{-1}_{\mathrm{tr}_{B}\left(\phi_{x}\right)}(X)|_{\lambda=0}=R^{-1}_{|\zeta_{0}(x)|^{2}\mathrm{tr}_{B}\left(\phi\right)}(X)=\frac{1}{|\zeta_{0}(x)|^{2}}R^{-1}_{\mathrm{tr}_{B}\left(\phi\right)}(X)=\frac{1}{|\zeta_{0}(x)|^{2}}\sum_{j,k}\frac{\ln\xi_{j}-\ln\xi_{k}}{\xi_{j}-\xi_{k}}\langle\xi_{j}|X|\xi_{k}\rangle|\xi_{j}\rangle\langle\xi_{k}|$
(S75)

with $\mathrm{tr}_{B}\left(\phi\right)|\xi_{k}\rangle=\xi_{k}|\xi_{k}\rangle$,
i.e. $\xi_{k}$’s and $|\xi_{k}\rangle$’s are the eigenvalues and the
eigenvectors of the reduced state of subsystem $A$, correspondingly. With this
we write the last term as

$\displaystyle-\frac{1}{2\ln
2}\sum_{x}|\zeta_{1}(x)|^{2}\sum_{j,k}\frac{\ln\xi_{j}-\ln\xi_{k}}{\xi_{j}-\xi_{k}}\bigg{|}e^{-i\mu_{01}(x)}\langle\xi_{j}|\mathrm{tr}_{B}\left(c\phi\right)|\xi_{k}\rangle+e^{+i\mu_{01}(x)}\langle\xi_{j}|\mathrm{tr}_{B}\left(\phi
c^{\dagger}\right)|\xi_{k}\rangle\bigg{|}^{2}.$ (S76)

In summary, we obtain for non-vacuum-resolving measurement the following
expression for the time-derivative of the EAEE:

$\displaystyle\dot{\overline{E}}_{\mathrm{NVR}}=\frac{1}{2\ln
2}\sum_{x}|\zeta_{1}(x)|^{2}\Bigg{[}$
$\displaystyle\bigg{|}e^{-i\mu_{01}(x)}\mathrm{tr}\left(c\phi\right)+e^{+i\mu_{01}(x)}\mathrm{tr}\left(\phi
c^{\dagger}\right)\bigg{|}^{2}-$ (S78)
$\displaystyle-\sum_{j,k}\frac{\ln\xi_{j}-\ln\xi_{k}}{\xi_{j}-\xi_{k}}\bigg{|}e^{-i\mu_{01}(x)}\langle\xi_{j}|\mathrm{tr}_{B}\left(c\phi\right)|\xi_{k}\rangle+e^{+i\mu_{01}(x)}\langle\xi_{j}|\mathrm{tr}_{B}\left(\phi
c^{\dagger}\right)|\xi_{k}\rangle\bigg{|}^{2}\Bigg{]}.$ (S79)

Both terms exhibit a dependence on the measurement basis, specifically on its
“phase” $\mu_{01}(x)$.

###  III.4 Vacuum-resolving measurement

In the vacuum-resolving case we need to consider the terms in the sum over $x$
where $\zeta_{0}(x)=0$ separately. Thus we split the sum in
$\dot{\overline{E}}_{j}$ into two parts, one part with $\zeta_{0}(x)\neq 0$
and another part with $\zeta_{0}(x)=0$:

$\displaystyle\dot{\overline{E}}_{j}=$ $\displaystyle\frac{1}{2\ln
2}\lim_{\lambda\rightarrow 0}\bigg{(}\sum_{x|\zeta_{0}(x)\neq
0}\mathrm{tr}\left(\partial_{\lambda}^{2}\phi_{x}\right)\ln\mathrm{tr}\left(\phi_{x}\right)-\mathrm{tr}\left[{\mathrm{tr}_{B}\left(\partial_{\lambda}^{2}\phi_{x}\right)\ln\mathrm{tr}_{B}\left(\phi_{x}\right)}\right]$
(S80)
$\displaystyle+\frac{\left[\mathrm{tr}\left(\partial_{\lambda}\phi_{x}\right)\right]^{2}}{\mathrm{tr}\left(\phi_{x}\right)}-\mathrm{tr}\left\\{{\mathrm{tr}_{B}\left(\partial_{\lambda}\phi_{x}\right)R^{-1}_{\mathrm{tr}_{B}\left(\phi_{x}\right)}\left[\mathrm{tr}_{B}\left(\partial_{\lambda}\phi_{x}\right)\right]}\right\\}\bigg{)}$
(S81) $\displaystyle+$ $\displaystyle\frac{1}{2\ln 2}\lim_{\lambda\rightarrow
0}\bigg{(}\sum_{x|\zeta_{0}(x)=0}\mathrm{tr}\left(\partial_{\lambda}^{2}\phi_{x}\right)\ln\mathrm{tr}\left(\phi_{x}\right)-\mathrm{tr}\left[{\mathrm{tr}_{B}\left(\partial_{\lambda}^{2}\phi_{x}\right)\ln\mathrm{tr}_{B}\left(\phi_{x}\right)}\right]$
(S82)
$\displaystyle+\frac{\left[\mathrm{tr}\left(\partial_{\lambda}\phi_{x}\right)\right]^{2}}{\mathrm{tr}\left(\phi_{x}\right)}-\mathrm{tr}\left\\{{\mathrm{tr}_{B}\left(\partial_{\lambda}\phi_{x}\right)R^{-1}_{\mathrm{tr}_{B}\left(\phi_{x}\right)}\left[\mathrm{tr}_{B}\left(\partial_{\lambda}\phi_{x}\right)\right]}\right\\}\bigg{)}.$
(S83)

For the first line (S80) we can use the above result to find

$\displaystyle\lim_{\lambda\rightarrow 0}\frac{1}{2\ln
2}\sum_{x|\zeta_{0}(x)\neq
0}\mathrm{tr}\left(\partial_{\lambda}^{2}\phi_{x}\right)\ln\mathrm{tr}\left(\phi_{x}\right)-\mathrm{tr}\left[{\mathrm{tr}_{B}\left(\partial_{\lambda}^{2}\phi_{x}\right)\ln\mathrm{tr}_{B}\left(\phi_{x}\right)}\right]$
(S84) $\displaystyle=-\frac{1}{2\ln 2}\sum_{x|\zeta_{0}(x)\neq
0}\mathrm{tr}\left\\{\left(|\zeta_{1}(x)|^{2}-|\zeta_{0}(x)|^{2}\right)\mathrm{tr}_{B}\left(2c\phi
c^{\dagger}\right)\ln{\left[\mathrm{tr}_{B}\left(\phi\right)\right]}\right\\}.$
(S85)

The contribution of the next line (S81) can be rewritten as

$\displaystyle\frac{1}{2\ln 2}\sum_{x|\zeta_{0}(x)\neq
0}|\zeta_{1}(x)|^{2}\Bigg{[}$
$\displaystyle\bigg{|}e^{-i\mu_{01}(x)}\mathrm{tr}\left(c\phi\right)+e^{+i\mu_{01}(x)}\mathrm{tr}\left(\phi
c^{\dagger}\right)\bigg{|}^{2}$ (S86)
$\displaystyle-\sum_{j,k}\frac{\ln\xi_{j}-\ln\xi_{k}}{\xi_{j}-\xi_{k}}\bigg{|}e^{-i\mu_{01}(x)}\langle\xi_{j}|\mathrm{tr}_{B}\left(c\phi\right)|\xi_{k}\rangle+e^{+i\mu_{01}(x)}\langle\xi_{j}|\mathrm{tr}_{B}\left(\phi
c^{\dagger}\right)|\xi_{k}\rangle\bigg{|}^{2}\Bigg{]}.$ (S87)

The terms in (S82) are non-trivial. To calculate them, first note that for
$\zeta_{0}(x_{0})=0$, in $\phi_{x_{0}}$ both the zeroth and first order
contribution in $\lambda$ vanish, i.e.
$\phi_{x_{0}}=|\zeta_{1}(x_{0})|^{2}2c\phi
c^{\dagger}\lambda^{2}/2+\mathcal{O}(\lambda^{3})$. We consider the case
$\zeta_{1}(x_{0})\neq 0$ (the case $\zeta_{1}(x_{0})=0$ gives a vanishing
result):

$\displaystyle\lim_{\lambda\rightarrow 0}$
$\displaystyle\bigg{[}\mathrm{tr}\left(\partial_{\lambda}^{2}\phi_{x_{0}}\right)\ln\mathrm{tr}\left(\phi_{x_{0}}\right)-\mathrm{tr}\left(\mathrm{tr}_{B}\left(\partial_{\lambda}^{2}\phi_{x_{0}}\right)\ln\mathrm{tr}_{B}\left(\phi_{x_{0}}\right)\right)\bigg{]}$
(S88) $\displaystyle=$ $\displaystyle\lim_{\lambda\rightarrow
0}|\zeta_{1}(x_{0})|^{2}\bigg{(}\mathrm{tr}\left(2c\phi
c^{\dagger}\right)\ln{\left[\lambda^{2}|\zeta_{1}(x_{0})|^{2}\mathrm{tr}\left(c\phi
c^{\dagger}\right)\right]}-\mathrm{tr}\left\\{\mathrm{tr}_{B}\left(2c\phi
c^{\dagger}\right)\ln{\left[\lambda^{2}|\zeta_{1}(x_{0})|^{2}\mathrm{tr}_{B}\left(c\phi
c^{\dagger}\right)\right]}\right\\}\bigg{)}$ (S89) $\displaystyle=$
$\displaystyle|\zeta_{1}(x_{0})|^{2}\bigg{(}\mathrm{tr}\left(2c\phi
c^{\dagger}\right)\ln{\left[\mathrm{tr}\left(c\phi
c^{\dagger}\right)\right]}-\mathrm{tr}\left\\{\mathrm{tr}_{B}\left(2c\phi
c^{\dagger}\right)\ln{\left[\mathrm{tr}_{B}\left(c\phi
c^{\dagger}\right)\right]}\right\\}\bigg{)}.$ (S90)

For the two terms in (S83) note that $\zeta_{0}(x_{0})=0$ implies
$\lim_{\lambda\rightarrow 0}\mathrm{tr}\left(\phi_{x_{0}}\right)=0$ and also
$\lim_{\lambda\rightarrow
0}\mathrm{tr}\left(\partial_{\lambda}\phi_{x_{0}}\right)=0$. From that we find
$\lim_{\lambda\rightarrow
0}\frac{\mathrm{tr}\left(\partial_{\lambda}\phi_{x_{0}}\right)^{2}}{\mathrm{tr}\left(\phi_{x_{0}}\right)}=2\mathrm{tr}\left(\partial_{\lambda}^{2}\phi_{x_{0}}\right)$.
This follows simply by
$\phi_{x}=\phi_{x}|_{\lambda=0}+\lambda\partial_{\lambda}\phi_{x}|_{\lambda=0}+\frac{1}{2}\lambda^{2}\partial^{2}_{\lambda}\phi_{x}|_{\lambda=0}+\mathcal{O}(\lambda^{3})$,
and
$\partial_{\lambda}\phi_{x}=\partial_{\lambda}\phi_{x}|_{\lambda=0}+\lambda\partial^{2}_{\lambda}\phi_{x}|_{\lambda=0}+\mathcal{O}(\lambda^{2})$,
and the L’Hospital’s rule. A similar issue arises in the second term of (S83),
since also $\lim_{\lambda\rightarrow
0}\mathrm{tr}_{B}\left(\phi_{x_{0}}\right)=0$ (as an operator). Then we have

$\displaystyle\lim_{\lambda\rightarrow
0}\mathrm{tr}\bigg{\\{}{\mathrm{tr}_{B}\left(\partial_{\lambda}\phi_{x_{0}}\right)R^{-1}_{\mathrm{tr}_{B}\left(\phi_{x_{0}}\right)}\left[\mathrm{tr}_{B}\left(\partial_{\lambda}\phi_{x_{0}}\right)\right]}\bigg{\\}}=2\mathrm{tr}\bigg{\\{}{\mathrm{tr}_{B}\left(\partial^{2}_{\lambda}\phi_{x_{0}}\right)R^{-1}_{\mathrm{tr}_{B}\left(\partial^{2}_{\lambda}\phi_{x_{0}}\right)}\left[\mathrm{tr}_{B}\left(\partial^{2}_{\lambda}\phi_{x_{0}}\right)\right]}\bigg{\\}}\bigg{|}_{\lambda=0}=2\mathrm{tr}\left(\partial^{2}_{\lambda}\phi_{x_{0}}\right)|_{\lambda=0}.$
(S91)

This can be proven as follows. Consider the expression with the logarithmic
derivative:

$\displaystyle\lim_{\lambda\rightarrow
0}\mathrm{tr}\bigg{\\{}{\mathrm{tr}_{B}\left(\partial_{\lambda}\phi_{x_{0}}\right)R^{-1}_{\mathrm{tr}_{B}\left(\phi_{x_{0}}\right)}\left[\mathrm{tr}_{B}\left(\partial_{\lambda}\phi_{x_{0}}\right)\right]}\bigg{\\}},$
(S92)

where $\phi_{x_{0}}=\lambda^{2}/2|\langle x|1\rangle|^{2}2c\phi c^{\dagger}$
and $\partial_{\lambda}\phi_{x_{0}}=\lambda|\langle x|1\rangle|^{2}2c\phi
c^{\dagger}$ for $|\langle x|0\rangle|^{2}=0$ (see definitions above). We can
set the matrix $\mathrm{tr}_{B}\left(2c\phi
c^{\dagger}\right)=\sum_{i}r_{i}|r_{i}\rangle\langle r_{i}|$ diagonalized in
basis $\\{|r_{i}\rangle\\}$, such that

$\displaystyle\mathrm{tr}_{B}\left(\phi_{x_{0}}\right)=\frac{\lambda^{2}}{2}|\langle
x|1\rangle|^{2}\sum_{i}r_{i}|r_{i}\rangle\langle r_{i}|,$ (S93)
$\displaystyle\mathrm{tr}_{B}\left(\partial_{\lambda}\phi_{x_{0}}\right)=\lambda|\langle
x|1\rangle|^{2}\sum_{i}r_{i}|r_{i}\rangle\langle r_{i}|.$ (S94)

Thus the limit simplifies as follows:

$\displaystyle\lim_{\lambda\rightarrow 0}$
$\displaystyle\mathrm{tr}\bigg{\\{}{\mathrm{tr}_{B}\left(\partial_{\lambda}\phi_{x_{0}}\right)R^{-1}_{\mathrm{tr}_{B}\left(\phi_{x_{0}}\right)}\left[\mathrm{tr}_{B}\left(\partial_{\lambda}\phi_{x_{0}}\right)\right]}\bigg{\\}}=$
(S95) $\displaystyle=$ $\displaystyle\lim_{\lambda\rightarrow
0}\mathrm{tr}\Bigg{\\{}{|\langle
x|1\rangle|^{2}\lambda\left(\sum_{i}r_{i}|r_{i}\rangle\langle
r_{i}|\right)\left[\sum_{j,k}\frac{\ln r_{j}-\ln
r_{k}}{\lambda^{2}/2(r_{j}-r_{k})}\langle
r_{j}|\lambda\left(\sum_{i}r_{i}|r_{i}\rangle\langle
r_{i}|\right)|r_{k}\rangle|r_{j}\rangle\langle r_{k}|\right]}\Bigg{\\}}=$
(S96) $\displaystyle=$ $\displaystyle\lim_{\lambda\rightarrow
0}\mathrm{tr}\left[{|\langle
x|1\rangle|^{2}\lambda\left(\sum_{i}r_{i}|r_{i}\rangle\langle
r_{i}|\right)\frac{2\lambda}{\lambda^{2}}\left(\sum_{i}|r_{i}\rangle\langle
r_{i}|\right)}\right]=2|\langle x|1\rangle|^{2}\sum_{i}r_{i}\equiv
2\mathrm{tr}\left(\partial^{2}_{\lambda}\phi_{x0}\right)|_{\lambda=0},$ (S97)

where again ${r_{i}}$ is the set of eigenvalues of matrix $2c\phi c^{\dagger}$
and we used the L’Hospital’s rule for the expression:

$\displaystyle\lim_{r_{j}\rightarrow
r_{i}}\frac{\ln{r_{j}}-\ln{r_{i}}}{r_{j}-r_{i}}=\frac{1}{r_{i}}.$ (S98)

Thus the terms in (S83) exactly cancel. In summary, we obtain for a vacuum-
resolving case the following expression:

$\displaystyle\dot{\overline{E}}_{\mathrm{VR}}=$ $\displaystyle-\frac{1}{2\ln
2}\sum_{x|\zeta_{0}(x)\neq
0}\left(|\zeta_{1}(x)|^{2}-|\zeta_{0}(x)|^{2}\right)\mathrm{tr}\left\\{\mathrm{tr}_{B}\left(2c\phi
c^{\dagger}\right)\ln{\left[\mathrm{tr}_{B}\left(\phi\right)\right]}\right\\}$
(S99) $\displaystyle+\frac{1}{2\ln 2}\sum_{x|\zeta_{0}(x)\neq
0}|\zeta_{1}(x)|^{2}\Bigg{[}\bigg{|}e^{-i\mu_{01}(x)}\mathrm{tr}\left(c\phi\right)+e^{+i\mu_{01}(x)}\mathrm{tr}\left(\phi
c^{\dagger}\right)\bigg{|}^{2}$ (S100)
$\displaystyle\qquad\qquad\qquad-\sum_{j,k}\frac{\ln\xi_{j}-\ln\xi_{k}}{\xi_{j}-\xi_{k}}\bigg{|}e^{-i\mu_{01}(x)}\langle\xi_{j}|\mathrm{tr}_{B}\left(c\phi\right)|\xi_{k}\rangle+e^{+i\mu_{01}(x)}\langle\xi_{j}|\mathrm{tr}_{B}\left(\phi
c^{\dagger}\right)|\xi_{k}\rangle\bigg{|}^{2}\Bigg{]}$ (S101)
$\displaystyle+\frac{1}{2\ln
2}\sum_{x|\zeta_{0}(x)=0}|\zeta_{1}(x)|^{2}\bigg{(}\mathrm{tr}\left(2c\phi
c^{\dagger}\right)\ln{\left[\mathrm{tr}\left(c\phi
c^{\dagger}\right)\right]}-\mathrm{tr}\left\\{\mathrm{tr}_{B}\left(2c\phi
c^{\dagger}\right)\ln{\left[\mathrm{tr}_{B}\left(c\phi
c^{\dagger}\right)\right]}\right\\}\bigg{)}.$ (S102)

We can further simplify this expression by noting that
$\sum_{x|\zeta_{0}(x)\neq 0}\dots=\sum_{x}\dots-\sum_{x|\zeta_{0}(x)=0}\dots$.
The first line above then simplifies to

$\displaystyle+\frac{1}{2\ln
2}\sum_{x|\zeta_{0}(x)=0}|\zeta_{1}(x)|^{2}\mathrm{tr}\left\\{\mathrm{tr}_{B}\left(2c\phi
c^{\dagger}\right)\ln{\left[\mathrm{tr}_{B}\left(\phi\right)\right]}\right\\}.$
(S103)

The expression (S103) seems to be problematic in the case when
$\mathrm{tr}_{B}\left(\phi\right)$ has zero eigenvalues. However one can
notice that, since the jump operator $c$ is local in the B-partition of the
system, the matrices $\mathrm{tr}_{B}\left(\phi\right)$ and
$\mathrm{tr}_{B}\left(2c\phi c^{\dagger}\right)$ share the same eigenspace
with zero eigenvalue, thus avoiding divergences.

Thus we finally obtain for the vacuum-resolving case:

$\displaystyle\dot{\overline{E}}_{\mathrm{VR}}=$ $\displaystyle\frac{1}{2\ln
2}\sum_{x|\zeta_{0}(x)\neq
0}|\zeta_{1}(x)|^{2}\bigg{(}\bigg{|}e^{-i\mu_{01}(x)}\mathrm{tr}\left(c\phi\right)+e^{+i\mu_{01}(x)}\mathrm{tr}\left(\phi
c^{\dagger}\right)\bigg{|}^{2}$ (S104)
$\displaystyle\qquad\qquad\qquad-\sum_{j,k}\frac{\ln\xi_{j}-\ln\xi_{k}}{\xi_{j}-\xi_{k}}\bigg{|}e^{-i\mu_{01}(x)}\langle\xi_{j}|\mathrm{tr}_{B}\left(c\phi\right)|\xi_{k}\rangle+e^{+i\mu_{01}(x)}\langle\xi_{j}|\mathrm{tr}_{B}\left(\phi
c^{\dagger}\right)|\xi_{k}\rangle\bigg{|}^{2}\bigg{)}$ (S105)
$\displaystyle+\frac{1}{2\ln
2}\sum_{x|\zeta_{0}(x)=0}|\zeta_{1}(x)|^{2}\Bigg{[}\mathrm{tr}\left(2c\phi
c^{\dagger}\right)\ln{\left[\mathrm{tr}\left(c\phi
c^{\dagger}\right)\right]}+\mathrm{tr}\left(\mathrm{tr}_{B}\left(2c\phi
c^{\dagger}\right)\left\\{\ln{\left[\mathrm{tr}_{B}\left(\phi\right)\right]}-\ln{\left[\mathrm{tr}_{B}\left(2c\phi
c^{\dagger}\right)\right]}\right\\}\right)\Bigg{]}.$ (S106)

###  III.5 Number measurement

Let us consider the special case of a number measurement, i.e. the measurement
basis $\\{|0\rangle,|1\rangle,|2\rangle,\dots\\}$. One can immediately see
that this corresponds to a vacuum-resolving measurement. We note that for this
measurement, when $|\zeta_{0}(x)|^{2}=0$, it immediately implies that
$|\zeta_{1}(x)|^{2}=1$, and vice versa. So that the first terms for
$\dot{\overline{E}}$, Eq. (S104) and Eq. (S105), vanish, as $\zeta_{1}(x)=0$,
thus the phase-sensitive part of the expression vanishes, as expected for a
number measurement setup. The last term, Eq. (S106), stays and gives the
expression for a number measurement case:

$\displaystyle\dot{\overline{E}}_{\mathrm{num}}=\mathrm{tr}\left(c\phi
c^{\dagger}\right)\log_{2}{\left[\mathrm{tr}\left(c\phi
c^{\dagger}\right)\right]}+\mathrm{tr}\bigg{(}\mathrm{tr}_{B}\left(c\phi
c^{\dagger}\right)\left\\{\log_{2}{\left[\mathrm{tr}_{B}\left(\phi\right)\right]}-\log_{2}{\left[\mathrm{tr}_{B}\left(c\phi
c^{\dagger}\right)\right]}\right\\}\bigg{)}.$ (S107)

Note that this expression only depends on the state itself and the jump
operator.

###  III.6 Homodyne measurement

Let us now consider a continuous measurement basis in case of an ideal
homodyne measurement, in which case the measurement basis is the eigenbasis of
$b_{j}(t)e^{i\varphi}+b_{j}^{\dagger}(t)e^{-i\varphi}$. Since this operator
has a continuous spectrum, we replace the summation over $x$ with an
integration, $\sum_{x}\rightarrow\int\,dx$, and the corresponding wave
functions of continuous measurement are the wave functions of a harmonic
oscillator:

$\displaystyle\zeta_{0}(x)$ $\displaystyle=\langle
x|0\rangle=\frac{1}{\pi^{1/4}}e^{-x^{2}/2},$ (S108)
$\displaystyle\zeta_{1}(x)$ $\displaystyle=\langle
x|1\rangle=\frac{\sqrt{2}xe^{-i\varphi}}{\pi^{1/4}}e^{-x^{2}/2},$ (S109)
$\displaystyle\mu_{01}(x)$ $\displaystyle=\mu_{01}=\varphi.$ (S110)

Clearly, this measurement falls into the non-vacuum-resolving category (see
Subsection III.3). Thus, the expression for the time derivative of
entanglement entropy is obtained as [Eqs. (S78-S79)]:

$\displaystyle\dot{\overline{E}}_{\mathrm{hom}}=\frac{1}{2\ln
2}\left(\int\frac{2x^{2}}{\sqrt{\pi}}e^{-x^{2}}\,dx\right)\Bigg{[}$
$\displaystyle\bigg{|}e^{-i\varphi}\mathrm{tr}\left(c\phi\right)+e^{+i\varphi}\mathrm{tr}\left(\phi
c^{\dagger}\right)\bigg{|}^{2}-$ (S111)
$\displaystyle-\sum_{j,k}\frac{\ln\xi_{j}-\ln\xi_{k}}{\xi_{j}-\xi_{k}}\bigg{|}e^{-i\varphi}\langle\xi_{j}|\mathrm{tr}_{B}\left(c\phi\right)|\xi_{k}\rangle+e^{+i\varphi}\langle\xi_{j}|\mathrm{tr}_{B}\left(\phi
c^{\dagger}\right)|\xi_{k}\rangle\bigg{|}^{2}\Bigg{]},$ (S112)

where the integral in brackets is equal to 1, so that we finally obtain:

$\displaystyle\dot{\overline{E}}_{\mathrm{hom}}=\frac{1}{2\ln 2}\Bigg{[}$
$\displaystyle\bigg{|}e^{-i\varphi}\mathrm{tr}\left(c\phi\right)+e^{+i\varphi}\mathrm{tr}\left(\phi
c^{\dagger}\right)\bigg{|}^{2}-$ (S113)
$\displaystyle-\sum_{j,k}\frac{\ln\xi_{j}-\ln\xi_{k}}{\xi_{j}-\xi_{k}}\bigg{|}e^{-i\varphi}\langle\xi_{j}|\mathrm{tr}_{B}\left(c\phi\right)|\xi_{k}\rangle+e^{+i\varphi}\langle\xi_{j}|\mathrm{tr}_{B}\left(\phi
c^{\dagger}\right)|\xi_{k}\rangle\bigg{|}^{2}\Bigg{]}.$ (S114)

Note that $\dot{\overline{E}}_{\mathrm{hom}}$ depends on $\varphi$. In the
following we write this dependence explicitly as
$\dot{\overline{E}}_{\mathrm{hom}}^{\varphi}$.

###  III.7 Optimal measurements

####  III.7.1 Non-vacuum-resolving case

Now let us consider an arbitrary non-vacuum-resolving measurement $\nu$. Note
that $\dot{\overline{E}}_{\nu}\leq 0$, if $c$ has support only in $A$ or $B$,
i.e. if $c$ is local. Let us define new variables $r$ and $\beta$:

$\displaystyle\sum_{x}|\zeta_{1}(x)|^{2}e^{-2i\mu_{01}(x)}=r^{2}e^{-2i\beta}$
(S115)

with a property $\sum_{x}|\zeta_{1}(x)|^{2}=1$. Here $0\leq r\leq 1$, and
$0\leq\beta<2\pi$. Note that $r$ and $\beta$ are determined by $\nu$. If $\nu$
is the the ideal homodyne measurement, then we get $r=1$ and $\beta=\varphi$.
With these definitions we have

$\displaystyle\dot{\overline{E}}_{\nu}$ $\displaystyle=\frac{1}{2\ln
2}\bigg{[}\left(r^{2}e^{-2i\beta}a_{0}^{2}+2|a_{0}|^{2}+r^{2}e^{2i\beta}{a_{0}^{\ast}}^{2}\right)-\sum_{i}\left(r^{2}e^{-2i\beta}a_{i}^{2}+2|a_{i}|^{2}+r^{2}e^{2i\beta}{a_{i}^{\ast}}^{2}\right)\bigg{]},$
(S116)

where $a_{0}=\mathrm{tr}\left(c\phi\right)$ and
$a_{i}=a_{jk}=\langle\xi_{j}|\mathrm{tr}_{B}\left(c\phi\right)|\xi_{k}\rangle$,
where the index $i$ refers to a double-index $j,k$. Simplifying the latter
expression, one can obtain:

$\displaystyle\dot{\overline{E}}_{\nu}$
$\displaystyle=r^{2}\dot{\overline{E}}_{\rm
hom}^{(\varphi=\beta)}+\frac{1-r^{2}}{\ln
2}\left(|a_{0}|^{2}-\sum_{i}|a_{i}|^{2}\right),$ (S117)

The expression above is a monotonic function of $r$ in $0\leq r\leq 1$. Thus
its minimum is obtained at one of the boundary points, either at $r=0$ or at
$r=1$. The options for the minimum value are thus:

  1. (i)

$\dot{\overline{E}}_{\nu}=(|a_{0}|^{2}-\sum_{i}|a_{i}|^{2})/{\ln 2}$, or

  2. (ii)

$\dot{\overline{E}}_{\nu}=\dot{\overline{E}}_{\rm hom}^{(\varphi=\beta)}$.

If $(|a_{0}|^{2}-\sum_{i}|a_{i}|^{2})/{\ln 2}>\dot{\overline{E}}_{\rm
hom}^{(\varphi=\beta)}$, then then homodyne measurement with $\varphi=\beta$
achieves a better $\dot{\overline{E}}$ then the measurement $\nu$. On the
other hand, if $(|a_{0}|^{2}-\sum_{i}|a_{i}|^{2})/{\ln
2}<\dot{\overline{E}}_{\rm hom}^{(\varphi=\beta)}$, then the homodyne
measurement with $\varphi=\beta+\pi/2$ achieves a better $\dot{\overline{E}}$
then the measurement $\nu$. The last fact can be explained by considering the
following relation:

$\displaystyle\frac{\dot{\overline{E}}_{\rm
hom}^{(\varphi)}+\dot{\overline{E}}_{\rm
hom}^{(\varphi+\pi/2)}}{2}=\frac{|a_{0}|^{2}-\sum_{i}|a_{i}|^{2}}{{\ln 2}}.$
(S118)

Since $\dot{\overline{E}}_{\rm hom}$ is a real value, either
$\dot{\overline{E}}_{\rm
hom}^{(\beta)}\leq(|a_{0}|^{2}-\sum_{i}|a_{i}|^{2})/{\ln 2}$ for either
$\beta=\varphi$ or $\beta=\varphi+\pi/2$. Thus there exists no non-vacuum-
resolving measurement that could give a lower $\dot{\overline{E}}$ than the
ideal homodyne measurement.

####  III.7.2 Vacuum-resolving case

Note that $\dot{\overline{E}}_{\rm VR}\leq 0$, if $c$ has support only in $A$
or $B$. Let us again define the variables $r$, $\beta$, and $s$:

$\displaystyle\sum_{x|\zeta_{0}(x)\neq
0}|\zeta_{1}(x)|^{2}e^{-2i\mu_{01}(x)}=r^{2}e^{-2i\beta},$ (S119)
$\displaystyle\sum_{x|\zeta_{0}(x)\neq 0}|\zeta_{1}(x)|^{2}=s^{2},$ (S120)
$\displaystyle\sum_{x|\zeta_{0}(x)=0}|\zeta_{1}(x)|^{2}=1-s^{2}.$ (S121)

Here $0\leq r\leq s\leq 1$, and $0\leq\beta<2\pi$. Again, introducing
$a_{0}=\mathrm{tr}\left(c\phi\right)$ and
$a_{i}=a_{jk}=\langle\xi_{j}|\mathrm{tr}_{B}\left(c\phi\right)|\xi_{k}\rangle$,
we can have:

$\displaystyle\dot{\overline{E}}_{\rm VR}$ $\displaystyle=\frac{1}{2\ln
2}\left[\left(r^{2}e^{-2i\beta}a_{0}^{2}+2|a_{0}|^{2}+r^{2}e^{2i\beta}{a_{0}^{\ast}}^{2}\right)-\sum_{i}\left(r^{2}e^{-2i\beta}a_{i}^{2}+2s^{2}|a_{i}|^{2}+r^{2}e^{2i\beta}{a_{i}^{\ast}}^{2}\right)\right]+\left(1-s^{2}\right)\dot{\overline{E}}_{\rm
num},$ (S122)

which can be simplified to:

$\displaystyle\dot{\overline{E}}_{\rm VR}$
$\displaystyle=r^{2}\dot{\overline{E}}_{\rm
hom}^{(\varphi=\beta)}+\frac{s^{2}-r^{2}}{\ln
2}\left(|a_{0}|^{2}-\sum_{i}|a_{i}|^{2}\right)+\left(1-s^{2}\right)\dot{\overline{E}}_{\rm
num}.$ (S123)

The above expression is a monotonic function of $r$ in $0\leq r\leq s$ and a
monotonic function of $s$ in $0\leq s\leq 1$. Thus the minimum has to be on
the boundary of the region $0\leq r\leq s\leq 1$. To proceed we have to
distinguish the two cases:

  1. (i)

$\dot{\overline{E}}_{\rm
hom}^{(\varphi=\beta)}\leq(|a_{0}|^{2}-\sum_{i}|a_{i}|^{2})/\ln 2$, and

  2. (ii)

$\dot{\overline{E}}_{\rm
hom}^{(\varphi=\beta)}\geq(|a_{0}|^{2}-\sum_{i}|a_{i}|^{2})/\ln 2$.

Let us start with the case (i). In this case $\dot{\overline{E}}$ is a
monotonically decreasing function of $r$. Thus the minimum is achieved at the
largest value of $r$, i.e. for $r=s$. Therefore we have:

$\displaystyle\dot{\overline{E}}=s^{2}\dot{\overline{E}}_{\rm
hom}^{(\varphi=\beta)}+(1-s^{2})\dot{\overline{E}}_{\rm num}$ (S124)

This is a monotonic function of $0\leq s\leq 1$. The minimum is thus obtained
at the boundary, where the expression is either $\dot{\overline{E}}_{\rm
hom}^{(\varphi=\beta)}$ or $\dot{\overline{E}}_{\rm num}$.

In the case (ii), $\dot{\overline{E}}$ is a monotonically increasing function
of $r$. That is, it’s minimum is obtained at $r=0$. Therefore we have:

$\displaystyle\dot{\overline{E}}=s^{2}\frac{|a_{0}|^{2}-\sum_{i}|a_{i}|^{2}}{\ln
2}+(1-s^{2})\dot{\overline{E}}_{\rm num}$ (S125)

This is again a monotonic function in $s$. It’s minimum on $0\leq s\leq 1$ is
the smallest of $2(|a_{0}|^{2}-\sum_{i}|a_{i}|^{2})$ and
$\dot{\overline{E}}_{\rm num}$. In case when $\dot{\overline{E}}_{\rm num}$ is
smaller than $2(|a_{0}|^{2}-\sum_{i}|a_{i}|^{2})$, we choose $s=0$. In the
opposite case, note that from the non-vacuum-resolving case above we have:

$\displaystyle\dot{\overline{E}}_{\rm
hom}^{(\varphi=\beta+\pi/2)}\leq\frac{|a_{0}|^{2}-\sum_{i}|a_{i}|^{2}}{\ln
2}.$ (S126)

Thus, no matter what measurement one makes, either the number or homodyne
measurement will give the minimal value for $\dot{\overline{E}}$.

##  IV Matrix product state representation and stochastic updates

In this section we consider in more detail the evolution of an open many-body
quantum system represented by an MPS. We describe the state of the system of
$n$ constituents via so-called $\Gamma\lambda$-form of the MPS (see Ref. [3]
in the main text):

$\displaystyle|\phi\rangle=\sum_{i_{1},i_{2},i_{3},...,i_{n}}c_{i_{1},i_{2},i_{3},...,i_{n}}|i_{1},i_{2},i_{3},...,i_{n}\rangle,$
(S127) $\displaystyle\textrm{where
}c_{i_{1},i_{2},i_{3},...,i_{n}}=\sum_{\alpha,\beta,\gamma,...,\chi,\zeta}\Gamma_{\alpha}^{i_{1}}[1]\lambda_{\alpha\alpha}[1]\Gamma_{\alpha\beta}^{i_{2}}[2]\lambda_{\beta\beta}[2]\textrm{
}...\textrm{
}\lambda_{\xi\xi}\Gamma_{\xi\zeta}^{i_{n-1}}[n-1]\lambda_{\zeta\zeta}[n-1]\Gamma_{\zeta}^{i_{n}}[n],$
(S128)

where $i_{j}$ are the state indices of $j^{\mathrm{th}}$ constituent,
$\Gamma[j]$ and $\lambda[j]$ are a three- and two-dimensional tensors (see
Fig. S2a). The index $i_{j}$ in the superscript of the $\Gamma[j]$ tensor is
called a physical index and has dimension $d$, indices in the subscript of
$\Gamma[j]$ and $\lambda[j]$ tensors correspond to bond indices and have
dimension $\chi$.

Figure S2: Main diagramatic MPS identities used in the paper. (a)
$\Gamma\lambda$-form of the state. (b) TEBD-update of two neighbouring sites
of the state. (c) Single-site (local) update of the state. (d) Example of MPS
contractions needed to calculate the entranglement entropy change rate
$\dot{\bar{E}}_{j}$.

As discussed in the main text, trajectories are propagated by Trotterizing the
coherent and incoherent evolution, leading to alternating layers of state
updates as depicted in Fig. 2a of the main text.

The coherent part is governed by (generally) nearest-neighbour Hamiltonian
$H_{\mathrm{sys}}$. The technique to update the state of the system with such
a Hamiltonian is called time-evolving block decimation (TEBD, see Ref. [4] in
the main text). This technique consists in contracting the physical indices of
two-local unitary tensors (reshaped unitary matrices) with the corresponding
physical indices of the state, as shown in Fig. 2a in the main text.
Mathematically one step of this technique is implemented as follows (Fig.
S2b):

  1. 1.

Apply the unitary tensor to two neighbouring $\Gamma[j]$ and $\Gamma[j+1]$
matrices (Fig. S2b, first step):

$\displaystyle U|\phi\rangle=$
$\displaystyle\sum_{i^{\prime}_{j},i^{\prime}_{j+1}}\sum_{i_{1},...,i_{j},i_{j+1},...,i_{n}}\sum_{\alpha,...,\zeta}\Gamma_{\alpha}^{i_{1}}[1]\textrm{
}...\textrm{
}U_{i_{j},i_{j+1}}^{i^{\prime}_{j},i^{\prime}_{j+1}}\lambda_{\theta\theta}[j-1]\Gamma_{\theta\iota}^{i_{j}}[j]\lambda_{\iota\iota}[j]\Gamma_{\iota\kappa}^{i_{j+1}}[j+1]\lambda_{\kappa\kappa}[j+1]\textrm{
}...\textrm{ }\Gamma_{\zeta}^{i_{n}}[n]|i_{1},i_{2},i_{3},...,i_{n}\rangle=$
(S129)
$\displaystyle\sum_{i_{1},...,i^{\prime}_{j},i^{\prime}_{j+1},...,i_{n}}\sum_{\alpha,...,\zeta}\Gamma_{\alpha}^{i_{1}}[1]\textrm{
}...\textrm{ }\Theta_{\theta\kappa}^{i^{\prime}_{j},i^{\prime}_{j+1}}\textrm{
}...\textrm{ }\Gamma_{\zeta}^{i_{n}}[n]|i_{1},i_{2},i_{3},...,i_{n}\rangle;$
(S130)

  2. 2.

Perform a singular value decomposition (SVD) on tensor $\Theta$ to obtain
updated versions of $\Gamma[j]$, $\lambda[j]$ and $\Gamma[j+1]$ tensors (Fig.
S2b, second and third steps);

  3. 3.

Go to the next pair of matrices $\Gamma[j+2]$ and $\Gamma[j+3]$ and repeat the
steps.

This routine is first applied to all the odd pairs of $\Gamma$-matrices and
then to all the even pairs.

The incoherent part of the each layer is generated by applying a sequence of
single-site operations (denoted by $K_{j}$ in the main text) to the state of
the system. In the MPS framework the application of the operator $K_{j}$ is
achieved by contracting the physical indices with the corresponding index of
the $\Gamma[j]$-matrix (Fig. S2c, first step), resulting in an update of the
$\Gamma[j]$-matrix (Fig. S2c, second step):

$\displaystyle O|\phi\rangle=$
$\displaystyle\sum_{i^{\prime}_{j}}\sum_{i_{1},...,i_{j},i_{j+1},...,i_{n}}\sum_{\alpha,...,\zeta}\Gamma_{\alpha}^{i_{1}}[1]\textrm{
}...\textrm{
}O_{i_{j}}^{i^{\prime}_{j}}\Gamma_{\theta\iota}^{i_{j}}[j]\textrm{
}...\textrm{ }\Gamma_{\zeta}^{i_{n}}[n]|i_{1},i_{2},i_{3},...,i_{n}\rangle=$
(S131)
$\displaystyle\sum_{i_{1},...,i^{\prime}_{j},...,i_{n}}\sum_{\alpha,...,\zeta}\Gamma_{\alpha}^{i_{1}}[1]\textrm{
}...\textrm{ }{\Gamma^{\prime}}_{\theta\iota}^{i^{\prime}_{j}}[j]\textrm{
}...\textrm{ }\Gamma_{\zeta}^{i_{n}}[n]|i_{1},i_{2},i_{3},...,i_{n}\rangle.$
(S132)

This procedure is applied to all sites, after which the MPS is not in
canonical form and not normalized. At the end of the stochastic layer, the MPS
is renormalized and the canonical $\Gamma\lambda$-form is reinstated by TEBD
technique with an identity unitary matrix (Fig. S2c, last step, see Ref. [3]
in the main text).

The choice of single-site operations $K_{j}$ depends on a specific type of
unravelling. In the EOQT method that we propose in the main text, the choices
of $K_{j}$ are time-, trajectory- and site-dependent, since they are based on
the minimization of the EAEE change rate $\dot{\bar{E}}_{j}$ [for example, see
Eqs. (10) and (11) in the main text]. Note that $\dot{\bar{E}}_{j}$ can be
efficiently calculated within the framework of MPSs. To see this, consider,
for example, the following term from Eq. (10) in the main text:

$\displaystyle\mathrm{tr}\bigg{\\{}\mathrm{tr}_{B}\left(c\phi
c^{\dagger}\right)\ln\left[\mathrm{tr}_{B}\left(c\phi
c^{\dagger}\right)\right]\bigg{\\}},$ (S133)

the diagram of which is shown in Fig. S2d. First, the MPS representation of a
density matrix $\phi$ is constructed, to which two operators are applied to
obtain $c\phi c^{\dagger}$ (Fig. S2d, step one). Then the partial trace over
$B$ is taken to obtain $\mathrm{tr}_{B}\left(c\phi c^{\dagger}\right)$ (Fig.
S2d, step two and three). Finally, the logarithm of eigenvalues of
$\mathrm{tr}_{B}\left(c\phi c^{\dagger}\right)$ is taken and tracing over $A$
of two objects, $\mathrm{tr}_{B}\left(c\phi c^{\dagger}\right)$ and
$\ln\left[\mathrm{tr}_{B}\left(c\phi c^{\dagger}\right)\right]$, is performed
to obtain Eq. (S133). All the other terms entering $\dot{\bar{E}}_{j}$ can be
calculated in a similar way.

##  V Numerical case studies

In this section we apply the EOQT method as well as other QT methods to
several examples in order to illustrate the convergence properties on concrete
physical models.

###  V.1 Convergence with number of trajectories

Figure S3: Convergence of various QT approaches to the ME solution on the
example of the nearest-neighbour Ising model with dephasing [see Eq. (S134)].
We show the evolution of the excited state population of the second atom in
the chain, $\rho^{ee}_{2}(t)$, and correlation
$\left<\sigma^{z}_{1}\sigma^{z}_{2}\right>-\left<\sigma^{z}_{1}\right>\left<\sigma^{z}_{2}\right>$,
in case of the quantum jump method [(a) and (d)], the quantum state diffusion
method [(b) and (e)] and the EOQT method [(c) and (f)], for $N=10$, $N=10^{2}$
and $N=10^{4}$ trajectories in each case. The black dashed line indicates the
exact ME solution, which coincides with the ensemble averages for $N=10^{4}$
trajectories. The initial state is the fully polarized state with all atoms in
the state $|1\rangle$, and the system parameters are $g/\gamma=2.5$,
$J/\gamma=0.5$, $h/\gamma=-0.5$, and $n=m=4$.  Figure S4: Comparison of noise
properties of existing QT methods. Stochastic fluctuations of mean excited
state population of the second atom in the chain, $\sigma^{2}_{ee}(t)$, in
case of quantum jump method (blue), quantum state diffusion method (red) and
EOQT method (green) calculated by sampling $N=10^{4}$ trajectories for each
method. One can notice that EOQT method inherits noise properties from both
number and homodyne measurement. The system with a number of constituents
$n=4$ is initialized in a product state $|\psi\rangle=|1111\rangle$.
Parameters of the system are: $g/\gamma=2.5$, $J/\gamma=0.5$, $h/\gamma=-0.5$.

We first apply various QT methods to the paradigmatic example of the open 1D
Ising models. This serves, first, to illustrate numerically that the all QT
methods indeed converge to the exact solution of the ME and, second, to
provide an analysis of the statistical convergence properties on a well-
studied example. Specifically, we consider an open 1D spin chain described by
the ME $\dot{\rho}=-i[H_{\rm sys},\rho]+\gamma\sum_{j}(c_{j}\rho
c_{j}^{\dagger}-\frac{1}{2}c_{j}^{\dagger}c_{j}\rho-\frac{1}{2}\rho
c_{j}^{\dagger}c_{j})$, whose coherent part is given by a nearest-neighbor
Ising interaction with transverse and longitudinal fields:

$\displaystyle
H_{\mathrm{sys}}=\sum_{j=1}^{n}\left[h\sigma_{j}^{z}-g\sigma_{j}^{x}\right]+\sum_{j=1}^{n-1}J\sigma_{j}^{z}\sigma_{j+1}^{z},$
(S134)

and the incoherent part is generated by the operators
$c_{j}=|0\rangle_{j}\,\langle 1|$.

In Fig. S3 we show exemplary results obtained from the exact integration of
the ME and compare it with the QT methods, where the trajectories are produced
by continuous number measurements [Fig. S3a and S3d], quantum state diffusion
method that corresponds to homodyne measurements [Fig. S3b and S3e], and the
EOQT method that utilizes the optimal propagator explained in the main text
[Fig. S3c and S3f]. As one can see, all of the possible trajectories converge
to the ME solution (black dashed lines in Fig. S3), in case when the number of
trajectories $N$ is sufficiently large. For all trajectory methods, the
statistical error of a linear observable $O$ decreases with the number of
trajectories as $\sim C_{O}/\sqrt{N}$. The prefactor $C_{O}$ depends in
general on the observable of interest, but, more importantly for the present
discussion, also on the QT method. This can be quantified by the variance of
expectation values of interest across trajectories,
$C_{O}^{2}\sim\sigma_{O}^{2}=\sum_{i}p^{(i)}\langle\phi^{(i)}|O|\phi^{(i)}\rangle^{2}-\left(\sum_{i}p^{(i)}\langle\phi^{(i)}|O|\phi^{(i)}\rangle\right)^{2}$,
shown in Fig. S4. This plot illustrates that, even though all types of the
trajectories converge to the ME solution, the convergence rates differ for
various environment monitoring methods.

###  V.2 Convergence with bond dimension

The main idea underlying the design of our EOQT method is the expectation that
it allows to solve the ME with trajectories that require a reduced bond
dimension (compared to other QT methods). To illustrate this feature on a
concrete physical model, we choose to consider in the following a model of an
open 1D chain of 3-level atoms (with states denoted
$|g_{1}\rangle,|g_{2}\rangle$ and $|r\rangle$). We already demonstrated this
in the main text on the example of an open RBC (Fig. 3 of the main text),
where some QT methods lead to volume-law-entangled trajectories (MPSs with
finite bond dimension will fail to represent these trajectories accurately),
while other methods, including EOQT, generate area-law-entangled trajectories,
allowing for an accurate representation by an MPS with finite bond dimension.
While the open RBC allows to elegantly demonstrate this feature semi-
analytically, other aspects regarding the convergence with bond dimension are
less transparent in this model: Since the open RBC quickly generates a
maximally mixed state, uncontrolled errors in the representation of individual
trajectories due to bond dimension limitations often do not result in errors
in physically relevant observables, since the state is maximally mixed. We
therefore choose to illustrate the difference of various QT schemes regarding
the convergence of physically relevant observables with the bond dimension on
a model that does not suffer from this shortcoming. The coherent part of the
evolution is given by

$\displaystyle
H_{\mathrm{EIT}}=-\sum_{j=1}^{n}\left[\frac{\Omega_{1}}{2}\left(\sigma^{+}_{1,j}+\sigma^{-}_{1,j}\right)+\frac{\Omega_{2}}{2}\left(\sigma^{+}_{2,j}+\sigma^{-}_{2,j}\right)\right]+\sum_{j=1}^{N-1}V\sigma^{z}_{1,j}\otimes\sigma^{z}_{1,j+1},$
(S135)

where $\sigma^{-}_{i,j}=|g_{i}\rangle_{j}\langle r|$ is the spin-lowering
operator of the transition $|g_{i}\rangle\leftrightarrow|r\rangle$ of atom
$j$, and $\sigma^{+}_{i,j}=(\sigma^{-}_{i,j})^{\dagger}$. The incoherent part
is induced by jump operators of the form $c_{j}=|r\rangle_{j}\langle r|$ (i.e.
dephasing channels). This is a physically relevant model describing a Rydberg
atom array driven under conditions of electromagnetically induced transparency
(EIT). In Fig. S5 we show the dynamics of this chain obtained from an exact
solution using the ME and compare it with various QT schemes using MPSs with
different values of the maximum bond dimension $\chi$. As one can see, all QT
schemes reproduce the result of the ME in the limit of large bond dimension
(as expected). However, when the maximum bond dimension is reduced, some QT
methods result in significant errors, while others, including the EOQT method,
still produce accurate results.

Figure S5: Comparison between different unravellings in terms of their ability
to solve the ME at a given bond dimension $\chi$. Population of
$|g_{1}\rangle$ state of the first atom in the chain, produced by homodyne
measurements with phases $\varphi_{j}=0$ (red solid lines) and
$\varphi_{j}=\pi/2$ (grey solid lines), number measurements (blue solid lines)
and EOQT-measurements (green solid lines) for bond dimension (a) $\chi=9$ and
(b) $\chi=1$. The dashed line corresponds to the ME solution. One can see that
trajectories with $\varphi_{j}=\pi/2$ fail to solve the ME, while trajectories
with $\varphi_{j}=0$ as well as optimal trajectories do not. (c) EAEE of the
ensembles produced by homodyne measurements with phases $\varphi_{j}=0$ (red
solid line) and $\varphi_{j}=\pi/2$ (grey solid line), number measurements
(blue solid line) and EOQT-measurements (green solid line) for $\chi=9$. (d)
The error in the population of $|g_{1}\rangle$ state of the first atom in the
chain, $\sigma^{2}_{\rho_{1}^{g_{1}g_{1}}}$, obtained for the time instant
$\gamma t=17.5$, produced by homodyne measurements with phases $\varphi_{j}=0$
(red solid line) and $\varphi_{j}=\pi/2$ (grey solid line), number
measurements (blue solid line) and EOQT-measurements (green solid line). In
this figure: Number of atoms and dephasing channels $n=m=4$ [for (a), (b) and
(c)] and $n=m=20$ [for (d)], number of trajectories $N=10^{4}$,
$\Omega_{1}=\Omega_{2}=0.5\gamma$ and $V=\gamma$, where $\gamma$ is the
decoherence rate.

##  VI Further details on open Random Brownian circuits

In this section we provide further details on some points in the analysis of
the RBC. To start, consider the stochastic propagator associated with the
homodyne measurement for a jump operator $c_{j}=\sigma^{z}_{j}$ and equal
decay rates $\gamma_{j}=\gamma$ [Eq. (5) in the main text]:

$\displaystyle K_{j}=$ $\displaystyle e^{-\gamma
dt/2}+\sqrt{\gamma}\sigma^{z}_{j}e^{i\varphi_{j}}\left[2\sqrt{\gamma}\left<\sigma^{z}_{j}\right>(t)\cos{\varphi_{j}}dt+dW_{j}(t)\right],$
(S136)

where
$\left<\sigma^{z}_{j}\right>(t)=\langle\phi_{j}^{(k)}(t)|\sigma^{z}_{j}|\phi_{j}^{(k)}(t)\rangle$
($k$ is the trajectory index) and $dW_{j}(t)$ is a Wiener increment. In the
discussion of the RBC in the main text we use a different form of the
propagator Eq. (S136), given by
$\exp\left\\{e^{i\varphi_{j}}\sqrt{\gamma}\sigma^{z}_{j}\left[2\sqrt{\gamma}\left<\sigma^{z}_{j}\right>(t)\cos{\varphi_{j}}dt+dW_{j}(t)\right]\right\\}$.
To show the equivalence between the two propagators, let us consider the
expansion:

$\displaystyle\exp\left\\{e^{i\varphi_{j}}\sqrt{\gamma}\sigma^{z}_{j}\left[2\sqrt{\gamma}\left<\sigma^{z}_{j}\right>(t)\cos{\varphi_{j}}dt+dW_{j}(t)\right]\right\\}$
(S137)
$\displaystyle=\mathbb{1}+e^{i\varphi_{j}}\sqrt{\gamma}\sigma^{z}_{j}\left[2\sqrt{\gamma}\left<\sigma^{z}_{j}\right>(t)\cos{\varphi_{j}}dt+dW_{j}(t)\right]+\frac{1}{2}e^{i2\varphi_{j}}dW^{2}_{j}(t)\gamma\mathbb{1}+\mathcal{O}\left(dt^{3/2}\right)$
(S138)
$\displaystyle=\mathbb{1}+e^{i\varphi_{j}}\sqrt{\gamma}\sigma^{z}_{j}\left[2\sqrt{\gamma}\left<\sigma^{z}_{j}\right>(t)\cos{\varphi_{j}}dt+dW_{j}(t)\right]-\frac{\gamma
dt}{2}\mathbb{1}+\frac{\gamma}{2}\left[dt+e^{i2\varphi_{j}}dW^{2}_{j}(t)\right]\mathbb{1}+\mathcal{O}\left(dt^{3/2}\right)$
(S139)
$\displaystyle=K_{j}+\frac{\gamma}{2}\left[dt+e^{i2\varphi_{j}}dW^{2}_{j}(t)\right]K_{j}+\mathcal{O}\left(dt^{3/2}\right)$
(S140)
$\displaystyle=K_{j}\exp{\left\\{\frac{\gamma}{2}\left[dt+e^{i2\varphi_{j}}dW^{2}_{j}(t)\right]\right\\}}+\mathcal{O}\left(dt^{3/2}\right).$
(S141)

The two propagators are thus equivalent up to some higher-order terms
vanishing in the $dt\rightarrow 0$ limit and a factor $f=\exp{\left[\gamma
dt/2+\gamma e^{i2\varphi_{j}}dW^{2}_{j}(t)/2\right]}$. The factor $f$ is a
complex number that only affects the normalization and the global phase of the
state. Since after each update the state in normalized, this factor is
irrelevant. Therefore the propagation with $K_{j}$ is equivalent to
propagation with
$\exp\left\\{e^{i\varphi_{j}}\sqrt{\gamma}\sigma^{z}_{j}\left[2\sqrt{\gamma}\left<\sigma^{z}_{j}\right>(t)\cos{\varphi_{j}}dt+dW_{j}(t)\right]\right\\}$.