Topological ordering in the Majorana Toric Code

TL;DR

This paper investigates the topological phases of a Majorana toric code system, characterizing topological order using lattice gauge theory and Wilson-loop order parameters, and analyzing the effects of Cooper pair tunneling.

Contribution

It introduces a lattice gauge theory framework and uses a generalized Wilson-loop order parameter to identify topological phases in the Majorana toric code.

Findings

Toric code phase exists in Mott insulator and charge-2e superconductor phases.

Toric code phase disappears in charge-e superconductor phase.

Cooper pair tunneling stabilizes the topological gap of the toric code.

Abstract

At zero temperature, a two-dimensional lattice of Majorana zero modes on mesoscopic superconducting islands exhibits a topologically-ordered toric code phase. Recently, a Landau field theory was used to describe the different phases of the aforementioned system and the phase-transitions separating them. While the field theory provides details on the properties of the system close to the phase-transitions, signatures of topological ordering in the different phases have not been computed. This is the primary goal of the current work. We describe a lattice gauge theory of the Majorana toric code in terms of $\mathrm{U}(1)$ matter fields coupled to an emergent $\mathbb{Z}_2$ gauge field. Subsequently, we use a generalized Wilson-loop order-parameter, the equal-time Fredenhagen-Marcu order parameter, to characterize the topological ordering in the different phases. Our computation provides evidence of the toric code phase both in the Mott insulator and the charge-$2e$ superconductor phases, while showing that the toric code phase disappears in the charge-$e$ superconductor phase. In addition, we perturbatively analyze the influence of Cooper pair tunneling on the topological gap of the toric code in the limit of strong charging energy and show that the toric code phase is, in fact, stabilized by the Cooper pair tunneling. Our results are relevant for experimental realizations of the Majorana toric code.

Equations

$$H_C$$

$$H_J$$

$$H_M$$

$$Q_j|\text{Phys}⟩$$

$$Q_j=-{\gamma }_a^j{\gamma }_b^j{\gamma }_c^j{\gamma }_d^j{\text{e}}^{2\pi i{n}_j}.$$

$$i{\gamma }^j{\gamma }^k\to {\sigma }_{jk}^z\text{and}-{\gamma }_a^j{\gamma }_b^j{\gamma }_c^j{\gamma }_d^j\to \prod _{+}{\sigma }^x,$$

$$H_{MTC}$$

$$=4E_C\sum _j{n}_j^2-E_J\sum _{⟨j,k⟩}\cos \left({\varphi }_j-{\varphi }_k\right)$$

$$-E_M\sum _{⟨j,k⟩}{\sigma }_{jk}^z\cos \left(\frac{{\varphi }_j-{\varphi }_k}{2}\right).$$

$$Q_j={\text{e}}^{2\pi i{n}_j}\prod _{+}{\sigma }^x,$$

$$W(L)$$

$$W_{1/2}(L)=⟨G|e^{i{\varphi }_s/2}\prod _{l\in {\mathcal{C}}_{1/2}}{\sigma }_l^z{e}^{-i{\varphi }_{s^{\prime }}/2}|G⟩.$$

$$R(L)=\frac{W_{1/2}(L)}{\sqrt{W(L)}}.$$

