Measuring von Neumann entanglement entropies without wave functions

TL;DR

This paper introduces a novel method to measure von Neumann entanglement entropy in quantum many-body systems without needing the wave function, using thermodynamic properties of entanglement Hamiltonians applicable to simulations and experiments.

Contribution

The authors develop a thermodynamic approach to determine entanglement entropy directly, bypassing the need for wave function access, and validate it on various quantum systems.

Findings

Successfully measured entanglement entropy in quantum spin chains and 2D magnets.

Identified the onset of area law behavior and number of Goldstone bosons.

Confirmed a conjecture on geometric entanglement at critical points.

Abstract

We present a method to measure the von Neumann entanglement entropy of ground states of quantum many-body systems which does not require access to the system wave function. The technique is based on a direct thermodynamic study of entanglement Hamiltonians, whose functional form is available from field theoretical insights. The method is applicable to classical simulations such as quantum Monte Carlo methods, and to experiments that allow for thermodynamic measurements such as the density of states, accessible via quantum quenches. We benchmark our technique on critical quantum spin chains, and apply it to several two-dimensional quantum magnets, where we are able to unambiguously determine the onset of area law in the entanglement entropy, the number of Goldstone bosons, and to check a recent conjecture on geometric entanglement contribution at critical points described by strongly coupled field theories.

Equations

$$S_A=-{\text{Tr}}_A{\rho }_A\ln {\rho }_A,{\rho }_A={\text{Tr}}_B|\Psi ⟩⟨\Psi |.$$

$${\tilde{H}}_A=2\pi {\int }_{\vec{x}\in A}d\vec{x}\left[x_{1}H(\vec{x})\right]+c^{\prime },$$

$${\tilde{H}}_{EH}={\beta }_{EH}\left[\sum _{n=1}^{L-1}\Gamma (n)h_{n,n+1}+\Gamma (n-1/2)l_n\right],$$

$${\rho }_{BW}=\frac{e^{-{\beta }_{EH}{H}_{EH}}}{Z_{EH}},S_{BW}=-\text{Tr}{\rho }_{BW}\ln {\rho }_{BW}$$

$$H_{BW}=\sum _i\Gamma (i){\vec{S}}_i{\vec{S}}_{i+1},$$

$$H=-\sum _i\Gamma (i)S_i^z{S}_{i+1}^z-\sum _i\Gamma (i-1/2)S_i^x.$$

$$\Gamma (x)=\frac{L}{\pi }\sin \left(\frac{\pi x}{L}\right),$$

$$\Gamma (x)=\frac{2L}{\pi }\sin \left(\frac{\pi x}{2L}\right).$$

$$H_{\text{BW}}$$

$$\Gamma \left(x\right)=x,$$

$$\Gamma \left(x\right)=\frac{x(L-x)}{L},$$

$$Z=\mathrm{Tr}{e}^{-\beta H}=\sum _{n=0}^{\infty }\frac{{\beta }^n}{n!}\mathrm{Tr}{\left(-H\right)}^n=\sum _{n=0}^{\infty }{\beta }^{n}g(n),$$

$$S(\beta )=\frac{1}{N_r}\sum _{i=1}^{N_r}{S}_i(\beta ).$$

Full text

Introduction. -

Over the last twenty years, entanglement has emerged as a paramount tool to characterize quantum wave-functions Amico et al. (2008); Fradkin (2013); Eisert et al. (2010); Laflorencie (2016). A striking example is ground states $|\Psi\rangle$ of many-body systems where, given a spatial bipartition dividing the system into regions $A$ and $B$, the entanglement between $A$ and $B$ is measured by the von Neumann entropy (VNE):

[TABLE]

