# A new multi-partite entanglement measure and its holographic dual

Abhijit Gadde  Vineeth Krishna  Trakshu Sharma
Department of Theoretical Physics,
Tata Institute of Fundamental Research

###### Abstract

In this letter we define a natural generalization of the von Neumann entropy
to multiple parties that is symmetric with respect to all the parties. We call
this measure multi-entropy. We show that for conformal field theories with
holographic duals, the multi-entropy is computed by the area of an appropriate
“soap-film” anchored on the boundary. We conjecture the quantum version of
this prescription that takes into account the sub-leading corrections in
$G_{N}$.

Entanglement, Multi-partite, Holography

††preprint: TIFR/TH/22-34

##  I Introduction

In recent years, quantum information theoretic notions such as entanglement
entropy [1, 2, 3, 4, 5, 6] have been immensely useful in shedding light on
some of the important questions pertaining to quantum gravity . In particular,
it has helped understand holographic encoding of the gravitational Hilbert
space and in turn, the Black-hole entropy [7, 8] and the information paradox.
However, with the exception of reflected entropy [9] and negativity [10, 11],
almost all the discussion has been centered on the entanglement between two
parties i.e. on bi-partite entanglement. It is expected that knowing multi-
partite entanglement structure in quantum gravity would refine our
understanding of the holographic encoding of gravitational Hilbert space in
conformal field theory Hilbert space. To this end, we introduce a multi-
partite entanglement measure and its holographic dual.

##  II The Measure

Let $|\Psi\rangle\in\otimes_{{\mathtt{a}}={1}}^{\mathtt{q}}{\cal
H}_{\mathtt{a}}$ be a quantum state of a ${\mathtt{q}}$-party system. Let
${\rm dim}\,{\mathcal{H}}_{\mathtt{a}}=d_{\mathtt{a}}$ and
$|\alpha_{\mathtt{a}}\rangle,\,\alpha_{\mathtt{a}}\in\\{1,\ldots,d_{\mathtt{a}}\\}$
be its orthonormal basis. In $|\alpha_{\mathtt{a}}\rangle$ basis, the state
$|\Psi\rangle$ is given as

$\displaystyle|\Psi\rangle=\sum_{\alpha_{1}=1}^{d_{1}}\ldots\sum_{\alpha_{\mathtt{q}}=1}^{d_{\mathtt{q}}}\,\psi_{\alpha_{1}\ldots\alpha_{\mathtt{q}}}|\alpha_{1}\rangle\otimes\ldots\otimes|\alpha_{\mathtt{q}}\rangle.$
(1)

A ${\mathtt{q}}$-party entanglement measure is the information contained in
the wavefunction $\psi$ (and its complex conjugate) that is invariant under
“local unitary transformations”. They are products of unitary operations
${\mathcal{U}}_{\mathtt{a}}\in U(d_{\mathtt{a}})$ performed in individual
system ${\mathtt{a}}$. Naturally, such a measure is obtained by contracting
$\alpha_{\mathtt{a}}$ indices for all ${\mathtt{a}}$. Now
$\alpha_{\mathtt{a}}$ is a fundamental index on $\psi$ and an anti-fundamental
index of $\bar{\psi}$. Hence, in order to construct invariants we need an
equal number of $\psi$’s and $\bar{\psi}$’s. We will call this number the
replica number $n$. We index the replicas by the superscript $(i)$. In
particular, the Hilbert space ${\mathcal{H}}_{a}$ of the $i$-th replica is
denoted as ${\mathcal{H}}_{\mathtt{a}}^{(i)}$ and the associated basis as
$|\alpha_{\mathtt{a}}^{(i)}\rangle$. The wavefunction of the $i$-th replica is
then $\psi_{\alpha_{1}^{(i)}\ldots\alpha_{\mathtt{q}}^{(i)}}$ and similarly
its conjugate is
$\bar{\psi}^{\alpha_{1}^{(i)}\ldots\alpha_{\mathtt{q}}^{(i)}}$.

###  II.1 Bi-partite

Let us consider the following bi-partite entanglement measure with a fixed
replica number $n$

$\displaystyle{\mathcal{Z}}_{n}\equiv{\rm Tr}(\rho^{n}),\quad{\rm
where}\quad\rho_{\alpha^{(1)}_{1}}^{\alpha^{(2)}_{1}}\equiv\psi_{\alpha^{(1)}_{1}\alpha^{(1)}_{2}}\bar{\psi}^{\alpha^{(2)}_{1}\alpha^{(1)}_{2}}.$
(2)

Here $\rho$ is the density matrix for party ${1}$. With this notation, the
squared norm of $|\Psi\rangle$ is ${\mathcal{Z}}_{1}$. This measure is related
to the familiar $n$-th Renyi entropy $S_{n}$ as

$\displaystyle S_{n}=\frac{1}{1-n}{\rm
Log}({\mathcal{Z}}_{n}/{\mathcal{Z}}_{1}^{n}).$ (3)

Let us reformulate this measure using permutations acting on replicas. This
will be useful while dealing with entanglement in quantum field theories using
twist operators.

$\displaystyle{\mathcal{Z}}_{n}=(\psi_{\alpha^{(1)}_{1}\alpha^{(1)}_{2}}\psi_{\alpha^{(2)}_{1}\alpha^{(2)}_{2}}\ldots)(\bar{\psi}^{\alpha^{(1)}_{1}\alpha^{(\sigma\cdot
1)}_{2}}\bar{\psi}^{\alpha^{(2)}_{1}\alpha^{(\sigma\cdot 2)}_{2}}\ldots)$ (4)

where $\sigma$ is a cyclic permutation element acting in the space of replicas
$i=1,\ldots,n$ of the second party i.e.
$\sigma(\alpha_{2}^{(i)})=\alpha_{2}^{(i+1)}$ with
$\alpha_{2}^{(n+1)}\equiv\alpha_{2}^{(1)}$. Because this permutation acts only
on the replica copies of party-${2}$, we denote it as $\sigma_{2}$. A priori,
it would seem that we should specify two permutation elements
$(\sigma_{1},\sigma_{2})$, the first one acting on replicas of party-${1}$ and
the other acting on replicas of party-${2}$ to get the most general index
contraction, i.e.

$\displaystyle(\psi_{\alpha^{(1)}_{1}\alpha^{(1)}_{2}}\psi_{\alpha^{(2)}_{1}\alpha^{(2)}_{2}}\ldots)(\bar{\psi}^{\alpha^{(\sigma_{1}\cdot
1)}_{1}\alpha^{(\sigma_{2}\cdot 1)}_{2}}\bar{\psi}^{\alpha^{(\sigma_{1}\cdot
2)}_{1}\alpha^{(\sigma_{2}\cdot 2)}_{2}}\ldots)$ (5)

However, we can always relabel replicas of $\bar{\psi}$ so that one of the
$\sigma$’s, say $\sigma_{1}$, is brought to ${\rm id}$ form. In other words,
$(\sigma_{1},\sigma_{2})$ and $(\sigma_{1}g,\sigma_{2}g),g\in{\mathbb{S}}_{n}$
yield the same measure. We express this as an equivalence relation

$\displaystyle(\sigma_{1},\sigma_{2})\sim(\sigma_{1},\sigma_{2})\,g$ (6)

There is yet another equivalence, namely under simultaneous relabeling of
replica copies of _both_ $\psi$’s and $\bar{\psi}$’s. Such a relabeling leads
to conjugation of both $\sigma_{\mathtt{a}}$ by some $h\in{\mathbb{S}}_{n}$.
We express this as the equivalence relation

$\displaystyle(\sigma_{1},\sigma_{2})\sim h^{-1}\,(\sigma_{1},\sigma_{2})\,h.$
(7)

If we use equation (6) to set $\sigma_{1}={\rm id}$, then thanks to equation
(7), the bi-partite entanglement measure depends only on the conjugacy class
of $\sigma_{2}$. Conjugacy class in ${\mathbb{S}}_{n}$ is specified by cycle
decomposition $\\{p_{k}\\}$ where $p_{k}$ is the number of $k$-cycles such
that $\sum_{k}k\,p_{k}=n$. For this class, it is easy to see that the
resulting measure (5) is

$\displaystyle\prod_{k}\,({\rm
Tr}(\rho^{k}))^{p_{k}}=\prod_{k}\,({\mathcal{Z}}_{k})^{p_{k}}.$ (8)

This shows that the ring of bi-partite measures is generated by
${\mathcal{Z}}_{k}$’s.

###  II.2 Multi-partite

Following this discussion, we can write a general multi-partite measure by
giving permutation elements $(\sigma_{1},\ldots,\sigma_{\mathtt{q}})$ with
equivalence relations

$\displaystyle(\sigma_{1},\ldots,\sigma_{\mathtt{q}})$
$\displaystyle\sim(\sigma_{1},\ldots,\sigma_{\mathtt{q}})\cdot g$ (9)
$\displaystyle(\sigma_{1},\ldots,\sigma_{\mathtt{q}})$ $\displaystyle\sim
h^{-1}\,(\sigma_{1},\ldots,\sigma_{\mathtt{q}})\,h.$ (10)

Unlike the bi-partite case, where Renyi entropies are the only independent
measures of entanglement, for the case of three or higher number of parties,
the number of measures increase exponentially or faster with the replica
number $n$. In this letter, we will not concern ourselves with general
${\mathtt{q}}$-party measures but rather focus on a particular family of
measures labeled by an integer that is symmetric under the exchange of all the
parties.

As in the bi-partite case, equation (9) can be used to set $\sigma_{1}={\rm
id}$. But this “gauge fixing” obscures the symmetry between all the parties as
it treats party-${1}$ differently from others. So we will not use the gauge
freedom (9) and (10) just yet and specify all $\sigma_{\mathtt{a}}$’s. We will
do so with $n^{{\mathtt{q}}}$ replicas. We index the replicas with a
${\mathtt{q}}$-dimensional index vector $(i_{1},\ldots,i_{\mathtt{q}})$ where
each $i_{\mathtt{a}}=\\{1,\ldots,n\\}$. The permutation element
$\sigma_{\mathtt{a}}$ is defined as the cyclic element acting only on
$i_{\mathtt{a}}$:

$\displaystyle\sigma_{\mathtt{a}}\cdot(\ldots,i_{\mathtt{a}},\ldots)=(\ldots,i_{\mathtt{a}}+1,\ldots),\quad{\mathtt{a}}=1,\ldots,{\mathtt{q}}.$
(11)

As an element of the permutation group ${\mathbb{S}}_{n^{{\mathtt{q}}}}$, its
conjugacy class is $p_{n}=n^{{\mathtt{q}}-1}$ with all other $p_{k}=0$. It is
also clear that this measure is symmetric in all the ${\mathtt{q}}$ parties.

Now we use the gauge freedom (9) with $g=\sigma_{1}^{-1}$. This gives us an
equivalent set of permutation elements

$\displaystyle({\rm
id},\sigma_{2}\sigma_{1}^{-1},\ldots,\sigma_{\mathtt{q}}\sigma_{1}^{-1}).$
(12)

Note that the action of any of these permutations on
$(i_{1},\ldots,i_{\mathtt{q}})$ keeps the sum
$\sum_{\mathtt{a}}\,i_{\mathtt{a}}$ invariant. As a result, their action on
$n^{{\mathtt{q}}}$ replicas splits into $n$ orbits of $n^{{\mathtt{q}}-1}$
elements with $\sum_{\mathtt{a}}i_{\mathtt{a}}={\rm constant}({\rm mod}\,n)$.
Each orbit gives rise to the same invariant of $\psi$’s and $\bar{\psi}$’s. It
is convenient to work with a single orbit, say with
$\sum_{\mathtt{a}}i_{\mathtt{a}}=0({\rm mod}\,n)$. For this orbit, we can use
the relabeling freedom (10) to set the replica index $i_{1}$ to $1$. With this
gauge fixing, we get a convenient presentation of the rest of the permutation
elements
$(\sigma_{2}\sigma_{1}^{-1},\ldots,\sigma_{\mathtt{q}}\sigma_{1}^{-1})\equiv(\hat{\sigma}_{2},\ldots,\hat{\sigma}_{\mathtt{q}})$.

$\displaystyle\hat{\sigma}_{\mathtt{a}}\cdot(\ldots,i_{\mathtt{a}},\ldots)=(\ldots,i_{\mathtt{a}}+1,\ldots),\quad{\mathtt{a}}={2},\ldots,{\mathtt{q}}.$
(13)

We have denoted the gauge fixed permutation elements as
$\hat{\sigma}_{\mathtt{a}}$ to avoid confusion with un-gauge fixed permutation
elements $\sigma_{\mathtt{a}}$. The index set in equation (13) starts from
$i_{2}$ because $i_{1}$ has already been set to $1$. We denote the measure
defined by $\hat{\sigma}_{\mathtt{a}}$ elements in (13) as
${\mathcal{Z}}_{n}^{({\mathtt{q}})}$ and define Renyi multi-entropy
$S_{n}^{({\mathtt{q}})}$ as

$\displaystyle S^{({\mathtt{q}})}_{n}\equiv\frac{1}{1-n}{\rm
Log}\,({\mathcal{Z}}_{n}^{({\mathtt{q}})}/({\mathcal{Z}}_{1}^{({\mathtt{q}})})^{n^{{\mathtt{q}}-1}}).$
(14)

Here the factor $({\mathcal{Z}}_{1}^{({\mathtt{q}})})^{n^{{\mathtt{q}}-1}}$
serves to normalize the state. This is the family of ${\mathtt{q}}$-partite
measures that we are interested in. Let us pause for a moment to consider
properties of the set
$\\{\hat{\sigma}_{\mathtt{a}}:{\mathtt{a}}={2},\ldots,{\mathtt{q}}\\}$.

  * •

$\hat{\sigma}_{\mathtt{a}}$’s have the same equivalence class given by the
cycle $p_{n}=n^{{\mathtt{q}}-2}$ with all other $p_{k}=0$.

  * •

$\hat{\sigma}_{\mathtt{a}}^{-1}\hat{\sigma}_{\mathtt{b}}$ for
${\mathtt{a}}\neq{\mathtt{b}}$ also have the same equivalence class.

  * •

Together, $\hat{\sigma}_{\mathtt{a}}$s generate the subgroup
${\mathbb{Z}}_{n}^{\otimes{{\mathtt{q}}-1}}$ of the permutation group.

These observations will be important in the future discussions.

