Detecting Topological Entanglement Entropy in a Lattice of Quantum Harmonic Oscillators

The Kitaev surface-code model is the most studied example of a topologically ordered phase and typically involves four-spin interactions on a two-dimensional surface. A universal signature of this phase is topological entanglement entropy (TEE), but due to low signal to noise, it is extremely difficult to observe in these systems, and one usually resorts to measuring anyonic statistics of excitations or non-local string operators to reveal the order. We describe a continuous-variable analog to the surface code using quantum harmonic oscillators on a two-dimensional lattice, which has the distinctive property of needing only two-body nearest-neighbor interactions for its creation. Though such a model is gapless, satisfies an area law, and the ground state can be simply prepared by measurements on a finitely squeezed and gapped two-dimensional cluster state, which does not have topological order. Asymptotically, the TEE grows linearly with the squeezing parameter, and we show that its mixed-state generalization, the topological mutual information, is robust to some forms of state preparation error and can be detected simply using single-mode quadrature measurements. Finally, we discuss scalable implementation of these methods using optical and circuit-QED technology.

Topological order describes a phase of matter whose correlations satisfy an area law while maintaining long-range entanglement and ground-state degeneracy impervious to all local perturbations. These properties make such systems attractive candidates for stable quantum memories or processors [1]. However, the lack of a local order parameter makes measuring topological order an experimentally onerous task. Some possibilities include measuring non-local string operators [2] or the statistics of anyonic excitations above the ground state, as has been demonstrated experimentally with small photonic networks [3,4]. However, due to finite correlation lengths of local operators [2,5], these methods suffer from low visibility if the system is not prepared in a pure phase with vanishing two-point correlations.
An alternative is to study properties of the state itself that are robust to small changes in the correlation length. For a topologically ordered phase, the entanglement entropy of a subsystem in state ρ A is the von-Neumann entropy S(ρ A ) ≡ − tr[ρ A log 2 (ρ A )] = α|∂A| − γ + ε , where α ∈ R, |∂A| is the boundary size of A, and ε → 0 for |∂A| → ∞, [6][7][8]. The parameter γ is termed the topological entanglement entropy (TEE) [1], which is an intrinsically non-local quantity that characterizes topological phases in a variety of systems, including spin lattices such as the qubit surface code [6,9], bosonic spin liquids [10], and fermionic Laughlin states [11].
While useful for numerics, actually measuring TEE in a physical system is a daunting task since extracting the von Neumann entropy requires knowledge of the complete spectrum of the reduced state. A different option is to instead measure the Renyi entropy S (α) (ρ A ) ≡ 1 1−α log 2 tr[ρ α A ] since it was shown [12] that γ is the same when replacing S(ρ A ) with S (α) (ρ A ) ∀ α. The value α = 2 is an attractive choice since the purity tr[ρ 2 ] is observable via a simple swap-test measurement on two copies of the state [13]. A pure topological phase, such as the qudit (d-level spins) surface-code state [14], has tr[ρ 2 A ] = d 1−|∂A| , meaning γ = log(d) [15]. In contrast, the purity of another area law state with no TEE, such as the qu-dit cluster state [16], is tr[(ρ A ) 2 ] = d −|∂A| . Thus, even using Renyi entropy one still requires a number of measurements exponential in the size |∂A| to distinguish the two phases.
In this work we study, for the first time, topological order in a continuous-variable (CV) Gaussian state [17] analogue of the discrete-variable surface-code state. We describe how to prepare it efficiently using an intermediate mapping to first an ideal (infinitely squeezed) and then a physical (finitely squeezed) CV cluster state [18][19][20][21][22][23][24][25]. Remarkably, the TEE of the CV surface code can be easily computed simply from quadrature measurements. We show that unlike their qubit (or qudit) counterparts, the CV surface-code state has a parent Hamiltonian that is gapless in the thermodynamic limit. Nonetheless, the state is topologically ordered with a TEE that asymptotically grows linearly with the squeezing parameter. Other gapless models with topological order have been investigated in different contexts [26]. We propose experimental realizations for this model that are accessible with today's technology, and we conclude by analyzing the stability of the CV topological order against two forms of noise: thermalization (in the case of preparation by cooling) and noisy input states (in the case of active construction).
The ideal CV cluster state [18,19,27] is the CV analog of its qubit-based cousin [28,29] and may be obtained by sending zero-momentum eigenstates through pairwise controlled-Z gates CẐ j,k = e iq jqk in accord with an undirected, unweighted graph with one qumode (quantum mode-i.e., harmonic oscillator) per vertex. A CV cluster state is described completely by its nullifiers, which are comparable to stabilizers but with eigenvalue 0 instead of 1. Ideal CV cluster states have a complete set of nullifiers given by {η j = p j − ∑ k∈N ( j)q j }, where j runs over all vertices, and N ( j) is the neighborhood of j. From these nullifiers, we can construct a HamiltonianĤ ideal CS = ∑ jη 2 j , which has the ideal CV cluster state as its (non-normalizable) ground state. Note that H ideal CS has a continuous spectrum and is therefore gapless. See [19] and Appendix A for details.  Inspired by the dynamical mapping of qubit cluster states to qubit surface codes introduced in [30], there is a simple scheme [31] that transforms the CV cluster state into the corresponding CV surface code. First, start with the CV cluster state, and label vertices by row and column. Then, measure inp (inq) those vertices on rows and columns that are both odd (both even), as in Fig. 1a. This scheme is equivalent up to translation and/or inversion of thep,q measurements. After the measurements, we are left with a CV surface code state with a new set of nullifiers.
It is convenient to borrow notation from the qudit version of surface codes [14] that describes the nature of the coupling involved in the nullifiers in terms of a surface-code graph Λ = {V , E,F }. This graph defines an oriented surface where oscillators reside on the edges e ∈ E, and the oriented edges meet at vertices v ∈ V and surround oriented faces f ∈ F . Edge orientations are chosen such that for each vertex, all incident edges point toward it or all away from it (Fig. 1b). The new nullifiers areâ v = ∑ e|v∈∂eqe andb f = ∑ e∈∂ f o(e, f )p e , where o(e, f ) = ±1 if e is oriented with (+) or against (−) f . The CV surface code is the non-normalizable ground state of the quadratic HamiltonianĤ ideal SC = ∑ vâ † vâv + ∑ fb † fb f , which has a fully continuous spectrum and is therefore gapless.
With experimental realization in mind, we need to consider the physical case of finite squeezing. Squeezing is a Gaussian transformation, performed by the unitary operator S(s) = e − i 2 (log s)(qp+pq) , with s > 0, where log s is traditionally known as the squeezing parameter. In the Heisenberg picture,Ŝ(s) †qŜ (s) = sq, andŜ(s) †pŜ (s) =p/s, such that the variance ofp (ofq) after squeezing is a factor of s −2 (of s 2 ) times its original value. To produce the CV cluster state in the finite-squeezing case by the canonical method [18,19], which is the easiest method to work with theoretically, one starts from N vacuum modes |0 ⊗ N . Sinceâ j |0 ⊗N = 0 ∀ j, the initial nullifier set is {â j = 1 √ 2 (q j + ip j )}. These states are then all squeezed byŜ(s), followed by the pairwise CẐ couplings, yielding the transformed nullifierŝ These operators satisfy the canonical commutation relations for normal-mode operators, [η s j ,η s k ] = 0 and [η s j , (η s k ) † ] = δ j,k , and therefore the CV cluster-state Hamiltonian is [23] For finite s, the system has a gap of 2s −2 . The prefactor provides for finite energy even in the limit of infinite squeezing where lim s→∞ĤCS (s) =Ĥ ideal CS . Using the same measurement pattern as in the ideal case ( Fig. 1), the finitely squeezed CV cluster state can be mapped to the finitely squeezed CV surface code. Taking linear combinations of neighboring cluster-state nullifiersspecifically sums of neighboring nullifiers around thepmeasured nodes-and alternating signed cyclic sums around q-measured modes, one finds the general form of the surface code nullifiers [32]. Since the finitely squeezed cluster state is Gaussian, and quadrature measurements are Gaussian operations [33], so is the finitely squeezed surface code state. In the case of a square lattice with toroidal boundary conditions, the nullifiers arê where s = √ 5s 2 + s −2 . See Fig. 1b. We can construct a Hamiltonian using these nullifiers: The squeezing dependence of the prefactors is done to ensure the Hamiltonian has finite energy for s → ∞. Unlike the discrete-variable case, this Hamiltonian is gapless in the thermodynamic limit. This arises because the nullifiers do not define normal modes. Rather, neighboring nullifiers have nontrivial commutation relations, which allow for low-energy mode excitations. For a square n × m (n or m odd) lattice with n ≤ m the gap is ∆(s) ≈ 4π 2 /s 2 n 2 , and generically the system is gapless, see Appendix B. Hence, in distinction to the cluster-state HamiltonianĤ CS (s), the surface-code Hamilto-nianĤ SC (s) is gapless in the thermodynamic limit, though for infinite squeezing both models are gapless. A zero-mean N-mode Gaussian state is completely and uniquely described [34] via its symmetrized covariance matrix Γ j,k = Re tr[ρr jrk ], wherer = (q 1 , ...,q N ,p 1 , ...,p N ) T is a 2N-dimensional column vector of quadrature operators [34] FIG. 2: Sections used for calculations of topological entanglement entropy by the methods of (a) Kitaev-Preskill [7] and (b) Levin-Wen [8]. The areas of the regions satisfy the equality depicted.
. A Gaussian pure state's entanglement entropy can be calculated [35] in terms of the symplectic eigenvalues of Γ , which are the positive elements of the N eigenvalue pairs {±σ j } of the matrix product iΓΩ, with Ω j,k = −i[r j ,r k ]. The entropy for an N A -mode Gaussian subsystem ρ A is calculated using the reduced symplectic spectrum {σ A 1 , . . . , σ A N A } obtained deleting all of B's rows and columns from the covariance matrix.
In addition to the covariance-matrix representation, every zero-mean N-mode Gaussian pure state can be uniquely represented by an N-node, undirected, complex-weighted graph whose adjacency matrix is called Z = V + iU [33,36]. When Z is purely imaginary (V = 0), the state's covariance matrix becomes Γ ≡ 1 2 (Γ x ⊕ Γ p ) = 1 2 (U −1 ⊕ U). As shown in Fig. 1c, this is the case for the CV surface-code state, for which Z = iU SC = i(2s 2 + s −2 )I + is −2 A SC , where A SC is the corresponding unweighted adjacency matrix without self-loops. From this, we see immediately thatp-p correlations (determined by U) have range at most 1, andp-q correlations are zero. Further, as shown in Appendix E, theq-q correlations are bounded above by q iq j ≤ Ce −(d(i, j)+1)/ξ where d(i, j) is the shortest distance between modes i and j on the surface code graph. The constant is C = (1+ √ 8s 4 + 1) 2 /4(8s 2 +s −2 ) and the correlation length ξ is: Therefore, because for Gaussian states all higher order correlations are generated by the linear and quadratic ones, the CV surface code state obeys an area law. Numerically we find for squeezing log s = 3.2 the correlation length is ξ = 2.42, see Fig. 4, and the scale factor of entropy with area is α = 4.68.
To study the topological properties of the CV surface-code state, we make use of two alternative (but equivalent) definitions of TEE , one introduced by Kitaev and Preskill (KP) [7], and another by Levin and Wen (LW) [8], with regions shown in Figs. 2a and 2b, respectively. If the system is not topologically ordered, these combinations are exactly zero. Thus, we say that a model is topologically ordered only when γ > 0. In our simulation we use the ground state of a 36×36 mode CV surface-code and the squeezing-dependent values of the TEE are calculated selecting from the covariance matrix the reductions corresponding to regions chosen as prescribed by Eq. (7) and Eq. (8). See Fig. 3 for numerical results. Further, we plot the Topological Log-Negativity (TLN) for the KP regions based on recent results that show this quantity to be a good witness of topological order for stabiliser states of Abelian anyon models [39,40]. The log-negativity of the reduced state with support on a subsystem A is [41] where µ A = P A ⊥ ⊕ (−P A ), with P X being the projector onto the modes in region X and λ i being the i-th eigenvalue of the matrix argument. Remarkably, we find that the TLN is a rather tight upper bound on the TEE with the same asymptotic slope.
To derive this slope we consider the entanglement entropy for one mode of the smallest meaningful portion of the surface code, specifically a three-mode correlated state (see Appendix F for the details of the argument). The squeezing-dependent symplectic eigenvalue of the reduced covariance matrix corresponding to the mode used is given by σ 1 = 1 2 (1 + 3s 4 + 2s 8 ) 1/2 (1 + 3s 4 ) −1/2 . The mode entropy grows linearly with an asymptotic slope of lim s→∞ dγ(s) d(log s) = 2/ ln(2) 2.8854, which matches the slope of TEE and TLN we find numerically for larger systems. Interestingly, since the controlled-Z gates introduce additional squeezing [42,43], even for s → 1 (starting with vacuum states instead of momentum-squeezed states), the TEE is very small but non-zero.
To model noise in the state preparation, we consider a thermal state with respect to the cluster-state Hamiltonian, ρ CS (β) = e −βĤ CS (s) / tr[·] as the pre-measurement initial state. This could be generated by engineering the Hamiltonian H CS (s), which is gapped for finite squeezing, and then waiting until the system reaches equilibrium with an environment at temperature β −1 . Alternatively, using e.g. networks of non-interacting photons, one could start with separable modes each in a thermal input state ρ(β) = ∏ j e −β 2 s 2 a † j a j / tr[·] and then generate the thermal cluster state as before.
For such mixed states we can detect topological order by making use of the Topological Mutual Information (TMI) [15]. The TMI is constructed replacing in Eq. (7) the von Neumann entropy S X with (half of) the mutual information I X = S X + S X c − S X∪X c between a region X and its complement X c : Assume that, under a given sequence of symplectic transformations and homodyne detections, the ground-state covari-  7); the crosses (+) label the TEE using Eq. (8); the squares ( ) depicts the TLN; and the circles (•) indicate the lower bound for the topological mutual information (TMI) γ l MI for noisy state preparation. The main graph shows the maximum single-mode squeezing of 12.7 dB achieved to date [53,54]. Finally, the dashed line corresponds to TEE for the CV cluster state, which is always zero. Inset: TEE/TMI for levels of multimode squeezing with 5 dB marked as achievable with current optical technology [52]. The TLN is above the scale shown here.
ance matrix Γ CS of the CV cluster-state HamiltonianĤ CS (s) maps to a CV surface-code state Γ. Because all normal modes ofĤ CS (s) have the same frequency, the thermal covariance matrix ofĤ CS (s) at temperature β −1 is just κΓ CS and, under the same evolution, maps to κΓ, where κ = coth(βε 0 /2), and ε 0 is the energy gap ofĤ CS (s). Note this is not a thermal state forĤ SC (s) because the spectrum ofĤ SC (s) is nonuniform.
The TMI for this class of mixed states is lower bounded by the value computed as κ → ∞, as shown in Appendix F. This is γ l MI = − 1 2 ∑ X ζ(X) ∑ i log 2 (eσ X i ), where X runs over all the regions in Fig. 2a and their complements, ζ(X) = ±1 in accordance with Eq. (10), and the prime on the sum indicates that we need only include the zero-temperature symplectic eigenvalues σ X i for which σ X i > 1/2. (see Fig. 3). The CV cluster states considered in this work are canonical CV cluster states [33], so named because they are the states that would result if one were to use the canonical method of constructing them [18,19,27]. This method generates Gaussian states with graphs of the form Z = V + is −2 I, where the entries of V are either 0 or 1. While this method is straightforward theoretically, the CẐ gates are experimentally difficult and inefficient in an optical setting, where the most progress has been made [42,43].
More efficient and scalable optical construction methods exist [20-22, 24, 42, 45] but produce cluster states with uniform non-unit edge weight. These optical methods can produce very large states [50,51], including a recently demonstrated 10,000-mode cluster state with linear topology [52]. Very large square lattices with toroidal [21], cylindrical [21,24], and planar topology can also be made, as well as higher-dimensional lattices [46]. Similar states might also be created by cooling a circuit-QED system to the ground state ofĤ CS (s) [Eq. (2)] [23]. In Ref. [47] it was shown that using superconducting co-planar waveguides coupled pairwise via dissipative Cooper pair boxes, one can engineer an effectiveq-q interaction between the microwave modes in neighbouring waveguides. By changing the location of the box in the waveguides, one can also generatep-p couplings. While the cluster-state Hamiltonian in Eq. (2) also hasq-p couplings, there exist parent Hamiltonians for CV cluster states (up to phase shifts) that consist only ofq-q andp-p couplings [25,33].
All of these methods can be used to efficiently produce square-lattice CV cluster states with uniform edge weight g, with 0 < g < 1 4 . In the optical case, these states can be very large (thousands of modes) [22] and all produce a Gaussian state with a graph of the form Z g = gV + iεI. By squeezing each mode inq by a factor of Despite being constructed by a completely different method, the resulting state is a canonical CV cluster state with effective initial squeezings = g/ε. Since entanglement measures are local-unitary invariant, all of our results apply to these states if we take s =s. Furthermore, we don't need to actively perform the single-mode squeezing before measuring the TEE/TMI. We can simply rescale the outcomes of measurements on the original (g-weighted) state [48].
We have presented a model of correlated quantum harmonic oscillators in a topologically ordered continuous-variable surface-code state constructed using only Gaussian operations and with the remarkable property that its topological entanglement entropy γ can be observed simply via quadrature measurements. In contrast to discrete-variable systems, now γ is a continuous function of a system parameter, specifically the squeezing; however, this is not surprising since the parent Hamiltonian for the system is gapless. The CV surface code state can be prepared beginning in a gapped CV cluster state phase, and the topological entanglement entropy is robust to preparation errors modeled as thermal input. This provides a practical path to observe topological order in bosonic systems using current technology.
Appendix A: Infinitely squeezed CV cluster states and surface-code states In this appendix, we provide additional information about ideal (i.e., infinitely squeezed) CV cluster states [19] and ideal CV surface-code states [31].