The VNE remarkably provides a systematic way to connect wave function properties to all operational definitions of entanglement, and is of pivotal importance to both quantum information purposes and as a diagnostic tool in quantum many-body theory. Examples of its relevance include the existence of area laws bounding entanglement in ground state of local Hamiltonians Eisert et al. (2010), the sharp characterization of conformal field theories (CFTs) in one-dimension (1D) Holzhey et al. (1994); Calabrese and Cardy (2004, 2009), topological order Kitaev and Preskill (2006); Levin and Wen (2006) and spontaneous symmetry breaking Metlitski and Grover (2011), and its importance in understanding the complexity of classical simulations Schollwöck (2005). However, differently from remarkable theoretical Moura Alves and Jaksch (2004); Daley et al. (2012); Alba (2017); Elben et al. (2018); Vermersch et al. (2018) and experimental Islam et al. (2015); Brydges et al. (2018) studies aimed at measuring Renyi entropies, the VNE has so far eluded a direct experimental verification, as it requires tomographic access to the system wave function - which becomes quickly impractical beyond few spins. The same limitations affect numerical methods such as quantum Monte Carlo (QMC) Kaul et al. (2013), which typically cannot sample wave functions, and, thus, cannot compute the VNE.

In this work, we describe an approach to measure the von Neumann entanglement entropy of ground states without relying on probing wave functions, that can be efficiently implemented in both experiments and QMC simulations. The backbone of the technique is the formulation of the entanglement measurement problem in terms of the thermodynamic properties of the entanglement (modular) Hamiltonian (EH) $\tilde{H}_{A}=-\ln\rho_{A}$ Bisognano and Wichmann (1975, 1976); Witten (2018), whose structure is schematically depicted in Fig. 1. As we show below, our method allows to perform accurate entanglement-based measurements of universal quantitites, such as the number of Nambu-Goldstone modes Metlitski and Grover (2011) and central charges Holzhey et al. (1994); Calabrese and Cardy (2004), at the percent level, even for modest system sizes. Most remarkably, it allows the calculation of the entanglement of many-body systems in a scalable manner (and well beyond what can be done with alternative numerical methods), thanks to its thermodynamic analogy. This is a key point when interested in universal quantities, as those are captured by subleading corrections to the entropy in dimensions $D>1$. In terms of experiments, our work is complementary to other techniques that have been proposed to measure either Renyi entropies Moura Alves and Jaksch (2004); Daley et al. (2012); Islam et al. (2015); Brydges et al. (2018); Elben et al. (2018); Vermersch et al. (2018); Abanin and Demler (2012) and entanglement spectra Pichler et al. (2016); Dalmonte et al. (2018), both in terms of observables and in terms of resources, as we detail below.

After benchmarking our method on 1D examples, we carry out QMC simulations on a series of two-dimensional lattice models and calculate their entanglement entropy for partitions of up to thousand sites - inaccessible to any other controlled numerical method. For the 2D Heisenberg and XY models, we provide direct evidence that (i) the VNE is constrained by the area law (confirming lower bounds based on Renyi entropies), and (ii) the number of Goldstone modes can be determined with percent accuracy solely from entanglement properties. For the bilayer Heisenberg model, we study the geometric contribution to the entanglement entropy at its strongly coupled critical point, and verify a recent conjecture on $O(N)$ models Metlitski et al. (2009).

Thermodynamics of entanglement Hamiltonians. -

The relation between entanglement and thermodynamic quantities has been widely exploited in the quantum mechanics and field theory literature: an epitome in this context is the Unruh effect Unruh (1976), that describes how the vacuum appears as an equilibrium finite temperature state from the point of view of an accelerating observer. In the context of axiomatic field theory, this relation is conveniently expressed by the Bisognano-Wichmann (BW) theorem Bisognano and Wichmann (1975, 1976); Witten (2018). For a Lorentz invariant theory with Hamiltonian density $H(\vec{x})$, $\vec{x}=(x_{1},...,x_{D})$, in $D$ spatial dimensions, the entanglement Hamiltonian of a half-plane bipartition A defined by $x_{1}>0$ reads:

[TABLE]