$$\underset{L\to \infty }{\lim }R(L)=\{\begin{array}{ll}0,& {\mathbb{Z}}_2\text{deconfined,}\\ const.,& {\mathbb{Z}}_2\text{confined.}\end{array}$$

$$\prod _{+}{\sigma }^x|G⟩=|G⟩.$$

$$|G⟩=w\sum _{s\in \mathcal{S}}|n_i=0,s⟩=w\sum _{s\in \mathcal{S}}|s⟩.$$

$$R(L)={\text{e}}^{-\ln \left(4\frac{E_C}{E_M}\right)L}\stackrel{L\to \infty }{⟶}0,$$

$$H$$

$$\text{and}{Q}_i$$

$$H=-\Delta \sum _i\prod _{\square }{\sigma }^z-E_M\sum _{<i,j>}{\sigma }_{ij}^x.$$

$$|G⟩=w\sum _{s\in \mathcal{S}}|{\tau }_i^x=1,s⟩=w\sum _{s\in \mathcal{S}}|s⟩,$$

$$R(L)={\text{e}}^{-\ln \left(4\frac{\Delta }{E_M}\right)L}\stackrel{L\to \infty }{⟶}0,$$

$$|G⟩=w\sum _{m\in \mathcal{M}}|m⟩,$$

$$R(L)={\text{e}}^{-\frac{1}{16}{(\frac{\Delta }{E_M})}^2}\stackrel{L\to \infty }{⟶}\text{const.}$$

$$R(L)={\text{e}}^{-\frac{1}{16}{(\frac{{\Delta }^{\prime }}{E_M})}^2}\stackrel{L\to \infty }{⟶}\text{const.}$$

$$H_{\text{eff}}^{(4)}=\text{const.}-\frac{5}{16}\frac{E_M^4}{E_C^3}\sum _i\prod _{\square }{\sigma }^z.$$

$$H_{\text{eff}}^{(5)}=\text{const.}-\frac{61}{144}\frac{E_M^4{E}_J}{E_C^4}\sum _i\prod _{\square }{\sigma }^z.$$

$$E_M\to -E_M\text{and}{\varphi }_k\to {\varphi }_k+2\pi ,\forall k\in K,$$

$$H_{\text{eff}}^{(6)}=\text{const.}+H_{\text{eff}}^{(6,1)}+H_{\text{eff}}^{(6,2)}.$$

$$H_{\text{eff}}^{(6,1)}=-\left[\frac{259}{1728}\frac{E_M^6}{E_C^5}+\frac{24509861}{71124480}\frac{E_M^4{E}_J^2}{E_C^5}\right]\sum _i\prod _{\square }{\sigma }^z.$$

$$H_{\text{eff}}^{(6,2)}=-\frac{63}{256}\frac{E_{M}^{6}}{E_{C}^{5}}\sum\limits_{i}~{}\biggl{(}~{}\prod\limits_{\Box\Box}\sigma^{z}+\prod\limits_{\begin{subarray}{c}\text{$\Box$}\\ \text{$\Box$}\end{subarray}}\sigma^{z}~{}\biggr{)}\,.$$

$$i{\gamma }^k{\gamma }^l,\forall (k,l)\in \mathcal{L}\text{and}{P}_i,\forall i\in \mathcal{I}.$$

$${\sigma }_{kl}^z,\forall (k,l)\in \mathcal{L}$$

Full text

I Introduction

A promising candidate for realizing fault-tolerant quantum computation is Kitaev’s two-dimensional toric code Kitaev (1997, 2003, 2006); Fowler et al. (2012). The ground (code) space is topologically ordered Kitaev (1997); Levin and Wen (2005, 2006) and is four-fold degenerate. Thus, it can encode two logical qubits. These encoded logical qubits are robust to local perturbations since the degeneracy of the ground space depends only on the topology of the embedding space where the code is implemented. There are several approaches to realize the toric code. The first approach involves tessellating a two-dimensional plane with a regular lattice, with physical qubits on the links of the lattice. The Hamiltonian of the system is the toric code Hamiltonian which comprises vertex terms and plaquette terms. Each vertex term is a product of $\sigma_{x}$ operators of the qubits residing on the links incident at the vertex, while each plaquette term is a product of $\sigma_{z}$ operators of the qubits residing on the links around a plaquette Kitaev (2003). A different approach to realize the toric code is to design Hamiltonians of interacting many-body systems which, in the low-energy sector, give rise to the toric code Hamiltonian. This is the case of Kitaev’s honeycomb model, which describes SU(2) spins, interacting with alternating $\sigma_{x}\sigma_{x}$, $\sigma_{y}\sigma_{y}$ and $\sigma_{z}\sigma_{z}$ interactions around a plaquette Kitaev (2006). In this work, we concentrate on a different implementation of the latter approach involving a lattice of mesoscopic superconducting islands with Majorana zero modes (MZM-s) on each island Xu and Fu (2010); Terhal et al. (2012); Nussinov et al. (2012); Vijay et al. (2015); Landau et al. (2016); Karzig et al. (2017); Litinski et al. (2017); Burrello et al. (2013); Litinski and von Oppen (2018); Wille et al. (2019); Plugge et al. (2017); Hoffman et al. (2016). Note that any physical implementation of either approaches eventually also involves active measurements which evacuate the entropy due to thermal fluctuations since the topological ordering of the toric code disappears at any nonzero temperature. A toric code built out of MZM-based physical qubits holds the promise of better error correction properties compared to its superconducting transmon qubit-based counterpart. This is due to the potentially superior coherence properties of the MZM-based physical qubits Plugge et al. (2016); Landau et al. (2016); Vijay and Fu (2016); Vijay et al. (2015); Karzig et al. (2017).

The physical system of the Majorana toric code comprises mesoscopic superconducting islands which have finite charging energy ($E_{C}$). These nearest-neighboring islands are separated by tunnel junctions (see below for details), through which Cooper-pairs can tunnel coherently between islands (at rate $E_{J}$). Furthermore, due to the presence of the MZM-s, single electrons can also coherently tunnel between neighboring islands (at rate $E_{M}$). It was shown that for $E_{J}\ll E_{M}\ll E_{C}$, the system is in a topologically ordered toric code phase, while being a Mott-insulator Xu and Fu (2010); Landau et al. (2016). The presence of the MZM-s lead to an emergent $\mathbb{Z}_{2}$ gauge field and the topological ordering of the toric code. Increase of $E_{M}/E_{C}$, while keeping $E_{J}$ much smaller than $E_{M},E_{C}$, causes the system to undergo a 3D-XY type quantum phase transition to a charge-$e$ superconductor phase, the latter property arising from the dominant tunneling of single electrons. In this phase, the topological ordering of the toric code disappears Xu and Fu (2010). On the other hand, for $E_{M}\ll E_{C}\ll E_{J}$, the system was again shown to be in a topologically ordered toric code phase, while being a charge-$2e$ superconductor Terhal et al. (2012), the latter property arising from the dominant tunneling of Cooper-pairs. Increase of $E_{M}/E_{C}$, keeping $E_{J}$ significantly larger than $E_{M},E_{C}$, causes the system to undergo a 3D-Ising type quantum phase-transition to a charge-$e$ superconductor phase (without toric code ordering). These two different limiting cases were shown to be connected by tricritical points and first order transitions Roy et al. (2017). Using a Landau-Ginzburg-Wilson field theory, the phase-diagram of the model for general $E_{C},E_{M},E_{J}$ was predicted Roy et al. (2017) and experimentally accessible charge transport characteristics were also computed Roy and Hassler (2018). However, the above field theory calculations focused on the charge characteristics of the system and provided only indirect evidence of the existence of the topological order.

Our goal, in this work, is to fill this void and compute non-local order parameters which provide direct evidence of the presence/absence of topological ordering in the different phases. In contrast to the coarse-grained field theory, we analyze the system at a microsopic lattice level, by transforming the problem to a lattice gauge theory one with $\mathrm{U}(1)$ matter fields and a $\mathbb{Z}_{2}$ gauge field. In this mapping the emergent $\mathbb{Z}_{2}$ gauge field of the toric code arises naturally. The coupling between gauge and matter fields causes additional complications in the detection of topological ordering. For a gauge theory which has local symmetry, in absence of matter fields, the Wilson loop is the relevant non-local order parameter Elitzur (1975); Kogut (1979). The characteristic area vs perimeter law decay of the latter is able to distinguish between the different phases of the theory. This should be contrasted to a theory with global symmetry, where the phases are distinguished by a local order parameter. However, the presence of matter fields in a gauge theory changes the situation since the fluctuations of the matter field can screen the interaction between the charges of the gauge field. This can cause the Wilson loop to decay with perimeter law in all the phases Fradkin and Shenker (1979). As will be shown in this work, this is indeed the case for the Majorana toric code. For such theories, the concept of Wilson loops has to be generalized and the nonlocal order parameters capable of distinguishing between the different phases of these theories are the Fredenhagen-Marcu operators Fredenhagen and Marcu (1986); Gregor et al. (2011). We compute the latter operators in the different phases using time-independent perturbation theory. We show that the topological ordering exists in the Mott-insulator and charge-$2e$ superconductor phases of the system, while it disappears in the charge-$e$ superconductor phase.

Recent experimental endeavours in mesoscopic superconducting systems have shown a lot of promise towards detection of the MZM-s Mourik et al. (2012); Albrecht et al. (2016); Nichele et al. (2017); Lutchyn et al. (2018); Zhang et al. (2018); Wang et al. (2018). Motivated by these developments, we also compute the topological gap using perturbation theory for small, but finite, $E_{J}$ and $E_{M}\ll E_{C}$. This computation is relevant for near-term experimental implementation of the Majorana toric code since the experiments are most likely to be done for finite Josephson tunneling rate.

The article is organized as follows. In section II, we describe the microscopic Hamiltonian of a two-dimensional square lattice of MZM-s on mesoscopic islands and map it to a lattice gauge theory of $\mathrm{U}(1)$ matter and $\mathbb{Z}_{2}$ gauge fields. In Sec. III, we introduce and compute the equal-time Fredenhagen-Marcu operator in the different phases of the system. In Sec. IV, we determine the influence of inter-island Cooper pair hopping on the toric code gap in the limit of dominant charging energy. In Sec. V, a concluding summary is provided. Details of calculations and a short review of lattice gauge order parameters are given in the appendices.

II The model

The microscopic model consists of a two dimensional lattice of mesoscopic superconducting islands, each carrying four MZM-s occuring as edge-modes of Kitaev wires Kitaev (2001). The nearest-neighbor islands are separated by tunnel junctions (see Fig. 1). The Hamiltonian of the system is given by $H=H_{C}+H_{J}+H_{M}$, where

[TABLE]

The first term, $H_{C}$, denotes the charging energy associated with each mesoscopic island. Here $E_{C}$ is the overall charging energy scale $e^{2}/2C$, where $C$ is the self-capacitance of each island. In contrast to the Cooper-pair box Koch et al. (2007); Shnirman et al. (1997), here $n_{i}$ denotes excess number of fermions on each island and thus, can take both integer or half-integer values. The four MZM-s on each island are denoted by Hermitian operators $\gamma^{j}_{\alpha}$, $\alpha=a,b,c,d$, satisfying anticommutation relations: $\{\gamma^{j}_{\alpha},\gamma^{k}_{\beta}\}=2\delta_{\alpha\beta}\delta_{jk}$. Even though they are not directly present in $H_{C}$, the MZM-s impose a gauge (parity) constraint satisfied by the physical state, $\ket{\text{Phys}}$, of the system, given by Fu (2010); Roy and Hassler (2018); van Heck et al. (2012)

[TABLE]

The second term, $H_{J}$, describes the coherent tunneling of Cooper pairs between neighboring islands, $E_{J}$ being the Josephson tunneling rate. Finally, the third term, $H_{M}$, describes the Majorana-assisted single-electron tunneling, $E_{M}$ being the relevant tunneling rate. The factor of $1/2$ in the argument of the cosine indicates that $1/2$ of the charge of the Cooper-pair is being transferred, while the fermionic operators keep track of the change in the fermion number parity Fu (2010); van Heck et al. (2011). To avoid clutter, we have dropped the subscripts of the MZM operators in Eq. (1c). Throughout this work, we consider the case of zero off-set charge.

In order to analyze the topological ordering in the system, we map the Hamiltonian in Eq. (1) to that of interacting spins coupled to quantum rotors. This is done by the transformation

[TABLE]

using the bond algebraic approach (Cobanera et al., 2011; Nussinov et al., 2012) (see Appendix A for details). This can be viewed as a Jordan-Wigner transformation followed by a duality transform and maps the product of the MZM-s located at the endpoints of a link to a spin placed centrally on it. This transformation keeps the number of degrees of freedom (DOF) invariant and the resulting Hamiltonian is given by (see Fig. 2)

[TABLE]

The gauge constraint is transformed to

[TABLE]

where in the last equation, $\prod_{+}\sigma^{x}$ indicates the product of the Pauli-X operators of the spins residing on the links incident at the $j^{th}$ island.

Equation (5) is the discretized Gauss’ law for the system Kogut (1979), but unlike conventional electrodynamics, a $\mathrm{U}(1)$ matter field ($e^{i\phi_{j}}$), residing on the vertices, is coupled through its parity to a $\mathbb{Z}_{2}$ gauge field ($\sigma_{jk}^{z}$), residing on the links. The gauge fields do not have dynamics of their own since only $\sigma_{jk}^{z}$ operators occur in the Hamiltonian [Eq. (4)]. However, they acquire dynamics through the gauge constraint which couples them to the parity of the $\mathrm{U}(1)$ matter fields.

The phase-diagram (see Fig. 4) of the system is rich and has been analyzed Terhal et al. (2012); Roy et al. (2017); Xu and Fu (2010); Roy and Hassler (2018). For small $E_{J}\ll E_{M}\ll E_{C}$, the system is in a toric code phase, while being a Mott insulator. Increasing $E_{M}/E_{C}$ keeping $E_{J}\ll E_{M},E_{C}$, causes the system to undergo a 3D-XY quantum phase transition to a charge-$e$ superconductor phase. For $E_{J}\gg E_{M}\gg E_{C}$, the system is again in a toric code phase, while being a conventional superconductor. Increasing $E_{M}/E_{C}$, while keeping $E_{J}\gg E_{M},E_{C}$, causes the system to undergo a 3D-Ising type quantum phase transition to a charge-$e$ superconductor state. These different phase-transitions are connected by a couple of tricritical points and first-order transitions. Finally, a 3D-XY type quantum phase transition separates the Mott-insulator and the conventional superconductor phases.

While charge signatures of the different phases have been analyzed in the earlier works, the existence of the toric code has been indirectly inferred. In the next section, we compute non-local order parameters to provide direct evidence of the presence/absence of toric code ordering.

III Analysis of topological ordering using The Fredenhagen-Marcu operator

For models with global symmetry, the phases are distinguished by a local order parameter (the celebrated symmetry-breaking paradigm of Landau). In contrast, for models with gauge symmetry, which cannot be broken Elitzur (1975), the different phases can be differentiated by non-local order parameters Nussinov and Ortiz (2009). For a pure gauge theory, like the Ising gauge theory, the Wilson loops have their characteristic area (perimeter) law decay in the confined (deconfined) phases (see Ref. Kogut, 1979 and references therein). However, the presence of dynamical matter fields, alongside gauge fields, introduces additional complications. The matter fields potentially screen the fluctuations of the gauge fields and can cause the Wilson loop of the gauge fields to decay with a perimeter law in both the confined and the deconfined phases (throughout we refer to confinemenet/deconfinement of the external charges of the gauge field) of the system, e.g., the Ising gauge theory in the presence of Ising matter fields Fradkin and Shenker (1979). The Majorana toric code, as was shown in the previous section, can be viewed as such a theory, where dynamical $\mathrm{U}(1)$ matter fields interact with $\mathbb{Z}_{2}$ gauge fields.

For these theories, the notion of the Wilson loop has to be generalized to distinguish between the different phases of the model and the relevant non-local order parameter is the Fredenhagen-Marcu (FM) operator Fredenhagen and Marcu (1986), proposed initially for models of particle physics and has been used for condensed matter systems Gregor et al. (2011). In contrast to the Wilson loop which has the same behavior in both confined and deconfined phases, the FM operator has the desired feature of an order parameter: it is zero (nonzero) in the confined (deconfined) phases of the system. There are three different operator formulations depending on the space-time orientation of the non-local operators Gregor et al. (2011) and we use the equal-time formulation (a short review of the Wilson loop and different FM order parameters is given in Appendix B).

The Wilson loop for the system is given by

[TABLE]

where $\sigma^{z}_{l}$ are the Pauli-Z operators of the spins residing on the contour $\mathcal{C}$, the latter being a square of side length $L$ [cf. Fig. 3(a)], and $|G\rangle$ is (one of) the ground state(s) of the system. In order to construct the FM operator, we need also the modified half Wilson loop. It comprises a product of a string of gauge field operators along the contour $\mathcal{C}_{1/2}$, terminated by two matter field operators at sites $s$ and $s^{\prime}$, given by

[TABLE]

The equal-time FM operator, $R(L)$, is given by

[TABLE]

In the limit of an infinitely large loop, one can show that (see Appendix  B for details)

[TABLE]

The construction of the FM operator can be understood as follows. In order to correctly diagnose the signatures of the confining/deconfining phases, a ‘correct string tension’, which determines the potential between the gauge charges, needs to be computed. To that end, consider half a Wilson loop [see Fig. 3(a)]. Just by itself, this quantity is not gauge-invariant, but can be made so by gluing matter fields at the end-points. In the confining phase, this half-Wilson loop with matter fields at the ends, $W_{1/2}(L)$, acquires a finite expectation value, whereas in the deconfining phase, it vanishes; see Appendix B. The division by a square-root of the Wilson loop is done to distill the dependence of the diagnostic on the coupling constant from the size of the loop.

Intuitively, both the Wilson loop and the FM operator can be best interpreted in the space-time formulation [see Fig. 3(b)]. In this formulation, loops extend in one spatial and one temporal direction. In absence of matter charges, the Wilson loop is a quantity purely involving gauge DOF and is related to the potential between external, static, charges of the gauge fields. In the presence of dynamical matter charges which can potentially screen the interaction between the gauge field charges, the Wilson loop is no longer capable of correctly diagnosing the potential between the gauge charges. This is remedied by the FM operator which is defined in terms of the dynamical matter charges to compensate for the screening effect. More on the FM operator can be found in Appendix B and Ref. Gregor et al., 2011.

The Majorana toric code describes interacting $\mathrm{U}(1)$ matter fields and $\mathbb{Z}_{2}$ gauge fields in two spatial and one (Euclidean) time direction. The toric code phase corresponds to the system being in the $\mathbb{Z}_{2}$ deconfined phase. Below, we compute $R(L)$ in the different phases of the Majorana toric code using time-independent perturbation theory. Since we evaluate the FM operator expectation value at equal times, it is sufficient to consider the expectation value in the ground state of the system.

Before providing the details of the computation, we provide a summary of the findings in Fig. 4. The Mott-insulator [phase (a)] and the charge-2$e$ superconductor [phase (b)] are in the toric code phase, indicated by $R(L)\rightarrow 0$, while the charge-$e$ superconductor [phase (c)], with $R(L)\rightarrow$ constant, does not exhibit toric code ordering. This result also indicates that there are no additional phase-transitions as $E_{J}$ is varied between [math] and $\infty$ for large $E_{M}$. Thus, the results of our computations of the FM operator support the earlier field theory computations.

III.1 Dominant charging energy

First, we consider the parameter regime where the Josephson as well as the single-electron tunneling rate is small compared to the charging energy: $E_{J},E_{M}\ll E_{C}$. It is evident that the charge (excess number of fermions) on each island is a good quantum number and the system is a Mott insulator Roy et al. (2017). The unperturbed Hamiltonian is $H_{C}$, given in Eq. (4), with the unperturbed ground state being $|G\rangle=|n_{i}=0\rangle$, $\forall i$. As a consequence, the gauge constraint reduces to

[TABLE]

Imposing this constraint allows eight configurations of gauge fields at each node, as shown in Fig. 5, where out of the four spins incident at each vertex, an even number of them are in $\sigma^{x}=+1$ (blue circles) and $\sigma^{x}=-1$ (red circles).

If we imagine coloring the links of the lattice with red whenever the spins are in the $\sigma^{x}=-1$ eigenstate, then all possible closed red loop configurations of the lattice are allowed and the ground state is a loop condensate Kitaev (2003); Fradkin (2013). The degeneracy of the ground state is four-fold, due to the presence of the four non-contractible loops on a torus. While calculating using perturbation theory 111For the perturbation theory calculation, it is sufficient to consider one of the four degenerate ground states since their splitting is exponentially suppressed in the system size., we impose the constraint in Eq. (10) directly on the unperturbed ground state and obtain equal weight ($w$) superpositions of all possible closed loop configurations $\mathcal{S}$

[TABLE]

Including the single and double charge hopping from Eq. (4) as perturbations $V=H_{M}+H_{J}$, time independent calculations yield

[TABLE]

indicating the presence of a $\mathbb{Z}_{2}$ deconfined (toric code) phase for dominant charging energy (see Appendix C.1 for details of the computation).

III.2 Dominant Josephson tunneling rate

Next, we analyze the case when the Josephson tunneling rate is the strongest: $E_{J}\gg E_{M},E_{C}$, when the system can be mapped to an effective Ising gauge theory, as explained below. In this limit, the phase differences on the islands are pinned to multiples of $2\pi$ up to a gauge choice. In this reduced subspace, the $\mathrm{U}(1)$ DOF behaves like an effective $\mathbb{Z}_{2}$ DOF. The operator $\text{e}^{\pm i\phi_{i}/2}$ acts like an effective Pauli-Z operator: $\tau_{i}^{z}$ whose eigenvalues can be $\pm 1$ depending on whether the $\mathrm{U}(1)$ DOF are pinned to even/odd multiples of $2\pi$. The charging energy term induces quantum fluctuations of the phase, leading to constant energy shifts and $2\pi$-phase slips on the nodes. Thus, the exponent of the number operator $n_{i}$ acts like a Pauli-X operator: $\text{e}^{2\pi in_{i}}~{}\rightarrow~{}\tau^{x}_{i}$, where $\commutator{\tau_{i}^{z}}{\tau_{j}^{x}}=\commutator{\text{e}^{\pm i\phi_{i}/2}}{\text{e}^{2\pi in_{j}}}$. Disregarding the constant contribution, the charging energy term can effectively be incorporated into the $\mathbb{Z}_{2}$ formulation as $H_{C,i}~{}\rightarrow~{}-\Delta\tau_{i}^{x}$ leading to the effective Hamiltonian

[TABLE]

To determine the dominant coupling dependence of $\Delta$, assume the slips to occur separately on each island. This approximation is analogous to the one dimensional case treated in Ref. van Heck et al., 2014. The coupling of independent slips is then calculated as the tunnelling amplitude in the cosine potential using semi-classical WKB approximation Landau and Lifshitz (1958); Terhal et al. (2012), leading to $\ln{\Delta}\propto-\sqrt{E_{J}/E_{C}}$. The obtained Hamiltonian is that of the 2D quantum $\mathbb{Z}_{2}$ gauge theory. This is seen by choosing the London gauge, where the matter fields are unity: $\tau^{z}_{i}=1$ and using the gauge constraint to replace $\tau^{x}_{i}=\prod_{+}\sigma^{x}$. Switching to the dual lattice and rotating the basis yields

[TABLE]

This directly verifies the limit considered in Ref. Terhal et al., 2012, where the toric code phase was predicted to exist for $\Delta\gg E_{M}$. Increasing $E_{M}/\Delta$ caused the system to undergo a 3D-Ising type phase-transition to a non-toric code phase. The argument of Ref. Terhal et al., 2012 is based on working in a dual picture and analyzing the phase-transition of the quantum Ising model. Below, we provide a more direct proof of the toric code ordering by computing the expectation value of the FM operator in the ground state of the gauge-invariant Hamiltonian given in Eq. (13).

III.2.1 Large charging energy: $E_{C}\gg E_{M}$

Since in this limit, $\Delta\gg E_{M}$, from Eq. (13), the matter fields are all pinned: $\tau^{x}_{i}=1$. From Eq. (14), only those configurations of the gauge field are allowed which satisfy the gauge constraint. Those are closed loops where loops are formed by the gauge field in $\sigma^{x}=-1$ state, exactly as for the case $E_{C}\gg E_{M},E_{J}$ analyzed in Sec. III.1. Thus, the ground state is given by

[TABLE]

which is again a four-fold degenerate loop-condensate. Time independent perturbation theory then yields

[TABLE]

again indicating a toric code phase.

III.2.2 Large single-electron tunneling rate: $E_{M}\gg E_{C}$

In this limit, $E_{M}\gg\Delta$ and from the gauge constraint in Eq. (14), the valid vertex configurations displayed in Fig. 6 allow for open (red) loops of $\sigma^{x}=-1$ terminated by flipped parity ($\tau_{i}^{x}=-1$). The unperturbed ground state is

[TABLE]

where $\mathcal{M}$ denotes the set of all open loop configurations. Note that closed loops are treated as open loops with coinciding start- and endpoint such that those are also contained in $\mathcal{M}$. In contrast to the closed loop condensate, where only an even number of non-contractible loops are generated, Eq. (18) contains also odd numbers of non-contractible loops in the equal weight superposition leading to a non-degenerate ground state. Computing the expectation value of the FM operator in the ground state, we find (see Appendix C.2 for details)

[TABLE]

This indicates that the system is in a $\mathbb{Z}_{2}$ confined phase [see Eq. (9)] and thus, does not have the toric code ordering.

III.3 Dominant single electron tunneling rate

We consider the limit of $E_{M}\gg E_{C},E_{J}$ and we restrict ourselves to the case of $E_{J}=0$ since the limit of large $E_{J}$ was already analyzed in the previous section. For $E_{M}\rightarrow\infty$ in Eq. (4), we find that $\phi_{i}-\phi_{j}=0$ or $2\pi$ depending on whether $\sigma_{ij}^{z}=1$ or $-1$. Assuming a fixed configuration $\{\sigma_{ij}^{z}\}$, the phases are again pinned to $\phi_{i}=2\pi m_{i}$, where $m\in\mathbb{Z}$ . Thus, we can apply the same mapping as in the previous seciton to map the $\mathrm{U}(1)$ DOF to the Ising variable $\tau_{i}^{z}$. Furthermore, the charging term induces quantum fluctuations such that $H_{C,i}\rightarrow\Delta^{\prime}\tau_{i}^{x}$. Assuming independent slips, $\Delta^{\prime 2}$ can be estimated by the tunnelling event in the $4\pi$ periodic cosine potential using the WKB method Landau and Lifshitz (1958). This yields $\ln{\Delta^{\prime}}\propto\-\sqrt{E_{M}/E_{C}}$. Thus, the Hamiltonian takes the same form as in Eq. (13) only differing by an adjusted tunnelling rate $\Delta^{\prime}$. It has the same non-degenerate unperturbed ground state resulting in the Fredenhagen-Marcu operator to be evaluated as

[TABLE]

Thus, the system is again in a $\mathbb{Z}_{2}$ confined phase and there is no toric code ordering. The summary of the our findings is given in Fig. 4.

IV Stability of the Toric Code in the dominant charging energy regime

In this section, we compute the topological gap of the toric code phase perturbatively in $E_{J}/E_{C}$ and $E_{M}/E_{C}$. We consider the Hamiltonian [Eq. (4)] of the system in the Mott insulator phase (a), where $E_{C}\gg E_{J},E_{M}$ and determine the effective Hamiltonian up to sixth order in perturbation theory, for details see Appendix D. The goal of this computation is to provide a quantitative prediction of the effect of Cooper pair tunneling on the toric code gap. The different virtual processes that contribute to the effective Hamiltonian are depicted in Fig. 7. The number of lines indicates the number of single charges that are transferred and the color green denotes the appearance of $\sigma^{z}$ for the link. Every valid contribution to $H_{\text{eff}}$ is given by a process that starts and ends in a charge less state [see the projectors $P_{-}$ in Eq. (68)].

The lowest order contribution to the effective Hamiltonian that is not a constant, is obtained at fourth order in the perturbation theory. This process is depicted in Fig. 7. This is due to the process which transfers a single charge around the smallest loop on the lattice, a plaquette. The $\mathbb{Z}_{2}$ (gauge) DOF are carried around by the $\mathrm{U}(1)$ charge transfer operators $e^{\pm i\phi/2}$ and lead to a non-constant contribution to the Hamiltonian [see Fig. 7]. The respective fourth order term in operator formulation is given by Landau et al. (2016)

[TABLE]

In combination with the gauge constraint [Eq. (5)], this gives the toric code Hamiltonian.

To determine the leading order contribution in $E_{J}$, we have to go to fifth order in perturbation and find the process shown in Fig. 7. A single charge effectively flows around a plaquette. However, in contrast to Fig. 7, a single-electron tunnels opposite to the direction of the overall charge-flow around the plaquette. This opposite movement of the single electron is compensated by a transfer of Cooper pair along the direction of the overall charge-flow around the plaquette. Since the gauge operators appear independent of the flow direction [cf. Eq. (4)] it yields a plaquette term. The contribution to the effective Hamiltonian is given by

[TABLE]

Thus, for small $E_{J}$, the hopping of Cooper pairs increases the gap of the toric code and stabilizes the toric code space. In contrast to the Cooper pair box Hamiltonian $H_{C}+H_{J}$, where the sign of $E_{J}$ leaves the spectrum invariant, it becomes important due to the introduction of $H_{M}$. Instead, Eq. (4) is only left invariant by the transformation

[TABLE]

where $K$ is a set of indices containing next nearest neighbouring islands, one sub lattice in a bipartite lattice. The phase diagram (cf. Fig. 4) is thus not expected to be symmetric along the $E_{J}=0$ axis, while the transition lines approach the $E_{M}=0$ at a right angle.

At sixth order, the non-constant contributions can be separated in two terms:

[TABLE]

First, there are next to leading order plaquette terms, which we compute to be

[TABLE]

Second, we find an additional process depicted in Fig. 7. A single charge can flow around the perimeter of two adjacent plaquettes, which is the second smallest cycle in the square lattice. The gauge fields on the perimeter get thereby flipped. The sum of the relevant diagrams for this contribution is calculated to be

[TABLE]

The $\sigma^{z}$ on the shared link can be added for clarity as they square to identity. This term introduces a nearest neighbour interaction between adjacent plaquettes. The interaction commutes with the single plaquette term and results in an additional stabilization of the ground state space Terhal et al. (2012).

V Conclusion

To summarize, we have analyzed the signatures of topological ordering in the different phases of the Majorana toric code. To that end, we mapped the system onto a model of interacting $\mathrm{U}(1)$ matter fields and $\mathbb{Z}_{2}$ gauge fields. We computed a non-local order parameter, the equal-time Fredenhagen-Marcu order parameter to capture the topological ordering in the system. Our calculations confirm that the Mott-insulator and the charge-2$e$ superconductor phases of the system show an emergent toric code, while the charge-$e$ superconductor phase does not exhibit toric code ordering. Our results are compatible with the earlier field theory predictions in the different parts of the phase-diagram. Furthermore, we computed the topological gap of the toric code in the presence of finite Cooper pair tunneling and showed that a small amount of Cooper pair tunneling actually stabilizes the toric code phase since the topological gap increases with a small amount of Cooper pair tunneling.

VI Acknowledgements

A.Z. and F.H. acknowledge funding by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany’s Excellence Strategy – Cluster of Excellence Matter and Light for Quantum Computing (ML4Q) EXC 2004/1 – 390534769. A.R. acknowledges the support of the Alexander von Humboldt foundation.

Appendix A From Majoranas to spins

In this section, we construct the duality transformation given in Eq. (3) based on the bond algebraic approach from Refs. Nussinov et al., 2012; Cobanera et al., 2011. First, we define $\mathcal{L}$ the set of all nearest neighbour bonds, $\mathcal{I}$ the set of all nodes and the operator $P_{i}=-\gamma_{a}^{i}\gamma_{b}^{i}\gamma_{c}^{i}\gamma_{d}^{i}$ to compactify the notation. The set of bonds generating the bond algebra is then given by

[TABLE]

Making use of the Clifford algebra one can show that the bond algebra $\mathcal{A}$ contains the intensive relations, i.e., independent of the lattice size,

• I)

