# Separability transitions in topological states induced by local decoherence

Yu-Hsueh Chen  Department of Physics, University of California at San Diego,
La Jolla, California 92093, USA  Tarun Grover  Department of Physics,
University of California at San Diego, La Jolla, California 92093, USA

###### Abstract

We study states with intrinsic topological order subjected to local
decoherence from the perspective of separability, i.e., whether a decohered
mixed state can be expressed as an ensemble of short-range entangled pure
states. We focus on toric codes and the X-cube fracton state and provide
evidence for the existence of decoherence-induced separability transitions
that precisely coincide with the error-recovery transitions. A key insight is
that local decoherence acting on the ‘parent’ cluster states of these models
results in a Gibbs state.

In this work we will explore aspects of many-body topological states subjected
to decoherence from the perspective of separability, i.e., whether the
resulting mixed-state can be expressed as a convex sum of short-range
entangled (SRE) states [1, 2, 3]. This criteria is central to the definition
of what constitutes a trivial or a non-trivial mixed-state, and various mixed-
state entanglement measures such as negativity[4, 5, 6, 7, 8] and entanglement
of formation [9] are defined so as to quantify non-separability. At the same
time, currently there exist no faithful measures of mixed-state entanglement
that are efficiently calculable for many-body states [3] (a faithful measure
is zero if and only if a mixed-state can be written as a convex sum of
unentangled pure states). This limitation makes it rather challenging to know
if a mixed state is separable. We will be particularly interested in
decoherence-induced “separability transitions”, i.e., transitions tuned by
decoherence such that the density matrix in one regime is expressible as a
convex sum of SRE states, and in the other regime, it is not. One salient
distinction between pure state Vs mixed-state dynamics is that although a
short-depth unitary evolution cannot change long-range entanglement encoded in
a pure state, a short-depth local channel can fundamentally alter long-range
mixed-state entanglement. Therefore, even the limited class of mixed states
that are obtained by the action of local short-depth channels on an entangled
pure state offer an opportunity to explore mixed-state phases and phase
transitions [10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21]. We will focus on
mixed-states that are obtained via subjecting several well-understood
topologically ordered phases of matter to short-depth quantum channels.

Figure 1:  (a) Topological orders under local decoherence can undergo a
separability transition, where only above a certain critical error rate, the
decohered mixed state $\rho_{\textrm{dec}}$ can be written as a convex sum of
SRE pure states. The bottom depicts the parent cluster states and their
offspring models obtained by appropriate measurements (indicated by an arrow)
(b) 2d cluster Hamiltonian and 2d toric code, (c) 3d cluster Hamiltonian and
3d toric code, and (d) “Cluster-X” Hamiltonian [22] and the X-cube
Hamiltonian.

Error-threshold theorems [23, 24, 25, 26, 27, 28] indicate that if a pure
state has intrinsic topological order, then it will be perturbatively stable
against decoherence from a short-depth, local quantum channel. Therefore, the
mixed state is expected to undergo a phase transition as a function of the
decoherence rate [29]. Such transitions were originally studied from the
perspective of quantum error correction (QEC) in Refs.[30, 31] and more
recently using mixed-state entanglement measures such as topological
negativity [14], and other non-linear functions of the density matrix
(Refs.[14, 15, 13]). These approaches clearly establish at least two different
mixed-state phases: one where the topological qubit can be decoded, and the
other where it can’t. However, it is not obvious if the density matrix in the
regime where decoding fails can be expressed as a convex sum of SRE pure
states, which, following Refs.[1, 2], we will take as the definition of an SRE
mixed state. Our main result is that for several topologically ordered phases
subjected to local decoherence, one can explicitly write down the decohered
mixed state as a convex sum of pure states which we argue all undergo a
topological phase transition, from being long-ranged entangled to short-ranged
entangled, at a threshold that precisely corresponds to the optimal threshold
for QEC. Therefore, in these examples, we argue that the error-recovery
transition does indeed coincide with a separability transition. A key
observation in our analysis is that when the parent cluster states of these
models [32, 33, 22, 34] are subjected to local decoherence, one obtains a
Gibbs state. This naturally leads to the aforementioned ‘canonical’
decomposition of the decohered density matrix in terms of pure states which
captures the separability transition. This decomposition also provides a
direct connection between the separability transitions and statistical
mechanics models that naturally appear in the context of classical and quantum
error-correcting codes [35, 36, 37, 38, 30, 31, 39, 14].

Let us begin by considering the ground state of 2d toric code (see Fig.1(b))
with Hamiltonian $H_{\textrm{2d toric}}=-\sum_{v}(\prod_{e\in
v}Z_{e})-\sum_{p}(\prod_{e\in p}X_{e})$ subjected to phase-flip errors. Here
the Hilbert space consists of qubits residing on the edges (denoted as ‘$e$’)
of a square lattice and we assume periodic boundary conditions. If one denotes
the ground state as $\rho_{0}$, then the Kraus map corresponding to the phase-
flip errors act on an edge $e$ as:
$\mathcal{E}_{e}[\rho]=pZ_{e}\rho_{0}Z_{e}+(1-p)\rho_{0}$, and the full map is
given by the composition of this single-edge map over all edges. The key first
step is to identify a ‘parent’ cluster Hamiltonian (in the sense of Refs.[32,
33, 22, 34]) such that the application of the aforementioned Kraus map to its
ground state results in a Gibbs state. For example, for the problem at hand,
consider $H_{\textrm{2d cluster}}=\sum_{v}h_{v}+\sum_{e}h_{e}$ where
$h_{v}=-X_{v}(\prod_{e\ni v}Z_{e})$ and $h_{e}=-X_{e}(\prod_{v\in e}Z_{v})$,
see Fig.1(b). Note that here the Hilbert space consists of qubits both on the
vertices and the edges of the square lattice. A simple calculation shows that
the action of the aforementioned Kraus map on the ground state
$\rho_{C,0}\,(\propto\prod_{e}(I-h_{e})\prod_{v}(I-h_{v}))$ of $H_{\textrm{2d
Cluster}}$ yields $\mathcal{E}_{e}[\rho_{C,0}]\propto
e^{-\beta\sum_{e}h_{e}}\prod_{v}(I-h_{v})$ where $\tanh(\beta)=1-2p$. The
ground state density matrix $\rho_{0}$ of the 2d toric code can be written as
$\rho_{0}\propto\langle x_{\mathbf{v}}=1|\rho_{C,0}|x_{\mathbf{v}}=1\rangle$,
where $|x_{\mathbf{v}}=1\rangle=\otimes_{v}|x_{v}=1\rangle$ is the product
state in the Pauli-$X$ basis. The projection selects one specific ground state
of the toric code that is an eigenvector of the non-contractible Wilson loops
$W_{\ell}=\prod_{e\in\mathcal{\ell}}X_{e}$ with eigenvalue +1 along both
cycles $\mathcal{\ell}$ of the torus. Using the above Gibbs form of
$\mathcal{E}_{e}[\rho_{C,0}]$, the decohered density matrix of the toric code
is

$\displaystyle\rho$ $\displaystyle\propto\langle
x_{\mathbf{v}}=1|e^{-\beta\sum_{e}h_{e}}|x_{\mathbf{v}}=1\rangle P_{Z},$ (1)

where $P_{Z}=\prod_{v}(I+\prod_{e\ni v}Z_{e})$. By inserting a complete set of
states $\\{|x_{\mathbf{e}},z_{\mathbf{v}}\rangle\\}$ between
$e^{-\beta\sum_{e}h_{e}}$ and $|x_{\mathbf{v}}=1\rangle$, one may simplify the
above expression to obtain $\rho\propto P_{Z}\rho_{e}P_{Z}$ where
$\rho_{e}=\sum_{x_{\mathbf{e}}}\mathcal{Z}_{\text{2d
Ising},x_{\mathbf{e}}}|x_{\mathbf{e}}\rangle\langle x_{\mathbf{e}}|$ and
$\mathcal{Z}_{\text{2d
Ising},x_{\mathbf{e}}}=\sum_{z_{\mathbf{v}}}e^{\beta\sum_{e}x_{e}\prod_{v\in
e}z_{v}}$ is the partition function of the 2d Ising model with Ising
interactions determined by $\\{x_{e}\\}$. Thus,
$\rho\propto\sum_{x_{\mathbf{e}}}\mathcal{Z}_{\text{2d
Ising},x_{\mathbf{e}}}|\Omega_{x_{\mathbf{e}}}\rangle\langle\Omega_{x_{\mathbf{e}}}|$,
where $|\Omega_{x_{\mathbf{e}}}\rangle\propto\prod_{v}(I+\prod_{e\ni
v}Z_{e})|x_{\mathbf{e}}\rangle$ are nothing but a subset of toric code
eigenstates. Note that in this derivation, the 2d Ising model emerges due to
the $h_{e}$ terms in the parent cluster Hamiltonian, and ultimately, this will
lead to the relation between the separability transition and the statistical
mechanics of the 2d random-bond Ising model (RBIM) that also describes the
error-recovery transition [30]. We note that the above spectral representation
of $\rho$ in terms of toric code eigenstates has also previously appeared in
Ref.[13], using a different derivation. Since non-contractible cycles of the
torus will play an important role below, let us note that distinct eigenstates
$|\Omega_{x_{\mathbf{e}}}\rangle$ can be uniquely specified by two labels: the
first label corresponds to the set of local $\mathbb{Z}_{2}$ fluxes
$f_{p}=\prod_{e\in p}x_{e}$ through elementary plaquettes $p$, while the
second label $\mathbf{L}=(L_{x}=\pm 1,L_{y}=\pm 1)$ with
$L_{x}=\prod_{e\in\ell,e\parallel\hat{x}}x_{e},L_{y}=\prod_{e\in\ell,e\parallel\hat{y}}x_{e}$
and $\ell$ a non-contractible loop along $\hat{x}/\hat{y}$ direction,
specifies the topological sector (‘Logical data’) in which a given eigenstate
$|\Omega_{x_{\mathbf{e}}}\rangle$ lives.

We now probe the mixed state $\rho$ using the separability criteria, i.e., we
ask whether it can be decomposed as a convex sum of SRE states. Clearly, the
aforementioned spectral representation is not a useful decomposition since it
is written in terms of toric code eigenstates which are all long-range
entangled (LRE). Taking cue from the argument for separability of the Gibbs
state of toric codes [40], we decompose $\rho$ as
$\rho=\sum_{z_{\mathbf{e}}}\rho^{1/2}|z_{\mathbf{e}}\rangle\langle
z_{\mathbf{e}}|\rho^{1/2}\equiv\sum_{m}|\psi_{m}\rangle\langle\psi_{m}|$ where
$\\{z_{\mathbf{e}}\\}$ are a complete set of product states in the Pauli-$Z$
basis, and $|\psi_{m}\rangle=\rho^{1/2}|z_{\mathbf{e}}\rangle$. Generically,
to determine whether $\rho$ is an SRE mixed state, one needs to determine
whether each $|\psi_{m}\rangle$ is SRE. However, for the current case of
interest, it suffices to consider only $|\psi\rangle=\rho^{1/2}|m_{0}\rangle$
with $|m_{0}\rangle=|z_{\mathbf{e}}=1\rangle$. The reason is as follows. The
Gauss’s law ($\prod_{e\ni v}Z_{e}=1$) implies that the Hilbert space only
contains states that are closed loops in the $Z$ basis. Therefore, one may
write $|m\rangle=g_{x}|m_{0}\rangle$ where $g_{x}$ is a product of single-site
Pauli-$X$s forming closed loops. Since $[g_{x},\rho]=0$, this implies that
$|\psi_{m}\rangle\equiv|\psi_{g_{x}}\rangle=g_{x}|\psi\rangle$, and therefore,
if $|\psi\rangle$ is SRE (LRE), so is $|\psi_{g_{x}}\rangle$. $\rho(\beta)$
may then be written as:
$\rho(\beta)=\sum_{g_{x}}|\psi_{g_{x}}(\beta)\rangle\langle\psi_{g_{x}}(\beta)|$.
Now, using the aforementioned spectral representation of $\rho$, the (non-
normalized) state $|\psi\rangle=\rho^{1/2}|z_{\mathbf{e}}=1\rangle$ is:

$\displaystyle|\psi(\beta)\rangle$
$\displaystyle\propto\sum_{x_{\mathbf{e}}}[\mathcal{Z}_{\text{2d
Ising},x_{\mathbf{e}}}(\beta)]^{1/2}|x_{\mathbf{e}}\rangle,$ (2)

It is easy to see that when $\beta=\infty$,
$|\psi\rangle\propto|\Omega_{0}\rangle$, the non-decohered toric code ground
state, while when $\beta=0$, $|\psi\rangle\propto|z_{\mathbf{e}}=1\rangle$ is
a product state. This hints at a phase transition for $|\psi(\beta)\rangle$
from being an LRE state to an SRE state as we increase the error rate $p$
(i.e. decrease $\beta$). It is worth emphasizing that if we succeed in showing
that $\rho$ is a convex sum of SRE pure states, then long-range mixed-state
entanglement must be zero (as quantified by any valid mixed-state entanglement
measure, including ones which require an optimization over all possible
decompositions). Therefore, thinking in terms of convex decomposition allows
one to potentially leverage well-understood diagnostics of pure state
entanglement to make general statements about the complexity/entanglement of a
mixed state.