where $c^{\prime}$ is a constant that ensures $\text{Tr}_{A}\rho_{A}=1$. The BW theorem has been extended to different geometrical partitions in the presence of conformal invariance Hislop and Longo (1982); Cardy and Tonni (2016); Casini et al. (2011).

These results can be cast on a discrete space-time lattice Dalmonte et al. (2018); Giudici et al. (2018) as follows. For the sake of simplicity, let us focus on 1D systems with nearest-neighbor interaction, $h_{n,n+1}$, and on-site terms, $l_{n}$; the 2D case is discussed in the SM. Up to $c^{\prime}$, the lattice BW-EH ansatz of a subsytem of length $L$ is

[TABLE]

The coefficients $\Gamma$ depend on the geometry of the partitionBisognano and Wichmann (1975); Hislop and Longo (1982); Cardy and Tonni (2016); Giudici et al. (2018): (i) for a half-infinite partition under open boundary conditions (OBC), $\Gamma(n)=n$ [see Fig. 1 (b2)]; (ii) for subsystem embedded in an infinite system, $\Gamma(n)=n\left(L-n\right)/L$ Hislop and Longo (1982), which corresponds to periodic boundary condition (infinite PBC); and for finite systems (iii) with both PBC, $\Gamma(n)=\frac{L}{\pi}\sin\left(\frac{\pi x}{L}\right)$, and (iv) OBC, $\Gamma(x)=\frac{2L}{\pi}\sin\left(\frac{\pi x}{2L}\right)$. It is straightforward to generalize the BW-EH ansatz for a N-dimensional lattice model Giudici et al. (2018): in Fig. 1b, we schematically illustrate it for the cylinder geometries discussed below; see the SM for more details sup .

The overall energy scale of Eq. (3) is related to the “speed of light“, $v$, in the corresponding low-energy field theory, $\beta_{EH}=\frac{2\pi}{v}$ and plays the role of an effective inverse temperature. The corresponding reduced density matrix on the lattice and the corresponding thermal entropy read:

[TABLE]

where the normalization factor is interpreted as a partition function $Z_{EH}=\text{Tr}e^{-\beta_{EH}H_{EH}}$ 111We note that, while in principle the entropy is nothing but the expectation value of the EH, the normalization factor $Z_{EH}$, which is not universal, makes this approach hardly applicable..

The predictive power of the BW theorem on the lattice goes well beyond the low-lying entanglement spectrum, for which analytical and numerical evidence is abundant Peschel et al. (1999); Peschel and Eisler (2009); Itoyama and Thacker (1987); Nienhuis et al. (2009); Kim et al. (2016); Eisler and Peschel (2017); Dalmonte et al. (2018); Eisler and Peschel (2018); Zhu et al. ; Parisen Toldin and Assaad (2018); Kosior et al. (2018); Giudici et al. (2018): in particular, the entanglement entropy of a lattice system $S_{A}$ can be evaluated via its thermal equivalent $S_{BW}$ so to determine universal quantities at the percent level even from modest system sizes. The accuracy progressively increases with system size, as expected for any field theoretical expectation on the lattice. We remark that the validity of Eq. (2) relies on the underlying field theory being Lorentz invariant: as such, it is applicable to a broad range of phenomena, including quantum critical phases and points with emergent relativistic description, topological phases, lattice gauge theories, just to name few examples; it cannot be applied to other situations, such as disordered systems or ferromagnets, where low-energy Lorentz invariance is absent.

Measuring entanglement entropy in numerical simulations and experiments. -

Before discussing the concrete validation examples, we illustrate how to measure VNE in both numerical and real experiments. Both strategies rely on the basic fact that, in analogy with thermodynamics, the VNE is uniquely determined once the density of states of the EH, $D_{EH}$, is known.

From the point of view of classical computations, this can be achieved using any algorithm based on metadynamics. Below, we illustrate this by applying the quantum version of the Wang-Landau method performed in the stochastic series expansion (SSE) QMC framework Troyer et al. (2003); Wessel et al. (2007). This method allows the direct calculation of the free-energy and the entropy of the BW-EH; see the SM for more details. One key point we emphasize is that the method straightforward to implement on a working Wang-Landau code [only requires to implement an inhomogeneous version of the system Hamiltonian, as Eq. (3)] cod .

