Relaxation dynamics of the toric code in contact with a thermal reservoir: finite-size scaling in a low temperature regime

The toric code provides a simple exactly soluble model with a topologically ordered ground state that may provide topologically protected qubits at $T=0$ Kitaev (2003); Dennis et al. (2002). The toric code is defined on a square lattice, where Ising spins sit on the links of the lattice. We define the linear dimension of the lattice as $L$ and the number of spins $N=2L^{2}$. The Hamiltonian involves four-spin interactions around the plaquettes and vertices of the lattice:

where the sums over $v$ and $p$ are over the vertices and plaquettes of the lattice, respectively (see Fig. 1). The ground states are the $+1$ eigenstate of all $A_{v}$ and $B_{p}$ operators, since all such operators commute. On a torus, there are four degenerate ground states that are distinguished by the expectation values of non-local winding operators $W_{1,2}^{x}$,$W_{1,2}^{z}$:

where $\Gamma_{1,2}$ and $\tilde{\Gamma}_{1,2}$ are topologically non-trivial loops along the links and plaquettes of the lattice, respectively, that wind around each of the two axes of the torus. There is a finite gap $\Delta_{e,m}=4J_{e,m}$ to excited states that are $-1$ eigenstates of some $A_{v}$ and/or $B_{p}$. These correspond to $e$-type and $m$-type anyonic quasiparticle excitations, respectively Kitaev (2003).

Consider the $W_{1,2}^{z}$ basis for the degenerate ground state subspace; we may label the four ground states by the eigenvalues of $W_{1,2}^{z}$:

where $+/-$ represent the $\pm 1$ eigenstates of $W_{1}^{z}$ and $W_{2}^{z}$, respectively. We choose this basis to be the logical basis for a two qubit quantum memory. To simplify the discussion of errors we will consider the limit $J_{m}\rightarrow\infty$, so that only $e$-type quasiparticles have finite energy. The lowest excited eigenstates have a single pair of localized $e$-type quasiparticles which are connected to a ground state by operation of an error string:

Here, $v$ and $v^{\prime}$ are the vertices where the quasiparticles are located and the string $\Gamma_{v,v^{\prime}}$ connects $v$ and $v^{\prime}$. Error strings that form topologically non-trivial loops generate the $W_{1,2}^{x}$ operators and drive transitions between ground states; such error strings create noncorrectable errors, i.e., errors in the logical subspace that cannot be corrected. Correctable errors, or self-correcting errors, are those error strings that close without causing a change in winding number.

Under local perturbations to $H_{\mathrm{TC}}$, such topologically nontrivial error strings only occur at order $L$ in perturbation theory; consequently both the splitting of the ground state degeneracy and transitions between ground states are suppressed exponentially, and thus this ground state subspace is “topologically protected” from such unitary perturbations Klich (2010); Bravyi et al. (2010); Bravyi and Hastings (2011); Trebst et al. (2007); Tupitsyn et al. (2010); Dusuel et al. (2011); Michalakis and Zwolak (2013); Kay (2011). The toric code can thus act as a self-correcting quantum memory at zero temperature.

Despite the topological protection at zero temperature, in the thermodynamic limit the topological order of the toric code is destroyed at any finite temperature Nussinov and Ortiz (2008); Castelnovo and Chamon (2007). Consequently, a topological qubit would be thermally fragile. While topologically nontrivial error strings due to unitary perturbations are exponentially suppressed in the system size, non-trivial error strings may also be generated by non-unitary perturbations, e.g., from contact with a thermal reservoir. Non-correctable errors due to non-unitary perturbations are not exponentially suppressed in the system size.

The thermal fragility of the 2D toric code can be understood from a simple picture of the dissipative dynamics generated from local interactions with an external bath. A local system-bath interaction can generate a trivial error string by flipping a single spin, thus creating a single pair of neighboring quasiparticles:

where $v$ and $v^{\prime}$ are the vertices on either end of the link $j$. The rate of such a process is suppressed exponentially in the inverse temperature, due to the energy gap $\Delta$ to such excited states. Such a trivial error is correctable by applying another $\sigma_{j}^{x}$. However, additional trivial error strings $\sigma_{j^{\prime}}^{x}$ with $j\neq j^{\prime}$ applied to $v$ or $v^{\prime}$ will generate a longer, non-trivial error string at no energy cost. Consequently, local coupling to a bath can drive a random walk of quasiparticle pairs around the lattice with a rate that is only suppressed by a single Boltzmann factor. Such random walks may generate topologically nontrivial error loops and return the system back to the ground state subspace. If an error loop has an odd winding number, this error loop has caused a non-correctable error by driving a transition between ground states. Alternatively, if the error loop has an even winding number, the error is self-correcting. Indeed, Alicki et al. have placed a system-size-independent upper bound on the relaxation time of a pure toric code ground state that explicitly demonstrates this thermal fragility of the toric code in the thermodynamic limit Alicki et al. (2007, 2009). Additionally, the analysis of Nussinov and Ortiz demonstrates the lack of “spontaneous topological symmetry breaking” at finite-temperature, as the autocorrelation time of the winding operators is sub-exponential in lattice dimension in the thermodynamic limit Nussinov and Ortiz (2008, 2009).

While the toric code is thermally fragile at all non-zero temperatures in the thermodynamic limit, we can also consider how this fragility is affected by the finite size of a lattice. As outlined above, the dissipative error processes which lead to thermal fragility occur when there is a single quasiparticle pair present. Since the number of excitations in equilibrium is suppressed by the Boltzmann factor at low temperatures, we expect that at sufficiently low temperatures the number of quasiparticles in equilibrium will be vanishingly small. We can then define an equilibrium crossover temperature $T_{eq}^{*}$ which distinguishes the thermally fragile regime from a low temperature regime with reduced dissipative error processes by:

Castelnovo and Chamon define an equilibrium crossover temperature $T_{eq}^{*}$ above which the topological entanglement entropy vanishes Castelnovo and Chamon (2007). They find that this equilibrium definition of $T^{*}$ scales inversely with the log of the system size, and that this becomes a zero temperature phase transition in the thermodynamic limit. Conversely, on a finite sized system the crossover temperature $T^{*}$ defines a low temperature regime where topological order persists as a finite-size effect. This opens the possibility of using finite-size effects to exploit the zero temperature topological order at finite temperatures. The usefulness of this low temperature regime for quantum information processing depends on the scaling of the relaxation time in this regime. Robustness to unitary perturbations requires a sufficiently large system size to minimize the splitting of the degeneracy and the matrix elements between ground states, while the thermal fragility increases with system size. Thus one may expect that there is an optimal size for computational performance.

Below, we directly address this question of the finite-size scaling of the relaxation time of a toric code ground state in contact with a thermal reservoir. This analysis complements the growing literature concerning active error correction on the toric code by use of the stabilizer space and associated stabilizer operations Dennis et al. (2002); Bravyi and Terhal (2009); Chesi et al. (2010a); Jouzdani et al. (2013). Usually, stabilizer error analysis is considered in the context of effectively infinite temperature thermal instabilityChesi et al. (2010b), in contrast to the finite-temperature dynamics presented here. A complete toolkit for understanding and controlling errors in physical implementations of the toric code would need a predictive low temperature model, as well as a recipe for understanding how the fidelity of stabilizer operations affects the finite temperature operation of the toric code. Here we focus on the robustness of the passive error correcting (i.e., self-correcting) dynamics under the action of the toric code Hamiltonian in contact with a thermal reservoir.

We present a microscopic model of the real-time non-equilibrium dynamics of the toric code in contact with a thermal reservoir. Due to the fact that the spectrum of $H_{\mathrm{TC}}$ has a finite gap to excited eigenstates with localized quasiparticle excitations, such dynamics may be described by a Lindblad master equationAlicki et al. (2007, 2009):

here $\rho$ is the toric code system density matrix and $\{c_{\omega}\}$ is a set of Lindblad operators generated by local system-bath interactions. We will consider the limit where $J_{m}\gg J_{e}$, such that at low temperatures the system will remain in the +1 eigensector of all $B_{p}$ operators and only $e$-type quasiparticles will be excited by the reservoir. We will only consider local system-bath couplings, for which the bath generates single spin flips in the system. The relevant Lindblad operators are:

where $E^{e\dagger}_{vv^{\prime}}$ ($E^{e}_{vv^{\prime}}$) creates (annihilates) a pair of quasiparticles at neighboring vertices and $T^{e}_{vv^{\prime}}$ translates a quasiparticles across a link. These operators are defined by:

Since the Lindblad form of the master equation only connects diagonal elements of the density matrix $\rho$ to other diagonal elements, expectation values of diagonal elements will evolve independently of off-diagonal density matrix elements. Correspondingly, the time evolution of diagonal matrix elements reduces to a classical master equation:

where $n$ labels an eigenstate of $H_{\mathrm{TC}}$, $P_{n}=\rho_{nn}$ are the diagonal matrix elements (probabilities) and $\{\gamma_{0},\gamma_{+},\gamma_{-}\}$ are the rates at which the operators $\{T^{e},E^{e\dagger},E^{e}\}$ act, respectively. Similarly, in (11), for a given $n$, the indices of $n_{0}$, $n_{+}$, and $n_{-}$ label the sets of eigenstates connected to $|n\rangle$ by the operators $T^{e}$, $E^{e\dagger}$, and $E^{e}$, respectively. The ratio $\gamma_{+}/\gamma_{-}$ is fixed by detailed balance to be

but the nature of the bath and the coupling strength determines $\gamma_{0}$ and the magnitude of $\gamma_{-}$ (or equivalently $\gamma_{+}$).

We consider here an Ohmic, Markovian bath with a power spectrum given by:

with ${\omega}_{c}$ a high frequency cutoff much larger than $J_{e}$. Taking $w_{c}$ to infinity gives rise to decay rates of the formChesi et al. (2010b):

where $\xi$ sets the strength of the phenomenological system-bath coupling. This leads to the following rates:

We are most interested in the dynamics deriving from the initial condition of a pure ground state. We characterize the relaxation from a pure ground state by considering the time evolution of the expectation value of the winding operators:

The population dynamics are governed by the creation of quasiparticle pairs that undergo random walks on the torus and then annihilate. Thermal transitions between ground states occur when the quasiparticle pair undergoes a topologically non-trivial random walk before annihilating. The decay of the expectation values in (16) from their values in a pure state with eigenvalue $\pm 1$ can be due to both topologically nontrivial random walks generating transitions between ground states, as well as transitions to excited states via propagating, open error strings. Consequently, the statistics of such topological random walks affect the scaling of the lifetime of a ground state.

Nussinov and Ortiz showed that one can take advantage of the equivalence of the partition function of the toric code and that of 1D classical Ising chains to compute equilibrium properties of the toric code Nussinov and Ortiz (2008, 2009). However, one can not readily take advantage of this mapping for the study of non-equilibrium properties. While the partition function is only a function of the spectrum of the system, non-equilibrium dynamics depend on the nature of the (local) coupling to the external reservoir. Since the mapping of the toric code to an Ising chain maps a 2D model onto a 1D model, local couplings of the toric code to an external reservoir in general can lead to non-local couplings in the corresponding Ising model. Thus, a simple model of the 1D non-equilibrium dynamics of the Ising chains locally coupled to an external bath cannot describe the non-equilibrium dynamics of the toric code with a local bath coupling. Fundamentally, the dynamics of each system at low temperatures are governed by the random walk of defects (anyons in the toric code, domain walls in the Ising models); the defects of the Ising model undergo 1D random walks, whereas those of the toric code undergo 2D random walks. Consequently, we cannot directly compute the finite-size relaxation times of the toric code ground states from an analysis of the dynamics of the Ising chain. It is nevertheless useful to discuss the nature of thermal relaxation in a finite-sized Ising chain to help inform our discussion of such dynamics in the toric code.

Consider a periodic 1D chain of $L$ classical Ising spins $s_{i}=\pm 1$ with energy

where $J>0$ is a ferromagnetic coupling constant. The ground state of (17) is a ferromagnet, but the long range order is destroyed at all nonzero temperatures. Excitations above the degenerate ground states are pairs of domain walls with energy cost $\Delta=4J$. If a pair of domain walls undergoes a topologically non-trivial 1D random walk, this drives a thermal transition between ground states, which is we will refer to as “ground state relaxation”. At sufficiently low temperatures on a finite-sized system, we can expect that these 1D topologically nontrivial walks will dominate the relaxation time of the magnetization. Such a low temperature regime must occur when there is less than a single pair of defects in equilibrium:

At low enough temperatures, there may be a separation of time scales such that $\gamma_{+}\ll\gamma_{0}\ll\gamma_{-}$. Intuitively, the domain wall production rate, $\gamma_{+}$, can be tuned much less than the domain wall annihilation rate, $\gamma_{-}$, simply by lowering the temperature (c.f. (12)). For certain choices of bath model, the domain wall hopping rate, $\gamma_{0}$, can be tuned between the latter two rates. On a finite size lattice the time scale for an extensive random walk of the defects can be estimated by the diffusion equation to be of the order of $L^{2}/\gamma_{0}$. We consider the low temperature regime of a finite size lattice where $\gamma_{+}\ll\gamma_{0}/L^{2}$. In this limit, the rate of topologically nontrivial walks occurring is determined by the rate of production of defect pairs and by the probability that such pairs undergo a topologically nontrivial walk before annihilating. This is because any domain wall pairs which proceed to an extensive random walk will carry out their walk and annihilate much faster than another pair of defects will be created.

We consider only the lowest order processes at first order in $\gamma_{+}$. Once a defect pair is created, the probability of the pair separating instead of trivially annihilating is of the order of $\gamma_{0}/\gamma_{-}$. The lowest order processes will annihilate upon their first return to neighboring links. Processes for which the defect pair do not annihilate after the first return to neighboring links will occur at higher order in $\gamma_{0}/\gamma_{-}$. Consequently the overall order of these lowest order processes is $\gamma_{+}\gamma_{0}/\gamma_{-}$. We may then introduce a phenomenological form of the relaxation rate from a ferromagnetic ground state:

The linear scaling in $L$ arises from the number of locations for domain wall pairs to be created. The factor $P^{\Omega}_{{\rm 1D}}(L)$ is the probability that any given domain wall that does not immediately annihilate eventually undergoes a topologically nontrivial random walk. A symmetry argument shows that $P^{\Omega}_{{\rm 1D}}(L)\sim L^{-1}$ (see Appendix B). This linear scaling of $P^{\Omega}_{{\rm 1D}}(L)$ suggests that the relaxation rate of the Ising model is size independent: ${\Gamma}_{\rm{Ising}}\sim\gamma_{0}\cdot e^{-\Delta/T}$.

In Ref. Glauber, 1963 Glauber solved the exact dynamics of the Ising chain in contact with a thermal reservoir. Glauber uses a bath where $\gamma_{0}$ is taken to be a constant and

The relaxation time of this model is found to be

At low temperatures, $\Gamma_{\rm{Glauber}}\approx\gamma_{0}e^{-\Delta/T}$, in agreement with the expected size independent form of the low temperature single defect pair model(19), despite the fact that $\gamma_{0}/\gamma_{-}\ll 1$ is not satisfied.

We can now make an analogous argument for the toric code. At sufficiently low temperatures on a finite-sized lattice, the relaxation rate from a ground state should be dominated by the dynamics of a single quasiparticle pair. Transitions between ground states are generated by pairs of excitations that annihilate after undergoing a topologically non-trival 2D random walk with an odd winding number. The low temperature regime dominated by single defect pairs occurs when

We consider a separation of time scales for the Ohmic bath defined by (15):

In this regime, the lowest-order processes are of the order of,

and we write a phenomenological ground state relaxation rate of the form,

Analogous to the phenomenological relaxation rate for the Ising model (i.e., (19)), the term $\xi\cdot T\cdot e^{-\Delta/T}\cdot L^{2}$ encodes the rate at which free quasiparticle pairs are produced. The number of spins where a defect pair can be created is $N=2L^{2}$. The scaling of the topological factor $P^{\Omega}_{{\rm 2D}}(L)$, which is the probability of a 2D topologically nontrivial walk (i.e., a walk with odd winding number) will control finite-size scaling of $\Gamma_{\rm{TC}}$; only if $P^{\Omega}_{{\rm 2D}}(L)\sim L^{-2}$ will the relaxation rate of the toric code be system size independent, as for the classical Ising chain. Note that $P^{\Omega}_{{\rm 2D}}\left(L\right)$ is, in general, a function of temperature (see Appendix A).

