# Unconventional topological mixed-state transition and critical phase induced
by self-dual coherent errors

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

A topological phase can undergo a phase transition driven by anyon
condensation. A potential obstruction to such a mechanism could arise if there
exists a symmetry between anyons that have non-trivial mutual statistics. Here
we consider toric code subjected to errors that tend to proliferate anyons
with non-trivial mutual statistics. Using triangle inequality, we show that in
the presence of electromagnetic duality and a partial-transpose symmetry, a
decoherence induced phase transition out of the topological phase must be
rather unconventional and lie beyond standard rules of anyon condensation. To
explore such physics, we first subject toric code to a self-dual quantum
channel where Kraus operators are proportional to $X+Z$. We find that the
topological phase is stable up to the maximal error rate, when viewing density
matrix as a pure state in the double Hilbert space. To access an
unconventional transition, we then consider a perturbed toric code subjected
to the self-dual channel, and find numerical evidence that beyond a critical
error rate, the topological phase is destroyed resulting in a critical phase
where anyons are only power-law condensed.

##  I Introduction

Topological phases of matter host anyons that can have non-trivial mutual
statistics with respect to each other [1]. The anyonic statistics is not only
responsible for the universal properties of a gapped, topological phase, it
also puts strong constraints on the nature of quantum criticality and the
proximate phases. For example, condensing an anyon confines all other anyons
that have non-trivial mutual statistics with respect to it [2, 3, 4, 5, 6].
Relatedly, one cannot condense two anyons simultaneously that have non-trivial
mutual statistics. This can lead to phases of matter or critical points that
are difficult to understand in terms of anyon condensation [7, 8, 9, 10, 11,
12]. Motivated by such considerations, in this paper we study mixed states
obtained by subjecting topologically ordered pure states to quantum channels
that have a tendency to condense anyons with non-trivial mutual statistics.
For a 2+1-D $\mathbb{Z}_{2}$ topologically ordered state, we derive certain
general constraints that imply the decoherence driven phase transition that
destroys the topological order must violate the standard rules of anyon
condensation in one way or another. Specifically for the toric code, we find
numerical evidence for a phase transition out of the topological state to a
critical phase of matter in which anyons are ‘algebraically condensed’ i.e.
the associated string operators decay as a power law.

Figure 1:  Schematic phase diagram of the mixed state studied in this work
based on the 2nd Renyi entropy. The mixed state is obtained by subjecting the
pure state $|\Psi(h)\rangle=\prod_{e}[I+h(X_{e}+Z_{e})/\sqrt{2}]|\text{Toric
code}\rangle$ to the coherent channel
$\mathcal{E}_{e}[\cdot]=(1-p)(\cdot)+p\sigma_{e}(\theta)(\cdot)\sigma_{e}(\theta)$,
where $\sigma_{e}(\theta)=\cos(\theta)Z_{e}+\sin(\theta)X_{e}$. We focus
primarily on the self-dual line (i.e. $\theta=\pi/4$). The entanglement along
the red line is consistent with a central charge $c=2$ conformal field theory.
The plot shown here corresponds to $h=0.2$.

To motivate our setup, let us recall the phase diagram of the toric code
Hamiltonian [13], $H_{TC}$, perturbed by uniform fields along $X$ and $Z$
directions in the spin-space:
$H=H_{\text{TC}}-h_{x}\sum_{i}Z_{i}-h_{z}\sum_{i}X_{i}$ [14, 15, 7, 8, 9, 10,
12]. When $h_{x}>h_{z}$ ($h_{x}<h_{z}$), increasing $h_{x}$ ($h_{z}$) leads to
a second-order phase transition that can be understood as a condensation of a
self-boson anyon, namely, either the electric charge $e$ or the magnetic flux
$m$ [15]. However, along the self-dual line $h_{x}=h_{z}$, the mass gap to the
$e$ and the $m$ anyon is equal, and due to their $\pi$ mutual statistics, the
fate of the ground state is not so obvious as one leaves the topological
phase. Numerically, one finds an intriguing multicritical point along the
$h_{x}=h_{z}$ line as one exits the topological phase [7, 8, 9, 10, 12] whose
universal exponents are approximately known, but the corresponding field
theory is not well-understood [12].

Let us now consider a different setup where the toric code ground state is
subjected to local decoherence, so that the system is described by a mixed
state. For example, let’s first subject the ground state to a channel with
local Kraus operators $\propto X$ whose action on the ground state generates a
pair of electric anyons $e$ (in the convention where the plaquette term in the
toric code Hamiltonian is $\prod_{e\in p}X_{e}$ while the star term is
$\prod_{e\in v}Z_{e}$). Increasing the strength of the decoherence leads to a
phase transition to a non-topological phase which can be characterized in
multiple ways [16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26]. One way to
understand this transition is via the ‘double-state formalism’ [27, 18, 20],
in which the density matrix is mapped to a pure state $|\rho\rangle$ that
lives in the tensor product of the original Hilbert space with its own copy
[28, 29, 30, 31]. One can show that quantities linear in density matrix cannot
detect local decoherence induced phase transitions, and one of the simplest
quantities that does capture such transitions is the normalization of the
double state $\langle\rho|\rho\rangle$, which equals the second Renyi entropy
of the mixed state [18, 19, 20]. Furthermore, if the double state
$|\rho\rangle$ is short-ranged-entangled (SRE), then the density matrix
$\propto\rho^{2}\otimes I$ can be expressed as a convex sum of SRE states
[32]. Although the universality class or the threshold for the transition
within the double state formalism is not identical to that for the intrinsic
mixed-state transition for a single copy of the system [16, 18, 19, 21, 22,
25], it still provides an intuitive understanding of the transitionIn
particular, the phase transition for the double state corresponds to
condensing the bound state $e\bar{e}$ where $e$ and $\bar{e}$ refer to the
(electric) charge anyons in the two copies respectively [18, 20]. Similarly,
subjecting the toric code ground state to local Kraus operators $\propto Z$
also destroys the double topologically ordered state and the corresponding
transition can now be understood as condensing the bound state $m\bar{m}$ of
anyons where $m$ and $\bar{m}$ refer to the (magnetic) flux anyons in the two
copies respectively [18, 20]. If one now simultaneously subjects the system to
the composition of the aforementioned two channels, the resulting transition
can also be understood in terms of anyon condensation, since $e\bar{e}$ and
$m\bar{m}$ are mutual bosons. What happens if one instead considers a quantum
channel with Kraus operators of the form $X+Z$? Now it is not obvious that one
can condense either $e\bar{e}$ or $m\bar{m}$ anyons, since the probability
amplitude for condensing $e$ or $m$ in either copy is equal due to the
symmetry of the channel. This is reminiscent to the aforementioned (pure)
toric code ground state along the self-dual line. Such Kraus operators will be
the focus of our study.

To illustrate the special properties of the Kraus operator $X+Z$, we first ask
a general question: what conditions must the Kraus operators satisfy such that
a conventional phase transition out of the topologically ordered double state
of the toric code is not even allowed? Here by ‘conventional’ we mean a
transition that corresponds to condensation of a self-boson (such as
$e\bar{e}$) and where the standard rules of anyon condensation apply, e.g.,
one expects that the string operator corresponding to the condensing anyon has
long-range order, while anyons that have non-trivial mutual statistics with
the condensed anyon confine. In Sec.II we show that if a channel possesses
electromagnetic duality (EMD) symmetry and a ‘partial-transpose symmetry’
(which essentially corresponds to sending $X\rightarrow X,Z\rightarrow
Z,Y\rightarrow-Y$ in only one copy of the double state), then any transition
out of the topological phase must be unconventional. We show this result using
Cauchy-Schwarz (i.e. triangle) inequality. The basic idea is to express anyon
condensation rules in terms of asymptotic behavior of certain matrix elements
between excited states that correspond to the action of anyon-creating string
operators on the ground state. Using EMD and partial-transpose symmetry,
matrix elements corresponding to physically distinct processes can be related
to each other. At the same time, one can use triangle inequality to bound the
product of matrix elements in terms of other matrix elements. If one now
assumes that the phase transition out of the double topologically ordered
state happens due to condensation of a self-boson, one can then show that EMD
symmetry along with the partial-transpose symmetry leads to a contradiction
with the triangle inequality. On that note, the density matrix obtained by
decoherence with local Kraus operators $X+Z$ possesses both EMD and partial-
transpose symmetry, which motivate us to study this specific channel in
detail.

To illustrate the implications of the aforementioned constraints on
topological phase transitions due to EMD and partial-transpose symmetry, we
start with a class of pure states of the form
$|\Psi(h)\rangle=\prod_{e}(I+h(X_{e}+Z_{e})/\sqrt{2})|\textrm{Toric code
ground state}\rangle$ where $h\in[0,1]$. Such a pure state was first studied
in Ref.[11] where an interesting critical phase for a range of $h$ was
identified. We subject these class of pure states to Kraus operators of the
form $X\cos(\theta)+Z\sin(\theta)$. When $\theta=\pi/4$, the density matrix
respects both EMD and partial-transpose symmetry. We then use the
aforementioned double state formalism to map the problem to a non-unitarily
evolved pure state, where the evolution time depends on the strength of the
decoherence. The normalization of the double state can be interpreted as a
partition function of a statistical mechanics problem whose phase diagram
reveals different phases of the mixed state, analogously to the case of $Z$ or
$X$ decoherence [18, 20]. By employing this mapping, we analyze the phase
diagram for our problem with a combination of numerical techniques and
analytical arguments. The result is schematically shown in Fig.1. We find that
along the self-dual line, i.e., $\theta=\pi/4$, the double state undergoes a
phase transition from a topological ordered phase to a phase with critical
(i.e. power-law) correlations. The extent of the critical phase depends on the
parameter $h$ in the non-decohered original state $|\Psi(h)\rangle$. When
$h=0$, i.e., when the original state is simply the toric code, the correlation
functions along the self-dual line can be exactly mapped to the self-dual
Ashkin-Teller (AT) model, and it turns out that in this special case the
critical phase shrinks to a single point whose properties correspond to a
single, gapless complex scalar in $1+1$-D i.e., central charge $c=1$ theory.
In contrast, when $h>0$, one finds a critical phase for a whole range of
decoherence strength whose entanglement is consistent with a central charge
$c=2$ theory. We also find a relation of the corresponding statistical
mechanics model with two copies of the AT model along the self-dual line which
suggests that the low-energy theory of this phase is simply two copies of a
gapless, complex scalar.

The remaining part of the paper is organized as follows. In Sec.II, we discuss
the conditions under which a channel guarantees that a mixed-state transition
(if it exists) can not be described by a standard anyon condensation scheme.
In Sec.III, we consider a self-dual coherent channel that satisfies the
conditions in Sec.II and study its effect on the toric code ground state. In
Sec.IV, we apply the self-dual channel to a non-fixed point toric code and
obtain the phase diagram schematically sketched in Fig.1. We conclude in Sec.V
with a summary and discussion of our results.

##  II Constraints on anyon condensation from self-duality and partial-
transpose symmetry

The standard anyon condensation scheme leads to various constraints for
possible gapped ground states proximate to a given topological phase [2, 3, 4,
5, 6]. For example, condensing a self-boson anyon confines all particles that
have non-trivial mutual statistics with it. In this section we will study
constraints from the EMD symmetry and a partial-transpose symmetry (introduced
below) on anyon condensation in a decohered toric code, which we will serve as
a guide to search for decoherence induced novel transitions that must lie
beyond the standard anyon condensation scheme.

Let us first introduce our notation for anyon-creating string operators that
we will employ extensively. The Hamiltonian of the toric code can be written
as

$H_{\text{TC}}=-\sum_{v}A_{v}-\sum_{p}B_{p},$ (1)

where $A_{v}=\prod_{e\in v}Z_{e}$, $B_{p}=\prod_{e\in p}X_{e}$, and the qubits
reside on the edges (denoted as $e$) of a square lattice (Fig.2). The ground
state $|\Psi_{0}\rangle$ satisfies
$A_{v}|\Psi_{0}\rangle=B_{p}|\Psi_{0}\rangle=|\Psi_{0}\rangle,\ \forall v,p$.
A pair of charge anyons $|e_{i}e_{j}\rangle$ located at site $i$ and $j$
correspond to the violation of $A_{v}$ on sites $i$ and $j$, and can be
created by applying $w_{e}(l)=\prod_{e\in l}X_{e}$, a Pauli-$X$ string $l$
with the end points located at sites $i$ and $j$, to $|\Psi_{0}\rangle$ (see
Fig.2). Similarly, a pair of flux anyons $|m_{\tilde{i}}m_{\tilde{j}}\rangle$
located on the dual lattice with site $\tilde{i}$ and $\tilde{j}$ corresponds
to the violation of $B_{p}$ on sites $\tilde{i}$ and $\tilde{j}$, and can be
created by applying $w_{m}(\tilde{l})=\prod_{e\in\tilde{l}}Z_{e}$, a Pauli-$Z$
string $\tilde{l}$ with the end points located at sites $\tilde{i}$ and
$\tilde{j}$, to $|\Psi_{0}\rangle$.

Figure 2:  (a) Hamiltonian of the toric code and the string operator
$w_{e}(l)[w_{m}(\tilde{l})]$ that creates a pair of charge (flux) anyons at
sites $i(\tilde{i})$ and $j(\tilde{j})$. The orange arrows depict the
translation along the diagonal direction. (b) The resulting Hamiltonian and
the string operators after conjugating with the EMD unitary
$U_{D}=T_{\delta}\prod_{e}[(X_{e}+Z_{e})/\sqrt{2}]$.

The toric code Hamiltonian possesses an electro-magnetic duality (EMD)
symmetry corresponding to the unitary transformation

$U_{D}=T_{\delta}\prod_{e}\Big{(}\frac{X_{e}+Z_{e}}{\sqrt{2}}\Big{)}=T_{\delta}\prod_{e}\sigma_{e},$
(2)

where $T_{\delta}$ is a translation operator along the diagonal direction [see
Fig.2(a)]. The involvement of the translation operator in $U_{D}$ leads to
some subtleties which we now discuss. First, conjugating $U_{D}$ with the
aforementioned charge and flux creation operators yields

$U_{D}w_{e}(l)U_{D}^{\dagger}=w_{m}(\tilde{l}),\
U_{D}w_{m}(\tilde{l})U_{D}^{\dagger}=w_{e}({l}+2\delta)$ (3)

[see Fig.2(b)]. Therefore, EMD transformation exchanges the charge and flux
operators up to a lattice translation. One may then naively think that the
fermion creation operator $w_{f}(l)\equiv w_{e}(l)w_{m}(\tilde{l})$ is
invariant under the EMD transformation up to a lattice translation (we define
fermions on the original lattice without loss of generality). However, a
straightforward algebra shows that

$\displaystyle U_{D}w_{f}(l)U_{D}^{\dagger}$
$\displaystyle=(U_{D}w_{e}(l)U_{D}^{\dagger})(U_{D}w_{m}(\tilde{l})U_{D}^{\dagger})$
(4) $\displaystyle=w_{e}(l+2\delta)w_{m}(\tilde{l}),$

which is not equal to
$w_{f}(l+2\delta)=w_{e}(l+2\delta)w_{m}(\tilde{l}+2\delta)$. Indeed,
conjugating $w_{f}(l)$ with $U_{D}$ not only shifts the locations of the
fermions by $2\delta$ but also rotates them by an angle $\pi$. Second, since
we will be interested in applying a local quantum channel $\mathcal{E}$ to the
toric code, one can ask whether $\mathcal{E}$ respects the “strong” or “weak”
symmetry of $U_{D}$ [33]. Following the definition in Ref.[33], a channel
$\mathcal{E}$ respects the strong symmetry of $U_{D}$ if each of its Kraus
operator individually commutes with $U_{D}$ up to a constant phase factor
while only respects the weak symmetry if
$U^{\dagger}_{D}\mathcal{E}[U_{D}(\cdot)U_{D}^{\dagger}]U_{D}=\mathcal{E}[\cdot]$
. Since $U_{D}$ involves translation and we require $\mathcal{E}$ to be the
composition of local channels (we assume that the decoherence is generated by
local errors), $\mathcal{E}$ can only respect the weak symmetry of $U_{D}$.

Based on the intuition that one can’t simultaneously condense anyons with non-
trivial mutual statistics, one might naively think that a weak EMD symmetry is
already enough to rule out standard anyon condensation transition for a toric
code subjected to local decoherence. This is not true however. A channel that
satisfies the weak EMD symmetry is the composition of bit-flip and phase-flip
channels on all edges with equal error rate. It has already been shown in
Ref.[19] that for toric code, the bit-flip and phase-flip channels act
independently on the density matrix, and the statistical mechanics model for
the transition is simply two decoupled random-bond Ising models along the
Nishimori line. Therefore, the weak EMD symmetry is not enough to induce novel
mixed-state transitions. In fact, this can be easily understood by mapping the
density matrix into a state vector in the double Hilbert space [18, 20].
Denoting the charge (flux) anyons in the two copies as $e(m)$ and
$\bar{e}(\bar{m})$, the effect of the simultaneous bit-flip and phase-flip
errors then corresponds to condensing mutual bosons, $m\bar{m}$ and $e\bar{e}$
simultaneously [20, 18], which does not violate any anyon condensation rule.
Based on this intuition, it is reasonable to search for additional conditions
that may lead to stronger constraints that prohibit standard anyon
condensation transition in the double Hilbert space.

To keep the discussion self-contained, we will first briefly review the
double-state formalism introduced in Ref.[27, 18, 20], along with an overview
of order parameters to detect anyon condensation and confinement [34, 35, 36].
By employing this notation, we will soon argue that if one imposes a ‘partial-
transpose symmetry’ defined below in addition to the EMD symmetry, then the
double state corresponding to the decohered toric code is not allowed to
undergo a conventional anyon condensation transition.