From an experimental viewpoint, $D_{EH}$ can be conveniently measured following the procedures proposed in Refs. Osborne ; Schrodi et al. (2017) in the context of conventional Hamiltonian operators. The main idea here is that the density of states of a given operator $A$ can be obtained by analysing the quench dynamics starting from a randomized set of product states, and evolving according to $A$. We remark that certain universal quantities are then simply accessible via the density of states itself, without the need to evaluate the VNE: an example is the central charge in one-dimensional critical systems Calabrese and Lefevre (2008).

A second way of measuring the VNE experimentally is via the specific heat with respect to the EH: in analogy with thermodynamics, this quantity, when integrated over the temperature range $[1/\beta_{EH},\infty]$ returns the von Neumann entropy. These temperatures are in units of the speed of sound, which can be measured independently via spectroscopy. However, this approach is limited to systems where a controlled exchange with the reservoir is possible, e.g., cold atoms in optical lattices Catani et al. (2009).

One-dimensional critical systems. -

We now benchmark our strategy for one-dimensional critical systems, where the calculation of the VNE is amenable to both exact and tensor network simulations. In this case, the VNE of a subsystem of size $L$ diverges logarithmically, $S(L)\propto c\ln L$, where $c$ is the central charge of the underlying CFT.

In Fig. 2, we plot the BW VNE of the one-dimensional Heisenberg model (HM) and the quantum Ising model (QIM) at its quantum critical point, and under PBCs. Throughout this work, we employ dimensionless energy units for the sake of convenience. For the two models, the exact value of the entropy (empty circles) is evaluated using density-matrix-renormalization-group White (1992) (HM) and exact diagonalization methods for a biparition of size $L$ embedded in systems of size $2L$. The calculations of the BW-EH thermal entropy are carried out with QMC with both local and SSE directed-loop updates Sandvik and Kurkijärvi (1991); Syljuåsen and Sandvik (2002) for the HM, and exact diagonalization for the QIM. In addition to the finite-size EH (red triangles), for the sake of comparison, we also compute the entropy obtained utilizing the EH of a finite partition in an infinite system (black circles) Giudici et al. (2018): the two are separated only by a constant shift that depends solely on the central charge sup .

In both cases, the VNE increases logarithmically as expected: the corresponding central charge considering systems up to $L=80$ $(100)$ is in within 1% (0.05%) level with the exact results for the HM (QIM) - see Figs. 2 (a) and (b). The difference $\Delta S(L)=S_{BW}-S_{exact}$ goes to zero as $L\to\infty$. Remarkably, the BW-EH captures also corrections to the CFT scaling due to marginal operators Laflorencie et al. (2006); Calabrese et al. (2010) under OBCs, as we discuss in the SM. These results strongly suggest that the EH-BW entropy provides the exact VNE in the limit $L\to\infty$ for one-dimensional critical systems, capturing both universal (leading contribution) and non-universal (constant) contributions.

Two-dimensional quantum magnets.

The VNE also describes universal properties of two-dimensional systems. For instance, the VNE of 2D ground states that break a continuous symmetry scales as $S(L)=AL+B\ln(L)+D$, where $L$ is the linear size of the boundary. The $A$ is the non-universal area law term Eisert et al. (2010), while, for a smooth boundary, the prefactor of the logarithmic term is a universal quantity related to the number of Nambu-Goldstone modes $n_{b}$, $B=n_{b}/2$, of the associated spontaneuosly-symmetry-broken (SSB) phase Kallin et al. (2011); Metlitski and Grover (2011). As examples of SSB, we consider the 2D XY model and the Heisenberg model. In both cases, we perform QMC simulations of the EH and extract the corresponding VNE as a function of the subsystem linear size, $L$. The entropy is evaluated at $\beta_{EH}=2\pi/v$, with $v_{Heis}=1.658J$ Sen et al. (2015) and $v_{XY}=1.134J$ Sandvik and Hamer (1999), using the SSE-Wang Landau algorithm.