Ideal CV cluster states
There are many ways to construct physical CV cluster states [18][19][20][21][22][23][24]42], all of which give slightly different states in the finitely squeezed case [33], Each has important differences that manifest when using them for measurement-based quantum computation [48]. In the infinitely squeezed case, however, these differences become largely irrelevant, and since infinitely squeezed states are unphysical anyway, we are free to choose for their analysis the method that is simplest from a theoretical perspective. For this reason, we choose the canonical method [18], which is the most straightforward generalization of the qubit cluster-state preparation procedure, despite its inefficiency when used in practice [42].
A CV cluster state is described in terms of a continuous generalization of the Pauli group, namely, the Weyl-Heisenberg group [19]. This generalization is most easily accomplished by thinking of the qubitσ x andσ z gates as implementing one-unit cyclic shifts in "position" (|0 ↔ |1 ) and "momentum" (|+ ↔ |− ), respectively. The CV analogs of the qubitσ x andσ z gates are the translation (position-shift) operatorsX(t) and boost (momentum-shift) operatorsẐ(u), respectively, with t, u ∈ R. While for qubits the shift is by an element of Z 2 , in the CV case, one may implement a shift in position or momentum by any real-valued amount. These shift operators are generated by the canonical self-adjoint quadrature operatorsp and −q, respectively, which satisfy [q,p] = i. Specifically, with the group commutator To construct a qubit cluster state [29], one starts with a collection of qubits in the |+ state (+1 eigenstate of the Pauli operatorσ x ) and applies controlled-σ z gates pairwise in accord with a chosen graph (traditionally, a square lattice), with one qubit per node of the graph. To construct an ideal CV cluster state by the canonical method [18], one begins with a collection of N quantum harmonic oscillators, each prepared in a 0-eigenstate of momentum, |0 ⊗N p , where the subscript means thatp j |0 p j = 0 for each mode. A CV cluster state |CS is generated from these initial states through the pairwise application of controlled-Z gates CẐ ( j,k) = e iq jqk upon all the nearest-neighbor modes j, k of the initial state, in accord with a given graph, with one mode per node of the graph. It is worth pointing that there is some additional freedom in this construction procedure, even at the ideal level. In particular, ideal CV cluster-state graphs can have a nonzero real-valued weight g ∈ R * associated with each edge. This modifies the strength of the CẐ gate represented by that edge: CẐ ( j,k) [g] := e igq jqk . These weights were first introduced in Ref. [20] as a way to enable new methods of constructing CV cluster states. They have shown themselves to be very important when considering the computational properties of these states [48] and when considering efficient construction of cluster states with very large graphs [22,24,25,33]. For the purposes of all derivations in this work, we set g = 1, but in the main text ("Experimental Implementation"), we showed how our results apply to CV cluster states constructed with non-unit-but still uniform-weight g.
Just like in the qubit case, ideal CV cluster states admit a description in terms of stabilizers. Given a stabilizer operatorK j for a state |ψ (i.e.,K j |ψ = |ψ ), under a unitary transformation of the state, |ψ =Û|ψ ,K j transforms asK j =ÛK jÛ † in order to preserve its status as a stabilizer:K j |ψ = |ψ . Note that the transformationK j −→K j under the action ofÛ is opposite from the Heisenberg evolution of observables under the same unitaryÛ. This is because we are not modeling the evolution of an observable when we evolve stabilizers. Instead, we are evolving the old stabilizer into a new stabilizer for the new state. Thus, the unitary evolution applied to the stabilizer must counteract that applied to the state in order to maintain the stabilizer's role as such.
To derive the Hamiltonian of the ideal CV cluster state, we construct the qubit case and then show the analogous steps in the CV case. Consider, for simplicity, a square-lattice graph with a qubit placed on each vertex. Each of the qubits is initially in a |+ state and hence is stabilized by aσ x operator:σ x j |+ j = |+ j . Next, the cluster state is created by applying a controlled-σ z gate on every pair of nearest-neighbor qubits. Consequently, the stabilizers are transformed into new elements equal tô with N (i) identifying the nearest neighbors of qubit j, and with N, S, E, and W indicating the qubit to the North, South, East, and West of qubit j, respectively. Finally, the qubit cluster state is the ground state of the Hamiltonian constructed by imposing an energy penalty for violating any of the stabilizer conditions:Ĥ To similarly define the ideal CV cluster state, we substitute the N qubits on the vertices of the square-lattice graph with N qumodes in the |0 p state as defined before. The global state is therefore |0 ⊗N p , and this is stabilized by the single-mode operatorsX j (s) in the sense thatX j (s)|0 ⊗N p = |0 ⊗N p . For CV state, there exists an equivalent way to express the stabilizer relations by using nullifiers [33]. An operatorη is called a nullifier for a state |ψ whenη|ψ = 0 holds. In this case, we havê and the state |0 ⊗N p is nullified by the set {p j } of generators of the stabilizer group. As in the qubit case, the ideal CV cluster state |CS is the result of the pairwise application of nearest-neighbor controlled-Z gates CẐ ( j,k) on |0 ⊗N p . Under the CẐ ( j,k) evolution, the quadrature operators transform as and thus the initial state stabilizers {X j (s)} are changed into the cluster state stabilizers {X j (s) ∏ k∈N ( j)Ẑk (s)}, which have the same form of the operators in Eq. (A4). We can also express the stabilizers by [33] equivalent to the elements of the CV cluster-state nullifier set {η j } being defined aŝ All theη j commute and are the elements of the algebra that generates the stabilizer group of |CS . Again, we can construct a HamiltonianĤ CS whose ground state is |CS by imposing an energy penalty for violating any of the nullifier con- Since all nullifiers {η j } commute and have a continuous spectrum of eigenvalues (R), this Hamiltonian also has a continuous spectrum, [0, ∞), and is therefore gapless.