In a way, fixing $\hat{\sigma}_{1}={\rm id}$ has a very natural
interpretation. If we consider a single $\psi$ and a single $\bar{\psi}$ with
$i_{1}$ index contraction, we get the density matrix on
${\mathcal{H}}_{2}\otimes\ldots\otimes{\mathcal{H}}_{\mathtt{q}}$. This
density matrix has the index structure

$\displaystyle\rho_{\alpha_{2}^{(i_{2})}\ldots\alpha_{\mathtt{q}}^{(i_{\mathtt{q}})}}^{\alpha_{2}^{(\hat{\sigma}_{2}\cdot
i_{2})}\ldots\alpha_{\mathtt{q}}^{(\hat{\sigma}_{\mathtt{q}}\cdot
i_{\mathtt{q}})}}.$ (15)

This means the gauge fixed permutations $\hat{\sigma}_{\mathtt{a}}$’s describe
index contractions of this density matrix. Let us see this with two examples,
for ${\mathtt{q}}=2$,

$\displaystyle{\mathcal{Z}}^{({2})}_{n}=\sum\prod_{i_{2}}\,\rho_{\alpha_{2}^{(i_{2})}}^{\alpha_{2}^{(\hat{\sigma}_{2}\cdot
i_{2})}}=\sum\rho_{\alpha_{2}^{(1)}}^{\alpha_{2}^{(\hat{\sigma}_{2}\cdot
1)}}\rho_{\alpha_{2}^{(2)}}^{\alpha_{2}^{(\hat{\sigma}_{2}\cdot
2)}}\ldots\rho_{\alpha_{2}^{(n)}}^{\alpha_{2}^{(\hat{\sigma}_{2}\cdot n)}}.$

Here sum is over all repeated indices. With $\hat{\sigma}_{2}$ being a cyclic
element as specified by equation (13), this is nothing but ${\mathcal{Z}}_{n}$
defined earlier. So we have showed $S_{n}^{({2})}=S_{n}$. Graphically,
denoting $\rho$ as in figure 1, we get the graphical representation of
$S^{({2})}_{n}$ as a circular lattice of $\rho$’s of length $n$ as in figure
1. Each directed link in this lattice represents index contraction from a
fundamental index to an anti-fundamental index. The advantage of this
graphical notation is that it admits straightforward extension to higher
number of parties.

Figure 1: Graphical notation for the density matrix for ${\mathtt{q}}=2$ and
the index contraction that gives ${\mathcal{Z}}_{4}^{(2)}$. Figure 2:
Graphical notation for the density matrix for ${\mathtt{q}}=3$ and the index
contraction that gives ${\mathcal{Z}}_{4}^{(3)}$

For instance, in the case of ${3}$-parties, we denote the density matrix as in
figure 2. Then ${\mathcal{Z}}^{({3})}_{n}$ is the index contraction given by
the toric lattice of $\rho$’s with both sides being of length $n$. This is
shown in figure 2. Similarly, ${\mathcal{Z}}^{({\mathtt{q}})}_{n}$ is the
index contraction given by $({\mathtt{q}}-1)$-dimensional toric lattice with
all sides being of length $n$.

Arguably the most important measure of bi-partite entanglement, the von-
Neumann entropy $S$, is obtained as a limit

$\displaystyle S={\rm lim}_{n\to{1}}\,S_{n}.$ (16)

The Renyi multi-entropy $S^{({\mathtt{q}})}_{n}$ is defined such that the
limit

$\displaystyle S^{({\mathtt{q}})}={\rm
lim}_{n\to{1}}\,S^{({\mathtt{q}})}_{n},$ (17)

called _the_ multi-entropy, has a number of nice properties.

  1. 1.

It is symmetric in all the parties.

  2. 2.

For ${\mathtt{q}}=2$, it reduces to the von-Neumann entropy.

  3. 3.

It admits a convenient holographic description.

Properties $1$ and $2$ are obvious from the discussion so far. In the rest of
the paper, we will discuss multi-entropy for holographic theories and describe
what we exactly mean by property $3$. The multi-entropy can also be computed
from ${{\mathcal{Z}}}_{n}^{({\mathtt{q}})}$ using the formula

$\displaystyle S^{({\mathtt{q}})}=-\partial_{n}{\rm
log}\Big{(}{\mathcal{Z}}_{n}^{({\mathtt{q}})}/({\mathcal{Z}}_{1}^{({\mathtt{q}})})^{n^{{\mathtt{q}}-1}}\Big{)}|_{n=1}.$
(18)

We will use this formula and equation (17) interchangeably.

###  II.3 A note on analytic continuation

We have defined the multi-entropy as a limit of the Renyi multi-entropy in
equation (17). This requires analytic continuation of $S_{n}^{({\mathtt{q}})}$
away from integer $n$. One may rightly wonder whether such an analytic
continuation exists and under what conditions is it unique. The issue of
existence was recently emphasized in [12]. The uniqueness issue is usually
tackled using the Carlson’s theorem. It says

  * •

A function $f(z)$ that is analytic for ${\rm Re}(z)>0$ and satisfies

$\displaystyle f(|z|)\leq Ae^{c|z|},\quad f(iy)\leq Ae^{\pi y}$ (19)

for some real constants $A$ and $c$ for all $z\in{\mathbb{C}}$ and
$y\in{\mathbb{R}}$ and vanishes for any non-negative integer must be
identically zero.

In order to argue for the uniqueness of the analytic continuation (if it
exists) we must look for analytic functions obeying the analyticity and
boundedness conditions required by the Carlson’s theorem, henceforth called
the Carlson conditions. The problem of constructing such a function from its
values at non-negative integers or “input data” was considered in a paper by
Regge and Viano [13]. There this problem was motivated from the theory of
complex angular momentum. The authors obtain a solution that is even amenable
to numerical analysis. However, it is not clear whether their proposed
analytic continuation obeys the Carlson conditions. This question was
expounded upon in the paper [14]. In this paper, authors propose an elaborate
set of conditions on the input data that is necessary for the Carlson
conditions to hold 111We would like to thank Geoff Penington to bring this
paper to our attention.. We are currently exploring if these conditions
satisfied Renyi multi-entropies of qubit states. If these conditions hold for
general quantum states then the multi-entropy would be unambiguously defined.
For the purposes of this paper, we will define multi-entropy only for those
states whose Renyi multi-entropy admit a manifest analytic continuation. As we
will see below, this is the case for holographic states.

###  II.4 Example

Although the multi-entropy is amenable to holographic computation as we will
see shortly, it is difficult to compute for finite dimensional quantum
systems. This is because multi-entropy is defined only via the replica trick
and, in general, it is difficult to compute
${\mathcal{Z}}_{n}^{({\mathtt{q}})}$ as an analytic function of $n$. However,
for a special class of states that we call generalized GHZ states, the multi-
entropy is readily computed. Consider
$|\Psi\rangle\in\otimes_{{\mathtt{a}}=1}^{\mathtt{q}}{\mathcal{H}}_{\mathtt{a}}$
with ${\rm dim}\,{\mathcal{H}}_{{\mathtt{a}}}=d$. The generalized GHZ state is
defined as

$\displaystyle|\Psi\rangle_{\rm
GHZ}=\sum_{i}\lambda_{i}|i\rangle\otimes\ldots\otimes|i\rangle.$ (20)

Here $|i\rangle$ is an orthonormal basis in ${\mathcal{H}}_{{\mathtt{a}}}$ for
all ${\mathtt{a}}$. The state is normalized so $\sum_{i}|\lambda_{i}|^{2}=1$.
It is easy to compute ${\mathcal{Z}}_{n}^{({\mathtt{q}})}$ because the same
index $i$ runs through all the contractions. It contributes
$|\lambda_{i}|^{2n^{{\mathtt{q}}-1}}$. Then the ${\mathtt{q}}$-Renyi entropy
is

$\displaystyle S_{n}^{({\mathtt{q}})}=\frac{1}{1-n}{\rm
log}\Big{(}\sum_{i}|\lambda_{i}|^{2n^{{\mathtt{q}}-1}}\Big{)}.$ (21)

Taking the $n\to 1$ limit,

$\displaystyle
S^{({\mathtt{q}})}=(1-{\mathtt{q}})\sum_{i}|\lambda_{i}|^{2}{\rm
log}|\lambda_{i}|^{2}.$ (22)

A bi-partite state can always be taken to the generalized GHZ form via Schmidt
decomposition. Then $|\lambda_{i}|^{2}$ are the eigenvalues of the density
matrix. It is clear that $S^{({\mathtt{q}})}$ agrees with the von-Neumann
entropy for ${\mathtt{q}}=2$. We believe that it is extremely important to
develop techniques to compute multi-entropy for general states to understand
its quantum information theoretic properties.

###  II.5 Distinguishing isospectral density matrices

In this section, we will discuss the effectiveness of Renyi multi-entropy in
distinguishing multi-partite entanglement in two mixed states. For this
purpose, we will present density matrices $\rho_{12}$ and $\sigma_{12}$ on two
qubits such that the spectrum of $\rho_{12},\rho_{1}$ and $\rho_{2}$ is
identical to that of $\sigma_{12},\sigma_{1}$ and $\sigma_{2}$ respectively.
Such density matrices are known as isospectral. This example of the
isospectral pair is borrowed from [16]. The isospectral property of these
density matrices guarantees that no bi-partite measure of entanglement can
distinguish between these states. However, as we will show, all the Renyi
multi-entropies for these two states are different.

$\displaystyle\rho_{12}=\frac{1}{3}\begin{pmatrix}1&0&0&0\\\ 0&1&1&0\\\
0&1&1&0\\\
0&0&0&0\end{pmatrix},\quad\sigma_{12}=\frac{1}{3}\begin{pmatrix}1&0&0&0\\\
0&0&0&0\\\ 0&0&0&0\\\ 0&0&0&2\end{pmatrix}$ (23)

The non-zero eigenvalues of $\rho_{12},\rho_{1},\rho_{2}$ (and also for
$\sigma_{12},\sigma_{1},\sigma_{2}$) are $\frac{1}{3},\frac{2}{3}$. The values
of the first few Renyi multi-entropies are

| $\rho$ | $\sigma$
---|---|---
$S^{(3)}_{2}$ | ${\rm Log}\left(9\right)$ | ${\rm Log}\left(\frac{81}{17}\right)$
$S^{(3)}_{3}$ | $\frac{1}{2}{\rm Log}\left(\frac{6561}{14}\right)$ | $\frac{1}{2}{\rm Log}\left(\frac{729}{19}\right)$
$S^{(3)}_{4}$ | $\frac{1}{3}{\rm Log}\left(\frac{14348907}{139}\right)$ | $\frac{1}{3}{\rm Log}\left(\frac{43046721}{65537}\right)$

In fact, it is easy to see that $\rho_{12}$ is obtained by tracing out one
party in the pure W-state $(|100\rangle+|010\rangle+|001\rangle)/\sqrt{3}$ and
$\sigma_{12}$ is obtained analogously from the pure generalized GHZ state
$(|000\rangle+\sqrt{2}|111\rangle)/\sqrt{3}$. The density matrix $\sigma_{12}$
is separable while the density matrix $\rho_{12}$ is not. The fact that the
(Renyi) multi-entropy distinguishes the two suggests that it can lead to a new
separability criterion. We are currently exploring this possibility.

##  III Multi-entropy from Holography

In this section, we will discuss computation of multi-entropy in a
$D$-dimensional conformal field theory ${\cal T}$. Let the state
$|\Psi\rangle$ be defined on a time symmetric Cauchy slice ${\mathcal{R}}$ of
a $D$-manifold ${\mathcal{M}}$. It is given by a Euclidean path integral on
the half-space ${\mathcal{M}}_{\Psi}$ such that
$\partial{\mathcal{M}}_{\Psi}={\mathcal{R}}$. The dual bra $\langle\Psi|$ is
constructed by Euclidean path integral on the other half
${\mathcal{M}}_{\bar{\Psi}}$. It is obtained from ${\mathcal{M}}_{\Psi}$ by
reflecting across ${\mathcal{R}}$. The squared norm of $\Psi$ is the partition
function ${\mathcal{Z}}_{\mathcal{M}}$ on ${\mathcal{M}}$. Let us decompose
${\mathcal{R}}$ into ${\mathtt{q}}$ number of disjoint regions
${\mathcal{R}}_{\mathtt{a}}$, such that
$\cup_{\mathtt{a}}{\mathcal{R}}_{\mathtt{a}}={\mathcal{R}}$. Let the Hilbert
space on region ${\mathcal{R}}_{\mathtt{a}}$ be ${\mathcal{H}}_{\mathtt{a}}$.
We are interested in computing multi-entropy of the state $|\Psi\rangle$ under
the decomposition $\otimes_{\mathtt{a}}{\mathcal{H}}_{\mathtt{a}}$. For
theories that admit a weakly coupled gravity dual, this problem can be
addressed holographically.

The replica trick [1] involves working with the tensor product theory ${\cal
T}^{\otimes n^{{\mathtt{q}}-1}}$ on ${\mathcal{M}}$. This theory has the
discrete symmetry ${\mathbb{S}}_{n^{{\mathtt{q}}-1}}$. It admits co-dimension
$2$ twist defects labeled by elements of this permutation group. For every
pair of regions $({\mathcal{R}}_{\mathtt{a}},{\mathcal{R}}_{\mathtt{b}})$ that
share a boundary, we insert the twist operator
${{\mathcal{O}}}_{\hat{\sigma}_{\mathtt{a}}^{-1}\hat{\sigma}_{{\mathtt{b}}}}$
on the common boundary. The measure ${\mathcal{Z}}_{n}^{({\mathtt{q}})}$ is
then given by the correlation function of these twist operators. Let us denote
the resulting replicated manifold as ${\cal M}_{n}$. The correlation function
of twist operators is the partition function
${\mathcal{Z}}_{{\mathcal{M}}_{n}}$ on ${\mathcal{M}}_{n}$. Using (14) we have

$\displaystyle S^{{\mathtt{q}}}_{n}=\frac{1}{1-n}{\rm
Log}({\mathcal{Z}}_{{\mathcal{M}}_{n}}/({\mathcal{Z}}_{\mathcal{M}})^{n^{{\mathtt{q}}-1}}).$
(24)