In Fig. 3a, we show the scaling of the BW VNE for both cylinder and torus geometries. The scaling is clearly linear. In the case of the HM on a torus, we extracted the coefficient $A$ by fitting these results to $S(L)=AL+B\ln(L)+D$, and obtain $A=0.372(6)$, which is in agreement with a prediction based on spin-wave approximation Song et al. (2011) (discrepancy $<3\%$).

In Fig. 3b, we extract the subleading logarithmic correction by considering the entropy difference $2S(L)-S(2L)\simeq\frac{n_{b}\log(L)}{2}$ in toroidal geometries of circumference $2L$. The number of Nambu-Goldstone modes obtained from the prefactor of this term is in perfect agreement with field theoretical expectations Metlitski and Grover (2011); Song et al. (2011); Humeniuk and Roscilde (2012); Kulchytskyy et al. (2015), with accuracy at the percent level or lower. The fact that the VNE returns a value which is considerably closer to the field theoretical prediction when compared to the one extracted from Renyi entropies Kallin et al. (2011); Humeniuk and Roscilde (2012) may signal the fact that the latter are more affected by irrelevant operators, as observed in 1D Laflorencie et al. (2006); Calabrese et al. (2010); Dalmonte et al. (2011), or may be due to the smoother continuity properties of the VNE.

Strongly coupled Quantum criticality. -

As a second example of 2D system, we consider the bilayer Heisenberg model Sandvik and Scalapino (1994); Wang et al. (2006). This model describes a quantum phase transition induced by the inter-coupling $g$ that belongs to the $O(3)$ universality class. We compute the BW-EH entropy at the QCP, $g_{c}=2.522$, considering $\beta_{EH}=2\pi/v$, with $v=1.9001(2)$ Sen et al. (2015). The details of how we cast the BW-EH in the bilayer HM are discussed in Ref. sup . For this universal class, it has been argued that there is a universal constant correction to the entanglement entropy that depends solely on the aspect ratio Metlitski et al. (2009): for PBC, this constant has been conjectured to vanish, in sharp contrast to anti-PBC. Verifying this conjecture requires accurate values of the entropy at large system sizes of several hundred sites.

Our results up to partition of size $L=18$ are depicted in Fig. 3c. Within error bars, our results show that $S(L)$ is independent of the aspect ratio of the subsytem, see Fig. 3 (c), and have no detectable logarithmic subleading term (the $S(L)=AL+B\ln(L)+D$ fitting, gives $B=-0.05(8)$). The $y$-intercept of $S(L)$ is $0.010(7)$, which confirms the conjecture for the universal constant contribution for the $O(N)$ model Metlitski et al. (2009).

Stability of BW-EH entropy. -

Finally, we discuss how stable are the BW-EH entropy to perturbations. The two main effects we consider are (i) errors in the value of $\beta_{EH}$, and (ii) random perturbations in the EH couplings $\Gamma(n)$. The first is relevant for both QMC simulations and experiments, where in general the value of $\beta_{EH}$ is not known exactly, while the second accounts for possible imperfect experimental realizations of the EH.

In Fig.4 (a,b1), we show the value of the extracted entropy obtained via Wang-Landau sampling as a function of $\beta$, for both 1D and 2D HM. The insets magnify the region in the vicinity of the exact value of $\beta_{EH}$, signaled by a dashed vertical line: in this regime, the entropy is linearly sensitive to $\beta$. This implies that the accuracy in estimating $S$ is ultimately limited by the accuracy on the sound velocity: this strengthen the applicability of our method to QMC simulations, where $v$ can be measured very accurately via a variety of techniques Kaul et al. (2013); Sen et al. (2015).