To locate the transition point $\beta_{c}$ for the state
$|\psi(\beta)\rangle$, we consider the expectation value of the ‘anyon
condensation operator’ (also known as ‘t Hooft loop) defined as [41, 42, 43,
15] $T_{\mathcal{\tilde{\ell}}}=\prod_{e\in\mathcal{\tilde{\ell}}}Z_{e}$,
where $\mathcal{\tilde{\ell}}$ denotes a homologically non-contractible loop
on the dual lattice (recall that in the language of $\mathbb{Z}_{2}$ gauge
theory [44, 43], $Z_{e}\sim e^{i\pi\textrm{(Electric field)}_{e}}$, and if one
starts with a ground state of toric code that is an eigenstate of Wilson loop
operators $W_{\ell}$ along non-contractible cycles $\ell$, e.g.
$|\Omega_{0}\rangle$, then the action of $T_{\mathcal{\tilde{\ell}}}$ takes
one to a different topological sector). Physically, $\langle
T_{\mathcal{\tilde{\ell}}}\rangle\equiv\langle\psi|T_{\mathcal{\tilde{\ell}}}|\psi\rangle/\langle\psi|\psi\rangle$
captures the amplitude of tunneling from one logical subspace to an orthogonal
one, and therefore it is zero in the $\mathbb{Z}_{2}$ topologically ordered
phase, and non-zero in the topologically trivial phase (=anyon condensed
phase). 111We note that the ‘t Hooft loop works as an order parameter for
detecting topological order in our pure state $|\psi\rangle$ because the
Gauss’s law ($\prod_{e\ni v}Z_{e}=1$) is satisfied exactly and therefore there
are no dynamical charges. Besides, note that $\tr(\rho
T_{\mathcal{\tilde{\ell}}})$, where $\rho$ is the decohered mixed state, is
identically zero, and clearly does not capture mixed-state entanglement. One
may easily verify that $\langle T_{\mathcal{\tilde{\ell}}}\rangle=0\,(1)$ when
$\beta=\infty\,(\beta=0)$. Using Eq.(2), the effect of
$T_{\mathcal{\tilde{\ell}}}$ on $|x_{\mathbf{e}}\rangle$ is to flip spins
along the curve $\mathcal{\tilde{\ell}}$ (i.e., $x_{e}\rightarrow-
x_{e},\forall e\in\mathcal{\tilde{\ell}}$), and we denote the corresponding
configuration as $x_{\mathcal{\tilde{\ell}},\mathbf{e}}$. Crucially, while
$x_{\mathcal{\tilde{\ell}},\mathbf{e}}$ and $x_{\mathbf{e}}$ have the same
flux through every elementary plaquette, they correspond to different values
of the logical label $\mathbf{L}$ and cannot be obtained from one to another
through a local transformation $x_{e}\rightarrow x_{e}\prod_{v\in
e}s_{v},s_{v}=\pm 1$. Therefore,
$T_{\mathcal{\tilde{\ell}}}|\psi\rangle\propto\sum_{x_{\mathbf{e}}}[\mathcal{Z}_{x_{\mathbf{e}}}]^{1/2}|x_{\mathcal{\tilde{\ell}},\mathbf{e}}\rangle=\sum_{x_{\mathbf{e}}}[\mathcal{Z}_{x_{\mathcal{\tilde{\ell}},\mathbf{e}}}]^{1/2}|x_{\mathbf{e}}\rangle$,
where we have changed the dummy variable $x_{\mathbf{e}}\rightarrow
x_{\mathcal{\tilde{\ell}},\mathbf{e}}$, and suppressed the subscript ‘2d
Ising’ under the partition function $\mathcal{Z}$ for notational convenience.
Thus,

$\displaystyle\langle T_{\mathcal{\tilde{\ell}}}\rangle$
$\displaystyle=\frac{\sum_{x_{\mathbf{e}}}\sqrt{\mathcal{Z}_{x_{\mathbf{e}}}\mathcal{Z}_{x_{\mathcal{\tilde{\ell}},\mathbf{e}}}}}{\sum_{x_{\mathbf{e}}}\mathcal{Z}_{x_{\mathbf{e}}}}=\frac{\sum_{x_{\mathbf{e}}}\mathcal{Z}_{x_{\mathbf{e}}}e^{-\Delta
F_{x_{\mathcal{\tilde{\ell}},\mathbf{e}}}/2}}{\sum_{x_{\mathbf{e}}}\mathcal{Z}_{x_{\mathbf{e}}}}$
(3) $\displaystyle=\langle e^{-\Delta F_{\mathcal{\tilde{\ell}}}/2}\rangle\geq
e^{-\langle\Delta F_{\mathcal{\tilde{\ell}}}\rangle/2}$

where $\Delta
F_{x_{\mathcal{\tilde{\ell}},\mathbf{e}}}=-\log(\mathcal{Z}_{x_{\mathcal{\tilde{\ell}},\mathbf{e}}}/\mathcal{Z}_{x_{\mathbf{e}}})$
is the free energy cost of inserting a domain wall of size $|\tilde{\ell}|\sim
L$ (= system’s linear size) in the RBIM along the Nishimori line [30], and we
have used Jensen’s inequality in the last sentence. Note that we are along the
Nishimori line because the probability of a given gauge invariant label
$\\{f_{\mathbf{p}},\mathbf{L}\\}$ along the Nishimori line is precisely the
partition function $\mathcal{Z}_{x_{\mathbf{e}}}$ [30]. Since $\langle\Delta
F_{\mathcal{\tilde{\ell}}}\rangle$, the disorder-averaged free energy cost,
diverges with $L$ in the ferromagnetic phase of the RBIM while converges to a
constant in the paramagnetic phase [30], Eq.(3) rigorously shows that for
$p>p_{c}=p_{\text{2d RBIM}}\approx 0.109$ [45], $\langle
T_{\mathcal{\tilde{\ell}}}\rangle$ saturates to a non-zero constant.
Therefore, $|\psi\rangle$ is a topologically trivial state when $p>p_{c}$, and
hence the mixed state is SRE for $p>p_{c}$. In contrast, for $p<p_{c}$, due to
non-vanishing ferromagnetic order (and associated domain wall cost) of the
RBIM, we expect that $\langle T_{\mathcal{\tilde{\ell}}}\rangle\sim
e^{-\langle\Delta F_{\mathcal{\tilde{\ell}}}\rangle/2}\sim e^{-cL}\rightarrow
0$ in the thermodynamic limit ($c>0$ is a constant), implying that
$|\psi\rangle$ is topologically ordered. This is also consistent with the non-
vanishing topological entanglement negativity for $p<p_{c}$ [46]. Therefore,
we conclude that for this model, error-recovery transition also corresponds to
a separability transition for the mixed-state.

Another diagnostic of topological order in pure states is the bipartite
(Renyi) topological entanglement entropy (TEE) [47, 48, 49]. Dividing the
whole system in real-space as $A\cup B$, we define the reduced density matrix
$\rho_{A}=\tr_{B}|\psi(\beta)\rangle\langle\psi(\beta)|$ for the state
$|\psi(\beta)\rangle$ of our interest (Eq.(2)). One finds (see Appendix A):

$\displaystyle\tr(\rho_{A}^{2})=$
$\displaystyle\frac{\sum_{\begin{subarray}{c}x_{\mathbf{e}},x^{\prime}_{\mathbf{e}}\end{subarray}}\mathcal{Z}_{x_{A},x_{B}}\mathcal{Z}_{x_{A}^{\prime},x_{B}^{\prime}}e^{-\Delta
F_{AB}(x_{\mathbf{e}},x^{\prime}_{\mathbf{e}})/2}}{\sum_{x_{\mathbf{e}},x^{\prime}_{\mathbf{e}}}\mathcal{Z}_{x_{A},x_{B}}\mathcal{Z}_{x^{\prime}_{A},x^{\prime}_{B}}}.$
(4)