Following [17], we will proceed to analyze this problem holographically. Let
${\mathcal{B}}_{n}$ be dominant the gravity solutions such that
$\partial{\mathcal{B}}_{n}={\mathcal{M}}_{n}$. The manifold
${\mathcal{B}}_{n}$ is a smooth manifold e.g. in the case of Einstein gravity
with negative cosmological constant, it is of constant negative curvature. The
holographic dictionary gives

$\displaystyle{\rm Log}\,{\mathcal{Z}}_{{\mathcal{M}}_{n}}=-{\mathcal{S}}_{\rm
grav}({\mathcal{B}}_{n}).$ (25)

Here ${\mathcal{S}}_{\rm grav}({\cal X})$ is the gravitational action
evaluated on the solution ${\cal X}$. For $n=1$, this gives ${\rm
Log}\,{\mathcal{Z}}_{{\mathcal{M}}}=-{\mathcal{S}}_{\rm
grav}({\mathcal{B}}),{\mathcal{B}}\equiv{\mathcal{B}}_{1}$.

The background fields on the manifold ${\mathcal{M}}_{n}$ enjoy a replica
symmetry. This is the symmetry generated by the permutation elements
$\hat{\sigma}_{\mathtt{a}}^{-1}\hat{\sigma}_{\mathtt{b}}$ associated to all
the twist operators. As remarked earlier, this group is
${\mathbb{Z}}_{n}^{\otimes{{\mathtt{q}}-1}}$. Following [17], we will assume
that the dominant bulk solution ${\mathcal{B}}_{n}$ enjoys this symmetry. The
solution ${\mathcal{B}}_{n}$ consists of co-dimension $2$ loci that are
invariant under certain subgroups of the replica symmetry group. Some of these
loci, called “external”, are anchored at the fixed points on the boundary
(these are locations of twist operator insertions on ${\mathcal{M}}$) while
the rest are “internal”. Let us denote the loci that are anchored at the fixed
points corresponding to the twist operator ${\mathcal{O}}_{g}$ as
${\mathcal{L}}_{g}$. They are invariant under the ${\mathbb{Z}}_{n}$ subgroup
generated by $g$. Generically, two ${\mathcal{L}}$s can merge to form a
different ${\mathcal{L}}$. Merging obeys the algebra
${\mathcal{L}}_{g_{1}}\cdot{\mathcal{L}}_{g_{2}}\to{\mathcal{L}}_{g_{1}g_{2}}$.
The internal loci come about because of such merging. Below we will assume
that

  1. 1.

Every fixed point locus is of the form
${\mathcal{L}}_{\hat{\sigma}_{\mathtt{a}}^{-1}\hat{\sigma}_{\mathtt{b}}}$.

We will justify this assumption shortly. Thanks to the special property of our
permutation elements $\hat{\sigma}_{\mathtt{a}}$ stated below equation (14),
this means that any locus is invariant under _some_ ${\mathbb{Z}}_{n}$
subgroup of the replica symmetry. Every locus ${\mathcal{L}}_{g}$ appears in
groups of $n^{{\mathtt{q}}-2}$. This is because, the orbit of the action of
replica group ${\mathbb{Z}}_{n}^{\otimes{{\mathtt{q}}-1}}$ consists of
$n^{{\mathtt{q}}-2}$ elements as its stabilizer is ${\mathbb{Z}}_{n}$. For
${\mathtt{q}}=2$, the number of elements in the orbit is $1$.

We now make use of the replica symmetry in the bulk to construct the orbifold
$\widetilde{\mathcal{B}}_{n}\equiv{\mathcal{B}}_{n}/({\mathbb{Z}}_{n})^{{{\mathtt{q}}-1}}$.
Due to symmetry, the classical gravitational actions on the two spaces are
related as

$\displaystyle{\mathcal{S}}_{\rm
grav}({\mathcal{B}}_{n})=n^{{\mathtt{q}}-1}\,{\mathcal{S}}_{\rm
grav}(\widetilde{\mathcal{B}}_{n}).$ (26)

The orbifold $\widetilde{\mathcal{B}}_{n}$ has a nice property that
$\partial\widetilde{\mathcal{B}}_{n}={\mathcal{M}}$. A group of
$n^{{\mathtt{q}}-2}$ number of ${\mathcal{L}}_{g}$ become a single conical
singularity of opening angle $2\pi/n$ in the orbifold
$\widetilde{\mathcal{B}}_{n}$. Let us denote this singularity as
$\widetilde{\mathcal{L}}_{g}$. Let us denote the web created by these
singularities as ${\mathcal{W}}$. Consider a co-dimension $1$ slice
${\mathcal{C}}\in\widetilde{\mathcal{B}}_{n}$ that contains ${\mathcal{W}}$
and $\partial{\mathcal{C}}={\mathcal{R}}$. There are multiple such slices and
the precise choice doesn’t matter for the following discussion. Every
singularity becomes a co-dimension $1$ wall in ${\mathcal{C}}$ and the web
${\mathcal{W}}$ yields its chamber decomposition. As we move from
${\mathcal{R}}_{\mathtt{a}}$ to ${\mathcal{R}}_{\mathtt{b}}$ through
${\mathcal{C}}$, we must encounter at least one wall because the permutation
elements $\hat{\sigma}_{\mathtt{a}}$ and $\hat{\sigma}_{\mathtt{b}}$ are
different. At this stage, we make an assumption about ${\mathcal{W}}$ that

  1. 2.

There is no chamber which lies completely in the interior of ${\mathcal{C}}$.

As a result, we get a one-to-one map between the chambers and boundary
regions. Let us denote the chamber adjacent to ${\mathcal{R}}_{\mathtt{a}}$ as
${\mathcal{C}}_{\mathtt{a}}$. It has the property
$\partial{\mathcal{C}}_{a}\cap{\mathcal{M}}={\mathcal{R}}_{\mathtt{a}}$. The
web ${\mathcal{W}}$ consists of only those walls that separate
${\mathcal{C}}_{\mathtt{a}}$ and ${\mathcal{C}}_{\mathtt{b}}$ for some
$({\mathtt{a}},{\mathtt{b}})$. Such a wall must be of the type
$\widetilde{\mathcal{L}}_{\hat{\sigma}_{\mathtt{a}}^{-1}\hat{\sigma}_{\mathtt{b}}}$.
The parent ${\mathcal{L}}\in{\mathcal{B}}_{n}$ must also be of the same type.
This justifies our assumption 1.

To compute the ${\mathtt{q}}$-Renyi entropy using equations (25) and (26), we
need to evaluate the gravitational action on the orbifold solution
$\widetilde{{\mathcal{B}}}_{n}$. In what follows, we will specialize to the
case of Einstein gravity. The orbifold solution is a smooth solution of
Einstein’s equations with constant negative curvature except at
$\widetilde{{\mathcal{L}}}$ where it has a conical singularity of opening
angle $2\pi/n$. The gravitational action for this solution is computed as
follows [17]. We excise a small tubular neighborhood, say of radius $a$,
around the conical singularity ${\widetilde{\mathcal{L}}}$. The gravitational
action comes from the Gibbons-Hawking-York (GHY) term evaluated on the
resulting boundary. This contribution is extensive in the co-dimension $2$
area of ${\widetilde{\mathcal{L}}}$. In the limit $a\to 0$, the action takes
the form [18], $\frac{1-n}{4nG_{N}}\int dy^{D-1}\sqrt{h}$. Here $y$ is a
coordinate along ${\widetilde{\mathcal{L}}}$ and $h$ is the induced metric on
${\widetilde{\mathcal{L}}}$. To compute the equations that are obeyed by the
singular locus, it is useful to introduce a cosmic brane with the same action
supported on the singularity.

$\displaystyle{\mathcal{S}}_{\rm br}^{(n)}=\frac{1-n}{4nG_{N}}\int
dy^{D-1}\sqrt{h}=\frac{1-n}{4nG_{N}}A.$ (27)

The $y$ integral is simply the area $A$ of the brane web ${\mathcal{W}}$. The
solution is then computed by solving the equations coming from the action
${\cal S}_{\rm grav}+{\cal S}^{(n)}_{\rm br}$. It is known [19, 20] that such
a cosmic brane action indeed gives rise to a conical singularity of opening
angle $2\pi/n$ in the two transverse directions as desired. Once the solution
is found by extremizing ${\mathcal{S}}_{\rm grav}+{\mathcal{S}}_{\rm
br}^{(n)}$, we need to evaluate only the action $S_{\rm grav}$ on the
solution. This is because the cosmic brane is not actually present at the
singularity but merely used as a trick to model the singularity. The novelty
in the multi-partite case is that multiple $\widetilde{{\mathcal{L}}}_{g}$’s
can meet each other at higher co-dimensional loci. To accommodate such a
meeting we simply let the corresponding cosmic branes meet.

A priori, there can be additional terms in the brane action supported only at
the meeting locus. To compute such terms we again excise the tubular
neighborhood of ${\widetilde{\mathcal{L}}}$’s and focus near their higher co-
dimensional junction. As multiple ${\widetilde{\mathcal{L}}}$’s meet, their
corresponding tubular neighborhoods also meet forming corners as shown in
figure 3.

Figure 3: Top view of the junction of tubular neighborhoods of three
${\widetilde{\mathcal{L}}}^{\prime}s$.

For the gravitational variational principle to be well-defined, we need to add
the so-called Hayward term at the corner just as the way we add the GHY term
on the boundary.

$\displaystyle{\cal S}_{\rm Hayward}=-\frac{1}{8\pi G_{N}}\int
d\xi^{D-2}\,(\theta-\pi)\sqrt{\gamma}.$ (28)

Here $\theta$ is the angle subtended by the two boundaries at the corner,
$\xi$ is the coordinate along the corner and $\gamma$ is the induced metric on
the corner. The contribution of the Hayward term is extensive along the co-
dimension $3$ meeting locus and the proportionality constant is computed by
integrating the Hayward term along the one-dimensional corner seam of the
tubular neighborhood as shown in figure 3. This integral is proportional to
$a$ and vanishes as we take $a\to 0$. This had to be the case on dimensional
grounds because the contribution at the meeting locus must take the form
$\sim\frac{1}{G_{N}}\ell^{(3)}_{\rm meeting}\ell_{\rm scale}$. Here
$\ell^{(3)}_{\rm meeting}$ is the length (co-dimension $3$) of the meeting
locus in ${\mathcal{W}}$ and $\ell_{\rm scale}$ is some length scale needed to
obtain a dimensionless answer. The only length scale that could serve this
purpose is $a$ 222Another length scale present in the problem is the AdS scale
$\ell_{\rm AdS}$, however it is irrelevant to this completely local
computation. which we take to $0$. This argument also shows that there is no
extra contribution to the action even at higher co-dimension meeting loci.

Finding the solution to the theory ${\mathcal{S}}_{\rm
grav}+{\mathcal{S}}_{\rm br}^{(n)}$ for general $n$ is still a daunting task
(see [22, 23, 18] for computation of bi-partite Renyi entropy). But to compute
multi-entropy, only the limit $n\to 1$ is relevant. As the only $n$ dependence
appears in the coefficient in ${\mathcal{S}}_{\rm br}^{(n)}$, the solution can
be analytically continued away from $n$ integer. Moreover, in the limit $n\to
1$, the tension of the brane goes to zero and the solution can be found in the
probe limit as brane web configuration that extremizes ${\mathcal{S}}_{\rm
br}$ in the fixed background ${\mathcal{B}}$. The solution obeys the equation
of motion $\delta_{g_{\mu\nu}}({\mathcal{S}}_{\rm grav}+{\mathcal{S}}_{\rm
br}^{(n)})=0$. In particular,

$\displaystyle\partial_{n}{\mathcal{S}}_{\rm
grav}=-\partial_{n}{\mathcal{S}}_{\rm br}^{(n)}=\frac{A}{4G_{N}}.$ (29)

This, along with equations (17), (24) and (25), shows that the multi-entropy
is given a simple formula

$\displaystyle S^{({\mathtt{q}})}=\frac{A({\mathcal{W}})}{4G_{N}}.$ (30)

Here $A({\mathcal{W}})$ is the area of the minimal brane web ${\mathcal{W}}$
in ${\mathcal{B}}$. The brane web ${\mathcal{W}}$ obeys the topological
conditions,

  1. 1.

${\mathcal{W}}$ is anchored at the boundaries of all the regions
${\mathcal{R}}_{\mathtt{a}}$’s.

  2. 2.

${\mathcal{W}}$ contains sub-webs that are homologous to all the regions
${\mathcal{R}}_{\mathtt{a}}$’s.

The second condition is the reformulation of the statement that between any
two chambers ${\mathcal{C}}_{\mathtt{a}}$ and ${\mathcal{C}}_{\mathtt{b}}$
there must be at least one wall. As the solution minimizes the area subject to
these conditions, it doesn’t allow any chamber that lies completely in the
interior of ${\mathcal{C}}$. This justifies our assumption 2.

The brane web in question is extremely familiar in $D=4$. In the three
dimensional Cauchy slice ${\mathcal{C}}$ it resembles a soap film in
hyperbolic space anchored on a given “wire frame” at infinity. Specifying the
anchor does not specify the soap film uniquely but when combined with the
homology condition and the global minimum condition, it does. The close
analogy to soap-films leads us to call our prescription, the soap-film
prescription. To describe the soap-film it is convenient to introduce the
label ${\widetilde{\mathcal{L}}}^{(k)}$ for the special locus of co-dimension
$k$. It is defined inductively as the meeting locus of
${\widetilde{\mathcal{L}}}^{(k-1)}$ with
${\widetilde{\mathcal{L}}}^{(2)}\equiv{\widetilde{\mathcal{L}}}$. For $D=4$ it
is known that

  * •

Three ${\widetilde{\mathcal{L}}}^{(2)}$’s meet at
${\widetilde{\mathcal{L}}}^{(3)}$ at an angle $2\pi/3=\cos^{-1}(-1/2)$.

  * •

Four ${\widetilde{\mathcal{L}}}^{(3)}$’s meet at
${\widetilde{\mathcal{L}}}^{(4)}$ at an angle $\cos^{-1}(-1/3)$.

These are known as Plateau’s laws of soap-film. In dimension $D$ we expect
$k+1$ of ${\widetilde{\mathcal{L}}}^{(k)}$’s to meet at
${\widetilde{\mathcal{L}}}^{(k+1)}$ at an angle $\cos^{-1}(-1/k)$ and so on
until we get to ${\widetilde{\mathcal{L}}}^{(D)}$.

We have verified our proposal for $2D$ CFTs with large central charge. These
calculations will appear in an accompanying paper [24].

