Entanglement in gapless resonating valence bond states

We study resonating-valence-bond (RVB) states on the square lattice of spins and of dimers, as well as SU(N)-invariant states that interpolate between the two. These states are ground states of gapless models, although the SU(2)-invariant spin RVB state is also believed to be a gapped liquid in its spinful sector. We show that the gapless behavior in spin and dimer RVB states is qualitatively similar by studying the R\'enyi entropy for splitting a torus into two cylinders, We compute this exactly for dimers, showing it behaves similarly to the familiar one-dimensional log term, although not identically. We extend the exact computation to an effective theory believed to interpolate among these states. By numerical calculations for the SU(2) RVB state and its SU(N)-invariant generalizations, we provide further support for this belief. We also show how the entanglement entropy behaves qualitatively differently for different values of the R\'enyi index $n$, with large values of $n$ proving a more sensitive probe here, by virtue of exhibiting a striking even/odd effect.


Introduction
A vast amount of recent theoretical and experimental effort has been devoted to the study of and search for spin liquids. A spin-liquid phase typically exhibits exotic behavior such as quasiparticles with fractionalized quantum numbers, or spin-charge separation. Such behavior is manifestly non-perturbative, and so requires advanced theoretical tools to analyze.
A vast amount of recent theoretical effort has been devoted to understanding how the entanglement properties of a particular state illuminate its physical properties. Since a major lesson of condensed matter physics in recent decades was that a phase can often be characterized by a model state, a feature of studying entanglement is that states can be studied directly, without needing to know a specific Hamiltonian. Since studying entanglement properties of a model state is often both analytically and numerically more tractable than analyzing a model Hamiltonian, much recent progress has been made in condensed matter physics from the former.
The behavior of quantities such as the Rényi entropy thus can be used to explore spinliquid states directly. The quantum spin liquids best characterized theoretically are gapped, where a gauge symmetry typically emerges. An important and useful concept applicable here is topological order, which relates the emergent gauge symmetry to a topological degeneracy, and the related concept of topological entropy [1][2][3]. However, many experimental systems currently identified as good candidates for spin liquids appear to be gapless. Compared to their gapped counterparts, these phases are less amenable to simple unifying attributes such as topological order. Although it is likely that other quantities, related perhaps to correlations or manifestations of long-range entanglement, can be used to characterize gapless phases, few viable examples of models-particularly those tractable by large-scale numerical simulation-are known to contain gapless spin liquids.
We study a particular family of states, resonating-valence-bond (RVB) states on the square lattice, that are gapless and exhibit spin-liquid-type behavior [4][5][6]. Our main (but not only) tool is to analyze a subleading term in the Rényi entropy particularly useful for characterizing gapless systems [7]. The Rényi entanglement entropies S n are defined as where ρ A is the reduced density matrix for a region A, and the von Neumann entropy S 1 is defined as the limit n → 1. The term we analyze describes the entanglement between two regions formed by cutting a torus into two cylinders. This geometry is useful for studying gapless systems because the length of the boundaries separating the cylinders is independent of the area of the cylinders. Because of the long-range correlations, the entanglement entropy of gapless systems will be sensitive to varying the size of the cylinders. This dependence was studied for three completely different gapless systems in [7], and argued to have universal behavior. One of the results of this paper is to provide additional evidence in support of this assertion.
One model with a gapless RVB ground state, the square-lattice quantum dimer model, is quite amenable to an analytic treatment. For example, all equal-time correlators diagonal in the dimer basis are those of a classical dimer model [4]. These can be computed exactly directly in the dimer model by using Pfaffian techniques [8][9][10], or in the continuum limit by utilizing a two-dimensional (2D) classical free-boson field theory [11,12]. Such exact computations have been done for entanglement quantities as well [13][14][15][16][17][18]. In this paper, we extend these methods to compute the two-cylinder entanglement entropy for the quantum dimer model exactly in the scaling limit. For Rényi index n greater than a specific (non-universal) value n c , it is simply given by a ratio of cylindrical partition functions. From this we find for example that the striking even-odd effect observed numerically in [7] for the SU (2)-invariant RVB state was not a finitesize artifact; it also occurs for the quantum dimer model, and our computation illuminates its origin.
One lesson from our results is that universal entanglement properties can be qualitatively different for values of the Rényi index n. For n > n c , the two-cylinder entanglement entropy depends strongly on whether the cylinders are of odd or even length. This is not shocking from the point of view of the quantum dimer state, because the correlators exhibit similar behavior. However, at the von Neumann point (n → 1), no such even-odd effect is possible. Indeed, the shape of the entanglement entropy has to be a concave-down function of y, as required by the strong subadditivity property [19]. Other values of n are not constrained by subadditivity, and we will derive this even/odd effect explicitly for the dimer RVB state. Thus only for n > n c is the Rényi entropy sensitive enough to include this effect. This effect, whereby a Rényi entropy may exhibit qualitatively different behavior for some higher n, is reminiscent of the manifestation of finite-temperature criticality in Rényi entropies for n > 1, but not n = 1 [20].
at the sites of some lattices. For the N -site square lattice, the SU (2) RVB state is where so each spin i on one sublattice is in a singlet with one of its nearest neighbors j i on the other sublattice. A VB between neighbors i and j can be labeled by a dimer on this link between i and j, and so each VB state can be viewed as a dimer configuration C with exactly one dimer touching each site. Our computations probe a state directly and so do not require a Hamiltonian. However, it is worth noting that this SU (2) RVB state is the exact ground state of a quantum spin Hamiltonian with local interactions [29]. This RVB state breaks no symmetry, and so is a natural candidate for a state exhibiting spin-liquid behavior. For example, if the 2D square lattice is periodic (i.e. is a torus) and has an even number of sites in both directions, there is a winding number for each direction. These are defined simply by counting the dimers around a closed path going straight around a cycle, weighting each dimer by +1 for an even number of steps along the path, −1 for an odd number of steps, with no contribution from empty links. An RVB state can be defined for each value of the winding numbers, and each will be the ground state of a local Hamiltonian. The number of ground states will then depend on the number of non-trivial cycles, a basic characteristic of topological order.
For the square lattice, extensive numerical analysis indicates that dimer-dimer correlators in the RVB state decay algebraically [5,6]. A theorem of Hastings guarantees that when any local operators have an algebraically decaying correlator, any local Hamiltonian with this as a ground state must be gapless [30]. Thus the square-lattice RVB state could only describe a gapless system. However, it has long been known, both from numerics and strong analytic arguments, that spin-spin correlators are exponentially decaying [28]. Coupled with the fact that no local order parameter has been found [5,6], it seems very possible that there is a spin gap, and that the state effectively behaves as a spin liquid. In this paper we do not address this question directly. Instead we focus on the gapless behavior of nearest-neighbor RVB states on the square lattice.

