Detecting topological entanglement entropy in a lattice of quantum harmonic oscillators

In the following discussion we only consider Gaussian states. A zero-mean, N-mode Gaussian state can be described conveniently and completely by an easy algebraic formalism that follows from the form of its characteristic function [32], which is solely a function of the vector of the first statistical moments and the matrix $\Gamma$ that carries the information about the second moments. Recall that $\Gamma$ is the covariance matrix of the Gaussian state defined as

$\begin{eqnarray}&&{{\Gamma }_{j,k}}={\rm Retr}[\rho \;{{\widehat{{\bar{r}}}}_{i}}{{\widehat{{\bar{r}}}}_{j}}],\end{eqnarray} \tag{ 7 }$

where $\widehat{{\bar{r}}}={{\left( {{{\hat{q}}}_{1}},...,{{{\hat{q}}}_{n}},{{{\hat{p}}}_{1}},...,{{{\hat{p}}}_{n}} \right)}^{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 [33]. Explicit calculations of Gaussian states entanglement entropy are performed making use of the symplectic spectrum of $\Gamma$. Let us introduce the symplectic form $\Omega$,

$\begin{eqnarray}&&{{\Omega }_{j,k}}=-i[{{\widehat{{\bar{r}}}}_{i}},{{\widehat{{\bar{r}}}}_{j}}],\end{eqnarray} \tag{ 8 }$

which is a skew-symmetric matrix that incapsulates the canonical commutation relations of the quadrature operators. For a Gaussian state ρ with covariance matrix $\Gamma$, the positive elements of the N pairs of eigenvalues $\{\pm {{\sigma }_{i}}\}$ of the matrix product $i\Gamma \Omega$ are called symplectic eigenvalues. The entropy for an N_A -mode Gaussian subsystem ${{\rho }_{A}}$ is

$\begin{eqnarray}&&S({{\rho }_{A}})=\mathop{\sum }\limits_{\{\sigma _{i}^{A}\}}\left[ (\sigma _{i}^{A}+\frac{1}{2}){{{\rm log} }_{2}}(\sigma _{i}^{A}+\frac{1}{2})-(\sigma _{i}^{A}-\frac{1}{2}){{{\rm log} }_{2}}(\sigma _{i}^{A}-\frac{1}{2}) \right],\end{eqnarray} \tag{ 9 }$

calculated using the reduced symplectic spectrum $\{\sigma _{1}^{A},\ldots ,\sigma _{{{N}_{A}}}^{A}\}$ obtained deleting the complementary rows and columns from the covariance matrix.

We define a Gaussian operation any unitary transformation that maps a Gaussian state onto a Gaussian state [33]. These are generated by Hamiltonians at most quadratic in the quadrature operators. Note that all CV Clifford operations are Gaussian operations. Each Gaussian operation ${{\hat{U}}_{G}}$ has a correspondent matrix representation ${\bf Y}$ that belongs to the Symplectic group Sp $(2n,\mathbb{R})$, while quadratures measurements have a well-defined action on the covariance matrix $\Gamma$ of the state, which is described below.

A symplectic transformation ${\bf Y}$ preserves the canonical commutation relations as follows

$\begin{eqnarray}&&{\bf Y}\Omega {{{\bf Y}}^{T}}=\Omega ,\qquad \forall \;{\bf Y}\in Sp(2n,\mathbb{R}),\end{eqnarray} \tag{ 10 }$

while the action of a Gaussian transformation ${{\hat{U}}_{G}}$ on the quadratures can be expressed by

$\begin{eqnarray}&&{{\widehat{{\bar{r}}}}^{\prime }}=\hat{U}_{G}^{\dagger }\widehat{{\bar{r}}}{{\hat{U}}_{G}}\longrightarrow {{\widehat{{\bar{r}}}}^{\prime }}={\bf Y}\widehat{{\bar{r}}},\end{eqnarray} \tag{ 11 }$

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

$\begin{eqnarray}\begin{array}{rcl} {{\Gamma }^{\prime }}=\operatorname{Re}\langle {{\widehat{{\bar{r}}}}^{\prime }}{{\widehat{{\bar{r}}}}^{\prime \,T}}\rangle \longrightarrow {{\Gamma }^{\prime }} & = & \operatorname{Re}\langle {\bf Y}\widehat{{\bar{r}}}{{({\bf Y}\widehat{{\bar{r}}})}^{T}}\rangle \\ {} & = & {\bf Y}\;\operatorname{Re}\langle \widehat{{\bar{r}}}{{\widehat{{\bar{r}}}}^{T}}\rangle {{{\bf Y}}^{T}} \\ {} & = & {\bf Y}\Gamma {{{\bf Y}}^{T}}. \\ \end{array}\end{eqnarray} \tag{ 12 }$

A state of N independent vacua is described by the covariance matrix

$\begin{eqnarray}&&{{\Gamma }_{0}}=\frac{1}{2}{{{\bf I}}_{2N}}.\end{eqnarray} \tag{ 13 }$

After a Gaussian unitary transformation ${{\hat{U}}_{{\bf Y}}}$ represented by the matrix ${\bf Y}$ is applied to Γ₀, the resulting Gaussian state is

$\begin{eqnarray}&&{{\Gamma }_{{\bf Y}}}=\frac{1}{2}{\bf Y}{{{\bf Y}}^{T}}.\end{eqnarray} \tag{ 14 }$

Exploiting the decomposition properties of symplectic matrices [31, 34], the product ${\bf Y}{{{\bf Y}}^{T}}$ is uniquely specified by ${{{\bf Y}}_{({\bf UV})}}{\bf Y}_{({\bf UV})}^{T}$, with

$\begin{eqnarray}{{{\bf Y}}_{({\bf UV})}}=\left( \begin{array}{ccccccccccccccc} {{{\bf U}}^{-1/2}} & 0 \\ {\bf V}{{{\bf U}}^{-1/2}} & {{{\bf U}}^{1/2}} \\ \end{array} \right),\end{eqnarray} \tag{ 15 }$

where, for an N-mode state, both ${\bf U}$ and ${\bf V}$ are $N\times N$ symmetric matrices, and ${\bf U}\gt 0$. Hence, the complex linear combination

$\begin{eqnarray}&&{\bf Z}:={\bf V}+i{\bf U},\end{eqnarray} \tag{ 16 }$

offers an alternative description for a pure Gaussian state. The graph ${\bf Z}$ shows up directly in the position-space wavefunction ${{\psi }_{{\bf Z}}}({\bf q})$ for an N-mode Gaussian state $|{{\psi }_{{\bf Z}}}\rangle$:

$\begin{eqnarray}&&{{\psi }_{{\bf Z}}}({\bf q})={{\pi }^{-N/4}}{{\left( {\rm det} {\bf U} \right)}^{1/4}}{\rm exp} \left( \frac{i}{2}{{{\bf q}}^{T}}{\bf Zq} \right),\end{eqnarray} \tag{ 17 }$

where ${\bf q}={{\left( {{q}_{1}},\ldots ,{{q}_{2}} \right)}^{T}}$ is a column vector of position-space variables. For this state, ${{\Gamma }_{{\bf Y}}}$ from equation (15) can be rewritten as

$\begin{eqnarray}{{\Gamma }_{{\bf Y}}}=\frac{1}{2}\left( \begin{array}{ccccccccccccccc} {{{\bf U}}^{-1}} & {{{\bf U}}^{-1}}{\bf V} \\ {\bf V}{{{\bf U}}^{-1}} & {\bf U}+{\bf V}{{{\bf U}}^{-1}}{\bf V} \\ \end{array} \right).\end{eqnarray} \tag{ 18 }$

The matrix ${\bf Z}$ has a simple transformation rule under a symplectic transformation ${\bf Y}$. If ${\bf Y}$ is decomposed into block form

$\begin{eqnarray}{\bf Y}=\left( \begin{array}{ccccccccccccccc} {\bf A} & {\bf B} \\ {\bf C} & {\bf D} \\ \end{array} \right),\end{eqnarray} \tag{ 19 }$

then the ${{{\bf Z}}^{\prime }}$ matrix associated to the transformed state is given by

$\begin{eqnarray}&&{{{\bf Z}}^{\prime }}=({\bf C}+{\bf DZ}){{({\bf A}+{\bf BZ})}^{-1}}.\end{eqnarray} \tag{ 20 }$

The usefulness of this approach lies in the simple transformation rules of ${\bf Z}$ for the most common laboratory procedures corresponding to Gaussian unitary transformations, as listed in [31]. This permits the study of Gaussian state evolution simply in terms of appropriate transformations on ${\bf Z}$.

Measurements also have a straightforward translation in the ${\bf Z}$ transformation rules language. A $\hat{q}$ measurement on the kth mode is equivalent to deleting the kth row and column of the ${\bf Z}$ matrix, while a $\hat{p}$ measurement is equivalent to applying a $\pi /2$ phase shift on the kth mode and then measuring the $\hat{q}$ quadrature. At the level of the ${\bf Z}$ matrix, any $\hat{p}$ measurement deletes the measured mode while generating new connections among its nearest neighbors.

The CV surface code can be easily described and efficiently prepared using an intermediate mapping to the CV cluster-state via a simple pattern of quadrature measurements. Intuitively, the CV cluster-state is the CV analog of its qubit-based cousin [36]. There are many ways to construct physical CV cluster-states [18–24], all of which give slightly different states in the finitely squeezed case [31]. Each has important differences that manifest when using them for measurement-based quantum computation. In the ideal, 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. For this reason we choose the canonical method [18]: despite its inefficiency when used in practice [37], this is the most straightforward generalization of the qubit cluster-state preparation [38].

Given a square lattice defined on a graph $G=(\mathcal{V},\mathcal{E})$, we first substitute the N qubits on the vertices with N qumodes (quantum modes—i.e., harmonic oscillators) initialized in the infinitely-squeezed zero-momentum eigenstate $|0{{\rangle }_{p}}$, which is equivalent to the $|+\rangle$ state for qubits. The global state is therefore $|0\rangle _{p}^{\otimes N}$, and this is stabilized by the single mode operators ${{\hat{X}}_{j}}(t)$ in the sense that

$\begin{eqnarray}&&{{\hat{X}}_{j}}(t)|0\rangle _{p}^{\otimes N}=|0\rangle _{p}^{\otimes N},\quad \forall t\in \mathbb{R}\;,\;\forall j\in (1,...,N).\end{eqnarray} \tag{ 21 }$

Alternatively, the same state is nullified, by the set $\{{{\hat{p}}_{j}}\}$ of generators of the stabilizer group since,

$\begin{eqnarray}&&{{\hat{X}}_{j}}(t)|0{{\rangle }_{{{p}_{j}}}}={{{\rm e}}^{-it{{{\hat{p}}}_{j}}}}|0{{\rangle }_{{{p}_{j}}}}=|0{{\rangle }_{{{p}_{j}}}}{{\hat{p}}_{j}}|0{{\rangle }_{{{p}_{j}}}}=0.\end{eqnarray} \tag{ 22 }$

As in the qubit case, the ideal CV cluster-state $|{\rm CS}\rangle$ is the result of the pairwise application of controlled-Z gates ${{{\rm \hat{C}}}_{\hat{Z}(j,k)}}={{{\rm e}}^{i{{{\hat{q}}}_{j}}{{{\hat{q}}}_{k}}}}$ upon all the nearest-neighbor modes $\langle j,k\rangle$ of the initial state $|0\rangle _{p}^{\otimes N}$ as depicted by the graph, explicitly

$\begin{eqnarray}&&\mathop{\prod }\limits_{\langle j,k\rangle }^{N}{{{\rm \hat{C}}}_{\hat{Z}(j,k)}}|0\rangle _{p}^{\otimes N}=|{\rm CS}\rangle .\end{eqnarray} \tag{ 23 }$

How does the application of these gates affect the stabilizer set? Under the ${{{\rm \hat{C}}}_{\hat{Z}(j,k)}}$ evolution, the quadrature operators transform as [39]

$\begin{eqnarray}&&{{{\rm \hat{C}}}_{\hat{Z}(j,k)}}{{\hat{q}}_{j}}{\rm \hat{C}}_{\hat{Z}(j,k)}^{\dagger }={{\hat{q}}_{j}},\end{eqnarray} \tag{ 24 }$

$\begin{eqnarray}&&{{{\rm \hat{C}}}_{\hat{Z}(j,k)}}{{\hat{p}}_{j}}{\rm \hat{C}}_{\hat{Z}(j,k)}^{\dagger }={{\hat{p}}_{j}}-{{\hat{q}}_{k}},\end{eqnarray} \tag{ 25 }$

and thus the initial state stabilizers $\{{{\hat{X}}_{j}}(t)\}$ are changed into the cluster-state stabilizers

$\begin{eqnarray}&&{{\hat{K}}_{{\rm CS}}}={{\hat{X}}_{j}}(t)\mathop{\prod }\limits_{k\in \mathcal{N}(j)}{{\hat{Z}}_{k}}(t),\end{eqnarray} \tag{ 26 }$

which have the same form of the qubit cluster-state stabilizers [40], where $\mathcal{N}(j)$ indicates the nearest neighbors of the jth mode. We can also express these stabilizers by [31]

$\begin{eqnarray}&&{{\hat{X}}_{j}}(t)\mathop{\prod }\limits_{k\in \mathcal{N}(j)}{{\hat{Z}}_{k}}(t)={{{\rm e}}^{-it({{{\hat{p}}}_{j}}-{{\sum }_{k\in \mathcal{N}(j)}}{{{\hat{q}}}_{k}})}},\end{eqnarray} \tag{ 27 }$

equivalent to define the elements of the CV cluster-state nullifier set $\{{{\hat{\eta }}_{j}}\}$ as

$\begin{eqnarray}&&{{\hat{\eta }}_{j}}={{\hat{p}}_{j}}-\mathop{\sum }\limits_{k\in \mathcal{N}(j)}{{\hat{q}}_{k}}\;.\end{eqnarray} \tag{ 28 }$

All the ${{\hat{\eta }}_{j}}$ (and linear combinations of them) commute and are the elements of the algebra that generates the stabilizer group of $|{\rm CS}\rangle$. In analogy with the qubit case, we can construct a Hamiltonian $\hat{H}_{{\rm CS}}^{{\rm ideal}}$ whose ground state is the CV cluster-state $|{\rm CS}\rangle$ by imposing an energy penalty for violating any of the nullifier conditions:

$\begin{eqnarray}&&\hat{H}_{{\rm CS}}^{{\rm ideal}}=\mathop{\sum }\limits_{j=1}^{N}\hat{\eta }_{j}^{2}.\end{eqnarray} \tag{ 29 }$

Since all nullifiers commute and have a continuous spectrum of eigenvalues $(\mathbb{R})$, this Hamiltonian also has a continuous spectrum $[0,\infty )$, and is therefore gapless [19].

Even at the ideal level there is some additional freedom in the CV cluster-state construction procedure. In particular, ideal CV cluster-state graphs can have any non-zero real-valued weight $g\in \mathbb{R}$ associated with each edge. This modifies the strength of the ${{{\rm \hat{C}}}_{\hat{Z}(j,k)}}$ gate represented by that edge: ${{{\rm \hat{C}}}_{\hat{Z}(j,k)}}[g]:={{{\rm e}}^{ig{{{\hat{q}}}_{j}}{{{\hat{q}}}_{k}}}}$. These weights were first introduced in [20] as a way to enable new methods of construction. They have shown themselves to be very important for the computational properties of these states [41] and when considering efficient construction of cluster-states with very large graphs [22, 23, 31]. For the purposes of all future derivations we set g = 1, but when describing possible experimental implementations of the CV codes at the end of this manuscript, we will show how the results obtained with this assumption apply indistinctly to CV cluster-states constructed with non-unit—but still uniform—weight g.

In this section we present the description of the ideal CV surface code. We start introducing the stabilizer description for a general lattice implementation and then look at the special case of a planar square lattice with defined orientations.

The ideal CV cluster-state on a square lattice can be transformed into the corresponding ideal CV surface code by a simple scheme of quadrature measurements [35] inspired by the dynamical mapping of qubit cluster-states to surface codes [42]. In short, start with the CV cluster-state and label vertices by row and column. Then measure in $\hat{p}$ those modes on rows and columns that are both odd and in $\hat{q}$ those that are both even, as in figure 1 for the case of a square lattice. This scheme is equivalent up to translation and/or inversion of the $\hat{p}$, $\hat{q}$ measurements. Since the CV cluster-state is Gaussian, and quadrature measurements are Gaussian operations [31], after the measurements, we are left with a new Gaussian state. The measurements change the form of the stabilizers [31] and therefore after the measurements we are left with a state described by a new set of nullifiers that corresponds exactly to the CV analog of Kitaevʼs surface codes [9].

Zoom In Zoom Out Reset image size

To derive the post-measurement nullifier set, 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 $\Lambda =\{\mathcal{V},\mathcal{E},\mathcal{F}\}$. We assume the graph is oriented and that every face f inherit this orientation. Each quantum mode reside on an edge ${{e}_{j}}\in \mathcal{E}$, with the orientation of any edge determined by $e=[v,v^{\prime} ]$ for the base of the edge starting at vertex v and the head at vertex $v^{\prime}$.

The derivation of the nullifier set for the CV surface code implies knowledge of the nullifiers transformation rules under quadrature measurements. The explicit form of the surface code nullifiers for a square lattice (while allowing for smooth or open boundaries) are obtained by taking linear combination of neighboring nullifiers of the CV cluster-state. A quadrature measurement on a qumode removes it from the cluster [19]. Given a set of exact nullifiers for a Gaussian state [31], one can obtain new nullifiers by a three-step process:

• (1)  
  Given a quadrature measurement ${{\hat{x}}_{j}}$ to be made on mode j, where ${{\hat{x}}_{j}}\in \{\hat{q},\hat{p}\}$, using linear combinations of the original nullifiers, write a new set of nullifiers (remember, they need to commute) such that the canonically conjugate local quadrature ${{\hat{y}}_{j}}$ (where $[{{\hat{x}}_{j}},{{\hat{y}}_{j}}]=\pm i)$ appears in only one nullifier in the new set.
• (2)  
  In each new nullifier, replace ${{\hat{x}}_{j}}$ with the real-valued measurement outcome.
• (3)  
  Eliminate the nullifier that contains ${{\hat{y}}_{j}}$.

We always assume that the outcome of the measurement is zero 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 [33].

Then, in this very general case, the stabilizers for the CV surface code are equal to

$\begin{eqnarray}&&{{\hat{A}}_{v}}(t)=\mathop{\prod }\limits_{e|v\in \partial e}{{\hat{Z}}_{e}}[o(e,v)t],\qquad {{\hat{B}}_{f}}(u)=\mathop{\prod }\limits_{e\in \partial f}{{\hat{X}}_{e}}[-o(e,f)u],\end{eqnarray} \tag{ 30 }$

where the symbols in these expressions denote

$\begin{eqnarray}\begin{array}{rcl} o(e,v) & = & \left\{ \begin{array}{ccccccccccccccc} +1 & {\rm if}\ e\in [v,\cdot ], \\ -1 & {\rm if}\ e\in [\cdot ,v], \\ \end{array} \right. \\ o(e,f) & = & \left\{ \begin{array}{ccccccccccccccc} +1 & {\rm if}\ e\ {\rm is}\;{\rm oriented}\;{\rm the}\;{\rm same}\;{\rm as}\ f, \\ -1 & {\rm otherwise}, \\ \end{array} \right. \\ \end{array}\end{eqnarray} \tag{ 31 }$

and the dot (·) stands for any vertex. By construction, the stabilizers commute:

$\begin{eqnarray}&&[{{\hat{A}}_{v}},{{\hat{A}}_{v^{\prime} }}]=[{{\hat{B}}_{f}},{{\hat{B}}_{f^{\prime} }}]=[{{\hat{A}}_{v}},{{\hat{B}}_{f}}]=0\;.\end{eqnarray} \tag{ 32 }$

The CV surface-code subspace is the $+1$ co-eigenspace of the stabilizers ${{\hat{A}}_{v}}$ and ${{\hat{B}}_{f}}$, similarly to Kitaevʼs surface codes [1].

To simplify the discussion, in the following we specialize to the case of the code graph $\Lambda$ being a toroidal square lattice. We then fix the edge orientations such that at any vertex v, all incident edges point toward v, or all point away from v and the faces inherit equal counterclockwise orientation; see figure 2. Under these assumptions, the symbol $o(e,v)$ is a constant $\pm 1$ for any vertex v, thus it only amounts to a sign flip on t in equation (30), which has no effect on a stabilizerʼs role as such, and we can ignore it. In this case, the CV stabilizers function as the CV analogs of the surface code stabilizers for qubits [9] (up to a $X\leftrightarrow Z$ swap)

$\begin{eqnarray}&&{{\hat{A}}_{v}}(t)\hat{A}_{v}^{q}=\mathop{\prod }\limits_{j\in v}\hat{\sigma }_{j}^{X}\;,\qquad {{\hat{B}}_{f}}(u)\hat{B}_{f}^{q}=\mathop{\prod }\limits_{j\in \partial f}\hat{\sigma }_{j}^{Z}\;,\end{eqnarray} \tag{ 33 }$

and explicitly become

$\begin{eqnarray}&&{{\hat{A}}_{v}}(t)\longrightarrow {{e}^{i{{{\hat{a}}}_{v}}t}},\qquad \;{{\hat{B}}_{f}}(u)\longrightarrow {{e}^{i{{{\hat{b}}}_{f}}u}},\end{eqnarray} \tag{ 34 }$

where the stabilizers generators, i.e the nullifiers, are:

$\begin{eqnarray}&&{{\hat{a}}_{v}}=\mathop{\sum }\limits_{e|v\in \partial e}{{\hat{q}}_{e}}\;,\quad \quad {{\hat{b}}_{f}}=\mathop{\sum }\limits_{e\in \partial f}o(e,f){{\hat{p}}_{e}}\;.\end{eqnarray} \tag{ 35 }$

Zoom In Zoom Out Reset image size

If, on the boundary, one of the edges is missing, then that mode is not included in the nullifiers. It should be clear that the choice of the orientation is purely conventional and does not change the state. Hence, if we were to flip the orientation of every edge, ${{\hat{b}}_{f}}$ simply becomes $-{{\hat{b}}_{f}}$, which is still a nullifier.

The ideal CV surface code is the non-normalizable ground state of the quadratic Hamiltonian [35]

$\begin{eqnarray}&&\hat{H}_{{\rm SC}}^{{\rm ideal}}=\mathop{\sum }\limits_{v\in \mathcal{V}}\hat{a}_{v}^{\dagger }{{\hat{a}}_{v}}+\mathop{\sum }\limits_{f\in \mathcal{F}}\hat{b}_{f}^{\dagger }{{\hat{b}}_{f}}\;.\end{eqnarray} \tag{ 36 }$

In analogy with the ideal cluster-state Hamiltonian, the spurious mode operators ${{\hat{a}}_{v}}$ and ${{\hat{b}}_{f}}$ are actually Hermitian quadrature operators that all commute and do not obey the canonical commutation relations. Consequently, because of

$\begin{eqnarray}\begin{array}{rcl} \left[ {{{\hat{a}}}_{v}},\hat{H}_{{\rm SC}}^{{\rm ideal}} \right] & = & 0\quad \forall v, \\ \left[ {{{\hat{b}}}_{f}},\hat{H}_{{\rm SC}}^{{\rm ideal}} \right] & = & 0\quad \forall f, \\ \end{array}\end{eqnarray} \tag{ 37 }$

an extensive number of gapless modes exist. Therefore this Hamiltonian has a fully continuous spectrum $[0,\infty )$ and is gapless for any number of systems, even on a square lattice with boundary. Anyonic braiding in this model was studied in [35, 43, 44].

Similarly to the previous section, we start considering the CV cluster-state first [18]. To produce the CV cluster-state in the finite-squeezing case by the canonical method one starts from N vacuum modes $|0{{\rangle }^{{{\otimes }_{N}}}}$. The elements of the initial nullifier set are simply the dimensionless single-mode annihilation operators ${{\hat{a}}_{j}}=\frac{1}{\sqrt{2}}({{\hat{q}}_{j}}+i{{\hat{p}}_{j}})$, since

$\begin{eqnarray}&&{{\hat{a}}_{j}}|0{{\rangle }^{\otimes N}}=0,\quad \forall \;j=1,...,N.\end{eqnarray} \tag{ 38 }$

The vacua are then all squeezed by $\hat{S}(s)$: squeezing is performed by the unitary operator $\hat{S}(s)={{e}^{-\frac{i}{2}({\rm log} s)(\hat{q}\hat{p}+\hat{p}\hat{q})}}$, with $s\gt 0$, where ${\rm log} s$ is traditionally known as the squeezing parameter. In the Heisenberg picture, $\hat{S}{{(s)}^{\dagger }}\hat{q}\hat{S}(s)=s\hat{q}$, and $\hat{S}{{(s)}^{\dagger }}\hat{p}\hat{S}(s)=\hat{p}/s$, such that the variance of $\hat{p}$ (of $\hat{q}$) after squeezing is a factor of ${{s}^{-2}}$ (of s²) times its original value. Analogously to the ideal case, the squeezed states are then sent the through pairwise controlled-Z gates ${{{\rm C}}_{{\hat{Z}}}}_{j,k}={{e}^{i{{{\hat{q}}}_{j}}{{{\hat{q}}}_{k}}}}$ in accord with an undirected, unweighted graph with one qumode per vertex. Under the action of the squeezing and the ${{{\rm C}}_{{\hat{Z}}}}_{j,k}$ couplings, the nullifier set for the CV cluster-state is $\{\hat{\eta }_{j}^{s}\}=\{{{{\rm C}}_{{\hat{Z}}}}_{j,k}{{\hat{S}}_{j}}(s){{\hat{a}}_{j}}{{\hat{S}}_{j}}{{(s)}^{\dagger }}{{{\rm C}}_{{\hat{Z}}}}_{j,k}^{\dagger }\}$, with

$\begin{eqnarray}&&\hat{\eta }_{j}^{s}=\frac{s}{\sqrt{2}}\left[ {{s}^{-2}}{{{\hat{q}}}_{j}}+i\left( {{{\hat{p}}}_{j}}-\mathop{\sum }\limits_{k\in \in (j)}{{{\hat{q}}}_{k}} \right) \right].\end{eqnarray} \tag{ 39 }$

These operators satisfy the canonical commutation relations for normal-mode operators, $[\hat{\eta }_{j}^{s},\hat{\eta }_{k}^{s}]=0$ and $[\hat{\eta }_{j}^{s},{{(\hat{\eta }_{k}^{s})}^{\dagger }}]={{\delta }_{j,k}}$, and therefore the CV cluster-state Hamiltonian is [24]

$\begin{eqnarray}&&{{\hat{H}}_{{\rm CS}}}(s)=\mathop{\sum }\limits_{j=1}^{N}\frac{2}{{{s}^{2}}}\left( {{\left( \hat{\eta }_{j}^{s} \right)}^{\dagger }}\hat{\eta }_{j}^{s}+\frac{1}{2} \right).\end{eqnarray} \tag{ 40 }$

For finite s, the system has a gap of $2{{s}^{-2}}$. The prefactor provides for finite energy even in the limit of infinite squeezing.

Using the same measurement pattern shown in figure 1, the finitely squeezed CV cluster-state can be mapped to the finitely squeezed CV surface code. The procedure to obtain the nullifier sets in this case is more elaborate: taking linear combinations of neighboring cluster-state nullifiers—specifically sums of neighboring nullifiers around the $\hat{p}$-measured nodes—and alternating signed cyclic sums around $\hat{q}$-measured modes, one finds the general form of the surface code nullifiers for a generic lattice with possibly incomplete vertices and faces [39]:

$\begin{eqnarray}\begin{array}{rcl} \hat{a}_{v}^{s} & = & \frac{{{s}_{v}}}{\sqrt{2V(v)(1+{{\left( s/{{s}_{v}} \right)}^{2}})}} \\ {} & {} & \times \left[ \mathop{\sum }\limits_{e|v\in \partial e}\left( {{{\hat{q}}}_{e}}+\frac{i}{s_{v}^{2}}{{{\hat{p}}}_{e}} \right)+\frac{{{s}^{2}}}{s_{v}^{2}}\mathop{\sum }\limits_{\frac{v^{\prime} |[v^{\prime} ,v]\in \mathcal{R}}{e|{{v}^{\prime }}\in \partial e\wedge v\notin \partial e}}{{{\hat{q}}}_{e}} \right], \\ \hat{b}_{f}^{s} & = & \frac{s}{\sqrt{2|\partial f|}}\mathop{\sum }\limits_{e\in \partial f}o(e,f)\left( {{{\hat{p}}}_{e}}-\frac{i}{{{s}^{2}}}{{{\hat{q}}}_{e}} \right), \\ \end{array}\end{eqnarray} \tag{ 41 }$

with ${{s}_{v}}=\sqrt{V(v){{s}^{2}}+{{s}^{-2}}}$, V(v) valence of vertex v and $|\partial f|$ boundary size of the lattice face. It is important to realize that now there is a dependence upon the position on the lattice of the vertex or face that we consider. The form of the commutation relations for these nullifiers is not as straightforward as in the previous cases. They are determined by the distance between the vertices and the boundary size of the faces as

$\begin{eqnarray}\begin{array}{rcl} \left[ \hat{a}_{v}^{s},\hat{a}_{v^{\prime} }^{s\,\dagger } \right] & = & \left\{ \begin{array}{ccccccccccccccc} 1 & {\rm if}\ d(v,v^{\prime} )=0, \\ \frac{(s_{v}^{2}+s_{v^{\prime} }^{2}+{{s}^{2}}(V(v)+V(v^{\prime} )))/2{{s}_{v}}{{s}_{v^{\prime} }}}{{{\left[ V(v)V(v^{\prime} )(1+{{\left( s/{{s}_{v}} \right)}^{2}})(1+{{\left( s/{{s}_{v^{\prime} }} \right)}^{2}}) \right]}^{1/2}}} & {\rm if}\ d(v,v^{\prime} )=1, \\ \frac{2{{s}^{2}}/{{s}_{v}}{{s}_{v^{\prime} }}}{{{\left[ V(v)V(v^{\prime} )(1+{{\left( s/{{s}_{v}} \right)}^{2}})(1+{{\left( s/{{s}_{v^{\prime} }} \right)}^{2}}) \right]}^{1/2}}} & {\rm if}\ d(v,v^{\prime} )=\sqrt{2}, \\ \frac{{{s}^{2}}/{{s}_{v}}{{s}_{v^{\prime} }}}{{{\left[ V(v)V(v^{\prime} )(1+{{\left( s/{{s}_{v}} \right)}^{2}})(1+{{\left( s/{{s}_{v^{\prime} }} \right)}^{2}}) \right]}^{1/2}}} & {\rm if}\ d(v,v^{\prime} )=2, \\ 0 & {\rm if}\ d(v,v^{\prime} )\gt 2, \\ \end{array} \right. \\ \left[ \hat{b}_{f}^{s},\hat{b}_{f^{\prime} }^{s\,\dagger } \right] & = & \left\{ \begin{array}{ccccccccccccccc} 1 & {\rm if}\ f=f^{\prime} , \\ \frac{1}{\sqrt{|\partial f||\partial f^{\prime} |}} & {\rm if}\ [f,f^{\prime} ]\in \mathcal{E}, \\ 0 & {\rm otherwise}, \\ \end{array} \right. \\ \left[ \hat{b}_{f}^{s},\hat{b}_{f^{\prime} }^{s} \right] & = & \left[ \hat{a}_{v}^{s},\hat{a}_{v^{\prime} }^{s} \right]=\left[ \hat{a}_{v}^{s},\hat{b}_{f}^{s} \right]=\left[ \hat{a}_{v}^{s},\hat{b}_{f}^{s\,\dagger } \right]=0. \\ \end{array}\end{eqnarray} \tag{ 42 }$

Here $d(v,v^{\prime} )$ is the Euclidean distance between the two vertices v and $v^{\prime}$, where the edge lengths of the graph are unit length. Then the most generic CV physical surface code Hamiltonian is given by

$\begin{eqnarray}&&{{\hat{H}}_{{\rm SC}}}(s)=\mathop{\sum }\limits_{v\in \mathcal{V}}\frac{2V(v)(1+{{s}^{2}}/s_{v}^{2})}{s_{v}^{2}}\hat{a}_{v}^{s\,\dagger }\hat{a}_{v}^{s}+\mathop{\sum }\limits_{f\in \mathcal{F}}\frac{2|\partial f|}{{{s}^{2}}}\hat{b}_{f}^{s\,\dagger }\hat{b}_{f}^{s},\end{eqnarray} \tag{ 43 }$

and the squeezing dependence of the prefactors for the vertex and face parts ensures the Hamiltonian has finite energy in the infinitely squeezed limit:

$\begin{eqnarray}&&\mathop{{\rm lim} }\limits_{s\to \infty }\;{{\hat{H}}_{{\rm SC}}}(s)=\mathop{\sum }\limits_{v\in \mathcal{V}}{{\left( \mathop{\sum }\limits_{e|v\in \partial e}{{{\hat{q}}}_{e}} \right)}^{2}}+\mathop{\sum }\limits_{f}{{\left( \mathop{\sum }\limits_{e\in \partial f}o(e,f){{{\hat{p}}}_{e}} \right)}^{2}}.\end{eqnarray} \tag{ 44 }$

Here we used the fact that for infinite squeezing, each vertex nullifier involves a sum of $\hat{q}$ʼs around that vertex and its four neighboring vertices, and since they all commute, the parent Hamiltonian is simply the squared sum of $\hat{q}$ʼs around each vertex.

In the simpler case of a square lattice with toroidal boundary conditions, the nullifiers become

$\begin{eqnarray}\begin{array}{rcl} {{{\hat{a}}}_{v}} & = & \frac{s^{\prime} }{\sqrt{8}}\left[ \mathop{\sum }\limits_{e|v\in \partial e}\left( {{{\hat{q}}}_{e}}+\frac{i}{s{{^{\prime} }^{2}}}{{{\hat{p}}}_{e}} \right)+\frac{{{s}^{2}}}{s{{^{\prime} }^{2}}}\mathop{\sum }\limits_{\frac{v^{\prime} |[v^{\prime} ,v]\in \mathcal{E}}{e|{{v}^{\prime }}\in \partial e\wedge v\notin \partial e}}{{{\hat{q}}}_{e}} \right], \\ {{{\hat{b}}}_{f}} & = & \frac{s}{\sqrt{8}}\mathop{\sum }\limits_{e\in \partial f}o(e,f)\left( {{{\hat{p}}}_{e}}-\frac{i}{{{s}^{2}}}{{{\hat{q}}}_{e}} \right), \\ \end{array}\end{eqnarray} \tag{ 45 }$

where $s^{\prime} =\sqrt{5{{s}^{2}}+{{s}^{-2}}}$. The state defined by these nullifiers is the CV surface code. Note that now the next-nearest neighbors contribute to the structure of the nullifiers, in contrast to the ideal case, as illustrated in figure 3. Without loss of generality we can consider a square $n\times m$ lattice: Then the nullifiers commutation relations from equation (42) can be rewritten as

$\begin{eqnarray}\begin{array}{rcl} \left[ {{{\hat{a}}}_{v}},\hat{a}_{v^{\prime} }^{\dagger } \right] & = & w\left( d(v,v^{\prime} ) \right), \\ \left[ {{{\hat{b}}}_{f}},\hat{b}_{f^{\prime} }^{\dagger } \right] & = & x(d(f,f^{\prime} )), \\ \left[ {{{\hat{a}}}_{v}},{{{\hat{a}}}_{v^{\prime} }} \right] & = & \left[ {{{\hat{b}}}_{f}},{{{\hat{b}}}_{f^{\prime} }} \right]=\left[ {{{\hat{a}}}_{v}},{{{\hat{b}}}_{f}} \right]=\left[ {{{\hat{a}}}_{v}},\hat{b}_{f}^{\dagger } \right]=0, \\ \end{array}\end{eqnarray} \tag{ 46 }$

where $d(v,v^{\prime} )$ and $d(f,f^{\prime} )$ are the Euclidean distance between vertices and faces, respectively on the unit-edge-length lattice and dual lattice. The functions w and x are given by:

$\begin{eqnarray}\begin{array}{rcl} w(0) & = & 1,\quad w(1)=\frac{(1+8{{s}^{4}})}{4(1+5{{s}^{4}})},\quad w(\sqrt{2})=\frac{{{s}^{4}}}{2(1+5{{s}^{4}})}, \\ w(2) & = & \frac{{{s}^{4}}}{4(1+5{{s}^{4}})},\qquad w(d\gt 2)=0, \\ \end{array}\end{eqnarray} \tag{ 47 }$

and

$\begin{eqnarray}&&x(0)=1,\quad x(1)=\frac{1}{4},\quad x(d\gt 1)=0.\end{eqnarray} \tag{ 48 }$

Using these nullifiers we can construct a Hamiltonian for the physical CV surface code:

$\begin{eqnarray}&&{{\hat{H}}_{{\rm SC}}}(s)=\mathop{\sum }\limits_{v}\frac{8}{s{{^{\prime} }^{2}}}\hat{a}_{v}^{s\,\dagger }\hat{a}_{v}^{s}+\mathop{\sum }\limits_{f}\frac{8}{{{s}^{2}}}\hat{b}_{f}^{s\,\dagger }\hat{b}_{f}^{s}.\end{eqnarray} \tag{ 49 }$

The squeezing dependence of the prefactors is done to ensure the Hamiltonian has finite energy for $s\to \infty$. Unlike the discrete-variable case [6], this Hamiltonian is gapless in the thermodynamic limit. This arises because the nullifiers do not define normal modes. Rather, neighboring nullifiers have non-trivial commutation relations, which allow for low-energy mode excitations.

Zoom In Zoom Out Reset image size

On the $n\times m$ torus, the specifics of the graph are $|\mathcal{E}|=2\;nm$, $|\mathcal{F}|=nm$, and $|\mathcal{V}|=nm$. We first focus on the case where $n\times m$ is odd, so that there are $|\mathcal{E}|$ independent nullifiers spanning the space of all the physical-mode annihilation operators. To diagonalize the Hamiltonian we introduce the normal-mode operators

$\begin{eqnarray}&&{{\hat{c}}_{j}}=\mathop{\sum }\limits_{r=0}^{n-1}\mathop{\sum }\limits_{s=0}^{m-1}\alpha _{r,s}^{(j)}\hat{a}_{{{v}_{r,s}}}^{s}\;,\qquad \;{{\hat{d}}_{j}}=\mathop{\sum }\limits_{i=0}^{n-1}\mathop{\sum }\limits_{l=0}^{m-1}\beta _{i,l}^{(j)}\hat{b}_{{{v}_{i,l}}}^{s},\end{eqnarray} \tag{ 50 }$

where the vertices at the lattice sites have coordinates $\{{{v}_{r,s}}\}$, and the faces at the dual lattice sites have coordinates $\{{{f}_{i,\;l}}\}$. In this base the Hamiltonian is

$\begin{eqnarray}&&{{\hat{H}}_{{\rm SC}}}(s)=\mathop{\sum }\limits_{j}\frac{8{{\omega }_{j}}}{s{{^{\prime} }^{2}}}\hat{c}_{j}^{\dagger }{{\hat{c}}_{j}}+\mathop{\sum }\limits_{j}\frac{8{{\delta }_{j}}}{{{s}^{2}}}\hat{d}_{j}^{\dagger }{{\hat{d}}_{j}}\;.\end{eqnarray} \tag{ 51 }$

To find the normal-mode frequencies, we need to solve the equations

$\begin{eqnarray}&&\left[ {{{\hat{c}}}_{j}},{{{\hat{H}}}_{{\rm SC}}}(s) \right]\equiv \left[ {{{\hat{c}}}_{j}},\mathop{\sum }\limits_{v}\hat{a}_{v}^{\dagger }{{{\hat{a}}}_{v}} \right]={{\omega }_{j}}{{\hat{c}}_{j}},\end{eqnarray} \tag{ 52 }$

where the equality is ensured by $[{{\hat{c}}_{j}},{{\sum }_{f}}\hat{b}_{f}^{s\,\dagger }\hat{b}_{f}^{s}]=0$, and

$\begin{eqnarray}&&\left[ {{{\hat{d}}}_{j}},{{{\hat{H}}}_{{\rm SC}}}(s) \right]\equiv \left[ {{{\hat{d}}}_{j}},\mathop{\sum }\limits_{f}\hat{b}_{f}^{\dagger }{{{\hat{b}}}_{f}} \right]={{\delta }_{j}}{{\hat{d}}_{j}}\;,\end{eqnarray} \tag{ 53 }$

with $[\widehat{{{d}_{j}}},{{\sum }_{v}}\hat{a}_{v}^{s\,\dagger }\hat{a}_{v}^{s}]=0$. If we introduce the generic vertex label state $\{|r,s\rangle \}$ and the face label state $\{|i,l\rangle \}$ the two linear equations above can be vectorized and rewritten in these two basis respectively as

$\begin{eqnarray}&&{{{\bf M}}_{v}}|{{\alpha }^{(j)}}\rangle ={{\omega }_{j}}|{{\alpha }^{(j)}}\rangle ,\quad {{{\bf M}}_{f}}|{{\beta }^{(j)}}\rangle ={{\delta }_{j}}|{{\beta }^{(j)}}\rangle ,\end{eqnarray} \tag{ 54 }$

where $|{{\alpha }^{(j)}}\rangle$ and $|{{\beta }^{(j)}}\rangle$ are the vectorized form of the operators ${{\hat{c}}_{j}}$ and ${{\hat{d}}_{j}}$ that result from equation (50). In the following we define operators for a periodic square lattice that make the problem solvable. First, define the shift operator ${{\hat{X}}_{x}}=\sum _{k=0}^{x-1}|k{{\oplus }_{r}}1\rangle \langle k|$, whose action on the generic vertex label state $|r,s\rangle$ is:

$\begin{eqnarray}\begin{array}{rcl} {{{\hat{X}}}_{n}}|r,s\rangle & = & |r+1,s\rangle ,\quad \hat{X}_{n}^{\dagger }|r,s\rangle =|r-1,s\rangle \;, \\ {{{\hat{X}}}_{m}}|r,s\rangle & = & |r,s+1\rangle ,\quad \hat{X}_{m}^{\dagger }|r,s\rangle =|r,s-1\rangle \;. \\ \end{array}\end{eqnarray} \tag{ 55 }$

The action of the shift operator on the face basis is completely analogous. This allows to rewrite the matrices ${{{\bf M}}_{v}}$ and ${{{\bf M}}_{f}}$, whose elements are non-zero according to the form of the commutation relations in equation (46), in the following elegant way:

$\begin{eqnarray}\begin{array}{rcl} {{{\bf M}}_{v}} & = & {{{\hat{I}}}_{nm}}+w(1)\left[ {{{\hat{I}}}_{n}}\otimes ({{{\hat{X}}}_{m}}+\hat{X}_{m}^{\dagger })+({{{\hat{X}}}_{n}}+\hat{X}_{n}^{\dagger })\otimes {{{\hat{I}}}_{m}} \right] \\ {} & {} & +w(\sqrt{2})\left[ {{{\hat{X}}}_{n}}\otimes {{{\hat{X}}}_{m}}+\hat{X}_{n}^{\dagger }\otimes \hat{X}_{m}^{\dagger }+{{{\hat{X}}}_{n}}\otimes \hat{X}_{m}^{\dagger }+\hat{X}_{n}^{\dagger }\otimes {{{\hat{X}}}_{m}} \right] \\ {} & {} & +w(2)\left[ {{{\hat{I}}}_{n}}\otimes (\hat{X}_{m}^{2}+\hat{X}_{m}^{2\dagger })+(\hat{X}_{n}^{2}+\hat{X}_{n}^{2\dagger })\otimes {{{\hat{I}}}_{m}} \right], \\ {{{\bf M}}_{f}} & = & {{{\hat{I}}}_{nm}}+x(1)\left[ {{{\hat{I}}}_{n}}\otimes ({{{\hat{X}}}_{m}}+\hat{X}_{m}^{\dagger })+({{{\hat{X}}}_{n}}+\hat{X}_{n}^{\dagger })\otimes {{{\hat{I}}}_{m}} \right]. \\ \end{array}\end{eqnarray} \tag{ 56 }$

The linear equations in equation (54) can be solved in the Fourier basis via ${{\hat{F}}_{n}}\otimes {{\hat{F}}_{m}}$, where

$\begin{eqnarray}&&{{\hat{F}}_{x}}=\frac{1}{\sqrt{x}}\mathop{\sum }\limits_{j,k=0}^{r-1}{{e}^{ijk\frac{2\pi }{r}}}|j\rangle \langle k|\;,\end{eqnarray} \tag{ 57 }$

and the nullifiers in the Fourier basis are decomposed as

$\begin{eqnarray}&&\hat{a}_{{{v}_{r,s}}}^{s}=\frac{1}{\sqrt{n}\sqrt{m}}\mathop{\sum }\limits_{{{k}_{x}},{{k}_{y}}}{{e}^{i(r{{k}_{r}}\frac{2\pi }{n}+s{{k}_{s}}\frac{2\pi }{m})}}\widehat{\tilde{a}}_{{{v}_{{{k}_{r}},{{k}_{s}}}}}^{s}.\end{eqnarray} \tag{ 58 }$

In this way, the solutions written in the basis $\{{{\hat{F}}_{n}}|r\rangle \otimes {{\hat{F}}_{m}}|s\rangle \}$ and $\{{{\hat{F}}_{n}}|i\rangle \otimes {{\hat{F}}_{m}}|l\rangle \}$ are

$\begin{eqnarray*}\begin{array}{rcl} \{{{\omega }_{j}}\} & = & \left\{ 1+2w(1)\left[ {\rm cos} \left( \frac{2\pi {{j}_{x}}}{n} \right)+{\rm cos} \left( \frac{2\pi {{j}_{y}}}{m} \right) \right] \right. \\ {} & {} & +2w(\sqrt{2})\left[ {\rm cos} \left( \frac{2\pi {{j}_{x}}}{n}+\frac{2\pi {{j}_{y}}}{m} \right)+{\rm cos} \left( \frac{2\pi {{j}_{x}}}{n}-\frac{2\pi {{j}_{y}}}{m} \right) \right] \\ {} & {} & +\left. 2w(2)\left[ {\rm cos} \left( \frac{4\pi {{j}_{x}}}{n} \right)+{\rm cos} \left( \frac{4\pi {{j}_{y}}}{m} \right) \right] \right\}_{{{j}_{x}}=0,{{j}_{y}}=0}^{n-1,m-1}, \\ \end{array}\end{eqnarray*}$

$\begin{eqnarray}\begin{array}{rcl} \{{{\delta }_{j}}\} & = & \left\{ 1+2x(1)\left[ {\rm cos} \left( \frac{2\pi {{j}_{x}}}{n} \right)+{\rm cos} \left( \frac{2\pi {{j}_{y}}}{m} \right) \right] \right\}_{{{j}_{x}}=0,{{j}_{y}}=0}^{n-1,m-1}, \\ \end{array}\end{eqnarray} \tag{ 59 }$

treating the normal mode index $j=(\;{{j}_{x}},{{j}_{y}})\in {{\mathbb{Z}}_{n}}\times {{\mathbb{Z}}_{m}}$ as a collective index. Then, the (squeezing-dependent) gap energy is the lowest-frequency mode energy:

$\begin{eqnarray}&&\Delta E(s)={{{\rm min} }_{{{j}_{x}},\;{{j}_{y}}}}\left\{ \frac{8{{s}^{2}}{{\omega }_{j}}}{1+5{{s}^{4}}},\frac{8{{\delta }_{j}}}{{{s}^{2}}} \right\}.\end{eqnarray} \tag{ 60 }$

For large system sizes, i.e. $n,m\gg 1$, and choosing freely that $n\leqslant m$, the gap is equal to

$\begin{eqnarray}&&\Delta E(s)\approx \frac{4{{\pi }^{2}}}{{{s}^{2}}{{n}^{2}}}\;,\end{eqnarray} \tag{ 61 }$

and in the thermodynamic limit ${\rm lim} n\to \infty$ the gap goes to zero.

If n and m are even, then not all the face nullifiers are independent. To see this simply bicolor all the lattice faces and assign a plus sign to face operators of one color and a minus sign to faces of the other, then add them to get zero. Thus, the Hamiltonian ${{\hat{H}}_{{\rm SC}}}(s)$ is underconstrained, and there exists an exact gapless zero mode. For a square lattice with planar boundaries (and not toroidal as we discussed so far), there are boundary effects, but these make only a small modification to the gap, which still scales like the inverse of the system size.

Hence, in distinction to the cluster-state Hamiltonian ${{\hat{H}}_{{\rm CS}}}(s)$, the surface code Hamiltonian ${{\hat{H}}_{{\rm SC}}}(s)$ is gapless in the thermodynamic limit, though for infinite squeezing both models are gapless. Note that neither ${{\hat{H}}_{{\rm SC}}}(s)$ nor ${{\hat{H}}_{{\rm CS}}}(s)$ is assumed to be a physical Hamiltonian governing the evolution of the system in question. One may choose to create the initial cluster-state by cooling a physical instance of ${{\hat{H}}_{{\rm CS}}}(s)$, but this is not assumed.

As a consequence of the CV surface code being Gaussian, to evaluate each subsystem entropy, S_X , one only needs the covariance matrix of the post measurement state. To obtain this we find it convenient to use the graphical calculus introduced in section 2.1. Thanks to this formalism it is possible to define the ${{{\bf Z}}_{{\rm CS}}}$ matrix for the CV cluster-state directly, solely making use of the adjacency matrix ${{{\bf A}}_{d}}$ that describes the square-lattice pattern of connections among the modes. Explicitly [31]:

$\begin{eqnarray}&&{{{\bf Z}}_{{\rm CS}}}(s):={{{\bf A}}_{d}}+i{{s}^{-2}}{{{\bf I}}_{N}},\end{eqnarray} \tag{ 64 }$

with squeezing parameter ${\rm log} s$ and ${{{\bf I}}_{N}}$ the $N\times N$ unit matrix. This is illustrated in figure 5. To derive the graph of the finitely squeezed CV cluster-state with open boundary conditions, we start with ${{{\bf Z}}_{{\rm CS}}}$ and perform the measurement scheme as shown in figure 1, using the graph transformation rules explained previously.

Zoom In Zoom Out Reset image size

Then the ${{{\bf Z}}_{{\rm CS}}}$ matrix transforms into the CV finitely squeezed surface code state ${{{\bf Z}}_{{\rm SC}}}$, shown in figure 6. In this case, it is a purely imaginary matrix whose entries are given by

$\begin{eqnarray}&&{{{\bf Z}}_{{\rm SC}}}(s)=i{{{\bf U}}_{{\rm SC}}}(s),\end{eqnarray} \tag{ 65 }$

where

$\begin{eqnarray}&&{{{\bf U}}_{{\rm SC}}}(s)={{s}^{2}}{{{\bf A}}_{{\rm SC}}}+({{s}^{-2}}+2{{s}^{2}}){{{\bf I}}_{N}},\end{eqnarray} \tag{ 66 }$

and ${{{\bf A}}_{{\rm SC}}}$ is the unweighted adjacency matrix of the surface code (without self-loops, see figure 6). Note that the ${{{\bf V}}_{{\rm SC}}}$ component of the ${{{\bf Z}}_{{\rm SC}}}$ matrix, see equation (16), is zero for the surface code. Using the connection between ${\bf Z}$ matrix and covariance matrix from equation (18), we derive the covariance matrix of the physical CV surface code state:

$\begin{eqnarray}{{\Gamma }_{{\rm SC}}}(s)=\frac{1}{2}\left( \begin{array}{ccccccccccccccc} {\bf U}_{{\rm SC}}^{-1}(s) & 0 \\ 0 & {{{\bf U}}_{{\rm SC}}}(s) \\ \end{array} \right).\end{eqnarray} \tag{ 67 }$

Thanks to the formula in equation (9), complete knowledge of the covariance matrix makes calculations of entanglement entropy for different regions of the system straightforward.

Zoom In Zoom Out Reset image size

Using the TEE to characterize the system is meaningful only if the CV surface code state obeys an area law for the entropy. For this assumption to be true, the quadrature correlations among the modes must decay exponentially with lattice separation [45]. From the form of the CV surface code covariance matrix, we see immediately that $\hat{p}$–$\hat{p}$ correlations (determined by ${\bf U}$) have range at most 1, and $\hat{p}$–$\hat{q}$ correlations are zero. On the other hand, the $\hat{q}$-correlations require a more elaborate analysis because they are determined by the matrix ${\bf U}_{{\rm SC}}^{-1}(s)$, which is more complicated. However, it is possible to provide analytical bounds on the $\hat{q}$-correlations for arbitrary squeezing. To prove this, first we show that the spectral range of ${{{\bf U}}_{{\rm SC}}}(s)$, denoted $\sigma ({{{\bf U}}_{{\rm SC}}}(s))$, satisfies $\sigma ({{{\bf U}}_{{\rm SC}}}(s))\subset [a,b]$ where $a={{s}^{-2}}$ and $b={{s}^{2}}(8+{{s}^{-4}})$. For the minimum eigenvalue, note that for finite squeezing the matrix ${{{\bf U}}_{{\rm SC}}}(s)/{{s}^{2}}$ is positive definite but for $s\to \infty$, some $\hat{q}$–$\hat{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 ${{{\bf U}}_{{\rm SC}}}(s)/{{s}^{2}}$ by ${{s}^{-4}}$ so the minimum eigenvalue of ${{{\bf U}}_{{\rm SC}}}(s)$ is $a={{s}^{-2}}$.

To derive the largest eigenvalue of ${{{\bf U}}_{{\rm SC}}}(s)$, observe that for the surface code on a lattice with periodic boundaries, the ${\bf Z}$ graph associated with the adjacency matrix ${{{\bf A}}_{{\rm 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 $\nu =1,...,1$ [46], therefore the maximum eigenvalue of ${{{\bf A}}_{{\rm SC}}}$ is 6. From these considerations it follows that the largest eigenvalue of ${{{\bf U}}_{{\rm SC}}}(s)$ is $b={{s}^{2}}(8+{{s}^{-4}})$. For lattices with open boundary conditions, the eigenvalues of ${{{\bf A}}_{{\rm SC}}}$ are upper bounded by the maximal degree with corrections that fall off with the system size, so the spectrum of ${{{\bf U}}_{{\rm SC}}}(s)$ lies in the interval $[a,b]$.

On a generic $n\times m$ lattice, the $nm\times nm$ matrix ${{{\bf U}}_{{\rm SC}}}(s)$ is block tridiagonal with n identical $m\times n$ matrices ${\bf A}$ on the diagonal and identical ${\bf B}$ on the immediate upper and lower blocks. Now the matrix coordinates (i, j) correspond to Euclidean coordinates $(({{i}_{x}},{{i}_{y}}),({{j}_{x}},{{j}_{y}}))$ on the lattice where $i=m{{i}_{x}}+{{i}_{y}}$ for ${{i}_{x}}\in \{0,...,n-1\}$ and ${{i}_{y}}\in \{0,...,m-1\}$, etc. It is convenient to define a graph distance $d(i,j)={\rm 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 bulk of the ${\bf Z}$ graph is the union of a square graph with a graph having two diagonal edges passing through every other face, 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

$\begin{eqnarray}&&ed(i,j)/\sqrt{2}\leqslant d(i,j)\leqslant ed(i,j),\end{eqnarray} \tag{ 68 }$

where $ed(i,j)=\sqrt{{{\left( {{i}_{x}}-{{j}_{x}} \right)}^{2}}+{{\left( {{i}_{y}}-{{j}_{y}} \right)}^{2}}}$ is the Euclidean distance. The matrix ${\bf A}$ is itself tridiagonal with elements $\alpha =2{{s}^{2}}+{{s}^{-2}}$ on the main diagonal and $\beta ={{s}^{2}}$ on the immediate upper and lower diagonal. The matrix ${\bf B}$ is also tridiagonal with diagonal elements $\delta ={{s}^{2}}$ and immediate upper and lower diagonal elements either equal to zero or δ. A theorem of Demko, Moss and Smith [47, 48] shows that banded matrices of a certain class have inverses with matrix elements that decay exponentially with the distance from the diagonal. Specifically, for matrices ${\bf M}$ of size $N\times N$ and spectral range $\sigma ({\bf M})\subset [a,b]$ with $a\gt 0$ ([47] proposition 5.1):

$\begin{eqnarray}&&{\rm sup} \{|{\bf M}_{i,j}^{-1}|:(i,j)\in {{D}_{n}}({\bf M})\}\leqslant {{C}_{0}}{{q}^{n+1}},\end{eqnarray} \tag{ 69 }$

where the decay sets are

$\begin{eqnarray*}&&{{D}_{n}}({\bf M})=(\{1,\ldots N\}\times \{1,\ldots N\})\backslash {{S}_{p}}({\bf M}),\end{eqnarray*}$

and the support sets are

$\begin{eqnarray*}&&{{S}_{p}}({\bf M})=\mathop{\bigcup }\limits_{k=0}^{p}\{(i,j):{\bf M}_{i,j}^{k}\ne 0\}.\end{eqnarray*}$

Here ${{C}_{0}}=\frac{{{\left( 1+\sqrt{b/a} \right)}^{2}}}{2b}$ and $q=\frac{\sqrt{b/a}-1}{\sqrt{b/a}+1}\lt 1$.

The matrices relevant to our problem are in this class. The matrix power ${\bf U}_{{\rm SC}}^{k}(s)$ is a banded block symmetric matrix with blocks of size $m\times m$ and block band width $2k+1$. Furthermore, each such block is banded with band width $2k+1$. Thus the support set ${{S}_{p}}({{{\bf U}}_{{\rm 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 equation (69) is that for nodes separated in graph distance $d(i,j)\gt 2p+1$ with associated matrix coordinates (i, j) the inverse matrix element ${\bf U}_{{\rm SC}}^{-1}{{(s)}_{i,j}}$ falls off exponentially with graph distance as ${{q}^{(d(i,j)+1)/2}}$.

This implies that the position correlations between nodes i and j separated in graph distance by $d(i,j)$ satisfy

$\begin{eqnarray}&&\langle {{\hat{q}}_{i}}{{\hat{q}}_{j}}\rangle \leqslant C{{e}^{-(d(i,j)+1)/\xi }},\end{eqnarray} \tag{ 70 }$

where the constant is given by

$\begin{eqnarray*}&&C=\frac{{{\left( 1+\sqrt{8{{s}^{4}}+1} \right)}^{2}}}{4(8{{s}^{2}}+{{s}^{-2}})},\end{eqnarray*}$

and the correlation length is

$\begin{eqnarray}&&\xi =\frac{2}{{\rm ln} \left[ \frac{\sqrt{8{{s}^{4}}+1}+1}{\sqrt{8{{s}^{4}}+1}-1} \right]}.\end{eqnarray} \tag{ 71 }$

For Gaussian states all higher order correlations are generated by the linear and quadratic ones. Thus the CV surface code state obeys an area law, and the use of the TEE formulas is appropriate. Numerically we find for squeezing ${\rm log} s=3.2$ the correlation length is $\xi =2.44$ and the scale factor of entropy with area is $\alpha =4.68$. At 5 dB squeezing $\xi =0.33$.

In our simulations we use the ground state of a $36\times 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 equations (62) and (63). See figure 7 for numerical results. Note that even for $s\to 1$, which corresponds to starting (in the CV cluster-state preparation) with vacuum states instead of momentum-squeezed states, the TEE is very small but non-zero. This is not a numerical artifact but a consequence of applying the controlled-Z gates, which introduce additional squeezing [37, 52] to input vacuum states. The values of the TEE calculated using the KP formula (62) and the LW formula (63) are extremely close for relatively large squeezing, and the difference accounts for numerics approximation, as expected. However, on the inset the difference is greater that the absolute values of $S_{{\rm topo}}^{{\rm KP}}$ and $S_{{\rm topo}}^{{\rm LW}}$: this is due to the different formulas used, which result in different values of γ for small squeezing. As a sanity check, we also calculated the TEE for a $36\times 36$ mode CV cluster-state, which proves to be zero for any value of the squeezing. This is in complete analogy with the qubit-based cluster-state, which does not exhibit topological order.

Zoom In Zoom Out Reset image size

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 [53, 54]. For a state with covariance matrix $\Gamma ={{\Gamma }_{q}}\oplus {{\Gamma }_{p}}$, as in our case, the log-negativity of the reduced state with support on a subsystem A is [55]

$\begin{eqnarray}&&\mathcal{N}({{\rho }_{A}})=-\frac{1}{2}\mathop{\sum }\limits_{i=1}^{N}{{{\rm log} }_{2}}[{\rm min} (1,{{\lambda }_{i}}({{\Gamma }_{q}}{{\mu }_{A}}{{\Gamma }_{p}}{{\mu }_{A}}))]\;,\end{eqnarray} \tag{ 72 }$

where ${{\mu }_{A}}={{P}_{{{A}^{\bot }}}}\oplus (-{{P}_{A}})$, with P_X being the projector onto the modes in region X and λ_i being the ith eigenvalue of the matrix argument. Then the TLN is constructed replacing the von Neumann entropy with the logarithmic negativity corresponding to the same region,

$\begin{eqnarray}&&{{\gamma }_{{\rm LN}}}\equiv -({{\mathcal{N}}_{A}}+{{\mathcal{N}}_{B}}+{{\mathcal{N}}_{C}}-{{\mathcal{N}}_{AB}}-{{\mathcal{N}}_{BC}}-{{\mathcal{N}}_{AC}}+{{\mathcal{N}}_{ABC}}).\end{eqnarray} \tag{ 73 }$

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 A. The squeezing-dependent symplectic eigenvalue of the reduced covariance matrix corresponding to the mode used is given by ${{\sigma }_{1}}=\frac{1}{2}{{\left( 1+3{{s}^{4}}+2{{s}^{8}} \right)}^{1/2}}{{\left( 1+3{{s}^{4}} \right)}^{-1/2}}$. The mode entropy grows linearly with an asymptotic slope of ${{{\rm lim} }_{s\to \infty }}\frac{d\gamma (s)}{d({\rm log} s)}=2/{\rm ln} (2)\simeq 2.8854$, which matches the slope of TEE and TLN we find numerically for larger systems.

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={{\mathbb{Z}}_{2}}$). As shown in [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 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({{\rho }_{A}})=(|\partial A|-1){{{\rm log} }_{2}}|G|$, implying $\gamma ={{{\rm log} }_{2}}|G|={{{\rm log} }_{2}}(\mathcal{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 figures A1 (a), (b). For the toric code the ground state would be the GHZ state $(|000\rangle +|111\rangle )/\sqrt{2}$ since both vertices implement the stabilizer $\hat{\sigma }_{1}^{X}\hat{\sigma }_{2}^{X}\hat{\sigma }_{3}^{X}$, and one face enforces $\hat{\sigma }_{1}^{Z}\hat{\sigma }_{2}^{Z}$, while the other face enforces $\hat{\sigma }_{2}^{Z}\hat{\sigma }_{3}^{Z}$. Identifying two qubits on one of the faces with subsystem B, the subsystem entropy is $S({{\rho }_{A}})=S({{\rho }_{B}})=2-\gamma =1$, where 2 comes form the size of the boundary of region B and therefore $\gamma =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 (figure A1(d)). 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 $\{\pm {{\sigma }_{1}}\}$ with

$\begin{eqnarray}&&{{\sigma }_{1}}=\frac{1}{2}\sqrt{\frac{1+3{{s}^{4}}+2{{s}^{8}}}{1+3{{s}^{4}}}},\end{eqnarray} \tag{ A.1 }$

thus, using the formula for the entanglement entropy in equation (9), S(A) is

$\begin{eqnarray}&&S(A)=\left[ \left( {{\sigma }_{1}}+\frac{1}{2} \right){{{\rm log} }_{2}}\left( {{\sigma }_{1}}+\frac{1}{2} \right)-\left( {{\sigma }_{1}}-\frac{1}{2} \right){{{\rm log} }_{2}}\left( {{\sigma }_{1}}-\frac{1}{2} \right) \right],\end{eqnarray} \tag{ A.2 }$

and the upper bound for the TEE can be expressed as

$\begin{eqnarray}&&\gamma \leqslant S(A).\end{eqnarray} \tag{ A.3 }$

Zoom In Zoom Out Reset image size

Here we derive a lower bound for the TMI of the (noisy) CV surface code state introduced in the main text, analyzing the limit $\kappa \to \infty$ (with $\kappa ={\rm coth} \beta {{\epsilon }_{0}}/2$, ${{\epsilon }_{0}}=2/{{s}^{2}}$ gap of the CV physical cluster-state), which corresponds to the strongest possible decrease of the TMI from the TEE for the noise model considered. Recall that for a reduction ${{\rho }_{A}}$ of a pure state ρ, the von Neumann entropy $S({{\rho }_{A}})$ determines the entanglement entropy of the subsystem with respect to its complement. For Gaussian states, one can use the formula

$\begin{eqnarray}&&S({{\rho }_{A}})=\mathop{\sum }\limits_{i=1}^{n_{A}^{\gt }}\left[ \left( {{\sigma }_{i}}+\frac{1}{2} \right){{{\rm log} }_{2}}\left( {{\sigma }_{I}}+\frac{1}{2} \right)-\left( {{\sigma }_{i}}-\frac{1}{2} \right){{{\rm log} }_{2}}\left( {{\sigma }_{I}}-\frac{1}{2} \right) \right],\end{eqnarray} \tag{ A.4 }$

where ${{\{{{\sigma }_{i}}\}}_{A}}$ is the collection of $n_{A}^{\gt }$ symplectic eigenvalues associated to the reduced covariance matrix ${{\Gamma }_{A}}$ of the subsystem ${{\rho }_{A}}$. About the notation, in the following $n_{X}^{\geqslant }={{n}_{X}}$ indicates the total number of symplectic eigenvalues $\geqslant \frac{1}{2}$ associated to a region X, $n_{X}^{\gt }$ indicates the number of symplectic eigenvalues $\gt \frac{1}{2}$, and $n_{X}^{=}$ denotes those $=\;\frac{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 (equation (40) 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 $\kappa ={\rm coth} (\beta /{{s}^{2}})$. To detect the topological order of the resulting mixed CV surface-code state, we use the TMI:

$\begin{eqnarray}&&{{\gamma }_{MI}}=-\frac{1}{2}({{I}_{A}}+{{I}_{B}}+{{I}_{C}}-{{I}_{AB}}-{{I}_{BC}}-{{I}_{AC}}+{{I}_{ABC}}).\end{eqnarray} \tag{ A.5 }$

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 (Γ₀) times κ, mathematically

$\begin{eqnarray}&&\Gamma =\kappa {{\Gamma }_{0}}.\end{eqnarray} \tag{ A.6 }$

As a consequence, the symplectic eigenvalues of $\Gamma$ (or any reduced section of it) are simply given by the'pure' symplectic eigenvalues multiplied by the overall κ factor, $\{{{\sigma }_{i}}\}=\{\kappa \sigma _{i}^{0}\}$. This simple transformation of $\Gamma$ is a special case that only arises due to the fact that all normal modes of the CV cluster-state Hamiltonian are identical, resulting in equal symplectic eigenvalues $\kappa /2$.

For a large value of κ, we find the following asymptotic expression for each eigenvalue contribution to equation (A.4):

$\begin{eqnarray*}&&\left( \kappa \sigma +\frac{1}{2} \right){{{\rm log} }_{2}}\left( \kappa \sigma +\frac{1}{2} \right)-\left( \kappa \sigma -\frac{1}{2} \right){{{\rm log} }_{2}}\left( \kappa \sigma -\frac{1}{2} \right)\approx {{{\rm log} }_{2}}(e\kappa \sigma ).\end{eqnarray*}$

We can use this to show the behavior of the TMI as $\kappa \to \infty$. To start, consider the first term of the TMI formula, equation (A.5):

$\begin{eqnarray}&&{{I}_{A}}={{S}_{A}}+{{S}_{BCD}}-{{S}_{ABCD}},\end{eqnarray} \tag{ A.7 }$

where the regions used are shown in (A2 ). If the total state ABCD has N modes, then, for region A, we have n_A modes and for its complement BCD we have n_BCD modes, such that ${{n}_{A}}+{{n}_{BCD}}=N$. Hence, the von Neumann entropy for A is given by

$\begin{eqnarray}&&{{S}_{A}}=\mathop{\sum }\limits_{i=1}^{{{n}_{A}}}\left[ \left( \kappa \sigma _{i}^{A}+\frac{1}{2} \right){{{\rm log} }_{2}}\left( \kappa \sigma _{i}^{A}+\frac{1}{2} \right) \right.-\left. \left( \kappa \sigma _{i}^{A}-\frac{1}{2} \right){{{\rm log} }_{2}}\left( \kappa \sigma _{i}^{A}-\frac{1}{2} \right) \right].\end{eqnarray} \tag{ A.8 }$

In the limit of very high temperature, corresponding to $\kappa \to \infty$, we can use the approximation (A.7) and rewrite the von Neumann entropy S_A as

$\begin{eqnarray}&&\mathop{{\rm lim} }\limits_{\kappa \to \infty }{{S}_{A}}\approx S_{A}^{l}=\mathop{\sum }\limits_{i=1}^{{{n}_{A}}}{{{\rm log} }_{2}}(e\kappa \sigma _{i}^{A})=\mathop{\sum }\limits_{i=1}^{{{n}_{A}}}{{{\rm log} }_{2}}(e\sigma _{i}^{A})+{{n}_{A}}{{{\rm log} }_{2}}\kappa .\end{eqnarray} \tag{ A.9 }$

Divide now the n_A symplectic eigenvalues into the two sets $n_{A}^{\gt }$ and $n_{A}^{=}$, such that $n_{A}^{\gt }+n_{A}^{=}{{n}_{A}}$. Hence, we can rewrite S_A ^l as:

$\begin{eqnarray}&&S_{A}^{l}=\mathop{\sum }\limits_{i=1}^{n_{A}^{\gt }}{{{\rm log} }_{2}}(e\sigma _{i}^{A})+\mathop{\sum }\limits_{i=1}^{n_{A}^{=}}{{{\rm log} }_{2}}\left( \frac{e}{2} \right)+{{n}_{A}}{{{\rm log} }_{2}}\kappa \;.\end{eqnarray} \tag{ A.10 }$

The same argument can be repeated for the region BCD and find that

$\begin{eqnarray}&&S_{BCD}^{l}=\mathop{\sum }\limits_{i=1}^{n_{BCD}^{\gt }}{{{\rm log} }_{2}}(e\sigma _{i}^{BCD})+\mathop{\sum }\limits_{i=1}^{n_{BCD}^{=}}{{{\rm log} }_{2}}\left( \frac{e}{2} \right)+{{n}_{BCD}}{{{\rm log} }_{2}}\kappa \;.\end{eqnarray} \tag{ A.11 }$

The N symplectic eigenvalues of ABCD are all equal to $\kappa /2$, consequently:

$\begin{eqnarray}&&S_{ABCD}^{l}=\mathop{\sum }\limits_{i=1}^{N}{{{\rm log} }_{2}}(\kappa \frac{e}{2})=\mathop{\sum }\limits_{i=1}^{N}{{{\rm log} }_{2}}\kappa +\mathop{\sum }\limits_{i=1}^{N}{{{\rm log} }_{2}}\frac{e}{2}=N{{{\rm log} }_{2}}\kappa +N{{{\rm log} }_{2}}\left( \frac{e}{2} \right),\end{eqnarray} \tag{ A.12 }$

and the value for the mutual information I_A ^l in the $\kappa \to \infty$ limit is given by

$\begin{eqnarray}\begin{array}{rcl} I_{A}^{l} & = & S_{A}^{l}+S_{BCD}^{l}-S_{ABCD}^{l} \\ {} & = & \mathop{\sum }\limits_{i=1}^{n_{A}^{\gt }}{{{\rm log} }_{2}}(e\sigma _{i}^{A})+\mathop{\sum }\limits_{i=1}^{n_{BCD}^{\gt }}{{{\rm log} }_{2}}(e\sigma _{i}^{BCD})+(n_{A}^{=}+n_{BCD}^{=}-N){{{\rm log} }_{2}}\left( \frac{e}{2} \right). \\ \end{array}\end{eqnarray} \tag{ A.13 }$

Notice that the κ-contributions cancel out exactly and, asymptotically, I_A ^l is not a function of κ. Using $n_{A}^{=}+n_{BCD}^{=}N-n_{A}^{\gt }-n_{BCD}^{\gt }$ and the area law behavior for the entropy, i.e., for regions sufficiently big, $n_{A}^{\gt }=n_{BCD}^{\gt }$ (although this does not mean that the sets $\{\sigma _{A}^{i}\}$ and $\{\sigma _{i}^{BCD}\}$ are the same), we can rewrite the mutual information as

$\begin{eqnarray}&&I_{A}^{l}=\mathop{\sum }\limits_{i=1}^{n_{A}^{\gt }}{{{\rm log} }_{2}}(e\sigma _{i}^{A})+\mathop{\sum }\limits_{i=1}^{n_{BCD}^{\gt }}{{{\rm log} }_{2}}(e\sigma _{i}^{BCD})-2n_{A}^{\gt }{{{\rm log} }_{2}}\left( \frac{e}{2} \right).\end{eqnarray} \tag{ A.14 }$

The same argument applies for each other contribution to the TMI, for example, for region C we have:

$\begin{eqnarray}&&I_{C}^{l}=\mathop{\sum }\limits_{i=1}^{n_{C}^{\gt }}{{{\rm log} }_{2}}(e\sigma _{i}^{C})+\mathop{\sum }\limits_{i=1}^{n_{ABD}^{\gt }}{{{\rm log} }_{2}}(e\sigma _{i}^{ABD})-2n_{C}^{\gt }{{{\rm log} }_{2}}\left( \frac{e}{2} \right).\end{eqnarray} \tag{ A.15 }$

When substituting these expressions for the mutual information into the TMI formula in equation (A.5), all the elements $2n_{X}^{\gt }{{{\rm log} }_{2}}(\frac{e}{2})$, which depend on the size of the region boundary, sum to zero for the same reason why all the area-dependent elements cancel out in the KP argument [7] for the TEE.

Zoom In Zoom Out Reset image size

Consequently, the lower limit for the TMI is simply given by

$\begin{eqnarray}&&\gamma _{MI}^{l}=-\frac{1}{2}\mathop{\sum }\limits_{X}\zeta (X)\mathop{\sum }\limits_{i=1}^{n_{X}^{\gt }}{{{\rm log} }_{2}}(e\sigma _{i}^{X}),\end{eqnarray} \tag{ A.16 }$

where X runs over all the possible combinations of regions, and the function $\zeta (X)$ is defined as

$\begin{eqnarray}\zeta (X)=\left\{ \begin{array}{ccccccccccccccc} +1 & {\rm if}\ X\in \{A,B,C,D,ABC,ABD,ACD,BCD\} \\ -1 & {\rm if}\ X\in \{AB,AC,AD,BC,BD,CD\} \\ \end{array} \right.\end{eqnarray} \tag{ A.17 }$

for a partitioning of the system as in figure A2. This formula confirms our numerics and prove that that the value of the TMI, for the particular case of the CV surface code state, is not affected by κ in the limit of thermal input states at high-temperature.