$(i\gamma^{k}\gamma^{l})^{2}=P_{j}^{2}=\bm{1}\,,\qquad\forall(k,l)\in\mathcal{L}\,,\forall j\in\mathcal{I}\,,$

• II)

$\anticommutator{P_{k}}{i\gamma^{k}\gamma^{l}}=\anticommutator{P_{l}}{i\gamma^{k}\gamma^{l}}=0\,,\qquad\forall(k,l)\in\mathcal{L}\,,$

• III)

all remaining bond combinations commute.

For periodic boundary conditions, there remains one extensive, i.e., lattice size dependent, relation

• IV)

$\prod\limits_{j\in\mathcal{I}}P_{j}=\alpha\prod\limits_{(k,l)\in\mathcal{L}}i\gamma^{k}\gamma^{l}$  ,

where $\alpha=\pm 1$ accounts for the different orderings of the Majorana operators on both sides of the equation. This last equation considers the global $\mathbb{Z}_{2}$ symmetry of the Majoranas. We want to reproduce the same bond algebra $\mathcal{A}$ with alternative bonds, where there is a spin-$1/2$ degree of freedom placed on the links between two islands (cf. Fig. 2). This DOF is represented by the regular Pauli operators $\sigma_{ij}^{z},\sigma_{ij}^{y},\sigma_{ij}^{x}$ obeying the Pauli algebra. The sizes of the Hilbert spaces agree as there are two Majorana modes per island which can be occupied or not in case of the original bonds and two links per island containing a spin-$1/2$ DOF which can be up or down in the alternative bond description. The new set of bonds can be defined as