Given a density matrix $\rho_{\mathcal{H}}$ acting on the Hilbert space
$\mathcal{H}$, one can define a state vector (which we refer to as “double
state”) in the doubled Hilbert space $\mathcal{H}\otimes\bar{\mathcal{H}}$
using Choi-Jamiolkowski (C-J) map [28, 29, 30, 31]:
$|\rho\rangle_{\mathcal{H}\otimes\bar{\mathcal{H}}}\equiv\rho_{\mathcal{H}}\otimes
I_{\bar{\mathcal{H}}}|\Phi\rangle_{\mathcal{H}\otimes\bar{\mathcal{H}}}$,
where $|\Phi\rangle_{\mathcal{H}\otimes\bar{\mathcal{H}}}$ is a product of
maximally entangled pairs connecting $\mathcal{H}$ and $\bar{\mathcal{H}}$.
Henceforth, for notational simplicity, we omit the subscripts labeling the
Hilbert space. After applying the C-J map, one finds that the initial pure
state density matrix $\rho_{0}=|\Psi_{0}\rangle\langle\Psi_{0}|$ is mapped to
$|\rho_{0}\rangle=|\Psi_{0}\rangle|\bar{\Psi}_{0}\rangle$ while the channel
$\mathcal{E}[\cdot]=\sum_{n}K_{n}(\cdot)K_{n}^{\dagger}$ is mapped to the
operator $\mathcal{E}=\sum_{n}K_{n}\otimes\bar{K}_{n}$. Therefore, the non-
decohered (initial) state $|\rho_{0}\rangle$ has double topological order
(e.g. the topological entanglement entropy will be $2\log(2)$), and the
decohered density matrix $\rho=\mathcal{E}[\rho]$ after the CJ map corresponds
to evolving $|\rho_{0}\rangle$ through the non-unitary operator $\mathcal{E}$,
i.e., $|\rho\rangle=\mathcal{E}|\rho_{0}\rangle$. It is also natural to
consider the error-corrupted anyon states
$|\rho_{\alpha\bar{\alpha}}\rangle=\mathcal{E}|\rho_{\alpha\bar{\alpha},0}\rangle=\mathcal{E}w_{\alpha}(l)\bar{w}_{\alpha}(l)|\rho_{0}\rangle$,
where $w_{\alpha}(l)$ is a string operator along the path $l$ that creates a
pair of anyons at the end points of $l$. To detect topological transitions in
the double Hilbert space, one can consider the following $\alpha\bar{\alpha}$
anyon condensation parameter [20]:

$\langle
I\bar{I}|\alpha\bar{\alpha}\rangle\equiv\lim_{|i-j|\rightarrow\infty}\frac{\langle\rho|\rho_{\alpha\bar{\alpha}}\rangle}{\langle\rho|\rho\rangle}=\lim_{|i-j|\rightarrow\infty}\frac{\langle\rho_{0}|\mathcal{E}^{\dagger}\mathcal{E}|\rho_{\alpha\bar{\alpha},0}\rangle}{\langle\rho_{0}|\mathcal{E}^{\dagger}\mathcal{E}|\rho_{0}\rangle},$
(5)

where $i$ and $j$ denote the end points of the string $w_{\alpha}(l)$.
Physically, $\langle I\bar{I}|\alpha\bar{\alpha}\rangle$ detects whether the
states $|\rho_{0}\rangle$ and $|\rho_{\alpha\bar{\alpha},0}\rangle$ are still
physically distinct after evolving both of them through the non-unitary
operator $\mathcal{E}$. Therefore, $\langle
I\bar{I}|\alpha\bar{\alpha}\rangle$ saturates to a finite constant when
$\alpha\bar{\alpha}$ is condensed. This will lead to the confinement of any
anyons $\beta\bar{\gamma}$ that have non-trivial mutual statistics with
$\alpha\bar{\alpha}$ [5, 2], which can be detected through the following anyon
confinement parameter [34, 35, 36]:

$\langle\beta\bar{\gamma}|\beta\bar{\gamma}\rangle\equiv\lim_{|i-j|\rightarrow\infty}\frac{\langle\rho_{\beta\bar{\gamma}}|\rho_{\beta\bar{\gamma}}\rangle}{\langle\rho|\rho\rangle},$
(6)

where
$|\rho_{\beta\bar{\gamma}}\rangle=\mathcal{E}|\rho_{\beta\bar{\gamma},0}\rangle=\mathcal{E}w_{\beta}(l)\bar{w}_{\gamma}(l)|\rho_{0}\rangle=\mathcal{E}|\Psi_{\beta},\bar{\Psi}_{\gamma}\rangle$.
Intuitively, $\langle\rho_{\beta\bar{\gamma}}|\rho_{\beta\bar{\gamma}}\rangle$
for a long string $l$ detects the probability for processes that lead to a
well-separated pair of the bound state $\beta\bar{\gamma}$. Confinement of
$\beta\bar{\gamma}$ implies that this amplitude decays exponentially with
$|l|=|i-j|$, and thus $\langle\beta\bar{\gamma}|\beta\bar{\gamma}\rangle=0$
when the anyon $\beta\bar{\gamma}$ is confined. From Eq.(5) and Eq.(6), one
observes that the most general quantities of interest are
$\langle{\alpha\bar{\beta}}|{\gamma\bar{\delta}}\rangle\equiv\langle\rho_{\alpha\bar{\beta}}|\rho_{\gamma\bar{\delta}}\rangle/\langle\rho|\rho\rangle$,
the “overlap” between the error-corrputed anyon eigenstates. In the original
Hilbert space, these quantities are all related to the second moment of the
density matrix. We also note that by definition, the order parameters
$\langle{\alpha\bar{\beta}}|{\gamma\bar{\delta}}\rangle$ are not conventional
physical observables of $|\rho\rangle$, i.e., they cannot be written as
$\langle\rho|O|\rho\rangle/\langle\rho|\rho\rangle$ where $O$ is a state-
independent operator. However, one may formulate
$\langle{\alpha\bar{\beta}}|{\gamma\bar{\delta}}\rangle$ as the expectation
value of a state-dependent operator
$O^{\alpha\bar{\beta}}_{\gamma\bar{\delta}}\equiv(\mathcal{E}^{\dagger})^{-1}w_{\alpha}^{\dagger}w_{\bar{\beta}}^{\dagger}\mathcal{E}^{\dagger}\mathcal{E}w_{\gamma}w_{\bar{\delta}}\mathcal{E}^{-1}+h.c.$
as long as the inverse of $\mathcal{E}$ exists. This implies that these order
parameters may not be a faithful measure of the underlying phase of matter
when $\mathcal{E}$ has non-trivial kernel. We will later encounter such a
situation in Sec.III.

We now propose the conditions under which a decoherence induced transition in
toric code cannot proceed via condensing a self-boson, and which therefore
guarantees a mixed-state transition (if any) beyond the standard anyon
condensation scheme. Consider the decohered double state
$|\rho\rangle=\mathcal{E}|\rho_{0}\rangle=\mathcal{E}|\Psi_{0},\bar{\Psi}_{0}\rangle$
with $|\Psi_{0}\rangle$ a $\mathbb{Z}_{2}$ topologically ordered state
respecting EMD and time-reversal (i.e., complex conjugation) symmetry. If the
C-J transformed channel $\mathcal{E}$ satisfies the following properties, then
any transition of $|\rho\rangle$ out of the double $\mathbb{Z}_{2}$
topologically ordered state can’t be described via standard anyon condensation
scheme.

  1. (i)

$\mathcal{E}$ respects the EMD symmetry:
${\mathbb{U}_{D}}\mathcal{E}{\mathbb{U}_{D}^{\dagger}}=\mathcal{E}$ where
$\mathbb{U}_{D}\equiv U_{\text{D}}\bar{U}_{\text{D}}$.

  2. (ii)

$\mathcal{E}^{\dagger}\mathcal{E}$ is invariant under partial transpose on
either $\mathcal{H}$ and $\bar{\mathcal{H}}$, i.e.,
$({\mathcal{E}}^{\dagger}\mathcal{E})^{T_{\mathcal{H}}}=({\mathcal{E}}^{\dagger}\mathcal{E})^{\bar{T}_{\bar{\mathcal{H}}}}={\mathcal{E}}^{\dagger}\mathcal{E}$,

where we denote $T_{\mathcal{H}}$ and $\bar{T}_{\bar{\mathcal{H}}}$ as the
partial transpose in the Hilbert space $\mathcal{H}$ and $\bar{\mathcal{H}}$,
respectively. We note that the condition (i) is nothing but the C-J
transformed version of the aforementioned weak EMD symmetry for $\mathcal{E}$.
To see the implication of condition (ii), let us introduce the notation
$|\Psi_{\alpha}\rangle\equiv w_{\alpha}(l)|\Psi_{0}\rangle$ for the $\alpha$
anyon state. Since $|\Psi_{\alpha}\rangle$ and $|\bar{\Psi}_{\beta}\rangle$
are time-reversal symmetric, condition (ii) implies

$\displaystyle\langle\Psi_{\alpha},\bar{\Psi}_{\beta}|\mathcal{E}^{\dagger}\mathcal{E}|\Psi_{\gamma},\bar{\Psi}_{\delta}\rangle$
$\displaystyle=\langle\Psi_{\gamma},\bar{\Psi}_{\beta}|\mathcal{E}^{\dagger}\mathcal{E}|\Psi_{\alpha},\bar{\Psi}_{\delta}\rangle$
(7)
$\displaystyle=\langle\Psi_{\alpha},\bar{\Psi}_{\delta}|\mathcal{E}^{\dagger}\mathcal{E}|\Psi_{\gamma},\bar{\Psi}_{\beta}\rangle.$

This means that the anyon overlap parameter
$\langle{\alpha\bar{\beta}}|{\gamma\bar{\delta}}\rangle$ is invaraint under
the exchange of the two copies in $\mathcal{H}/\bar{\mathcal{H}}$, i.e.,
$\langle{\alpha\bar{\beta}}|{\gamma\bar{\delta}}\rangle=\langle{\gamma\bar{\beta}}|{\alpha\bar{\delta}}\rangle=\langle\alpha\bar{\delta}|{\gamma\bar{\beta}}\rangle$.
To see that the EMD and partial-transpose symmetry together are inconsistent
with any anyon condensed phase, we will first use the EMD symmetry to rule out
all possibilities except for the following two: (a) both $e\bar{e}$ and
$m\bar{m}$ are condensed. (b) both $e\bar{m}$ and $m\bar{e}$ are condensed. We
will later show that these two situations are impossible if one further
insists on the partial-transpose symmetry.

Since $\mathbb{U}_{D}=U_{D}\bar{U}_{D}$ exchanges charges and fluxes on
$\mathcal{H}$ and $\bar{\mathcal{H}}$ simultaneously up to a lattice
translation and $|\rho_{0}\rangle$ is the eigenstate of $\mathbb{U}_{D}$, any
anyon overlap parameters
$\langle{\alpha\bar{\beta}}|{\gamma\bar{\delta}}\rangle$ without involving
fermions are invariant under the exchange of charges and fluxes on
$\mathcal{H}$ and $\bar{\mathcal{H}}$ simultaneously. Taking $\langle
I\bar{I}|m\bar{m}\rangle\propto\langle\rho_{0}|\mathcal{E}^{\dagger}\mathcal{E}w_{m}(\tilde{l})\bar{w}_{m}(\tilde{l})|\rho_{0}\rangle$
as an example, one finds

$\displaystyle\langle$
$\displaystyle\rho_{0}|\mathcal{E}^{\dagger}\mathcal{E}w_{m}(\tilde{l})\bar{w}_{m}(\tilde{l})|\rho_{0}\rangle$
(8)
$\displaystyle=\langle\rho_{0}|\mathbb{U}_{D}^{\dagger}\mathbb{U}_{D}\mathcal{E}^{\dagger}\mathcal{E}\mathbb{U}_{D}^{\dagger}\mathbb{U}_{D}w_{m}(\tilde{l})\bar{w}_{m}(\tilde{l})\mathbb{U}_{D}^{\dagger}\mathbb{U}_{D}|\rho_{0}\rangle$
$\displaystyle=\langle\rho_{0}|(\mathbb{U}_{D}\mathcal{E}^{\dagger}\mathcal{E}\mathbb{U}_{D}^{\dagger})(\mathbb{U}_{D}w_{m}(\tilde{l})\bar{w}_{m}(\tilde{l})\mathbb{U}_{D}^{\dagger})|\rho_{0}\rangle$
$\displaystyle=\langle\rho_{0}|\mathcal{E}^{\dagger}\mathcal{E}w_{e}({l}+2\delta)\bar{w}_{e}({l}+2\delta)|\rho_{0}\rangle$
$\displaystyle=\langle\rho_{0}|\mathcal{E}^{\dagger}\mathcal{E}w_{e}({l})\bar{w}_{e}({l})|\rho_{0}\rangle,$

which implies $\langle I\bar{I}|m\bar{m}\rangle=\langle
I\bar{I}|e\bar{e}\rangle$. In the final line, we have used the fact that both
$\mathcal{E}$ and $|\rho_{0}\rangle$ are translationally invariant. It is then
straightforward to see that the weak EMD symmetry rules out the possibility of
condensing self-bosons merely in $\mathcal{H}$ and $\bar{H}$. For example,
condensing $e$ implies $m$ is also condensed (since $\langle
I\bar{I}|e\bar{I}\rangle=\langle I\bar{I}|m\bar{I}\rangle$ by EMD symmetry),
which is in contradiction with the standard anyon condensation rule that the
anyons ($m$) with non-trivial mutual statistics with the condensed self-boson
($e$) should be confined. Therefore, if one requires $|\rho\rangle$ to
encounter a transition through anyon condensation in the presence of weak EMD
symmetry, only bound anyons in both $\mathcal{H}$ and $\bar{\mathcal{H}}$ can
condense, which leaves us with three possibilities distinguished by their
condensate generators: (a) condensation of $e\bar{e}$ and $m\bar{m}$, (b)
condensation of $e\bar{m}$ and $m\bar{e}$, or (c) condensation only of
$f\bar{f}=em\bar{e}\bar{m}$. It turns out that the situation (c) can also be
ruled out by the EMD symmetry, as we now argue that condensing $f\bar{f}$ in
fact implies that both $e\bar{e}$ and $m\bar{m}$ must also condense. To see
this, let’s use Eq.(4) to find how the $f\bar{f}$ creation operators
$w_{f}(l)\bar{w}_{f}(l)=w_{e}(l)\bar{w}_{e}(l)w_{m}(\tilde{l})\bar{w}_{m}(\tilde{l})$
transform under $\mathbb{U}_{D}$:
$\mathbb{U}_{D}w_{f}(l)\bar{w}_{f}(l)\mathbb{U}_{D}^{\dagger}=[w_{e}(l+2\delta)w_{m}(\tilde{l})][\bar{w}_{e}(l+2\delta)\bar{w}_{m}(\tilde{l})]=w_{f,\pi}(l)\bar{w}_{f,\pi}$,
where we denote $w_{f,\pi}\equiv w_{e}(l+2\delta)w_{m}(\tilde{l})$ as the
fermion creation operators rotated by angle $\pi$. Therefore, the same
argument in Eq.(8) shows that the EMD symmetry requires
$w_{f,\pi}(l)\bar{w}_{f,\pi}(l)$ to be condensed whenever
$w_{f}(l)\bar{w}_{f}(l)$ is condensed. However, fusing
$w_{f}(l)\bar{w}_{f}(l)$ and $w_{f,\pi}(l)\bar{w}_{f,\pi}(l)$ together yields
$w_{e}(l)\bar{w}_{e}(l)w_{e}(l+2\delta)\bar{w}_{e}(l+2\delta)$, which
corresponds to creating four $e\bar{e}$ anyons at the endpoints of $l$ and
$l+2\delta$. It is then obvious that repeatedly fusing $w_{f}\bar{w}_{f}$ and
$w_{f,\pi}\bar{w}_{f,\pi}$ corresponds to creating four $e\bar{e}$ anyons far
apart. Therefore, condensing $f\bar{f}$ along with the EMD symmetry implies
that $e\bar{e}$ is also condensed. One can use the similar argument to show
that $m\bar{m}$ is condensed as well.

The remaining two possibilities that are compatible with the EMD symmetry are:
(a) condensing $e\bar{e}$ and $m\bar{m}$, and (b) condensing $e\bar{m}$ and
$m\bar{e}$. We now use triangle (Cauchy-Schwarz) inequality to show that that
both of these possibilities are also disallowed if one further imposes the
aforementioned partial-transpose symmetry. Let’s assume $|\rho\rangle$
encounters a transition through condensing $e\bar{e}$ and $m\bar{m}$. This
implies that $\langle I\bar{I}|e\bar{e}\rangle=\langle
I\bar{I}|m\bar{m}\rangle$ saturates to a finite constant. Since $e$ and
$\bar{e}$ anyons have a non-trivial mutual statistics with $m\bar{m}$, they
must be confined, which implies $\langle e\bar{I}|e\bar{I}\rangle=\langle
I\bar{e}|I\bar{e}\rangle=0$. On the other hand, the partial-transpose symmetry
requires $\langle I\bar{I}|e\bar{e}\rangle=\langle I\bar{e}|e\bar{I}\rangle$
[due to Eq.(7)] to be a finite constant. However, a simple use of Cauchy-
Schwarz inequality $\langle A|B\rangle\leq\sqrt{\langle A|A\rangle\langle
B|B\rangle}$ yields

$\displaystyle\big{(}\langle\rho_{0}|$
$\displaystyle\bar{w}_{e}\mathcal{E}^{\dagger}\big{)}\big{(}\mathcal{E}w_{e}|\rho_{0}\rangle\big{)}$
(9)
$\displaystyle\leq\sqrt{\langle\rho_{0}|\bar{w}_{e}\mathcal{E}^{\dagger}\mathcal{E}\bar{w}_{e}|\rho_{0}\rangle\langle\rho_{0}|w_{e}\mathcal{E}^{\dagger}\mathcal{E}w_{e}|\rho_{0}\rangle}.$

After dividing
$\langle\rho_{0}|\mathcal{E}\mathcal{E}^{\dagger}|\rho_{0}\rangle$ on both
sides, Eq.(9) implies $\langle I\bar{e}|e\bar{I}\rangle\leq\sqrt{\langle
e\bar{I}|e\bar{I}\rangle\langle I\bar{e}|I\bar{e}\rangle}=0$, which is in
contradiction with $\langle I\bar{e}|e\bar{I}\rangle$ being a finite constant.
Therefore, the assumption of $|\rho\rangle$ encountering a phase transiton
thorough condensing both $e\bar{e}$ and $m\bar{m}$ must be wrong. One can
similarly show that condensing $e\bar{m}$ and $m\bar{e}$ is disallowed:
partial-transpose symmetry requires $\langle e\bar{I}|I\bar{m}\rangle=\langle
I\bar{I}|e\bar{m}\rangle$ being a finite constant. However, condensing
$e\bar{m}$ and $m\bar{e}$ implies confining $e$ and $\bar{m}$, i.e., $\langle
e\bar{I}|e\bar{I}\rangle=\langle I\bar{m}|I\bar{m}\rangle=0$, which violates
the Cauchy-Schwarz inequality $\langle
e\bar{I}|I\bar{m}\rangle\leq\sqrt{\langle e\bar{I}|e\bar{I}\rangle\langle
I\bar{m}|I\bar{m}\rangle}$. As a result, we reach the conclustion that the C-J
transformed channel $\mathcal{E}$ respects both EMD and partial-transpose
symmetry, any transition out of the topologically ordered toric code phase
must lie beyond the standard ayon condensation scheme. In the next two
sections we consider channels that do indeed respect EMD and partial-transpose
symmetries, and consider the fate of the toric code ground state and its close
relative due to the action of such channels.