Next, we consider the effect of disordered couplings in the EH, $\Gamma(n)\to\Gamma(n)(1+\delta_{n})$, where $\delta_{n}$ is a random number in the interval $[-\delta,\delta]$ 222We remind that this imperfect EH corresponds to the GS of a clean system, and is not related to the entanglement properties of disordered systems.. In In Fig.4 (b2), we plot the average value of the entropy as a function of $\delta$ for the 2D HM (the 1D case is discussed in Ref. sup ). The average value of the entropy is very robust to disorder strengths of magnitude 10%: this remarkable stability is in contrast to what is typically found when studying the effects of disorder in the Hamiltonian couplings, which have a quantitatively larger effect on entropies. A possible element in support of this unexpected resilience is the fact that the VNE is endowed with particularly robust continuity properties with respect to changes in the entanglement spectrum (which is instead expected to be directly affected by the random couplings).

Conclusions.

We have presented a method to measure the ground state von Neumann entropy of a broad class of lattice models via direct thermodynamic probe of the correspondent entanglement Hamiltonian. The method is straightforward to implement in quantum Monte Carlo codes, and is of immediate applicability to experiments capable of measuring the density of states. It enables accurate predictions of universal quantities solely based on entanglement, thanks in particular to its immediate scalability in numerical simulations. Future perspectives include the application of the method to other entanglement related quantities, such as the negativity, its extension to lattice gauge theories, and its integration with methods to determine the EH at finite temperature Parisen Toldin and Assaad (2018); Frerot (2017).

I Entanglement entropy of one-dimensional critical systems

I.1 Entanglement entropy of Heisenberg and Ising models on open chains

In this section, we provide additional details about the calculation of the von Neumann entropy (VNE) in one-dimensional critical systems. As discussed in the main text, the VNE is related to the thermal entropy of the lattice Bisognano-Wichmann entanglement Hamiltonian (BW-EH). For the Heisenberg model (HM) one has

[TABLE]

and for the quantum Ising model (QIM) at the critical point

[TABLE]

The low energy properties of these lattice models are described by Lorentz-invariant quantum field theories, a prerequisite for the applicability of our approach. Additionally, the emergence of conformal symmetries in these systems allows us to consider the extensions of the BW theorem to different geometrical partitions Hislop and Longo (1982); Cardy and Tonni (2016), see Fig.5 (a). For the half-partition of a ring with length $2L$, for instance, one has

[TABLE]

while for the open chain

[TABLE]

The definition of the BW-EH for finite systems allows us to directly compare the BW VNE with the exact results of the VNE obtained with exact diagonalization or the density-matrix-renormalization-group. In the PBC case, for instance, the BW entropy is in perfect agreement with exact results, see Fig. 2 of the main text. The discrepancy $\Delta S(L)=|S_{BW}-S_{exact}|$ goes to zero in the limit $L\to\infty$. Furthermore, the corresponding central charge considering system sizes up to $L=80$ ($L=100$) is only $1\%$ ($0.05\%$) away from the exact results for the HM (QIM). It is interesting to note that, when we directly compare the absolute value of the BW-EH entropy for the infinite PBC with the exact VNE (subsystem with size $L$ embedded in a system of size $2L$), we observe a shift $\epsilon$. This shift is almost equal to $\epsilon=\left(c/3\right)\ln(\pi/2)$, as expected from the finite-size corrections given by the CFT expression of the VNE Calabrese and Cardy (2004, 2009).