[TABLE]

so that $\tilde{P}_{i}$ is given by the product of $\sigma^{x}$ of all four links emanating from island $i$. With the Pauli algebra it can be verified that the alternative bonds satisfy the intensive relations

• I)

$(\sigma_{kl}^{z})^{2}=\tilde{P}_{j}^{2}=\bm{1}\,,\qquad\forall(k,l)\in\mathcal{L}\,,\forall j\in\mathcal{I}\,,$

• II)

$\anticommutator{\tilde{P}_{k}}{\sigma_{kl}^{z}}=\anticommutator{\tilde{P}_{l}}{\sigma_{kl}^{z}}=0\,,\qquad\forall(k,l)\in\mathcal{L}\,,$

• III)

all remaining bond combinations commute.

In order to fulfill the extensive relation, we choose one arbitrary, but fixed island and set $\tilde{P}_{0}=\prod\nolimits_{+_{0}}\sigma^{x}\cdot(\alpha\prod_{(k,l)\in\mathcal{L}}\sigma_{kl}^{z})$. With this choice, the intensive relations still hold, as the additional prefactor commutes with all remaining bonds and we obtain

• IV)

$\prod\limits_{j\in\mathcal{I}}\tilde{P}_{j}=\alpha\prod\limits_{(k,l)\in\mathcal{L}}\sigma_{kl}^{z}\,,$