##  III Example 1: self-dual coherent errors in toric code

We now consider one simple example that satisfies the aforementioned
conditions for a phase transition that can’t be described within the standard
anyon condensation scheme. Specifically, we subject the toric code ground
state to the composition of the following local channel on all edges $e$:

$\mathcal{E}_{e}[\rho]=(1-p)\rho+p\sigma_{e}\rho\sigma_{e},\
\sigma_{e}=\frac{Z_{e}+X_{e}}{\sqrt{2}}.$ (10)

Since $\sigma_{e}$ creates coherent superposition between states with
different types of anyons, we will use the adjective ‘coherent’ for the
channel in Eq.(10). After C-J map,
$\mathcal{E}=\prod_{e}[(1-p)+p\sigma_{e}\bar{\sigma}_{e}]$ and one can easily
show that
${U_{\text{D}}\bar{U}_{\text{D}}\mathcal{E}(U_{\text{D}}\bar{U}_{\text{D}})}^{\dagger}=\mathcal{E}$.
Besides, using the property $\sigma^{T}_{e}=\sigma_{e}$ and
$\mathcal{E}_{p}^{\dagger}\mathcal{E}_{p}=\mathcal{E}_{p^{\prime}}$ with
$p^{\prime}=2p(1-p)$, one finds
$({\mathcal{E}}^{\dagger}\mathcal{E})^{T_{\mathcal{H}}}=({\mathcal{E}}^{\dagger}\mathcal{E})^{{T}_{\bar{\mathcal{H}}}}={\mathcal{E}}^{\dagger}\mathcal{E}$.
Therefore, the transition for $|\rho\rangle$ out of the double topologically
ordered state should lie beyond the anyon condensation scheme as discussed in
the previous section.

In the following, we will map the transition induced by the channel in Eq.(10)
to several observables of the isotropic Ashkin-Teller (AT) model along the
self-dual line, which explicitly shows that $|\rho\rangle$ encounters a
Berezinskii–Kosterlitz–Thouless (BKT) transition [37, 38]. Interestingly, this
transition happens at the maximum error rate (i.e. $p=1/2$), and this is
suggestive that the mixed state is topologically ordered up to $p=1/2$. Later,
in Sec.IV, we will apply the channel in Eq.IV to a non-fixed point toric code
which will allow us to access an unconventional phase transition out of the
topologically ordered phase into a mixed state with critical correlations that
is stable over a whole range of parameters, as long as the EMD symmetry is
preserved.

###  III.1 Statistical mechanical mapping

The method of mapping the double state of toric code under decoherence to a
statistical mechanical model follows closely to the idea in Ref.[11], where a
pure-state transition for a single toric code was studied. Since the overall
scheme is fairly general and independent of the channel, we will first discuss
the mapping without specifying the explicit form of the channel. We will later
focus on the channel in Eq.(10) and show that the corresponding classical
model is the AT model along the self-dual line.

To begin with, we make use of the duality [39, 14, 40] by writing toric code
as the wave function $|\Psi_{0}\rangle$ obtained from performing forced
measurement of the 2d cluster state [41, 42]. The Hilbert state of the cluster
state consists of qubits both on the edges and vertices of a square lattice.
The cluster Hamiltonian takes the form
$H_{\textrm{cluster}}=\sum_{v}h_{v}+\sum_{e}h_{e}$, where
$h_{v}=-X_{v}A_{v}=-X_{v}\prod_{e\in v}Z_{e}$ and
$h_{e}=-X_{e}Z_{v}Z_{v^{\prime}}$ with $v$ and $v^{\prime}$ denote the two
vertices connected to the edge $e$. It follows that the toric code wave
function can be written as

$\displaystyle|\Psi_{0}\rangle$ $\displaystyle\propto\langle
x_{\mathbf{v}}=1|\Psi_{\text{cluster}}\rangle$ (11)
$\displaystyle\propto\langle
x_{\mathbf{v}}=1|\prod_{e}(I+X_{e}Z_{v}Z_{v^{\prime}})|x_{\mathbf{v}}=1,z_{\mathbf{e}}=1\rangle$
$\displaystyle=\sum_{x_{\mathbf{e}}}\sum_{z_{v}}\prod_{e}(1+x_{e}z_{v}z_{v^{\prime}})|x_{\mathbf{e}}\rangle.$

In the second line, we use
$|\Psi_{\text{cluster}}\rangle=\prod_{e}(1+h_{e})|x_{\mathbf{v}}=1,z_{\mathbf{e}}=1\rangle$
with $|x_{\mathbf{v}}=1,z_{\mathbf{e}}=1\rangle$ the product state satisfying
$h_{v}|x_{\mathbf{v}}=1,z_{\mathbf{e}}=1\rangle=1,\forall e$. On the other
hand, the final line can be obtained by inserting the complete basis
$\\{|z_{\mathbf{v}},x_{\mathbf{e}}\rangle\\}$ between
$(I+X_{e}Z_{v}Z_{v^{\prime}})$ and
$|x_{\mathbf{v}}=1,z_{\mathbf{e}}=1\rangle$.

Let’s now investigate how Eq.(11) assists our understanding of the effect of
decoherence and the aforementioned anyon condensation and confinement
parameters. We first compute the normalization of
$|\rho_{0}\rangle=|\Psi_{0},\bar{\Psi}_{0}\rangle$. Denoting
$\langle\Psi_{0}|=\sum_{x_{\mathbf{e}}}\sum_{t_{\mathbf{v}}}\prod_{e}(1+x_{e}t_{v}t_{v^{\prime}})\langle
x_{\mathbf{e}}|$, one finds the wavefunction normalization
$\langle\Psi_{0}|\Psi_{0}\rangle=\sum_{z_{\mathbf{v}},t_{\mathbf{v}}}\prod_{e}(1+z_{v}z_{v^{\prime}}t_{v}t_{v^{\prime}})=\sum_{z_{\mathbf{v}},t_{\mathbf{v}}}e^{\infty\sum_{e}z_{v}z_{v^{\prime}}t_{v}t_{v^{\prime}}}$.
It follows that

$\langle\rho_{0}|\rho_{0}\rangle=\sum_{z_{\mathbf{v}},\bar{z}_{\mathbf{v}},t_{\mathbf{v}},\bar{t}_{\mathbf{v}}}e^{\infty\sum_{e}(z_{v}z_{v^{\prime}}t_{v}t_{v^{\prime}}+\bar{z}_{v}\bar{z}_{v^{\prime}}\bar{t}_{v}\bar{t}_{v^{\prime}})}.$
(12)

Therefore, in the absence of decoherence, the Ising matter fields
corresponding to the state $|\Psi_{0}\rangle$, i.e., $z_{\mathbf{v}}$ and
$t_{\mathbf{v}}$ are strongly coupled with each other, and similarly those
originating from the state $|\bar{\Psi}_{0}\rangle$, i.e.,
$\bar{z}_{\mathbf{v}}$) and $\bar{t}_{\mathbf{v}}$) are also strongly coupled.
On the other hand, when the decoherence is introduced, the normalization of
the decohered double state $|\rho\rangle=\mathcal{E}|\rho_{0}\rangle$ takes
the form

$\displaystyle\langle\rho|\rho\rangle$
$\displaystyle=\sum_{z_{\mathbf{v}},\bar{z}_{\mathbf{v}},t_{\mathbf{v}},\bar{t}_{\mathbf{v}}}e^{-H(z_{\mathbf{v}},\bar{z}_{\mathbf{v}},t_{\mathbf{v}},\bar{t}_{\mathbf{v}})},$
(13)

where

$\displaystyle
e^{-H(z_{\mathbf{v}},\bar{z}_{\mathbf{v}},t_{\mathbf{v}},\bar{t}_{\mathbf{v}})}=\sum_{\begin{subarray}{c}x_{\mathbf{e}},\bar{x}_{\mathbf{e}},\\\
x^{\prime}_{\mathbf{e}},\bar{x}^{\prime}_{\mathbf{e}}\end{subarray}}\Big{[}\langle
x^{\prime}_{\mathbf{e}},\bar{x}^{\prime}_{\mathbf{e}}|{\mathcal{E}}^{\dagger}{\mathcal{E}}|x_{\mathbf{e}},\bar{x}_{\mathbf{e}}\rangle\times$
(14)
$\displaystyle\prod_{e}(1+x_{e}z_{v}z_{v^{\prime}})(1+\bar{x}_{e}\bar{z}_{v}\bar{z}_{v^{\prime}})(1+x^{\prime}_{e}t_{v}t_{v^{\prime}})(1+\bar{x}^{\prime}_{e}\bar{t}_{v}\bar{t}_{v^{\prime}})\Big{]}.$

While we haven’t specified the explicit form of the operator $\mathcal{E}$,
there is a constraint that $\mathcal{E}$ is a finite-depth local operator
since $\mathcal{E}[\cdot]$ is a finite-depth composition of local channels.
This implies that $H$ should be a local 2d Hamiltonian, as the four Ising
matter fields only interact with one another through $\langle
x^{\prime}_{\mathbf{e}},\bar{x}^{\prime}_{\mathbf{e}}|{\mathcal{E}}^{\dagger}{\mathcal{E}}|x_{\mathbf{e}},\bar{x}_{\mathbf{e}}\rangle$
in Eq.(14). Another essential feature of
$H(z_{\mathbf{v}},\bar{z}_{\mathbf{v}},t_{\mathbf{v}},\bar{t}_{\mathbf{v}})$
is that it will be invariant under the individual global spin-flip operation
$s_{\mathbf{v}}\rightarrow-s_{\mathbf{v}}$ for any
$s_{\mathbf{v}}=z_{\mathbf{v}},\bar{z}_{\mathbf{v}},t_{\mathbf{v}},\bar{t}_{\mathbf{v}}$.
This can be seen by noting that each term in the summation of Eq.(14) is
invariant under the individual global spin-flip operation. Therefore, it is
natural to expect that different quantum phases of $|\rho\rangle$ are mapped
to different classical phases of $H$ depending on what symmetries are
spontaneously broken. For example, from Eq.(12), one finds that the initial
double toric code $|\Psi_{0},\bar{\Psi}_{0}\rangle$ is in the double partial
ordered phase characterized by $\langle
z_{i}t_{i}z_{j}t_{j}\rangle=\langle\bar{z}_{i}\bar{t}_{i}\bar{z}_{j}\bar{t}_{j}\rangle$
saturating to a non-zero constant as $|i-j|\rightarrow\infty$ while $\langle
z_{i}z_{j}\rangle=\langle
t_{i}t_{j}\rangle=\langle\bar{z}_{i}\bar{z}_{j}\rangle=\langle\bar{t}_{i}\bar{t}_{j}\rangle$
decay exponentially as a function of $|i-j|$. Remarkably, these order
parameters precisely corrrespond to the $e$ condensation and confinement
parameters no matter what the operator $\mathcal{E}$ is (see Appendix A.1 for
details):

$\displaystyle\langle I\bar{I}|e\bar{I}\rangle$
$\displaystyle=\lim_{|i-j|\rightarrow\infty}\langle z_{i}z_{j}\rangle,$ (15)
$\displaystyle\langle e\bar{I}|e\bar{I}\rangle$
$\displaystyle=\lim_{|i-j|\rightarrow\infty}\langle
z_{i}t_{i}z_{j}t_{j}\rangle.$

One can similarly show that the $e\bar{e}$ condensation parameter $\langle
I\bar{I}|e\bar{e}\rangle$ is mapped to $\lim_{|i-j|\rightarrow\infty}\langle
z_{i}\bar{z}_{i}z_{j}\bar{z}_{j}\rangle$ and also vanishes for the intial
double toric code. On the other hand, the $m\bar{m}$ condensation parameter is
mapped to the following disorder parameter:

$\displaystyle\langle I\bar{I}|m\bar{m}\rangle$
$\displaystyle=\lim_{|\tilde{i}-\tilde{j}|\rightarrow\infty}\frac{\sum_{z_{\mathbf{v}},\bar{z}_{\mathbf{v}},t_{\mathbf{v}},\bar{t}_{\mathbf{v}}}e^{-H_{z_{\tilde{l}}\bar{z}_{\tilde{l}}}(z_{\mathbf{v}},\bar{z}_{\mathbf{v}},t_{\mathbf{v}},\bar{t}_{\mathbf{v}})}}{\sum_{z_{\mathbf{v}},\bar{z}_{\mathbf{v}},t_{\mathbf{v}},\bar{t}_{\mathbf{v}}}e^{-H(z_{\mathbf{v}},\bar{z}_{\mathbf{v}},t_{\mathbf{v}},\bar{t}_{\mathbf{v}})}}$
(16)
$\displaystyle\equiv\lim_{|\tilde{i}-\tilde{j}|\rightarrow\infty}\langle\mu^{z}_{\tilde{i}}\mu^{\bar{z}}_{\tilde{i}}\mu^{z}_{\tilde{j}}\mu^{\bar{z}}_{\tilde{j}}\rangle.$

Here $\tilde{l}$ denotes any string on the dual lattice with the endpoints
located at sites $\tilde{i}$ and $\tilde{j}$, and

$\displaystyle
e^{-H_{z_{\tilde{l}}\bar{z}_{\tilde{l}}}(z_{\mathbf{v}},\bar{z}_{\mathbf{v}},t_{\mathbf{v}},\bar{t}_{\mathbf{v}})}=\sum_{\begin{subarray}{c}x_{\mathbf{e}},\bar{x}_{\mathbf{e}},\\\
x^{\prime}_{\mathbf{e}},\bar{x}^{\prime}_{\mathbf{e}}\end{subarray}}\Big{[}\langle
x^{\prime}_{\mathbf{e}},\bar{x}^{\prime}_{\mathbf{e}}|{\mathcal{E}}^{\dagger}{\mathcal{E}}|x_{\mathbf{e}},\bar{x}_{\mathbf{e}}\rangle\times$
(17) $\displaystyle\prod_{e}$
$\displaystyle(1+\eta_{e}x_{e}z_{v}z_{v^{\prime}})(1+\eta_{e}\bar{x}_{e}\bar{z}_{v}\bar{z}_{v^{\prime}})(1+x^{\prime}_{e}t_{v}t_{v^{\prime}})(1+\bar{x}^{\prime}_{e}\bar{t}_{v}\bar{t}_{v^{\prime}})\Big{]}$

with $\eta_{e}=-1(+1)$ if $e\in\tilde{l}(\ni\tilde{l})$. One similarly finds
that the $m$ confinement parameter $\langle
m\bar{I}|m\bar{I}\rangle=\lim_{|\tilde{i}-\tilde{j}|\rightarrow\infty}\langle\mu^{z}_{\tilde{i}}\mu^{t}_{\tilde{i}}\mu^{z}_{\tilde{j}}\mu^{t}_{\tilde{j}}\rangle$.
Here
$\langle\mu^{z}_{\tilde{i}}\mu^{t}_{\tilde{i}}\mu^{z}_{\tilde{j}}\mu^{t}_{\tilde{j}}\rangle$
is defined analogously to
$\langle\mu^{z}_{\tilde{i}}\mu^{\bar{z}}_{\tilde{i}}\mu^{z}_{\tilde{j}}\mu^{\bar{z}}_{\tilde{j}}\rangle$
(see Eq.(16)). In fact, we can identify anyon overlap parameters
$\langle{\alpha\bar{\beta}}|{\gamma\bar{\delta}}\rangle\equiv\langle\rho_{\alpha\bar{\beta}}|\rho_{\gamma\bar{\delta}}\rangle/\langle\rho|\rho\rangle$
as certain combinations of order and disorder parameters in the statistical
model (see Appendix A.1).

Another interesting feature of the mapping is that the anyon condensation
rules now correspond to the constraints on order and disorder parameters of
any $\mathbb{Z}_{2}$ symmetric system with finite correlation length [43]. For
example, one finds $\langle I\bar{I}|e\bar{e}\rangle\propto\langle
z_{i}\bar{z}_{i}z_{j}\bar{z}_{j}\rangle$ detects the spontaneous breaking of
any symmetries that $z_{i}\bar{z}_{j}$ is charged under while $\langle
m\bar{I}|m\bar{I}\rangle\propto\langle\mu^{z}_{\tilde{i}}\mu^{{t}}_{\tilde{i}}\mu^{z}_{\tilde{j}}\mu^{{t}}_{\tilde{j}}\rangle$
detects the disorder of $z_{\mathbf{v}}t_{\mathbf{v}}$. The anyon condensation
rule that the condensation of $e\bar{e}$ implies the confinement of $m$ then
corresponds to the constraint that the order and disorder parameters of the
spin-flip operation
$(z_{\mathbf{v}},t_{\mathbf{v}})\rightarrow-(z_{\mathbf{v}},t_{\mathbf{v}})$
cannot be nonzero at the same time for any finite correlation length system
[43], and thus
$\lim_{|\tilde{l}|\rightarrow\infty}\langle\mu^{z}_{\tilde{i}}\mu^{{t}}_{\tilde{i}}\mu^{z}_{\tilde{j}}\mu^{{t}}_{\tilde{j}}\rangle$
must vanish when $\lim_{|i-j|\rightarrow\infty}\langle
z_{i}\bar{z}_{i}z_{j}\bar{z}_{j}\rangle$ is a finite constant. This statement
requires an assumption that the lowest eigenvalue of the corresponding
transfer matrix within the even subspace is non-degenerate, see Ref.[44] for
an example where this assumption is violated. Given this connection, it’s
illuminating to investigate how the weak EMD and partial-transpose symmetry
imply $|\rho\rangle$ must induce transition (if any) beyond anyon condensation
scheme using the constraints on order and disorder parameters. Let’s assume
that $e\bar{e}$ and $m\bar{m}$ condense simultaneously, i.e.,
$\lim_{|i-j|\rightarrow\infty}\langle
z_{i}\bar{z}_{i}z_{j}\bar{z}_{j}\rangle=\lim_{|\tilde{i}-\tilde{j}|\rightarrow\infty}\langle\mu^{z}_{i}\mu^{\bar{z}}_{i}\mu^{z}_{j}\mu^{\bar{z}}_{j}\rangle$
is a finite constant (the situation of condensing $e\bar{m}$ and $m\bar{e}$
can be analyzed using similar arguments). On the other hand, the partial-
transpose symmetry means $\lim_{|i-j|\rightarrow\infty}\langle
z_{i}\bar{t}_{i}z_{j}\bar{t}_{j}\rangle=\lim_{|i-j|\rightarrow\infty}\langle
z_{i}\bar{z}_{i}z_{j}\bar{z}_{j}\rangle$ is also nonzero. However, since
$\langle\mu^{z}_{i}\mu^{\bar{z}}_{i}\mu^{z}_{j}\mu^{\bar{z}}_{j}\rangle$ is
the disorder parameter of the symmetry transformation
$(z_{\mathbf{v}},\bar{z}_{\mathbf{v}})\rightarrow-(z_{\mathbf{v}},\bar{z}_{\mathbf{v}})$
and $z_{i}\bar{t}_{i}$ is charged under this symmetry,
$\lim_{|i-j|\rightarrow\infty}\langle z_{i}\bar{t}_{i}z_{j}\bar{t}_{j}\rangle$
and
$\lim_{|\tilde{i}-\tilde{j}|\rightarrow\infty}\langle\mu^{z}_{i}\mu^{\bar{z}}_{i}\mu^{z}_{j}\mu^{\bar{z}}_{j}\rangle$
cannot be nonzero at the same time for any system with finite correlation
length. We therefore reach a contradiction with the initial assumption that
$e\bar{e}$ and $m\bar{m}$ are condensed at the same time.