At high temperatures, the relaxation of a toric code ground state will be dominated by the growing population of quasiparticles, rather than the dynamics of a single quasiparticle pair. In this regime, $\gamma_{+}\approx\gamma_{-}\approx\gamma_{0}$ and the relaxation rate, i.e., the rate of decay of the expectation values $\langle W_{1,2}^{z}\rangle$, is due to error strings created across the length of the $W_{1,2}^{z}$ operator. Consequently we expect the decay to be linear in system size:

This linear scaling arises from the short time dynamics of the master equation (8) Viyuela et al. (2012) and is independent of the topological processes that contribute to $P^{\Omega}_{{\rm 2D}}(L)$ and dominate the low temperature regime.

We will now use the form of the scaling of non-trivial annihilation probabilities discussed above to construct an effective low temperature minimal Markov model of the toric code ground state subspace. The ground state Markov model is defined by the master equation

where $\boldsymbol{P}=(P_{++},P_{+-},P_{-+},P_{--})$ is the vector of all ground state probabilities and $\boldsymbol{\Gamma}$ is the matrix of transition rates between ground states.

Here we assume the low temperature form of the transition rate between ground states to take the form of (26); correspondingly the matrix elements of $\boldsymbol{\Gamma}$ take the form

where $\Gamma_{ij}$ corresponds to the transition $i\rightarrow j$, $\lambda$ is a rate of production of anyon pairs that undergo a nontrivial random walk, and $P_{ij}$ is the probability that a given anyon pair will undergo a topologically nontrivial walk causing the transition $i\rightarrow j$. We take the form of $\lambda$ to be:

where $2L^{2}\gamma_{+}$ is the rate of pair creation for the entire lattice. The remaining factor in brackets is exactly the probability that an adjacent pair of quasiparticles on an otherwise empty lattice does not annihilate. The numerical factors (i.e., $6,(2L^{2}-7),1$) simply index the number of edges that can be acted upon by the different Lindblad operators for the lattice configuration with a single pair of adjacent quasiparticles. By detailed balance, the probability of a given Lindblad operator acting on the system (e.g., $E^{e}$) is then just the ratio of the rate of that operator (e.g., $\gamma_{-}$) to the weighted sum of the rates of the other available operators, weighted by the number of edges available to each operator (e.g., $6\gamma_{0}+\left(2L^{2}-7\right)\gamma_{+}+\gamma_{-}$). Thus, $\lambda$ accounts for the creation rate of quasiparticle pairs that do not immediately annihilate, or those pairs which can generate nontrivial random walks.

The form of $P_{ij}^{\Omega}$ is determined by whether the matrix element is relating ground states that differ by a winding on one axis of the torus or on both. We define $P^{\Omega}_{\delta 1}$ and $P^{\Omega}_{\delta 2}$ to be the probabilities of a topologically nontrivial annihilation that has an odd winding about one axis and both axes of the torus, respectively. Then the form of $P_{ij}$ is:

We may solve this Markov model exactly by integrating (28): for $P_{++}(t=0)=1$ we obtain:

When $P^{\Omega}_{\delta 1}\approx P^{\Omega}_{\delta 2}\approx P^{\Omega}_{\rm{2D}}$, we see that all $P(t)$ are well described by an exponential decay with rate $4P^{\Omega}_{\rm{2D}}(L)\lambda$ at a finite temperature, $T$ (see Appendix A).

As discussed above, at sufficiently low temperatures on a finite size lattice, we expect the relaxation time of a toric code ground state to depend on the statistics of topologically nontrivial random walks on the torus. In this section we present a numerical study of discrete random walks on a square lattice on a torus using Monte Carlo simulations. Without loss of generality, we may map the processes of pair creation, two-particle random walk, and annihilation to a single random walker undergoing a random walk that starts and ends at the origin. We can compute the probability of a quasiparticle pair generating a transition between ground states after annihilation from the statistics of topologically non-trivial walks of the single walker with odd winding.

To estimate the scaling of $P^{\Omega}_{{\rm 2D}}(L)$, we consider a related quantity: the probability that two random walkers will annihilate after $n$ steps $p(n)$. Topological walks must have radius of $L$; given that the radius of a 2D random walk scales as $\sqrt{n}$, we may assume that topological walks hve a minimum number of steps that scales as $n_{\rm{topo}}\sim L^{2}$. $P^{\Omega}_{{\rm 2D}}(L)$ may then be estimated as