Ideal CV surface-code states
Now we present a general description of the ideal CV surface code using a representation reminiscent of case for qudits [14]. As stated in the main text, consider a graph Λ = {V , E,F } where V is a vertex set, E is an edge set, and F is a face set. We assume the graph is oriented and that the faces inherit this orientation. Each quantum mode is located on an edge e j ∈ E, with the orientation of any edge e determined by e = [v, v ] for the base starting at vertex v and the head a vertex v .
The construction for the ideal CV surface-code state in the main text focussed on creating the state by starting with an ideal CV cluster state and then performing quadrature measurements on selected modes of the system in accord with a simple pattern [31] in analogy to the qubit case [30]. This was involved specializing to the simple case of Λ being a planar square lattice with edge orientations chosen such that at any vertex v, all incident edges point toward v, or all point away from v, and nullifiers were derived under these assumptions. This Appendix will diverge from the main text and make no such assumptions about the edge orientations.
In this more general case, the stabilizers for the surface code areÂ where the dot (·) stands for any vertex. By construction, the stabilizers commute: Oriented edges are useful when making the connection to qudits, but they are unnecessary in the special (and original) case of the qubit surface code [9]. In this case, the CV stabilizers in Eq. (A11) function once again as the CV analogs of the surface-code stabilizers for qubits: In the special case considered in the main text-i.e., when Λ is a planar square lattice with edge orientations chosen such that at any vertex v, all incident edges point toward v, or all point away from v, o(e, v) is a constant (±1) for any vertex v. Thus, it amounts to only a sign flip on t in Eq. (A11), which has no effect on a stabilizer's role as such, and in the main text we ignore it for simplicity. Under these conditions, the stabilizers in Eq. (A11) becomê where, as noted in the main text, the nullifiersâ v ,b f arê If, on the boundary, one of the edges is missing, then that mode is not included in the nullifiers. The Hamiltonian whose ground state is the ideal CV surface code is given bŷ While we have written this Hamiltonian as if it were that of a system of coupled oscillators, this connection is spurious in the ideal case [although it will be correct in the finitely squeezed case; see Eq. (B3)]. Since the "mode operators"â v andb f are actually (Hermitian) quadrature operators that all commute, this Hamiltonian has a continuous spectrum [0, ∞) and is therefore gapless for any number of systems.
Appendix B: Hamiltonian of the finitely squeezed CV surface code In this section we provide details on the form of the finitely squeezed CV surface-code Hamiltonian and its spectrum. The explicit form of the surface-code nullifiers for a square lattice (while allowing for smooth or rough open boundaries) are obtained by taking linear combination of neighboring nullifiers of the finitely squeezed CV cluster state. Given a set of exact nullifiers for a Gaussian state [33], one can obtain new nullifiers by a three-step process: 1. Given a quadrature measurementx j to be made on mode j, wherex ∈ {q,p}, using linear combinations of the original nullifiers, write a new set of nullifiers such that the canonically conjugate local quadratureŷ j (where [x j ,ŷ j ] = ±i) appears in only one nullifier in the new set.
2. In each new nullifier, replacex j with the real-valued measurement outcome.
We always assume that the outcome of the measurement is 0 because any other outcome would merely result in the same state up to displacements in phase space. These displacements can always be undone by local unitaries and therefore do not change any entanglement measure we might want to calculate. To obtain the face nullifiersb f , one sums the cluster-state nullifiers immediately adjacent to the node in question (along the cardinal directions) with orientation-dependent signs (e.g., η s N +η s S −η s E −η s W ). The vertex nullifiersâ v are more complex and require next-nearest-neighbor nullifiers to be added to the sum in order to achieve step 1 in the procedure above. The result for a lattice with possibly incomplete vertices and faces iŝ with s v = V (v)s 2 + s −2 , V (v) valence of vertex v and |∂ f | boundary size of the lattice face.
Note that now there is a dependence on the lattice position of the vertex or plaquette considered. The nullifiers satisfy the following commutation relations: Here d(v, v ) is the Euclidean distance between the two vertices where the edge lengths of the graph are unit length. The Hamiltonian is given bŷ and the squeezing dependence of the prefactors forĥ V and h F ensures the Hamiltonian has finite energy in the infinitely squeezed limit: Here we have used the fact that in the infinite-squeezing limit, each vertex nullifier involves a sum ofq's around that vertex and its four neighboring vertices, and since they all commute, the parent Hamiltonian is simply the squared sum ofq's around each vertex. Now we compute the gap of the surface code Hamiltonian, Eq. (B3). We first consider a square n × m lattice with toroidal boundary conditions: where d(v, v ) (respectively, d( f , f )) is the Euclidean distance between vertices (faces) on the unit-edge-length lattice (dual lattice): and On a torus |E| = 2nm,|F | = nm, and |V | = nm. We first focus on the case where n × m is odd such that there are |E| independent nullifiers which therefore span the space of all the physical-mode annihilation operators. Introducing normalmode operatorŝ where the vertices at the lattice sites have vertex coordinates {v r,s } and the sites of the dual lattice have face coordinates { f k,l }, the Hamiltonian iŝ To find the normal-mode frequencies, we need to solve the equations These two linear equations can be vectorized and rewritten as where |α ( j) is the vectorized form of the operatorĉ j , and |α ( j) is that ford j . Defining the shift operatorX r = ∑ r−1 k=0 |k ⊕ r 1 k|, we have The linear equations Eq. (B11) can be solved in the Fourier basis viaF n ⊗F m whereF r = 1 √ r ∑ r−1 j,k=0 e i2π jk/r | j k|.
treating j = ( j x , j y ) ∈ Z n × Z m as a collective index. The gap energy is the lowest-frequency mode energy: For large systems sizes, n, m 1, and choosing without loss of generality n ≤ m, the gap is If n, m even on the torus, then not all face nullifiers are independent (simply bicolor the faces, and assign plus signs to face operators on one color and minus signs on the other, and then add to get zero). Thus, the HamiltonianĤ SC (s) is underconstrained, and there is an exact gapless zero mode. For a square lattice with planar boundaries, there are boundary effects, but this makes only a small modification to the gap, which still scales like the inverse of the system size for large lattices.
Appendix C: Topological S-matrix for symmetric CV surface code One of the defining characteristics of topologically ordered matter is the existence of non-local interactions between particle-like excitations, which are described in terms of a scattering matrix. The scattering matrix can be evaluated in terms of expectation values of products of ribbon operators around a torus [1]. We show how these arise in the context of continuous-variable topologically ordered states.
For simplicity, consider the symmetric Hamiltonian mentioned in the text: For the case n, m even, because of the simple structure of the vertex and face nullifiers, there are only n − 1 independent vertex and m − 1 independent face nullifiers, so in fact these operators do not span the space of physical modes. Rather, there are two gapless modes. To see this, repeat the argument above of bicoloring the lattice (respectively, the dual lattice) and assigning +1 weight to vertex (face) nullifiers of one color, assigning−1 to the other, and adding to get zero. On the torus, there are non-local string symmetries of the ground-state subspace: where P andP are oriented paths on the lattice and dual lattice, respectively. Here o(e) = ±1 if the edge e is oriented in the same (opposite) direction as P , and f (e) = ±1 if the edge e has the same (opposite) framing to the pathP . The framing of the pathP is to the right, normal to its direction. Since each string touches an even number of modes of each nullifier with opposite signs due to edge orientations, then we have by construction: However, [Ẑ P (t),b † f ] = 0, and [XP (t),â † v ] = 0, so these string operators are not symmetries of the HamiltonianĤ SC (s).
However they are symmetries of the ground subspace H GS .
To see this notice that the normal mode operators for this symmetric model are defined as in Eq. (B8), call themĉ j andd j . Sinceĉ j andd j are linear combinations of the nullifiersâ v and b f operators respectively, we havê For contractable paths P andP , the generators ofẐ P (t) and XP (r) are linear combinations ofb f orâ v operators inside the loops. We can compute the topological S-matrix for this model by computing the monondromy [56]. We loosely refer to the string operators as worldlines for particles that are created out of the vacuum, make an excursion around one noncontractible loop on the torus, and then annihilate. There are two types of particles: electric charges with charge t associated with loopsẐ P (t) and magnetic fluxes with flux r associated with loopsXP (r). The scattering-matrix element associated with a t electric charge braiding around an r magnetic flux is where P 1,2 are two loops along the vertical direction, andP 1,2 are two loops along the horizontal direction. This can be evaluated by using the symmetry of the loop operators: which is the expected statistics based on the group commutator for the Weyl representation of the Heisenberg group: e −itq e irp e itq e −irp = e irt .
Here we have used the facts thatẐ −1 P (t) =Ẑ P (−t) and X −1 P (r) =XP (−r) and that the only nontrivial action is at four intersection mode locations {e a , e b , e c , e d }. Note there is a factor (−1) f (e a ) f (e c )+ f (e d )+o(e a )+o(e c )+o(e d ) also in the exponent which we have assumed is equal to one. If for our chosen edge orientations this is −1, we simply redefine either r → −r or t → −t. Note that in the limit of infinite squeezing, s → ∞, the string operators become unitary.