Having discussed the general mapping, we now restrict ourselves to the self-
dual coherent channel in Eq.(10) and show the it is mapped to the AT model
along the self-dual line. We calculate the normalization of the double state
in Appendix A.2, and find that

$\displaystyle\langle\rho|\rho\rangle$
$\displaystyle=\sum_{z_{\mathbf{v}},\bar{z}_{\mathbf{v}},t_{\mathbf{v}}}e^{-H(z_{\mathbf{v}},\bar{z}_{\mathbf{v}},t_{\mathbf{v}})},$
(18) $\displaystyle-H$
$\displaystyle=K\sum_{e}(z_{v}z_{v^{\prime}}\bar{z}_{v}\bar{z}_{v^{\prime}}+t_{v}t_{v^{\prime}}\bar{z}_{v}\bar{z}_{v^{\prime}})+K_{4}\sum_{e}z_{v}z_{v^{\prime}}t_{v}t_{v^{\prime}},$

where $K$ and $K_{4}$ are determined by $p$ through

$\displaystyle\tanh(K)=\frac{2-\sqrt{4-{\lambda}^{2}}}{\lambda},\
\tanh(K_{4})=\frac{2-\lambda\sqrt{4-{\lambda}^{2}}}{2-{\lambda}^{2}},$ (19)

with $\lambda=2p(1-p)/(1-2p+2p^{2})$. In particular, $p=0$ corresponds to
$(K,K_{4})=(0,\infty)$ while $p=1/2$ corresponds to
$K=K_{4}=\tanh^{-1}(2-\sqrt{3})\approx 0.275$. Therefore, increasing
decoherence effectively decreases the intra-copy interaction, while increasing
the inter-copy interaction. Note that we have used the fact that
$z_{v}z_{v^{\prime}}\bar{z}_{v}\bar{z}_{v^{\prime}}t_{v}t_{v^{\prime}}\bar{t}_{v}\bar{t}_{v^{\prime}}=1,\forall
e$ to eliminate $\bar{t}_{\mathbf{v}}$ in Eq.(18) (see Appendix A.2 for
details). It is also worth mentioning that the coupling stengths of
$z_{v}z_{v^{\prime}}\bar{z}_{v}\bar{z}_{v^{\prime}}$ and
$t_{v}t_{v^{\prime}}\bar{z}_{v}\bar{z}_{v^{\prime}}$ being equal is precisely
due to the partial-transpose symmetry of $\mathcal{E}^{\dagger}\mathcal{E}$.
To show the equivalence between $\langle\rho|\rho\rangle$ and the partition
function of the isotropic Ashkin-Teller model, we define new variables
$s_{v}=z_{v}\bar{z}_{v},\ \tau_{v}=t_{v}\bar{z}_{v}$ and write Eq.(18) as

$-H=K\sum_{e}(s_{v}s_{v^{\prime}}+\tau_{v}\tau_{v^{\prime}})+K_{4}\sum_{e}s_{v}s_{v^{\prime}}\tau_{v}\tau_{v^{\prime}},$
(20)

which is the desired form of the isotropic AT model. Besides, Eq.(19) implies
$K$ and $K_{4}$ are related to each other through $e^{-2K_{4}}=\sinh(2K)$,
which is the self-duality condition of the AT model. This is not a
coincidence, as the self-duality condition essentially corresponds to the
average EMD duality. To make their connection more precise, we rewrite the
$e\bar{e}$ condensation parameter $\langle I\bar{I}|e\bar{e}\rangle$ and the
$e$ confinement parameter $\langle e\bar{I}|e\bar{I}\rangle$ in terms of the
new variables, and they take the form

$\displaystyle\langle I\bar{I}|e\bar{e}\rangle$
$\displaystyle=\lim_{|i-j|\rightarrow\infty}\langle
z_{i}\bar{z}_{i}z_{j}\bar{z}_{j}\rangle=\lim_{|i-j|\rightarrow\infty}\langle
s_{i}s_{j}\rangle,$ (21) $\displaystyle\langle e\bar{I}|e\bar{I}\rangle$
$\displaystyle=\lim_{|i-j|\rightarrow\infty}\langle
z_{i}t_{i}z_{j}t_{j}\rangle=\lim_{|i-j|\rightarrow\infty}\langle
s_{i}\tau_{i}s_{j}\tau_{j}\rangle.$

Similarly, one finds

$\displaystyle\langle I\bar{I}|m\bar{m}\rangle$
$\displaystyle=\lim_{|\tilde{i}-\tilde{j}|\rightarrow\infty}\langle\mu^{z}_{\tilde{i}}\mu^{\bar{z}}_{\tilde{i}}\mu^{z}_{\tilde{j}}\mu^{\bar{z}}_{\tilde{j}}\rangle$
(22)
$\displaystyle=\lim_{|\tilde{i}-\tilde{j}|\rightarrow\infty}\langle\mu^{\tau}_{\tilde{i}}\mu^{\tau}_{\tilde{j}}\rangle=\lim_{|\tilde{i}-\tilde{j}|\rightarrow\infty}\langle\mu^{s}_{\tilde{i}}\mu^{s}_{\tilde{j}}\rangle,$
$\displaystyle\langle m\bar{I}|m\bar{I}\rangle$
$\displaystyle=\lim_{|\tilde{i}-\tilde{j}|\rightarrow\infty}\langle\mu^{z}_{\tilde{i}}\mu^{{t}}_{\tilde{i}}\mu^{z}_{\tilde{j}}\mu^{{t}}_{\tilde{j}}\rangle$
$\displaystyle=\lim_{|\tilde{i}-\tilde{j}|\rightarrow\infty}\langle\mu^{s}_{\tilde{i}}\mu^{\tau}_{\tilde{i}}\mu^{s}_{\tilde{j}}\mu^{\tau}_{\tilde{j}}\rangle.$

It is then obvious that requiring $\langle I\bar{I}|e\bar{e}\rangle=\langle
I\bar{I}|m\bar{m}\rangle$ and $\langle e\bar{I}|e\bar{I}\rangle=\langle
m\bar{I}|m\bar{I}\rangle$ forces the classical model to lie along the self-
dual line.

The exact mapping allows us to borrow the knowledge of the AT model [45] to
understand the fate of the decohered mixed state [see Fig.3(a)]. Along the
self-dual line of the isotropic AT model and when $K_{4}>K(p<1/2)$, it belongs
to the partial ordered phase with $\langle
s_{j}\rangle=\langle\tau_{j}\rangle=0$ while $\langle s_{j}\tau_{j}\rangle\neq
0$. Therefore, $\langle s_{i}s_{j}\rangle$ decays exponentially as a function
of $|i-j|$ while $\lim_{|i-j|\rightarrow\infty}\langle
s_{i}\tau_{i}s_{j}\tau_{j}\rangle$ saturates to a finite constant. This
implies $e\bar{e}$ is not condensed and $e$ is deconfined, consistent with the
double topologically ordered state. On the other hand, when
$K=K_{4}=\tanh^{-1}(2-\sqrt{3})$(corresponding to $p=1/2$), the classical
model hits the Berezinskii-Kosterlitz-Thouless (BKT) transition point and then
enters into a critical line with continuously varying exponents for $K\geq
K_{4}$. However, for the decohered state, $K>K_{4}$ is not possible since
$p=1/2$ is maximum allowed error rate. In Sec.IV we will access the BKT regime
by starting with a perturbed toric code.

###  III.2 Convex decomposition at the maximum error rate

The above mapping of the double state normalization to the AT model shows that
the double topological order of $|\rho_{0}\rangle$ is stable as long as the
error rate is less than the maximum value of $p=1/2$. This is suggestive that
perhaps the decohered mixed state $\rho$ is long-range entangled for $p<1/2$.
To establish this, one would need to calculate a true measure of mixed-state
entanglement such as entanglement negativity, or perhaps more ambitiously,
show that the mixed state is not separable for $p<1/2$, i.e., there exists no
convex decomposition of the mixed state in terms of SRE states for $p<1/2$.
Here we would content ourselves with studying the separability aspects just at
$p=1/2$. When $p=1/2$,
$\mathcal{E}\propto\prod_{e}(I+\sigma_{e}\bar{\sigma}_{e})$ is a projector to
the even sector of $\sigma_{e}\bar{\sigma}_{e}$ and hence is non-invertible
(i.e. has a non-trivial kernel). Therefore, as mentioned in Sec.II, the power-
law decay of several anyon condensation parameters does not guarantee that
$|\rho\rangle$ has long-range correlations. Furthermore, even if the double
state $|\rho\rangle$ has long-range correlations, they may originate from
classical correlations of the original density matrix $\rho$ and do not
necessarily imply long-range quantum entanglement for $\rho$ itself.

We now explictly show that $\rho$ at $p=1/2$ can be written as a convex sum of
product states, which guarantees that $\rho$ can at most have long-range
classical correlations. This can be done by noting that the Kraus operators
for a given channel are non-unique, and the self-dual coherent channel in
Eq.(10) can be equivalently written as $\mathcal{E}_{e}[\rho]=K_{e,+}\rho
K_{e,+}+K_{e,-}\rho K_{e,-}$, where the Kraus operators
$K_{e,\pm}=(\sqrt{1-p}I\pm\sqrt{p}\sigma_{e})/\sqrt{2}$. Given such a
representation, one can express $\rho$ as the following convex decomposition:

$\displaystyle\rho$
$\displaystyle=\sum_{m_{\mathbf{e}}}\Big{(}\prod_{e}K_{e,m_{e}}|\Psi_{0}\rangle\Big{)}\Big{(}\langle\Psi_{0}|\prod_{e}K_{e,m_{e}}\Big{)}$
(23)
$\displaystyle=\sum_{m_{\mathbf{e}}}|\psi_{m_{\mathbf{e}}}\rangle\langle\psi_{m_{\mathbf{e}}}|=\sum_{m_{\mathbf{e}}}P_{m_{\mathbf{e}}}|\tilde{\psi}_{m_{\mathbf{e}}}\rangle\langle\tilde{\psi}_{m_{\mathbf{e}}}|.$

Here $m_{\mathbf{e}}=\\{m_{e}=\pm 1\\}$,
$|\psi_{m_{\mathbf{e}}}\rangle=\prod_{e}K_{e,m_{e}}|\Psi_{0}\rangle$ is the
non-normalized wave function,
$P_{m_{\mathbf{e}}}=\langle\psi_{m_{\mathbf{e}}}|\psi_{m_{\mathbf{e}}}\rangle$
is the probability for $|\psi_{m_{\mathbf{e}}}\rangle$ to occur, and
$|\tilde{\psi}_{m_{\mathbf{e}}}\rangle$ is the normalized version of
$|\psi_{m_{\mathbf{e}}}\rangle$. Since $K_{e,\pm}$ becomes a projector to the
$\pm$ sector of $\sigma_{e}$ at $p=1/2$, one finds
$\sigma_{e}|\tilde{\psi}_{m_{\mathbf{e}}}\rangle=m_{e}|\tilde{\psi}_{m_{\mathbf{e}}}\rangle$,
and thus $|\tilde{\psi}_{m_{\mathbf{e}}}\rangle$ is a product state in the
$\sigma_{e}$ basis.

It is also worth mentioning that the probability for each
$|\psi_{m_{\mathbf{e}}}\rangle$ to occur at a given error rate $p$ takes the
following interesting form:

$\displaystyle P_{m_{\mathbf{e}}}=$ (24)
$\displaystyle\prod_{e}\Big{[}1+\frac{m_{e}\mu}{\sqrt{2}+m_{e}}(z_{v}z_{v^{\prime}}+t_{v}t_{v^{\prime}})+\frac{\sqrt{2}-m_{e}\mu}{\sqrt{2}+m_{e}\mu}z_{v}z_{v^{\prime}}t_{v}t_{v^{\prime}}\Big{]},$

where $\mu=2\sqrt{p(1-p)}$. At $p=1/2$, one can use
$\langle\tilde{\psi}_{m_{\mathbf{e}}^{\prime}}|\tilde{\psi}_{m_{\mathbf{e}}}\rangle=\delta_{m_{\mathbf{e}}^{\prime},m_{\mathbf{e}}}$
to show that $\tr(\rho^{2})=\sum_{m_{\mathbf{e}}}P^{2}_{m_{\mathbf{e}}}$,
which is the classical 2nd Renyi entropy of the convex decomposition.
Therefore, the density matrix at $p=1/2$ possesses classical long-range
correlations characterized by the classical 2nd Renyi entropy.

Figure 3: (a) The phase diagram for the double state of toric code subjected
to the self-dual coherent channel with Kraus operator
$\sigma_{e}=(Z_{e}+X_{e})/\sqrt{2}$. The BKT transition happens at the maximum
error rate $p=1/2$, and before that $|\rho\rangle$ possesses double
topologically ordered. (b) The phase diagram of the toric code subjected to a
more general coherent channel with with Kraus operator
$\sigma_{e}=\cos(\theta)Z_{e}+\sin(\theta)X_{e}$.

###  III.3 More general coherent-error channel

Let us briefly comment on the fate of toric code under a more general coherent
channel with Kraus operator
$\sigma_{e}(\theta)=\cos(\theta)Z_{e}+\sin(\theta)X_{e}$ which generically
does not respect the EMD symmetry. We find that $\langle\rho|\rho\rangle$ is
now mapped to the general isotropic AT model, with $K$ and $K_{4}$ related to
$\lambda$ and $\theta$ through (see Appendix A.2 for details)

$\displaystyle\frac{\tanh(K)[1+\tanh(K_{4})]}{1+\tanh^{2}(K)\tanh(K_{4})}$
$\displaystyle=\frac{\lambda\sin^{2}(\theta)}{1+\lambda\cos^{2}(\theta)},$
(25)
$\displaystyle\frac{\tanh^{2}(K)+\tanh(K_{4})}{1+\tanh^{2}(K)\tanh(K_{4})}$
$\displaystyle=\frac{1-\lambda\cos^{2}(\theta)}{1+\lambda\cos^{2}(\theta)}.$

One can now use our knowledge about the phase diagram of the general isotropic
ferromagnetic AT model to understand the phase diagram of the double state
$|\rho\rangle$ subjected to the channel with Kraus operator
$\sigma_{e}(\theta)$. In particular, the ferromagnetic AT model has three
phases [45]:

(i) The partial ordered phase with $\langle
s_{j}\rangle=\langle\tau_{j}\rangle=0$ while $\langle s_{j}\tau_{j}\rangle\neq
0$. In this regime, $\lim_{|i-j|\rightarrow\infty}\langle s_{i}s_{j}\rangle=0$
and thus $\langle I\bar{I}|e\bar{e}\rangle=0$. On the other hand, the ordering
of $s_{\mathbf{v}}\tau_{\mathbf{v}}$ implies that the free energy cost of
creating the domain wall for $s_{\mathbf{v}}\tau_{\mathbf{v}}$ must diverge
with $|\tilde{l}|$. As a result,
$\lim_{|\tilde{l}|\rightarrow\infty}\langle\mu^{s}_{\tilde{i}}\mu^{s}_{\tilde{j}}\rangle=0$
and thus $\langle I\bar{I}|m\bar{m}\rangle=0$. One can similarly show that
both $\langle e\bar{I}|e\bar{I}\rangle$ and $\langle m\bar{I}|m\bar{I}\rangle$
are finite. Therefore, $|\rho\rangle$ is a double $\mathbb{Z}_{2}$
topologically ordered state, suggesting that the quantum memory may be
retained for the mixed state $\rho$.

(ii) The Baxter (ordered) phase with $\langle
s_{j}\rangle=\langle\tau_{j}\rangle\neq 0$ and $\langle
s_{j}\tau_{j}\rangle\neq 0$. In this regime,
$\lim_{|i-j|\rightarrow\infty}\langle s_{i}s_{j}\rangle$ is a finite constant
and thus $\langle I\bar{I}|e\bar{e}\rangle$ is finite. Similar to the partial
ordered phase, the ordering of $\langle s_{j}\sigma_{j}\rangle$ implies
$\langle I\bar{I}|m\bar{m}\rangle=0$. One can also show that both $\langle
e\bar{I}|e\bar{I}\rangle$ and $\langle m\bar{I}|m\bar{I}\rangle$ are finite
using similar arguments. Therefore, $|\rho\rangle$ only has a single
topological order due to the condensation of $e\bar{e}$. This suggests that
the mixed state $\rho$ only possesses classical memory.

(iii) The paramagnetic phase with $\langle
s_{j}\rangle=\langle\tau_{j}\rangle=0$ and $\langle s_{j}\tau_{j}\rangle=0$.
Due to the disordering of $\langle s_{j}\tau_{j}\rangle$, the free energy cost
of creating the domain wall for $s_{\mathbf{v}}\tau_{\mathbf{v}}$ along the
path $\tilde{l}$ now becomes finite independent of the system size. Therefore,
$\langle I\bar{I}|m\bar{m}\rangle$ saturates to a finite value. On the other
hand, $\langle I\bar{I}|e\bar{e}\rangle=0$ similar to the partial ordered
phase. As a result, $|\rho\rangle$ only has a single topological order due to
the condensation of $m\bar{m}$, suggesting that the mixed state $\rho$ only
possesses classical memory.

We summarize these observations schematically in the phase diagram shown in
Fig.3(b).

##  IV Example 2: self-dual coherent errors in a perturbed toric code

Figure 4:  (a) Phase diagram for the double state of the perturbed toric code
(Eq. (26)) subjected to the self-dual coherent channel. Here the locations of
BKT, $\text{BKT}^{2}$, and $\text{Ising}^{4}$ are known from exact mappings
while the other green dots are determined by calculations similar to (c). (b)
Finite-entanglement scaling of the approximate MPS at different locations. The
legend shows the numerically extracted central charges at various locations in
the $(h,p)$-plane. (c) Inverse of the correlation length $1/\xi$ as a function
of the error rate $p$ along the lines $h=0.1$, $0.2$, and $0.3$. The solid
vertical lines mark the positions where $\xi$ starts to diverge.