##  IV Generalizations

In this section we will comment on the generalizations of our soap-film
prescription in various directions.

###  IV.1 Covariant prescription

The covariant generalization of the Ryu-Takayanagi formula was presented in
[3]. It was argued there that the entanglement entropy is given by area of the
extremal surface rather than the minimal one. This prescription was later
rephrased as a maximin prescription in [25], where the surface in question is
obtained first by minimizing the area on some achronal slice $\Sigma$ and then
by maximizing the area with respect to the variation of $\Sigma$. In case
there are multiple extremal solutions, we pick the one with the minimum area
as that corresponds to the most dominant solution.

The soap-film prescription proposed here admits a natural covariant
generalization along the lines of [3] and [25]. We conjecture that the
covariant multi-entropy is obtained as the area of the minimal extremal soap-
film rather than the globally minimal one and that this choice is equivalent
to the one obtained by the maximin prescription of [25].

###  IV.2 Higher derivative gravity

For higher derivative gravitational theories, we do expect a non-trivial
contribution localized at the higher co-dimensional meeting loci
${\widetilde{\mathcal{L}}}^{(k)}$’s because the dimensional analysis argument
presented in section III ceases to be valid as $\ell_{\rm scale}$ can be
provided by the inverse mass-scale of the higher derivative corrections. It
would be interesting to compute these terms in general higher derivative
theory of gravity.

###  IV.3 Quantum corrections

Quantum correction to the multi-entropy can be found using the replica trick
in the bulk. We consider $n^{{\mathtt{q}}-1}$ replica copies of the orbifold
$\widetilde{\mathcal{B}}_{n}$ and insert the bulk twist operators
${\mathcal{V}}_{g}$ at $\widetilde{\mathcal{L}}_{g}$. This has the effect of
reversing the orbifold and give back the geometry ${\mathcal{B}}_{n}$. The
partition function on the replicated manifold is then simply the partition
function ${\mathcal{Z}}_{n}^{({\mathtt{q}})}$. This is the partition function
that goes into the calculation of the multi-entropy. If the replica trick was
performed on the bulk geometry ${\mathcal{B}}$, it would have given the bulk
multi-entropy $S^{({\mathtt{q}})}_{\rm bulk}({\mathcal{W}})$ corresponding to
the chamber decomposition ${\mathcal{C}}_{a}$ directly, however, because the
replica trick was performed on $\widetilde{\mathcal{B}}_{n}$, it is not
obvious that what we get is $S^{({\mathtt{q}})}_{\rm bulk}({\mathcal{W}})$.
This situation is similar to the bi-partite case [26] (see also [27]). There,
the difference between the two quantities is captured by changing the
classical solution by ${\mathcal{O}}(G_{N})$ to account for the one loop
expectation value of the stress tensor. This changes the area and hence the
entanglement entropy by ${\mathcal{O}}(1)$. We expect a similar formula to
give the sub-leading correction to the multi-entropy

$\displaystyle
S^{({\mathtt{q}})}=\frac{\langle\hat{A}({\mathcal{W}})\rangle}{4G_{N}}+S_{\rm
bulk}^{({\mathtt{q}})}({\mathcal{W}})+{\rm c.t.}.$ (31)

Here $\hat{A}({\mathcal{W}})$ is the area operator of the soap film and ${\rm
c.t.}$ are the counter-terms that render $S_{\rm
bulk}^{({\mathtt{q}})}({\mathcal{W}})$ finite.

Following [6], [28] we conjecture a formula that is valid to all orders in
$1/G_{N}$ perturbation theory: Multi-entropy is given by the above formula but
${\mathcal{W}}$ is not the ordinary area extremizing soap-film but rather the
“quantum extremal soap-film” i.e. the soap film that extremizes the
combination

$\displaystyle\frac{A({\mathcal{W}})}{4G_{N}}+S_{\rm
bulk}^{({\mathtt{q}})}({\mathcal{W}}).$ (32)

In analogy with the bi-partite case, we call this prescription the quantum
extremal soap-film prescription.

## Acknowledgements

We would like to thank Akshay Gadde, Gautam Mandal, Shiraz Minwalla, Piyush
Shrivastava, Pushkal Shrivastava and Sandip Trivedi for interesting
discussions. We are particularly indebted to Onkar Parrikar for stimulating
discussions and insightful comments. We would also like to thank Shraiyance
Jain for collaboration on related projects. This work is supported by the
Infosys Endowment for the study of the Quantum Structure of Spacetime and by
the SERB Ramanujan fellowship. We acknowledge support of the Department of
Atomic Energy, Government of India, under Project Identification No. RTI 4002.
We would also like to acknowledge our debt to the people of India for their
steady support to the study of the basic sciences.

## References

  * Calabrese and Cardy [2004] P. Calabrese and J. L. Cardy, Entanglement entropy and quantum field theory, J. Stat. Mech. 0406, P06002 (2004), arXiv:hep-th/0405152 .
  * Ryu and Takayanagi [2006] S. Ryu and T. Takayanagi, Holographic derivation of entanglement entropy from AdS/CFT, Phys. Rev. Lett. 96, 181602 (2006), arXiv:hep-th/0603001 .
  * Hubeny  _et al._ [2007] V. E. Hubeny, M. Rangamani, and T. Takayanagi, A Covariant holographic entanglement entropy proposal, JHEP 07, 062, arXiv:0705.0016 [hep-th] .
  * Casini  _et al._ [2011] H. Casini, M. Huerta, and R. C. Myers, Towards a derivation of holographic entanglement entropy, JHEP 05, 036, arXiv:1102.0440 [hep-th] .
  * Dong [2014] X. Dong, Holographic Entanglement Entropy for General Higher Derivative Gravity, JHEP 01, 044, arXiv:1310.5713 [hep-th] .
  * Engelhardt and Wall [2015] N. Engelhardt and A. C. Wall, Quantum Extremal Surfaces: Holographic Entanglement Entropy beyond the Classical Regime, JHEP 01, 073, arXiv:1408.3203 [hep-th] .
  * Bekenstein [1973] J. D. Bekenstein, Black holes and entropy, Phys. Rev. D 7, 2333 (1973).
  * Hawking [1975] S. W. Hawking, Particle Creation by Black Holes, Commun. Math. Phys. 43, 199 (1975).
  * Dutta and Faulkner [2021] S. Dutta and T. Faulkner, A canonical purification for the entanglement wedge cross-section, JHEP 03, 178, arXiv:1905.00577 [hep-th] .
  * Calabrese  _et al._ [2012] P. Calabrese, J. Cardy, and E. Tonni, Entanglement negativity in quantum field theory, Phys. Rev. Lett. 109, 130502 (2012), arXiv:1206.3092 [cond-mat.stat-mech] .
  * Kusuki  _et al._ [2019] Y. Kusuki, J. Kudler-Flam, and S. Ryu, Derivation of holographic negativity in AdS3/CFT2, Phys. Rev. Lett. 123, 131603 (2019), arXiv:1907.07824 [hep-th] .
  * Penington  _et al._ [2022] G. Penington, M. Walter, and F. Witteveen, Fun with replicas: tripartitions in tensor networks and gravity, (2022), arXiv:2211.16045 [hep-th] .
  * Regge and Viano [1962] T. Regge and G. A. Viano, The interpolation problem in the theory of complex angular momentum, Nuovo Cimento 10.1007/BF02733142 (1962).
  * Albino  _et al._ [1964] E. Albino, M. Bertero, and G. A. Viano, An analyticity test in the theory of complex angular momentum, Nuovo Cimento 10.1007/BF02726068 (1964).
  * Note [1] We would like to thank Geoff Penington to bring this paper to our attention.
  * Nielsen and Kempe [2001] M. A. Nielsen and J. Kempe, Separable states are more disordered globally than locally, Physical Review Letters 86, 5184 (2001).
  * Lewkowycz and Maldacena [2013] A. Lewkowycz and J. Maldacena, Generalized gravitational entropy, JHEP 08, 090, arXiv:1304.4926 [hep-th] .
  * Dong [2016] X. Dong, The Gravity Dual of Renyi Entropy, Nature Commun. 7, 12472 (2016), arXiv:1601.06788 [hep-th] .
  * Unruh  _et al._ [1989] W. G. Unruh, G. Hayward, W. Israel, and D. Mcmanus, Cosmic string loops are straight, Phys. Rev. Lett. 62, 2897 (1989).
  * Boisseau  _et al._ [1997] B. Boisseau, C. Charmousis, and B. Linet, Dynamics of a selfgravitating thin cosmic string, Phys. Rev. D 55, 616 (1997), arXiv:gr-qc/9607029 .
  * Note [2] Another length scale present in the problem is the AdS scale $\ell_{\rm AdS}$, however it is irrelevant to this completely local computation.
  * Headrick [2010] M. Headrick, Entanglement Renyi entropies in holographic theories, Phys. Rev. D 82, 126010 (2010), arXiv:1006.0047 [hep-th] .
  * Hung  _et al._ [2011] L.-Y. Hung, R. C. Myers, M. Smolkin, and A. Yale, Holographic Calculations of Renyi Entropy, JHEP 12, 047, arXiv:1110.1084 [hep-th] .
  * [24] A. Gadde, V. Krishna, and T. Sharma, Multi-partite entanglement measures and their holographic duals, To appear .
  * Wall [2014] A. C. Wall, Maximin surfaces, and the strong subadditivity of the covariant holographic entanglement entropy, Classical and Quantum Gravity 31, 225007 (2014).
  * Faulkner  _et al._ [2013] T. Faulkner, A. Lewkowycz, and J. Maldacena, Quantum corrections to holographic entanglement entropy, JHEP 11, 074, arXiv:1307.2892 [hep-th] .
  * Barrella  _et al._ [2013] T. Barrella, X. Dong, S. A. Hartnoll, and V. L. Martin, Holographic entanglement beyond classical gravity, JHEP 09, 109, arXiv:1306.4682 [hep-th] .
  * Dong and Lewkowycz [2018] X. Dong and A. Lewkowycz, Entropy, Extremality, Euclidean Variations, and the Equations of Motion, JHEP 01, 081, arXiv:1705.08453 [hep-th] .


# Positivity-preserving non-Markovian Master Equation for Open Quantum System
Dynamics: Stochastic Schrödinger Equation Approach

Wufu Shi1,3<EMAIL_ADDRESS>Yusui Chen2<EMAIL_ADDRESS>Quanzhen
Ding3 Jin Wang1<EMAIL_ADDRESS>Ting Yu3<EMAIL_ADDRESS>1Department of Chemistry, Stony Brook University, Stony Brook, New York 11794,
USA  2Physics Department, New York Institute of Technology, Old Westbury, NY
11568, USA  3Physics Department, Stevens Institute of Technology, Hoboken, NJ
07030, USA

###### Abstract

Positivity preservation is naturally guaranteed in exact non-Markovian master
equations for open quantum system dynamics. However, in many approximated non-
Markovian master equations, the positivity of the reduced density matrix is
not guaranteed. In this paper, we provide a general class of time-local,
perturbative and positivity-preserving non-Markovian master equations
generated from stochastic Schrödinger equations, particularly quantum-state-
diffusion equations. Our method has an expanded range of applicability for
accommodating a variety of non-Markovian environments. We show the positivity-
preserving master equation for a three-level system coupled to a dissipative
bosonic environment as a particular example to exemplify our general approach.
We illustrate the numerical simulations with an analysis explaining why the
previous approximated non-Markovian master equations cannot guarantee
positivity. Our work provides a consistent master equation for studying the
non-Markovian dynamics in ultrafast quantum processes and strong-coupling
systems.

Open quantum system, non-Markovian dynamics, positivity-preserving, master
equation

††preprint: APS/123-QED

##  I Introduction

A density matrix of a quantum system is positive semidefinite, as its
eigenvalues naturally are the probabilities of the associated eigenstates. For
a closed system, the positivity of the density matrix is always preserved in
the dynamical equations, e.g., the von Neumann equation. However, no quantum
systems can be isolated from the surrounding environment. In the context of
the theory of open quantum systems (OQSs), the state of the central quantum
system is characterized by the reduced density matrix whose time evolution
equation is the master equation (ME) instead [1, 2]. Generally, it is
extremely difficult to obtain the exact ME due to the infinite number of
degrees of freedom of the environment. Several perturbation strategies have
been developed in the past decades to obtain approximated MEs [3, 4]. For
instance, Lindblad-type [5] and Redfield-type [6] MEs based on the Born-Markov
approximation effectively describe the Markovian dynamics of many physical
processes [7, 8, 9]. However, among the two typical Markovian MEs, the former
can preserve positivity, while the latter cannot [10]. It is a dilemma to
preserve the positivity of MEs with perturbative methods beyond the original
Lindblad MEs.

Moreover, the theory of non-Markovian OQSs recently attracted great interests
because Markovian approximations are not valid in certain ultrafast processes
[11, 12, 13, 14, 15, 16]. A comprised solution is to use various _weaker_
approximations to maintain certain non-Markovian features beyond the Lindblad
ME. Usually, such changes would lead to the new ME which cannot guarantee
positivity preservation. Another feasible solution is the hierarchical
equations of motion (HEOM) technique [17, 18], which is a numerically exact
approach to studying the evolution of a density matrix without the typical
assumptions that conventional Lindblad or Redfield MEs use. HEOM technique is
applicable in computing expectation values of quantum observables, at both
zero and finite temperature. But HEOM is not a conventional ME, which is
supposed to be a homogeneous equation of the reduced density matrix only. In
studying the detailed balance breaking in open quantum systems [19], the
numerically generated density matrix, in chronological order,
$\hat{\rho}_{r}(t)$ is often insufficient to compute the probability flow and
analyze the flow’s detailed components, while a conventional ME can interpret
transition processes between arbitrary two states explicitly. As the result,
it is crucial and necessary to obtain a self-consistent ME. And the consequent
challenge is twofold: (1) obtaining exact non-Markovian MEs is extremely
difficult because of the lack of mathematical tools [20, 21, 22, 23, 24, 25,
26, 27, 28, 29, 30, 31, 32, 33]; (2) positivity preservation is not guaranteed
in perturbative MEs when certain approximations are applied [34]. The purpose
of this paper is to solve this long-standing problem.