since all $\sigma^{x}$ square to identity. Now that the Hilbert space sizes are equal and both the intensive as well as the extensive relations are satisfied, we can write down the new model description as in Eq. (4). Both Hamiltonians are linked by a unitary transformation, see Ref. Nussinov et al., 2012.

Appendix B Lattice gauge order parameters

In this section, we give a short review on the Wilson loop as well as an in-depth discussion of the Fredenhagen-Marcu operator to make the article self-contained.

B.1 The Wilson loop operator

The Wilson loop was introduced by Wegner Wegner (1971) for the isotropic Ising gauge model, where it distinguishes the area from perimeter law Kogut (1979). For anisotropic models, e.g. with continuous (imaginary) time direction as in Eq. (4), distinct operator formulations for the Wilson loop are possible (cf. Fig. 8).

The space-space operator formulation of the Wilson loop is given by

[TABLE]

To discuss the contour $\mathcal{C}$, we have to look at Figure 8(a). The contour is a square of side length $L$, thus the $L$ dependence on the left hand side of Eq. (29) stems from this implicit dependence in $\mathcal{C}$. The space-space formulation yields the ground state expectation of a square loop of side length $L$. It is a gauge invariant measure of the correlation between gauge fields, that are a distance $L$ apart.

To construct the space-time formulation, we define the imaginary time dependent gauge spin operators per link $l$ as Gregor et al. (2011)