In the previous section we showed the double topologically ordered toric code
ground state is stable under the self-dual coherent channel [Eq.(10)] as long
as $p<1/2$. One may wonder whether this statement is true for any
$\mathbb{Z}_{2}$ topologically ordered ground state subjected to self-dual
coherent channel. To make progress on this question, we consider applying the
same channel as in the previous section to a perturbed toric code that has
non-zero correlation length, and numerically study the corresponding
statistical mechanics model. We find that the transition induced by the
quantum channel can indeed happen below the maximum error rate, and as the
topological order is destroyed, the system enters a critical mixed-state phase
protected by the EMD symmetry (at least based on the quantities related to
$\tr(\rho^{2})$).

We start with the following pure state corresponding to a toric code ground
state $|\Psi_{0}\rangle$ perturbed by a non-unitary self-dual operator [11]:

$|\Psi(h)\rangle=\mathcal{N}|\Psi_{0}\rangle=\prod_{e}(I+h\sigma_{e})|\Psi_{0}\rangle,$
(26)

where $h\in[0,1]$ and we again denote $\sigma_{e}=(X_{e}+Z_{e})/\sqrt{2}$. We
emphasize that since $\mathcal{N}$ is a product of single-site non-unitary
operators, it acts as a temporal defect when considering the normalization of
$|\Psi(h)\rangle$ in the path intergral representation. Therefore,
$\mathcal{N}$ can only induce a (2+0)-D transition instead of the generic
(2+1)-D transition encountered when tuning parameters in a Hamiltonian. In
particular, it has been shown in Ref.[11] that $|\Psi(h)\rangle$ remains
topologically ordered when $h<h_{c}=\sqrt{2}-1\approx 0.414$ while for
$h>h_{c}$, it becomes a novel gapless state with continuously varying
exponents. The correlation functions of the anyon string operators in this
critical regime can be related to the conformal field theory corresponding to
a free, massless, complex scalar. $h=1$ corresponds to two decoupled critical
Ising models, while $h=h_{c}=\sqrt{2}-1$ corresponds to the BKT point at which
an operator becomes relevant and opens a mass gap resulting in the
topologically ordered phase for $h<h_{c}$.

Next, we subject $|\Psi(h)\rangle$ to the self-dual coherent channel in
Eq.(10). The corresponding double state can be written as
$|\rho\rangle=\mathcal{E}(\mathcal{N}\otimes\bar{\mathcal{N}})|\Psi_{0},\bar{\Psi}_{0}\rangle$.
This amounts to evolving the double toric code
$|\Psi_{0},\bar{\Psi}_{0}\rangle$ through the non-unitary operator
$\tilde{\mathcal{E}}=\mathcal{E}(\mathcal{N}\otimes\bar{\mathcal{N}})$.
Similar to the statistical mechanical mapping in Sec.III.1, we can write
$\langle\rho|\rho\rangle=\langle\Psi_{0},\bar{\Psi}_{0}|\tilde{\mathcal{E}}^{\dagger}\tilde{\mathcal{E}}|\Psi_{0},\bar{\Psi}_{0}\rangle$
as the partition function of 4-flavored Ising model that takes the form in
Eq.(14). One can also define the corrupted anyon state
$|\rho_{\alpha\bar{\beta}}\rangle=\tilde{\mathcal{E}}|\rho_{\alpha\bar{\beta},0}\rangle=\tilde{\mathcal{E}}w_{\alpha}(l)\bar{w}_{\beta}(l)|\rho_{0}\rangle$
and evaluate the anyon overlap parameters
$\langle{\alpha\bar{\beta}}|{\gamma\bar{\delta}}\rangle\equiv\langle\rho_{\alpha\bar{\beta}}|\rho_{\gamma\bar{\delta}}\rangle/\langle\rho|\rho\rangle$.
They again correspond to the combinations of order and disorder parameters.

Let us first briefly discuss the analytic structure of the statistical
mechanics model corresponding to $\langle\rho|\rho\rangle$. Recall that the
Hamiltonian of the 2D AT model along the self-dual line depends only on one
parameter, which we denote as $h$. Let’s write the partition function for a
single copy of translationally invariant AT model as
$\mathcal{Z}_{\text{AT}}=\sum_{t,z}\prod_{e}\omega_{e}(t,z;h)$ where
explicitly, $\omega_{e}$ on edge $e$ of a square lattice is given by
$\omega_{e}(t,z;h)=\left[\sqrt{2}+h+h(z_{v}z_{v^{\prime}}+t_{v}t_{v^{\prime}})+(\sqrt{2}-h)z_{v}z_{v^{\prime}}t_{v}t_{v^{\prime}}\right]$.
We find (see Appendix B)

$\displaystyle\langle\rho|\rho\rangle$ $\displaystyle=$
$\displaystyle\sum_{t,z,\bar{t},\bar{z}}\prod_{e}\left[\omega_{e}(t,z;h^{\prime})\omega_{e}(\bar{t},\bar{z};h^{\prime})\right.$
(27) $\displaystyle\left.\hskip
25.6073pt+f(p)\,\omega_{e}(1/h^{\prime},t,z)\omega_{e}(\bar{t},\bar{z};1/h^{\prime})\right]$

where $h’=2h/(1+h^{2})$ and $f(p)=2p(1-p)h^{\prime 2}/(1-2p+2p^{2})$. The
edges $e$ belong to a square lattice. By expanding the product in Eq.(27), one
can then write $\langle\rho|\rho\rangle$ as a sum of $2^{N_{e}}$ number of
partition functions ($N_{e}$ = total number of edges on the underlying square
lattice):

$\langle\rho|\rho\rangle=\sum_{\\{S\\}}|f(p)|^{|\bar{S}|}\mathcal{Z}_{S}\bar{\mathcal{Z}}_{S}$
(28)

where $\\{S\\}$ denotes the set of all distinct collection of edges, while
$\mathcal{Z}_{S}$ and $\bar{\mathcal{Z}}_{S}$ are partition functions for the
variables $\\{z_{\mathbf{v}},t_{\mathbf{v}}\\}$ and
$\\{\bar{z}_{\mathbf{v}},\bar{t}_{\mathbf{v}}\\}$ respectively.
$\mathcal{Z}_{S}$ is explicitly given by

$\mathcal{Z}_{S}=\sum_{t_{\mathbf{v}},z_{\mathbf{v}}}\prod_{e\in
S}\omega_{e}(z_{\mathbf{v}},t_{\mathbf{v}};h^{\prime})\prod_{e\notin{S}}\omega_{e^{\prime}}(z_{\mathbf{v}},t_{\mathbf{v}};1/h^{\prime})$
(29)

and $\bar{\mathcal{Z}}_{S}$ has the identical form with
$\\{t_{\mathbf{v}},z_{\mathbf{v}}\\}\rightarrow\\{\bar{t}_{\mathbf{v}},\bar{z}_{\mathbf{v}}\\}$.
Therefore, each $\mathcal{Z}_{S}$ corresponds to the partition function of a
self-dual AT model with binary disorder (the translational invariance is
restored due to the sum over $S$). More precisely, one can define a dual set
of variables $\mu^{z},\mu^{t}$ on the dual lattice and the self-duality
implies that the partition function of the dual variables is identical to
those of the original variables $z,t$ [45, 46]. This suggests independent
Kramers-Wannier duality for the fields $t_{\mathbf{v}},z_{\mathbf{v}}$ and
$\bar{t}_{\mathbf{v}},\bar{z}_{\mathbf{v}}$ in the low energy theory, whose
potential implications we will consider later. Incidentally, from Eq.(27), one
notices that when $h=1$, $\langle\rho|\rho\rangle$ is completely independent
of $p$ up to a proportionality constant, and therefore, corresponds to four
decoupled, critical Ising models for all $p$. We do not know an analytically
tractable way to solve the partition function in Eqs.(27),(29). Therefore, we
will instead use standard matrix product state (MPS) algorithms to carve out
the global phase diagram in the $(h,p)$-plane for this problem.

Fig.4(a) shows the global phase diagram in the $(h,p)$ plane. Here BKT refers
to the complex scalar CFT at the BKT point, while $\text{BKT}^{2}$ refers to
two copies of the BKT theory i.e. two complex scalars at the BKT point.
$\text{Ising}^{4}$ refers to four copies of the critical Ising theory. The
position of the BKT point, i.e., $(h=0,p=0.5)$, is known from the exact
mapping discussed in Sec.III.1 while the positions of $\text{BKT}^{2}$
$(h=\sqrt{2}-1,p=0)$ and $\text{Ising}^{4}$ $(h=1,p=0)$ are exactly known from
Ref.[11]. We now discuss the details of how this phase diagram is obtained.

The basic idea behind the numerics is as follows. One can write
$\langle\rho|\rho\rangle=\tr(T^{L_{y}})$ with $T$ the $2^{4L_{x}}\times
2^{4L_{x}}$ transfer matrix for the system of size $L_{x}\times L_{y}$. The
dominant eigenvector of $T$, which we denote as $|\psi\rangle$, is equivalent
to the ground state of the corresponding $(1+1)$-D Hamiltonian. Assuming
translational invariance, one can directly consider the thermodynamic limit
$L_{x},L_{y}\rightarrow\infty$ and then approximate the ground state
$|\psi\rangle$ using a translationally invariant MPS ansatz
$|\psi_{\chi}\rangle$ with a finite bond dimension $\chi$. For a system with
finite correlation length, the approximations of $|\psi_{\chi}\rangle$ as
$|\psi\rangle$ converge rapidly as one increases $\chi$ [47]. On the other
hand, $|\psi_{\chi}\rangle$ can never faithfully represent the entanglement of
true critical ground states, as its entanglement entropy $S_{\text{vN}}$ of
the bipartition on the infinite chain is bounded by $\log\chi$. However, it
turns out that the error of the approximation can be used to extract universal
behaviors of the critical system [48]. In particular, the von-Neumann entropy
$S_{\text{vN}}$, the bond dimension $\chi$, and the coorrelation length $\xi$
of $|\psi_{\chi}\rangle$ are related to one another through
$S_{\text{vN}}\approx c\kappa\log(\chi)/6\approx c\log(\xi)/6$ with $c$ the
central charge and $\kappa=6/(\sqrt{12c}+c)$. One can then consider a class of
$|\psi_{\chi}\rangle$ with different $\chi$ to extract the central charge of
the critical system, and such a method is called finite-entanglement scaling
[48]. In the following, we will first fix $\chi=48$ to compute the correlation
length of system in the $(h,p)$ plane [see Fig.4(c)]. The locations where
$\xi$ starts to diverge signals the possibility of critical phase (note that
for any finite $\chi$, the correlation length of $|\psi_{\chi}\rangle$ is
always finite). We will later use finite-entanglement scaling to provide more
evidence of the existence of the critical phase [Fig.4(b)].

Fig.4(c) shows the inverse of correlation length $1/\xi$ as a function the
error rate $p$ along the lines $h=0.1$, $0.2$, and $0.3$. When $p=0$, we find
that $1/\xi$ is some nonzero constant far from zero, indicating the system has
finite correlation length. This is consistent with the exact mapping that the
system is gapped when $h<h_{c}=\sqrt{2}-1\approx 0.414$ [11]. As one increases
the error rate, $1/\xi$ for a given $h$ decreases and eventually becomes close
to zero at a particular point $p_{c}(h)$, indicating the occurrence of a phase
transition. Besides, $1/\xi$ remains extremely small as one further increases
the error rate, suggesting that the system remains critical above $p_{c}(h)$.
We also checked that for any $p$ at $h>h_{c}$, the system always has very
large correlation length, suggesting that the system remains critical as one
increases the error rate. Further, we confirmed that the gapped phase below
the phase boundary is indeed a double topologically ordered state
characterized by the uncondenesd and deconfined charge/flux anyon parameters.
As mentioned above, for Fig.4(c) we used $\chi=48$. We also verified that
similar results are obtained when one further increases the bond dimension.

To confirm the presence of a critical phase outside the double topologically
ordered phase and to identify its universal features, we employ finite-
entanglement scaling [48] to extract the central charge $c$. We first
benchmark the positions $(h,p)=(0.0,0.5)$ and $(0.5,0.0)$, whose central
charges are exactly known to be $c=1$ and $c=2$, respectively. As shown in
Fig.4(b), we find the numerically extracted central charges are $c\approx
0.87$ at $(0.0,0.5)$ and $c\approx 1.86$ at $(0.5,0.0)$, which are close to
the exact values. We then proceed to compute central charges at various points
in the putative critical phase. We find that as long as the system is not too
close to the BKT point $(h,p)=(0.0,0.5)$, the numerically extracted central
charge is always around $c=2$ [see Fig.4(b)]. Therefore, we expect the regime
outside the double topologically ordered phase is described by a central
charge $c=2$ theory except for the special point at $(h,p)=(0.0,0.5)$ which
has central charge $c=1$. Since the system is known to be described by two
copies of the AT model at $p=0$ for $\sqrt{2}-1<h<1$, which has a marginal
deformation (the boson radius/Luttinger parameter), c-theorem [49] strongly
indicates that the whole critical regime is described by two complex, free,
massless scalars.

We also numerically considered the stability of the gapless phase against
perturbations that break the partial-transpose symmetry while retaining the
EMD symmetry. Such perturbations can be introduced by adding either the
channel with Kraus operator a Pauli-$Y$ matrix or the combinations of phase-
flip (Pauli-$Z$) and bit-flip (Pauli-$X$) channels with equal error rate. We
applied these perturbations simultaneously with equal error rate of
$p^{\prime}=0.1$ to the critical state at the location $(h,p)=(0.5,0.3)$, and
find the numerically extracted central charge remains close to $2$. This
suggests that our critical phase is also stable under the perturbations that
respect the EMD symmetry but break the partial-transpose symmetry.

Figure 5: (a) The phase diagram for the perturbed toric code in Eq.(30) under
phase-flip channel based on the quantities related to
$\log(\tr\rho^{n})/(1-n)$ with $n=2$. The central charge along the red line is
$c=1$, while that along the blue lines is $c=1/2$. (b) Conjectured phase
diagram for the limit $n\rightarrow 1$. The location of the green circle along
the $h=0$ line is known from Refs. [16, 50], the location of the blue circle
along the $p=0$ line is known from Ref.[51], while the location of the other
blue circle is determined in Appendix C.

Given that we are able to explore new mixed-state phases and phase transitions
by starting with a finite correlation length pure state, it is interesting to
consider subjecting finite correlation length pure states to other quantum
channels. To explore this direction, we also studied a differently perturbed
toric code subjected to phase flip errors. In short, we consider the state

$|\Psi(h)\rangle=\prod_{e}(1+hZ_{e})|\Psi_{0}\rangle\propto
e^{K\sum_{e}Z_{e}}|\Psi_{0}\rangle$ (30)

subjected to the phase-flip channel,
$\mathcal{E}[\cdot]=\prod_{e}\mathcal{E}_{e}[\cdot]$, where
$\mathcal{E}_{e}[\cdot]=(1-p)\rho+pZ_{e}\rho Z_{e}$. Interestingly, this case
is analytically tractable and we provide the details in Appendix C. The
resulting phase diagram is schematically shown in Fig.5. One finds a rather
rich phase diagram for the double state, with three phases, namely, a phase
where both copies are separately topologically ordered (so that topological
entanglement entropy (TEE) is $2\log(2)$), a phase where $e\bar{e}$ is
condensed (so that TEE is $\log(2)$), and a phase without any topological
order.

Besides the unconventional phase transition out of the topologically ordered
phase, perhaps the most striking feature in our phase diagram is the presence
of a critical phase for a whole range of parameters (Fig.4). It is worth
examining this phase and related aspects from the perspective of an effective
field theory. Let us start with the pure-state limit in the absence of
decoherence (i.e. $p=0$). The normalization
$\langle\rho|\rho\rangle=\langle\Psi(h)|\Psi(h)\rangle\langle\bar{\Psi}(h)|\bar{\Psi}(h)\rangle$
can be written as a path integral whose Lagrangian corresponds to two
decoupled sine-Gordon fields [11]:

$\mathcal{L}_{\text{sG}}=\Big{[}\frac{1}{2\pi
K}((\partial_{\tau}\phi)^{2}+(\partial_{x}\phi)^{2})+\lambda\sin(4\phi)\Big{]}+\Big{[}\phi\leftrightarrow\bar{\phi}\Big{]},$
(31)

where $K$ is the Luttinger parameter that depends on $h$, and
$\phi(\bar{\phi})$ denotes the bosonic field in
$\mathcal{H}(\bar{\mathcal{H}})$. The connection between $\phi$ and the
classical Ising spins $z_{\mathbf{v}},t_{\mathbf{v}}$ that enter the partition
function $\langle\rho|\rho\rangle$ is through the Majorana spinors
$\gamma_{R}+i\eta_{R}\sim e^{i(\phi+\theta)}$, $\gamma_{L}+i\eta_{L}\sim
e^{-i(\phi-\theta)}$, where $\theta$ is the field dual to $\phi$. Here
$\gamma=(\gamma_{L},\gamma_{R})$ [$\eta=(\eta_{L},\eta_{R})$] correspond to
the Ising spins $z_{\mathbf{v}}$($t_{\mathbf{v}}$) in the standard way [52,
53]: the transfer matrix for the critical Ising model corresponds to a (1+1)-D
transverse field Ising model at criticality, which can be fermionized via the
Jordan-Wigner transformation. Similar relations hold for $\bar{\phi}$ and
$\bar{z}_{\mathbf{v}},\bar{t}_{\mathbf{v}}$. EMD symmetry corresponds to the
lattice translation for Majorana fermions, which in the continuum limit acts
as
$(\gamma_{L},\eta_{L},\bar{\gamma}_{L},\bar{\eta}_{L})\rightarrow-(\gamma_{L},\eta_{L},\bar{\gamma}_{L},\bar{\eta}_{L})$
while the right-moving Majorana fermions remain invariant. The gapless regime
$h\in[\sqrt{2}-1,1]$ corresponds to $K\in[1/2,1]$ such that the scaling
dimension of $\sin(4\phi)=4K<2$, and thus this term is irrelevant. In
particular, the free-fermion limit $K=1$ corresponds to $h=1$, as
$\mathcal{N}(h)=\prod_{e}(I+h\sigma_{e})$ becomes a projector that decouples
$\gamma(\bar{\gamma})$ and $\eta(\bar{\eta})$. On the other hand, the double
topologically ordered phase corresponds to $K<1/2$ (i.e. $0\leq h<\sqrt{2}-1$)
and $\lambda>0$ such that $\sin(2\phi)\sim
i(\gamma_{L}\eta_{R}-\eta_{L}\gamma_{R})$ and $\sim
i(\bar{\gamma}_{L}\bar{\eta}_{R}-\bar{\eta}_{L}\bar{\gamma}_{R})$ are pinned
to $\pm 1$. This can be seen by recalling that the double topological phase
corresponds to the double partial order phase in the statistical model (see
Sec.III), and the partial order phase can be detected by the operator that is
even under the simultaneous spin-flip action
$\gamma,\eta\rightarrow-\gamma,-\eta$ but odd under the single spin-flip
operation $\gamma\rightarrow-\gamma$ or $\eta\rightarrow-\eta$.