We demonstrate our solution in Fig. 1. Due to the failure to guarantee
positivity, the path of obtaining the consistent approximated ME from the
exact ME is blocked. However, in contrast to the dynamics of mixed states of
OQSs, the positivity preservation is always guaranteed in the pure state
dynamics, even when approximations are applied ($|\psi\rangle\langle\psi|$ is
always positive-semidefinite, here the state $|\psi\rangle$ does not have to
be normalized). Additionally, the stochastic Schrödinger equation (SSE)
approach has provided a rigorous method to obtain the associated ME by
averaging an infinite number of stochastic pure-state trajectories in the
Markov limit [35, 36, 37, 38, 39, 40, 41]. Thus, we propose a strategy to
generate positivity-preserving non-Markovian MEs: (1) start with a formal
exact non-Markovian SSE; (2) apply truncation and obtain the approximated SSE;
(3) recover the approximated ME from the approximated SSE rigorously. In the
end, we assure that the generated approximated ME using the strategy can
guarantee positivity as expected.

$\textstyle{\rm{ME}\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces}$$\scriptstyle{unraveling}$$\scriptstyle{approximation}$$\textstyle{\rm{SSE}\ignorespaces\ignorespaces\ignorespaces\ignorespaces}$$\scriptstyle{approximation}$$\textstyle{\times\ignorespaces\ignorespaces\ignorespaces\ignorespaces}$$\begin{array}[]{c}negative\\\
probability\end{array}$$\textstyle{\rm{ME_{\rm app}}}$$\textstyle{\rm{SSE_{\rm
app}}\ignorespaces\ignorespaces\ignorespaces\ignorespaces}$$\scriptstyle{ensemble\;\;average}$

Figure 1: The mind map of obtaining positivity-preserving MEs.

If we restrict the measurement on the environment to Bargmann coherent states,
the evolution of the stochastic pure state can be characterized by the
quantum-state-diffusion (QSD) equation [42, 43, 44, 45, 46, 47, 48], a time-
evolution equation of a stochastic quantum trajectory
$|\psi_{z}\rangle\equiv\langle z||\Psi_{{\rm tot}}\rangle$, where $||z\rangle$
is the Bargmann coherent state of the environment and $z$ represents a large
number of complex Gaussian random variables. Consequently, the reduced density
matrix can be recovered using $\hat{\rho}_{r}={\rm Tr}_{\rm E}[|\Psi_{{\rm
tot}}\rangle\langle\Psi_{{\rm tot}}|]$, which is equivalent to the ensemble
average over all stochastic quantum trajectories
$\hat{\rho}_{r}={\mathcal{M}}(|\psi_{z}\rangle\langle\psi_{z}|)$. As shown in
Fig.1, the exact ME and SSE are rigorously equivalent. When some
approximations have to be employed, the approximated SSE ($\rm{SSE_{\rm
app}}$) can numerically generate the positivity-preserving reduced density
matrix $\hat{\rho}_{r}(t)$ at any time. But the approximated ME, if derived in
the same manner, is not positivity-preserving guaranteed. In this paper, we
develop a method to generate an approximated ME ($\rm{ME_{\rm app}}$), which
is positivity-preserving and has the same reduced density matrix as the ones
numerically recovered from the corresponding approximated SSE at any time.

This paper is organized as follows. In Sec. II, we briefly review the QSD
approach and introduce our method to derive the positivity-preserving MEs. In
Sec. III, we study a dissipative three-level system and demonstrate how to
derive the positivity-preserving ME explicitly. We close with a conclusion in
Sec. IV.

##  II General Methods

The theory of OQSs studies the dynamics of a quantum system coupled with an
external quantum system or an environment. Generally, the system’s dynamics
are significantly influenced by the environment, e.g., the quantum decoherence
process and the quantum entanglement regeneration process. The total
Hamiltonian of the combined system and environment is usually written as:

$\hat{H}_{{\rm tot}}=\hat{H}_{\rm S}+\hat{H}_{{\rm int}}+\hat{H}_{\rm E}.$ (1)

Here, the environment Hamiltonian $\hat{H}_{\rm E}$ contains an infinite
number of bosonic modes. And the interaction Hamiltonian of the coupling
between system and environment can be assumed and formally written as:

$\displaystyle\hat{H}_{\rm E}$ $\displaystyle=$
$\displaystyle\sum_{k}\omega_{k}\hat{b}_{k}^{\dagger}\hat{b}_{k},$
$\displaystyle\hat{H}_{{\rm int}}$ $\displaystyle=$
$\displaystyle\hat{L}\sum_{k}g_{k}\hat{b}_{k}^{\dagger}+{\rm H.c.},$

where $\hat{L}$ is the system’s operator linearly coupled to the environment.
We assume that the environment is at zero temperature and the initial state of
the combined system and environment is a product state, $|\Psi_{{\rm
tot}}(t=0)\rangle=|\psi_{\rm S}(t=0)\rangle\otimes|0_{\rm E}\rangle$. As
mentioned above, having restricted the measurement on the environment to
Bargmann coherent states will lead to the quantum trajectory in the form of
$|\psi_{z}\rangle\equiv\langle z||\Psi_{{\rm tot}}\rangle$, where
$||z\rangle\equiv\otimes_{k}||z_{k}\rangle$ is the Bargmann coherent state of
the entire environment, satisfying $\hat{b}_{k}||z\rangle=z_{k}||z\rangle$.
The evolution of $|\psi_{z}\rangle$ is governed by the QSD equation.

###  II.1 Quantum-state-diffusion approach

In the environmental interaction picture, the interaction Hamiltonian reads:

$\displaystyle\hat{H}_{{\rm int}}^{\rm
I}(t)=\hat{L}\sum_{k}g_{k}\hat{b}_{k}^{\dagger}e^{i\omega_{k}t}+\hat{L}^{\dagger}\sum_{k}g_{k}^{*}\hat{b}_{k}e^{-i\omega_{k}t}.$
(2)

Using the identity resolution of Bargmann coherent states

$\displaystyle\hat{I}_{\rm
E}=\int\frac{d^{2}z}{\pi}e^{-|z|^{2}}||z\rangle\langle z||,$ (3)

the Schrödinger equation, regarding the evolution of the pure state of the
composite system $|\Psi_{{\rm tot}}\rangle$, can be rewritten in the Bargmann
space representation (setting $\hbar=1$),

$\displaystyle\partial_{t}\langle z||\Psi_{{\rm tot}}\rangle$ $\displaystyle=$
$\displaystyle-i\langle z||(\hat{H}_{\rm S}+\hat{H}_{{\rm int}}^{\rm
I}(t))|\Psi_{{\rm tot}}\rangle$ (4) $\displaystyle=$
$\displaystyle-i(\hat{H}_{\rm
S}+\hat{L}\sum_{k}g_{k}z_{k}^{*}e^{i\omega_{k}t})\langle z||\Psi_{{\rm
tot}}\rangle$
$\displaystyle-i\hat{L}^{\dagger}\sum_{k}g_{k}^{*}e^{-i\omega_{k}t}\frac{\partial}{\partial
z_{k}^{*}}\langle z||\Psi_{{\rm tot}}\rangle.$

Here, we define a complex Gaussian process
$z_{t}^{*}\equiv-i\sum_{k}g_{k}z_{k}^{*}e^{i\omega_{k}t}$, which satisfies the
following relations: ${\mathcal{M}}(z_{t})={\mathcal{M}}(z_{t}z_{s})=0$, and
$\alpha(t,s)\equiv{\mathcal{M}}(z_{t}z_{s}^{*})=\sum_{k}|g_{k}|^{2}e^{-i\omega_{k}(t-s)}$,
where $\alpha(t,s)$ is the correlation function of the complex Gaussian
process $z_{t}^{*}$. Then, using the chain rule
$\frac{\partial(\cdot)}{\partial z_{k}^{*}}=\int_{0}^{t}ds\frac{\partial
z_{s}^{*}}{\partial z_{k}^{*}}\frac{\delta(\cdot)}{\delta z_{s}^{*}}$, Eq. (4)
can lead to a formal linear non-Markovian QSD equation:

$\displaystyle\partial_{t}|\psi_{z}\rangle=(-i\hat{H}_{\rm
S}+z_{t}^{*}\hat{L}-\hat{L}^{\dagger}\int_{0}^{t}ds\alpha(t,s)\delta_{z_{s}^{*}})|\psi_{z}\rangle.$
(5)

By taking the statistical mean over all trajectories, the reduced density
matrix of the system can be recovered

$\displaystyle\hat{\rho}_{r}={\mathcal{M}}(|\psi_{z}\rangle\langle\psi_{z}|).$
(6)

The functional derivative term in Eq. (5) can be formally written using a to-
be-determined operator $\hat{O}$, defined as

$\displaystyle\hat{O}(t,s)|\psi_{z}\rangle\equiv\delta_{z^{*}_{s}}|\psi_{z}\rangle,$
(7)

which can be solved through an operator evolution equation:

$\displaystyle\partial_{t}\hat{O}(t,s)=-i[\hat{H}_{\rm
eff},\hat{O}(t,s)]-i\delta_{z^{*}_{s}}\hat{H}_{\rm eff},$ (8)

where $\hat{H}_{\rm eff}$ is the effective Hamiltonian:

$\displaystyle\hat{H}_{\rm eff}\equiv\hat{H}_{\rm
S}+iz_{t}^{*}\hat{L}-i\hat{L}^{\dagger}\bar{O},$ (9)

and $\bar{O}(t)\equiv\int_{0}^{t}ds\,\alpha(t,s)\hat{O}(t,s)$. Therefore, the
formal linear QSD equation reads:

$\displaystyle\partial_{t}|\psi_{z}\rangle=(-i\hat{H}_{\rm
S}+z_{t}^{*}\hat{L}-\hat{L}^{\dagger}\bar{O}(t))|\psi_{z}\rangle.$ (10)

Generally, the structure of the exact O-operator is complicated. Only a few
models can be solved with the exact O-operator. A compromised solution to this
difficulty is to replace the exact O-operator with an approximated one. One
can drop certain terms of the O-operator to simplify the calculation, called a
truncation operation. How to truncate the O-operator depends on the type of
interaction and the size of the system. Without the loss of generality, the
O-operator can be written as a sum of all component operators [49]:

$\displaystyle\hat{O}(t,s,z^{*})=\sum_{n}\hat{O}_{n}(t,s,z^{*}).$ (11)

And the approximated O-operator after the truncation is defined as:

$\displaystyle\hat{O}^{N}(t,s,z^{*})\equiv\sum_{n=0}^{N}\hat{O}_{n}(t,s,z^{*}).$
(12)

Consequently, the approximated QSD equation after the truncation reads:

$\displaystyle\partial_{t}|\psi^{N}_{z}\rangle=(-i\hat{H}_{\rm
S}+z_{t}^{*}\hat{L}-\hat{L}^{\dagger}\int_{0}^{t}ds\alpha(t,s)\hat{O}^{N})|\psi^{N}_{z}\rangle.$
(13)

where the trajectory $|\psi_{z}^{N}\rangle$ is the associated approximated
trajectory.