The dimer RVB state
One key tool is to study the dimer RVB state [4]. Here one essentially forgets the underlying spins: the only degrees of freedom are the nearest-neighbor dimers. Each configuration of dimers is a linearly independent basis element of the Hilbert space of the quantum dimer model, and these dimer 'states' span the Hilbert space. As above, a dimer configuration C has exactly one dimer touching each site of the lattice, i.e. the dimers are close packed and hard-core. The RVB state for dimers is the equal-amplitude sum over all these states: Different dimer configurations are defined to be orthogonal: The 'RK' Hamiltonian [4] is a simple local Hamiltonian with this dimer RVB state as its ground state.
One very convenient feature of the dimer RVB state is that any correlator diagonal in the dimer basis can be computed exactly, because they are exactly those of the 2D classical dimer model. This follows immediately from the definitions of the state and the inner product. For example, the normalization of the dimer RVB state is the partition function of the classical dimer model, the number of distinct dimer configurations: which behaves asymptotically on the square lattice as ∼ (1.791 62 · · ·) N [8]. A classic result of classical statistical mechanics is that any correlator in the classical dimer model on any planar lattice can be expressed as a Pfaffian [8,9]. This is equivalent to saying that the model can be rewritten in terms of free fermions [31]. For the dimer-dimer correlation function, the asymptotic behavior can easily be obtained using this expression. For the square lattice long ago this was shown to be algebraically decaying: with exponent α = 2 [10]. The coefficient C depends on both the direction and the parity of the number of steps along any path connecting the two dimers; in some special cases it vanishes and so subleading terms dominate. In fermionic language, the dimer creation operator is a fermion bilinear, and so is of dimension 1. Algebraic decay occurs only for bipartite lattices; for non-bipartite lattices, correlators are exponentially decaying, a fact crucial in showing that the triangular-lattice quantum dimer model is in a gapped RVB phase, i.e. has topological order [32]. It is also possible to generalize this construction, with 'enforced' orthogonality (5), to other models. For example starting from a classical six-vertex model, we get a quantum sixvertex wavefunction [33]. This wavefunction also has algebraically decaying correlations, with a continuously varying exponent α as a function of the vertex weights on its critical line. The classical model, although in general more complicated than the dimer model, is nevertheless integrable. It is also constructed with a similar set of constraints, the fully packed and hardcore constraints being replaced by the ice rule. All the analytical results we will obtain for the entropy of quantum dimers can be straightforwardly generalized to the quantum six-vertex case.
The inner product is an important difference between the dimer and SU (2) RVB states. In the RVB state (2,3), the inner product is the usual one for spins: each different spin state in the S z diagonal basis is orthogonal to the others. As is easily shown [22], this means that different VB states in the spin model are not orthogonal. There is a convenient loop-gas description to represent their inner product. One places the two dimer graphs on top of each other, obtaining the 'transition' graph. Because of the requirement that there be exactly one dimer per site, this means that the transition graph is comprised of closed loops. The shortest loops are of length two, when both graphs have a dimer on the same link. Using the spin inner product then gives [34], n l being the number of loops in the transition graph. Despite this overlap between dimer configurations making analytic computations much more difficult, its presence is convenient in the algorithm used for computing the second Rényi entropy using quantum Monte Carlo simulations [35], since 'swapped' basis configurations (described in section 6.1) will always contribute a non-zero expectation value, even with a naive QMC updating algorithm. Despite this important difference between the dimer and SU (2) RVB states, their obvious similarity makes it plausible that the dimer-dimer correlator will behave similarly in the two cases. Further evidence for this fact is that rewritten in loop language, the spin-spin correlator corresponds to the probability that the two spins are on the same loop. This decays exponentially for loop models with weight per loop greater than 2 [36], and the inner product (8) amounts to a weight per loop of 4 [28]. This exponential decay arises because the classical partition function is dominated by configurations with many short loops. Thus correlators in the SU (2) RVB state are dominated by short-loop configurations, i.e. the dimers. Whereas the presence of longer loops renormalizes the exponent α, it seems very likely that dimer-dimer correlators in the RVB state will decay algebraically like they do for the dimer RVB state.
This expectation was convincingly demonstrated numerically [5,6]. The dimer-dimer correlator for dimers was found to be of the form (7) with exponent close to α ≈ 1.2. This is an even weaker falloff than in the dimer RVB state, so the presence of the loops in the inner product of the SU (2) RVB state only makes this correlator longer-ranged. By the theorem of [30], algebraic decay of local operators in a 2D ground state requires that any local Hamiltonian with this ground state be gapless. Thus not only is the quantum dimer model with RK Hamiltonian gapless, the spin Hamiltonian of [29], whose exact ground state is the SU (2) RVB state, must be as well.

The SU (N ) RVB state
It is possible to interpolate between the dimer and SU (2) RVB states, either by dressing the lattice itself [37], or by modifying the SU (2) symmetry to SU (N ). In the SU (N ) RVB state, the 'spins' have N = 2S + 1 components, and an SU (N ) singlet between sites i and j is given by When i and j are nearest neighbors, such a singlet can be viewed as a dimer occupying the link between i and j, as in the SU (2) case. The nearest-neighbor RVB wavefunction is then defined as an equally weighted superposition of these VB states (2) as before, with For SU (N ) VB states, the orthogonality relation becomes where n l is again the number of loops in the transition graph formed by superimposing the two dimer configurations C and C . Notice that n l can be at most N /2, and this happens only when all the loops are trivially flat, i.e. when the two dimer configurations are identical. We have however in the limit S → ∞ (or equivalently N → ∞). As a result, any correlation calculated in the limit N → ∞ will be identical to that in the corresponding dimer model. It is however important to stress that this argument strictly only holds at the level of correlation functions. We will present in this paper a variety of results for the SU (N ) case, and show that there is no sign of a phase transition as N is varied, thus giving further evidence for the similarity in behavior between the SU (2) and dimer RVB states. We will see, for example, that for N = 2, 3, 4, 5 the dimer-dimer correlator is algebraically decaying, with exponent depending on N .
Given the changing exponents, the SU (N ) RVB states are obviously not identical to those for dimers, even in the scaling limit. In fact, we will present strong evidence that equal-time correlators for the dimer and SU (N ) RVB states can all be written in terms of a free boson, or Coulomb gas.