Next, let’s discuss the physics of the critical phase (Fig.4). As discussed
above, numerically we find that the central charge of the statistical
mechanics model for the second Renyi entropy
$\tr\rho^{2}=\langle\rho|\rho\rangle$ throughout the critical phase equals two
within our numerical accuracy. Our expectation is that the low-energy theory
in this phase is again described by Eq.(31) where the Luttinger parameter will
generically be a function of the decoherence rate $p$ and the parameter $h$
that defines the pure state (Eq.(26)). Here we encounter a puzzle whose
resolution we don’t fully understand. The decohered density matrix only
respects the EMD symmetry that acts on both copies simultaneously, i.e.,
$\mathbb{U}_{D}^{\dagger}\rho\mathbb{U}_{D}=\rho$ where
$\mathbb{U}_{D}=U_{D}\bar{U}_{D}$ (see Eq.(2) for the definition of $U_{D}$).
Therefore, one might expect that the following operator, which is invariant
under $\mathbb{U}_{D}$ but odd under $U_{D}$ and $\bar{U}_{D}$, is allowed:

$\Delta\mathcal{L}=\lambda^{\prime}(\gamma_{L}\gamma_{R}+\eta_{L}\eta_{R})(\bar{\gamma}_{L}\bar{\gamma}_{R}+\bar{\eta}_{L}\bar{\eta}_{R}).$
(32)

After bosonizing Eq.(32),
$\Delta\mathcal{L}=\lambda^{\prime}\cos(2\phi)\cos(2\bar{\phi})$, which has
scaling dimension $\text{dim}[\cos(2\phi)\cos(2\bar{\phi})]=2K$. When $p=0$
and $h\in[\sqrt{2}-1,1)$, one finds that $K\in[1/2,1)$ and thus
$\Delta\mathcal{L}$ is relevant. This will lead to spontaneous symmetry
breaking resulting in $\cos(2\phi)=\cos(2\bar{\phi})=\pm 1$. This is at
variance with our numerical observation that the critical phase is in fact
stable everywhere in the phase diagram outside the topological phase. First,
let’s briefly consider the consequences for the low-energy physics if
$\Delta\mathcal{L}$ were indeed present in the effective field theory. One
expects that $\Delta\mathcal{L}$ being relevant will pin the system to a state
that breaks the EMD symmetry spontaneously, and satisfies
$\cos(2\phi)=\cos(2\bar{\phi})=\pm 1$ where $+1$ and $-1$ correspond to the
condensation of flux and charge respectively. This implies that the true
ground state will be a superposition of states corresponding to flux and
charge condensates. Such a state will violate the standard intuition that
anyons with non-trivial mutual statistics can’t condense at the same time.
Relatedly, the ground state of the transfer matrix of the corresponding
statistical mechanics model will have simultaneous existence of order and
disorder parameters, which is rather unusual but still allowed if certain
conditions are satisfied [44, 43]. We did not find evidence for such a state
in our phase diagram (Fig.4). However, it is indeed possible that in the
thermodynamic limit, the ground state is gapped with a very large correlation
length so that it appears critical for the maximum bond dimension we used.

Given the stability of the critical phase observed in our numerics, we suspect
that $\Delta\mathcal{L}$ is in fact not allowed at low-energies. The absence
of $\Delta\mathcal{L}$ in the low energy theory is suggested by Eqs.(27),(29),
where we find $\langle\rho|\rho\rangle$ can be expressed as sum of partition
functions that, despite being non-translationally invariant, are self-dual
under separate Kramers-Wannier duality for the fields
$\\{t_{\mathbf{v}},z_{\mathbf{v}}\\}$ and
$\\{\bar{t}_{\mathbf{v}},\bar{z}_{\mathbf{v}}\\}$. This suggests that one may
write the continuum version of the partition function as

$\displaystyle\mathcal{Z}$ $\displaystyle=$ $\displaystyle\int D\phi\int
DK\,D\lambda\,P(K,\lambda)e^{-\int dxd\tau\mathcal{L}}$
$\displaystyle\mathcal{L}$ $\displaystyle=$ $\displaystyle\Big{[}\frac{1}{2\pi
K(x,\tau)}\left((\partial_{\tau}\phi(x,\tau))^{2}+(\partial_{x}\phi(x,\tau))^{2}+\Big{[}\phi\leftrightarrow\bar{\phi}\Big{]}\right)\Big{.}$
(33)
$\displaystyle+\Big{.}\lambda(x,\tau)\left(\sin(4\phi(x,\tau))+\Big{[}\phi\leftrightarrow\bar{\phi}\Big{]}\right)\Big{]}$

where $P(K,\lambda)$ is the probability distribution for the couplings $K$ and
$\lambda$ and is the continuum analog of the factor $|f(p)|^{|\bar{S}|}$ in
the lattice partition function (see Eq.(27)), and integrating over $K,\lambda$
will generate interactions between the fields $\phi$ and $\bar{\phi}$. The
partition function in Eq.(IV) clearly possesses a ‘doubled Kramers-Wannier’
symmetry corresponding to independently sending $\phi\rightarrow\phi+\pi/2$,
or $\bar{\phi}\rightarrow\bar{\phi}+\pi/2$. We leave further investigation of
this interesting direction to the future.

Assuming that $\Delta\mathcal{L}$ is not allowed in the effective action, one
can then understand the phase transition between the topologically ordered
phase and the critical phase in terms of the most relevant operator allowed by
symmetries, namely, $\sin(4\phi)$. Approaching this transition from the
critical phase, as decoherence decreases, one expects that the Luttinger
parameter $K$ will decrease and eventually approach $K_{c}=1/2$ so that
$\sin(4\phi)$ becomes relevant. This is because as discussed in Sec.III.1, the
effect of decoherence will not only introduce the coupling between the fields
in $\mathcal{H}$ and $\bar{\mathcal{H}}$ (corresponding to a new term that
couples $\phi$ and $\bar{\phi}$) but also decrease the coupling between the
two copies. This implies decreasing the error rate is expected to decrease the
Luttinger parameter $K$. This qualitatively explains the phase boundary
between the topological phase and the critical phase in our phase diagram
(Fig.4).

##  V Summary and discussion

In this work we obtained constraints on anyon condensation in mixed states of
toric code by imposing weak EMD and partial-transpose symmetries. Throughout,
we focussed on the double state corresponding to the mixed state. Our
conclusion is that in the presence of EMD and partial-transpose symmetry,
conventional anyon condensation mechanism, where string operators
corresponding to anyons take non-zero expectation value, is ruled out. We
derived this result using a combination of symmetry arguments and the triangle
inequality. To illustrate the implications of this constraint, we provided an
example where a phase transition out of the double topologically ordered state
is necessarily unconventional in the sense that it cannot occur via
conventional anyon condensation. In particular, we studied a finite-
correlation version of toric code (Eq.(26)), and subjected it to a channel
that preserves EMD and partial-transpose symmetry. Above a critical strength,
the topological phase is destroyed via an unconventional transition, and we
find numerical evidence for a critical phase where anyons are only power-law
condensed, i.e., the string operators that characterize anyon condensation
have a power-law scaling. Our numerical results and analytical arguments are
consistent with the conclusion that these correlations are described by a CFT
corresponding to two free, massless, complex scalars. The phase diagram is
schematically shown in Fig.1 with more detailed aspects shown in Fig.4. We
numerically find that the critical phase is stable against breaking the
partial-transpose symmetry and only requires the EMD symmetry. For the special
case corresponding to the zero correlation length toric code subjected to the
same channel, we find that the double topologically ordered state is stable up
to the maximal decoherence strength, and right at the maximal decoherence
strength, the topological ordered phase is destroyed and the system possesses
critical correlations corresponding to a single free, massless, complex scalar
(Figs.3,4). For this case, the statistical mechanics model corresponding to
the normalization of the double state is mapped to a single copy of the
Ashkin-Teller model along the self-dual line in the critical regime. For the
more general case of the aforementioned, decohered perturbed toric code, the
statistical mechanics model can be expressed in terms of two coupled AT models
along the self-dual line.

In this work we only focussed on the double state subjected to $X+Z$ Kraus
operators. It will be very interesting to study the intrinsic properties of
the mixed state corresponding to a single copy, such as aspects related to
error-correction, mixed-state entanglement and separability [16, 17, 18, 19,
20, 21, 25]. For example, in the case of toric code decohered by phase-flip
errors, the double state undergoes a transition in the translationally
invariant Ising universality, while the universal aspects of error correction
and entanglement for the single copy are related to the random bond Ising
model along the Nishimori line. For our problem, the universality for the
double state transition is a free, complex scalar, and it is interesting to
speculate if the single copy transition belongs to a random version of this
problem along the Nishimori line. We leave this problem for future
investigation.

Another interesting direction is to subject the ground state of
$H=H_{\text{TC}}-h\sum(X_{i}+Z_{i})$ to the channel with Kraus operators
$X_{i}+Z_{i}$. Would the corresponding double state exhibit a phase diagram
similar to that in Fig.4? It will also be interesting to consider self-dual
decoherence in models where the EMD symmetry is onsite and does not involve
translations, such as those studied in [54]. Such models will allow one to
enforce strong EMD symmetry and may have a richer phase diagram. Another
potential direction is to generalize the arguments in Sec.II to more general
topologically ordered states with EMD symmetry such as $S_{3}$ gauge theory
[55].

Note added: While this manuscript was being finalized, Ref.[56] appeared which
has overlap with our Sec.III where we study the fixed-point toric code
decohered by $X+Z$ Kraus operators and its connection with the self-dual
Ashkin-Teller model. Although our focus and observables studied are different
from Ref.[56], the qualitative features do seem to match where they overlap.

###### Acknowledgements.

We thank John McGreevy for helpful comments on the manuscript. TG is supported
by the National Science Foundation under Grant No. DMR-1752417.

## References

  * Wen [2004] X.-G. Wen, Quantum field theory of many-body systems: from the origin of sound to an origin of light and electrons, oup oxford, (2004).
  * Burnell [2018] F. J. Burnell, Anyon condensation and its applications, Annual Review of Condensed Matter Physics 9, 307 (2018).
  * Bais  _et al._ [2002] F. A. Bais, B. J. Schroers, and J. K. Slingerland, Broken quantum symmetry and confinement phases in planar physics, Physical review letters 89, 181601 (2002).
  * Bais and Mathy [2007] F. Bais and C. Mathy, The breaking of quantum double symmetries by defect condensation, Annals of Physics 322, 552 (2007).
  * Bais and Slingerland [2009] F. A. Bais and J. K. Slingerland, Condensate-induced transitions between topologically ordered phases, Phys. Rev. B 79, 045316 (2009).
  * Neupert  _et al._ [2016] T. Neupert, H. He, C. Von Keyserlingk, G. Sierra, and B. A. Bernevig, Boson condensation in topologically ordered quantum liquids, Physical Review B 93, 115103 (2016).
  * Jongeward  _et al._ [1980] G. A. Jongeward, J. D. Stack, and C. Jayaprakash, Monte carlo calculations on z 2 gauge-higgs theories, Physical Review D 21, 3360 (1980).
  * Tupitsyn  _et al._ [2010] I. S. Tupitsyn, A. Kitaev, N. V. Prokof’ev, and P. C. E. Stamp, Topological multicritical point in the phase diagram of the toric code model and three-dimensional lattice gauge higgs model, Phys. Rev. B 82, 085114 (2010).
  * Vidal  _et al._ [2009] J. Vidal, S. Dusuel, and K. P. Schmidt, Low-energy effective theory of the toric code model in a parallel magnetic field, Phys. Rev. B 79, 033109 (2009).
  * Wu  _et al._ [2012] F. Wu, Y. Deng, and N. Prokof’ev, Phase diagram of the toric code model in a parallel magnetic field, Physical Review B 85, 195104 (2012).
  * Zhu and Zhang [2019] G.-Y. Zhu and G.-M. Zhang, Gapless coulomb state emerging from a self-dual topological tensor-network state, Phys. Rev. Lett. 122, 176401 (2019).
  * Somoza  _et al._ [2021] A. M. Somoza, P. Serna, and A. Nahum, Self-dual criticality in three-dimensional z 2 gauge theory with matter, Physical Review X 11, 041008 (2021).
  * Kitaev [2003] A. Y. Kitaev, Fault-tolerant quantum computation by anyons, Annals of Physics 303, 2 (2003).
  * Wegner [1971] F. J. Wegner, Duality in generalized ising models and phase transitions without local order parameters, Journal of Mathematical Physics 12, 2259 (1971).
  * Fradkin and Shenker [1979] E. Fradkin and S. H. Shenker, Phase diagrams of lattice gauge theories with higgs fields, Physical Review D 19, 3682 (1979).
  * Dennis  _et al._ [2002] E. Dennis, A. Kitaev, A. Landahl, and J. Preskill, Topological quantum memory, Journal of Mathematical Physics 43, 4452 (2002).
  * Wang  _et al._ [2003] C. Wang, J. Harrington, and J. Preskill, Confinement-higgs transition in a disordered gauge theory and the accuracy threshold for quantum memory, Annals of Physics 303, 31 (2003).
  * Lee  _et al._ [2023] J. Y. Lee, C.-M. Jian, and C. Xu, Quantum criticality under decoherence or weak measurement, PRX Quantum 4, 030317 (2023).
  * Fan  _et al._ [2023] R. Fan, Y. Bao, E. Altman, and A. Vishwanath, Diagnostics of mixed-state topological order and breakdown of quantum memory, arXiv preprint arXiv:2301.05689 (2023).
  * Bao  _et al._ [2023] Y. Bao, R. Fan, A. Vishwanath, and E. Altman, Mixed-state topological order and the errorfield double formulation of decoherence-induced transitions, arXiv preprint arXiv:2301.05687 (2023).
  * Chen and Grover [2023a] Y.-H. Chen and T. Grover, Separability transitions in topological states induced by local decoherence, arXiv preprint arXiv:2309.11879 (2023a).
  * Sang  _et al._ [2023] S. Sang, Y. Zou, and T. H. Hsieh, Mixed-state quantum phases: Renormalization and quantum error correction, arXiv preprint arXiv:2310.08639 (2023).
  * Li and Mong [2024] Z. Li and R. S. Mong, Replica topological order in quantum mixed states and quantum error correction, arXiv preprint arXiv:2402.09516 (2024).
  * Su  _et al._ [2024] K. Su, Z. Yang, and C.-M. Jian, Tapestry of dualities in decohered quantum error correction codes, arXiv preprint arXiv:2401.17359 (2024).
  * Lee [2024] J. Y. Lee, Exact calculations of coherent information for toric codes under decoherence: Identifying the fundamental error threshold, arXiv preprint arXiv:2402.16937 (2024).
  * Lyons [2024] A. Lyons, Understanding stabilizer codes under local decoherence through a general statistical mechanics mapping, (2024), arXiv:2403.03955 [quant-ph] .
  * Lee  _et al._ [2022] J. Y. Lee, Y.-Z. You, and C. Xu, Symmetry protected topological phases under decoherence, arXiv preprint arXiv:2210.16323 (2022).
  * Schmutz [1978] M. Schmutz, Real-time green’s functions in many body problems, Zeitschrift für Physik B Condensed Matter 30, 97 (1978).
  * Prosen [2008] T. Prosen, Third quantization: a general method to solve master equations for quadratic open fermi systems, New Journal of Physics 10, 043026 (2008).
  * Jamiołkowski [1972] A. Jamiołkowski, Linear transformations which preserve trace and positive semidefiniteness of operators, Reports on Mathematical Physics 3, 275 (1972).
  * Choi [1975] M.-D. Choi, Completely positive linear maps on complex matrices, Linear algebra and its applications 10, 285 (1975).
  * Chen and Grover [2023b] Y.-H. Chen and T. Grover, Symmetry-enforced many-body separability transitions, arXiv preprint arXiv:2310.07286 (2023b).
  * de Groot  _et al._ [2022] C. de Groot, A. Turz illo, and N. Schuch, Symmetry Protected Topological Order in Open Quantum Systems, Quantum 6, 856 (2022).
  * Schuch  _et al._ [2013] N. Schuch, D. Poilblanc, J. I. Cirac, and D. Pérez-Garc\́text{i}a, Topological order in the projected entangled-pair states formalism: Transfer operator and boundary hamiltonians, Phys. Rev. Lett. 111, 090501 (2013).
  * Haegeman  _et al._ [2015] J. Haegeman, V. Zauner, N. Schuch, and F. Verstraete, Shadows of anyons and the entanglement structure of topological phases, Nature communications 6, 8284 (2015).
  * Duivenvoorden  _et al._ [2017] K. Duivenvoorden, M. Iqbal, J. Haegeman, F. Verstraete, and N. Schuch, Entanglement phases as holographic duals of anyon condensates, Phys. Rev. B 95, 235119 (2017).
  * Berezinskii [1971] V. Berezinskii, Destruction of long-range order in one-dimensional and two-dimensional systems having a continuous symmetry group i. classical systems, Sov. Phys. JETP 32, 493 (1971).
  * Kosterlitz and Thouless [1973] J. M. Kosterlitz and D. J. Thouless, Ordering, metastability and phase transitions in two-dimensional systems, Journal of Physics C: Solid State Physics 6, 1181 (1973).
  * Kramers and Wannier [1941] H. A. Kramers and G. H. Wannier, Statistics of the two-dimensional ferromagnet. part i, Phys. Rev. 60, 252 (1941).
  * Kubica and Yoshida [2018] A. Kubica and B. Yoshida, Ungauging quantum error-correcting codes, arXiv preprint arXiv:1805.01836 (2018).
  * Raussendorf  _et al._ [2005] R. Raussendorf, S. Bravyi, and J. Harrington, Long-range quantum entanglement in noisy cluster states, Phys. Rev. A 71, 062313 (2005).
  * Aguado  _et al._ [2008] M. Aguado, G. K. Brennen, F. Verstraete, and J. I. Cirac, Creation, manipulation, and detection of abelian and non-abelian anyons in optical lattices, Phys. Rev. Lett. 101, 260501 (2008).
  * Levin [2020] M. Levin, Constraints on order and disorder parameters in quantum spin chains, Communications in Mathematical Physics 378, 1081 (2020).
  * O’Brien and Fendley [2018] E. O’Brien and P. Fendley, Lattice supersymmetry and order-disorder coexistence in the tricritical ising model, Phys. Rev. Lett. 120, 206403 (2018).
  * Baxter [2016] R. J. Baxter,  _Exactly solved models in statistical mechanics_ (Elsevier, 2016).
  * Wiseman and Domany [1995] S. Wiseman and E. Domany, Critical behavior of the random-bond ashkin-teller model: A monte carlo study, Physical Review E 51, 3074 (1995).
  * Eisert  _et al._ [2010] J. Eisert, M. Cramer, and M. B. Plenio, Colloquium: Area laws for the entanglement entropy, Rev. Mod. Phys. 82, 277 (2010).
  * Pollmann  _et al._ [2009] F. Pollmann, S. Mukerjee, A. M. Turner, and J. E. Moore, Theory of finite-entanglement scaling at one-dimensional quantum critical points, Phys. Rev. Lett. 102, 255701 (2009).
  * Zamolodchikov [1986] A. B. Zamolodchikov, Irreversibility of the flux of the renormalization group in a 2d field theory, JETP lett 43, 730 (1986).
  * Honecker  _et al._ [2001] A. Honecker, M. Picco, and P. Pujol, Universality class of the nishimori point in the 2d $\pm{}\mathit{J}$ random-bond ising model, Phys. Rev. Lett. 87, 047201 (2001).
  * Castelnovo and Chamon [2008] C. Castelnovo and C. Chamon, Quantum topological phase transition at the microscopic level, Physical Review B 77, 054433 (2008).
  * Schultz  _et al._ [1964] T. D. Schultz, D. C. Mattis, and E. H. Lieb, Two-dimensional ising model as a soluble problem of many fermions, Reviews of Modern Physics 36, 856 (1964).
  * Jordan and Wigner [1993] P. Jordan and E. P. Wigner,  _Über das paulische äquivalenzverbot_ (Springer, 1993).
  * Cheng  _et al._ [2017] M. Cheng, Z.-C. Gu, S. Jiang, and Y. Qi, Exactly solvable models for symmetry-enriched topological phases, Physical Review B 96, 115107 (2017).
  * Beigi  _et al._ [2011] S. Beigi, P. W. Shor, and D. Whalen, The quantum double model with boundary: condensations and symmetries, Communications in mathematical physics 306, 663 (2011).
  * Eckstein  _et al._ [2024] F. Eckstein, B. Han, S. Trebst, and G.-Y. Zhu, Robust teleportation of a surface code and cascade of topological quantum phase transitions, (2024), arXiv:2403.04767 [quant-ph] .
  * Grest and Widom [1981] G. S. Grest and M. Widom, $n$-color ashkin-teller model, Phys. Rev. B 24, 6508 (1981).

