Measuring von Neumann entanglement entropies without wave functions

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 lattice entanglement Hamiltonians—recently proposed in the paper [Dalmonte et al 2018 Nat. Phys. 14 827] via field theoretical insights—and can be performed by quantum Monte Carlo methods. 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. The protocol can also be adapted to measure entanglement in experiments via quantum quenches.


Introduction
Over the last twenty years, entanglement has emerged as a paramount tool to characterize quantum wave functions [1][2][3][4]. A striking example is ground states |Yñ 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): The VNE remarkably provides a systematic way to connect wave function properties to 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 [3], the sharp characterization of conformal field theories (CFTs) in one-dimension (1D) [5][6][7], topological order [8,9] and spontaneous symmetry breaking [10], and its importance in understanding the complexity of classical simulations [11]. Despite its pivotal importance, the current understanding of entanglement measures in many-body systems is essentially limited to non-interacting theories or to lower bounds provided by Renyi entropies, due to the lack of methods to calculate the VNE in any dimension D>1. This represents a key obstacle in determining both the capabilities of many-body systems in terms of quantum information processing (e.g. how much entanglement can be distilled from a given partition), and the generic relation between universal field theoretical descriptions and entanglement.
In this work, we describe an approach to compute the von Neumann entanglement entropy of ground states without relying on probing wave functions, that (i) is applicable in any dimension, and to a broad class of physical phenomena, including quantum critical and topological matter; and (ii) it is amenable to simulations based on Quantum Monte Carlo, and thus scalable to systems sizes order of magnitudes larger than attainable with other methods. 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)˜r = -H ln A A [12][13][14], whose structure is given by the lattice version of the Bisognano-Wichmann (BW) theorem, see the schematic explanation in figure 1. Concerning low-lying entanglement spectra and correlation functions, the applicability of the BW theorem on the lattice has been verified in recent works [15,16]. Here, we take a considerable step forward, and show that it can also be applied to do accurate entanglement-based measurements of universal quantities, such as the number of Nambu-Goldstone modes [10] and central charges [5,6], 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 (well beyond what can be done with alternative numerical methods), thanks to its thermodynamic analogy: this allows us to verify the onset of area law in two-dimensional quantum magnets, up to system sizes including ( )  10 3 spins. Such scalability is a key point when interested in universal quantities, as those are captured by subleading corrections to the entropy in dimensions D>1. We remark here that, in general, entanglement measures are unrelated to the low-lying entanglement spectrum, and do instead critically depend on the distribution of all eigenvalues, a much more delicate quantity to deal with that was never addressed in the context of lattice adaption of the BW theorem. In terms of techniques, our work complements the already successful QMC toolbox to lower-bound many-body entanglement via Renyi entropies [4,[17][18][19][20][21].
After benchmarking our method on 1D examples (an extension of those results is shown in [22]) we carry out QMC simulations on a series of 2D lattice models. For the 2D Heisenberg and XY models, we provide direct evidence that (i) the VNE is constrained by the area law (in agreement with 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 [23].

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 [24], 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 BW theorem [12][13][14]. For a Lorentz invariant theory with Hamiltonian density ( ) , in D spatial dimensions, the entanglement Hamiltonian of a half-plane bipartition A defined by x 1 >0 reads: where c′ is a constant that ensures r = Tr 1 A A . The BW theorem has been extended to different geometrical partitions in the presence of conformal invariance [25][26][27].
These results can be cast on a discrete space-time lattice [15,16] 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  [12,13]. A plane is divided into two half-planes at x=0. The reduced density matrix obtained from the vacuum of the field theory upon tracing the x<0 region can be interpreted as a thermal equilibrium state with inverse temperature β increasing as a function of the distance from the boundary. Hot region (red) are typically more entangled then the cold (blue) ones. Panel (b): adaption to cylinder geometries. In analogy with the infinite plane case, the inverse entanglement temperature is constant at fixed x (path γ y ), while it increases at fixed y (path γ x ). This picture can immediately be adapted to lattices [panel b1)] [15,16]: as depicted in b2), the couplings of the corresponding lattice entanglement Hamiltonian are constant along y (Γ y ), while they increase along x (Γ x ).
latter. Up to ¢ c , the lattice BW-EH ansatz of a subsystem of length L is ⎡ where n is the indices of sites within the subsystem A. The coefficients Γ depend on the geometry of the partition [12,16,25,26] It is straightforward to generalize the BW-EH ansatz for a N-dimensional lattice model [16]: in figure 1(b), we schematically illustrate it for the cylinder geometries discussed below.
The overall energy scale of equation (3) is related to the 'speed of light', v, in the corresponding low-energy field theory, b = p v EH 2 and plays the role of an effective inverse temperature. The lattice density matrix corresponding to the subsystem A and its thermal entropy read: . In the present framework of thermodynamic studies of EH, one needs to understand the predictive power of the BW theorem in relation to entropies, and to identify a proper theoretical framework to compute the VNE from thermodynamics.
The predictive power of the BW theorem on the lattice has been broadly verified regarding the low-lying entanglement spectrum and correlation functions [15,16,[28][29][30][31][32][33][34][35][36][37]. While these results represent a promising first step, they are not informative on the capability of the BW-EH to capture entanglement measures, since: (1) they are limited to some observables and do not shed light on the exact structure of the EH (which has been discussed only for free theories and for some gapped phases under specific conditions), and (2) they cannot be extended (in a scalable way) to calculate the full entanglement spectrum for interacting theories, and they are in fact limited to only few dozens eigenvalues for symmetry sector. Our approach, however, is based on the relation between the VNE and the thermodynamic entropy of the EH (equation (6)), which can be computed in a scalable way by quantum Monte Carlo methods in D>1, as we describe in the section 3. Below, we report several systematic checks on the validity of our approach, via exact benchmarking to one-dimensional spin systems-for additional checks, see [22], and demonstrate that it allows to make numerical predictions on entanglement properties in two-dimensions that are not accessible by any other method.

Measuring entanglement entropy in numerical simulations
Before discussing the concrete validation examples, we illustrate how to measure VNE in numerical simulations which do not have access to the system wave function. The strategy relies on any numerical method that is able to compute the thermodynamic entropy of the BW-EH at the entanglement temperature, β EH . This can be achieved using QMC algorithms based on Wang-Landau (WL) sampling [38]. Below, we illustrate this by applying the quantum version of the WL method performed in the stochastic series expansion (SSE) QMC framework [39,40]. Compared with the conventional 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' β 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 [38]. The key idea of the method is to calculate the density of states, ρ(E) by considering a non-Markovian sampling. 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 where the nth order series coefficient g(n) plays the role of the density of states in the classical algorithm. We use both local and loop updates (directed loop updates) in the WL-SSE sampling. The use of loop updates is particularly important to avoid problems with the slowing down of the configuration selection process in a inhomogeneous Hamiltonian such as the BW-EH [15] and at critical points [47]. We refer to [39,40] for the general details of the computation of g(n), and in the appendix we discuss the technical aspects of the simulation that are relevant to reproduce our results.
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)∼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)→g old (n)f. This process 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 )→ln( f old )/2 [41]. This process is repeated until convergence is achieved. Here we use the convergence condition proposed in [41,42].
In addition to the aforementioned algorithm, we consider the optimized-broad-histogram algorithm proposed in [40] for the 2D Heisenberg model, see figure 3(a). These results were obtained with the ALPS code [43] 7 . In this case, we confirm that the two methods give the same results (within error bars). One point that is worth to emphasize is that the method is straightforward to implement on a working WL code (only requires to implement an inhomogeneous version of the system Hamiltonian, as equation (3)) (see footnote 6).

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)∝cln L, where c is the central charge of the underlying CFT.
We consider the BW-EH for the one-dimensional Heisenberg model (HM) at its quantum critical point g=1. In figure 2, we plot the BW VNE under both PBCs and OBCs. 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 [44] (HM) and exact diagonalization methods (QIM) 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 [45,46] 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) [16]: the two are separated only by a constant shift that depends solely on the central charge.
For the PBC case 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 figures 2(b1) and (c1). For the OBC case, we observe an alternating term of the BW-EH entropy for the HM, but not for the QIM, see figures 2(b2) and (c2). These results is in agreement with the exact VNE. As discussed in [47,48], those oscillations are universal and due to the antiferromagnetic nature of the interactions, not appearing in the QIM [49] (in the latter, the effective Fermi momentum is either 0 or π). From the CFT perspective, the oscillations can be viewed as lattice corrections of scaling dimension Δ p : their decay as a function of the bipartition size is a power law whose exponent is related to Δ p [48,50].
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 is reported elsewhere [22]). 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 nonuniversal contributions are accurately captured: for instance, Δ S(L) goes to zero in the limit  ¥ L 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 [16]. 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.

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 [3], 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 spontaneously-symmetry-broken (SSB) phase [10,17]. 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 β EH =2π/v, with v Heis =1.658J [51] and v XY =1.134J [52], using the WL-SSE algorithm.
In order to illustrate how to cast the BW-EH on two-dimensional lattices [16], we consider the 2D Heisenberg model in a square lattice L x ×L y . In this case, the BW-EH is , , x y x y x y x y 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 figure 3(a). The function Γ(x) is given by the BW theorem (equation (4)) which represents the EH of a half-bipartition; we call this subsystem-geometry of cylinder. We also consider the CFT expression (equation (5)), which corresponds to the generalization of the BW to a subsystem that is embedded in a infinite system; we call this subsystem-geometry of toroid. We remind the reader that, as discussed in [16], for finite values of L y , formula equation (5) is in principle only applicable to conformal field theories. 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 ( ) P 1 model describing the emergent Nambu-Goldstone modes, and the gapped part of the theory is described by the massive Goldstone mode.
In figure 3(b1), 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 spinwave approximation [53] (discrepancy <3%). In figure 3(b2), we extract the subleading logarithmic correction by considering the entropy difference ( ) ( ) 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 [10,19,53,54], 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 [17,19] may signal the fact that the latter are more affected by irrelevant operators, as observed in 1D [47,48,50], or may be due to the smoother continuity properties of the VNE.  In panel (a) we present the sketch of the two-dimensional system considered in this work. The BW-EH is defined in the half-bipartition A. Panel (b) shows results for the HM and XY model. The x-axis of (b1) represents the linear size of the boundary, L y =L, and the subsystem aspect ratio for the HM (torus) is a.r.=L y /L x =1, while for the XY (torus) and the HM (cylinder), a.r.=2. In panel (b2), we remove the area law terms of S, and plot the subleading term of S as function of ln L. The number of Goldstone modes, n b =2b, extracted with a linear fit, is in agreement with expected results.

Strongly coupled quantum criticality
As a second example of 2D system, we consider the BW-EH for the bilayer Heisenberg model [55,56] x y x y x y x y where ( )  = i i i , x y labels the sites within the planes (square lattice), and l are the label of the planes. 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 β EH =2π/v, with v=1.9001(2) [51]. The function Γ(i x ) is given by the BW theorem (equation (4)), which represents the EH of a half-bipartition with a cylinder geometry; see figure 3(a). For this universal class and geometry, it has been argued that there is a universal constant correction to the entanglement entropy that depends solely on the aspect ratio [23,57]: for a cylinder geometry with PBC in the y direction, 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 figure 4. Within error bars, our results show that S (L) is independent of the aspect ratio of the subsystem, see figure 4, have no detectable logarithmic subleading term (the S(L)=AL+B ln(L)+D fitting, gives B=−0.05(8)), and the y-intercept of S(L) is D=0.010 (7). These results confirm that the scaling of the VNE is given by a pure area law behavior as predicted by the large N calculations for O(N) critical models [23,57]

Stability of BW-EH entropy
We now discuss the stability of the approach to measure the BW-EH utilizing QMC simulations. The most critical step are uncertainty due to errors in determining β EH . Since the density of states of the EH has qualitatively distinct properties from conventional density of states, it is of key importance to understand the sensitivity of the approach proposed here to such errors.
In figures 5(a1, b1), we show the value of the extracted entropy obtained via Wang-Landau sampling as a function of β, for both 1D and 2D HM. The insets magnify the region in the vicinity of the exact value of β EH , signaled by a dashed vertical line: in this regime, the entropy is linearly sensitive to β. 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 [51,58].

Stability with respect to inhomogeneous couplings
The second class of imperfections we discuss is the presence of inhomogeneities in the BW-EH couplings. This is motivated by potential experimental realizations of the EH: indeed, the generic approach described above can be extended to formulate protocols to measure the von Neumann entropy in experiments (complementing previous approaches based on Renyi entropies [59][60][61][62][63][64][65] and entanglement spectra [15,66]) as follows.
The key ingredient here is to obtain the density of states of the EH, whose microscopic implementation has been discussed in [15]. Such density of states can be obtained via quantum quenches, adapting to experiments the approach presented in [67,68]: the main idea is that, starting initially from the ground state, one can recover the eigenvalue density of a given operator (in our case, the BW-EH) by suddenly switching on a coherent dynamics given by the operator itself. This technique is a quantum quench analog of spectral decomposition, and is thus fully general and applicable to the EH.
The main challenges to be overcome in this direction are three: errors in the initial state preparation, finite quench time due to decoherence, and proper realization of the microscopic dynamics. Regarding the first two elements, the analysis of the EH case goes along the same lines of conventional Hamiltonians [67,68]. We thus focus here on the last element, which is unique to the present case due to the spatially modulated couplings. We thus address the effects of random perturbations in the EH couplings Γ(n), which accounts for possible imperfect experimental realizations of the EH. Such random perturbations are one of the possible cases discussed in  [15].
We consider disordered couplings, Γ(n)→Γ(n)(1+δ n ) 8 , where δ n =[−δ, δ], in the BW-EH of the HM (in 1D and 2D). Specifically, we are interested in understanding how the BW VNE is affected by a small amount of disorder.
In figure 5(a2, b2) we show that the mean value of the BW VNE is not appreciably affected by disorder up to strength of the order of 10%. For larger values of δ, we observe a considerable dependence on the disorder realization, as signaled by the visually large spreading of the values of S. Surprisingly, the mean value of the entropy is not dramatically affected. 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 simulations, 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 [36,69].

Acknowledgments
We acknowledge useful discussions with C Groß, N Laflorencie, A Laio, M Lukin, H Pichler, S Sachdev, E Tonni, B Vermersch, and P Zoller, and thank P Calabrese for feedback on the manuscript. This work has been carried out within the activities of the Trieste Institute for Theoretical Quantum Technologies (TQT). This work is partly supported by the ERC under grant number 758329 (AGEnTh), and has received funding from the European Union's Horizon 2020 research and innovation programme under grant agreement No 817482. TMS and MD acknowledge computing resources at Cineca Supercomputing Centre through the Italian SuperComputing Resource Allocation via the ISCRA grants QMCofEH. We used the ALPS to benchmark the results of our code [43]. TMS and GG equally contributed to this work.

Appendix. Quantum Wang-Landau sampling of the entanglement Hamiltonian
In this section we discuss some relevant technical aspects of the quantum Wang-Landau (WL) simulations used to obtain the thermodynamic entropy of the BW-EH. The orders of the series expansion shown in equation (7) that are relevant for a given β are sharply peaked around and the constant C is defined to guarantee that the QMC weights are always positive ( > á ñ C H BW ) [45]. Hence, one can truncate the expansion at an order Λ (i.e. n=0,1, ..., Λ ) without introducing systematic errors in the simulation. The choice of Λ is performed using the same algorithm of the conventional SSE simulations, see [45,46], which gives as a result ( ) The effect of introducing the cutoff Λ(β) is that the range of temperature that can be accessed is restricted to In order to obtain the results of figures 4(a, b1) of the main text, for instance, we simulate the BW-EH using Λ(3β EH ) as the cutoff. Instead, the computation of the BW VNE at β EH are obtained utilizing a cutoff Λ(αβ EH ). The results shown in figures 2 and 3 of the main text are obtained with α=1.3. We check that these results do not change upon increasing α.
The required numerical resources to compute the thermodynamic entropy of the BW-EH depend on the cutoff, Λ, introduced in the SSE series expansion. We can estimate the system-size dependence of the required numerical resources, by noting that the leading term of the size scaling of the SSE cutoff are Λ(β EH )∼L 3 in 2D and Λ(β EH )∼L 2 in 1D; where L is the linear size of the systems considered here. In order to establish these results we first note that the correlation functions that appear in the expression of á ñ H BW (see equation (10) In addition we note that β EH ∼O(1), and using the aforementioned argument one have that C∼L 3 . As the prefactor associated to C is larger then the one of á ñ H BW , we conclude that Λ(β EH )∼L 3 . Remarkably, the cutoff that one needs to introduce to compute the ground state properties of the original Hamiltonian H (the one from which the BW-EH ansatz is built) also scales as Λ∼L 3 . In this case, however, the ground state energy scale as ( ) b  ¥Ẽ L 2 , while one needs to consider β∼L in order to access ground state properties.
Finally, it is important to mention that the results of the BW entropy are obtained by doing an average of N r independent WL simulations, i.e. 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.