The Coulomb-gas approach to classical dimers
A well established and very useful fact about the classical dimer model is that long-distance properties of the dimers are described by a free bosonic field theory, or Coulomb gas. This can be seen either by taking the continuum limit of the fermionic description (see e.g. [38]) and then bosonizing the fermions, or by rewriting the dimers in terms of a discrete 'height' degree of freedom [12,39,40]. We follow the latter approach, because a careful treatment of the boundary conditions is necessary to obtain our results for the two-cylinder Rényi entropy.
The first step in the mapping of classical dimers on to the Coulomb gas is to associate an integer with each face of the lattice, as illustrated in figure 1. These integers are defined by the following construction. Set value on one face (say the lower leftmost one) to be zero. Then, turning counterclockwise around a site of the even (resp. odd) sublattice represented by black squares (resp. white squares), the integer changes by −3 (resp. +3) if it crosses a dimer, +1 (resp. −1) otherwise. Since by construction these integers are different on adjacent sites, in the bulk one defines a smoother height variable h(i, j) on the sites of the lattice, by averaging the values of the integers on four faces surrounding each point (and then for convenience's sake multiplying by π/2). Then the states where the dimers are arranged in columns correspond to a constant value of the height.
An important subtlety in this mapping is the boundary conditions, which will play a crucial role in understanding the even/odd effect in the two-cylinder Rényi entropy. With periodic boundary conditions and an even number of sites around this cycle (as in the left and right sides of figure 1), note that this means that heights can jump by 2π times an integer on from one side to the other. With open boundary conditions, as we have at the ends of the cylinder (the top and bottom of the figure), the natural definition of the boundary height comes from averaging the integers along the boundary. It is then easy to check that for a cylinder of even L x circumference and length L y , the height difference between top and bottom satisfies w ∈ Z, L y even, Even/odd effects have long been studied in the dimer model [41][42][43], but to our knowledge they have never been interpreted in the context of the height mapping. The standard Coulomb-gas arguments [36] then imply upon coarse-graining the height variable h(i, j) into a continuous field, the scaling limit of the classical dimer model is described by a free scalar field h(x, y) with Euclidean action The coupling κ is often known as the 'stiffness'; we show below that for the dimer case it is set to κ D = 1/2. This free-boson field theory is one of the simplest examples of a 2D conformal field theory (CFT) [44]. It is ubiquitous in 2D statistical models, as well as in condensed matter, where it describes Luttinger liquids. In order to be consistent with periodic boundary conditions, field values shifted by 2π must be identified; in conventional language it is said that the field h is compactified on a circle of radius 1: h ∼ h + 2π. The dimers are viewed as elementary electric charges, i.e. are created by the field cos(h) (the simplest function of h that satisfies the compactification h ∼ h + 2π ). The standard computation [36] gives for the dimer-dimer correlator decay exponent: This yields the stiffness for the dimer model κ D = 1/2. Elementary magnetic charges can also be identified with monomers, and the corresponding exponent describing their two-point function is given by β = κ, so β D = 1/2. Once the stiffness is fixed to the appropriate value, equation (14) describes all the universal properties of the classical dimer model. By virtue of the correspondence between correlators in the classical model and those in the dimer RVB state, the same exponents apply to the correlators in the latter as well. The same Coulomb-gas techniques can be used on the six-vertex model too, yielding a value of κ continuously varying with the couplings [36]. In this framework the quantum dimer and quantum six-vertex models can also be written in terms of a free scalar field in the 2 + 1-dimensional 'quantum Lifshitz' model [12,33]. This is a field-theory version of an RK Hamiltonian, where the equal-time correlators of the 2D quantum model are those of a 2D CFT. This model is critical but not Lorentz-invariant, with dynamical critical exponent z = 2.
Because of the arguments above, it is reasonable to hope that all the equal-time correlators in the SU (N ) RVB state for any N can be described by this Coulomb gas, albeit with stiffness depending on N . Substantial evidence for this was provided in the SU (2) case in [6], where several universal quantities, including the dimer-dimer and monomer-monomer decay exponents, were found numerically. With appropriate identification of the operators, each yields an independent value for the stiffness. The numerics gave approximately the same value of the stiffness κ 2 ≈ 0.83 for each, providing strong evidence that the Coulomb-gas description applies to correlators of spin singlets in the SU (2) RVB wavefunction.

Two-cylinder Rényi entropy from partition functions
In this section, we discuss the shape-dependent subleading contribution 6 to the entanglement entropy of 2D gapless systems discussed in [7]. For Rényi index n larger than a certain critical value n c , we derive it exactly for the scaling limit of the square-lattice dimer RVB state, the quantum six-vertex state, and the entire quantum critical line of the quantum Lifshitz theory. Thus we expect it to apply to the scaling limit of the SU (N ) RVB state as well, with an appropriate parameterization of κ in terms of N .
In particular, we show that when a torus is split into two cylinders labeled A and B, the Rényi entropy describing the entanglement between the two cylinders can be rewritten in terms of conformal field theory partition functions Z as: Below we detail the precise definition of the CFT partition functions, but the important fact is that for the cases at hand, these are known. Thus this expression gives an explicit formula for this entropy, exact in the scaling limit. The value of n c depends on the specifics of the model, but n c < 2 in all the examples we will consider, so that equation (16) applies to the second Rényi entropy to be studied numerically in later sections. Although our derivation applies only to the quantum Lifshitz ground state (i.e. RVB-type states that can be written in terms of a 2D classical free boson), we believe it possible that an analogous formula will apply more generally, in particular to any theory with an RK-type Hamiltonian.

Entanglement entropy as a Shannon entropy
We will now detail how to derive equation (16) in the simplest case of quantum dimers on the square lattice. This requires a slight generalization of two ingredients: the mapping to a classical Shannon entropy of [15] (see also [45]) and the boundary phase transition argument of [18] (see section 3.2). Although this derivation does not apply to the SU (N ) RVB case, we will later explain why we think equation (16) should still apply. The orthogonality of the different dimer configurations in the dimer RVB state allows for huge technical simplifications. As was shown in [15], the entanglement entropy can be expressed as a classical Shannon entropy. We give here only a brief summary of the results, and refer to [15] for more details.
Let us cut our torus into two subsystems A and B, as is shown for example in figure 2. We label the dimer configurations along the boundaries as |σ and |µ . The main difference with [15] is that there are two boundaries here, but the same arguments apply. The von Neumann entropy can be recast as a classical Shannon entropy where the p σ,µ are the probabilities of a given boundary configuration |σ ,|µ between A and B. The p σ,µ are precisely the eigenvalues of the reduced density matrix. Said differently, the Schmidt decomposition of the ground state reads where |ψ A σ,µ is a superposition of all dimer states in A compatible with the boundary conditions σ and µ, and the same goes for B. All these vectors are mutually orthogonal: This mapping therefore allows the Rényi entanglement entropy to be rewritten as We name this quantity the Rényi-Shannon entropy in the following. The probabilities can be conveniently rewritten in terms of dimer partition functions as where Z σ,µ is the number of dimer configurations compatible with the boundary configuration |σ , |µ , while Z is the number of dimer coverings on the whole torus.
It is important to emphasize the main reasons why such a remarkable quantum-classical correspondence holds. While (18) remains true for all the states we consider in this paper, the  orthogonality (19) is guaranteed only if three other conditions are satisfied: (i) The orthogonality of quantum states corresponding to different classical configurations.
(ii) A certain type of local constraint [15], which for Z 2 variables on the bonds of the squarelattice (e.g. dimer occupancies) amounts to: if three of the variables around each site are specified, then the value of the fourth is uniquely determined. (iii) Interactions are only between bond variables sharing a common site.
Quantum dimers obviously satisfy all these requirements. The same goes for the quantum sixvertex state, (ii) being guaranteed by the ice rule. The SU (N ) RVB states satisfy (ii) and (iii) but not (i) (see (8)). An example of a Coulomb-gas wavefunction that satisfies (i) and (ii) but not (iii) is a dimer wavefunction built out of an interacting dimer model [23,24]. Let us finally mention that when equation (19) does not hold, it is sometimes possible to nevertheless obtain the Schmidt eigenvalues by numerically diagonalizing the matrix whose elements are the ψ σ,µ |ψ σ ,µ [46,47].
Returning to the dimer RVB state, these expressions make it possible to use the computer to find essentially exact values for the Rényi entropy. Using free-fermion techniques described in appendix A, each probability can be expressed as a product of determinants of two matrices of size ∼ L/2, and computed numerically in a time of order ≈ L 3 . Using translational symmetry along the x-axis as well as the conservation of winding numbers, we are able to compute any Rényi entropy up to machine precision in a time of order L 3/2 × 4 L . The results for a 22 × 22 torus are summarized in figure 3, where the two cylinders are of length y and L y − y . The figure makes it obvious that the Rényi entropy exhibits completely different behavior for different values of n. In the 'replica' phase n 1, the curves are relatively flat. In the thermodynamic limit we expect them to go to universal functions, which in principle may be computed within boundary CFT, although we will not study this phase here. For n large enough, we observe an even-odd effect like that observed for the SU (2) RVB state in [7]. We explain the origin of this transition in dimer RVB state in the next subsection.

The entropy in the locked phase n > n c
Here we derive the relation of the two-cylinder entanglement entropy of the dimer RVB state to CFT partition functions in (16) for n > n c . This follows fairly simply from the description of the long-distance behavior of classical dimers in terms of the free-boson action (14).
When one utilizes an effective field-theory action to describe the scaling limit of a lattice model, one must include all terms in the action consistent with the symmetries of the field theory. As is familiar from the derivation of the Kosterlitz-Thouless transition in field theory, operators of the form with d integer are consistent with the compactification of the boson h → h + 2π . If such a term is relevant, then the action flows to a configuration where the height h is locked to one of the minimal values h = h 0 mod 2π/d. For dimers on the square lattice, it has long been known that any such allowed term is irrelevant in the bulk; this is why the action (14) applies and the correlators are algebraically decaying. However, we showed in the preceding section 3.1 that the probabilities used to compute the Rényi entropy depend on the boundary conditions. Thus the crucial observation of [18] (see also [15]) is that one needs to include terms of the form (22) localized at the boundary.
To be precise, in the continuum limit the probabilities p σ,µ depend now on the configurations of the field p(φ). The probability of observing field configuration φ at the boundary is still Gaussian and raising this probability to the n-th power yields Therefore, the system near the boundary effectively feels a stiffness κ = nκ. Standard Coulomb-gas/field-theory computations (see e.g. [36,44]) then give the boundary scaling dimension of (22) to be d 2 /(2nκ). As a consequence, V d is irrelevant as long as d 2 > 2nκ. Thinking of n as a continuous parameter, this defines two regions, separated by a (boundary) Kosterlitz-Thouless phase transition. The critical value of n in the computation of S n is therefore where d min is the smallest d allowed by the lattice symmetries. For dimers on the square lattice d min = 1 and κ = 1/2, so that the critical value is The two qualitatively different types of behavior seen in figure 3 indeed occur on opposite sides of n c . It is important to note that even though S n should take a universal form in both phases, the value itself of the critical Rényi parameter depends on a degeneracy d min , which is non-universal. For example, d min = 2 for the quantum six-vertex wavefunction, and d min = 3 for dimers on the hexagonal lattice. This observation allows us to derive the expression for S n valid for n > n c . In this phase, the boundary operator is relevant. This means that here the boundary value of the field locks on to a minimum of V 1 . This minimum must be fixed along the entire boundary; this is commonly known as Dirichlet boundary conditions on the boson. We refer therefore to the phase with n > n c as the 'locked' phase.
A constant value of the field along the boundary means that the action is at its minimum there, so this is the maximum value of the probability p max . Universal contributions to the Rényi entropy coming from the sum (20) are therefore dominated by the configuration: This result implies that in the scaling limit there is an 'entanglement gap' between p max and the other values. We have confirmed directly the presence of the entanglement gap for the dimer RVB state by numerical evaluation of the probabilities using the free-fermion techniques detailed in appendix A. This is distinct from the quantum Hall case [48][49][50][51][52], where there is no entanglement gap, in correspondence with the gapless edge excitations there.
In terms of dimers, the configuration with maximum probability corresponds to no dimers crossing the boundaries between A and B (see figure 2). This is a huge technical simplification, since the probability of having no dimers across the cut factorizes into pieces coming from region A and from region B. Namely, define Z cyl (L x , L y ) to be the partition function (the number of dimer configurations) for a cylindrical region of size L x in the periodic direction and L y in the other, with boundary conditions corresponding to no dimers leaving the region. Likewise, Z torus (L x , L y ) is the number of dimer configurations with periodic boundary conditions. Then it follows from (21) that when an L x × L y torus is split into cylinders of lengths y and L y − y , The analogous formula still holds if space is a cylinder instead of a torus, simply by replacing the denominator with the appropriate cylinder partition function. Note that the non-universal bulk parts of the partition functions cancel in this expression. The non-universal boundary parts do not, and contribute to the boundary law in (27). Thus putting (28) together with (27) and filtering out the boundary law term gives (16). For the dimer RVB state, the partition functions on the cylinder can be computed exactly both on the lattice and in the continuum. The lattice dimer result is detailed in appendix C. For the free-boson field theory, these partition functions are known explicitly [53]; we will discuss their application to the two-cylinder Rényi entropy in section 4. In particular, we will show how the even-odd effect is a signature of this locked phase n > n c , arising from the different boundary conditions necessary for the even and odd sectors.
It should also be possible to extend these results to the region n < n c , using boundary CFT methods [15][16][17][18]. Here, however, the boundary conditions will not allow the factorization of the probability into pieces coming from the two regions. Oshikawa [16] has already treated the y = y /L y = 1/2 case for the closely related cylinder geometry, using the replica approach n ∈ N. There will possibly be subtleties involving analytic continuation in n, similar to the twointerval 1D calculation [54,55]. The study of the 'marginal' case n = n c would presumably be even more challenging [18,56].

Entanglement in the dimer RVB state
This section is devoted to the exploration of the consequences of equation (16) for the twocylinder Rényi entropy of the dimer RVB state. We will calculate explicitly the partition functions, and so find results exact in the scaling limit.
The main subtlety in this computation is how to account in the field theory for the distinct results occurring when the cylinders are of even and of odd length, as apparent in the numerical results in figure 3 for the dimer RVB state and those in [7] for the SU (2) RVB state. The reason this is possible in the field theory is that a fixed/Dirichlet boundary condition is in fact a family of boundary conditions, corresponding to the specific value h 0 that the field takes at a boundary. This value can be always be shifted overall, since the bulk action used to compute the partition functions is independent of it. However, on a cylinder, there are two separate boundaries, and what cannot be shifted away is the difference of the values on the two ends [57,58]. The general arguments given above in section (3.2) do not specify this difference; it follows rather from the particular microscopic model. For dimers, we showed in section 2.4 that where a = 0 for even L y and a = 1/2 for odd L y . The computation of the free-boson partition function is completely standard; see e.g. [53,59,60]. Since the action is quadratic in the bosonic field, it can be decoupled into an oscillator part and a classical part. Only the classical part is affected by the compactification h ∼ h + 2π , and by the boundary conditions. The partition function for Dirichlet boundary conditions on both ends of the cylinder, with the height differing by 2πa mod 2π, is given by It is convenient to rewrite such sums in terms of the Jacobi-theta and Dedekind-eta functions defined in appendix B, to take advantage of many elegant identities. The even case a = 0 can be written, using (B.2), (B.7), (B.11), (B.12) and (B.15), as where we adopt the conventional variable τ = iL y /L x , and have set κ = 1/2, the dimer value.
For the odd case a = 1/2, we have, using this time (B.2), (B.5), (B.11), (B.14) and (B.15): Plugging these results into equation (16) and using the torus partition function [44,53] gives for the two-cylinder entanglement entropies Figure 4 shows a numerical test of the universal shape equation (33) for the dimer RVB states with two different aspect ratios L y /L x = 1 and L y /L x = 2. Obviously, these n > n c results are  Numerical results for the second Rényi entropy S 2 as a function of the subsystem-ratio y = y /L y for L x = 20. Black dots are the data for squarelattice dimer RVB states, while the CFT predictions are given by equation (33) in the even case and equation (34) in the odd case. Left: torus with aspect ratio L y /L x = 1. Right: torus with aspect ratio L y /L x = 2.
quite different in the even and odd sectors. The difference is proportional to n/(n − 1), and therefore slightly diminishes with increasing n, but does not go away. This even/odd effect disappears as L y /L x is increased; it follows from the CFT expressions that this occurs exponentially quickly. It does not approach the celebrated 1D result [61] ∝ ln sin(π y) in the effective 1D limit L y /L x → ∞, because the RK Hamiltonian becomes gapped with correlation length ξ ∝ L x in this limit.
The relatively large finite-size effects apparent in the plots should disappear in the thermodynamic limit. To illustrate the slowness of the convergence, we study the exactly solvable case n → ∞, where the agreement between CFT and the lattice can be made rigorous. It follows from (20) that only the largest probability (28) contributes: S ∞ = −ln p max , and so numerical analysis is much easier. We checked system sizes up to L x = 704, and found that while indeed there are strong finite-size effects, the data converges in the end to the CFT result. The slow convergence is illustrated by the numerical and analytical results performed in figure 5, for the aspect ratio L y /L x = 1 and the two system sizes L x = 16 and L x = 96. We also computed s ∞ in the W = (0, 0) winding sector directly for the dimer model using the techniques described in the appendices, finding Even though the winding numbers are zero, the result still depends on the compactification radius because of winding fluctuations at the boundaries. Interestingly, the even-odd effect is enhanced in the W = (0, 0) winding sector, in agreement with CFT. We also observe that the odd curves in both sectors almost coincide on the lattice, as well as in the scaling limit. This coincidence can be explained from the CFT result at y = y /L y = 1/2: This is a consequence of cancelations at the free fermion point α = 2; in general

Correlations in the SU(N) RVB wavefunction
We now start to explore the question of universality, with the motivation of understanding whether the results obtained for the dimer RVB state can be adapted to describe the SU (N ) RVB states. In this section, we use quantum Monte Carlo techniques to find the spin-spin and dimer-dimer correlators for N = 2, 3, 4, 5, generalizing some of the results of [5,6]. This shows that there is no evidence for a transition as N is varied, a strong piece of evidence for universality. It also allows us to make contact with the results of [21], where an expansion around the dimer RVB state is developed.

Monte Carlo methods
To study SU (N )-RVB wavefunctions, we use an unbiased Monte Carlo technique that computes expectation values in an equal-amplitude superposition of states using a Metropolis importance sampling scheme [6]. Numerically, the RVB wavefunction is represented in a combined VB-spin (or dimer-spin) basis, |V C = |D C |Z C , where |D C is a list of site-pairs, [i, j], specifying the ends of the VBs, and |Z C is the spin state for each site on the lattice. In order to obtain expectation values (as described below), two RVB wavefunctions are sampled simultaneously, V C | and |V C . The two VB states D C | and |D C can initially be chosen at random; spin states compatible with a non-zero overlap V C |V C are then constructed. Since the spin basis is orthogonal, we set |Z C = |Z C to ensure a non-zero overlap. Then, the Monte Carlo sampling scheme proceeds in two distinct steps.
After initializing the wavefunction in one configuration |D C |Z C , the first step is to sample a new spin state |Z C = |Z C . This is done at random, subject to the condition that the overlap V C |V C remains non-zero. The second step is to sample the VB states |D C and |D C . This must be done subject to the constraint imposed by the spin states. As in the SU (2) case [6,7], one has the choice of 'local' or 'non-local' updating schemes to sample all possible VB states [62]. In a local update, the simulation chooses two parallel VBs on a square plaquette and moves each to the empty bond, subject to the condition that the spin-basis is compatible with this new configuration: figure 6(a). However, if the spin-basis is not compatible with this proposed change, as in figures 6(b) and (c) where the right pointing arrow represents a different flavor spin for an SU (N > 2) system (e.g. m i = 0 for SU (3)), the move is rejected, and a new 'plaquette flip' is proposed by the simulation. While simple to perform, this move is not ergodic in that it cannot change the winding number of the system. Therefore, this update is most useful if one wants to fix the winding number to do measurements (e.g. for the correlation measurements below, we use the winding number W = (0, 0)). In contrast, a non-local loop update moves through the ensemble of states by creating a defect at some lattice site and propagating it through the system (subject to the spin constraints) until the defect reaches its initial point, thus closing the loop [62].
The important point to note is that the Monte Carlo update for the RVB wavefunction proceeds in two steps: sampling the spin state |Z C , and sampling the VB states |D C . Each basis configuration is constrained by the other such that the VB remains a quantum SU (N ) singlet. In contrast, for the classical dimer model, no spin state exists, and dimer wavefunctions are orthogonal. We now discuss the methods we have used to generalize the SU (2) updates used in [6,7] to account for the complexities introduced by the SU (N ) RVB wavefunctions. We do this by considering the similarities and the differences between the SU (2) and SU (3) wavefunctions. Given equation (9), the singlets on a pair of sites [i, j] can be written as While sampling the spins, there are N different spin flavors to choose for a site on sublattice A. Once this spin is chosen, all other spins in the same loop are determined by equation (9). For instance, consider the SU (2) wavefunction in figure 7(a). There are two types of loops: one with ↑ (loop 1) and another with ↓ (loop 2) on sublattice A. Once the spin basis is chosen in the transition graphs, the VB configurations in the bra and ket are independently updated as prescribed above according to the same spin basis. Updates, however, are only possible between loops of the same type (1 or 2). Therefore, instead of keeping track of the spins on each site, it suffices to keep track of the loop labels at each site. Then, generalizing this a little further, consider figure 7(b). One can easily see that there are three loop types: ↑, ↓ and 0 on sublattice A. So, for an SU (N ) system, there are N different loop types, and we quickly run into ergodicity issues as updates are more readily rejected. We find that as N > 4, the number of updates that one can successfully make becomes much smaller, making the statistical error bars much slower to converge in Monte Carlo. As shown by our data, some estimators for this large N show signs of ergodicity loss; this might be cured by a more sophisticated sampling scheme.
In our VB-basis Monte Carlo simulations, we begin by computing the spin-spin and the dimer-dimer correlations. While the two-point spin correlation is well-known [63,64], the fourpoint spin interaction was recently derived in [65]. The expectation values were then generalized in [66] for SU (N ) spins. Using the notation (i j) L for two sites i, j in the same loop in the transition graph and (i) L ( j) L for sites not in the same loop [6], Here, i j = 1 if (i, j) are on the same sublattice and i j = −1 if they are on different sublattices. We also measure parallel nearest-neighbor dimer-dimer correlations, where four sites i, k, j, l are related by k = i + e α , l = j + e α and j = i + me β , where α, β = x, y, α = β, and m ∈ {0, 1, . . .}. This simplifies the estimator, which is as follows and 0 otherwise. In the next section we present data for these correlation functions for SU (N ) with N = 2, 3, 4, 5, and fit the putative universal exponents associated with their decay.

Equal-time correlators
The connected dimer-dimer correlation function has been recently shown to decay algebraically in the SU (2) case in [5,6]: at distances large compared to the lattice spacing, but small compared to the torus size L = L x = L y . Using several different methods, the exponent was found to be approximately α 1.2, independent of the direction. We extend some of these results to the SU (N ) case, mainly focusing on the correlation between two parallel dimers in the transverse direction (r 1 − r 2 = x e x ). Our quantum Monte Carlo results are shown in figure 8 for an L × L torus, with L = 64. The data clearly indicate that the algebraic decay remains for N = 2, 3, 4, 5, with α increasing with N . Since SU (N ) dimer correlations reduce to pure dimer correlations in the limit N → ∞, we know that lim N →∞ α N = α D = 2, and our data is consistent with the smooth interpolation between N = 2 and ∞. For each N we extract the leading exponent by a fit (solid blue lines in figure 8) to equation (43) in the regime 1 x L x and for odd x . We find α 2 = 1.21 (7), α 3 = 1.43 (7), α 4 = 1.54 (8) and α 5 = 1.67 (9). Notice that as N gets larger, it becomes more and more difficult to converge the Monte Carlo data in a finite time, resulting in a loss of precision on the corresponding exponent.
We also observe the presence of a dipolar term which behaves as (−1) x −2 x for all N . This direction-dependent contribution is well known in the dimer model [24], and also appears for the SU (2) RVB state [6]. In the dimer model and for the transverse correlations considered here, it cancels exactly the universal term (43) at even distances, so that the dimer-dimer correlation behaves as −4 x [10]. In the SU (N ) RVB case the exponent α is different and this does not happen. The dipolar is then a subleading term, which introduces finite-size effect when trying to extract the exponents. For example we get α 2 1.14 at even distances, which is slightly smaller than the previous computed value 1.20 obtained for larger system sizes [5,6].  These values are in good agreement with the expansion around the dimer RVB state developed in [21]. This is a cluster expansion of the loop model, describing the non-trivial inner product in (11) perturbatively in 1/N in terms of a classical interacting dimer model. To first non-trivial order, this model is exactly the one studied in [23,24], and so the Monte Carlo results in figure 26 of [24] can be utilized; in their notation X 2 = α/2 and W = 1 + 1/N . This gives α 3 1.4, α 4 1.52 and α 5 1.6, in addition to the already reported α 2 1.22. Our numerical results therefore confirm that the leading-order term in this expansion gives accurate results. Note that higher order terms could slightly increase these values, as is also shown in the supplementary material of [21].
A presumably more-accurate value for α can be found by relaxing the constraint L, and fitting the data to the full curve predicted by conformal field theory, not just the straightline segment. While the two-point function on a torus is not fixed uniquely by conformal invariance, the appropriate form is known for the free-boson field theory for any κ [53]. Thus if we assume that the free-boson/Coulomb-gas results still apply with a different κ = 1/α, the SU (2) dimer-dimer correlators are modified to with x 1. f is a universal function of the two dimensionless ratios x /L x and τ = iL y /L x and can be expressed in terms of a Jacobi-theta function defined in appendix B. We have Notice that while this function reduces to f ( x /L x , τ ) = sin(π x /L x ) in the limit of a thin torus (L y /L x 1), the approximation remains excellent for L y /L x = 1. Performing the fits to equation (44) in the same region 1 x L x as before (red solid curve in figure 8), we observe  Figure 9. Spin-spin correlation function on a 64 × 64 torus. We look at spins separated by a distance x along the x-direction: that the data reproduces very well the free-boson CFT result, even the upturn when x is of order L x . The fact that the CFT result seems to still hold here provides additional evidence in favor of universality. However, the exponents determined this way are slightly different; for example we find α 2 1.26 in the odd case, and α 2 1.18 in the even case. This small discrepancy with the previous SU (2) results [5,6], also suggested by the large N calculations [21], could possibly be resolved by studying larger systems. We also studied spin correlations, which as discussed in section 2 are known to decay exponentially in the SU (2) case: S(0)S( x ) ∼ exp (− x /ξ ) [5,28]. We find the same behavior in the SU (N ) case, with our results displayed in figure 9. We also observe that ξ decreases with N , compatible with the intuition that the RVB state becomes more and more spin-disordered as N increases, approaching zero correlation length in the dimer limit N → ∞. Fitting the data gives ξ(N = 2) = 1.30(5) (compatible with [5]) and ξ(N = 3) = 0.68 (5), ξ(N = 4) = 0.56 (3) and ξ(N = 5) = 0.46 (2). This lack of long-range Néel order is consistent with quantum spinliquid behavior.

Entanglement in the SU(N) RVB wavefunction
Throughout this paper we have emphasized the similarities between the SU (N ) and dimer RVB states. These make it plausible that the two-cylinder Rényi entropy for the SU (2) RVB state in the scaling limit can be obtained from the free-boson result (16). If this is true, this is a much stronger statement than the results for the critical exponent α as a function of N in the previous section. The reason is that even though correlators in the classical dimer model can be obtained by taking the N → ∞ limit of the loop model defined by the SU (N ) RVB inner product, the Hilbert space of the two is completely different; the latter being much larger. Thus a dimer crossing the boundary would carry much more entropy. Moreover, the reduced density matrix for SU (N ) has blocks with non-zero total spin, while these are by construction forbidden in the quantum dimer model. Thus it is very unlikely that non-universal quantities will be the same in the two cases. Figure 10. The Swap A operator. In (a), we show the original state, |V C , which is composed of a physical (or 'real') system of eight lattice sites, and a noninteracting replica (or 'ancillary'). In (b), we show the result of applying the Swap A operator, Swap A |V C = |V C . When acting the operator on this state, it is important to note that both the spins and bond bases in region A get swapped between the real and ancillary copies.
All our results are consistent with the subleading part of the two-cylinder Rényi entropy being universal. Thus it is possible that in the scaling limit, this piece still could be given by the conformal field theory result, with the identification of κ coming from the dimer correlation functions results. In this section, we provide evidence for such universality by finding the curves for arbitrary α, and extending previous numerical results [7] for SU (2) to the SU (N ) case.

Monte Carlo evaluation of the Rényi entropies
In order to calculate the Rényi entropies, the Monte Carlo sampling algorithm of section 5.1 must be modified to incorporate a replica trick procedure, first described in [35]. Namely, in order to calculate the n-th order Rényi entropy equation (1), one measures the expectation value of a Swap A (or permutation) operator acting between n copies of the system. For the SU (N ) wavefunctions studied here, we restrict ourselves to n = 2.
A naive sampling of the Swap A operator, shown graphically in figure 10, is straightforward. One simply identifies the spatial region A and its complement B, then applies the operator to a sampled RVB state, Swap A |V C = |V C (where each state |V C contains n = 2 replicas). The second Rényi entropy can be calculated from [35] where n l is the number of loops in the transition graph of the two un-swapped RVB states (the denominator), n swap is the number of loops in the transition graph using one swapped state and N denotes the number of spin flavors from SU (N ). Thus, the reconfiguration of the VB endpoints is the important ingredient of the Swap A operator for this measurement. Note, since the Swap A operator serves to swap all basis states in the spatial region A, spin states |Z C can effectively be ignored when performing this measurement; they are only used as a condition to sample the original un-swapped VB configurations. However, one discovers quickly as the lattice size (and particularly the size of region A) grows, sampling statistics become exponentially poor when using the naive Swap A operator to calculate S 2 . One must therefore employ some variant of the ratio trick described in [35]. As discussed in detail there, this involves sampling where X is a spatial region smaller than A, i.e. X A. In order to implement this measurement in the RVB wavefunction, a completely different simulation must be run with the overlap V L |Swap X |V R occurring in the weight. The essential point of the ratio trick then is that simulations must be run with the basis states corresponding to region X already swapped before measurement of the Swap A operator. Since we work in a combined VB-spin basis, one imagines the simulation beginning with two states, |V L = |D L |Z L and |V R = |D R |Z R . However, the simulation is performed with V L |Swap X |V R in the weight, meaning that |D R |Z R = Swap X |D R |Z R . Therefore, it is important to sample spin states such that the overlap Z L |Z R is non-zero.
Note that, when performing the Monte Carlo update on the VB state |D R , one does not have to consider the overlap of V L |V R since the overcompleteness condition ensures that the inner product of new VB configurations (with the same spin configuration) are always nonzero. Hence, it may be convenient to consider the 'un-swapped' state |V R for the purposes of sampling the VB configuration, where one can straightforwardly perform local (plaquette-flip) or non-local (loop) updates on the VB state as described in section 5.1.
We employ ratio-trick evaluations of the Swap A operator to calculate the scaling of the second Rényi entropy in the SU (N ) RVB wavefunction. In the next subsection, we compare our results from these Monte Carlo simulations to those obtained from conformal field theory.

Comparison of numerical results and conformal field theory (CFT)
Here we compare our quantum Monte Carlo results with curves found by using equation (16) with partition functions coming from conformal field theory for arbitrary stiffness κ = 1/α. This two-cylinder entropy is that of the quantum Lifshitz state, the continuum limit of the quantum six-vertex state. If the universality discussed above holds true, they will also apply to the continuum limit of the square-lattice SU (N ) RVB states, with appropriate identification of α in terms of N .
Let us begin with the SU (2) case, and present our results for S 2 both for unrestricted windings and for the W = (0, 0) winding sector. The quantum Monte Carlo data is shown in figure 11. We observe a strong similarity to the dimer data. In particular, there is a strong even-odd effect that becomes even bigger in the W = (0, 0) sector. This suggests S 2 for the SU (2) RVB state lies in a locked phase, just as its dimer counterpart does. This scenario is also supported by exact diagonalizations on small systems, which show that the biggest eigenvalue of the reduced density matrix is non-degenerate: SU (2) ≡ −ln( p 1 / p max ) ≈ 2.67 for a 4 × 4 torus cut in two halves, and therefore d min = 1. (The exact value for dimers is dimers = 2 ln π 2.29 in the scaling limit [67].) Such a large entanglement gap has also been observed recently for the SU (2) RVB state, in the closely related infinite-cylinder geometry [68]. We also refer to [68] for a detailed study of the entanglement spectrum in this geometry. With no restrictions on the winding numbers, the CFT result for the two-cylinder Rényi entropy at general α is The prediction of equation (35) for the W = (0, 0) winding sector is generalized to: We have made the assumption that the boundary conditions are the same as those for dimers (i.e. a = 0, 1/2 for even and odd respectively for all α). The curves (48) and (49) with α = 1.2 are plotted in figure 11 along with the SU (2) Monte Carlo data. The trend in the data is clearly toward these curves, but as with the dimer case, the finite-size effects are very strong. Several additional interesting features become apparent by plotting the CFT curves (48) and (49) at α = 2 and 1.2, corresponding to the dimer and SU (2) RVB exponents respectively. As is apparent from figure 12(left), the curve for the even branch is essentially independent of the α over a large range, and resembles the 1D result S 2 ∝ ln sin(π y). The shape of the odd curve only weakly depends on α, but the curve does appreciably shift in this range. We also observe- figure 12(right)-that contrary to the dimer case (α = 2), the two odd curves corresponding to α = 1.2 are not nearly identical in the W = (0, 0) and W = all sectors. This feature is also in qualitative agreement with the data shown in figure 11. Also plotted in figure 12 is the curve a + b ln(sin π y), with a and b fitted to make the curves as close as possible. This curve is the 1D Rényi entropy for a conformal field theory with space a circle of circumference L split into segments of length y L and (1 − y)L. The logarithmic term is the leading term in one dimension, with its coefficient b proportional to the central charge [61]. One might expect that our 2D curve would resemble the 1D curve when the 2D system is obtained by coupling together critical 1D systems while maintaining the criticality. This for example is the case for fermions with π flux discussed in [7]. However, the squarelattice quantum dimer model is gapless only in the 2D limit (just like the Goldstone modes for the Heisenberg model also discussed in [7]), so we see no particular reason for the curves to be so similar. As can be seen from the figure, the two curves are quite close, especially if one considers the numerical accuracy with which they can usually be determined. They are, however, clearly not identical.
The same picture still seems to hold for the SU (N ) RVB state. We plot S 2 data for L = 16 and SU (2,3,4) is shown in figure 13. As can be seen, the error bars get bigger and bigger for larger N , but we still observe a clear even-odd effect. This strongly suggests that we are still in a locked phase, and that equation (48) applies in the thermodynamic limit.
Extracting an accurate exponent α N from the S 2 data is difficult, due to the combined finitesize effects and statistical errors. A possible try would be to look at the even-odd difference for y = 1/2. On the lattice it is defined as δ n (L x , L y ) = S n (L y /2 + 1) − S n (L y /2) for even L y . This difference should go to a positive constant in the thermodynamic limit. From equation (48) we predict Using our biggest system-size even-odd data in (52) yields α ≈ 2.27 for dimers and α ≈ 1.57 for SU (2) RVB. These numbers can be improved by extrapolation, but the strong finite-size effects described above prevent an accurate estimate of these numbers.  Figure 13. Numerical data for the second Rényi entropy S 2 in the SU (2,3,4) case. As can be seen, even on relatively small system sizes the numerical data is not totally converged.
With such strong finite-size effects, we cannot make a definitive conclusion. Our results, however, provide qualitative evidence for equation (16) applying to the SU (2) RVB state. They also support the convention of [7] that the subleading piece of the two-cylinder Rényi entropy is universal.
Unfortunately, we do not know how to derive these results as we did for the dimers. The mapping of [15] to a classical Rényi-Shannon entropy is not valid for the SU (N ) RVB states; because of the non-orthogonality of different VB configurations, the Schmidt decomposition of the RVB state is much more complicated. Moreover, a naive extension of the result of [18] would then yield a critical value n c = α/2 0.6 for SU (2). This would imply an even-odd effect even for the von Neumann entropy, in violation of strong subadditivity [19]. The formula n c = d 2 min /(2κ) = d 2 min α/2 therefore cannot hold in general for models that cannot be mapped onto a classical Rényi-Shannon entropy. The fact that it does apply to quantum six-vertex states does not contradict strong subadditivity, because in this case d min = 2 and α 1, so n c 1 in this case 7 .
A possible way to test our CFT results while sidestepping this difficulty would be to look at another quantity, defined as follows. Consider another SU (2) RVB wavefunction | living on the torus, but where all VB states |V C with singlets crossing the boundaries between A and B have been removed. Such a wavefunction can be rewritten as the tensor product of two RVB-wavefunctions living respectively on cylinders A and B: Consider then the logarithmic (bipartite) fidelity [69] defined by taking the scalar product with the original (torus) RVB wavefunction: This fidelity enjoys several properties very similar to the entanglement entropy, like a generic area law, and a universal logarithmic divergence for 1D quantum critical points. Using equation (8), it can be expressed using classical partition functions for the RVB loop gas. In this case loops longer than the dimers can cross the boundary between A and B, but since the dimers dominate, we expect our main result (16) to apply to this quantity, upon removing the n/(n − 1) factor (or equivalently setting n = ∞). Using the arguments of section 3.2, one can also check that the equality holds exactly on the lattice for the dimer RVB state. While this is not true anymore for SU (N ) RVB states, we expect (55) to remain true at the level of universality. Therefore, equations (48) and (49) should describe the shape of both S ∞ and F, when varying the respective sizes of the cylinders A and B. Studying this fidelity, for example using quantum Monte Carlo methods similar to the ones presented here, could probably shed some more light on the phase transition scenario.

Conclusion
In this paper, we studied a family of RVB states on the square lattice, focusing in particular on the Rényi entropy S n for dividing a torus into two cylinders. This provides an elegant probe into the physics of such gapless states. We found several interesting characteristics of their quantum entanglement, by deriving exact results for the dimer and quantum Lifshitz RVB states, and by numerical work on the SU (2) and the SU (N ) RVB states on the square lattice. One of our main results is confirmation that the Rényi entropy for RVB states exhibits a two-branch structure, when the length of the entangled region varies between an even and odd number of lattice sites. We computed this exactly for the dimer RVB functions by mapping the Rényi entanglement entropy to a classically-computable Rényi-Shannon entropy, and demonstrated that the two-branch structure exists for n > 1 (but not for n 1). By measuring the Monte Carlo expectation value of a 'swap' operator on a replicated lattice, we find S 2 for the SU (N ) RVB state for N = 2, 3, 4, and find the same two-branch structure. It would be interesting to extend the Monte Carlo calculations to higher n, where the replica trick can still be used, and to study the fidelity of [69] in this context.
The two-branch shape of the Rényi entropies has an interesting consequence: strong subadditivity is violated for the 'odd' branch, which is a concave-up function of the linear size of the entangled region. We demonstrate that this is far from a numerical artifact even away from the dimer case, by calculating the shape of these branches in an effective conformal field theory. This effective theory involves partition functions that can be evaluated exactly in the continuum, and exhibit the same two-branch structure observed in lattice calculations. This emphasizes the universality of the result and confirms that the two-branch Rényi entropies are not the result of finite-size lattice effects, contrary to what generically happens in 1D systems [70,71].
From this analysis, we have elucidated two important points regarding the shape of the Rényi entropy as a function of linear size of the entangled region. Firstly, the two-branch structure does not occur for every Rényi index n, but only for values strictly larger than some (non-universal) n c , which must be greater than unity since the von Neumann entropy S 1 is strictly constrained by strong-subadditivity (while the n > 1 Rényi entropies are not). We believe that in the square-lattice RVB, 1 n c < 2; since our Monte Carlo methods require n 2, they will always find a two-branch structure as observed here. It would be useful to change d min by studying the RVB wavefunction on other bipartite lattices. For example on the honeycomb lattice d min = 3 (and therefore n c = 3 2 = 9) for dimers, and S 2 for RVB would likely lie in the replica phase, and would not give rise to an even-odd effect. This prediction should be tested by future Monte Carlo studies.
Secondly, we have found an exact expression for the shape of the Rényi entropy curves in a geometry where a linear size of the entangled regions grows as . The fact that it is consistent with our numerical data for the SU (2) and dimer RVB wavefunctions is a strong indication that it is universal. We see no particular reason why this universality should be a characteristic only of RVB wavefunctions, so this conjecture should be tested on future gapless wavefunctions. Indeed, DMRG results for the 2D transverse-field Ising model critical point have already exhibited the same type of behavior [72]. As observed in previous studies of gapless systems [7], these curves are well approximated by the celebrated universal 1D result ln(sin(π /L)); however, deviations were observed. Our exact 2D results from effective field theory indeed confirm that the 1D result is only approximate in the RVB case, but also that it is a very good approximation.
We have also explored the crossover between the SU (2) and dimer RVB states by studying the SU (N ) states. For equal-time correlators, it is possible to interpolate between the SU (2) and dimer states by taking N → ∞, but there is no guarantee that there is no transition at some value of N . Previous studies demonstrated that the SU (2) RVB state exhibited algebraic decay with an exponent α 2 ≈ 1.2, while the dimer RVB has α ∞ = 2. We showed in this work that this exponent evolves smoothly as N is varied, with no evidence for a transition. This is strong evidence of universality, allowing us to link various physical aspects of the SU (2) case with exact results on the dimer RVB. Perhaps even more importantly, we also showed that this correspondence at least qualitatively holds for the putatively universal part of the two-cylinder Rényi entropy as well. This to us is evidence both for the universality and the effectiveness of using dimer models to probe at least some highly non-trivial aspects of the SU (2) RVB state.
It is worth recalling that the SU (N ) and dimer RVB states are the simplest wavefunctions exhibiting features of gapless spin-liquid behavior. Since gapless spin liquids are long-range entangled states not amenable to classification by simple quantities such as the topological entanglement entropy, we hope that some of the universal scaling properties uncovered here will be of use in future efforts to classify and characterize them in theory, numerics and experiment.

Appendix A. Rényi-Shannon entropy for dimers on the torus
We wish to compute the following classical entropy for any real n. The probabilities p σ,µ we are interested in are given by a ratio of partition functions, which can handily be expressed using a transfer matrix T is the transfer matrix of the dimer model, and acts on the vector space generated by dimer occupancies on vertical edges along a horizontal line: a|T |b = 1 if configuration |a and |b are compatible, 0 otherwise. The denominator is calculated in appendix C. Each factor on the numerator of equation (A.2) can be evaluated using free fermion techniques [15,24,73]. The mapping goes as follows. Let us first consider a particular configuration with staggered dimers, which we call the reference configuration (see figure A.1 on the left). Superimposing the reference on any other dimer configuration generates a collection of non-intersecting lattice paths, represented by black lines in figure A.1 (right).
There is a one to one correspondence between the dimer configurations and the lattice paths, provided the latter obey a certain set of rules. To clarify them, it is most convenient to reformulate everything in terms of fermions, because their statistics will naturally encode the non-intersection constraint. Define a fermion as a vertical link occupied by a reference (red) dimer only, or a real (blue) dimer only. These are represented by black zigzag lines in figure A.1. |σ and |µ can be rewritten using fermion creation operators as where the x i and y j label the ordered positions of the N fermions. Introducing of shift of +1 lattice spacing from one row to the other, it is easy to check from figure A.1 that the transfer matrix satisfies Equations (A.5)-(A.7) conserve the number of fermions, and encode the dimer close-packed and hard-core constraints on the square lattice. To account for periodic boundary conditions we also have to impose c † L x +1 = (−1) N +1 c † 1 . Rewriting this in matrix form using Einstein's summation convention we successively get and σ |T |µ = 0|c x 1 c x 2 · · · c x N T c † y 1 c † y 2 · · · c † y N |0 (A.10) = 0|c x 1 c x 2 · · · c x N (M ) y 1 z 1 c † z 1 (M ) y 2 z 2 c † z 2 · · · (M ) y N z N c † z N |0 . (A.11) By applying the Wick's theorem and using 0|c x c † y |0 = δ x y , we finally get σ |T |µ = det This type of result is referred to as the Karlin-McGregor [74] or Lindström-Gessel-Viennot [75,76] lemma in the mathematical literature. In the end, equation (A.2) reduces to a product of two determinants, which may be computed numerically.

Appendix B. Dedekind-eta and Jacobi-theta functions
For a given modulus τ , we introduce the squared nome: The b's provide information about W x , while the number of fermions does that for W y . Defining as before T c † i T −1 = M i j c † j , diagonalizing T amounts to diagonalizing M, which can be done in Fourier space, carefully taking into account the boundary condition on fermions We write the number of dimer coverings in the (W x , W y ) sector as Z W x ,W y . Setting b = e −π u x /L x , the winding generating function is given by Z torus (u x , u y ) = W x ,W y Z W x ,W y e −π(W x u x +W y u y ) (C.5) = Tr e −πu y (N −L x /2)/2 T L y . (C.6) If we now denote by d † k the set of fermionic operators that diagonalize the one-particle transfer matrix the λ k are then the eigenvalues of M. Using this, the transfer matrix can be expressed as It is also possible to rewrite Z ± ν as e πu x + e −πu x ± 2T L y i cos (2m + ν + iu y )π L x , (C. 13) where T m (x) = cos(m arccos x) is the m-th Chebyshev polynomial of the first kind. This expression allows for convenient numerical evaluations on large system sizes.
We now wish to extract the universal CFT partition function Z, i.e. the constant term in the asymptotic expansion of Z : Z torus (u x , u y ) ∼ A L x L y B L x C L y Z torus (u x , u y ). (C.14) A possible way to do so would be to combine Euler-Maclaurin expansions and majoration techniques, as has already been done for the honeycomb lattice [77]. Here we obtain Z in a more heuristic (and non-rigorous) manner, noticing that universal properties have to come from low energy excitation near the Fermi momenta k F = ±π/2. Linearizing ln(λ k ) around the two k F and extending the products in equation (C.9) over all integers, allows one to get, after a long calculation, very similar to that in [77]: Z torus (u x , u y ) = θ 3 (iπu x /2|τ /2) θ 3 (iπ u y /2|τ/2) (2 m τ ) 1/2 η(τ ) 2 , (C. 15) where τ = iL y /L x andτ = −1/τ = iL x /L y .

C.2. Cylinder
The cylinder geometry is slightly simpler than its torus counterpart, because only the windings along y are allowed. To still express the partition function as a trace, we look at the transfer matrixT acting on the configuration of rows of length L y with open boundary conditions. This way, the number of fermions keeps track of the winding number W y . The winding generating function can be expressed, after some algebra, as Z cyl (u y ) = e −π(L y /2)u y L y m=1 1 + e π u y µ L x m , (C. 16) µ m = cos mπ L y + 1 + 1 + cos 2 mπ L y + 1 . (C.17) Here the Fermi momentum is k F = π/2, and the CFT partition function can be accessed using the linearization trick. We however have to distinguish between the even and odd L y .
Odd case. We use the same method, but care must be taken because of the presence of a zeromode for m = (L y + 1)/2. We have Z cyl (u y ) = 2 cosh(π u y /2) In the end we obtain: Z cyl (u y ) = θ 2 (iπ u y /2|τ /2) η(τ /2) (C. 23) and once again recover the CFT result (32) after setting u y = 0 and modular transformation.

C.3. Zero-winding sectors
With the winding generating functions at hand, the calculation of s n (y, τ ) in the W = (0, 0) winding sector becomes straightforward. For example we have in the odd case (recall τ = iL y /L x andτ = −1/τ ) where we have selected the constant term in the product of the two theta functions. This allows us to recover equations (50) and (51) after once again performing a modular transformation.