[TABLE]

Furthermore, we apply the temporal gauge $\sigma_{l^{\prime}}^{z}=1$ for all links $l^{\prime}$ in the imaginary time direction. The length of the loop in imaginary time is given by $L_{\tau}=v\mathcal{T}$, where we assume the velocity $v$ to be unity. This results in a space-time operator formulation of

[TABLE]

The strings are defined as $S(L,\tau)=\prod\nolimits_{l\in S_{L}}\sigma_{l}^{z}(\tau)$, where $S_{L}$ is the curve of links depicted in Fig. 8(b). Those are the remaining spatial lines of the loop at different time slices. For convenience, we choose the strings to lie at times $\tau=-\mathcal{T}$ and $\tau=0$. This operator is connected to the potential $V$ between two external static, infinitely heavy, charges of the gauge field placed a distance $L$ apart Kogut (1979) with

[TABLE]

The latter relation leads to the labelling of phases in the Ising gauge model as

[TABLE]

where $A$ $(P)$ denotes the area (perimeter) of the Wilson loop and $\alpha,\beta$ are model dependent prefactors. In the presence of dynamical matter this diagnostic breaks down, if the matter fluctuations introduce screening and thus perimeter law everywhere Fradkin and Shenker (1979).

B.2 The Fredenhagen-Marcu operator

The interpretation of the equal-time Fredenhagen-Marcu operator defined in Eqs. (6)–(8) is analyzed using the lattice gauge theory of Ising gauge coupled to Ising matter, where we follow the treatment in Ref. Gregor et al., 2011. We consider a system of spins $\tau/\sigma$ on the nodes/links of a two dimensional square lattice. The respective Hamiltonian is given by

[TABLE]

We refer to nearest neighbour hopping of matter DOF [proportional to $\lambda$ in Eq. (34)] as edges and aligned gauge DOF around plaquettes [proportional to $\Delta$ in Eq. (34)] as surfaces. Since the perimeter of a square of side length $L$ is considerably smaller than its area for large $L$, the system will generate the (half) Wilson loop both by mainly constructing its perimeter with edges. More generally, we say that large loops will always be dominantly covered by edges as opposed to surfaces, as long as $\lambda\gtrapprox\Delta$ [cf. Fig. 9(a)]. Therefore, the numerator and denominator in Eq. (8) have the same scaling behaviour, except when surfaces are much cheaper than edges. This is the case, when the matter field is considerably heavy, while the gauge field is ordered $\Delta\gg\lambda$. Then the closed loop gets filled with surfaces. The open loop is equivalently covered with surfaces, but the remaining ($L$ dependent) line has to be covered by expensive edges [cf. Fig. 9(b)]. Thus, the numerator decays much faster than the denominator. Central to these scaling behaviours is the necessity of matter flow (edges) in the generation of the half Wilson loop. In conclusion, we have that

[TABLE]

Either the numerator decays faster or the scaling is equal up to a constant. By taking the $\lambda\rightarrow 0$ limit, we can rewrite Eq. (34) as

[TABLE]

where we replaced $\tau^{x}_{i}=\prod\nolimits_{+}\sigma^{x}$ via the gauge constraint. If also $\Delta\rightarrow\infty$, we find an emergent toric code. Using this we can infer

[TABLE]

This connection can be made more explicit by considering the original Fredenhagen-Marcu order parameter, where one spatial direction is exchanged for the imaginary time direction [cf. Fig. 3(b)].

We construct the operator formulation along the lines of Ref. Gregor et al., 2011. We assume that the two matter spins $\tau_{s}$ and $\tau_{s^{\prime}}$ lie at the same time slice, which we set to $\tau=0$. Choosing the temporal gauge and performing the imaginary time propagation analogously to Eq. (30), we find that

[TABLE]

The gauge flux generated from a charge at $s$ and anti-charge at $s^{\prime}$ gets introduced at time $-\mathcal{T}/2$ and the respective charges are created at time zero. Using a unit velocity $v$, we choose $v\mathcal{T}=L_{\tau}=L$ such that the loop only has $L$ dependence. For the Wilson loop we have

[TABLE]

We explicitly add and remove the charges for convenience. We define a state where we added a charge and anti-charge, that are a distance $L$ apart, to the ground state as

[TABLE]

If we assume the distance to go to infinity, we can focus on the charge and neglect the anti-charge at infinity. Thus, the state $\ket{S(L\rightarrow\infty)}$ can be interpreted as a free charge state. Applying the reformulation to the (half) Wilson loop gives

[TABLE]

The Fredenhagen-Marcu operator is rewritten to

[TABLE]

i.e., the overlap between the ground state and the normalized charge anti-charge state. Taking the $L\rightarrow\infty$ limit indicates a phase containing free (deconfined) charge states $\ket{S(L\rightarrow\infty)}$ only if $R(L)$ decays to zero, since the state is orthogonal to the ground state. If in turn there is a finite overlap, the state $\ket{S(L)}$ would decay into the ground state signalling only confined charges Gregor et al. (2011); Fredenhagen and Marcu (1986). Therefore, we write

[TABLE]

This diagnostic was first proposed in Ref. Fredenhagen and Marcu, 1986 to test whether a theory contains free or confined quarks. It generally detects a confinement-deconfinement transition in a lattice gauge theory containing dynamical matter charges.

Appendix C Evaluation of the Fredenhagen-Marcu operator

In this section, we discuss the perturbative evaluation of the Fredenhagen-Marcu operator in more detail.

C.1 Dominant charge energy

For the regime of dominant charging energy, $E_{C}\gg E_{J},E_{M}$, we found the unperturbed ground state in Eq. (11). For clarity, we expand the cosine functions and set $d_{j}=\text{e}^{-i\phi_{j}/2}$ to write the perturbation as

[TABLE]

The ground state to first order follows to

[TABLE]

We define $N^{\prime}$ to be the normalization factor. For the second term, the Majorana term, the two single charge creation/annihilation operators each lead to an $n_{i}^{2}=1/4$ on the respective island and thus result in an energy change of $2\cdot E_{C}$. This is the charging energy of two single charges. For the third term, the Josephson term, we analogously obtain $2\cdot 4E_{C}$, the charging energy of two Cooper pairs. With this perturbative ground state one can calculate the expectation value of the (half) Wilson loop, both to first order. In the following, we denote the state with respect to which we calculate the expectation value of the Wilson loop explicitly by redefining $W(L)=\expectationvalue{W(L)}{G}$. As the mixed terms vanish, this yields

[TABLE]