For the OBC case, we observe an alternating term of the BW-EH entropy for the HM, but not for the QIM, see Figs. 5 (b) and (c). These results is in agreement with the exact VNE. As discussed in Ref. Laflorencie et al. (2006); Calabrese et al. (2010), those oscillations are universal and due to the antiferromagnetic nature of the interactions, not appearing in the QIM Zhou et al. (2006) (in the latter, the effective Fermi momentum is either 0 or $\pi$). From the CFT perspective, the oscillations can be viewed as lattice corrections of scaling dimension $\Delta_{p}$: their decay as a function of the bipartition size is a power law whose exponent is related to $\Delta_{p}$  Calabrese et al. (2010); Dalmonte et al. (2011). The fact that the BW-EH faithfully reproduces not only the leading, but also the dominant subleading correction testifies its predictive power on generic universal quantities captured by the VNE (a CFT-specific analysis will be reported elsewhere Mendes-Santos et al. ). While, for instance, non-universal contributions such as additive constants in 1D shall not be immediately reproduced due to the field theoretical origin of the relation we employ, in all examples where a comparison to exact results is possible (essentially, 1D systems), we observe that even non-universal contributions are accurately captured: for instance, $\Delta S(L)$ goes to zero in the limit $L\to\infty$ both in the OBC and PBC cases. We attribute this to the fact that the BW-EH is actually able to reproduce a “partition function” whose corresponding Hamiltonian has the correct density of states, and whose generic correlation functions are correct Giudici et al. (2018). In case only the first element was true, and, for instance, the overall scaling correction was wrong, one would have generically expected incorrect correlation functions. From a methodological viewpoint, this implies that our method may be used to check convergence of tensor network states in conformal phases, especially for large values of the central charge.

I.2 Resilience to errors: the 1D Heisenberg model

We now consider the effect of disordered couplings, $\Gamma(n)\to\Gamma(n)(1+\delta_{n})$, where $\delta_{n}=[-\delta,\delta]$, in the BW-EH of the 1D HM. Specifically, we are interested in understanding how the BW VNE is affected by a small amount of disorder, which is an important issue in an experimental context.

In Fig. 6 we show that, similarly to the 2D case [see Fig. 4 (b1) of the main text], the BW VNE is not appreciably affected by disorder up to strength of the order of $10\%$. For larger values of $\delta$, we observe a considerable dependence on the disorder realization, as signalled by the visually large spreading of the values of $S$. Surprisingly, the mean value of the entropy is not dramatically affected.

II Bisognano-Wichmann entanglement Hamiltonian for 2D systems

We now review how to cast the BW-EH on two dimensional lattices Giudici et al. (2018). As a concrete example, we consider the 2D Heisenberg model in a square lattice $L_{x}\times L_{y}$. In this case, the BW-EH is

[TABLE]

where the lattice spacing has been set to 1 without loss of generality. The simulation of the subsystem BW-EH is performed considering periodic boundary condition in the $y$ direction, and open boundary condition in the $x$ direction, see Fig. 7.

The function $\Gamma(x)$ is given by the BW theorem

[TABLE]

which represents the EH of a half-bipartition; we call this subsystem-geometry of cylinder. Furthermore, we consider

[TABLE]

which corresponds to the generalization of the BW to a subsytem that is embbeded in a infinite system; we call this subsystem-geometry of toroid.

We remind the reader that, as discussed in Ref. Giudici et al. (2018), for finite values of $L_{y}$, formula Eq. (11) is in principle only applicable to conformal field theories. In order to mitigate this effect, we perform systematically extrapolations to large $L_{y}$ when determining universal properties. Let us illustrate here a simple, non-rigorous argument that partly justifies the applicability of this approach to generic (i.e., non conformal) 2D models. Typically, the low energy theory will be made of gapless and gapped sectors. The description of the former will be scale invariant and relativistic invariant: while this does not guarantee emergent conformal invariance, exceptions are rare. The gapped part of the theory will (at most) contribute to the entanglement properties only in the very vicinity of the edge of the partition, where it would actually behave like a gapless theory. Far from the boundary, the reduced density matrix with respect to these degrees of freedom will be an identity operator (up to degeneracies). This indicates that the CFT formulas used above shall be applicable also to more general cases where some low-energy degrees of freedom are actually gapped. In the context of the 2D HM, the role of gapless degrees of freedom is played by the $\mathbb{C}P(1)$ model describing the emergent Nambu-Goldstone modes, and the gapped part of the theory is described by the massive Goldstone mode.