##  Appendix A Statistical mechanical mapping for toric code under local
decoherence

The goal of this Appendix is to discuss the statistical mechanical mapping for
the toric code under local decoherence. We will show in Appendix A.1 that the
anyon overlap parameters discussed in Sec.II can always be mapped to some
combinations of ordered and disordered parameters in the corresponding
classical model. In Appendix A.2, we will focus on the mapping for the
coherent channel discussed in the main text. In particular, we will map the
second moment of the density matrix (equivalent to the normalization of the
double state) to the partition function of the isotropic Ashkin-Teller model.

###  A.1 General mapping for anyon overlap parameters

We now show that the anyon overlap parameters
$\langle{\alpha\bar{\beta}}|{\gamma\bar{\delta}}\rangle\equiv\langle\rho_{\alpha\bar{\beta}}|\rho_{\gamma\bar{\delta}}\rangle/\langle\rho|\rho\rangle$
can be mapped to various combinations of order and disorder parameters of the
spin variables $(z,\bar{z})$ and $(t,\bar{t})$. To begin with, we first use
Eq.(11) to write down the initial double state:

$\displaystyle|\rho_{0}\rangle=|\Psi_{0}\rangle|\bar{\Psi}_{0}\rangle=\sum_{x_{\mathbf{e}},\bar{x}_{\mathbf{e}}}[\sum_{z_{v},\bar{z}_{v}}\prod_{e}(1+x_{e}z_{v}z_{v^{\prime}})(1+\bar{x}_{e}\bar{z}_{v}\bar{z}_{v^{\prime}})]|x_{\mathbf{e}},\bar{x}_{\mathbf{e}}\rangle.$
(34)

One can also easily write down the initial anyon state
$|\rho_{\alpha\beta,0}\rangle=w_{\alpha}\bar{w}_{\beta}|\rho_{0}\rangle$. Take
$|\rho_{e\bar{e},0}\rangle$ as an example. Since a pair of $e\bar{e}$ anyon
state located at sites $i$ and $j$ can be created by any string operator
$w_{e}(l)\bar{w}_{e}(l)=\prod_{e\in l}X_{e}\bar{X}_{e}$ with the end points
located at sites $i$ and $j$, one finds

$\displaystyle|\rho_{e\bar{e},0}\rangle$ $\displaystyle=\prod_{e\in
l}X_{e}\bar{X}_{e}|\rho_{0}\rangle=\sum_{x_{\mathbf{e}},\bar{x}_{\mathbf{e}}}\sum_{z_{\mathbf{v}},\bar{z}_{\mathbf{v}}}\prod_{e\in
l}(x_{e}+z_{v}z_{v^{\prime}})(\bar{x}_{e}+\bar{z}_{v}\bar{z}_{v^{\prime}})\prod_{e\notin
l}(1+x_{e}z_{v}z_{v^{\prime}})(1+\bar{x}_{e}\bar{z}_{v}\bar{z}_{v^{\prime}})|x_{e},\bar{x}_{e}\rangle$
(35)
$\displaystyle=\sum_{x_{\mathbf{e}},\bar{x}_{\mathbf{e}}e}[\sum_{z_{\mathbf{v}},\bar{z}_{\mathbf{v}}}\prod_{e}z_{i}z_{j}\bar{z}_{i}\bar{z}_{j}(1+x_{e}z_{v}z_{v^{\prime}})(1+\bar{x}_{e}\bar{z}_{v}\bar{z}_{v^{\prime}})]|x_{e},\bar{x}_{e}\rangle$
$\displaystyle=\sum_{z_{\mathbf{v}},\bar{z}_{\mathbf{v}}}z_{i}z_{j}\bar{z}_{i}\bar{z}_{j}\sum_{x_{\mathbf{e}},\bar{x}_{\mathbf{e}}}[\prod_{e}(1+x_{e}z_{v}z_{v^{\prime}})(1+\bar{x}_{e}\bar{z}_{v}\bar{z}_{v^{\prime}})]|x_{e},\bar{x}_{e}\rangle.$

In the second line, we use
$(x_{e}+z_{v}z_{v^{\prime}})(\bar{x}_{e}+\bar{z}_{v}\bar{z}_{v^{\prime}})=z_{v}z_{v^{\prime}}\bar{z}_{v}\bar{z}_{v^{\prime}}(1+x_{e}z_{v}z_{v^{\prime}})(1+\bar{x}_{e}\bar{z}_{v}\bar{z}_{v^{\prime}})$
and the property that all the spins $z_{v}$ with $v\notin i,j$ appear even
number of times. In the third line, we simply change the order of the
summations $\sum_{x_{\mathbf{e}},\bar{x}_{\mathbf{e}}}$ and
$\sum_{z_{\mathbf{v}},\bar{z}_{\mathbf{v}}}$ so that it becomes transparent
that the effect of the operator $\prod_{e\in l}X_{e}\bar{X}_{e}$ on
$|\rho_{0}\rangle$ is transformed to $z_{i}z_{j}\bar{z}_{i}\bar{z}_{j}$, the
two-point function of the matter field. Now, let’s consider the $e\bar{e}$
condensation parameters
${\langle\rho_{0}|\mathcal{E}^{\dagger}\mathcal{E}|\rho_{e\bar{e},0}\rangle}/{\langle\rho_{0}|\mathcal{E}^{\dagger}\mathcal{E}|\rho_{0}\rangle}$.
As already shown in the main text,
$\langle\rho_{0}|\mathcal{E}^{\dagger}\mathcal{E}|\rho_{0}\rangle=\sum_{z_{\mathbf{v}},\bar{z}_{\mathbf{v}},t_{\mathbf{v}},\bar{t}_{\mathbf{v}}}e^{-H(z_{\mathbf{v}},\bar{z}_{\mathbf{v}},t_{\mathbf{v}},\bar{t}_{\mathbf{v}})}$,
where

$\displaystyle
e^{-H(z_{\mathbf{v}},\bar{z}_{\mathbf{v}},t_{\mathbf{v}},\bar{t}_{\mathbf{v}})}$
$\displaystyle=\sum_{\begin{subarray}{c}x_{\mathbf{e}},\bar{x}_{\mathbf{e}},\\\
x^{\prime}_{\mathbf{e}},\bar{x}^{\prime}_{\mathbf{e}}\end{subarray}}\Big{[}\langle
x^{\prime}_{\mathbf{e}},\bar{x}^{\prime}_{\mathbf{e}}|{\mathcal{E}}^{\dagger}{\mathcal{E}}|x_{\mathbf{e}},\bar{x}_{\mathbf{e}}\rangle\prod_{e}(1+x_{e}z_{v}z_{v^{\prime}})(1+\bar{x}_{e}\bar{z}_{v}\bar{z}_{v^{\prime}})(1+x^{\prime}_{e}t_{v}t_{v^{\prime}})(1+\bar{x}^{\prime}_{e}\bar{t}_{v}\bar{t}_{v^{\prime}})\Big{]}.$
(36)

Since $z_{i}z_{j}\bar{z}_{i}\bar{z}_{j}$ is outside the summation of
$\sum_{x_{\mathbf{e}},\bar{x}_{\mathbf{e}}}$, one finds
$\langle\rho_{0}|\mathcal{E}^{\dagger}\mathcal{E}|\rho_{e\bar{e},0}\rangle=\sum_{z_{\mathbf{v}},\bar{z}_{\mathbf{v}},t_{\mathbf{v}},\bar{t}_{\mathbf{v}}}z_{i}z_{j}\bar{z}_{i}\bar{z}e^{-H(z_{\mathbf{v}},\bar{z}_{\mathbf{v}},t_{\mathbf{v}},\bar{t}_{\mathbf{v}})}$,
and thus $\langle I\bar{I}|e\bar{e}\rangle$ corresponds to the two-point
correlation function of $z\bar{z}$ no matter what the channel $\mathcal{E}$
is. A straightforward generalization to other anyon states related to charges
then show that we can identify $(e,\bar{e})$ in
$|\Psi_{0}\rangle|\bar{\Psi}_{0}\rangle$ as $(z,\bar{z})$ and $(e,\bar{e})$ in
$\langle\Psi_{0}|\langle\bar{\Psi}_{0}|$ as $(t,\bar{t})$.

On the other hand, a pair of $m\bar{m}$ ayons located at sites $\tilde{i}$ and
$\tilde{j}$ on the dual lattice can be created by
$w_{m}(\tilde{l})\bar{w}_{m}(\tilde{l})=\prod_{e\in\tilde{l}}Z_{e}\bar{Z}_{e}$,
where $\tilde{l}$ can be any string with end points located at $\tilde{i}$ and
$\tilde{j}$. It follows that

$\displaystyle|\rho_{m\bar{m},0}\rangle$
$\displaystyle=\prod_{e\in\tilde{l}}Z_{e}\bar{Z}_{e}|\rho_{0}\rangle=\sum_{x_{\mathbf{e}},\bar{x}_{\mathbf{e}}}[\sum_{z_{\mathbf{v}},\bar{z}_{\mathbf{v}}}\prod_{e}(1+x_{e}z_{v}z_{v^{\prime}})(1+\bar{x}_{e}\bar{z}_{v}\bar{z}_{v^{\prime}})]|\eta_{\mathbf{e}}x_{\mathbf{e}},\eta_{\mathbf{e}}\bar{x}_{\mathbf{e}}\rangle$
(37)
$\displaystyle=\sum_{x_{\mathbf{e}},\bar{x}_{\mathbf{e}}}[\sum_{z_{\mathbf{v}},\bar{z}_{\mathbf{v}}}\prod_{e}(1+\eta_{e}x_{e}z_{v}z_{v^{\prime}})(1+\eta_{e}\bar{x}_{e}\bar{z}_{v}\bar{z}_{v^{\prime}})]|x_{\mathbf{e}},\bar{x}_{\mathbf{e}}\rangle,$

where $\eta_{e}=1(-1)$ if $e\notin\tilde{l}(\in\tilde{l})$. In the second
line, we relabel the variables $x_{\mathbf{e}}(\bar{x}_{\mathbf{e}})$ as
$\eta_{\mathbf{e}}x_{\mathbf{e}}(\eta_{\mathbf{e}}\bar{x}_{\mathbf{e}})$.
Therefore, the effect of $\prod_{e\in\tilde{l}}Z_{e}\bar{Z}_{e}$ on
$|\rho_{0}\rangle$ corresponds to changing the sign of the Ising interaction
for $z_{\mathbf{v}}$ and $\bar{z}_{\mathbf{v}}$. It is then obvious that
$\langle\rho_{0}|\mathcal{E}^{\dagger}\mathcal{E}|\rho_{m\bar{m},0}\rangle=\sum_{z_{\mathbf{v}},\bar{z}_{\mathbf{v}},t_{\mathbf{v}},\bar{t}_{\mathbf{v}}}e^{-H_{z_{\tilde{l}}\bar{z}_{\tilde{l}}}(z_{\mathbf{v}},\bar{z}_{\mathbf{v}},t_{\mathbf{v}},\bar{t}_{\mathbf{v}})}$,
where

$\displaystyle
e^{-H_{z_{\tilde{l}}\bar{z}_{\tilde{l}}}(z_{\mathbf{v}},\bar{z}_{\mathbf{v}},t_{\mathbf{v}},\bar{t}_{\mathbf{v}})}=\sum_{\begin{subarray}{c}x_{\mathbf{e}},\bar{x}_{\mathbf{e}},\\\
x^{\prime}_{\mathbf{e}},\bar{x}^{\prime}_{\mathbf{e}}\end{subarray}}\Big{[}\langle
x^{\prime}_{\mathbf{e}},\bar{x}^{\prime}_{\mathbf{e}}|{\mathcal{E}}^{\dagger}{\mathcal{E}}|x_{\mathbf{e}},\bar{x}_{\mathbf{e}}\rangle\times\prod_{e}(1+\eta_{e}x_{e}z_{v}z_{v^{\prime}})(1+\eta_{e}\bar{x}_{e}\bar{z}_{v}\bar{z}_{v^{\prime}})(1+x^{\prime}_{e}t_{v}t_{v^{\prime}})(1+\bar{x}^{\prime}_{e}\bar{t}_{v}\bar{t}_{v^{\prime}})\Big{]}.$
(38)

Therefore,
$\langle\rho_{0}|\mathcal{E}^{\dagger}\mathcal{E}|\rho_{m\bar{m},0}\rangle/\langle\rho_{0}|\mathcal{E}^{\dagger}\mathcal{E}|\rho_{0}\rangle$
is nothing but the disorder parameter
$\langle\mu^{z}_{\tilde{i}}\mu^{\bar{z}}_{\tilde{i}}\mu^{z}_{\tilde{j}}\mu^{\bar{z}}_{\tilde{j}}\rangle$.
A straightforward generalization to other anyon states related to flux then
show that we can identify $(m,\bar{m})$ in
$|\Psi_{0}\rangle|\bar{\Psi}_{0}\rangle$ as $(\mu^{z},\mu^{\bar{z}})$ and
$(m,\bar{m})$ in $\langle\Psi_{0}|\langle\bar{\Psi}_{0}|$ as
$(\mu^{t},\mu^{\bar{t}})$.

###  A.2 Coherent channel

We now focus on the specific coherent error channel
$\mathcal{E}[\cdot]=\prod_{e}\mathcal{E}_{e}[\cdot]$, where
$\mathcal{E}_{e}[\cdot]=(1-p)(\cdot)+p\sigma_{e}(\cdot)\sigma_{e},\
\sigma_{e}=\cos(\theta)Z_{e}+\sin(\theta)X_{e}$. After C-J map, one finds
$\mathcal{E}=\mathcal{E}(p,\theta)=\prod_{e}[(1-p)+p\sigma_{e}(\theta)\bar{\sigma}_{e}(\theta)]$.
To simplify the calculation of the normalization
$\langle\rho|\rho\rangle=\langle\rho_{0}|\mathcal{E}^{\dagger}\mathcal{E}|\rho_{0}\rangle$,
we note that that since
$\mathcal{E}(p,\theta)^{\dagger}\mathcal{E}(p,\theta)=\mathcal{E}(p^{\prime},\theta)$
with $p^{\prime}=2p(1-p)$, one finds
$\langle\rho_{0}|\rho(p^{\prime},\theta)\rangle=\langle\rho(p,\theta)|\rho(p,\
\theta)\rangle$. Therefore, instead of directly evaluating
$\langle\rho(p,\theta)|\rho(p,\ \theta)\rangle$, we will compute
$\langle\rho_{0}|\rho(p^{\prime},\theta)\rangle$, which is relatively simpler
to work with. Now, a direct substitution of $\mathcal{E}(p^{\prime},\theta)$
gives

$\langle\rho_{0}|\rho(p^{\prime},\theta)\rangle=\langle\rho_{0}|\prod_{e}[[1-2p(1-p)]I+2p(1-p)\sigma_{e}\bar{\sigma}_{e}]|\rho_{0}\rangle\propto\langle\rho_{0}|\prod_{e}(I+\lambda\sigma_{e}\bar{\sigma}_{e})|\rho_{0}\rangle,$
(39)

where $\lambda=2p(1-p)/(1-2p+2p^{2})$ and

$\sigma_{e}(\theta)\bar{\sigma}_{e}(\theta)=\sin^{2}(\theta)X_{e}\bar{X}_{e}+\sin(\theta)\cos(\theta)(X_{e}\bar{Z}_{e}+\bar{Z}_{e}X_{e})+\cos^{2}(\theta)Z_{e}\bar{Z}_{e}.$
(40)

Using the property that $\langle\Psi_{0}|\prod_{e\in
l}X_{e}|\Psi_{0}\rangle=0$
($\langle\Psi_{0}|\prod_{e\in\tilde{l}}Z_{e}|\Psi_{0}\rangle=0$) unless $l$
($\tilde{l}$) forms closed loops in the original (dual) lattice, one finds
that when expanding $\prod_{e}(I+\lambda\sigma_{e}\bar{\sigma}_{e})$, any
terms that include $(X_{e}\bar{Z}_{e}+\bar{Z}_{e}X_{e})$ take zero expectation
value with $|\rho_{0}\rangle=|\Psi_{0}\rangle|\Psi_{0}\rangle$. Therefore, we
only need to consider
$|\rho^{\prime}(\lambda,\theta)\rangle=\prod_{e}[1+\lambda(\sin^{2}(\theta)X_{e}\bar{X}_{e}+\cos^{2}(\theta)Z_{e}\bar{Z}_{e})]|\rho_{0}\rangle$,
and a straightforward algebra gives
$|\rho^{\prime}(\lambda,\theta)\rangle=\sum_{x_{e},\bar{x}_{e}}[\sum_{z_{v},\bar{z}_{v}}\prod_{e}\nu_{e}(z_{v},\bar{z}_{v})]|x_{e},\bar{x}_{e}\rangle$,
where

$\displaystyle\nu_{e}=(1+\lambda\sin^{2}(\theta)x_{e}\bar{x}_{e})(1+x_{e}z_{v}z_{v^{\prime}})(1+\bar{x}_{e}\bar{z}_{v}\bar{z}_{v^{\prime}})+\lambda\cos^{2}(\theta)(1-x_{e}z_{v}z_{v^{\prime}})(1-\bar{x}_{e}\bar{z}_{v}\bar{z}_{v^{\prime}}).$
(41)

Now, using
$\langle\rho_{0}|x_{e},\bar{x}_{e}\rangle=[\sum_{t_{v},\bar{t}_{v}}\prod_{e}(1+x_{e}t_{v}t_{v^{\prime}})(1+\bar{x}_{e}\bar{t}_{v}\bar{t}_{v^{\prime}})]$,
one finds
$\langle\rho_{0}|\rho(\lambda,\theta)\rangle=\sum_{x_{e},\bar{x}_{e}}\langle\rho_{0}|x_{e},\bar{x}_{e}\rangle\langle
x_{e},\bar{x}_{e}|\rho^{\prime}(\lambda,\theta)\rangle\propto\sum_{z_{\mathbf{v}},\bar{z}_{\mathbf{v}},t_{\mathbf{v}},\bar{t}_{\mathbf{v}}}\prod_{e}\omega_{e}(z_{\mathbf{v}},t_{\mathbf{v}})$,
where

$\displaystyle\omega_{e}=$ $\displaystyle
1+z_{v}z_{v^{\prime}}\bar{z}_{v}\bar{z}_{v^{\prime}}t_{v}t_{v^{\prime}}\bar{t}_{v}\bar{t}_{v^{\prime}}+\frac{1-\lambda\cos^{2}(\theta)}{1+\lambda\cos^{2}(\theta)}(t_{v}t_{v^{\prime}}z_{v}z_{v^{\prime}}+\bar{t}_{v}\bar{t}_{v^{\prime}}\bar{z}_{v}\bar{z}_{v^{\prime}})$
(42)
$\displaystyle+\frac{\lambda\sin^{2}(\theta)}{1+\lambda\cos^{2}(\theta)}(t_{v}t_{v^{\prime}}\bar{t}_{v}\bar{t}_{v^{\prime}}+z_{v}z_{v^{\prime}}\bar{z}_{v}\bar{z}_{v^{\prime}}+t_{v}t_{v^{\prime}}\bar{z}_{v}\bar{z}_{v^{\prime}}+\bar{t}_{v}\bar{t}_{v^{\prime}}z_{v}z_{v^{\prime}})$

One can simplify Eq.(42) by noting that
$\omega_{e}z_{v}z_{v^{\prime}}\bar{z}_{v}\bar{z}_{v^{\prime}}t_{v}t_{v^{\prime}}\bar{t}_{v}\bar{t}_{v^{\prime}}=\omega_{e}$
so that $\langle
z_{v}z_{v^{\prime}}\bar{z}_{v}\bar{z}_{v^{\prime}}t_{v}t_{v^{\prime}}\bar{t}_{v}\bar{t}_{v^{\prime}}\rangle=1,\forall
e$. Therefore, we can replace $\bar{t}_{v}\bar{t}_{v^{\prime}}$ as
$z_{v}z_{v^{\prime}}\bar{z}_{v}\bar{z}_{v^{\prime}}t_{v}t_{v^{\prime}}$, and
Eq.(42) becomes

$\omega_{e}\propto
1+\frac{\lambda\sin^{2}(\theta)}{1+\lambda\cos^{2}(\theta)}(z_{v}z_{v^{\prime}}\bar{z}_{v}\bar{z}_{v^{\prime}}+t_{v}t_{v^{\prime}}\bar{z}_{v}\bar{z}_{v^{\prime}})+\frac{1-\lambda\cos^{2}(\theta)}{1+\lambda\cos^{2}(\theta)}t_{v}t_{v^{\prime}}z_{v}z_{v^{\prime}}.$
(43)

Defining the variables $s_{v}=z_{v}\bar{z}_{v},\ \tau_{v}=t_{v}\bar{z}_{v}$,
$\omega_{e}$ can be further simplified to

$\displaystyle\omega_{e}$ $\displaystyle\propto
1+\frac{\lambda\sin^{2}(\theta)}{1+\lambda\cos^{2}(\theta)}(s_{v}s_{v^{\prime}}+\tau_{v}\tau_{v^{\prime}})+\frac{1-\lambda\cos^{2}(\theta)}{1+\lambda\cos^{2}(\theta)}s_{v}s_{v^{\prime}}\tau_{v}\tau_{v^{\prime}}$
(44)

By reexponentiating Eq.(44), one finds

$\omega_{e}(z_{\mathbf{v}},t_{\mathbf{v}})\propto
e^{K(s_{v}s_{v^{\prime}}+\tau_{v}\tau_{v^{\prime}})+K_{4}(s_{v}s_{v^{\prime}}\tau_{v}\tau_{v^{\prime}})}.$
(45)

which is nothing but the Ashkin-Teller model with the coeffecient $K$ and
$K_{4}$ related to $h^{\prime}$ and $\theta$ through

$\displaystyle\frac{\tanh(K)[1+\tanh(K_{4})]}{1+\tanh^{2}(K)\tanh(K_{4})}=\frac{\lambda\sin^{2}(\theta)}{1+\lambda\cos^{2}(\theta)},\
\frac{\tanh^{2}(K)+\tanh(K_{4})}{1+\tanh^{2}(K)\tanh(K_{4})}=\frac{1-\lambda\cos^{2}(\theta)}{1+\lambda\cos^{2}(\theta)}.$
(46)

In particular, when $\theta=\pi/4$, one finds

$\displaystyle\tanh(K)=\frac{2-\sqrt{4-{\lambda}^{2}}}{\lambda},\
\tanh(K_{4})=\frac{2-h\sqrt{4-{\lambda}^{2}}}{2-{\lambda}^{2}}.$ (47)

##  Appendix B Derivation of Eq.(27)

We here derive the compact form of the normalization of the decohered double
state
$|\rho(h,p)\rangle=\mathcal{E}(p)\mathcal{N}(h)\bar{\mathcal{N}}(h)|\Psi_{0},\bar{\Psi}_{0}\rangle$
in Eq.(27), where
$\mathcal{E}(p)=\prod_{e}[(1-p)+p\sigma_{e}\bar{\sigma}_{e}]$,
$\mathcal{N}(p)=\prod_{e}(1+h\sigma_{e})$, and
$\sigma_{e}=(X_{e}+Z_{e})/\sqrt{2}$. We first make use of the property that
$[\mathcal{E},\mathcal{N}]=0$,
$\mathcal{E}(p)^{\dagger}\mathcal{E}(p)=\mathcal{E}(p^{\prime})$ with
$p^{\prime}=2p(1-p)$, and
$\mathcal{N}(h)^{\dagger}\mathcal{N}(h)\propto\mathcal{N}(h^{\prime})$ with
$h^{\prime}=2h/(1+h^{2})$ to write
$\mathcal{N}(h)^{\dagger}\bar{\mathcal{N}}(h)^{\dagger}\mathcal{E}(p)^{\dagger}\mathcal{E}(p)\mathcal{N}(h)\bar{\mathcal{N}}(h)\propto\mathcal{E}(p^{\prime})\mathcal{N}(h^{\prime})\bar{\mathcal{N}}(h^{\prime})$,
where

$\displaystyle\mathcal{E}(p^{\prime})\mathcal{N}(h^{\prime})\bar{\mathcal{N}}(h^{\prime})$
$\displaystyle=\prod_{e}[(1-p^{\prime})+p^{\prime}\sigma_{e}\bar{\sigma}_{e}](1+h^{\prime}\sigma_{e})(1+h^{\prime}\bar{\sigma}_{e})$
(48)
$\displaystyle\propto\prod_{e}\Big{[}(1+h^{\prime}\sigma_{e})(1+h^{\prime}\bar{\sigma}_{e})+\frac{p^{\prime}}{1-p^{\prime}}(\sigma_{e}+h^{\prime})(\bar{\sigma}_{e}+h^{\prime})\Big{]}$
$\displaystyle=\prod_{e}\Big{[}(1+h^{\prime}\sigma_{e})(1+h^{\prime}\bar{\sigma}_{e})+\frac{p^{\prime}}{1-p^{\prime}}h^{\prime
2}(1+\frac{1}{h^{\prime}}\sigma_{e})(1+\frac{1}{h^{\prime}}\bar{\sigma}_{e})\Big{]}$

Here $f(p)=2p(1-p)h^{\prime 2}/(1-2p+2p^{2})$ and $\\{S\\}$ denotes the set of
all distinct collection of edges. To derive Eq.(27), we first note that
$\langle\Psi_{0},\bar{\Psi}_{0}|\prod_{e}(1+h^{\prime}\sigma_{e})(1+h^{\prime}\sigma_{e})|\Psi_{0},\bar{\Psi}_{0}\rangle=\sum_{z_{\mathbf{v}},t_{\mathbf{v}},\bar{z}_{\mathbf{v}},\bar{t}_{\mathbf{v}}}\prod_{e}\omega_{e}(z_{\mathbf{v}},t_{\mathbf{v}};h)\bar{\omega}_{e}(\bar{z}_{\mathbf{v}},\bar{t}_{\mathbf{v}};h)$,
where
$\omega_{e}(z_{\mathbf{v}},t_{\mathbf{v}};h)=\left[\sqrt{2}+h+h(z_{v}z_{v^{\prime}}+t_{v}t_{v^{\prime}})+(\sqrt{2}-h)z_{v}z_{v^{\prime}}t_{v}t_{v^{\prime}}\right]$.
It follows that

$\displaystyle\langle\rho|\rho\rangle$
$\displaystyle\propto\langle\Psi_{0},\bar{\Psi}_{0}|\mathcal{E}(p^{\prime})\mathcal{N}(h^{\prime})\bar{\mathcal{N}}(h^{\prime})|\Psi_{0},\bar{\Psi}_{0}\rangle$
(49)
$\displaystyle=\sum_{z_{\mathbf{v}},t_{\mathbf{v}},\bar{z}_{\mathbf{v}},\bar{t}_{\mathbf{v}}}\prod_{e}\left[\omega_{e}(z_{\mathbf{v}},t_{\mathbf{v}};h)\omega_{e}(\bar{z}_{\mathbf{v}},\bar{t}_{\mathbf{v}};h)+f(p)\,\omega_{e}(z_{\mathbf{v}},t_{\mathbf{v}};h)\omega_{e}(\bar{z}_{\mathbf{v}},\bar{t}_{\mathbf{v}};h)\right].$

##  Appendix C Statistical mechanical mapping for non-fixed-point toric code
under phase-flip channel

In this Appendix, we are interested in subjecting non-fixed point toric code
wave functions to the phase-flip channel
$\mathcal{E}[\cdot]=\prod_{e}\mathcal{E}_{e}[\cdot]$, where
$\mathcal{E}_{e}[\cdot]=(1-p)\rho+pZ_{e}\rho Z_{e}$. Specifically, we choose
the initial pure state as the toric code wavefunction evolved by a single-site
non-unitary operator with a tuning parameter $h\in[0,1]$:

$|\Psi(h)\rangle=\prod_{e}(1+hZ_{e})|\Psi_{0}\rangle\propto
e^{K\sum_{e}Z_{e}}|\Psi_{0}\rangle,$ (50)

where $K=\log(1+h)-\log(1-h)$. We note that Ref.[51] has shown that
$|\Psi(h)\rangle$ is a topologically ordered state when
$h<h_{c}=1+\sqrt{2}+\sqrt{2(1+\sqrt{2})}\approx 0.217$ while is a short-range-
entangled state when $h>h_{c}$. On the other hand, Ref.[19] has shown that
when $h=0$, the $n$-th moment of the decohered density matrix
$\rho=\mathcal{E}[|\Psi(0)\rangle\langle\Psi(0)|]$ can be mapped to $(n-1)$-th
flavored Ising model. In this Appendix, we will show that for any $p$ and $h$,
the $n$-th moment of the density matrix can be mapped to the $n$-flavored
Ising model with the $s$-th flavored spin coupled to the $(s+1)$-th (mod $n$)
flavored spin through four-spin interactions. We will also comment about what
we expect when taking the limit $n\rightarrow 1$.

To see how the initial pure state density matrix
$\rho_{h}=|\Psi(h)\rangle\langle\Psi(h)|$ evolves under the phase-flip
channel, it is more convenient to express $|\Psi(h)\rangle$ in terms of the
loop configurations using the fact that
$|\Psi_{0}\rangle\propto\sum_{g}g|0\rangle\propto\sum_{g}|g\rangle$ with $g$
denotes the closed-loop Pauli-X operators. A repeated use of
$e^{KZ_{e}}X_{e}|0\rangle=X_{e}e^{-KZ_{e}}|0\rangle=e^{-K}X_{e}|0\rangle$
shows that $|\Psi(h)\rangle=\sum_{g}e^{-K|g|}|g\rangle/\sqrt{\mathcal{Z}}$
where $\mathcal{Z}=\sum_{g}e^{-2K|g|}$ and $|g|$ denotes the total length of
the given loops configuration. Therefore, the initial pure state density
matrix takes the form

$\rho=\frac{1}{\mathcal{Z}}\sum_{g,\bar{g}}e^{-K(|g|+|\bar{g}|)}|g\rangle\langle\bar{g}|.$
(51)

The decohered density matrix after subjecting to channel
$\mathcal{E}_{e}[\cdot]=(1-p)\rho+pZ_{e}\rho Z_{e}$ can then be computed using
the similar method as $|\Psi(h)\rangle$, and the result is

$\rho=\frac{1}{\mathcal{Z}N_{G}}\sum_{g,\bar{g}}e^{-K(|g|+|\bar{g}|)-K_{4}|g\bar{g}|}|g\rangle\langle\bar{g}|,$
(52)

where $K_{4}=-\log(1-2p)$ and $N_{G}$ denotes the number of all possible loops
configurations. One can now write the $n$-th moment density matrix as

$\tr(\rho^{n})=\frac{1}{(\mathcal{Z}N_{G})^{n}}\sum_{\\{g^{(s)}\\}}e^{-2K\sum_{s=1}^{n}|g^{(s)}|-K_{4}\sum_{s=1}^{n}|g^{(s)}g^{((s+1)_{n})}|},$
(53)

where $(s+1)_{n}$ denotes the sum modulo $n$. Eq.(53) can be regarded as $n$
different flavors of loop models, and the given flavor $s$ coupled with the
flavor $(s+1)_{n}$ through the relative loop configurations. One can also map
it to the spin model by identifying the flavor $s$ loop as the domain walls of
the flavor $s$ Ising spins, and Eq.(53) can be written as
$\tr(\rho^{n})=\sum_{\\{\sigma^{(s)}_{\mathbf{v}}\\}}e^{-H_{n}}/\mathcal{Z}^{n}$,
where

$-H_{n}=\sum_{\langle
v,v^{\prime}\rangle}\sum_{s=1}^{n}\big{(}2K\sigma^{(s)}_{v}\sigma^{(s)}_{v^{\prime}}+K_{4}\sigma^{(s)}_{v}\sigma^{(s)}_{v^{\prime}}\sigma^{(s+1)_{n}}_{v}\sigma^{(s+1)_{n}}_{v^{\prime}}\big{)}.$
(54)

Therefore, we have mapped the $n$-th moment density matrix as the n-flavored
Ising model with flavor $s$ spin coupled with the flavor $(s-1)_{n}$ and
$(s+1)_{n}$ through the four-spin interactions. We note that for $n=2$ and
$3$, $H_{n}$ is equivalent to the 2-flavored and 3-flavored [57] Ashkin-Teller
(AT) model. However, we are not aware of any previous work that has studied
$H_{n}$ for $n\geq 4$.

We now focus on the case where $n=2$, whose error threshold is expected to be
the closest to the $n\rightarrow 1$ limit compared to other integer Rényi
indices. The exact mapping to the isotropic Ashkin-Teller model allows us to
determine the phase boundary and the universality of the transition as shown
in Fig.5(a). One can similarly use the anyon condensation and confinement
paremeters mentioned in the main text to identify three phases: the mixed-
state topologically ordered phase, the classical topologically ordered phase,
and the trivial phase. We note that the topologically ordered phase is mapped
to the paramagnetic phase of the AT model, which is different from the mapping
for the fixed-point toric code subjected to coherent channel. Another
interesting feature of the phase diagram is that when $h$ is large but below
the critical point of the pure-state transition, one can increase $p$ and
drive the system from topological to the trivial phase through the BKT
transition. However, we believe this phenomena is due to the replica index and
will disappear when $n\rightarrow 1$, since the Ising2 transtion along the
$p=0$ line (i.e., the pure-state transition line) becomes the usual Ising
transtion when $n\rightarrow 1$.

Finally, we breifly comment on the expected phase diagram in the $n\rightarrow
1$ limit. In this limit, both the $p=0$ and $p=0.5$ lines can be exactly
mapped to the 2d classical Ising model, and crucially, they have the same
critical field strength $h_{c}=1+\sqrt{2}+\sqrt{2(1+\sqrt{2})}\approx 0.217$.
Therefore, we expect a vertical boundary with 2d Ising transition at
$h=h_{c}$. On the other hand, the $h=0$ limit can be exactly mapped to the
random-bond (RB) Ising model along the Nishimori line, where the transtion
occurs at $[=p_{c}\approx 0.109$. Therefore, one expects a smooth line with RB
Ising transition connecting $(h,p)\approx(0,0.109)$ and
$(h,p)\approx(0.217,0)$. We summaize the expected $n\rightarrow 1$ phase
diagram in Fig.5 in the main text.