Appendix D: General properties of Gaussian states, Z matrix and entanglement entropy
Both qubit cluster states and the qubit surface code can be represented by graphs. Strictly speaking, however, their definitions in terms of graphs are incompatible. In particular, qubit cluster states (or qubit "graph states," as they are often confusingly called in the literature) are represented by graphs that act as a well-defined recipe for creating the state: nodes represent qubits prepared in |+ , and edges represent CPHASE gates between them. While one does not have to make cluster states this way, any cluster state can be so constructed. If we try to apply this recipe to the graphs for qubit surface-code states, however, we fail because for the surface code, qubits live on edges, while plaquettes and vertices define stabilizers in terms of allσ z 's or allσ x 's. There is no analog for this in the graph recipe used to define qubit cluster states. While one can define by fiat a connection between the two types of graphs, it is at best a patch-up job.
In contrast, the graphical calculus for Gaussian pure states [33] provides a unified graphical representation of any Gaussian pure state in terms of a recipe for its creation. Again, one does not have to follow the recipe to make a Gaussian state, but any Gaussian state can be made by the recipe defined (uniquely) by its graph. As such, for CV cluster states and the CV surface code, only one type of graph is needed. Furthermore, the measurement patterns that connect the two types of states have definite graph transformation rules, so the connection between the two graphs can be derived using the graphical calculus instead of invented ad hoc as a patch between two inequivalent types of graphs, as in the qubit case. This elegant connection is described in detail in this Appendix, where we also review the main features of Gaussian states.
Entanglement entropy of a generic state can be rather hard to compute since it generally requires knowledge of the spectrum of the density matrix. On the other hand, a zero-mean, Nmode Gaussian state can be described conveniently and completely by an easy algebraic formalism that follows from the form of its characteristic function [57], which is solely a function of the vector of the first statistical moments and the matrix Γ that carries the information about the second moments. Recall that Γ is the covariance matrix of the Gaussian state defined as wherer = (q 1 , ...,q n ,p 1 , ...,p n ) T is the 2N-dimensional column vector of the hermitian quadrature operators of the N modes. The information contained in the covariance matrix completely determines the entanglement properties of a Gaussian state [35]. Explicit calculations of Gaussian states entanglement entropy are performed making use of the symplectic spectrum of Γ. Let us introduce the symplectic form Ω, which is a skew-symmetric matrix that incapsulates the canonical commutation relations of the quadrature operators. For a Gaussian state ρ with covariance matrix Γ, the positive elements of the N pairs of eigenvalues {±σ i } of the matrix product iΓΩ are called symplectic eigenvalues. The connection between entanglement entropy S(ρ) and symplectic eigenvalues of ρ is given by the formula The Gaussianity of the state is preserved under Gaussian unitary transformationsÛ G . Each of these operations has a correspondent matrix representation Y that belongs to the Symplectic group Sp(2n, R), while quadratures measurements have a well-defined action on the covariance matrix Γ of the state, which is described below. A symplectic transformation Y preserves the canonical commutation relations as follows while the action of a Gaussian transformationÛ G on the quadratures can be expressed bŷ where the right-hand side corresponds to a matrix multiplication on the quadratures vector. At the level of the covariance matrix, this is reflected in the transformation rule A state of N independent vacua is described by the covariance matrix After a Gaussian unitary transformationÛ Y represented by the matrix Y is applied to Γ 0 , the resulting Gaussian state is Exploiting the decomposition properties of symplectic matrices [33,36], the product YY T is uniquely specified by where, for an N-mode state, both U and V are N × N symmetric matrices, and U > 0. Hence, the complex linear combination offers an alternative description for a pure Gaussian state. The graph Z shows up directly in the position-space wavefunction ψ Z (q) for an N-mode Gaussian state |ψ Z : where q = (q 1 , . . . , q 2 ) T is a column vector of c-number position variables. For this state, Γ Y from Eq. (D9) can be rewritten as The matrix Z has a simple transformation rule under a symplectic transformation Y. If Y is decomposed into block form then the Z matrix associated to the transformed state is given by The usefulness of this approach lies in the simple transformation rules of Z for the most common laboratory procedures corresponding to Gaussian unitary transformations, as listed in [33]. This permits the study of Gaussian state evolution simply in terms of appropriate transformations on Z. Furthermore, it is possible to define the Z CS matrix for the CV cluster state directly, solely making use of the adjacency matrix A d that describes the square-lattice pattern of connections among the modes. Explicitly: Applying the measurement pattern prescribed in Fig. 1a in the main text to the state represented by the graph in (a) produces a state with the graph in (b). In both graphs, log s is called the squeezing parameter as used throughout this document. The darkness of each line is proportional to its magnitude, with color indicating phase (red = positive real, cyan = positive imaginary).
with squeezing parameter log s and I N the N × N unit matrix. This is illustrated in Fig. 5a. For the finitely squeezed cluster state with open boundary conditions, we start with Z CS (s) and perform the measurement scheme displayed in Fig. 1 in the main text. Measurements have a straightforward translation in the Z transformation rules language. Aq measurement on the kth mode is equivalent to deleting the kth row and column of the Z matrix, while ap measurement is equivalent to applying a π/2 phase shift on the kth mode and then measuring theq quadrature. At the level of the Z matrix, anyp measurement deletes the measured mode while generating new connections among its nearest neighbors. The Z CS (s) transforms into the CV finitely squeezed surface-code state Z SC (s), shown in Fig. 5b. In this case, it is a matrix having only imaginary entries given by with and A SC adjacency matrix of the surface code. This immediately provides a simple expression for the covariance matrix of the newly evolved state: Complete knowledge of the covariance matrix makes straightforward calculations of entanglement entropy for different regions of the system. It is enough to derive the proper reduced covariance matrices, calculate the symplectic eigenvalues, and apply the formula from Eq. (D3). The covariance matrix Γ SC (s) allows us to study theq andp-correlations for the modes of the surface code. From Eq. (D17) one can immediately notice that thep-correlations are non-zero exclusively in the nearest neighborhood of each mode, due to the particular form of the matrix U SC (s). The graph in Fig. 4 shows the behavior of theq-correlations along the main axes of the system described by Eq. (D18), which are not simple because U −1 SC (s) is complicated. These correlations decay exponentially with the distance.
Appendix E: Correlations on the lattice A necessary assumption to use the topological entropy formulas is the exponential decay of the quadrature correlations. The covariance matrix Γ SC (s) allows us to study thê q andp-correlations for the modes of the surface code. From Eq. (D17) one can immediately notice that thep-correlations are non-zero exclusively in the nearest neighborhood of each mode, due to the particular form of the matrix U SC (s). This corresponds to a correlation length ξ p = 1. On the other hand, theq-correlations require a more elaborate analysis because the matrix U −1 SC (s) is more complicated. It is however possible to prove analytically that also thê q-correlations decay quickly enough. First we show that the spectral range of U SC (s), denoted σ(U SC (s)), satisfies σ(U SC (s)) ⊂ [a, b] where a = s −2 and b = s 2 (8 + s −4 ). For the minimum eigenvalue, note that for finite squeezing the matrix U SC (s)/s 2 is positive definite but for s → ∞, someq −q correlations become infinite. This indicates that the matrix is singular in that limit and the smallest eigenvalue is therefore zero. Adding finite squeezing shifts the spectrum of U SC (s)/s 2 by s −4 so the minimum eigenvalue of U SC (s) is a = s −2 .
To derive the largest eigenvalue of U SC (s), observe that for the surface code on a lattice with periodic boundaries, the Z graph associated with the adjacency matrix A SC is regular with degree 6, i.e. each node connects to six others. The largest eigenvalue of the adjacency matrix of a regular graph is equal to the degree with an associated eigenvector ν = 1, ..., 1 [58], therefore the maximum eigenvalue of A SC is 6. From these considerations it follows that the largest eigenvalue of U SC (s) is b = s 2 (8 + s −4 ). For lattices with open boundary conditions, the eigenvalues of A SC are upper bounded by the maximal degree with corrections that fall off with the system size, so the spectrum of U SC (s) lies in the interval [a, b].
On a n × m lattice, the nm × nm matrix U SC (s) is block tridiagonal with n identical m × n matrices A on the diagonal and identical m × m matrices B on the immediate upper and lower blocks. Here the matrix coordinates (i, j) correspond to Euclidean coordinates ((i x , i y ), ( j x , j y )) on the lattice where i = mi x + i y for i x ∈ {0, ..., n − 1} and i y ∈ {0, ..., m − 1}, etc. It is convenient to define a graph distance d(i, j) = max{|i x − j x |, |i y − j x |} between coordinates (i x , i y ) and ( j x , j y ). Since away from the edges the Z graph is a union of square graph with a graph having two diagonal edges passing though every other plaquette, the graph distance is the number of edges on the shortest path between (i x , i y ) and ( j x , j y ) and it satisfies ed(i, j)/ √ 2 ≤ d(i, j) ≤ ed(i, j) where ed(i, j) = (i x − j x ) 2 + (i y − j y ) 2 is the Euclidean distance. The matrix A is itself tridiagonal with elements α = 2s 2 + s −2 on the main diagonal and β = s 2 on the immediate upper and lower diagonal. The matrix B is also tridiagonal with diagonal elements γ = s 2 and immediate upper and lower diagonal elements either equal to zero or γ. A theorem of Demko, Moss, and Smith [37,38] shows that banded matrices of a certain class have inverses with matrix elements that decay exponentially with the distance from the diagonal. Specifically they show for matrices M of size N × N and spectral range σ(M) ⊂ [a, b] with a > 0 [Ref. [37] Proposition 5.1]: where the decay sets are The matrices relevant to our problem are in this class. The matrix power U k SC (s) is a banded block symmetric matrix with blocks of size m × m and block band width 2k + 1. Furthermore, each such block is banded with band width 2k + 1. Thus the support set S p (U SC (s)) is the set of those matrix coordinates (i, j) such that the graph distance d(i, j), is no more that 2p + 1. Similarly, the decay set is all matrix coordinates outside the support set. The statement in Eq. E1 is that for nodes separated in graph distance d(i, j) > 2p + 1 with associated matrix coordinates (i, j) the inverse matrix element U −1 SC (s) i, j falls off exponentially with graph distance as q (d(i, j)+1)/2 . This implies that theq-correlations between nodes i and j separated in graph distance by d(i, j) satisfy where the constant is and the correlation length is Appendix F: Topological entanglement entropy Let us give a physical interpretation of the different formulas used to identify the TEE. Topological order is a phase of matter resulting from non-local long-range correlations of topological origin that extend along all the modes of such systems. For this reason, local parameters cannot detect it, and therefore more effort is required in order to identify it. The first hint regarding the nature of the TEE was given in [6], where considerations about the toric code led the authors to note that the entanglement entropy of a region of the system is proportional to the perimeter of that region minus a constant. This idea was then broadened, and interesting relationships among entanglement entropy, total quantum dimension of the model, and TEE were introduced. The topology of the surface where a topologically ordered model is defined determines the species of anyons that are supported by the model. If we assign to each anyonic species a quantum dimension d, and we sum over all the quasiparticles species, we call the total quantum dimension the quantity Both Kitaev and Preskill (KP) [7] and Levin and Wen (LW) [8] had the intuition that for a topologically ordered phase the entanglement entropy of a region scales as where α ∈ R, |∂A| is the size of the boundary of A, ε is a contribution that goes to zero in the limit of |∂A| → ∞, and is the TEE. KP focus on constructing a linear combination of entanglement entropies such that all the terms proportional to the length of the boundaries of the regions cancel out, and only the topological term remains. Specifically, this corresponds to: Once the regions are chosen to be reasonably bigger than the correlation length of the system, it is also possible to demonstrate that deformations of the boundaries of the regions and smooth deformations of the Hamiltonian do not change the value of γ. This means that γ is a quantity that is both invariant (i.e., only determined by the underlying topology) and universal (i.e., local modifications of the Hamiltonian do not affect the global topological properties). LW moved along a different line of thought, considering partitions of the system with different topologies rather than more general regions as KP, see Fig. 2. The partitions are constructed such that their pairwise differences are equal, and chosen to be large enough, so that short-range correlations decay inside the boundaries. In analogy to Eq. (F4), an explicit expression for the TEE can be derived: Again, if the system is not topologically ordered, and if the entanglement entropies only depend on terms proportional to the boundaries, then the difference in Eq. (F5) is zero. On the other hand, if topological long-range correlations are present on the lattice, LW argue that non-local closed string operators with non-vanishing expectation value exist on the lattice. Operators that wind around the region A must then contribute to a nonzero value of Eq. (F5). They indeed call this contribution TEE, which is a witness for topological order in the following sense: the value of γ characterizes the global anyonic and entanglement properties of the topological state. When γ = 0, then it follows from Eq. (F3) that D = 1, which physically means that no anyons are supported by the system, and no long-range topological correlations contribute to the entanglement. The idea behind both the KP and LW derivations of TEE formulas is to suppress the non-topological correlations, extracting only the topological information.
For mixed states, the appropriate signature of topological order is given by a quantity called topological mutual information (TMI), introduced in [15]: using the same regions used by KP for the calculations of the TEE. Now the correlations are measured using (one half of) the mutual information I A between a region A of a system and its complement A c , i.e., I A = S A + S A c − S A∪A c , which replaces the von Neumann entropy S A (up to a factor of 2).
In [44] a modified version of the TMI was proposed, to avoid treacherous ambiguities in the recognition of topological order. In fact, certain mixed states (in particular mixture of degenerate ground states of gapped local Hamiltonians) can exhibit long-range non topological correlations and therefore affect the value of the TMI as it was defined in [7,8]. Together with a new formulation for the TMI, they also introduced lower and upper bounds that precisely determine the TMI whenever they coincide. Consider now the LW regions as defined in Fig. 2 and label E = A \ B, set difference of A and B, F 1 = A \C, and F 2 = D (from which follows F = F 1 ∩ F 2 = C). Then the value of the TMI lies between these two bounds: The circles (•) in Fig. 3 in the main text, are a plot of the lower bound in Eq. (F7) in the worst-case scenario of very high initial temperature (κ → ∞); see Appendix G. This value for the bound is the lowest possible bound when creating the surfacecode state from a thermal cluster state. It therefore illustrates the maximum extent to which the TMI can sink below the TEE for any given value of the squeezing parameter for the particular construction procedure described in the main text.