III Quantum Wang-Landau sampling of the entanglement Hamiltonian

In this section we discuss some relevant details of the quantum Wang-Landau (WL) simulations used in this work. Compared with the convetional quantum Monte Carlo (QMC) simulations, that is performed at a fixed temperature, the WL method features two main advantages for the study of the thermodynamic properties of the EH: (i) it allows to directly compute the thermal entropy at the “entanglement temperature” $\beta_{EH}$, and (ii) the thermodynamic properties of the EH are obtained for a broad range of temperature with a single run of the simulation.

The WL method was originally proposed for classical systems in Ref. Wang and Landau (2001). The key idea of the method is to calculate the density of states, $\rho(E)$. For a quantum Hamiltonian, such as the BW-EH, however, one must map the system to a classical one. This is done, for instance, using the SSE framework, which considers the following form for the partition function

[TABLE]

where the nth order series coefficient $g(n)$ plays the role of the density of states in the classical algorithm. We refer to Ref. Troyer et al. (2003); Wessel et al. (2007) for the general details of the computation of $g(n)$. Below we mention the technical aspects of the simulation that are relevant to reproduce our results.

The series expansion Eq.(12) can be truncated at an order $\Lambda$, i.e., $n=0,1,...,\Lambda$, without introducing systematic errors in the simulation. The choice of $\Lambda$ is performed as in the conventional SSE simulations, see Refs. Sandvik and Kurkijärvi (1991); Syljuåsen and Sandvik (2002), which gives as a result $\Lambda(\beta)\approx\beta|E(\beta)|$; where $E(\beta)$ is the expectation value of the total energy at inverse temperature $\beta$. The effect of the cutoff $\Lambda(\beta)$ is that the range of temperature that can be accessed is restricted to $\beta<\Lambda(\beta)/E(\beta)$. In order to obtain the results of Fig. 4(a,b1) of the main text, for instance, we simulate the BW-EH using $\Lambda(3\beta_{EH})$ as the cutoff. Instead, the computation of the BW VNE at $\beta_{EH}$ are obtained utilizing a cutoff $\Lambda(\alpha\beta_{EH})$. The results shown in Figs. 2 and 3 of the main text are obtained with $\alpha=1.3$. We check that these results don’t change upon increasing $\alpha$.

In the WL-SSE algorithm, the sampling of the SSE configurations with different n is performed with a probability function that is proportional to the inverse of the “density of states“, $1/g(n)$. The WL sampling generates a histogram for the distribution of n that is flat, i.e., $H(n)\sim const$; the histogram $H(n)$ is obtained counting the number of times a configuration with n is observed. The key point of the algorithm is that $g(n)$ can be computed by iteratively flattening $H(n)$. More specifically, one start with the guess $g(n)=1$. Further, each time the configuration n is accepted $g(n)$ is multiplied by a factor $f$, i.e., $g(n)\to g_{old}(n)f$. This proccess is repeated until $H(n)$ is flat. In practice, we consider as a condition for the flatness of $H(n)$ a maximum deviation of 20$\%$ from the mean value. Once $H(n)$ is flat, it is reset to zero, and $f$ is decreased by $\ln(f)\to\ln(f_{old})/2$ Inglis and Melko (2013). This proccess is repeated until convergence is achieved. Here we use the convergence condition proposed in Refs. Inglis and Melko (2013); Belardinelli et al. (2008).

In addition to the aforementioned algorithm, we consider the optimized-broad-histogram algorithm proposed in Ref. Wessel et al. (2007) for the 2D Heisenberg model, see Fig.3 (a) of the main text. These results were obtained with the ALPS code Bauer et al. (2011); cod In this case, we confirm that the two methods give the same results (within error bars).

Finally, it is important to mention that the results of the BW entropy are obained by doing an average of $N_{r}$ independent WL simulations, i.e.,

[TABLE]

The error bars are the standard deviation of the distribution $\{S_{i}\}$, and for all the results presented, we consider at least $N_{r}>200$.