The first term is the result for the unperturbed ground state. For the remaining terms it is important to note that the creation/annihilation operators $d_{j}^{\dagger}/d_{j}$ commute with the Pauli operators and that the perturbations are Hermitian. In the second term, the Majorana term, one finds that it only yields a finite contribution if the charge transfer occurs on the same link. This happens for $2N$ terms of the sum, where $N$ is the number of nodes in the system. For the third term, the Josephson term, we analogously find those $2N$ summands, but with different prefactors.

Most importantly we note that the expectation value for the Wilson loop is finite to zeroth order already. If we now perform the same analysis for the half Wilson loop with a length $L$ larger than two links, we find a vanishing expectation value

[TABLE]

From Eq. (47) we generalize to even higher orders. We find that contributions purely from the Josephson term cannot be finite even in higher orders of perturbation theory, since the Pauli operators $\sigma^{z}$ can never be squared to unity or form a closed contour. The first term yields a finite contribution if the open loop $C_{1/2}$ is bridged by consecutive single electron transfers. One can infer that even for mixed terms in higher orders, the lowest order contribution is still given by the pure single charge tunnelling. Those are at least of order $\mathcal{O}(\left[E_{M}/(4E_{C})\right]^{L})$. Now we are in a position to calculate the expectation value of the Fredenhagen-Marcu operator to leading order in the perturbation theory. With the ground state $\ket{G^{*}}$ and normalisation $N^{*}$ of at least $\order{L/2}$ in the perturbation, this reads

[TABLE]

The Wilson loop as well as the normalization are approximated to leading order as $1$. Thus, the expectation value of the Fredenhagen-Marcu operator vanishes for $L\rightarrow\infty$ in the regime of dominant charging energy $E_{C}\gg E_{J},E_{M}$.

C.2 Dominant Josephson energy and large single-electron tunneling rate

For the case of dominant Josephson energy and large single-electron tunneling rate, $E_{J}\gg E_{M}\gg E_{C}$, we found the ground state in Eq. (18). Furthermore, the perturbation is given by

[TABLE]

We compute the first order perturbative ground state to

[TABLE]

We define the parameter $\alpha$ to simplify the notation. Since we are interested in the ground state expectation value of $R(L)$ for the limit of an infinite lattice we can, to a reasonable degree, approximate higher order corrections to the ground state as independent flips of matter DOF ($\tau^{x}_{i}$). Under this assumption of locally separate perturbations, we find for the ground state at order $N_{\text{max}}$

[TABLE]

At $n$-th order, we thus assume $n$ independent flips to occur, where the $n!$ is necessary to avoid double counting. The normalization is then calculated to

[TABLE]

Under the assumption of independent flips, the mixed order contributions have to vanish, as otherwise at least one $\tau_{s}^{x}$ does not square to unity. In particular, this means that the two sums in the calculation of the normalization in Eq. (C.2) simplify with a Kronecker delta $\delta_{nm}$ to

[TABLE]

We used the fact that there are $n!$ ways to arrange the sums over $j$ such that all $\tau_{i}^{x}$ square to unity. The empty sums run over all nodes resulting in the total number of nodes $N$. In the last step, we added the remaining terms from $N_{\text{max}}$ to infinity, which is a negligible difference for large enough $N_{\text{max}}$. Moving on to the expectation value of the Wilson loop, we obtain to first order with Eq. (50)

[TABLE]

Since the $\tau^{x}$ commute with the $\sigma^{z}$, we again obtain a Kronecker delta because the $\tau^{x}$ have to square to unity to result in a finite contribution. Thus, we repeat the above higher order calculation for the normalization factor for the Wilson loop to find

[TABLE]

In contrast to that, the half Wilson loop expectation value yields to first order

[TABLE]

In total $N-2$ of the $\tau_{j}^{x}$ commute with $\tau^{z}_{s}$ and $\tau_{s^{\prime}}^{z}$. In two terms the $\tau$ operators anticommute so that the expectation value differs from the Wilson loop. Assuming independent perturbations analogously to Eq. (55) results in

[TABLE]

Collecting Eq. (C.2), (55) and (57), we obtain for the Fredenhagen-Marcu operator

[TABLE]

Resubstitution of the parameter gives

[TABLE]

Thus, for $E_{C}\gg E_{M}\gg E_{C}$ the Fredenhagen-Marcu operator takes a constant value and the phase is confining.

Appendix D Perturbation analysis

In this section we give a more detailed description of the perturbation analysis for the case of $E_{C}\gg E_{M},E_{J}$. Let us start by giving some additional notation. We define the projectors

[TABLE]

The $P_{-}$ projects onto the ground state space and $P_{+}$ projects onto the remaining space. With

[TABLE]

$V_{\pm\mp}=P_{\pm}VP_{\mp}\qquad\text{and}\qquad H_{0+}=P_{+}H_{0}P_{+}$, the effective Hamiltonian is given by

[TABLE]

We neglected the constant contribution and linear corrections vanish. Note that $\ket{G}$ is in general highly degenerate, if we only take $H_{0}$ into account. In the presence of the perturbation, effectively introducing the plaquette terms and the gauge constraint, this degeneracy is lifted for the most part. The ground state then remains $4$-fold degenerate due to the topology of the system under consideration. Switching between these four states will only be possible in $L_{\text{min}}/2$-th order, where $L_{\text{min}}$ denotes the total length of the smallest lattice dimension. Thus, for reasonably low orders of perturbation theory there is no mixing and we can concentrate on one of the ground states. It was proven, that the ground states remain degenerate even in the perturbed case Gregor et al. (2011). To do this, they defined two ground states $\ket{G}$ and $\ket{G^{\prime}}=\Gamma\ket{G}$, where $\Gamma=\prod_{l\in\mathcal{C}_{nc}}\sigma^{z}_{l}$ denotes a non-contractible loop. Then one can show that

[TABLE]

by expressing $\ket{G^{\prime}}$ in terms of $\ket{G}$ and commuting the latter $\Gamma$ to the front such that it squares to identity. Therefore, we only perform the calculation for one representative ground state.

Since the Green’s function for the ground state energy corrections is negative, we define $\tilde{G}_{0+}(E)=-G_{0+}(E)$ and $\tilde{V}=-V$ to find

[TABLE]

This allows us to keep better track of the signs and we omit the $\sim$ again for convenience. In contrast to earlier works, where the perturbation was only performed up to fourth order, we have to expand $G_{0+}$ also in terms of $E$ around $E_{0}=0$. We can compute this with the operator identity

[TABLE]

as

[TABLE]

To capture all important terms up to sixth order, we now have to set

[TABLE]

since there is no linear contribution. The third order is furthermore neglected as it cannot lead to a plaquette contribution up to sixth order. The expansion has to be inserted back into Eq. (62) to find all the necessary terms. Collecting the results, we find

[TABLE]

where

[TABLE]

and

[TABLE]

Here, we set $G_{0+}=1/H_{0+}$. Without the energy expansion, we would obtain only negative contributions to the effective Hamiltonian.