One of the advantages of the non-Markovian QSD approach is that any reduced
density operator recovered from quantum trajectories
$\hat{\rho}_{r}={\mathcal{M}}(|\psi_{z}\rangle\langle\psi_{z}|)$ is always
positivity preserved, even if quantum trajectories are numerically generated
by the approximated QSD equation (13),
$\hat{\rho}_{r}^{N}={\mathcal{M}}(|\psi^{N}_{z}\rangle\langle\psi^{N}_{z}|)$.
(For the single trajectory, we know that
$|\psi_{z}^{N}\rangle\langle\psi_{z}^{N}|$ must be positive semidefinite. The
ensemble average ${\mathcal{M}}(|\psi_{z}^{N}\rangle\langle\psi_{z}^{N}|)$ can
be considered as a convex combination of
$|\psi_{z}^{N}\rangle\langle\psi_{z}^{N}|$, therefore,
${\mathcal{M}}(|\psi_{z}^{N}\rangle\langle\psi_{z}^{N}|)$ is also positive
semidefinite. Principally, this is how we derive the positivity-preserving ME
from the approximated QSD equation.

###  II.2 Positivity-preserving ME

For a given exact QSD equation, the associated ME reads:

$\displaystyle\partial_{t}\hat{\rho}_{r}$ $\displaystyle=$
$\displaystyle{\mathcal{M}}(\frac{\partial|\psi_{z}\rangle}{\partial
t}\langle\psi_{z}|+|\psi_{z}\rangle\frac{\partial\langle\psi_{z}|}{\partial
t})$ (14) $\displaystyle=$ $\displaystyle{\mathcal{M}}(-iH_{\rm
eff}|\psi_{z}\rangle\langle\psi_{z}|+i|\psi_{z}\rangle\langle\psi_{z}|H_{\rm
eff}^{\dagger})$ $\displaystyle=$ $\displaystyle-i[\hat{H}_{\rm
S},\hat{\rho}_{r}]+\hat{L}{\mathcal{M}}(z_{t}^{*}\hat{P})+{\mathcal{M}}(z_{t}\hat{P})\hat{L}^{\dagger}$
$\displaystyle-\hat{L}^{\dagger}{\mathcal{M}}(\bar{O}\hat{P})-{\mathcal{M}}(\hat{P}\bar{O}^{\dagger})\hat{L},$

where $\hat{P}$ denotes the stochastic operator
$\hat{P}\equiv|\psi_{z}\rangle\langle\psi_{z}|$. Using the Novikov theorem
(see Appendix A), the two terms, ${\mathcal{M}}(z_{t}^{*}\hat{P})$ and
${\mathcal{M}}(z_{t}\hat{P})$, in the above equation, can be estimated,

$\displaystyle{\mathcal{M}}(z_{t}^{*}\hat{P})=\int_{0}^{t}ds\alpha^{*}(t,s){\mathcal{M}}(\delta_{z_{s}}\hat{P})={\mathcal{M}}(\hat{P}\bar{O}^{\dagger}).$
(15)

As a result, the formal ME is obtained,

$\displaystyle\partial_{t}\hat{\rho}_{r}=-i[\hat{H}_{\rm
S},\hat{\rho}_{r}]+[\hat{L},{\mathcal{M}}(\hat{P}\bar{O}^{\dagger})]-[\hat{L}^{\dagger},{\mathcal{M}}(\bar{O}\hat{P})].$
(16)

The above-derived ME is positivity-preserving since the reduced density matrix
$\hat{\rho}_{r}$ is rigorously equivalent to the exact stochastic quantum
trajectory governed by the QSD equation (5).

Next, we will demonstrate why the ME can not guarantee positivity if all the
four exact O-operators in Eq. (16) are replaced by the approximated one
$\hat{O}^{N}$. Following the similar method of obtaining Eq. (16), the
approximated ME for the perturbative reduced density matrix
$\hat{\rho}_{r}^{\prime}$ reads:

$\displaystyle\partial_{t}\hat{\rho}_{r}^{\prime}=-i[\hat{H}_{\rm
S},\hat{\rho}_{r}^{\prime}]+[\hat{L},{\mathcal{M}}(\hat{P}^{\prime}(\bar{O}^{N\dagger})]-[\hat{L}^{\dagger},{\mathcal{M}}(\bar{O}^{N}\hat{P}^{\prime})].$
(17)

However, it is worth pointing out that the reduced density matrix
$\hat{\rho}_{r}^{\prime}$ can violate positivity because the approximated ME
can not be unraveled by the QSD equation (13). To clarify the difference
between Eq. (17) and the approximated ME which is rigorously equivalent to the
approximated QSD equation (13), we recover the ME starting from the identity
$\hat{\rho}_{r}^{N}={\mathcal{M}}(|\psi_{z}^{N}\rangle\langle\psi_{z}^{N}|)$
and the QSD equation (13). The approximated ME reads

$\displaystyle\partial_{t}\hat{\rho}_{r}^{N}$ $\displaystyle=$
$\displaystyle-i[\hat{H}_{\rm
S},\hat{\rho}_{r}^{N}]+\hat{L}{\mathcal{M}}(z_{t}^{*}\hat{P}^{N})+{\mathcal{M}}(z_{t}\hat{P}^{N})\hat{L}^{\dagger}$
(18)
$\displaystyle-\hat{L}^{\dagger}{\mathcal{M}}(\bar{O}^{N}\hat{P}^{N})-{\mathcal{M}}(\hat{P}^{N}\bar{O}^{N}{}^{\dagger})\hat{L},$

where $\hat{P}^{N}\equiv|\psi_{z}^{N}\rangle\langle\psi_{z}^{N}|$. After
applying the Novikov theorem to simplify the term of
${\mathcal{M}}(z_{t}^{*}\hat{P}^{N})$, it is easy to verify that in the
general case

$\displaystyle{\mathcal{M}}(z_{t}^{*}\hat{P}^{N})=\int_{0}^{t}ds\,\alpha^{*}(t,s){\mathcal{M}}(\delta_{z_{s}}\hat{P}^{N})\neq{\mathcal{M}}(\hat{P}^{N}\bar{O}^{N\dagger}).$

This is why the above-mentioned approximated ME (17) cannot be unraveled by
the QSD equation (13). To solve this problem, we need to know the exact value
of $\delta_{z_{s}}\hat{P}^{N}$, therefore a new O-operator has to be
introduced

$\displaystyle\hat{O}_{d}(t,s,z^{*})|\psi_{z}^{N}\rangle\equiv\delta_{z_{s}^{*}}|\psi_{z}^{N}\rangle,$
(19)

where the new operator $\hat{O}_{d}(t,s,z^{*})$ is determined by an evolution
equation,

$\displaystyle\partial_{t}\hat{O}_{d}(t,s,z^{*})$ $\displaystyle=$
$\displaystyle[-i\hat{H}_{\rm
S}+z_{t}^{*}\hat{L}-\hat{L}^{\dagger}\bar{O}^{N},\hat{O}_{d}(t,s,z^{*})]$ (20)
$\displaystyle-\hat{L}^{\dagger}\delta_{z_{s}^{*}}\bar{O}^{N},$

together with the initial condition,

$\displaystyle\hat{O}_{d}(t=s,s,z^{*})=\hat{L}.$ (21)

The subtle difference between $\hat{O}^{N}$ and $\hat{O}_{d}$ is just the
reason of positivity violation in the ME (17). Consequently, the result of
applying the Novikov theorem is revised to

$\displaystyle{\mathcal{M}}(z_{t}^{*}\hat{P}^{N})$ $\displaystyle=$
$\displaystyle\int_{0}^{t}ds\,\alpha^{*}(t,s){\mathcal{M}}(\delta_{z_{s}}\hat{P}^{N})$
(22) $\displaystyle=$
$\displaystyle{\mathcal{M}}(\hat{P}^{N}\bar{O}_{d}^{\dagger}),$

where
$\bar{O}_{d}^{\dagger}(t,z^{*})\equiv\int_{0}^{t}ds\alpha^{*}(t,s)\hat{O}_{d}^{\dagger}(t,s,z^{*})$.
By substituting Eq. (22) into Eq. (18), we obtain the formal positivity-
preserving approximated ME,

$\displaystyle\partial_{t}\hat{\rho}_{r}^{N}$ $\displaystyle=$
$\displaystyle-i[\hat{H}_{\rm
S},\hat{\rho}_{r}^{N}]+\hat{L}{\mathcal{M}}(\hat{P}^{N}\bar{O}_{d}^{\dagger})-{\mathcal{M}}(\hat{P}^{N}\bar{O}^{N\dagger})\hat{L}$
(23)
$\displaystyle-\hat{L}^{\dagger}{\mathcal{M}}(\bar{O}^{N}\hat{P}^{N})+{\mathcal{M}}(\bar{O}_{d}\hat{P}^{N})\hat{L}^{\dagger}.$

##  III Models and Results

Here, we consider a ladder-type three-level system coupled with a dissipative
zero-temperature reservoir and use it to demonstrate how to obtain the
positivity-preserving approximated ME. The total Hamiltonian reads

$\displaystyle\hat{H}_{\rm
tot}=\omega\hat{J}_{z}+\sum_{k}g_{k}(\hat{J}_{+}\hat{b}_{k}+\hat{b}_{k}^{\dagger}\hat{J}_{-})+\sum_{k}\omega_{k}\hat{b}_{k}^{\dagger}\hat{b}_{k},$
(24)

where $g_{k}$ is the real coupling strength of the $k$th mode. $\hat{J}_{+}$
($\hat{J}_{-}$) is the raising (lowering) operator of the three-level system,
satisfying the commutation relation
$\hat{J}_{z}=\frac{1}{2}[\hat{J}_{+},\hat{J}_{-}]$. The three operators have
the matrix form

$\displaystyle\hat{J}_{z}=\begin{bmatrix}1&0&0\\\ 0&0&0\\\
0&0&-1\end{bmatrix},\hat{J}_{+}=\sqrt{2}\begin{bmatrix}0&1&0\\\ 0&0&1\\\
0&0&0\end{bmatrix},\hat{J}_{-}=\sqrt{2}\begin{bmatrix}0&0&0\\\ 1&0&0\\\
0&1&0\end{bmatrix}.$

###  III.1 Master equation for the three-level system

In Ref. [50, 26], it has been proved that the O-operator of a dissipative
three-level system contains noise up to the first order. We use a noise-free
operator $\hat{O}^{\\{0\\}}$ to replace the exact $\hat{O}$ in the QSD
equation:

$\displaystyle\partial_{t}|\psi_{z}^{\\{0\\}}\rangle=(-i\hat{H}_{\rm
S}+z^{*}_{t}\hat{L}-\hat{L}^{\dagger}\bar{O}^{\\{0\\}})|\psi_{z}^{\\{0\\}}\rangle,$
(25)

where
$\bar{O}^{\\{0\\}}\equiv\int_{0}^{t}ds\alpha(t,s)\hat{O}^{\\{0\\}}(t,s)$, and
the Lindblad operator $\hat{L}=\hat{J}_{-}$. The operator $\hat{O}^{\\{0\\}}$
is governed by its evolution equation

$\displaystyle\partial_{t}\hat{O}^{\\{0\\}}(t,s)=[-i\hat{H}_{\rm
S}-\hat{L}^{\dagger}\bar{O}^{\\{0\\}},\hat{O}^{\\{0\\}}(t,s)],$ (26)

with its initial condition

$\displaystyle\quad\quad\hat{O}^{\\{0\\}}(t=s,s)=\hat{L}.$ (27)

According to Eq. (26), the equation is valid only when the operator
$\hat{O}^{\\{0\\}}$ has the form of

$\displaystyle\hat{O}^{\\{0\\}}(t,s)\equiv
f_{1}(t,s)\hat{J}_{-}+g_{1}(t,s)\hat{J}_{z}\hat{J}_{-},$ (28)

where $f_{1}(t,s)$ and $g_{1}(t,s)$ are two to-be-determined evolution
coefficients. By substituting the ansatz (28) into Eq. (26), the coefficients
$f_{1}$ and $g_{1}$ can be determined by the following differential equations

$\displaystyle\partial_{t}f_{1}$ $\displaystyle=$
$\displaystyle(i\omega+2G_{1})f_{1},$ $\displaystyle\partial_{t}g_{1}$
$\displaystyle=$
$\displaystyle(-2F_{1}+4G_{1})f_{1}+(i\omega+2F_{1}-2G_{1})g_{1},\;\;$ (29)

associated with the initial conditions,

$\displaystyle f_{1}(t=s,s)=1,\quad\quad g_{1}(t=s,s)=0,$ (30)

where $F_{1}(t)\equiv\int_{0}^{t}ds\alpha(t,s)f_{1}(t,s)$, and
$G_{1}(t)\equiv\int_{0}^{t}ds\alpha(t,s)g_{1}(t,s)$. Subsequently, the time-
evolution equation of the operator $\hat{O}_{d}(t,s,z^{*})$ reads

$\displaystyle\partial_{t}\hat{O}_{d}$ $\displaystyle=$
$\displaystyle[-i\hat{H}_{\rm
S}+z_{t}^{*}\hat{L}-\hat{L}^{\dagger}\bar{O}^{\\{0\\}},\hat{O}_{d}].$ (31)

Note that the last functional derivative term in Eq. (20) has been eliminated
because $\hat{O}^{\\{0\\}}$ is noise-free. Similarly, the ansatz of the
operator $\hat{O}_{d}$ reads

$\displaystyle\hat{O}_{d}(t,s,{\bm{z}}^{*})$ $\displaystyle\equiv$
$\displaystyle f_{2}(t,s)\hat{J}_{-}+g_{2}(t,s)\hat{J}_{z}\hat{J}_{-}$ (32)
$\displaystyle+\int_{0}^{t}ds^{\prime}p_{2}(t,s,s^{\prime})z^{*}_{s^{\prime}}\hat{J}_{-}^{2}.$

By substituting the ansatz (32) into Eq. (31), the new set of coefficients,
$f_{2}(t,s),g_{2}(t,s)$ and $p_{2}(t,s,s^{\prime})$, are determined by

$\displaystyle\partial_{t}f_{2}$ $\displaystyle=$
$\displaystyle(i\omega+2G_{1})f_{2},$ $\displaystyle\partial_{t}g_{2}$
$\displaystyle=$
$\displaystyle(-2F_{1}+4G_{1})f_{2}+(i\omega+2F_{1}-2G_{1})g_{2},\;\;$
$\displaystyle\partial_{t}p_{2}$ $\displaystyle=$
$\displaystyle(2i\omega+2F_{1})p_{2},$ (33)

with the initial conditions

$\displaystyle\quad f_{2}(t=s,s)$ $\displaystyle=$ $\displaystyle 1,$
$\displaystyle\quad g_{2}(t=s,s)$ $\displaystyle=$ $\displaystyle 0,$
$\displaystyle p_{2}(t=s^{\prime},s,s^{\prime})$ $\displaystyle=$
$\displaystyle g_{2}(s^{\prime},s).$

By comparing Eq. (29) and Eq. (33), it is clear that $f_{1}=f_{2}$ and
$g_{1}=g_{2}$, since they have the same evolution equations and the same
initial conditions. As a result, $f_{1}$ and $g_{1}$, in the rest of the
paper, will be replaced by $f_{2}$ and $g_{2}$, respectively.

After obtaining operators $\hat{O}_{d}$ and $\hat{O}^{\\{0\\}}$, the formal ME
(23) for the dissipative three-level model reads

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\partial_{t}\hat{\rho}_{r}^{\\{0\\}}\\!\\!\\!$
$\displaystyle=$ $\displaystyle\\!\\!\\!-i[\hat{H}_{\rm
S},\hat{\rho}_{r}^{\\{0\\}}]+\\{[(F_{2}\hat{J}_{-}+G_{2}\hat{J}_{z}\hat{J}_{-})\hat{\rho}_{r}^{\\{0\\}},\hat{J}_{+}]$
(34)
$\displaystyle\\!\\!\\!+\hat{J}_{-}^{2}\int_{0}^{t}dsP_{2}(t,s){\mathcal{M}}(z_{s}^{*}\hat{P}^{\\{0\\}})\hat{J}_{+}\\}+{\rm
H.c.},$

where $F_{2}(t)$, $G_{2}(t)$, and $P_{2}(t,s^{\prime})$ are defined as

$\displaystyle F_{2}(t)$ $\displaystyle\equiv$
$\displaystyle\int_{0}^{t}ds\alpha(t,s)f_{2}(t,s),$ $\displaystyle G_{2}(t)$
$\displaystyle\equiv$ $\displaystyle\int_{0}^{t}ds\alpha(t,s)g_{2}(t,s),$
$\displaystyle P_{2}(t,s^{\prime})$ $\displaystyle\equiv$
$\displaystyle\int_{0}^{t}ds\alpha(t,s)p_{2}(t,s,s^{\prime}).$

Applying the Novikov theorem and the termination condition in Ref. [22], the
term
$\hat{J}_{-}^{2}\int_{0}^{t}dsP_{2}(t,s){\mathcal{M}}(z_{s}^{*}\hat{P}^{\\{0\\}})\hat{J}_{+}$
in the above ME can be further simplified to

$\displaystyle\hat{J}_{-}^{2}{\mathcal{M}}(z_{s}^{*}\hat{P}^{\\{0\\}})\hat{J}_{+}$
(35) $\displaystyle=$
$\displaystyle\int_{0}^{t}ds^{\prime}\alpha^{*}(s,s^{\prime})\hat{J}_{-}^{2}{\mathcal{M}}(\hat{P}^{\\{0\\}}\hat{O}_{d}^{\dagger}(t,s^{\prime},{\bm{z}}))\hat{J}_{+}$
$\displaystyle=$
$\displaystyle\int_{0}^{t}ds^{\prime}\alpha^{*}(s,s^{\prime})f^{*}_{2}(t,s^{\prime})\hat{J}_{-}^{2}\hat{\rho}_{r}^{\\{0\\}}\hat{J}_{+}\hat{J}_{+}.$

Subsequently, the approximated positivity-preserving ME reads,

$\displaystyle\\!\\!\\!\partial_{t}\hat{\rho}_{r}^{\\{0\\}}$ $\displaystyle=$
$\displaystyle-i[\hat{H}_{\rm
S},\hat{\rho}_{r}^{\\{0\\}}]+\\{[(F_{2}\hat{J}_{-}+G_{2}\hat{J}_{z}\hat{J}_{-})\hat{\rho}_{r}^{\\{0\\}},\hat{J}_{+}]$
(36)
$\displaystyle+P_{f^{*}}\hat{J}_{-}^{2}\hat{\rho}_{r}^{\\{0\\}}\hat{J}_{+}^{2}\\}+{\rm
H.c.},$

where the coefficient
$P_{f^{*}}(t)\equiv\int_{0}^{t}\int_{0}^{t}dsds^{\prime}P_{2}(t,s)\alpha^{*}(s,s^{\prime})f^{*}_{2}(t,s^{\prime})$.
Note that $\alpha(t,s)$ is the time correlation function, corresponding to a
variety of spectra, white or colored. For simplicity, we assume the
environment is described by an Ornstein-Uhlenbeck process. So, its correlation
function is $\alpha(t,s)=a\gamma e^{-\gamma|t-s|}e^{-i\Omega(t-s)}$, where
$1/\gamma$ is the scale of memory time and $\Omega$ is the central frequency
of the environment. As a result, the coefficients’ evolution equations can be
simplified from integrodifferential equations to differential equations that

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\partial_{t}F_{2}\\!\\!\\!$
$\displaystyle=$
$\displaystyle\\!\\!\\!a\gamma+(i\omega-\gamma-i\Omega+2G_{2})F_{2},$
$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\partial_{t}G_{2}\\!\\!\\!$
$\displaystyle=$
$\displaystyle\\!\\!\\!-2F_{2}^{2}+(i\omega-\gamma-i\Omega+6F_{2}-2G_{2})G_{2},$
$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\partial_{t}\tilde{P}_{2}\\!\\!\\!$
$\displaystyle=$ $\displaystyle\\!\\!\\!a\gamma
G_{2}+(2i\omega-2\gamma-2i\Omega+2F_{2})\tilde{P}_{2},$
$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\partial_{t}P_{f^{*}}\\!\\!\\!$
$\displaystyle=$
$\displaystyle\\!\\!\\!(i\omega-\gamma-i\Omega+2F_{2}+2G^{*}_{2})P_{f^{*}}+\tilde{P}_{2}+G_{2}F_{2}^{*},$

with the initial conditions,

$\displaystyle\quad
F_{2}(t=0)=G_{2}(t=0)=\tilde{P}_{2}(t=0)=P_{f^{*}}(t=0)=0,$

where $\tilde{P}_{2}(t)\equiv\int_{0}^{t}ds\alpha(t,s)P_{2}(t,s)$.

###  III.2 Numerical results

In this section, we compare the simulation results of the population of states
of the three-level system using four different methods: (1) the exact ME,
$\hat{\rho}_{r}={\mathcal{M}}(|\psi_{z}\rangle\langle\psi_{z}|)$ in Eq. (16);
(2) the approximated positivity-preserving ME,
$\hat{\rho}_{r}^{\\{0\\}}={\mathcal{M}}(|\psi_{z}^{\\{0\\}}\rangle\langle\psi_{z}^{\\{0\\}}|)$
in Eq. (23); (3) the approximated non-positivity-preserving ME,
$\hat{\rho}_{r}^{\prime}$ in Eq. (17); (4) the approximated QSD,
$|\psi_{z}^{\\{0\\}}\rangle$ in Eq. (25).

First of all, we plot the time evolution of the population of the dissipative
three-level system, $\rho_{00}$, $\rho_{11}$, and $\rho_{22}$, generated by
approximated QSD equation approach and the approximated positivity-preserving
ME approach. The initial state of the system is prepared at an excited state:
$|\psi_{z}(t=0)\rangle=|2\rangle$. We choose the frequency, $\omega=1$, and
the central frequency of the environment $\Omega=0$. In a strong non-Markovian
regime, $\gamma=0.05$, the simulation results from two methods, as shown in
Fig. 2, overlap each other. Since the reduced density matrix generated from
the approximated QSD approach is naturally positivity-preserving, the matched
dynamics prove that our approximated ME gives rise to the same degree of
accuracy as the approximated QSD equation. It indicates that the approximated
ME can guarantee positivity.

Figure 2: Time evolution of the population of the dissipative three-level
system, $\rho_{00}$, $\rho_{11}$ and $\rho_{22}$ are generated by two
different methods: approximated linear QSD equation approach (red circle), and
the derived positivity-preserving approximated ME (black cross). The
parameters are $\omega=1$, $a=0.8$, $\gamma=0.05$, and $\Omega=0$. (The result
of the approximated QSD equation approach is obtained by averaging over 5000
quantum trajectories.)

Next, we plot the time evolution of the population of the ground state using
three different ME approaches in Fig. 3. Using the same parameters as Fig. 2,
we observe that the exact ME and our approximated positivity-preserving ME
both preserve the positivity. Meanwhile, the simulation result of the non-
positivity-preserving ME leads to failure due to the appearance of negative
probabilities in some time intervals. Furthermore, the magnitude of the
negative probability increases with time up to infinity. Consequently, the
probabilities of the other two levels also increase to infinity
simultaneously. If simply replaces the exact O-operator in the exact ME with
the truncated operator $\hat{O}^{\\{0\\}}$, then the Eq. (17) can be
explicitly written as

$\displaystyle\partial_{t}\hat{\rho}_{r}^{\prime}=-i[\hat{H}_{\rm
S},\hat{\rho}_{r}^{\prime}]+\big{(}[(F_{2}\hat{J}_{-}+G_{2}\hat{J}_{z}\hat{J}_{-})\hat{\rho}_{r}^{\prime},\hat{J}_{+}]+{\rm
H.c.}\big{)}.$

It is clear that the above approximated ME does not preserve positivity in
some parameter regions, and may lead to meaningless physics.

Figure 3:  Comparison of the time evolution of the population of the ground
state, generated by three different ME approaches: exact ME ($\rho_{r,00}$,
red circle), positivity-preserving approximated ME ($\rho^{\\{0\\}}_{r,00}$,
black cross), and the non-positivity-preserving approximated ME
($\rho^{\prime}_{r,00}$, green dash-dotted line). The parameters are
$\omega=1$, $a=0.8$, $\gamma=0.05$, and $\Omega=0$.  Figure 4:  Comparison of
the time evolution of the population of the ground state, generated by three
different ME approaches: exact ME ($\rho_{r,00}$, red circle), positivity-
preserving approximated ME ($\rho^{\\{0\\}}_{r,00}$, black cross), and the
non-positivity-preserving approximated ME ($\rho^{\prime}_{r,00}$, green dash-
dotted line), in a _moderate_ non-Markovian regime. The parameters are
$\omega=1$, $a=0.2$, $\gamma=0.2$, and $\Omega=0$.

To further demonstrate the importance of our method in studying non-Markovian
dynamics, we plot Fig. 4, the time evolution of the population of the ground
state in a _moderate_ non-Markovian regime. When $\gamma=0.2$, a shorter
memory time compared with the parameter $\gamma=0.05$ used in Fig. 3, the
dynamics simulated by the approximated non-positivity-preserving ME
$\rho_{r,00}^{\prime}$ do not contain any negative probabilities. However, its
distance from the results of the exact ME approach is significantly larger
compared with the results of our positivity-preserving ME. Comparing Figs. 3
and 4, we show that the reduced density matrix $\rho_{r}^{\\{0\\}}$ can
guarantee positivity preservation from the Markovian to the strong non-
Markovian regime. In contrast, the reduced density matrix $\rho_{r}^{\prime}$
cannot offer such confidence. Moreover, for different models and interested
parameter spaces, our method is flexible for different approximations. It
provides a robust method to obtain positivity-preserving MEs for the analysis
of non-Markovian dynamics.

##  IV Conclusion

We have addressed a long-standing issue in OQSs on how to construct
positivity-preserving approximated MEs in a general situation. Although
several developed MEs, such as Lindblad and Redfield MEs, can provide powerful
and efficient mathematical tools, these approaches have common shortcomings
since they are rooted in the Born-Markov approximation. In this study, we
start from the fact that a reduced density matrix must carry over positivity
if recovered from the ensemble average over the stochastic pure states. Then
we consider a class of linear approximated QSD equations, exploring the
possibility of constructing MEs rigorously equivalent to QSD equations. We
find that previous works mistakenly used the approximation relation
$\delta_{z_{s}^{*}}|\psi_{z}^{N}\rangle\approx\hat{O}|\psi_{z}^{N}\rangle\approx\hat{O}^{N}|\psi_{z}^{N}\rangle$.
In fact, $\delta_{z_{s}^{*}}|\psi_{z}^{N}\rangle$ should be equal to
$\hat{O}_{d}|\psi_{z}^{N}\rangle$, where $\hat{O}_{d}$ is a newly defined
operator. No matter how small the difference between the $\hat{O}_{d}$ and the
approximated $\hat{O}^{N}$ is, replacing $\hat{O}_{d}$ by $\hat{O}^{N}$ in the
derivation may lead to a violation of positivity of MEs. Consequently, it is
necessary to introduce two different approximated $\hat{O}$ to generate the
positivity-preserving approximated ME. Generally, we explain why applying an
approximation directly on the exact ME may violate positivity, while applying
the same approximation on the exact SSE will not.

In the paper, we propose a general class of positivity-preserving non-
Markovian MEs generated from SSEs, in particular, QSD equations. We explicitly
derive the approximated positivity-preserving ME for a dissipative three-level
system as a specific example of our general results. Moreover, our simulations
also show that the negative probability generated by non-positivity-preserving
MEs sometimes ends up with negative infinity, which is definitely not a
trivial issue.

In summary, we have developed a systematic method to obtain a class of
approximated but positivity-preserving non-Markovian MEs originating from
approximated linear QSD equations. With such MEs, it is feasible to study
nonequilibrium dynamics in living or biological systems, perform reliable
error analysis for quantum engineering, and investigate dynamics and phase
transitions in many-body systems.

## Acknowledgments

We thank Lajos Diósi and Water Strunz for useful conversations on the topic
over a long period of time.

##  Appendix A Novikov Theorem

Firstly, we recall several notations:

$\displaystyle{\mathcal{M}}(\cdot)$ $\displaystyle=$
$\displaystyle\int\frac{d^{2}z}{\pi}e^{-|z|^{2}}(\cdot),$ (38)

where
$\frac{d^{2}z}{\pi}\equiv\frac{d^{2}z_{1}}{\pi}\frac{d^{2}z_{2}}{\pi}...\frac{d^{2}z_{k}}{\pi}...$,
$|z|^{2}\equiv\sum_{k}|z_{k}|^{2}$, $\int
d^{2}z_{k}\equiv\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}dx_{k}dy_{k}$,
and $x_{k}\equiv{\rm Re}(z_{k})$, $y_{k}\equiv{\rm Im}(z_{k})$. In addition,
the complex Gaussian process is defined as

$\displaystyle z_{t}^{*}\equiv-i\sum_{k}g_{k}z_{k}^{*}e^{i\omega_{k}t}.$ (39)

Consequently, the left-hand side of the Eq. (15) can be explicitly expanded as

$\displaystyle{\mathcal{M}}(z_{t}^{*}\hat{P})$ $\displaystyle=$
$\displaystyle\int\frac{d^{2}z}{\pi}e^{-|z|^{2}}z_{t}^{*}\hat{P}$
$\displaystyle=$
$\displaystyle\int\frac{d^{2}z}{\pi}e^{-|z|^{2}}\big{(}-i\sum_{k}g_{k}z_{k}^{*}e^{i\omega_{k}t}\big{)}\hat{P}$
$\displaystyle=$
$\displaystyle\int\frac{d^{2}z}{\pi}e^{-|z|^{2}}\big{(}-i\sum_{k}g_{k}(x_{k}-iy_{k})e^{i\omega_{k}t}\big{)}\hat{P}.\;\;\;$

Using integration by parts for the $k$th mode, we have

$\displaystyle\int d^{2}z_{k}e^{-|z_{k}|^{2}}x_{k}\hat{P}$ $\displaystyle=$
$\displaystyle-\frac{1}{2}\int d^{2}z_{k}\frac{\partial}{\partial
x_{k}}(e^{-|z_{k}|^{2}})\hat{P}$ $\displaystyle=$
$\displaystyle-\frac{1}{2}\int
dy_{k}(e^{-|z_{k}|^{2}}\hat{P})\big{|}^{x_{k}=\infty}_{x_{k}=-\infty}+\frac{1}{2}\int
d^{2}z_{k}e^{-|z_{k}|^{2}}\partial_{x_{k}}\hat{P}.$

Usually, the boundary terms
$(e^{-|z_{k}|^{2}}\hat{P})|_{x_{k}\rightarrow\infty}$ and
$(e^{-|z_{k}|^{2}}\hat{P})|_{x_{k}\rightarrow-\infty}$ in the above row
converge to zero rapidly. Thus, we have

$\displaystyle\int d^{2}z_{k}e^{-|z_{k}|^{2}}x_{k}\hat{P}=\frac{1}{2}\int
d^{2}z_{k}e^{-|z_{k}|^{2}}\partial_{x_{k}}\hat{P}.$ (42)

Similarly, we also have

$\displaystyle\int d^{2}z_{k}e^{-|z_{k}|^{2}}y_{k}\hat{P}=\frac{1}{2}\int
d^{2}z_{k}e^{-|z_{k}|^{2}}\partial_{y_{k}}\hat{P}.$ (43)

According to the chain rule, we have

$\displaystyle\begin{bmatrix}\partial_{x_{k}}\\\
\partial_{y_{k}}\end{bmatrix}$ $\displaystyle=$
$\displaystyle\begin{bmatrix}\partial z_{k}/\partial x_{k}&\partial
z_{k}^{*}/\partial x_{k}\\\ \partial z_{k}/\partial y_{k}&\partial
z_{k}^{*}/\partial y_{k}\end{bmatrix}\begin{bmatrix}\partial_{z_{k}}\\\
\partial_{z_{k}^{*}}\end{bmatrix}$ (44) $\displaystyle=$
$\displaystyle\begin{bmatrix}1&1\\\
i&-i\end{bmatrix}\begin{bmatrix}\partial_{z_{k}}\\\
\partial_{z_{k}^{*}}\end{bmatrix}.$

By employing Eqs. (42), (43), and (44), we obtain the conclusion:

$\displaystyle\int d^{2}z_{k}e^{-|z_{k}|^{2}}z_{k}^{*}\hat{P}$ (45)
$\displaystyle=$ $\displaystyle\int
d^{2}z_{k}e^{-|z_{k}|^{2}}(x_{k}-iy_{k})\hat{P}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\int
d^{2}z_{k}e^{-|z_{k}|^{2}}\begin{bmatrix}1&-i\end{bmatrix}\begin{bmatrix}1&1\\\
i&-i\end{bmatrix}\begin{bmatrix}\partial_{z_{k}}\\\
\partial_{z_{k}^{*}}\end{bmatrix}\hat{P}$ $\displaystyle=$ $\displaystyle\int
d^{2}z_{k}e^{-|z_{k}|^{2}}\partial_{z_{k}}\hat{P}.$

Substituting Eq. (45) into Eq. (LABEL:zt), we obtain

$\displaystyle{\mathcal{M}}(z_{t}^{*}\hat{P})$ $\displaystyle=$
$\displaystyle-i\sum_{k}g_{k}e^{i\omega_{k}t}\int\frac{d^{2}z}{\pi}e^{-|z|^{2}}\frac{\partial\hat{P}}{\partial
z_{k}}$ (46) $\displaystyle=$ $\displaystyle\frac{\partial z_{t}^{*}}{\partial
z_{k}^{*}}\int\frac{d^{2}z}{\pi}e^{-|z|^{2}}\frac{\partial\hat{P}}{\partial
z_{k}}.$

Then, we apply the chain rule:

$\displaystyle\frac{\partial(\cdot)}{\partial
z_{k}}=\int_{0}^{t}ds\frac{\partial z_{s}}{\partial
z_{k}}\frac{\delta(\cdot)}{\delta z_{s}},$ (47)

and obtain the Novikov theorem[51, 46],

$\displaystyle{\mathcal{M}}(z_{t}^{*}\hat{P})$ $\displaystyle=$
$\displaystyle\frac{\partial z_{t}^{*}}{\partial
z_{k}^{*}}\int\frac{d^{2}z}{\pi}e^{-|z|^{2}}\int_{0}^{t}ds\frac{\partial
z_{s}}{\partial z_{k}}\frac{\delta\hat{P}}{\delta z_{s}}$ (48)
$\displaystyle=$
$\displaystyle\int\frac{d^{2}z}{\pi}e^{-|z|^{2}}\int_{0}^{t}ds\Big{[}\frac{\partial
z_{t}^{*}}{\partial z_{k}^{*}}\frac{\partial z_{s}}{\partial
z_{k}}\Big{]}\frac{\delta\hat{P}}{\delta z_{s}}$ $\displaystyle=$
$\displaystyle\int_{0}^{t}ds{\mathcal{M}}(z_{t}^{*}z_{s}){\mathcal{M}}(\frac{\delta\hat{P}}{\delta
z_{s}}).$

## References

  * Breuer and Petruccione [2002] H.-P. Breuer and F. Petruccione,  _The theory of open quantum systems_ (Oxford University Press, Oxford ; New York, 2002).
  * Gardiner and Zoller [2000] C. W. Gardiner and P. Zoller,  _Quantum Noise_ , 2nd ed., edited by H. Haken (Springer, 2000).
  * Diósi [1993a] L. Diósi, On high-temperature markovian equation for quantum brownian motion, Europhysics Letters 22, 1 (1993a).
  * Diósi [1993b] L. Diósi, Calderia-leggett master equation and medium temperatures, Physica A: Statistical Mechanics and its Applications 199, 517 (1993b).
  * Lindblad [1976] G. Lindblad, On the generators of quantum dynamical semigroups, Communications in Mathematical Physics 48, 119 (1976).
  * Redfield [1957] A. G. Redfield, On the theory of relaxation processes, IBM Journal of Research and Development 1, 19 (1957).
  * Jean  _et al._ [1992] J. M. Jean, R. A. Friesner, and G. R. Fleming, Application of a multilevel redfield theory to electron transfer in condensed phases, The Journal of Chemical Physics 96, 5827 (1992).
  * Egorova  _et al._ [2003] D. Egorova, M. Thoss, W. Domcke, and H. Wang, Modeling of ultrafast electron-transfer processes: Validity of multilevel redfield theory, The Journal of Chemical Physics 119, 2761 (2003).
  * Breuer [2007] H.-P. Breuer, Non-markovian generalization of the lindblad theory of open quantum systems, Phys. Rev. A 75, 022103 (2007).
  * Hartmann and Strunz [2020] R. Hartmann and W. T. Strunz, Accuracy assessment of perturbative master equations: Embracing nonpositivity, Phys. Rev. A 101, 012103 (2020).
  * Wang  _et al._ [2017] Y. H. Wang, F. R. Liu, W. Q. Li, T. Fan, J. F. Yang, Z. M. Wang, F. Liu, and N. X. Sun, Study of non-equilibrium thermal transport in $\text{G}e_{2}\text{S}b_{2}\text{T}e_{5}$ thin films under ultrafast laser excitation using a photo-excited carrier integrated semiconductor model, Journal of Applied Physics 122, 043104 (2017).
  * Walowski and Münzenberg [2016] J. Walowski and M. Münzenberg, Perspective: Ultrafast magnetism and thz spintronics, Journal of Applied Physics 120, 140901 (2016).
  * Zwaschka  _et al._ [2021] G. Zwaschka, F. Lapointe, R. K. Campen, and Y. Tong, Characterization of ultrafast processes at metal/solution interfaces: Towards femtoelectrochemistry, Current Opinion in Electrochemistry 29, 100813 (2021).
  * O’Donnell  _et al._ [2020] L. O’Donnell, K. Dholakia, and G. D. Bruce, High speed determination of laser wavelength using Poincaré descriptors of speckle, Optics Communications 459, 124906 (2020).
  * Metzger  _et al._ [2017] N. K. Metzger, R. Spesyvtsev, G. D. Bruce, B. Miller, G. T. Maker, G. Malcolm, M. Mazilu, and K. Dholakia, Harnessing speckle for a sub-femtometre resolved broadband wavemeter and laser stabilization, Nat Commun 8, 15610 (2017).
  * Link  _et al._ [2022] V. Link, K. Muller, R. G. Lena, K. Luoma, F. Damanet, W. T. Strunz, and A. J. Daley, Non-Markovian Quantum Dynamics in Strongly Coupled Multimode Cavities Conditioned on Continuous Measurement, PRX Quantum 3, 020348 (2022).
  * Tanimura and Wolynes [1991] Y. Tanimura and P. G. Wolynes, Quantum and classical fokker-planck equations for a gaussian-markovian noise bath, Phys. Rev. A 43, 4131 (1991).
  * Tanimura [2020] Y. Tanimura, Numerically “exact” approach to open quantum dynamics: The hierarchical equations of motion (HEOM), J. Chem. Phys. 153, 020901 (2020).
  * Gnesotto  _et al._ [2018] F. S. Gnesotto, F. Mura, J. Gladrow, and C. P. Broedersz, Broken detailed balance and non-equilibrium dynamics in living systems: a review, Reports on Progress in Physics 81, 066601 (2018).
  * Ferialdi [2016] L. Ferialdi, Exact closed master equation for gaussian non-markovian dynamics, Phys. Rev. Lett. 116, 120402 (2016).
  * Pomyalov and Tannor [2005] A. Pomyalov and D. J. Tannor, The non-markovian quantum master equation in the collective-mode representation: Application to barrier crossing in the intermediate friction regime, The Journal of Chemical Physics 123, 204111 (2005).
  * Chen  _et al._ [2014] Y. Chen, J. Q. You, and T. Yu, Exact non-Markovian master equations for multiple qubit systems: Quantum-trajectory approach, Phys. Rev. A 90, 052104 (2014).
  * Megier  _et al._ [2018] N. Megier, W. T. Strunz, C. Viviescas, and K. Luoma, Parametrization and Optimization of Gaussian Non-Markovian Unravelings for Open Quantum Dynamics, Phys. Rev. Lett. 120, 150402 (2018).
  * Chen  _et al._ [2018] Y. Chen, J. Q. You, and T. Yu, Non-markovian quantum interference in multilevel quantum systems: Exact master equation approach, Quantum Information and Computation 18, 11 (2018).
  * Lally  _et al._ [2022] S. Lally, N. Werren, J. Al-Khalili, and A. Rocco, Master equation for non-markovian quantum brownian motion: The emergence of lateral coherences, Phys. Rev. A 105, 012209 (2022).
  * Chen  _et al._ [2020] Y. Chen, Q. Ding, W. Shi, J. Jun, and T. Yu, Exact entanglement dynamics mediated by leaky optical cavities, J. Phys. B: At. Mol. Opt. Phys. 53, 125501 (2020).
  * Ban  _et al._ [2017] M. Ban, S. Kitajima, T. Arimitsu, and F. Shibata, Linear response theory for open systems: Quantum master equation approach, Phys. Rev. A 95, 022126 (2017).
  * Carballeira  _et al._ [2021] R. Carballeira, D. Dolgitzer, P. Zhao, D. Zeng, and Y. Chen, Stochastic Schrödinger equation derivation of non-Markovian two-time correlation functions, Sci Rep 11, 11828 (2021).
  * Ma  _et al._ [2020a] Y.-H. Ma, Q.-Z. Ding, and T. Yu, Persistent spin squeezing of a dissipative one-axis twisting model embedded in a general thermal environment, Phys. Rev. A 101, 022327 (2020a).
  * Ma  _et al._ [2020b] Y.-H. Ma, Q.-Z. Ding, and E. Wu, Entanglement of two nitrogen-vacancy ensembles via a nanotube, Phys. Rev. A 101, 022311 (2020b).
  * Ma  _et al._ [2021] Y.-H. Ma, Y. Xu, Q.-Z. Ding, and Y.-S. Chen, Preparation of spin squeezed state in siv centers coupled by diamond waveguide*, Chinese Physics B 30, 100311 (2021).
  * Qiao  _et al._ [2020] Y. Qiao, J. Zhang, Y. Chen, J. Jing, and S. Zhu, Quantumness protection for open systems in a double-layer environment, Sci. China Phys. Mech. Astron. 63, 250312 (2020).
  * Ding  _et al._ [2021] Q. Ding, P. Zhao, Y. Ma, and Y. Chen, Impact of the central frequency of environment on non-markovian dynamics in piezoelectric optomechanical devices, Scientific Reports 11, 1814 (2021).
  * Vacchini [2016] B. Vacchini, Generalized Master Equations Leading to Completely Positive Dynamics, Phys. Rev. Lett. 117, 230401 (2016).
  * Mølmer  _et al._ [1993] K. Mølmer, Y. Castin, and J. Dalibard, Monte carlo wave-function method in quantum optics, J. Opt. Soc. Am. B 10, 524 (1993).
  * Dalibard  _et al._ [1992] J. Dalibard, Y. Castin, and K. Mølmer, Wave-function approach to dissipative processes in quantum optics, Phys. Rev. Lett. 68, 580 (1992).
  * Carmichael [1993] H. Carmichael,  _An Open Systems Approach to Quantum Optics: Lectures Presented at the Université Libre de Bruxelles, October 28 to November 4, 1991_, v. 18 (Springer Berlin Heidelberg, 1993).
  * Dum  _et al._ [1992] R. Dum, P. Zoller, and H. Ritsch, Monte carlo simulation of the atomic master equation for spontaneous emission, Phys. Rev. A 45, 4879 (1992).
  * Link and Strunz [2017] V. Link and W. T. Strunz, Stochastic Feshbach Projection for the Dynamics of Open Quantum Systems, Phys. Rev. Lett. 119, 180401 (2017).
  * Barchielli and Pellegrini [2010] A. Barchielli and C. Pellegrini, Jump-diffusion unravelling of a non-markovian generalized lindblad master equation, Journal of Mathematical Physics 51, 112104 (2010), https://doi.org/10.1063/1.3514539 .
  * Piilo  _et al._ [2008] J. Piilo, S. Maniscalco, K. Härkönen, and K.-A. Suominen, Non-markovian quantum jumps, Phys. Rev. Lett. 100, 180402 (2008).
  * Gisin [1984] N. Gisin, Quantum measurements and stochastic processes, Phys. Rev. Lett. 52, 1657 (1984).
  * Gisin and Percival [1992] N. Gisin and I. C. Percival, The quantum-state diffusion model applied to open systems, Journal of Physics A: Mathematical and General 25, 5677 (1992).
  * Strunz [1996] W. T. Strunz, Linear quantum state diffusion for non-markovian open quantum systems, Physics Letters A 224, 25 (1996).
  * Strunz  _et al._ [1999] W. T. Strunz, L. Diósi, and N. Gisin, Open system dynamics with non-markovian quantum trajectories, Phys. Rev. Lett. 82, 1801 (1999).
  * Yu  _et al._ [1999] T. Yu, L. Diósi, N. Gisin, and W. T. Strunz, Non-markovian quantum-state diffusion: Perturbation approach, Phys. Rev. A 60, 91 (1999).
  * Fleming [2012] C. H. Fleming, Non-Markovian dynamics of open quantum systems: Stochastic equations and their perturbative solutions, Annals of Physics , 39 (2012).
  * Luo  _et al._ [2015] D.-W. Luo, P. V. Pyshkin, C.-H. Lam, T. Yu, H.-Q. Lin, J. Q. You, and L.-A. Wu, Dynamical invariants in non-Markovian quantum state diffusion equation, Phys. Rev. A 92, 062127 (2015), arXiv: 1510.00518.
  * Jing  _et al._ [2013] J. Jing, X. Zhao, J. Q. You, W. T. Strunz, and T. Yu, Many-body quantum trajectories of non-markovian open systems, Phys. Rev. A 88, 052122 (2013).
  * Jing and Yu [2010] J. Jing and T. Yu, Non-markovian relaxation of a three-level system: Quantum trajectory approach, Phys. Rev. Lett. 105, 240403 (2010).
  * Novikov [1965] E. A. Novikov, Functionals and the random-force method in turbulence theory, Soviet Journal of Experimental and Theoretical Physics 20, 1290 (1965).