Upper bound on TEE
Here we compute an upper bound on the TEE for the CV surface-code state by calculating the subsystem entropy of a simpler network of entangled modes. To do this, we invoke a the calculation of subsystem entropy appropriate to stabilizer states, which have vanishing two-point correlation functions. Consider stabilizer states that are quantum doubles of a finite group G (such as the toric code with group G = Z 2 ). As shown in Ref. [6], the TEE is calculated by dividing the system into two subsystems A, B and identifying the redundant gauge transformations defined on the boundary between FIG. 6: Simplified network of modes with which to derive an upper bound on the TEE for the CV surface code: (a) This shows a surfacecode state for a quantum double model with discrete variables. There are six vertices, two faces, and seven edges where physical modes reside. Not all the vertex stabilizers are independent. Rather, their product is the identity, so the number of independent stabilizers is 5 + 2, equal to the number of physical modes. We can deform the lattice while preserving the topological order by replacing three modes of each face, excluding the mode on the shared boundary, with one mode as shown in (b). This network has one independent vertex stabilizer, and two independent face stabilizers which equals the number of physical modes. For the toric code, the network represents a ground state that is a GHZ state. (c) A finitely squeezed cluster-state graph with seven modes maps to the three-mode network by measuring the grey modes in theq basis and the black modes in thep basis.
In fact, the nullifiers on the top and bottom black modes of the graph act equivalently on the white modes meaning one of them is redundant. Therefore, an even simpler CV cluster-state graph suffices as depicted in (d). Upon measurements, the reduced network has three modes with a correlation matrix that can be computed exactly to yield an upper bound to the CV surface code TEE.
the two regions. The entanglement entropy of subsystem A is the logarithm of the number of the (all equivalent) Schmidt coefficients of the state. Exploiting the group properties of G allows one to write the entropy as S(ρ A ) = (|∂A| − 1) log 2 |G|, implying γ = log 2 |G| = log 2 (D).
For the CV surface codes, it is complicated to extract an analogous exact expression for the entropy of a subsystem because the Schmidt coefficients are not equal as in the discrete case. Furthermore, the TEE is infinite for infinitely squeezed CV surface code states, and the definition of quantum dimension is not so clear for finitely squeezed CV surface code states since we do not yet have a description of this model in terms of a quantum double of a group. Nevertheless, we can go ahead and compute the subsystem entropy in the same way that would be done for the discrete case and treat this as a bound for the TEE of the CV surface code state. It is simply an upper bound because we are ignoring longer-range correlations that degrade the topological order, but since the correlation length is bounded for any finite amount of squeezing, this should be a reasonably tight bound.
A simple configuration to start with is a quantum double model with a discrete group on a lattice with two faces. This can be realised using a graph with just three edges (physical modes) and two vertices, as shown in Figs. 6a,b. For the toric code the ground state would be the GHZ state (|000 + |111 )/ √ 2 since both vertices implement the stabilizerσ x 1σ x 2σ x 3 , and one face enforcesσ z 1σ z 2 , while the other face enforcesσ z 2σ z 3 . Identifying two qubits on one of the faces with subsystem B, the subsystem entropy is S(ρ A ) = S(ρ B ) = 2 − γ = 1, where 2 comes form the size of the boundary of re-gion B and therefore γ = 1. This simplified surface-code network of three modes can be obtained from a finitely squeezed CV cluster state with six modes after measuring out three of the modes (Fig. 6d). The resultant CV network has a correlation matrix that can be computed exactly. The symplectic spectrum for one subsystem consisting of one mode has two eigenvalues {±σ 1 } with σ 1 = 1 2 1 + 3s 4 + 2s 8 1 + 3s 4 .
Appendix G: Surface-code TMI from high-temperature CV cluster states In this appendix, we derive a bound for the topological mutual information (TMI), analyzing the limit κ → ∞, which corresponds to the strongest possible decrease of the TMI from the TEE for the noise model considered. Recall that for a reduction ρ A of a pure state ρ, the von Neumann entropy S(ρ A ) determines the entanglement entropy of the subsystem with respect to its complement. For Gaussian states, one can use the formula where {σ i } A is the collection of n > A symplectic eigenvalues associated to the reduced covariance matrix Γ A of the subsystem ρ A . About the notation, in the following n ≥ X = n X indicates the total number of symplectic eigenvalues ≥ 1 2 associated to a region X, n > X indicates the number of eigenvalues > 1 2 , and n = X denotes those = 1 2 . Noise in our scheme has been modeled as beginning with a thermalized CV cluster state rather than a pure one. Since the normal-mode energies of the CV cluster-state Hamiltonian (Eq. (2) in the main text) are all equal, this is equivalent to squeezing identical thermal states instead of vacuum states in the canonical construction procedure [18]. The inverse temperature β of the state defines a useful parameter κ = coth(β/s 2 ). To detect the topological order of the resulting mixed CV surface-code state, we use the TMI: The numerics show that the TMI does not decrease significantly with an increment of the initial value of κ. This is interesting and rather unintuitive, hence an analytical expression is required to confirm our findings. First of all, recall that the covariance matrix of the resulting mixed CV surface-code state is equivalent to the pure one (Γ 0 ) times κ, mathematically Γ = κΓ 0 . (G3) As a consequence, the symplectic eigenvalues of Γ (or any reduced section of it) are simply given by the "pure" symplectic eigenvalues multiplied by the overall κ factor, {σ i } = {κσ 0 i }. This simple transformation of Γ is a special case that only arises due to the fact that all normal modes of the CV clusterstate Hamiltonian are identical, resulting in equal symplectic eigenvalues κ/2.
To start, consider the first term of Eq. (G2): where the regions are shown in Fig. 7. If the total state ABCD has N modes, then, for region A, we have n A modes and for its complement n BCD modes, such that n A + n BCD = N. Hence, the von Neumann entropy for A is given by In the limit of very high temperature, corresponding to κ → ∞, we can use the approximation (G4) and rewrite the von Neumann entropy S A as Divide now the n A symplectic eigenvalues into the two sets n > A , n = A , such that n > A + n = A = n A . Hence, we can rewrite the expression in Eq. (G7) as: We can repeat the same argument for the region BCD and find that The N symplectic eigenvalues of ABCD are all equal to κ/2. Consequently, and the value for the mutual information I l A in the κ → ∞ limit is given by Notice that the κ-contributions cancel out exactly. Using n = A + n = BCD = N − n > A − n > BCD and the area law behavior for the entropy, i.e., for regions sufficiently big, n > A = n > BCD (although this does not mean that the sets {σ A i } and {σ BCD i } are the same), we can rewrite the mutual information as (G12) The same argument applies for each other contribution to the TMI-for example, for region C: