Solvable models for 2+1D quantum critical points: Loop soups of 1+1D conformal field theories

We construct a class of solvable models for 2+1D quantum critical points by attaching 1+1D conformal field theories (CFTs) to fluctuating domain walls forming a ``loop soup''. Specifically, our local Hamiltonian attaches gapless spin chains to the domain walls of a triangular lattice Ising antiferromagnet. The macroscopic degeneracy between antiferromagnetic configurations is split by the Casimir energy of each decorating CFT, which is usually negative and thus favors a short loop phase with a finite gap. However, we found a set of 1D CFT Hamiltonians for which the Casimir energy is effectively positive, making it favorable for domain walls to coalesce into a single ``snake'' which is macroscopically long and thus hosts a CFT with a vanishing gap. The snake configurations are geometrical objects also known as fully-packed self-avoiding walks or Hamiltonian walks which are described by an $\mathrm{O}(n=0)$ loop ensemble with a non-unitary 2+0D CFT description. Combining this description with the 1+1D decoration CFT, we obtain a 2+1D theory with unusual critical exponents and entanglement properties. Regarding the latter, we show that the $\log$ contributions from the decoration CFTs conspire with the spatial distribution of loops crossing the entanglement cut to generate a ``non-local area law''. Our predictions are verified by Monte Carlo simulations.

For gapped topological systems like SPTs, the decorated domain wall construction provides a way of constructing a D-dim SPT starting from a (D − 1)-dim SPT [25].The idea is to decorate domain walls of a (typically G = Z 2 ) symmetry with a (D − 1)-dim SPT, and to let the domain walls fluctuate in order to restore G.Is it possible to generalize this dimensional "bootstrapping" approach to gapless systems?In this work, we answer this question affirmatively by providing a decorated domain wall construction in which 1+1D CFTs are attached on domain walls, leading to a gapless 2+1D theory with remarkable properties.
The idea of attaching gapless 1D theories on fluctuating domain walls is motivated by a number of physical systems.First, the combination of frustrated magnetism and itinerant degrees of freedom appears in so-called "charge-ice" systems [26,27].Second, decorated domain wall models were shown to be relevant to certain models combining spin and orbital degrees of freedom [28], including the extended Hubbard model on the Kagome [29] and checkerboard [30] lattice.Third, when a two-dimensional antiferromagnetic insulator is doped, the charge concentrates into domain walls forming "metallic rivers" and effectively behaving as Luttinger liquids in an active environment [31][32][33].Fourth, if we reintepret the domains as distinct gapped topological phases, the gapless degrees of freedom appearing at the domain walls would have a natural intepretation as edge modes, like in the quantum Hall Figure 1: Example of a decorated domain wall configuration.σ z = ±1 spins are represented by filled and empty blue dots, respectively.τ spins (red dots) are located on every site of the hexagonal lattice (dashed blue lines), including the domain walls of the σ z spins (green lines).As an example, τ 1 and τ 2 are located on the domain wall between sites p and q.  and c) illustrate hexagon solid and stripe configurations, consisting of contractible and non-contractible domain walls of the shortest length, respectively.plateau transition [34].Finally, since Ising domain walls have been proposed as a model for strings [35][36][37][38], the decoration described here is analogous to the fermionic degrees of freedom which are added to form super-strings [39].
From a spin-liquid perspective, we will see that our construction starts with a familiar route, in which local frustration (in our case Ising antiferromagnetism on a triangular lattice) leads to an exponentially large number of classical ground states which is described by a familiar ensemble of fully packed loops (FPL) [40][41][42].However, the source of the exotic criticality we will describe below -fully packed self-avoiding walks described by a non-unitary c = −1 theory -lies elsewhere: It is due to an additional kind of non-local frustration whereby all allowed loops on the honeycomb lattice realize an effectively twisted boundary condition [43] for the CFT which lives on it.

Model and solution
We introduce a local model which attaches any translation-invariant 1D Hamiltonian H 1D to Ising domain walls.The model contains σ Ising spins living on a triangular lattice, and τ decoration degrees of freedom (DOFs) living on the vertices of the dual honeycomb lattice (see Fig. 1).The τ operators could be anything (spins, bosons, fermions,. . .), but we will take them to be spins for concreteness.The Hamiltonian couples the τ DOFs only along the domain walls of the σ spins.This is easily implemented on this lattice, as we now show for an example which involves nearest-neighbor coupling terms H dec (τ 1 , τ 2 ) between τ DOFs (e.g. where 〈pq〉 denotes bonds of the triangular lattice and τ pq 1 and τ pq 2 are the spins attached to the honeycomb bond crossing 〈pq〉(see Fig. 1).One can easily generalize this construction to the case of H dec with longer-range terms, as shown in Appendix A.
The Hamiltonian 1 is block diagonalized in the σ z basis.There is a two to one mapping between {σ z } configurations and their domain wall configurations, which form a set of noncrossing loops L. This enables us to rewrite the Hamiltonian as a sum of 1D Hamiltonians living on each loop l: where H 1D = i H dec (τ i , τ i+1 ) is the 1D Hamiltonian given by the sum of the H dec terms along the loop.Each H eigenstate, denoted |Ψ〉, is a tensor product of an H 1D eigenstate |ψ 1D 〉 on each loop l.We will be mostly interested in the case for which |ψ 1D 〉 is the ground state of H 1D , given by:2 The total energy of this state is the sum of the 1D ground state energies, and where L l is the length of loop l.Which configurations L minimize the total energy E[L]?Let us for now neglect finitelength effects in E GS (L) and only keep the leading-order term in L: E GS (L) ≃ ε 0 L, with ε 0 the ground state energy per site.If ε 0 < 0, minimizing the energy simply amounts to maximizing the total length of all domain walls, which corresponds to a "fully packed" loop configuration for which each honeycomb site is visited by a loop.The corresponding {σ z } are the ground states of the classical Ising AFM on a triangular lattice, first studied by Wannier.These configurations lead to a residual entropy per spin of 0.323066 [44].Note that ε 0 can always be made negative by adding a term −J1 to H dec with a large enough J.For the sake of simplicity, we will work in the limit of J → ∞ for the remainder of this work, which means the only allowed loop configurations are fully packed.
At this point, we have an exponential degeneracy between FPL configurations.However, this degeneracy will in general be split by the finite-length corrections to the ground state energy of each chain.For the rest of the letter, we will focus on the case when H 1D realizes a conformal field theory (CFT).For a CFT with periodic boundary conditions, the finite-size correction to the ground state energy -also called the Casimir energy -is expected to be universal and proportional to the central charge c: E GS /L = ε 0 − πc/3L 2 + . . . 3 Naively, this should put a damp on our hopes of realizing a gapless theory: since c > 0 for unitary theories, this means the Casimir energy is negative and the energy per site is thus minimized for short loops.In this scenario, the lowest-energy loop configurations are thus the ones which pave the plane with the shortest possible loops (which are hexagons of length 6), leading to a 3-fold breaking of lattice translation in a "hexagon solid" phase (see Fig. 1.b) [29].Since each loop has finite size, the 1D Hamiltonian living on it has a finite gap, and the theory is thus gapped.
However, if c was negative, the energy per site would be minimized by taking L → ∞, leading to a loop which is macroscopically long and which can thus host a 1+1D theory with a vanishing gap.Interestingly, this scenario can be realized by choosing CFTs which are "frustrated" for certain chain lengths [43] and which can effectively behave as if c < 0, as we now show.A key insight is that, for certain CFTs, the value of c which appears in the Casimir energy is not always equal to the actual central charge and may depend on the chain length modulo some integer.We will focus on cases where that integer is 4, leading to a more general formula Figure 2: The energy density of the ground state of the H X-ZXZ chain as a function of its length L (markers) along with a fit to Eq. 4 for the c 0 = 1 and c 2 = −2 branches (curves).Specifically, the figure shows the graph of E GS /L + ε 0 vs 1/L 2 , where the thermodynamic limit value for the energy density of ε 0 = 2/π was used.Although higher-order terms in 1/L could in principle become sizable for short loops, we find that the −πc/3L 2 formula (curves) for the energy density agrees extremely well with the exact values (markers) down to L = 6, which corresponds to the shortest loop allowed by the honeycomb lattice.
for the finite-size ground state energy density:4 For example, consider the following 1D Hamiltonian, which is the decoration model we will focus on for the rest of the article:5 It describes the quantum phase transition between a 1D trivial paramagnet and an SPT protected by a Z 2 × Z 2 symmetry which is generated by i τ x 2i and i τ x 2i+1 [45]. 6,7The mod 4 effect in the energy density of Eq. 5 can be calculated analytically (see Appendix B) and gives c 0 = 1 and c 2 = −2 (see Fig. 2). 8The most relevant aspect for us will be that c 2 < 0. Another example of a CFT with c 2 < 0 is the doubled version of the XX chain with where Q is at the corner of the Brillouin zone (see also Appendix.E.1 for more details).The linear fit gives η = 0.65 ± 0.01.The calculations have been performed at T ℓ 2 ≡ T = 4.
) which has c 0 = 4 and c 2 = −11 9 (see Appendix C for a detailed calculation and for more examples of decoration Hamiltonians).
Another crucial insight is that, on the honeycomb lattice, all contractible loops forming an FPL configuration have a length given by L = 4k + 2 with k an integer (see Appendix D for a proof).Putting aside non-contractible loops for now, this means the low energy properties of the system are only determined by the sign of c 2 .For c 2 > 0, the system forms a gapped solid of short loops, as explained before.However, for c 2 < 0, the Casimir energy is positive and thus the minimal energy per site is obtained for L → ∞.For a 2D system of linear size ℓ, the maximal loop length scales like ℓ 2 and is obtained for a single loop visiting every site of the lattice exactly once (also called a Hamiltonian walk, or fully-packed self-avoiding walk (FPSAW), or "snake" in the rest of the paper).Since the gap ∆ of H 1D scales like the inverse length L of the chain on which it lives, one finds ∆ ∼ 1/L ∼ 1/ℓ 2 , which indicates a 2D gapless theory with dynamical critical exponent z = 2.
Since there is an extensive number of degenerate "snake" configurations (with entropy per site s ≡ S/N ≈ 0.130812) [49], the low-energy physics is described by a statistical average over them, which is obtained in the zero-T limit of the following thermal density matrix: where and where ).In the zero-T limit, 10 the ensemble of {σ z } described by ρ maps to the O(n → 0) fully packed loop model, which is described by a c = −1 non-unitary CFT [49].This leads to unusual power laws for correlation functions.For example, the antiferromagnetic correlations captured by have a power law envelope C Z Z (x) ∼ |x| −η , which we extract from the structure factor (see Fig. 3 and Appendix.E.1).We find η = 0.65 ± 1, which agrees well with the prediction of η = 2/3 based on a Coulomb gas description of the c = −1 CFT (see Appendix E.2).This also shows conclusively that the snake phase is in a different universality class from the standard Ising triangular lattice antiferromagnet (also known as the O(n = 1) fully-packed loop model) which has η = 0.5 [50].
Figure 4: The average number of domain walls 〈n DW 〉 vs rescaled temperature T = T ℓ 2 , for an ℓ by ℓ torus.(Due to an even-odd effect, the data for even and odd ℓ was separated into two panels for clarity.)The plot indicates a phase transition between a snake phase at high T for which 〈n DW 〉 does not scale with ℓ, and a stripe phase at low T for which 〈n DW 〉 grows with ℓ.
We now discuss correlation functions of τ DOFs.Denoting φ(x, t) a scaling operator in the decoration CFT with scaling dimension ∆ φ , purely temporal correlations have the same value for any snake configuration, leading to: By contrast, for spatial correlations the average over snake configurations leads to an averaging over the 1D distance d measured along the snake which appears in the correlator ∼ 1/|d| 2∆ φ .Assuming φ is chosen so that it has no lattice-scale oscillations, we expect it is safe to replace |d| by its average value, which scales like |d| ∼ |x| 1/ν , with ν = 1/2 the known geometrical exponent of the FPSAWs.This leads to the following prediction for spatial correlations: To summarize, temporal scaling dimensions in the 2+1D theory are the same as for the underlying decoration 1+1D CFT, whereas spatial ones are multiplied by 2. This behavior is manifestly consistent with z = 2. (It is instructive to contrast this behavior with conformal quantum critical points (CQCPs) [51][52][53], which are 2+1D critical points constructed starting from a 2+0D CFT, and which often have z = 2.For CQCPs, the spatial scaling dimensions are the same as the underlying 2+0D CFT, whereas the dynamical ones are divided by 2.)

Finite temperature and finite size
How do the previous results survive at finite temperature?In order to study this, we have developed a worm Monte Carlo algorithm [54][55][56] (see Appendix F for more details) which probes the density matrix in Eq. 6.We calculated the average number of domain walls 〈n DW 〉 as a function of temperature (see Fig. 4), from which we can also infer the typical length of each domain wall L ∼ ℓ 2 / 〈n DW 〉.We find that, in the snake phase, the average number of domain walls scales as 〈n DW 〉 ∼ ℓ 2 /β, and the typical length of each domain wall therefore scales as L ∼ β.Finite temperature thus provides a spatial infrared cutoff of size β for the chain lengths on which the 1+1D CFTs live, along with the usual temporal cutoff β in imaginary time.This means that, as T approaches zero, both the spatial and temporal dimensions of the 1+1D tori hosting the CFTs diverge like β, in contrast to usual quantum critical scaling for which only the imaginary time dimension scales with β whereas the spatial dimension is independent of β.
While the discussion so far was done in the limit of ℓ → ∞ before T → 0, let us now consider a finite-size system, focusing on an ℓ by ℓ torus for concreteness.The main new ingredient on the torus is the presence of non-contractible loops.Indeed, the geometrical constraint that loop lengths have L = 4k + 2 for FPL configurations only applies to contractible loops.This means non-contractible loops can have L = 4k, and can thus be on the c 0 > 0 branch, which is lower in energy (see Fig. 2).Filling the space with such loops amounts to forming a stripe configuration for the σ z spins (see Fig. 1.c).The difference in total energy between the snake and the stripe configuration is calculated easily: ∆E = πc 0 6 +O(1/ℓ 2 ).Since c 0 > 0, the stripe phase is lower in energy, and is actually the true ground state for a finite-size torus.However, the snake phase has a zero-temperature finite entropy density s, whereas the stripe phase does not.Any small temperature should thus be enough to stabilize the snake phase at the expense of the stripe phase.Let us consider the difference in free energy between the two phases, given by ∆F = ∆E − 2T ℓ 2 s.Since ∆E = O(1), the entropy term dominates already at a temperature which is parametrically small in system size, motivating the definition of a rescaled temperature T ≡ T ℓ 2 .Setting ∆F = 0, we predict a first-order phase transition 11at Tc = πc 0 /12s between a stripe phase at low T and a snake phase at higher T .This gives Tc ≃ 2.0013 for the decoration Hamiltonian H X-ZXZ with c 0 = 1, which is confirmed by our numerics (see Fig. 4).For T > Tc , we observe that the average number of domain walls is independent of system size, consistently with the snake phase.For T < Tc , 〈n DW 〉 shows an increase with ℓ, consistently with a stripe phase.
In practice, since numerics is done for finite ℓ, we have to work at T > Tc in order to probe the snake phase, for which 〈n DW 〉 is not strictly equal to one, but remains of order O(1) (for example, we used T = 4 to generate data for Fig. 3, for which 〈n DW 〉 is slightly above 2).At any rate, local properties should not be able to distinguish between a single or an O(1) number of snakes.5 vs. x for different system sizes ℓ.The black curves are fitted to x ≥ 4 numerical data points according to Eq. 9.The snake configurations were generated at T = 4. Details on how S was calculated numerically, and about the non-universal contributions which were dropped, can be found in Appendix G.

Entanglement
Should we expect an ℓ A log(ℓ A ) term (with ℓ A the linear size of subsystem A) in the bipartite entanglement of a 2D soup of 1+1D CFTs such as the one we have constructed here (as suggested in Ref. [57] for example)?Surprisingly, we find that this is not the case, due to a subtle property of the spatial distribution of loops crossing the entanglement cut which restores the area law in 2D. 12 (It is however a non-local form of the area law, as we will explain).
Let us first give a physical argument for the restoration of the area law, followed by numerical results.As shown in Fig. 5, for any given snake configuration L, unfolding the snake to a straight line results in a periodic chain for the τ DOFs in which the subsystem A is mapped into several disjoint intervals.Hence, to determine the entanglement entropy S[L] of a given state |Ψ[L]〉, we only need to apply the known formula for the entanglement of disjoint intervals in a 1+1D CFT [58][59][60] (refer to Appendix G.1 for more details).In order to derive the entanglement scaling, we assume it is sufficient to approximate this formula by S ∼ (c/3) where t i is the length of interval i and N int is the number of intervals.We then average S[L] over L which, assuming a translation-invariant entanglement cut like that of Fig. 5 for simplicity, effectively corresponds to an average over interval lengths denoted by f (t) ≡ d t p(t) f (t) with a distribution p(t i ) ≡ p(t).We also know that N int = 2ℓ A /3 since each honeycomb edge crossed by the entanglement cut has a 2/3 probability of being occupied by a loop strand. 13verall, this leads to a simple prediction for the entanglement: The remaining task is thus to find how log(t) scales with ℓ A .Although p(t) is a priori unknown, we already know that t ∝ ℓ A since the total length of all intervals should scale like ℓ 2 A because the snake visits every site inside subsystem A. If we assumed that log(t) ∼ log t , we would thus find a ℓ A log(ℓ A ) term for the entanglement.This is incorrect however because, as we argue in Appendix G.3, p(t) effectively describes the distribution of the exit time at which a random walker first exits subsystem A, having started inside A one lattice spacing away from the entanglement cut.Such a process is expected to follow a Lévy-type distribution, for which the typical value t typ is O(1) even though the average value t diverges with ℓ A due to a long-time tail.One thus finds log(t) ∼ log t typ = O(1), and the area law is restored: S ∼ ℓ A .
We now focus on a strip subsystem of width x in a ℓ by ℓ torus (see Fig 5) [61].Using the random walker model mentioned above, we derive a distribution p(t) and calculate log(t) as a function of x (see Appendix G.3 for details), leading to the following prediction for the entanglement: with x = ℓ π sin πx ℓ .This formula gives a good agreement with our numerical results, as shown in Fig. 6 (see also Appendix G.4).In Eq. 9, the dependence of the area law prefactor on x due to the term proportional to B is unusual: it does not appear in the standard scaling forms used for gapless 2D systems, and leads to a much larger dependence on x of S strip than usually observed [61,62].The B term depends crucially on contributions from parametrically long intervals and thus reveals the non-local character of the area law: the distribution p(t) has a long-time tail which extends until the "Thouless time" t Th = x 2 /D, with D the effective diffusion constant of the random walker.(Another example of a non-local area law was recently proposed in Ref. [57]).

Discussion
By attaching CFTs with positive Casimir energy on domain walls, we have shown how to realize a 2+1D quantum critical point featuring a single domain wall visiting every site of the system, whose statistical fluctuations are described by the O(n = 0) fully-packed loop ensemble.
On the border of which phases does this QCP exist?There are at least two types of relevant perturbations.The first type is obtained by adding terms to the decoration Hamiltonian.In that case, any relevant perturbation of the 1+1D CFT would also be relevant for the 2+1D QCP.For example, for the case of H X-ZXZ considered in this work, we can use the following interpolation Hamiltonian as decoration: with 0 ≤ α ≤ 1.In this notation, the unperturbed decoration Hamiltonian H X-ZXZ corresponds to α = 1/2, at which the QCP occurs.For α < 1/2 (resp.α > 1/2), H 1D (α) flows to a trivial (resp.non-trivial) Z 2 × Z 2 1D gapped bosonic SPT [45].
When the decoration Hamiltonian flows to a gapped phase, we expect the spin degrees of freedom σ z to flow to the conventional triangular lattice Ising antiferromagnet, also known as the O(n = 1) fully-packed loop ensemble [40,41].Indeed, in the limit of vanishing correlation length for the decoration DOFs, the ground state energy of H 1D becomes independent of the domain wall length, and all fully-packed loop configurations become equally likely.For the example at hand, the QCP thus separates two phases of fully-packed loops with n = 1 fugacity which are decorated by a trivial (resp.non-trivial) 1D gapped Z 2 ×Z 2 SPT.(The 2D Hamiltonian obtained by using H 1D (α) as decoration corresponds to the same interpolation Hamiltonian between 2D Z 3  2 SPTs introduced in Ref. [63] (See also Refs.[64][65][66][67][68]), but with the addition of an infinitely large nearest-neighbor antiferromagnetic coupling on one of the three triangular sublattices in order to enforce the fully-packed loop constraint on the σ z spins).The class of QCPs we have proposed appear thus naturally at the transition between different 2D gapped SPTs with an underlying domain wall structure.Further, the kind of statistical average over snake configurations we constructed could describe a phase transition between average SPTs, which were introduced recently [69,70].
A second type of relevant perturbation is obtained by adding terms to the Hamiltonian for the σ degrees of freedom.The most natural term to add would be a transverse field term σ x , which would generate quantum fluctuations between the snake configurations.In fact, quantum fluctuations between decorated loops were studied in the context of the Kagome Hubbard model in Ref. [29].In that work, the relevant classical loop configurations formed a solid of short loops (the "hexagon solid" pictured in Fig. 1b), and the quantum fluctuations thus naturally stabilized a "plaquette" phase of resonant short loops.In our case, the relevant classical configurations feature a single snake, and the impact of local quantum fluctuations on such a non-local object seems difficult to predict a priori.Nonetheless, the possibility of stabilizing a quantum superposition of snake configurations, in analogy with earlier work work [51,71,72], is interesting and left for future work.
We note that the snake phase could be amenable to a fermionic description after performing a Jordan-Wigner transformation on the τ DOFs along the snake.A similar construction was used to study spin liquid models like the Kitaev honeycomb model [73,74], but for a fixed snake configuration with a simple geometry which necessarily breaks certain spatial symmetries, whereas in our case the symmetry is restored by averaging over snake configurations.Once expressed in terms of fermions, the transition from the stripe phase to the snake phase would be reminiscent of the smectic and/or nematic transitions in electronic systems [75].
Our construction generates a new kind of non-local frustration by combining two primary ingredients: (1) a non-trivial dependence of the Casimir energy of the decoration CFT on the chain length modulo some integer, and (2) a geometrical constraint which forces all loops on a given lattice to have a certain length modulo some integer.The first ingredient was recently understood as arising from effectively twisted boundary conditions for certain chain lengths [43].More generally, our work demonstrates the importance of the coefficients c r (which encode the Casimir energy dependence on the chain length modulo some integer) as an additional property of a CFT Hamiltonian with crucial physical consequences.It is remarkable that two CFTs with the same central charge can lead to completely different physics when used as decoration in our model due to their different sign for c 2 (e.g. the XX chain versus H X-ZXZ ).One natural generalization of our construction is to consider other lattices for which domain wall lengths are constrained in some other way, which could magnify the effect of the c r coefficients for r ̸ = 2. Another generalization is to consider other decoration Hamiltonians beyond the ones we have proposed here.This motivates further work on classifying the possible c r sequences for known CFTs, especially beyond c = 1 theories [43].
Finally, our calculation of the entanglement, which combines known results about entanglement in 1+1D CFTs with properties of exit time distributions, allowed us to explain the somewhat surprising presence of an area law and could be useful in other contexts, including other kinds of constrained models like that of Ref. [57].Since the bipartite entanglement turned out to follow an area law, one wonders whether other measures of entanglement could provide a more direct probe of the non-locality arising from parametrically long intervals.

A Generalization of the decoration Hamiltonian to further neighbor interactions
In the main text, we showed how to construct a 2D local Hamiltonian which attaches to domain walls any 1D decoration Hamiltonian composed of nearest-neighbor terms.In this appendix we generalize this construction to terms involving three consecutive sites, using for concreteness the example of the H X-ZXZ Hamiltonian we consider in the main text.The same procedure can be used for terms involving any number of sites.
For the X-ZXZ model given in Eq. 5, ) should be attached to domain walls of σ spins.Since H dec involves three honeycomb vertices and correspondingly two honeycomb bonds, we need two domain wall projectors to modify Eq. 1: where the first sum is over triangular plaquettes of the triangular lattice, the second sum is a sum over the three choices of pairs of edges b 1 ̸ = b 2 belonging to the triangle △ pqr , and where 1,2,3 are the three τ spins on the honeycomb vertices which are connected by the honeycomb bonds dual to b 1 and b 2 (see Fig. 7 for an example).

B Solution of the X-ZXZ model
In this section, we review the mapping of the H = X − Z X Z Hamiltonian to free fermions.We begin with a Kramers-Wannier duality transformation, for which it is convenient to split the Hilbert space into two sectors: i τ x i = ±1.The constraint for the positive sector can be satisfied with a standard duality transformation: τz i = τ z i τ z i+1 and τ x i = τx i−1 τx i for all i.Then it follows that Eq. 5 maps to For the negative sector, we employ the same mapping with the modification to satisfy the i τ x i = −1 constraint.Then Eq. 5 maps to H + and H − can then both be mapped to free fermions using a Jordan-Wigner transformation where s = (−1) L−N +1 for L sites and N particles.In this free fermion language, ground state energies can be easily computed.In the case of even L, s = −1 for both sectors, since the number of domain walls N must always be even.Hence, for even L, the positive and negative sectors possess antiperiodic and periodic boundary conditions, respectively.For a chain of length L = 0 mod 4, the ground state occurs in the positive sector with half-filling N = L/2.The corresponding energy density is given by since the momentum is quantized as k = 2π n + 1 2 /L in the case where L − N is even.Eq.B.4 can be expressed in closed form as E GS /L = − 2 L 1 sin(π/L) .For a chain of length L = 2 mod 4, we are restricted from N = L/2 since the number of particles must be even.Hence, the ground state occurs in the positive sector with N = L/2 ± 1.These two degenerate ground states have an energy density given by which can be expressed in closed form as E GS /L = −2 L cot π L .The Casimir energy then can be obtained from these closed form expressions in the limit of large L: According to Eq. 4, this consequently leads to c 0 = 1 and c 2 = −2.

C Other CFT Hamiltonians with a c r < 0 branch
Like the X-ZXZ model, there exist other 1D chains whose ground state energies, as a function of their length, have positive and negative c branches.Here we provide two other examples of such models.

C.1 Spin-1 chain
The model is a spin-1 chain with the Hamiltonian                              It can be numerically confirmed that the system exhibits the even-odd effect at any values of −1 < ∆ 1, where the system is in its gapless phase.However, the values of c 0 and c 1 are not universal and depend on ∆.
This model undergoes a phase transition at D c ≈ 1 between a trivial paramagnetic phase for large D and the Haldane phase for small D [48,77], which exhibits a 2 × 2 SPT [78][79][80].
The critical point belongs to the same universality class as the X -ZXZ chain and is described by a c = 1 CFT [79].As shown in Fig. S2 (a), the ground state energy of the system with periodic boundary conditions, has branches of c 0 > 0 for even system sizes and c 1 < 0 for odd ones.
In order to realize the snake physics with this Hamiltonian, it is necessary to have a mod 4 effect instead of an even-odd effect.This can be accomplished by using two independent copies of the spin-1 chain.More precisely, the ground state energy of at the quantum critical point D = D c has c 0 > 0 for L = 0 mod 4 and c 2 < 0 for L = 2 mod 4.

C.2 XXZ model
The spin-1/2 XXZ model the Hamiltonian The X X Z model in its critical phase exhibits an even-odd effect when it has periodic boundary conditions.The plot of energy density is provided in Fig. S2 (b-d) for ∆ = −0.5, ∆ = 0.5 and ∆ = 0 (X X chain).Similar to the previous examples, the doubled Hamiltonian, leads to a mod 4 effect.In the special case of ∆ = 0 (XX model), the model can be solved exactly using the Jordan-Wigner transformation, allowing us to analytically calculate c 0 and c 2 as explained below.First, we determine the relationship between the Casimir coefficients of the single and double chain models.The doubled model consists of two independent chains, each of length L/2, so we have: Now, for a single X X chain, the Jordan-Wigner transformation yields the ground state energy as Therefore, according to Eq. C.6, for the doubled X X model, we find c 0 = 4 and c 2 = −11.

D Proof of L = 2 mod 4 for contractible loops in a fully packed configuration
In this appendix, we prove that contractible loops in a fully packed configuration cannot be of length L = 0 mod 4, and must thus be of length L = 2 mod 4 (since loop lengths are always even).The proof is in two steps.First, we show that contractible loops with L = 0 mod 4 necessarily have an odd number of honeycomb vertices in their interior.Second, we show that an area of the honeycomb lattice with an odd number of vertices cannot host a fully packed configuration.By combining the two results, we conclude that a contractible domain wall with L = 0 mod 4 is not consistent with a fully packed configuration.

D.1 First result
In this first result, we show that contractible loops with L = 0 mod 4 necessarily have an odd number of honeycomb vertices in their interior.A simple proof can be found in Ref. [82], which we reproduce here for convenience.
Let L be the length of the loop, let k be the number of hexagons inside the loop, and let x be the number of vertices which are strictly inside the loop.Let z denote the number of obtuse loop vertices (i.e.vertices at which the interior angle drawn by the loop is 120 degrees) and y the number of reflex loop vertices (i.e.vertices at which the interior angle drawn by the loop is 240 degrees).If we orient the loop counterclockwise, the obtuse vertices correspond to left turns and the reflex vertices correspond to right turns.We know that z + y = L since each loop vertex is either obtuse of reflex.We also know that z − y = 6 since the loop needs to close: it needs to do 6 more left turns than right turns.
Consider cutting each interior hexagon into 12 right triangles by cutting along all its axes of symmetry.There are two ways of counting the number t of such triangles: t = 12k, but

D.2 Second result
A fully packed configuration must have a domain wall passing through each vertex of the honeycomb lattice.If x is the number of vertices inside an area of the honeycomb lattice, a fully packed configuration in that area is a configuration of loops such that x = i L i , where i is the loop index.Since each L i is even, there can be no fully packed configuration if x is odd.

E More details on the σ z σ z correlator E.1 Numerical results
As usual for triangular lattice antiferromagnets, the spin correlation function C Z Z (x) = 〈σ z (0)σ z (x)〉 has lattice scale oscillations corresponding to the K and K' wavevectors at the corner of the Brillouin zone.On top of these lattice scale oscillations, in the long distance limit, we expect the spin correlation C Z Z (x) = 〈σ z (0)σ z (x)〉 to have a power law decaying envelope: C Z Z (x) ∼ |x| −η .Assuming the structure factor defined by x has a sharp peak at k = Q, we find C Z Z (Q) ∼ ℓ −η .Fig. 9 shows that Q is indeed at the corners of the Brillouin zone.The value of η can be obtained using the graph of C Z Z (Q) vs. ℓ (Fig. 3), which gives η ≈ 0.65.In Appendix E.2, we propose a CFT argument based on Coulomb gas which predicts η = 2/3.

F Worm Update for Monte Carlo algorithm
Fully-packed loop configurations on a hexagonal lattice can be ergodically sampled using a classical Monte Carlo algorithm with worm updates [54][55][56].In such loop configurations, all vertices will generally touch an even number of occupied bonds (links).If B is the complete set of occupied bonds constituting a given loop configuration, then the worm update is initiated by attaching a single bond to B with vertices a and b: B ′ = B ∪ a b.In this case, a and b are referred to as vertex defects since they touch an odd number of occupied bonds, hence removing the bond configuration from the subspace of valid loop configurations.To return back to this subspace, one random bond whose boundary includes either a or b is added or removed to B ′ , thereby updating the pair of vertex defects, and this process is repeated until the two vertex defects coincide: a = b.In our case, this cluster update can then be accepted or rejected with standard Metropolis sampling according to the Boltzamann weight given by e −F [L]/T .This cluster update procedure is summarized in page 12 of [54].
Here, F [L] is the free energy of the domain wall configuration L, which is the sum of the free energies of each loop: F [L] = l∈L F 1D (L l , T ).In the low-T limit, F 1D (L, T ) can be expanded as which according to B.6, for the X − Z X Z chain will be The log(2) term appears due to the fact that the L = 2 mod 4 chain is doubly degenerate.By dropping unimportant constant terms and rewriting the above relations in terms of T , we obtain

G Calculations of the entanglement G.1 Entanglement on the 1D chain
In order to calculate the entanglement of a given eigenstate |Ψ[L]〉, we need to calculate the entanglement of the τ dofs on the snake.As show in Fig. 5, after unfolding the snake, the subsystem A is mapped to a union of disjoint intervals on the 1D chain.Following [58][59][60], Renyi entropies of disjoint intervals in a CFT can be computed using the following analytical formula where each interval (i running from 1 to N ) is between u i and v i , C n are non-universal constants, and F N ,n (x) is a model-dependent scaling function of all the invariant ratios which can be constructed out of all the 2N endpoints of the intervals.This formula is interpreted as the correlation function of twist operators inserted at each interval endpoint, where the twist operators have the universal scaling dimension In this work, we will limit ourselves to the dominant contributions to the entanglement in the limit of large perimeter ℓ A of the subsystem, which will turn out to follow a "non-local area law".For this reason, we omit the contribution from F N ,n (x) because it only gives a term of order O(1).Further, for the strip geometry we will consider below, we will show in G.3 that the contribution from the C N n factor leads to a simple, local contribution to the area law which model-dependent and not particularly interesting.We thus also omit it in the following.
The n th Renyi entropy, S n = − 1 n−1 log Tr[ρ n A ] , is then given by where the dots represent the non-universal contributions from F N ,n (x) and C N n which we dropped.In particular, the von Neumann entanglement entropy is obtained by taking n → 1, leading to Also, since we will always work with a finite-length snake with periodic boundary conditions, one should use the following replacement in Eqs.G.2 and G.3: and similarly for |v i − v j | and |u i − v j |, where L is the length of the chain.
To test the validity of Eq.G.2 for the X − Z X Z chain in the case of multiple intervals, we used the DMRG method described in Section 6.2 of Ref. [60], which computes the second Renyi entropy by performing "twist" operations between two MPS states at the endpoints of each interval.Specifically, we compute S 2 for four equal partitions: u 1 , u 2 = 0, L/2 and v 1 , v 2 = L/4, 3L/4 where L is the linear size of the chain.Fig. 10 compares the results from DMRG to Eq. G.2 with n = 2 and with the modification given in Eq.G.4.Indeed, we find good agreement between DMRG and the CFT prediction for sufficiently large L.This justifies the use of Eq.G.2 when calculating the entanglement for the snake phase.

G.2 Averaging over snake configurations
Since the σ degrees of freedom are classical, we are only interested in the quantum entanglement between τ degrees of freedom.We decide to average the value of the quantum entanglement over the thermally-generated σ configurations.This provides us with the "quantum" contribution to the von Neumann entropy of the thermal density matrix after performing a partial trace over the τ degrees of freedom over the complement of the A subsystem, as we now show.
The density matrix of the whole system is block diagonalized in the basis of domain wall configurations.
where Z is the partition function, is the thermal probability of being in the i th domain wall configuration and is the density matrix of the system in that configuration.By taking the partial trace over the τ spins in the subsystem B (complement of A), the density matrix of the snake in subsystem A can be expressed as The von Neumann entropy is obtained by the relation S A = −Tr(ρ A log ρ A ) which gives where S th = − i p i log p i is the thermal entropy of the domain wall configurations (σ spins) and S qu = − i p i Tr ρ A i log ρ A i is the quantum entanglement of the τ spins, thermally averaged over snake configurations.The latter term is the quantum contribution which we are interested in.We simply call it the entanglement entropy and denote it by S in the main text.
Practically, we calculate the average entanglement entropy over an ensemble of domain wall configurations generated by the Monte Carlo method discussed in Appendix.F. In order to find the entanglement entropy of each configuration, we add contributions of different domain walls calculated via Eq.G.3 and Eq.G.4.The only subtlety is that the definition of u i and v i is ambiguous on a lattice and needs to be regularized.We conventionally choose the following UV-regularization: first we number the sites on the loop from 1 to L. Then, for the i th interval, u i is defined as the first site of the interval placed inside the region A, and v i is defined as the first site which lies outside A. We also note that although in our calculations the temperature is not exactly zero, we use Eq.G.2 and Eq.G.4 and drop the finite-T contribution which is parametrically small in system size and is not desired since we seek quantum contributions.

G.3 Analytical prediction for the entanglement scaling
In this section, we derive a formula for the entanglement of a strip based on a Brownian motion model for the interval distribution (see Fig. 5 for the geometry of the strip and the definition of intervals t i ).As mentioned in the main text, since we only aim to reproduce the dominant scaling behavior of the entanglement with the partition size ℓ A , we assume it is sufficient to replace Eq.G.3 by S = (c/3) N int i=1 log(t i ), where t i is the length of interval i.The average over snake configurations is then replaced by an average over interval length f (t) ≡ d t p(t) f (t) with distribution p(t), and the entanglement reads S = (c/3)N int log(t) = (c/3)(2ℓ A /3)log(t).
There remains to find p(t).In order to do this, let us neglect the correlations inherent to an n = 0 fully-packed self-avoiding walk and assume we are dealing with a "plain" random walk instead.The justification is that the geometrical exponent ν = 1/2 for the n = 0 fully-packed self-avoiding walk is the same as that of a plain random walk.It might therefore be safe to neglect correlations and work with a random walk, as far as the qualitative behavior of p(t) is concerned.As a reminder, the exponent ν relates the average end-to-end distance |x| of a walk after t steps according to |x| ∼ t ν .
Using this random walk approximation, t becomes the number of steps a random walk spends in partition A before exiting, having started inside A one lattice site away from the entanglement cut.This is known as a first passage time and was studied in the literature in a variety of geometries [84].Let us consider an infinite cylinder of circumference ℓ and calculate the bipartite entanglement for a strip of width x.Taking the continuum limit, the random walk is described by a Brownian particle with probability distribution U(X , Y, t):  This means the number of loop strands crossing the entanglement cut is independent of x for any given loop configuration.We also know that the number of intervals is simply the number of loop strands crossing the entanglement cuts divided by two.The snake configurations can thus be divided into different sectors based on their value for N int .(For a discussion of these sectors, see [56]).Finally, we know that N int = 2ℓ/3 (up to fluctuations which vanish in the thermodynamic limit) since the N int = 2ℓ/3 sector has the largest entropy [56].All in all, this means we have dropped a contribution to S which is equal to (2ℓ/3)C ′ 1 , which can be absorbed in the constant A in Eq.G.15.

G.4 More numerical results for the entanglement
More numerical results on the entanglement are provided in this appendix.First, Fig. S6 (a) explicitly illustrates the dominant area law scaling of the entanglement for the strip geometry.The dependence of the slope on x reveals the non-local nature of the entanglement, as discussed in the main text.
Secondly, we show in Fig. S6 (b) the fitting parameters obtained by fitting Eq.G.15 to our numerical results, at fixed ℓ (see Fig. 6 in the main text for the fits).We denote the potentially ℓ-dependent values obtained by these fits A(ℓ) and B(ℓ).We find a small drift of the fitting parameter B with ℓ, and almost no drift for A. Based on Fig. S6 (b), we propose a fit of the form x , (G. 16) Figure1: Example of a decorated domain wall configuration.σ z = ±1 spins are represented by filled and empty blue dots, respectively.τ spins (red dots) are located on every site of the hexagonal lattice (dashed blue lines), including the domain walls of the σ z spins (green lines).As an example, τ 1 and τ 2 are located on the domain wall between sites p and q. Figure a) illustrates a snake configuration, consisting of a single domain wall traversing the entire lattice.Figures b)and c) illustrate hexagon solid and stripe configurations, consisting of contractible and non-contractible domain walls of the shortest length, respectively.

Figure 5 :
Figure 5: Top: Partition of an ℓ by ℓ torus for which we calculate entanglement.Bottom: After unfolding the snake, the partition becomes a union of disjoint intervals.Intervals belonging to the subsystem A are colored while those outside A are black and labeled B.

Figure 6 :
Figure6: Universal contributions to the entanglement entropy of the strip shown in Fig.5vs.x for different system sizes ℓ.The black curves are fitted to x ≥ 4 numerical data points according to Eq. 9.The snake configurations were generated at T = 4. Details on how S was calculated numerically, and about the non-universal contributions which were dropped, can be found in Appendix G.

Figure 7 :
Figure 7: Example of a decorated domain wall configuration.σ z = ±1 spins live on a triangular lattice and are represented by filled and empty blue dots, respectively.τ spins (red dots) are located on the vertices of the hexagonal lattice dual to the triangular lattice of the σ spins.Domain walls of σ z spins are depicted as green lines.An example of a triangle of σ z spins △pqr and the corresponding τ spins is also shown in reference to Eq. A.1.

2 Figure 8 :
Figure 8: Energy density of different 1D CFT chains with an even-odd effect.The energies are calculated through exact diagonalization.(a) Energy density of the spin-1 chain introduced in C.1 at D = 1 ≈ D c , (b-d) energy density of the X X Z chain at different values of ∆, namely ∆ = −0.5, 0, 0.5.It can be numerically confirmed that the system exhibits the even-odd effect at any values of −1 < ∆ 1, where the system is in its gapless phase.However, the values of c 0 and c 1 are not universal and depend on ∆.

Figure 9 :
Figure 9: Color map of C Z Z (k) for an ℓ = 48 system at T = 4 which peaks at the corners of the Brillouin zone.

Figure 10 :
Figure10: Second Renyi entropy for four equal partitions of the X − Z X Z chain with linear size L, computed according to the DMRG method in Ref.[60] (with maximum bond dimension χ = 100) and the analytical prediction for a CFT given in Eq.G.2.The entropies computed using DMRG are shifted by a constant value of −1.889 to account for the difference generated by the non-universal constant c 2 and scaling function F N ,2 (x) given in Equation G.1.

Figure 12 :
Figure 12: (a) Entanglement entropy of the strip shown in Fig. 5 vs. ℓ for fixed values of its width x.The graph shows that the entanglement entropy grows almost linearly, indicating that the leading contribution to the entropy follows the area law.(b) A(ℓ) and B(ℓ) along with their respective linear fits vs. 1/ℓ.