This rough estimate assumes that all walks larger than a certain length are necessarily topologically nontrivial.

Restricting to a planar square lattice with trivial topology, the annihilation probability $p^{\rm{p}}(n)$ can be computed to give the asymptotic behavior for large $n$ as Montalenti and Ferrando (2000):

For small $n$, the exact result may be computed numerically via a recursion relation Montalenti and Ferrando (2000). We then can estimate the scaling of $P^{\Omega}_{{\rm 2D}}$ by integrating the planar result:

We expect then that the poly-log scaling of $p^{\rm{p}}(n)$ will lead to an inverse logarithmic finite-size scaling of $P^{\Omega}_{{\rm 2D}}$. Below we compute $P^{\Omega}_{{\rm 2D}}$ numerically and demonstrate this finite-size scaling empirically.

In the low temperature limit, quasiparticle dynamics are dominated by trivial events where pairs that are created and immediately annihilate as $\gamma_{-}\gg\gamma_{0}$. Additionally, for nontrivial walks, once the anyon pair returns to occupy nearest neighbor sites on the lattice, the pair will annihilate with high probability. To model this low temperature regime, we consider an annihilation to occur as soon as a single random walker returns to a site adjacent to the origin Fig. 4. This can be understood as a “zero temperature” limit to the true quasiparticle statistics, as it ignores higher order processes that occur at finite temperature involving quasiparticle trajectories that meet in annihilation geometries, but then do not annihilate. Explicitly, this approximation amounts to taking $\gamma_{-}\rightarrow\infty$. Additionally, to improve efficiency of the Monte Carlo simulations, we start the walker at one of eight starting positions away from the origin; we account for the relative probabilities for reaching these starting positions via exact enumeration of the combinatorics of short topologically trivial walks (see Fig. 4). The random walker undergoes a discrete time random walk on the square lattice on a torus until it returns to one of the four vertices adjacent to the origin.

Using this approach, we compute the probability that two random walkers will annihilate after $n$ steps, $p^{\rm{t}}(n)$, and the average number of steps before annihilation, $\langle n\rangle$. For the true finite temperature toric code, the annihilation probability for nearest neighbor quasiparticles is less than one, as it is a function of $\gamma_{0}/\gamma_{-}$. The finite temperature probabilities $P^{\Omega}_{\rm{2D}}(L)$ may be computed from the zero temperature limit via a “resummation” method described in Appendix A.

Fig. 5 shows the annihilation probability on a torus $p^{\rm{t}}(n)$ as a function of the number of steps $n$ for the zero temperature model with several system sizes $L$, as computed via Monte Carlo. We see that $p^{\rm{t}}(n)$ agrees with $p^{\rm{p}}(n)$ up to a certain value of $n$ for each lattice size. We can therefore define a characteristic “departure time” $n_{d}(L)$, as:

Random walks that annihilate at small $n$ are not sensitive to the topology of the finite-size lattice. Thus, $n_{d}$ reflects the characteristic number of steps at which the random walk distribution is affected by the finite size torus topology. For $n>n_{d}$, we see that $p^{\rm{t}}(n)>p^{\rm{p}}(n)$ up to a characteristic number of steps. We therefore define the “crossing time” $n_{c}$ where $p^{\rm{t}}(n)$ crosses $p^{\rm{p}}(n)$ and then drops significantly. Fig. 6 shows the scaling of both dynamical quantities, $n_{c}$ and $n_{d}$ as a function of system size $L$. Both are seen to be well described by power laws:

with $\alpha_{c}=2.343\pm 0.001$ and $\alpha_{d}=1.66\pm 0.04$.

We also compute the initial and final topological sectors of each walk; this allows us to compute the probability of topologically nontrivial annihilation, $P^{\Omega}_{{\rm 2D}}(L)$, where the walk generates a topologically non-trivial path with odd winding. Fig. 7 shows the finite size scaling of $P^{\Omega}_{{\rm 2D}}(L)$; for larger system sizes, we find

where $c^{\Omega}_{{\rm 2D}}=0.472\pm 0.003$.

We may understand the scaling of $P^{\Omega}_{{\rm 2D}}(L)$ from an analysis of $p^{\rm{t}}(n)$. Since the deviations of $p^{\rm{t}}(n)$ from $p^{\rm{p}}(n)$ for $n>n_{d}$ are due to topologically nontrivial walks which contribute to $P^{\Omega}_{{\rm 2D}}(L)$, we can place approximate bounds on $P^{\Omega}_{{\rm 2D}}(L)$ from $p^{\rm{t}}(n)$. As an upper bound to $P^{\Omega}_{{\rm 2D}}(L)$, we assume that all walks for $n<n_{d}$ generate topologically nontrivial windings that contribute to $P^{\Omega}_{{\rm 2D}}(L)$; therefore we define the integrated probability:

As both $p^{\rm{p}}$ and $p^{\rm{t}}$ integrate to unity, $P_{d}$ is approximately equal to the same integral over $p^{\rm{t}}$. $P_{d}$ includes both topologically trivial walks and topologically nontrivial walks that have even winding; consequently we expect $P_{d}$ to provide an upper bound to $P^{\Omega}_{{\rm 2D}}(L)$.

Alternately, we can make the approximation that topologically nontrivial walks are only the excess probability for $n_{d}\leq n\leq n_{c}$. By making the assumption that $p^{\rm{t}}(n>n_{c})\approx 0$ (see Fig. 5), we may approximate this excess by:

again relying on the normalization of $p^{\rm{t}}$ and $p^{\rm{p}}$. While $P_{c}$ should over-count the events that contribute to $P^{\Omega}_{{\rm 2D}}(L)$ in the region $n_{d}\leq n\leq n_{c}$ (as only odd winding topological events contribute), the approximation $p^{\rm{t}}(n>n_{c})\approx 0$ leads to an underestimation of $P^{\Omega}_{{\rm 2D}}(L)$, i.e. $P_{c}$ provides a lower bound on $P^{\Omega}_{{\rm 2D}}(L)$.

Fig. 7 shows the finite size scaling of $P_{c}\left(L\right)$ and $P_{d}\left(L\right)$ and confirms that these provide a lower and upper bound to $P^{\Omega}_{{\rm 2D}}(L)$, respectively. Thus as with $P^{\Omega}_{{\rm 2D}}(L)$, we find that $P_{c}\left(L\right)$ and $P_{d}\left(L\right)$ scale as $(\ln L)^{-1}$; the lines in Fig. 7 represent a fit to $(\ln L)^{-1}$.

To see the origin of this $(\ln L)^{-1}$ scaling, we can use the asymptotic form of $p^{\rm{p}}(n)$ and the power law scaling of $n_{c,d}$ from (37), to approximate both $P_{c}$ and $P_{d}$:

Consequently, we can understand the origin of the $(\ln L)^{-1}$ scaling which we predicted for $P^{\Omega}_{{\rm 2D}}(L)$ to fundamentally be due to the particular form of the polylog scaling of $p^{\rm{p}}(n)$. We note that this inverse logarithmic scaling of $P^{\Omega}_{{\rm 2D}}(L)$ implies a nontrivial finite-size scaling of ${\Gamma}_{\rm{TC}}$; indeed for the phenomenological form from Eq. (26) we have ${\Gamma}_{\rm{TC}}\sim L^{2}/\ln L$.

We now present Monte Carlo simulations of the real time dynamics of the toric code in contact with an Ohmic bath, as described in section III.1. We use a continuous real time Monte Carlo method Chesi et al. (2010a) to numerically solve the master equation given in (11). We focus on the relaxation dynamics of the system when prepared initially in a pure ground state. We define the operator:

which is one for the $|\Psi_{0}^{++}\rangle$ ground state and vanishes for all other ground states. We then compute the expectation value $\langle\Pi_{++}\left(t\right)\rangle$ to study the decay from $|\Psi_{0}^{++}\rangle$.

The exponential nature of the decay of $\Pi_{++}(t)$ is displayed in Fig. 8, where the lines represent exponential fits to the Monte Carlo data. We find such decays are well described by exponential decays at all but intermediate temperatures (see Fig. 8 c.), where short time deviations lead to a stretched exponential decay. We fit $\Pi_{++}(t)$ to an exponential:

to extract the relaxation rate $\Gamma_{++}$. The additional systematic uncertainty of the rate $\Gamma_{++}$ due to the stretched exponential behavior in intermediate regimes does not appreciably affect the analysis.

Fig. 9 shows $\Gamma_{++}/e^{-\Delta/T}$ computed for four system sizes over a range of temperatures. We see three temperature regimes for each system size: a low temperature regime where $\Gamma_{++}\sim Te^{-\Delta/T}$, a high temperature regime where $\Gamma_{++}\sim e^{-\Delta/T}$ and an intermediate temperature regime smoothly connecting these two forms of the temperature scaling.

Fig. 9 shows the low temperature model predictions of (26) as well as the Monte Carlo data. The regime of linear behavior of $\Gamma_{++}/e^{-\Delta/T}$ and corresponding agreement with the effective model at low temperatures allows us to identify this regime as the low temperature regime where the finite-size scaling of the relaxation time is determined by the scaling of topologically non-trivial random walks. Fig. 10 shows the finite-size scaling of $\Gamma_{++}$ in this low temperature regime, where we have performed a data collapse to remove the leading temperature dependence of (28). We find good agreement with the parameter free low temperature model (28) (with temperature dependent $P^{\Omega}_{\rm{2D}}(L)$ obtained by the resummation procedure in appendix A) which has an approximate $L^{2}/\ln L$ scaling. This non-trivial finite-size scaling is a key feature of this low temperature regime.

On increasing temperature, we can understand the transition out of this low temperature regime as follows. The separation of timescales breaks down as the ratio $L^{2}\cdot\gamma_{+}/\gamma_{0}$ grows larger. As the temperature increases, the decay rate is significantly affected by multi-pair processes which are not accounted for in (28). Interactions between multiple pairs of quasiparticles modify the annihilation probability distributions used in the low temperature model. Additionally, at higher temperatures, as the lifetime of quasiparticle pairs increases, the decay of $\Pi_{++}$ is sensitive to trivial error strings (i.e., single applications of $E^{e\dagger}$) acting across the edges shared with the winding operators. At higher temperatures, these trivial error strings dominate the decay rate.

At higher temperatures in Fig. 9, we see that $\Gamma_{++}\sim e^{-\Delta/T}$; this is the high temperature regime described by (27) where we expect a linear scaling in system size. In Fig. 11 we show a fit to the linear finite-size scaling for several different temperatures, where we have scaled the $\Gamma_{++}$ by the rate of formation of quasiparticle pairs, $\gamma_{+}$. This one parameter linear fit to the scaled data gives a single functional form for all system sizes and temperatures

where we find the constant $c_{T_{H}}=2.5\pm 0.1$. If only the lowest order process trivial anyon pairs contributed to the decay across both winding operators, we would have $c_{T_{H}}=2$. Obtaining a fit to the finite-size scaling with $c_{T_{H}}>2$ suggests that higher order processes are also providing significant contributions. The solid red line in Fig. 9 shows that this one parameter high temperature fit is in good agreement with the Monte Carlo data in the high temperature regime. We note that the linear finite-size scaling of $\Gamma_{++}$ is distinct from the $L^{2}/\ln{L}$ scaling in the low temperature regime (see Fig. 10).

Analysis of the results shown in Fig. 9 strongly suggests that we can identify two distinct regimes where the relaxation rate is dominated by distinct physical processes. We can therefore define a dynamical crossover temperature $T^{*}_{\rm{dyn}}$ which signifies the crossover between these regimes. We define $T^{*}_{\rm{dyn}}$ as the local maxima on Fig. 9 where the linear temperature scaling breaks down; this does not correspond to a maximum of $\Gamma_{++}$ itself, which is monotonically increasing as a function of temperature, since we have removed the temperature dependence of the Boltzmann factor by rescaling. Clearly $T^{*}_{\rm{dyn}}$ is a function of system size, since the low-temperature regime shrinks as $L$ increases; Fig. 12 displays the finite-size scaling of $T^{*}_{\rm{dyn}}$ as well as the equilibrium crossover temperature $T^{*}_{\rm{eq}}$ computed in [Castelnovo and Chamon, 2007] from the topological entanglement entropy. We find an inverse logarithmic scaling of $T^{*}_{\rm{dyn}}$ with system-size, in agreement with the scaling of $T^{*}_{\rm{eq}}$.