Here $x_{A}(x_{B})$ denotes all the edges belonging to the region $A(B)$,
$\mathcal{Z}_{x_{A},x_{B}}$ denotes the partition function of the 2d Ising
model with the sign of Ising interactions determined by $x_{A}$ and $x_{B}$,
and $\Delta
F_{AB}(x_{\mathbf{e}},x^{\prime}_{\mathbf{e}})=-\log[\mathcal{Z}_{x_{A}^{\prime},x_{B}}\mathcal{Z}_{x_{A},x_{B}^{\prime}}/(\mathcal{Z}_{x_{A},x_{B}}\mathcal{Z}_{x_{A}^{\prime},x^{\prime}_{B}})]$
can be thought of as the free energy cost of swapping bonds between two copies
of RBIM in region $A$. Using Eq.(4), we provide a heuristic argument in
Appendix A that the TEE (= the system-size-independent subleading term for
$S_{2}=-\log\left(\tr(\rho^{2}_{A})\right)$ jumps from $\log(2)$ to zero at
$p_{c}=p_{\text{2d RBIM}}$. The main idea is that in the ferromagnetic phase
of the RBIM, the free energy penalty of creating a single Ising vortex leads
to a specific non-local constraint on the allowed configurations that
contribute to the sum in Eq.(4). The constraint is essentially that one needs
to minimize the free energy cost $\Delta
F_{AB}(x_{\mathbf{e}},x^{\prime}_{\mathbf{e}})$ for each fixed flux
configuration $\\{f_{p}\\}$ and $\\{f^{\prime}_{p}\\}$ corresponding to
$\\{x_{e}\\}$ and $\\{x^{\prime}_{e}\\}$ in Eq.(4). One finds that there
always exists a pair of configurations that contribute equally to
$\tr(\rho^{2}_{A})$ while satisfying the aforementioned constraint. This
results in a subleading contribution of $-\log(2)$ in the entanglement
entropy, that we identify as TEE. In the paramagnetic phase, the
aforementioned non-local constraint does not exist, and one therefore does not
expect a non-zero TEE. 222Note that the argument for TEE is reminiscent of the
argument in Refs.[50, 14] where non-zero TEE/topological negativity results
from a non-local constraint on the entanglement boundary. However, the precise
origin of the constraint is a bit different, as the calculation in Refs.[50,
14] is in a dual picture where the topological (trivial) phase corresponds to
a paramagnetic (ferromagnetic) phase. Therefore, we arrive at the same
conclusion as the one obtained using the anyon condensation operator.

Incidentally, Eq.(1) also allows us to construct an alternative convex
decomposition of the decohered mixed state $\rho$ that shows a phase
transition at a certain (non-optimal) threshold $p_{\textrm{non-optimal}}$
which is also related to the 2d RBIM. Despite being non-optimal, it is
conceptually interesting since the resulting separability threshold is related
via a Kramers-Wannier duality [51] to the RBIM. The main outcome is that
$\tanh(\beta_{\textrm{non-optimal}})=1-2p_{\textrm{non-optimal}}$ satisfies
$\tanh^{2}(\beta_{\textrm{non-optimal}}/2)=p_{\textrm{2d
RBIM}}/(1-p_{\textrm{2d RBIM}})$ which yields $p_{\textrm{non-optimal}}\approx
0.188$. See Appendix B for details where we also discuss an analogous non-
optimal decomposition for the 3d toric code.

Let us next consider 3d toric code with $H_{\textrm{3d
toric}}=-\sum_{f}(\prod_{e\in f}Z_{e})-\sum_{v}(\prod_{e\in v}X_{e})$ (see
Fig.1(c)). We again assume periodic boundary conditions. We will be interested
in subjecting the ground state of $H_{\textrm{3d toric}}$ to phase (bit)-flip
errors (non-trivial Kraus operators $\sim Z_{e}\,(X_{e})$). Let us first
consider phase-flip errors. We chose the parent cluster state as the ground
state $\rho_{C,0}$ of $H_{\textrm{3d cluster}}=-\sum_{e}X_{e}\prod_{f\ni
e}Z_{f}-\sum_{f}X_{f}\prod_{e\in f}Z_{e}=\sum_{e}h_{e}+\sum_{f}h_{f}$. The
corresponding ground state density matrix of the 3d toric code $\rho_{0}$ is
$\rho_{0}\propto\langle x_{\mathbf{f}}=1|\rho_{C,0}|x_{\mathbf{f}}=1\rangle$,
which is an eigenstate of the non-contractible ‘t Hooft membrane operators
$T_{xy},T_{yz},T_{zx}$ along the three planes with eigenvalue +1
($T_{xy}=\prod_{e\parallel z}X_{e}$ where the product is taken over all edges
parallel to the $z$ axis in any $xy$ plane. $T_{yz}$ and $T_{zx}$ are defined
analogously). Following essentially the same steps as in 2d toric code, one
obtains the decohered density matrix in the form
$\rho=\sum_{g_{x}}|\psi_{g_{x}}(\beta)\rangle\langle\psi_{g_{x}}(\beta)|$ with
$|\psi_{g_{x}}\rangle=g_{x}|\psi(\beta)\rangle$ and $g_{x}$ a product of
single-site Pauli-$X$s forming closed membranes. Therefore, we again only need
to analyze whether $|\psi\rangle\equiv\rho^{1/2}|z_{\mathbf{e}}=1\rangle$ is
SRE or LRE. Again, one may rewrite
$|\psi(\beta)\rangle\propto\sum_{x_{\mathbf{e}}}[\mathcal{Z}_{\text{3d
gauge},x_{\mathbf{e}}}(\beta)]^{1/2}|x_{\mathbf{e}}\rangle$ where
$\mathcal{Z}_{\text{3d
gauge},x_{\mathbf{e}}}=\sum_{z_{\mathbf{f}}}e^{\beta\sum_{e}x_{e}\prod_{f\ni
e}z_{f}}$ is now the partition function of a classical 3d Ising gauge theory
with the sign of each plaquette term determined by $\\{x_{e}\\}$. To probe the
topological transition as a function of $\beta$, we now consider the Wilson
loop operator $W_{\ell}=\prod_{e\in\ell}Z_{e}$, where $\mathcal{\ell}$ denotes
a homologically nontrivial cycle on the original lattice, say, along $z$ axis
(so that it pierces and anti-commutes with the aforementioned ‘t Hooft
operator $T_{xy}$). Essentially the same derivation as that for the 2d toric
code shows that $\langle W_{\ell}\rangle=\langle e^{-\Delta
F_{\ell}/2}\rangle\geq e^{-\langle\Delta F_{\ell}\rangle/2}$, where $\Delta
F_{\ell}$ now denotes the free energy cost of inserting a domain wall along
the non-contractible loop for the 3d random-plaquette gauge model (RPGM) along
the Nishimori line. Since $\langle\Delta F_{\ell}\rangle$ diverges as the
length of $|\mathcal{C}|\sim L$ (= system-size) in the Higgs (ordered) phase,
while converges to a constant in the confinement (disordered) phase [31], one
finds $\langle W_{\ell}\rangle$ saturates to a non-zero constant when
$p>p_{\text{3d RPGM}}\approx 0.029$ [31], while it vanishes for $p<p_{c}$.
Therefore, we expect that $|\psi\rangle$, and correspondingly the decohered
state $\rho$, is SRE when $p>p_{\text{3d RPGM}}$.

Let us next consider bit-flip errors. The aforementioned $H_{\textrm{3d
cluster}}$ is not a good starting point for our approach since we would like
the corresponding ground state to take a Gibbs form when subjected to Kraus
operators $\sim X_{e}$. Therefore, one instead considers
$H^{\prime}_{\textrm{3d Cluster}}=-\sum_{v}Z_{v}(\prod_{e\ni
v}X_{e})-\sum_{e}Z_{e}(\prod_{v\in e}X_{v})$, which has previously also
appeared in Ref.[52]. The rest of the analysis is quite similar to that for
the 2d toric code (with $X\leftrightarrow Z$ everywhere). After writing the
non-decohered density matrix of toric code as $\rho_{0}\propto\langle
z_{\mathbf{v}}=1|\rho_{C,0}|z_{\mathbf{v}}=1\rangle$, where $\rho_{C,0}$ is
the ground state of $H^{\prime}_{\textrm{3d Cluster}}$, and subjecting it to
bit-flip errors, the decohered state is schematically given by
$\rho\propto\sum_{z_{\mathbf{e}}}\mathcal{Z}_{\text{3d
Ising},z_{\mathbf{e}}}|\Omega_{z_{\mathbf{e}}}\rangle\langle\Omega_{z_{\mathbf{e}}}|$,
where $|\Omega_{z_{\mathbf{e}}}\rangle$ are toric code eigenstates and
$\mathcal{Z}_{\text{3d Ising},z_{\mathbf{e}}}$ is the partition function of
the 3d classical Ising model with interactions determined by $\\{z_{e}\\}$.
The analog of the state $|\psi(\beta)\rangle$ is
$|\psi(\beta)\rangle\propto\sum_{x_{\mathbf{e}}}[\mathcal{Z}_{\text{3d
Ising},z_{\mathbf{e}}}(\beta)]^{1/2}|z_{\mathbf{e}}\rangle$ and its
topological transition is indicated by the non-analyticity of the ‘t Hooft
operator similar to Eq.(3). The corresponding $p_{c}$ for the separability
transition is then determined by the transition out of the ferromagnetic phase
in the 3d RBIM along the Nishimori line, which matches the optimal error-
recovery threshold, $p_{c}\approx 0.233$ [53, 54].

Finally, let us briefly consider the 3d X-cube model (Ref.[55]), where the
Hilbert space consists of qubits residing on the edges(e) of a cubic lattice,
and the Hamiltonian is $H_{\textrm{X-cube}}=-\sum_{c}\prod_{e\in
c}Z_{e}-\sum_{v}(\prod_{e\in v_{x}}X_{e}+\prod_{e\in v_{y}}X_{e}+\prod_{e\in
v_{z}}X_{e})=-\sum_{c}A_{c}-\sum_{v}(B_{v_{x}}+B_{v_{y}}+B_{v_{z}})$ where
$e\in v_{\gamma},\ \gamma=x,y,z$ denotes all the edges emanating from the
vertex $v$ that are normal to the $\gamma$-direction (see Fig.1(d)). It was
shown in Ref.[22] that the ground state density matrix $\rho_{0}$ of the
X-cube model can be written as $\rho_{0}\propto\langle
x_{\mathbf{c}}=1|\rho_{C,0}|x_{\mathbf{c}}=1\rangle$ where
$\rho_{C,0}=\prod_{c}(I-h_{c})\prod_{e}(I-h_{e})$ ($h_{c}=-X_{c}\prod_{e\in
c}Z_{e}$ and $h_{e}=-X_{e}\prod_{c\ni e}Z_{c}$ denotes the ground state
density matrix of the parent cluster state, and
$|x_{\mathbf{c}}=1\rangle=\otimes_{c}|x_{c}=1\rangle$ is the product state in
the Pauli-$X$ basis. Note that the qubits in the parent cluster state live at
the edges and the centers of the cubes so that $h_{c}$ involves 13-qubit
interactions, and $h_{e}$ involves 5-qubit interactions (Fig.1(d)). Similar to
the previous cases, the density matrix after subjecting $\rho_{0}$ to the
phase-flip channel (Kraus operators $\sim Z_{e}$) can be written as
$\rho\propto\sum_{x_{\mathbf{e}}}\mathcal{Z}_{\text{3d
plaquette},x_{\mathbf{e}}}|\Omega_{x_{\mathbf{e}}}\rangle\langle\Omega_{x_{\mathbf{e}}}|$,
where $|\Omega_{x_{\mathbf{e}}}\rangle\propto\prod_{e}(I+\prod_{e\in
c}Z_{e})|x_{\mathbf{e}}\rangle$ and $\mathcal{Z}_{\text{3d
plaquette},x_{\mathbf{e}}}=\sum_{z_{\mathbf{c}}}e^{\beta\sum_{e}x_{e}\prod_{c\ni
e}z_{c}}$ is the partition function of the 3d plaquette Ising model [56] with
the sign of interaction on each plaquette determined by $\\{x_{e}\\}$. The
appearance of this model is not a coincidence since the statistical mechanics
of the error-recovery transition involves precisely this model [57]. One again
only needs to analyze the state
$|\psi\rangle=\rho^{1/2}|z_{\mathbf{e}}=1\rangle$ to study the separability
transition for $\rho$. Now there exists exponentially many topological sectors
[55], and in the non-decohered ground state $\rho_{0}$, the membrane operators
defined as $\prod_{e\parallel\hat{a}}X_{e}$ with $a=x,y,z$ for any plane have
expectation value one. To detect the presence/absence of topological order in
$|\psi\rangle$, one therefore considers non-contractible Wilson loop operators
$W_{\ell}=\prod_{e\in\ell}Z_{e}$ that anti-commute with the membrane operators
orthogonal to $\ell$. The expectation value of any such Wilson loop takes a
form similar to Eq.(3) where the partition function
$\mathcal{Z}_{x_{\mathbf{e}}}=\mathcal{Z}_{\text{3d
plaquette},x_{\mathbf{e}}}$ and one is again along the Nishimori line. This
again indicates that the pure state $|\psi\rangle$ undergoes a transition at
the error threshold $p_{c}=p_{\text{3d plaquette}}\approx 0.152$ [57].

We note that for all the examples considered, one may define a more general
class of wavefunctions
$|\psi^{(\alpha)}\rangle\propto\rho^{\alpha/2}|z_{\mathbf{e}}=1\rangle$.
Interestingly, the wavefunction corresponding to $\alpha=2$ has been studied
in detail previously in Ref.[50], with $\rho$ identical to our decohered
density matrix for the 2d toric code. One finds that $|\psi^{(2)}\rangle$
undergoes a phase transition from being topologically ordered to an SRE state
at a temperature that precisely corresponds to the critical temperature of the
2d translationally invariant classical Ising model (with all couplings set to
unity) [50]. Alternatively, one may locate the critical point from the non-
analyticity of the wavefunction overlap of the non-normalized wavefunction
$|\psi(\beta)\rangle$. Specifically, one finds that
$\log(\langle\psi^{(2)}(\beta)|\psi^{(2)}(\beta)\rangle)$ is proportional to
the free energy of the classical 2d Ising model, and therefore, is non-
analytic precisely at the same point where the TEE jumps from $\log(2)$ to
zero. One may think of the wavefunction overlap as a kind of ‘partition
function’ since it enters in the denominator for expectation value of various
observables. This motivates a generalization of this overlap to general
$\alpha$ for any of the models considered above by defining
$F_{\alpha}(\beta)=\frac{1}{1-\alpha}\log(\langle\psi^{(\alpha)}(\beta)|\psi^{(\alpha)}(\beta)\rangle)$.
Taking the limit $\alpha\to 1$, which corresponds to the wavefunction of our
main interest (Eq.(2) for the 2d toric code, and analogous states for the
other two models), one finds that $F_{1}(\beta)$ precisely corresponds to the
free energy of the corresponding statistical mechanics model along the
Nishimori line, which indeed shows a singularity at the optimal error-recovery
threshold $p_{c}$. This provides an alternative, albeit heuristic, approach to
locate the phase transition.

To summarize, we provided evidence that decoherence-induced separability
transitions in several topological states coincide with the optimal threshold
for QEC [30, 14, 15, 13, 57]. This implies that in these models, the inability
to correct logical errors implies an ability to prepare the mixed state using
an ensemble of short-depth unitary circuits. The convex decomposition we
constructed captures the universal aspects of the phase diagram, as well as
the threshold correctly, and it is ‘canonical’ or optimal in this sense
(although it need not be unique). It will be interesting to consider applying
our method to other states that can be obtained from cluster states by
measurements [22, 34], including non-CSS codes, e.g., the double semion
topological order. Relatedly, using cluster states to obtain decohered density
matrices is seemingly an alternative method to generate statistical mechanics
models for error-correcting codes and it will be interesting to probe its
generality. It will also be worthwhile to numerically study the entanglement
structure of the state in Eq.(2) (and its analogs for the other models) using
Quantum Monte Carlo. A more ambitious task would be to construct an explicit
short-depth circuit for this state when $p>p_{c}$.

###### Acknowledgements.

The authors thank Dan Arovas, Tim Hsieh and John McGreevy for helpful
discussions. TG is supported by the National Science Foundation under Grant
No. DMR-1752417. This research was supported in part by the National Science
Foundation under Grant No. NSF PHY-1748958.

## References

  * Werner [1989] R. F. Werner, Phys. Rev. A 40, 4277 (1989).
  * Hastings [2011] M. B. Hastings, Physical review letters 107, 210501 (2011).
  * Horodecki  _et al._ [2009] R. Horodecki, P. Horodecki, M. Horodecki, and K. Horodecki, Reviews of modern physics 81, 865 (2009).
  * Peres [1996] A. Peres, Phys. Rev. Lett. 77, 1413 (1996).
  * Horodecki  _et al._ [1996] M. Horodecki, P. Horodecki, and R. Horodecki, Physics Letters A 223, 1 (1996).
  * Eisert and Plenio [1999] J. Eisert and M. B. Plenio, Journal of Modern Optics 46, 145 (1999).
  * Vidal and Werner [2002] G. Vidal and R. F. Werner, Phys. Rev. A 65, 032314 (2002).
  * Plenio [2005] M. B. Plenio, Physical review letters 95, 090503 (2005).
  * Bennett  _et al._ [1996] C. H. Bennett, D. P. DiVincenzo, J. A. Smolin, and W. K. Wootters, Phys. Rev. A 54, 3824 (1996).
  * de Groot  _et al._ [2022] C. de Groot, A. Turzillo, and N. Schuch, Quantum 6, 856 (2022).
  * Ma and Wang [2022] R. Ma and C. Wang, arXiv preprint arXiv:2209.02723 (2022).
  * Lee  _et al._ [2022a] J. Y. Lee, Y.-Z. You, and C. Xu, arXiv preprint arXiv:2210.16323 (2022a).
  * Lee  _et al._ [2023] J. Y. Lee, C.-M. Jian, and C. Xu, arXiv preprint arXiv:2301.05238 (2023).
  * Fan  _et al._ [2023] R. Fan, Y. Bao, E. Altman, and A. Vishwanath, arXiv preprint arXiv:2301.05689 (2023).
  * Bao  _et al._ [2023] Y. Bao, R. Fan, A. Vishwanath, and E. Altman, arXiv preprint arXiv:2301.05687 (2023).
  * Zou  _et al._ [2023] Y. Zou, S. Sang, and T. H. Hsieh, Physical Review Letters 130, 250403 (2023).
  * Lu  _et al._ [2023] T.-C. Lu, Z. Zhang, S. Vijay, and T. H. Hsieh, PRX Quantum 4, 030318 (2023).
  * Ma  _et al._ [2023] R. Ma, J.-H. Zhang, Z. Bi, M. Cheng, and C. Wang, arXiv preprint arXiv:2305.16399 (2023).
  * Ma [2023] R. Ma, arXiv preprint arXiv:2304.08277 (2023).
  * Su  _et al._ [2023] K. Su, N. Myerson-Jain, and C. Xu, arXiv preprint arXiv:2305.13410 (2023).
  * Wang  _et al._ [2023] Z. Wang, Z. Wu, and Z. Wang, arXiv preprint arXiv:2307.13758 (2023).
  * Verresen  _et al._ [2021] R. Verresen, N. Tantivasadakarn, and A. Vishwanath, arXiv preprint arXiv:2112.03061 (2021).
  * Shor [1996] P. W. Shor, Proceedings of 37th conference on foundations of computer science , 56 (1996).
  * Aharonov and Ben-Or [1997] D. Aharonov and M. Ben-Or, Proceedings of the twenty-ninth annual ACM symposium on Theory of computing , 176 (1997).
  * Kitaev [2003] A. Y. Kitaev, Annals of Physics 303, 2 (2003).
  * Knill  _et al._ [1998] E. Knill, R. Laflamme, and W. H. Zurek, Science 279, 342 (1998).
  * Preskill [1998] J. Preskill, Proceedings of the Royal Society of London. Series A: Mathematical, Physical and Engineering Sciences 454, 385 (1998).
  * Terhal [2015] B. M. Terhal, Reviews of Modern Physics 87, 307 (2015).
  * Aharonov [2000] D. Aharonov, Physical Review A 62, 062311 (2000).
  * Dennis  _et al._ [2002] E. Dennis, A. Kitaev, A. Landahl, and J. Preskill, Journal of Mathematical Physics 43, 4452 (2002).
  * Wang  _et al._ [2003] C. Wang, J. Harrington, and J. Preskill, Annals of Physics 303, 31 (2003).
  * Raussendorf  _et al._ [2005] R. Raussendorf, S. Bravyi, and J. Harrington, Physical Review A 71, 062313 (2005).
  * Aguado  _et al._ [2008] M. Aguado, G. Brennen, F. Verstraete, and J. I. Cirac, Physical review letters 101, 260501 (2008).
  * Tantivasadakarn  _et al._ [2021] N. Tantivasadakarn, R. Thorngren, A. Vishwanath, and R. Verresen, arXiv preprint arXiv:2112.01519 (2021).
  * Nishimori [1981] H. Nishimori, Progress of Theoretical Physics 66, 1169 (1981).
  * Sourlas [1989] N. Sourlas, Nature 339, 693 (1989).
  * Ruján [1993] P. Ruján, Physical review letters 70, 2968 (1993).
  * Nishimori [1993] H. Nishimori, Journal of the Physical Society of Japan 62, 2973 (1993).
  * Mezard and Montanari [2009] M. Mezard and A. Montanari,  _Information, physics, and computation_ (Oxford University Press, 2009).
  * Lu  _et al._ [2020] T.-C. Lu, T. H. Hsieh, and T. Grover, Physical Review Letters 125, 116801 (2020).
  * Hooft [1978] G. Hooft, Nuclear Physics: B 138, 1 (1978).
  * Kogut [1979] J. B. Kogut, Reviews of Modern Physics 51, 659 (1979).
  * Fradkin [2013] E. Fradkin,  _Field theories of condensed matter physics_ (Cambridge University Press, 2013).
  * Kogut and Susskind [1975] J. Kogut and L. Susskind, Physical Review D 11, 395 (1975).
  * Honecker  _et al._ [2001] A. Honecker, M. Picco, and P. Pujol, Physical review letters 87, 047201 (2001).
  * Fan [2022] R. Fan, Physical Review Letters 129, 260403 (2022).
  * Levin and Wen [2006] M. Levin and X.-G. Wen, Phys. Rev. Lett. 96, 110405 (2006).
  * Kitaev and Preskill [2006] A. Kitaev and J. Preskill, Phys. Rev. Lett. 96, 110404 (2006).
  * Flammia  _et al._ [2009] S. T. Flammia, A. Hamma, T. L. Hughes, and X.-G. Wen, Physical review letters 103, 261601 (2009).
  * Castelnovo and Chamon [2008] C. Castelnovo and C. Chamon, Physical Review B 77, 054433 (2008).
  * Kramers and Wannier [1941] H. A. Kramers and G. H. Wannier, Physical Review 60, 252 (1941).
  * Lee  _et al._ [2022b] J. Y. Lee, W. Ji, Z. Bi, and M. Fisher, arXiv preprint arXiv:2208.11699 (2022b).
  * Ozeki and Ito [1998] Y. Ozeki and N. Ito, Journal of Physics A: Mathematical and General 31, 5451 (1998).
  * Kubica  _et al._ [2018] A. Kubica, M. E. Beverland, F. Brandão, J. Preskill, and K. M. Svore, Physical review letters 120, 180501 (2018).
  * Vijay  _et al._ [2016] S. Vijay, J. Haah, and L. Fu, Physical Review B 94, 235157 (2016).
  * Savvidy and Wegner [1994] G. Savvidy and F. Wegner, Nuclear Physics B 413, 605 (1994).
  * Song  _et al._ [2022] H. Song, J. Schönmeier-Kromer, K. Liu, O. Viyuela, L. Pollet, M. A. Martin-Delgado, _et al._ , Physical Review Letters 129, 230502 (2022).
  * Dong  _et al._ [2008] S. Dong, E. Fradkin, R. G. Leigh, and S. Nowling, Journal of High Energy Physics 2008, 016 (2008).
  * Zhang  _et al._ [2012] Y. Zhang, T. Grover, A. Turner, M. Oshikawa, and A. Vishwanath, Physical Review B 85, 235151 (2012).

##  Appendix A Entanglement entropy of the proposed optimal pure state for 2d
toric code

In this section, we compute the 2nd Renyi entanglement entropy of the
(normalized) state in Eq.(2), i.e.,
$|\tilde{\psi}\rangle={\sum_{x_{\mathbf{e}}}\mathcal{\sqrt{Z}}_{x_{\mathbf{e}}}|x_{\mathbf{e}}\rangle}/{\sqrt{\sum_{x_{\mathbf{e}}}\mathcal{Z}_{x_{\mathbf{e}}}}}$
with
$\mathcal{Z}_{x_{\mathbf{e}}}=\sum_{z_{\mathbf{v}}}e^{\beta\sum_{e}x_{e}\prod_{v\in
e}z_{v}}$, the partition funcion of the 2d Ising model with the Ising
interactions determined by $\\{x_{\mathbf{e}}\\}$. We will later provide a
heuristic argument that as one increases the error rate $p$ (which corresponds
to increasing the temperature $\beta^{-1}$ in the Ising model via the relation
$\tanh(\beta)=1-2p$), the topological Renyi entropy undergoes a transition
from $\log 2$ to $0$ at $p_{c}=p_{\text{2d RBIM}}$. We will only consider an
entanglement bipartition such that the entangling boundary is contractible.
Since the topological Renyi entropy for a contractible region does not depend
on the topological sector $\mathbf{L}$ in which the wavefunction resides [47,
48, 58, 59], below we will ignore the dependence of the partition function on
$\mathbf{L}$ (alternatively, one may assume open boundary conditions so that
there is a unique ground state).

The calculation of the 2nd Renyi entanglement entropy for a bipartition of the
total system $A\cup B$ can be divided into two steps. First, we compute the
reduced density matrix
$\rho_{A}=\tr_{B}(|\tilde{\psi}\rangle\langle\tilde{\psi}|)$. Next, we compute
$\tr\left(\rho^{2}_{A}\right)$ to get the 2nd Renyi entropy
$S_{2}=-\log[\tr(\rho_{A}^{2})]$. Denoting $x_{A}(x_{B})$ as the set of the
links belonging to region $A(B)$, the pure state density matrix can be written
as

$|\tilde{\psi}\rangle\langle\tilde{\psi}|=\frac{\sum_{x_{\mathbf{e}},x^{\prime}_{\mathbf{e}}}\sqrt{\mathcal{Z}_{x_{A},x_{B}}}\sqrt{\mathcal{Z}_{x_{A}^{\prime},x_{B}^{\prime}}}|x_{A},x_{B}\rangle\langle
x_{A}^{\prime},x_{B}^{\prime}|}{(\sum_{x_{e}}\mathcal{Z}_{x_{\mathbf{e}}})}$
(5)

Here $x_{A}(x_{B})$ denotes all the edges belonging to the region $A(B)$ and
$\mathcal{Z}_{x_{A},x_{B}}$ denotes the partition of the 2d Ising model with
the sign of Ising interactions determined by $x_{A}$ and $x_{B}$. Note that
any given bond either belongs only to $A$ or only to $B$, see Fig.2. One may
then compute $\rho_{A}=\tr_{B}(|\tilde{\psi}\rangle\langle\tilde{\psi}|)$ to
obtain

$\displaystyle\rho_{A}$
$\displaystyle=\frac{\sum_{x_{\mathbf{e}},x^{\prime}_{\mathbf{e}}}\sqrt{\mathcal{Z}_{x_{A},x_{B}}}\sqrt{\mathcal{Z}_{x_{A}^{\prime},x_{B}}}|x_{A}\rangle\langle
x_{A}^{\prime}|}{(\sum_{x_{e}}\mathcal{Z}_{x_{\mathbf{e}}})}.$ (6)

Therefore,

$\displaystyle\begin{aligned}
\tr(\rho_{A}^{2})&=\frac{\sum_{\begin{subarray}{c}x_{\mathbf{e}},x^{\prime}_{\mathbf{e}}\end{subarray}}\sqrt{\mathcal{Z}_{x_{A},x_{B}}\mathcal{Z}_{x_{A}^{\prime},x_{B}}\mathcal{Z}_{x_{A},x_{B}^{\prime}}\mathcal{Z}_{x_{A}^{\prime},x_{B}^{\prime}}}}{\sum_{x_{\mathbf{e}},x^{\prime}_{\mathbf{e}}}\mathcal{Z}_{x_{A},x_{B}}\mathcal{Z}_{x^{\prime}_{A},x^{\prime}_{B}}}\\\
&=\frac{\sum_{\begin{subarray}{c}x_{\mathbf{e}},x^{\prime}_{\mathbf{e}}\end{subarray}}\mathcal{Z}_{x_{A},x_{B}}\mathcal{Z}_{x_{A}^{\prime},x_{B}^{\prime}}e^{-\Delta
F_{AB}(x_{\mathbf{e}},x^{\prime}_{\mathbf{e}})/2}}{\sum_{x_{\mathbf{e}},x^{\prime}_{\mathbf{e}}}\mathcal{Z}_{x_{A},x_{B}}\mathcal{Z}_{x^{\prime}_{A},x^{\prime}_{B}}}\\\
\end{aligned}$ (7)

where

$\Delta
F_{AB}(x_{\mathbf{e}},x^{\prime}_{\mathbf{e}})=-\log(\frac{\mathcal{Z}_{x_{A}^{\prime},x_{B}}\mathcal{Z}_{x_{A},x_{B}^{\prime}}}{\mathcal{Z}_{x_{A}^{\prime},x^{\prime}_{B}}\mathcal{Z}_{x_{A},x_{B}}})$
(8)

is the free energy cost of swapping the bonds of two copies of RBIM (whose
respective probabilities are given by $\mathcal{Z}_{x_{A},x_{B}}$ and
$\mathcal{Z}_{x_{A}^{\prime},x^{\prime}_{B}}$) in region $A$. Therefore, one
finds that the second Renyi entropy $S_{2}=-\log[\tr(\rho_{A}^{2})]$ is
related to the averaged free energy cost along the Nishimori line of swapping
the bonds of two copies of RBIM in region $A$. The 2d RBIM along the Nishimori
line exhibits two phases: a ferromagnetic phase when $p<p_{\text{2d RBIM}}$
and a paramagnetic phase when $p>p_{\text{2d RBIM}}$. We now argue that the
topological Renyi entropy is $\log 2$ ($0$) in the ferromagnetic
(paramagnetic) phase.

To begin with, let’s first consider the extreme limit $\beta=\infty$. Since
$|\psi(\beta=\infty)\rangle$ is just the toric code ground state, we already
know that the TEE is $\log 2$. However, it is illuminating to see how this
subleading correction emerges from Eq.(7) in this limit. We first note that
the probability for a given configuration
$(x_{\mathbf{e}},x_{\mathbf{e}}^{\prime})$ to occur is completely determined
by $\mathcal{Z}_{x_{A},x_{B}}\mathcal{Z}_{x^{\prime}_{A},x^{\prime}_{B}}$,
which is merely a function of fluxes $f_{\mathbf{p}}(x_{\mathbf{e}})$ and
$f_{\mathbf{p}}(x_{\mathbf{e}}^{\prime})$, where we denote
$f_{\mathbf{p}}(x_{\mathbf{e}})=\\{f_{p}(x_{\mathbf{e}})\\}$ as the flux
configuration determined by $x_{\mathbf{e}}$ (recall that $f_{p}=\prod_{e\in
p}x_{e}$). As noted in the first paragraph of this appendix, if the system is
on a torus, then although $\mathcal{Z}$ technically also depends on the
logical sector $\mathbf{L}$, such dependence is not important for our
discussion here, as we are working with a region $A$ whose boundary is
contractible. When $\beta=\infty$, the partition function belonging to the
vacuum sector (i.e., $f_{p}=1,\forall p$) dominates. Therefore, the dominant
terms in the denominator of Eq.(7) are the configurations
$(x_{\mathbf{e}},x_{\mathbf{e}}^{\prime})$ that ensure
$f_{p}(x_{\mathbf{e}})=f_{p}(x_{\mathbf{e}}^{\prime})=1,\ \forall p$. On the
other hand, the dominant terms in the numerator of Eq.(7) should satisfy the
extra condition that $\Delta
F_{AB}(x_{\mathbf{e}},x_{\mathbf{e}}^{\prime})=0$. This forces the swapped
configurations $(x_{A}^{\prime},x_{B})$ and $(x_{A},x_{B}^{\prime})$ to also
lie in the vacuum sector, i.e.,
$f_{p}(x^{\prime}_{A},x_{B})=f_{p}(x_{A},x_{B}^{\prime})=1,\forall p$, where
we denote $f_{p}(x^{\prime}_{A},x_{B})$ as the flux determined by the bonds
$(x^{\prime}_{A},x_{B})$. Therefore, Eq.(7) in the limit $\beta=\infty$ simply
counts the degrees of freedom lost due to the additional constraint $\Delta
F_{AB}(x_{\mathbf{e}},x_{\mathbf{e}}^{\prime})=0$ in the vacuum sectors of
$x_{\mathbf{e}}$ and $x_{\mathbf{e}}^{\prime}$.

Figure 2: Given the boundary transformation $l_{\mathbf{v}}$ on
$x_{\mathbf{v}}$ in (a), one can choose the boundary transformation
$l^{\prime}_{\mathbf{v}}$ on $x^{\prime}_{\mathbf{e}}$ in (b) to be either
$l^{\prime}_{\mathbf{v}}=l_{\mathbf{v}}$(not shown) or
$l^{\prime}_{\mathbf{v}}=\partial A+l_{\mathbf{v}}\ (\text{mod}\ 2)$ to make
the boundary fluxes $f_{\mathbf{p}}{(x^{\prime}_{A},x_{B})}$ and
$f_{\mathbf{p}}{(x_{A},x^{\prime}_{B})}$ remain invariant. Recall that the
transformation is defined on vertices through $x_{e}\rightarrow
x_{e}\prod_{v\in e}\sigma_{v},\sigma_{v}=\pm 1$. Here we use blue circles to
denote the vertices $\sigma_{v}=-1$ and blue squares to denote the edges where
$\prod_{v\in e}\sigma_{v}=-1$.

Based on above considerations, to evaluate Eq.(7) for the case $\beta=\infty$,
we recall that different configurations $x_{\mathbf{e}}$ in the vacuum sector
can be obtained by applying the transformation (defined on the vertices)
$x_{e}\rightarrow x_{e}\prod_{v\in e}\sigma_{v},\sigma_{v}=\pm 1$ to the
disorder-free sector, i.e., $x_{e}=1,\forall e$. Therefore, without the
constraint $\Delta F_{AB}(x_{\mathbf{e}},x_{\mathbf{e}}^{\prime})=0$, there
are $2^{2|\partial A|}$ different transformations on the boundary vertices
(i.e., the vertices living on the boundary between region $A$ and $B$) that
make $x_{\mathbf{e}}$ and $x_{\mathbf{e}}^{\prime}$ stay in the vacuum
sectors, where $|\partial A|$ denotes the number of boundary vertices (for
example, $|\partial A|=10$ in Fig.2). On the other hand, after enforcing the
constraint $\Delta F_{AB}(x_{\mathbf{e}},x_{\mathbf{e}}^{\prime})=0$, the
transformation on the boundary vertices for $x_{\mathbf{e}}$ and
$x_{\mathbf{e}}^{\prime}$ are no longer independent, and can be counted as
follows (note that $\Delta F_{AB}(x_{\mathbf{e}},x_{\mathbf{e}}^{\prime})$
will only be modified by the transformations on the boundary, so the
transformations in the bulk for $x_{\mathbf{e}}$ and $x_{\mathbf{e}}^{\prime}$
are not subjected to any constraint). We first note that the disorder-free
configuration for both $x_{\mathbf{e}}$ and $x_{\mathbf{e}}^{\prime}$ of
course satisfy $\Delta F_{AB}(x_{\mathbf{e}},x_{\mathbf{e}}^{\prime})=0$.
Starting from the disorder-free configuration, there are $2^{|\partial A|}$
transformations (including the trivial transformation) that can be applied to
$x_{\mathbf{e}}$ so that $f_{p}(x_{\mathbf{e}})=1,\forall p$. Crucially, given
a specific boundary transformation $l_{\mathbf{v}}$ applied to
$x_{\mathbf{e}}$, there are two boundary transformations
$l^{\prime}_{\mathbf{v}}$ for $x^{\prime}_{\mathbf{e}}$ that make the swapped
configurations $(x_{A}^{\prime},x_{B})$ and $(x_{A},x_{B}^{\prime})$ remain in
the vacuum sector: either $l^{\prime}_{\mathbf{v}}=l_{\mathbf{v}}$ or
$l^{\prime}_{\mathbf{v}}=\partial A+l_{\mathbf{v}}\ (\text{mod}\ 2)$. We
provide one such example in Fig.2. As a consequence, there are $2^{|\partial
A|}\times 2$ transformations on the boundary that satisfy $\Delta
F_{AB}(x_{\mathbf{e}},x_{\mathbf{e}}^{\prime})=0$. It follows that

$\displaystyle S_{2}$ $\displaystyle=-\log(\frac{2^{|\partial
A|+1}}{2^{2|\partial A|}})=-\log(\frac{1}{2^{|\partial A|-1}})$ (9)
$\displaystyle=(|\partial A|-1)\log 2,$

which is precisely the entanglement entropy of the 2d toric code ground state,
as expected.

Given the aforementioned understanding in the limit $\beta=\infty$, let’s now
consider the whole ferromagnetic phase $\beta>\beta_{c}$, where $\beta_{c}$ is
the critical temperature of the RBIM along the Nishimori line. Since now
$\beta$ is no longer infinite, all flux configurations
$f_{\mathbf{p}}{(x_{\mathbf{e}})}$ and
$f_{\mathbf{p}}{(x_{\mathbf{e}}^{\prime})}$ contribute in Eq.(7). However,
crucially, the free energy cost of changing a single flux (i.e.,
$f_{p}\rightarrow-f_{p}$ for a given $p$) scales with the system size in the
ferromagnetic phase [30]. Therefore, for a given
$(f_{\mathbf{p}},f^{\prime}_{\mathbf{p}})$, it is reasonable to expect that
the numerator in Eq.(7) is dominated by the configurations that minimize
$\Delta F_{AB}(x_{\mathbf{e}},x_{\mathbf{e}}^{\prime})$, where
$(x_{\mathbf{e}},x_{\mathbf{e}}^{\prime})\in(f_{\mathbf{p}},f^{\prime}_{\mathbf{p}})$.
Similar to the counting of the vacuum sector in the limit $\beta=\infty$, all
such configurations can be obtained by applying transformations akin to Fig.2
to a configuration
$(x_{\mathbf{e},0},x_{\mathbf{e},0}^{\prime})\in(f_{\mathbf{p}},f^{\prime}_{\mathbf{p}})$
that minimizes $\Delta F_{AB}$. It follows that there are again $2^{|\partial
A|+1}$ boundary transformations for each
$(f_{\mathbf{p}},f^{\prime}_{\mathbf{p}})$, leading to a $-\log 2$ subleading
correction. We note that since the minimized free energy $\Delta
F_{AB}(x_{\mathbf{e}}\in f_{\mathbf{p}},x_{\mathbf{e}}^{\prime}\in
f^{\prime}_{\mathbf{p}})$ depends on the flux sector
$(f_{\mathbf{p}},f^{\prime}_{\mathbf{p}})$, the leading boundary-law
coefficient will be modified. Therefore, we expect the second Renyi entropy
$S_{2}=\alpha|\partial A|-\log 2$ in the ferromagnetic phase, with $\alpha$ a
non-universal constant. On the other hand, when $\beta<\beta_{c}$, so that the
system is in the paramagnetic phase, the free energy cost of changing a single
flux is finite. Therefore, there is no non-local constraint on the
configurations in Eq.(7), and therefore, one expects a purely boundary-law
scaling of the Renyi entropy without any subleading TEE.

##  Appendix B A non-optimal decomposition

Here we provide details of the non-optimal decompositions for the 2d and 3d
toric codes subjected to the phase-flip errors whose results are briefly
mentioned in the main text.

## 2d toric code

Figure 3:  The correspondence between the components of the (non-optimal)
convex decomposition for toric codes, i.e.,
$[\tanh(\beta/2)]^{|E^{\prime}|}\prod_{e\in
E^{\prime}}X_{e}|z_{\mathbf{e}}\rangle$ (see Eq.(14)) and the domain wall
configurations of the Ising model with the Ising interaction determined by
$\\{z_{e}\\}$ in (a) 2d and (b) 3d. Note that the classical spins reside on
plaquettes for case (a) while they reside on vertices for case (b). The edges
with $z_{e}=-1$ of the initial state $|z_{\mathbf{e}}\rangle$ label the chain
$E$ of antiferromagnetic bonds (red). The edges where Pauli-$X$s are applied
to correspond to the chain $E^{\prime}$ of excited bonds (blue), which is
determined by $E^{\prime}=D+E$ ($\textrm{mod}\ 2$) with $D$ labeling the
domain walls (the boundary of the shaded regions). The squared norm of the
amplitude, namely $[\tanh^{2}(\beta/2)]^{|E^{\prime}|}$, corresponds to the
probability for such a domain wall configuration. This notation is directly
borrowed from Ref.[30].

To obtain an alternative convex decomposition of the decohered toric code
density matrix in Eq.(1), we again return to the bigger Hilbert space of the
2d cluster state, and write the decohered state $\rho$ as
$\rho\propto\Gamma\Gamma^{\dagger}$, where

$\Gamma=\langle x_{\mathbf{v}}=1|e^{-\beta\sum_{e}h_{e}/2}\prod_{v}(I-h_{v})$
(10)

(recall $h_{v}=-X_{v}\prod_{e\ni v}Z_{e}$ and $h_{e}=-X_{e}\prod_{v\in
e}Z_{v}$ for the 2d cluster Hamiltonian). We note that unlike the operator
$\rho^{1/2}$ that was used for the optimal decomposition in the main text,
$\Gamma$ in Eq.(10) is not a square matrix, but instead is of size
$2^{N_{e}}\times 2^{N_{e}+N_{v}}$, where $N_{e}\,(N_{v})$ denotes the number
of edges (vertices) of the square lattice on which the cluster state lives. To
proceed, we insert the complete basis
$\\{|m\rangle=|z_{\mathbf{e}},x_{\mathbf{v}}\rangle\\}$ between $\Gamma$ and
$\Gamma^{\dagger}$, where $|z_{\mathbf{e}}\rangle=\otimes_{e}|z_{e}\rangle$
($|x_{\mathbf{v}}\rangle=\otimes_{v}|x_{v}\rangle$) denotes the eigenvector of
Pauli-$Z_{e}$ (Pauli-$X_{v}$) with eigenvalue $z_{e}\,(x_{v})=\pm 1$. The
corresponding $|\psi_{m}\rangle$ for a gvien $m$ can then be written as
$|\psi_{m}\rangle\propto\Gamma|z_{\mathbf{e}},x_{\mathbf{v}}\rangle$. Since
$\prod_{v}(I-h_{v})$ is the projector to the even sector of $X_{v}(\prod_{e\ni
v}Z_{e})$, $|\psi_{m}\rangle$ is non-vanishing only when $x_{v}=\prod_{e\in
v}z_{e}$. Therefore, $x_{\mathbf{v}}$ is completely determined by
$z_{\mathbf{e}}$. It is then more natural to label $|\psi_{m}\rangle$ using
$z_{\mathbf{e}}$, and each $|\psi_{z_{\mathbf{e}}}\rangle$ can be written as

$\displaystyle|\psi_{z_{\mathbf{e}}}\rangle\propto\langle
x_{\mathbf{v}}=1|e^{-\beta\sum_{e}h_{e}/2}|z_{\mathbf{e}},x_{\mathbf{v}}=\prod_{e\in\mathbf{v}}z_{e}\rangle.$
(11)

We emphasize that unlike the optimal decomposition mentioned in the main text
(namely, $\rho^{1/2}|z_{\mathbf{e}}\rangle,\forall z_{\mathbf{e}}$), the
states $|\psi_{z_{\mathbf{e}}}\rangle$ that enter the convex decomposition
here are not related to each other via a product of single-site unitaries.
This implies that one should analyze all the $|\psi_{z_{\mathbf{e}}}\rangle$
to determine whether
$\rho=\sum_{z_{\mathbf{e}}}p_{z_{\mathbf{e}}}|\psi_{z_{\mathbf{e}}}\rangle\langle\psi_{z_{\mathbf{e}}}|$
is a convex sum of SRE states. Besides, since the number of states
$|\psi_{z_{\mathbf{e}}}\rangle$ grows exponentially as a function of the
system size, the probability for each $|\psi_{z_{\mathbf{e}}}\rangle$ is
exponentially small in general. It follows that even if there exists a subset
of $|\psi_{z_{\mathbf{e}}}\rangle$ that are not SRE, the decohered state
$\rho$ may still be well approximated by an SRE mixed state, as long as the
total probability of the LRE pure states is vanishingly small. Therefore, the
notion of $\rho$ being a SRE mixed state must take into account the
probability for each $|\psi_{z_{\mathbf{e}}}\rangle$, and can only be made
precise in a statistical sense. To make progress, we characterize the ensemble
$\\{|\psi_{z_{\mathbf{e}}}\rangle\\}$ by the following average anyon
condensation parameter:

$[\langle
T_{\tilde{\mathcal{C}}}(\beta)\rangle^{2}]\equiv\sum_{z_{\mathbf{e}}}p_{z_{\mathbf{e}}}(\beta)\langle
T_{\tilde{\mathcal{C}}}(\beta)\rangle_{z_{\mathbf{e}}}^{2},$ (12)

where $T_{\tilde{\mathcal{C}}}=\prod_{e\in\tilde{\mathcal{C}}}Z_{e}$ with
$\tilde{\mathcal{C}}$ an open curve on the dual lattice, and $\langle
T_{\tilde{\mathcal{C}}}(\beta)\rangle_{z_{\mathbf{e}}}\equiv\langle\psi_{z_{\mathbf{e}}}|T_{\tilde{\mathcal{C}}}|\psi_{z_{\mathbf{e}}}\rangle/\langle\psi_{z_{\mathbf{e}}}|\psi_{z_{\mathbf{e}}}\rangle$.
Intuitively, $T_{\tilde{\mathcal{C}}}$ can be understood as an operator that
creates a pair of anyons on the plaquettes $p$ and $p^{\prime}$, the end
points of the curve $\tilde{\mathcal{C}}$. Therefore, the expectation value of
$T_{\tilde{\mathcal{C}}}$ vanishes for a topologically ordered (LRE) state as
$|\tilde{\mathcal{C}}|\rightarrow\infty$ and is nonzero for a non-
topologically ordered/anyon condensed (SRE) state. We note that a similar
order parameter has also been employed in Ref.[15] to understand the decohered
mixed state using C-J map/double state formalism. Our main result is the
following equality:

$[\langle T_{\tilde{\mathcal{C}}}(\beta)\rangle^{2}]=[\langle
s_{p}s_{p^{\prime}}(K)\rangle]_{\text{2d RBIM}},\ K=-\log\tanh(\beta/2),$ (13)

where $[\langle s_{p}s_{p^{\prime}}(K)\rangle]_{\text{2d RBIM}}$ denotes the
disorder-averaged spin-spin correlation function of the 2d RBIM along the
Nishimori line at inverse temperature $K$. Crucially, due to the relation
$e^{-K}=\tanh(\beta/2)$, the low $K$ limit (paramagnetic phase) now
corresponds to the large $\beta$ (low error rate) limit, in contrast to the
discussion of the optimal decomposition in the main text where the low error
rate corresponds to the ferromagnetic phase. This is a consequence of the
Kramers-Wannier duality. Therefore, $[\langle
T_{\tilde{\mathcal{C}}}\rangle^{2}]$ vanishes in the low error rate limit,
consistent with our expectation. To determine the critical error rate
$p_{\textrm{non-optimal}}$ i.e. the error rate beyond which the decohered
density matrix is (statistically) SRE, we use $\tanh(\beta_{\textrm{non-
optimal}})=1-2p_{\textrm{non-optimal}}$ and $\tanh^{2}(\beta_{\textrm{non-
optimal}}/2)=e^{-2K_{\textrm{2d RBIM}}}=p_{\textrm{2d RBIM}}/(1-p_{\textrm{2d
RBIM}})$ with $p_{\textrm{2d RBIM}}\approx 0.109$ to obtain $p_{\textrm{non-
optimal}}\approx 0.188$. As expected, $p_{\textrm{non-
optimal}}>p_{\textrm{optimal}}=p_{\textrm{2d RBIM}}$.

The sketch of our calculations that verify Eq.(13) is as follows. We will
first compute the probabilty $p_{z_{\mathbf{e}}}$ for each
$|\psi_{z_{\mathbf{e}}}\rangle$ to occur, which is simply proportional to the
normalization of $|\psi_{z_{\mathbf{e}}}\rangle$. We then evaluate the
expectation value of anyon condensation operator $\langle
T_{\tilde{\mathcal{C}}}\rangle_{z_{\mathbf{e}}}$. It will then become
transparent that Eq.(12) is equivalent to the 2d RBIM along the Nishimori
line.

To compute the probability
$p_{z_{\mathbf{e}}}\propto\langle\psi_{z_{\mathbf{e}}}|\psi_{z_{\mathbf{e}}}\rangle$,
we first simplify the expression of $|\psi_{z_{\mathbf{e}}}\rangle$. It is
illuminating to denote the edges where $z_{e}=-1$ as the string configuration
$E$ so that $x_{v}=\prod_{e\in v}z_{e}=-1$ corresponds to the end points of
$E$ (see Fig.3(a)). Now, expand
$e^{-\beta\sum_{e}h_{e}/2}\propto\prod_{e}(I-\tanh(\beta/2)h_{e})$ as a sum of
several string operators generated by product of $h_{e}$s, which we denote as
$E^{\prime}$. For $E^{\prime}$ to survive after the projection $\langle
x_{\mathbf{v}}=1|$, its end points should coincide with the end points of $E$
($x_{v}=-1$) so that $h_{e}$ flip those spins to $x_{v}=+1$. This implies that
$E^{\prime}+E$ ($\textrm{mod}\ 2$) should form closed loop configurations $D$,
and $|\psi_{z_{\mathbf{e}}}\rangle$ for a given ${z_{\mathbf{e}}}$ (which
fully specifies $E$) can be written as

$\displaystyle|\psi_{z_{\mathbf{e}}}\rangle\propto\sum_{E^{\prime}\textrm{s.t.}E^{\prime}+E=D}[\tanh(\beta/2)]^{|E^{\prime}|}\prod_{e\in
E^{\prime}}X_{e}|z_{\mathbf{e}}\rangle.$ (14)

Since each $\prod_{e\in E^{\prime}}X_{e}|z_{\mathbf{e}}\rangle$ is orthogonal
to one another, the probability
$p_{z_{\mathbf{e}}}\propto\langle\psi_{z_{\mathbf{e}}}|\psi_{z_{\mathbf{e}}}\rangle$
can be evaluated straightforwardly as

$p_{z_{\mathbf{e}}}\propto\sum_{E^{\prime}\textrm{s.t.}E^{\prime}+E=D}[\tanh^{2}(\beta/2)]^{|E^{\prime}|},$
(15)

which is nothing but the partition function of 2d Ising model
$\mathcal{Z}_{\textrm{2d Ising},{z_{\mathbf{e}}}}(K)$ at inverse temperature
$K=-\log\tanh(\beta/2)$ and the negative Ising interactions labeled by
$E^{\prime}$ (determined by ${z_{\mathbf{e}}}$). Besides, there is a one-to-
one mapping between the components of $|\psi_{z_{\mathbf{e}}}\rangle$, i.e.,
$[\tanh(\beta/2)]^{|E^{\prime}|}\prod_{e\in
E^{\prime}}X_{e}|z_{\mathbf{e}}\rangle$, and the domain wall configurations of
a statistical mechanics model. Specifically, the classical spins reside on the
plaquettes of the square lattice, and the bonds where $z_{e}=-1$ in the
initial state $|z_{\mathbf{e}}\rangle$ label the chain $E$ of the
antiferromagnetic bonds. The end points of $E$ labeled by $x_{v}=\prod_{e\in
v}z_{e}=-1$ then correspond to the Ising vortices. The edges where Pauli-$X$s
are applied correspond to the chain $E^{\prime}$ of excited bonds, which is
determined by $E^{\prime}=D+E$ ($\textrm{mod}\ 2$). The squared norm of the
amplitude $[\tanh^{2}(\beta/2)]^{|E^{\prime}|}=e^{-2K|E^{\prime}|}$
corresponds to the probability for such a domain wall configuration to appear.
We summarize this correspondence in Fig.3(a).

Now, we are ready to compute $\langle
T_{\tilde{\mathcal{C}}}\rangle_{z_{\mathbf{e}}}$ for a given $z_{\mathbf{e}}$.
Using $T_{\tilde{\mathcal{C}}}(\prod_{e\in
E^{\prime}}X_{e})=(-1)^{\\#\textrm{cross}(E^{\prime},\tilde{\mathcal{C}})}(\prod_{e\in
E^{\prime}}X_{e})T_{\tilde{\mathcal{C}}}$ and
$T_{\tilde{\mathcal{C}}}|z_{\mathbf{e}}\rangle=(-1)^{\\#\textrm{cross}(E,\tilde{\mathcal{C}})}|z_{\mathbf{e}}\rangle$
with $\\#\textrm{cross}(A,B)$ the number of crossings between $A$ and $B$, we
get

$\displaystyle
T_{\tilde{\mathcal{C}}}|\psi_{z_{\mathbf{e}}}\rangle=\sum_{E^{\prime}}(-1)^{\\#\textrm{cross}(D,\tilde{\mathcal{C}})}e^{-K|E^{\prime}|}\prod_{e\in
E^{\prime}}X_{e}|z_{\mathbf{e}}\rangle,$ (16)

where we have used $D=E+E^{\prime}$. It follows that

$\langle
T_{\tilde{\mathcal{C}}}(\beta)\rangle_{z_{\mathbf{e}}}=\frac{\sum_{E^{\prime}}(-1)^{\\#\textrm{cross}(D,\tilde{\mathcal{C}})}e^{-2K|E^{\prime}|}}{\mathcal{Z}_{\textrm{2d
Ising},{z_{\mathbf{e}}}}(K)},$ (17)

which corresponds to the spin-spin correlation function $\langle
s_{p}s_{p^{\prime}}(K)\rangle$ of the 2d Ising model with the Ising
interaction determined by $\\{z_{e}\\}$, where $p$ and $p^{\prime}$ denote the
end points of $\tilde{\mathcal{C}}$. It follows that we can characterize the
separability transition for the mixed state
$\rho=\sum_{z_{\mathbf{e}}}p_{z_{\mathbf{e}}}|\psi_{z_{\mathbf{e}}}\rangle\langle\psi_{z_{\mathbf{e}}}|$
using Eq.(12):

$\displaystyle(\beta)$
$\displaystyle=\sum_{z_{\mathbf{e}}}p_{z_{\mathbf{e}}}\langle
T_{\tilde{\mathcal{C}}}\rangle^{2}_{z_{\mathbf{e}}}$ (18)
$\displaystyle=\frac{\sum_{z_{\mathbf{e}}}\mathcal{Z}_{\text{2d
Ising},{z_{\mathbf{e}}}}(K)\langle
s_{p}s_{p}^{\prime}(K)\rangle}{\sum_{z_{\mathbf{e}}}\mathcal{Z}_{\text{2d
Ising},{z_{\mathbf{e}}}}(K)}.$

Therefore, as claimed in Eq.(13), $[\langle
T_{\tilde{\mathcal{C}}}\rangle^{2}](\beta)$ precisely corresponds to the
disorder-averaged spin-spin correlation function of the RBIM along the
Nishimori line with inverse temperature $K=-\log\tanh(\beta/2)$.

## 3d toric code

Here we derive a non-optimal decomposition for the density matrix of the 3d
toric code subjected to phase-flip errors. Similar to the 2d toric code, we
write $\rho\propto\Gamma\Gamma^{\dagger}$ with

$\Gamma=\langle x_{\mathbf{f}}=1|e^{-\beta\sum_{e}h_{e}/2}\prod_{f}(I-h_{f})$
(19)

(recall $h_{e}=-X_{e}\prod_{f\ni e}Z_{f}$ and $h_{f}=-X_{f}\prod_{e\in
f}Z_{e}$ for the 3d cluster Hamiltonian). This equation motivates the
decomposition $\rho=\sum_{m}|\psi_{m}\rangle\langle\psi_{m}|$, where
$|\psi_{m}\rangle\propto\Gamma|z_{\mathbf{e}},x_{\mathbf{f}}\rangle$ and
$|z_{\mathbf{e}}\rangle=\otimes_{e}|z_{e}\rangle$
($|x_{\mathbf{f}}\rangle=\otimes_{f}|x_{f}\rangle$) denotes the eigenvector of
Pauli-$Z_{e}$ (Pauli-$X_{f}$) with eigenvalue $z_{e}(x_{f})=\pm 1$. Since
$\prod_{f}(I-h_{f})$ is the projector to the even sector of $X_{f}(\prod_{e\in
f}Z_{e})$, $|\psi_{m}\rangle$ is non-vanishing only when $x_{f}=\prod_{e\in
f}z_{e}$. Therefore, $x_{\mathbf{f}}$ is completely determined by
$z_{\mathbf{e}}$, and it is more natural to label $|\psi_{m}\rangle$ using
$z_{\mathbf{e}}$, which can be written as

$\displaystyle|\psi_{z_{\mathbf{e}}}\rangle\sim\langle
x_{\mathbf{f}}=1|e^{-\beta\sum_{e}h_{e}/2}|z_{\mathbf{e}},x_{\mathbf{f}}=\prod_{e\in\mathbf{f}}z_{e}\rangle.$
(20)

To statistically characterize the separability condition for $\rho$, one can
use the anyon condensation operator which in this case is a ‘Wilson string
operator’ $W_{\mathcal{C}}$ where the open curve $\mathcal{C}$ lives on the
original lattice (note that unlike a proper $\mathbb{Z}_{2}$ theory,
$W_{\mathcal{C}}$ here does not vanish identically since the Gauss law is
being enforced only energetically, and not exactly). Specifically, our main
result is:

$[\langle W_{\mathcal{C}}(\beta)\rangle^{2}]=[\langle
s_{p}s_{p^{\prime}}(K)\rangle]_{\text{3d RBIM}},\ K=-\log\tanh(\beta/2),$ (21)

where $[\langle s_{p}s_{p^{\prime}}(K)\rangle]_{\text{3d RBIM}}$ denotes the
disorder-averaged spin-spin correlation function of the 3d RBIM along the
Nishimori line with inverse temperature $K$. We note that the 3d Ising model
is dual to the 3d Ising plaquette gauge model, consistent with our expectation
that the non-optimal state is related to the optimal one via a Kramers-Wannier
duality. To determine the bound on the error rate $p_{\textrm{non-optimal}}$
for separability of $\rho$ resulting from this calculation, we use
$\tanh(\beta_{\textrm{non-optimal}})=1-2p_{\textrm{non-optimal}}$ and
$\tanh^{2}(\beta_{\textrm{non-optimal}}/2)=e^{-2K_{\textrm{3d
RBIM}}}=p_{\textrm{3d RBIM}}/(1-p_{\textrm{3d RBIM}})$ with $p_{\textrm{3d
RBIM}}\approx 0.233$ [53] to obtain $p_{\textrm{non-optimal}}\approx 0.077$.
Therefore, the outcome of this analysis is that the decohered density matrix
is SRE for $p>p_{\textrm{non-optimal}}$. As expected, this non-optimal bound
is higher than the optimal threshold $p_{\text{3d RPGM}}\approx 0.029$ for
separability argued in the main text.

To verify Eq.(21), we will follow the same derivation as for the 2d toric code
by first computing the probabilty $p_{z_{\mathbf{e}}}$ and then evaluating
$\langle W_{{\mathcal{C}}}\rangle_{z_{\mathbf{e}}}$. Similar to the 2d toric
code, we denote the edges with $z_{e}=-1$ as $E$. Expanding
$e^{-\beta\sum_{e}h_{e}/2}|z_{\mathbf{e}},x_{\mathbf{f}}=\prod_{e\in\mathbf{f}}z_{e}\rangle\propto\prod_{f}(I-\tanh(\beta/2)h_{e})|z_{\mathbf{e}},x_{\mathbf{f}}=\prod_{e\in\mathbf{f}}z_{e}\rangle$
yields a superposition of terms in the $Z$-basis, and we denote the
corresponding configurations as $E^{\prime}$. To see what subset of
$E^{\prime}$ survives after the projection $\langle x_{f}=1|$, we note that
when one applies $h_{e}$, all the spins $x_{f}$ with $f\ni e$ will be flipped.
Therefore, $h_{e}$ acts non-trivially only on faces $f$, where $x_{f}=-1$ with
$e\in f$. The condition $x_{f}=-1$ is satisfied precisely on faces that are
bounded by edges $e$ with an odd number of $z_{e}=-1$ (since
$x_{f}=\prod_{e\in f}z_{e}$). This implies that $E+E^{\prime}$ should form
closed surfaces (denoted as $D$). Therefore, $|\psi_{z_{\mathbf{e}}}\rangle$
for a given $z_{\mathbf{e}}$ (which fully specifies $E$) has the same form as
Eq.(14), though here $D$ represents closed surfaces instead of closed loops.
Since each $\prod_{e\in E^{\prime}}X_{e}|z_{\mathbf{e}}\rangle$ is orthogonal
to one another,
$p_{z_{\mathbf{e}}}\propto\langle\psi_{z_{\mathbf{e}}}|\psi_{z_{\mathbf{e}}}\rangle$
can be computed as
$\sum_{E^{\prime}\textrm{s.t.}E^{\prime}+E=D}[\tanh^{2}(\beta/2)]^{|E^{\prime}|}$,
which is nothing but the partition function of 3d Ising model
$\mathcal{Z}_{\textrm{3d Ising}}(K)$ with inverse temperature
$K=-\log\tanh(\beta/2)$ and the negative Ising interactions labeled by
$E^{\prime}$ (determined by ${z_{\mathbf{e}}}$). Therefore, similar to the 2d
toric code, each component of the non-optimal $|\psi_{z_{\mathbf{e}}}\rangle$
has a one-to-one correspondence with the domain-wall configuration of the 3d
Ising model, which we summarize in Fig.3(b). A derivation similar to that for
the 2d toric code then shows that the average Wilson string operator (Eq.(21))
is equivalent to the disorder-averaged correlation function of the 3d RBIM
along the Nishimori line with inverse temperature $K=-\log\tanh(\beta/2)$.


# A Survey on Offline Model-Based Reinforcement Learning

Haoyang He
Electrical and Computer Engineering
Carnegie Mellon University
Pittsburgh, PA 15213
<EMAIL_ADDRESS>

###### Abstract

Model-based approaches are becoming increasingly popular in the field of
offline reinforcement learning, with high potential in real-world applications
due to the model’s capability of thoroughly utilizing the large historical
datasets available with supervised learning techniques. This paper presents a
literature review of recent work in offline model-based reinforcement
learning, a field that utilizes model-based approaches in offline
reinforcement learning. The survey provides a brief overview of the concepts
and recent developments in both offline reinforcement learning and model-based
reinforcement learning, and discuss the intersection of the two fields. We
then presents key relevant papers in the field of offline model-based
reinforcement learning and discuss their methods, particularly their
approaches in solving the issue of distributional shift, the main problem
faced by all current offline model-based reinforcement learning methods. We
further discuss key challenges faced by the field, and suggest possible
directions for future work.

##  1 Introduction

Reinforcement learning (RL) constitutes a subfield of machine learning,
concentrating on the development of agents capable of making optimal decisions
within uncertain environments. In an RL setting, agents engage with their
surroundings by performing actions and receiving feedback in the form of
rewards or penalties. The primary objective is for the agent to ascertain a
policy, a mapping from states to actions, which maximizes the accumulated
reward over a specified time period [14].

###  1.1 Offline Reinforcement Learning

Offline reinforcement learning, also referred to as batch reinforcement
learning, entails a form of reinforcement learning wherein the policy is
exclusively learned based on pre-existing historical data. In practical
applications, the available historical dataset is typically extensive, and
offline reinforcement learning seeks to maximize the utilization of this data,
learning the policy as a decision-making engine, contrasting with online
reinforcement learning’s emphasis on exploration [11]. In practice, the
preceding data collection policy, responsible for generating the historical
data, is customarily the current method or model in deployment. Recent
research endeavors often concentrate on optimizing the algorithm’s outcome in
relation to the data collection policy, with the aim of devising new
algorithms possessing practical deployment value [2, 10]. The algorithm
grapples with the challenge of distributional shift, signifying that the
deployment scenario exhibits a distinct distribution in comparison to the
historical dataset, with respect to state and action space [7]. The presence
of distributional shift constitutes a key obstacle in offline reinforcement
learning, warranting further in-depth examination.

###  1.2 Model-Based Reinforcement Learning

Model-based reinforcement learning encompasses a category of reinforcement
learning methodologies that focus on constructing a model to comprehend the
environment by approximating state transition dynamics and reward functions.
By leveraging the estimated model, these algorithms typically generate
simulated data to facilitate the learning of an optimal policy, in conjunction
with the actual data, also known as experienced data. Enhancing the
performance of the learned policy by capitalizing on the model necessitates
acquiring a model that more precisely characterizes the actual environment
[9]. Numerous techniques for augmenting model learning have been investigated
in the literature over time, including the use of approximation functions
trained with prediction loss, which is susceptible to horizontal compounding
error during updates. In order to mitigate this error, research has proposed
implementing constraints on the model, such as the Lipschitz Continuity
Constraint [1]. An alternative approach to diminishing this error involves
distribution matching via adversarial learning against a discriminator [17].

##  2 Offline Model-Based Reinforcement Learning

The intersection of the two previously discussed concepts defines offline
model-based reinforcement learning. Analogous to online model-based RL
methods, offline model-based RL approaches initially estimate state transition
dynamics and reward functions, subsequently employing them for environment
simulation and planning based on both simulated data and experienced data.
Owing to the characteristics of model-based learning, offline model-based
learning may particularly benefit from the robust model established by
supervised learning, which takes full advantage of the extensive historical
dataset when compared to offline model-free methods [11]. Nevertheless, due to
the intrinsic properties of offline learning, offline model-based RL also
grapples with distributional shifts in both state and action distributions.
Unlike online methods, no exploration is available for the model to rectify
itself. This challenge is exacerbated by model exploitation since the model is
derived from data collected by another policy with a distributional shift
relative to the deployment environment. Consequently, simple policy
optimization over the model may result in an optimized outcome in the
historical distribution, which is out of distribution for the deployment
environment, potentially leading to inferior results. Recent advancements in
offline model-based reinforcement learning have concentrated on addressing the
issue of distributional shift, typically by imposing certain constraints on
the model through modifications in state transition dynamics, reward
functions, or value functions [8, 6, 19, 10, 18, 13, 2].

###  2.1 MOReL

Kidambi et al. proposed MOReL [6], a framework for model-based offline
reinforcement learning that penalizes the reward for uncertainty, allowing the
use of classical planning methods with altered reward function and state
transition dynamics to address the issue of uncertainty and distributional
shift [7]. They constructed an Unknown State Action Detector based on a
threshold over the total variational divergence between the learned model and
the behavior model, which is then used to determine if a certain policy goes
into the absorbing state HALT to act as a control for the policy to enter the
above-threshold uncertain state, whose reward is set to a certain negative
value, addressing the model exploitation problem. In their practical
implementation, they constructs an ensemble of state transition dynamics
models for uncertainty estimation, using a measurement of disagreement between
the models given by eq.1,

$dist_{i,j}(s_{t},a_{t})=\mathbb{E}_{s_{t+1}^{i}\sim
T_{\psi_{i}}(\cdot|s_{t},a_{t}),s_{t+1}^{i}\sim
T_{\psi_{i}}(\cdot|s_{t},a_{t})}[||s^{i}_{t+1}-s^{j}_{t+1}||].$ (1)

in order to define the uncertainty of the state action pair as eq.2.

$U_{r}(s,a)=\left\\{\begin{aligned}
r_{max},\max_{i,j}dist_{i,j}(s,a)>threshold,\\\
0,\max_{i,j}dist_{i,j}(s,a)\leq threshold.\end{aligned}\right.$ (2)

###  2.2 MOPO

Yu et al. proposed MOPO [19], model-based offline policy optimization, a
similar attempt to address the problem of distributional shift by penalizing
the reward. In the practical implementation of MOPO, the ensemble of models
are estimated using the multivariate Gaussian distributions modeling the $i$th
transition dynamics of the model ensembles, denotes in eq.3,

$T_{\psi_{i}}(s^{\prime}|s,a)=\mathcal{N}(\mu_{i}(s,a),\Sigma_{i}(s,a)).$ (3)

while the uncertainty is estimated by the Frobenius norm of the maximum
distributional variance in among the ensemble of models, as shown in eq.4.

$U_{r}(s,a)=\max_{i}||\Sigma_{i}(s,a)||_{F}.$ (4)

Both methods are designed to build upon previously existing online model-based
RL methods, providing a penalized reward function, which can in turn be
directly used instead of the original reward function in classical planning
methods [11] like linear-quadratic regulator [15] or Monte Carlo tree search
[3].

###  2.3 BREMEN

An alternative approach to the distributional shift problem proposed by
Matsushima et al. is BREMEN [10], behavior-regularized model ensemble. Instead
of constructing a pessimistic or reward-penalized MDP model, BREMEN learns an
ensemble of dynamic models in conjunction with a policy using imaginary
rollouts from the ensemble and behavior regularization via conservative trust-
region updates. Each model of the ensemble is trained to optimize the
objective of eq.5, the mean squared error between the prediction of the next
state and the true next state, and imaginary model rollouts are generated
sequentially.

$\min_{\phi_{i}}\frac{1}{|\mathcal{D}|}\sum_{(s_{t},a_{t},s_{t+1})\in\mathcal{D}}\frac{1}{2}||s_{t+1}-\hat{f}_{\phi_{i}}(s_{t},a_{t})||^{2}_{2}$
(5)

To address the distributional shift problem, BREMEN uses an iterative policy
update via a trust-region constraint, re-initialized with the behavior cloned
policy after every deployment, with the behavior cloning objective function of
eq.6.

$\min_{\beta}\frac{1}{|\mathcal{D}|}\sum_{(s_{t},a_{t})\in\mathcal{D}}\frac{1}{2}||a_{t}-\hat{\pi}_{\beta}(s_{t})||^{2}_{2}$
(6)

After obtaining the estimated behavior policy, the target policy
$\pi_{\theta}$ is initialized as a Gaussian policy with mean from the behavior
policy and standard deviation of $1$. This behavior cloning initialization may
be seen as implicitly biasing the optimized target policy to be close to the
data-collection policy, and thus reduces the distributional shift. To further
bias the learned policy to be close to the data-collection policy, KL-
divergence trust region under a certain threshold is added as an additional
constraint, yielding the optimization problem of eq.7.

$\displaystyle\theta_{k+1}=\underset{\theta}{\arg\max}\underset{s,a\sim\pi_{\theta_{k}},\hat{f}_{\phi_{i}}}{\mathbb{E}}\left[\frac{\pi_{\theta}(a\mid
s)}{\pi_{\theta_{k}}(a\mid s)}A^{\pi_{\theta_{k}}}(s,a)\right],$ (7)
$\displaystyle\text{ s.t.
}\quad\underset{s,a\sim\pi_{\theta_{k}},\hat{f}_{\phi_{i}}}{\mathbb{E}}\left[D_{\mathrm{KL}}\left(\pi_{\theta}(\cdot\mid
s)\|\pi_{\theta_{k}}(\cdot\mid
s)\right)\right]\leq\delta,\quad\pi_{\theta_{0}}=\operatorname{Normal}\left(\hat{\pi}_{\beta},1\right).$

The authors claim an empirical improvement over direct KL-divergence penalties
such as previous methods, and showed that KL-divergence is controlled
implicitly in the optimization. The author also highlighted that BREMEN may
outperform similar methods in offline settings with limited data, attributing
to its sequential imaginary behavioral model rollouts, hence having higher
sample efficiency.

###  2.4 VI-LCB

VI-LCB [12], offline value iteration with lower confidence bound, is an
alternative approach to addressing the uncertainty created by distributional
shift proposed by Rashidinejad et al. With a similar idea of penalizing state-
action pairs under poor or partial coverage, VI-LCB builds upon the classic
value iteration algorithm, and instead of iterating upon the original values
given by the learned model, VI-LCB employs the concept of pessimism [12] and
subtracts a penalty function to the value function. The value function for
value iteration in VI-LCB is defined by eq.8, 9.

$Q(s,a)\mapsfrom
r(s,a)-b(s,a)+\gamma\sum_{s^{\prime}}P(s^{\prime}|s,a)V(s^{\prime}),\forall(s,a),$
(8) $V(s)\mapsfrom\max_{a}Q(s,a),\forall s.$ (9)

VI-LCB is developed shortly after PEVI [5], pessimistic value iteration
proposed by Jin et al. It is a simple framework to utilize pessimism in
offline reinforcement learning, which simply changes the addition to
subtraction for the term encouraging exploration in online RL algorithms to
use the term to penalize less-explored state-action pairs.

Li et al. proposed a version of VI-LCB with Bernstein-style penalties, denoted
by eq.10.

$b(s,a;V):=\min\\{\max\\{\sqrt{\frac{c_{b}\log\frac{N}{(1-\gamma)\delta}}{N(s,a)}Var_{\hat{P}_{s,a}}(V)},\frac{2c_{b}\log\frac{N}{(1-\gamma)\delta}}{(1-\gamma)N(s,a)}\\},\frac{1}{1-\gamma}\\}+\frac{5}{N}.$
(10)

Based on this version of VI-LCB, Li et al. further proved that model-based
approaches in offline reinforcement learning achieve minimax-optimal sample
complexity without burn-in cost for tabular Markov Decision Process (MDP) [8].
Specifically, given a single-policy clipped concentrability coefficient
$C^{*}_{clipped}$, model-based offline RL yields $\epsilon$-accuracy with
sample complexity given by

$\left\\{\begin{aligned}
\frac{SC^{*}_{clipped}}{(1-\gamma)^{3}\epsilon^{2}},\text{infinite-horizon
MDPs},\\\ \frac{H^{4}SC^{*}_{clipped}}{\epsilon^{2}},\text{finite-horizon
MDPs}.\end{aligned}\right.$ (11)

for $\gamma$-discounted infinite-horizon MDPs with effective horizon
$\frac{1}{1-\gamma}$ or finite-horizon MDPs with horizon $H$.

###  2.5 COMBO

COMBO [18], conservative offline model-based policy optimization, is an
approach proposed by Yu et al. to prevent distributional shift by directly
optimizing the model learning without explicitly providing an uncertainty
estimation in the planning process. In the construction of the model rollouts,
each policy is repeatedly evaluated on the conservative policy estimation
function 12,

$\hat{Q}^{k+1}\leftarrow\arg\min_{Q}\beta\left(\mathbb{E}_{\mathbf{s},\mathbf{a}\sim\rho(\mathbf{s},\mathbf{a})}[Q(\mathbf{s},\mathbf{a})]-\mathbb{E}_{\mathbf{s},\mathbf{a}\sim\mathcal{D}}[Q(\mathbf{s},\mathbf{a})]\right)+\frac{1}{2}\mathbb{E}_{\mathbf{s},\mathbf{a},\mathbf{s}^{\prime}\sim
d_{f}}\left[\left(Q(\mathbf{s},\mathbf{a})-\widehat{\mathcal{B}}^{\pi}\hat{Q}^{k}(\mathbf{s},\mathbf{a})\right)^{2}\right].$
(12)

which penalizes rewards that are more likely to be out of distribution, based
on the sampling distribution, a hyperparameter of choice to determine the
desired conservativeness. The policy is then improved using the conservative
critic Q obtained from previous conservative estimation 13.

$\hat{\pi}\leftarrow\arg\max_{\pi}\mathbb{E}_{s\sim\rho,a\sim\pi(\cdot|s)}\left[\hat{Q}^{\pi}(\mathbf{s},\mathbf{a})\right].$
(13)

The authors also showed that the obtained policy has a safe improvement
guarantee over the behavior policy. A key improvement of COMBO over similar
previous methods like MOReL and MOPO is that it no longer requires an
uncertainty estimation in order to perform optimization, which can in itself
be inaccurate and hard to obtain [11], but instead utilizes an adversarial
training process to optimize itself.

###  2.6 RAMBO

Rigter et al. proposed RAMBO-RL [13], Robust Adversarial Model-Based Offline
Reinforcement Learning, a self-adversarial model. Instead of previous
adversarial RL models that train an adversary, RAMBO directly trains the model
adversarially, while enforcing the conservatism in the training of the
adversarial dynamics model. The algorithm is basically an iterative update of
the optimization of the adversarial model updating the state transition
dynamics, where the optimization objective of the adversarial model is based
on a formulation of optimization objective to reduce distributional shift of a
previous work Uehara et al. [16], modified to eq.14,

$\min_{\widehat{T}_{\phi}}V_{\phi}^{\pi},\quad\text{ s.t.
}\mathbb{E}_{\mathcal{D}}\left[\mathrm{TV}\left(\widehat{T}_{\mathrm{MLE}}(\cdot\mid
s,a),\widehat{T}_{\phi}(\cdot\mid s,a)\right)^{2}\right]\leq\xi,$ (14)

which can be rewritten with Lagrangian relaxation to eq.15.

$\max_{\lambda\geq
0}\min_{\widehat{T}_{\phi}}\left(L(\hat{T},\lambda):=V_{\phi}^{\pi}+\lambda\left(\mathbb{E}_{\mathcal{D}}\left[\mathrm{TV}\left(\widehat{T}_{\mathrm{MLE}}(\cdot\mid
s,a),\widehat{T}_{\phi}(\cdot\mid s,a)\right)^{2}\right]-\xi\right)\right).$
(15)

Rather than optimizing the Lagrangian multiplier, the authors claimed it was
easier to fix the multiplier as a hyperparameter for tuning, reducing it as a
scaling factor, which can then be reduced to eq.16.

$\min_{\widehat{T}_{\phi}}\left(\lambda
V_{\phi}^{\pi}+\mathbb{E}_{\mathcal{D}}\left[\operatorname{TV}\left(\widehat{T}_{\mathrm{MLE}}(\cdot\mid
s,a),\widehat{T}_{\phi}(\cdot\mid s,a)\right)^{2}\right]\right).$ (16)

To make the optimization more simple, the authors changed from optimizing the
total variational divergence to simple MLE loss, which leads to the final
optimization objective eq.17.

$\mathcal{L}_{\phi}=\lambda
V_{\phi}^{\pi}-\mathbb{E}_{\left(s,a,r,s^{\prime}\right)\sim\mathcal{D}}\left[\log\widehat{T}_{\phi}\left(s^{\prime},r\mid
s,a\right)\right].$ (17)

The authors explicitly compared RAMBO to COMBO, and claimed that the
pessimistic value function updates used by COMBO create local maxima in the
Q-function which are present throughout training, making COMBO likely to stuck
in a local maxima, whereas RAMBO’s adversarial modifications resulted in a
less likelihood.

###  2.7 ARMOR

Bhardwaj et al. proposed ARMOR [2], adversarial model for offline
reinforcement learning, a novel framework capable of improving upon an
arbitrary reference policy regardless of data quality. The two-player game
consists of a learner policy and an adversary MDP model, defined as eq.18.

$\hat{\pi}=\arg\max_{\pi\in\Pi}\min_{M\in\mathcal{M}_{\alpha}}J_{M}(\pi)-J_{M}(\pi_{ref}).$
(18)

In the algorithm, the adversary MDP model is updated in the process 19 for
$i=1,2$, where $\mathcal{D}_{M}$ is the transition tuples using model
predictions based on minibatches of dataset samples from real and model data,
$f_{1},f_{2}$ are critic networks used by the double Q residual algorithm loss
from [4].

$\displaystyle l^{\text{adversary }}(f,M)$
$\displaystyle:=\mathcal{L}_{\mathcal{D}_{M}}\left(f,\pi,\pi_{\text{ref
}}\right)+\beta\left(\mathcal{E}_{\mathcal{D}_{M}}^{w}(f,M,\pi)+\lambda\mathcal{E}_{\mathcal{D}_{\text{real
}}^{\text{min }}}(M)\right),$ (19) $\displaystyle M$ $\displaystyle\leftarrow
M-\eta_{\text{fast }}\left(\nabla_{M}l^{\text{adversary
}}\left(f_{1},M\right)+\nabla_{M}l^{\text{adversary
}}\left(f_{2},M\right)\right),$ $\displaystyle f_{i}$
$\displaystyle\leftarrow\operatorname{Proj}_{\mathcal{F}}\left(f_{i}-\eta_{\text{fast
}}\nabla_{f_{i}}l^{\text{adversary }}\left(f_{i},M\right)\right),$
$\displaystyle\bar{f}_{i}$ $\displaystyle\leftarrow(1-\tau)\bar{f}_{i}+\tau
f_{i}.$

Then the actor network is updated with respect to the first critic $f_{1}$ and
the reference policy in the process 20.

$\displaystyle l^{\text{actor }}(\pi)$
$\displaystyle:=-\mathcal{L}_{\mathcal{D}_{M}}\left(f_{1},\pi,\pi_{\text{ref
}}\right),$ (20) $\displaystyle\pi$
$\displaystyle\leftarrow\operatorname{Proj}_{\Pi}\left(\pi-\eta_{\text{slow
}}\nabla_{\pi}l^{\text{actor }}(\pi)\right).$

The model state is updated by query of the MDP model after the updates to the
two players. Inspired by the idea of relative pessimism proposed by Cheng et
al. [4], ARMOR optimizes for the worst-case relative performance over
uncertainty, instead of a predetermined penalty term based on the state-action
itself. The authors showed that with relative pessimism, ARMOR policy does not
degrades the performance of a reference policy when given a fixed set of
hyperparameters, and is capable of competing against any policy covered by the
data when the right hyperparameter is chosen.

##  3 Challenges and Future Work

Distributional shift remains the most influential problem in offline model-
based reinforcement learning, and warrants further research. Even though the
discussed relevant work have proposed various methods to mitigate this
problem, none of the approaches have succeeded in providing a theoretically
provable data non-reliant absolute performance guarantee due to this problem.

Another promising direction is the absolute relative performance of offline
model-based reinforcement learning. As highlighted by [2], many theoretically
proven methods are not deployed in real-world scenarios because they failed to
achieve absolute more optimal performance over the current policy, making
policy replacement infeasible. Future work on the absolute relative
performance of learned policies can greatly improve the deployment of novel
offline model-based reinforcement learning algorithms.

## References

  * [1] Kavosh Asadi, Dipendra Misra, and Michael Littman.  Lipschitz continuity in model-based reinforcement learning.  In International Conference on Machine Learning, pages 264–273. PMLR, 2018.
  * [2] Mohak Bhardwaj, Tengyang Xie, Byron Boots, Nan Jiang, and Ching-An Cheng.  Adversarial model for offline reinforcement learning, 2023.
  * [3] Cameron B Browne, Edward Powley, Daniel Whitehouse, Simon M Lucas, Peter I Cowling, Philipp Rohlfshagen, Stephen Tavener, Diego Perez, Spyridon Samothrakis, and Simon Colton.  A survey of monte carlo tree search methods.  IEEE Transactions on Computational Intelligence and AI in games, 4(1):1–43, 2012.
  * [4] Ching-An Cheng, Tengyang Xie, Nan Jiang, and Alekh Agarwal.  Adversarially trained actor critic for offline reinforcement learning.  In International Conference on Machine Learning, pages 3852–3878. PMLR, 2022.
  * [5] Ying Jin, Zhuoran Yang, and Zhaoran Wang.  Is pessimism provably efficient for offline rl?  In Marina Meila and Tong Zhang, editors, Proceedings of the 38th International Conference on Machine Learning, volume 139 of Proceedings of Machine Learning Research, pages 5084–5096. PMLR, 18–24 Jul 2021.
  * [6] Rahul Kidambi, Aravind Rajeswaran, Praneeth Netrapalli, and Thorsten Joachims.  Morel: Model-based offline reinforcement learning.  In H. Larochelle, M. Ranzato, R. Hadsell, M.F. Balcan, and H. Lin, editors, Advances in Neural Information Processing Systems, volume 33, pages 21810–21823. Curran Associates, Inc., 2020.
  * [7] Sergey Levine, Aviral Kumar, George Tucker, and Justin Fu.  Offline reinforcement learning: Tutorial, review, and perspectives on open problems, 2020.
  * [8] Gen Li, Laixi Shi, Yuxin Chen, Yuejie Chi, and Yuting Wei.  Settling the sample complexity of model-based offline reinforcement learning.  arXiv preprint arXiv:2204.05275, 2022.
  * [9] Fan-Ming Luo, Tian Xu, Hang Lai, Xiong-Hui Chen, Weinan Zhang, and Yang Yu.  A survey on model-based reinforcement learning, 2022.
  * [10] Tatsuya Matsushima, Hiroki Furuta, Yutaka Matsuo, Ofir Nachum, and Shixiang Gu.  Deployment-efficient reinforcement learning via model-based offline optimization, 2020.
  * [11] Rafael Figueiredo Prudencio, Marcos R. O. A. Maximo, and Esther Luna Colombini.  A survey on offline reinforcement learning: Taxonomy, review, and open problems.  IEEE Transactions on Neural Networks and Learning Systems, pages 1–0, 2023.
  * [12] Paria Rashidinejad, Banghua Zhu, Cong Ma, Jiantao Jiao, and Stuart Russell.  Bridging offline reinforcement learning and imitation learning: A tale of pessimism.  In M. Ranzato, A. Beygelzimer, Y. Dauphin, P.S. Liang, and J. Wortman Vaughan, editors, Advances in Neural Information Processing Systems, volume 34, pages 11702–11716. Curran Associates, Inc., 2021.
  * [13] Marc Rigter, Bruno Lacerda, and Nick Hawes.  Rambo-rl: Robust adversarial model-based offline reinforcement learning, 2022.
  * [14] Richard S Sutton and Andrew G Barto.  Reinforcement learning: An introduction.  MIT press, 2018.
  * [15] Yuval Tassa, Tom Erez, and Emanuel Todorov.  Synthesis and stabilization of complex behaviors through online trajectory optimization.  In 2012 IEEE/RSJ International Conference on Intelligent Robots and Systems, pages 4906–4913. IEEE, 2012.
  * [16] Masatoshi Uehara and Wen Sun.  Pessimistic model-based offline reinforcement learning under partial coverage.  arXiv preprint arXiv:2107.06226, 2021.
  * [17] Tian Xu, Ziniu Li, and Yang Yu.  Error bounds of imitating policies and environments for reinforcement learning.  IEEE Transactions on Pattern Analysis and Machine Intelligence, 44(10):6968–6980, 2021.
  * [18] Tianhe Yu, Aviral Kumar, Rafael Rafailov, Aravind Rajeswaran, Sergey Levine, and Chelsea Finn.  Combo: Conservative offline model-based policy optimization.  In M. Ranzato, A. Beygelzimer, Y. Dauphin, P.S. Liang, and J. Wortman Vaughan, editors, Advances in Neural Information Processing Systems, volume 34, pages 28954–28967. Curran Associates, Inc., 2021.
  * [19] Tianhe Yu, Garrett Thomas, Lantao Yu, Stefano Ermon, James Y Zou, Sergey Levine, Chelsea Finn, and Tengyu Ma.  Mopo: Model-based offline policy optimization.  In H. Larochelle, M. Ranzato, R. Hadsell, M.F. Balcan, and H. Lin, editors, Advances in Neural Information Processing Systems, volume 33, pages 14129–14142. Curran Associates, Inc., 2020.