the notation here, referring the reader to Refs. [10–12] for details of implementation. The method is a T = 0 projector QMC, where high powers of a Hamiltonian H are used to project out the ground state | Ψ$_{0}$ 〉 from a trial wavefunction Ψ: ( − H ) $^{m}$| Ψ 〉 ∝ | Ψ$_{0}$ 〉 . In Sandvik’s QMC scheme, ( − H ) m is written as a sum of all products of m bond operators, basis. One begins by writing the expansion of a singlet wavefunction Ψ$_{0}$ of N spins as N/ 2 valence bonds: where ( ρ$_{A}$ )$_{α}$$_{1}$$_{,α}$$_{2}$ = ∑ $_{β$_{1}$}$C$_{α}$$_{1}$$_{,β}$$_{1}$ C$_{α}$$_{2}$$_{,β}$$_{1}$ denotes a matrix element of ρ$_{A}$ . Therefore, 〈 Swap$_{A}$ 〉 = 〈 〈 V ( l ) | Swap$_{A}$ | V ( r ) 〉 〈 V ( l ) | V ( r ) 〉 〉 , (9) 〈 Ψ$_{0}$ ⊗ Ψ$_{0}$ | Swap$_{A}$ | Ψ$_{0}$ ⊗ Ψ$_{0}$ 〉 = ∑ α$_{1}$,α$_{2}$,β$_{1}$,β$_{2}$ C$_{α}$$_{1}$$_{,β}$$_{1}$C$_{α}$$_{2}$$_{,β}$$_{1}$C$_{α}$$_{2}$$_{,β}$$_{2}$C$_{α}$$_{1}$$_{,β}$$_{2}$ = ∑ α$_{1}$,α$_{2}$ ( ρ$_{A}$ )$_{α}$$_{1}$$_{,α}$$_{2}$( ρ$_{A}$ )$_{α}$$_{2}$$_{,α}$$_{1}$ = Tr( ρ 2 $_{A}$) , (4) S$_{2}$ ( ρ$_{A}$ ) = − ln(Tr( ρ 2 $_{A}$)) = − ln( 〈 Swap$_{A}$ 〉 ) . (5) 〈O〉 = ∑ rl 〈 V | P ∗ l O P$_{r}$ | V 〉 rl 〈 V | P ∗ l P$_{r}$ | V 〉 = ∑ $_{rl}$W$_{l}$ W$_{r}$ 〈 V ( l ) |O| V ( r ) 〉 $_{rl}$W$_{l}$ W$_{r}$ 〈 V ( l ) | V ( r ) 〉 . (8) Swap$_{A}$ ( ∑ α$_{1}$,β$_{1}$ C$_{α}$$_{1}$$_{,β}$$_{1}$ | α$_{1}$ 〉| β$_{1}$ 〉 ) ⊗ ( ∑ α$_{2}$,β$_{2}$ D$_{α}$$_{2}$$_{,β}$$_{2}$ | α$_{2}$ 〉| β$_{2}$ 〉 ) = ∑ α$_{1}$,β$_{1}$ C$_{α}$$_{1}$$_{,β}$$_{1}$ ∑ α$_{2}$,β$_{2}$ D$_{α}$$_{2}$$_{,β}$$_{2}$ ( | α$_{2}$ 〉| β$_{1}$ 〉 ) ⊗ ( | α$_{1}$ 〉| β$_{2}$ 〉 ) . (3) P$_{r}$ | V 〉 = W$_{r}$ | V ( r ) 〉 . (7) ( − H ) m = ∑ r m ∏ i =1 H$_{a}$ r $_{i}$b r i = ∑ r P$_{r}$, (6) | Ψ$_{0}$ 〉 = ∑ r f$_{r}$ | ( a r $_{1}$, b r $_{1}$) . . . ( a r N/ $_{2}$, b r N/ $_{2}$) 〉 = ∑ r f$_{r}$ | V$_{r}$ 〉 (2) To derive an estimator for S$_{2}$ , we construct two copies of the system (higher entropies S$_{n}$ require n copies); one “real”, one “ancillary” (Fig. 1). Then we define an oper- ator Swap$_{A}$ which acts between the copies of the system, swapping the configurations within region A . To define this operator, we give its matrix elements in a product basis, e.g. the S z basis. Let | α 〉 be a complete basis of states in the region A and let | β 〉 be a complete ba- sis of states in the complement region B . The state of each copy can be decomposed in this product basis as | Ψ 〉 = ∑ $_{α,β}$C$_{α,β}$ | α 〉| β 〉 for some amplitudes C$_{α,β}$ . We define the unitary operator Swap$_{A}$ by where for concreteness we use the spin 1/2 Heisenberg model, defining H = − ∑ 〈 a,b $_{〉}$H$_{ab}$ and H$_{ab}$ = − ( S a · S b − 1 / 4). The “operator string” P$_{r}$ is sampled according to its weight, W$_{r}$ , accrued upon evolution of a trial valence bond state | V 〉 under projection: As shown in Ref. [10], for the Heisenberg model, W$_{r}$ is simply related to the number of off-diagonal operations m$_{off}$ in the projection, W$_{r}$ = 2 − $^{m$_{off}$}$. To sample the Swap$_{A}$ operator, one requires a double- projector valence bond QMC scheme [10], where a gen- eral expectation value for any observable O is given by In this case, the two operator strings, P$_{l}$ and P$_{r}$ , are ap- plied to two copies of the system; the expectation value of the Swap$_{A}$ operator as illustrated in Fig. 1 can then be calculated directly. Specifically, one performs im- portance sampling of operator strings according to the weight W$_{l}$W$_{r}$ 〈 V ( l ) | V ( r ) 〉 , and measures the QMC aver- age expectation value The two copies of the system will not interact with each other, so the ground state of the combined system will be the product of the ground states on each copy: | Ψ$_{0}$ ⊗ Ψ$_{0}$ 〉 . Thus, the expectation value of Swap$_{A}$ is calculating S$_{2}$ from Eq. (5). Results for a 1D chain are illustrated in Fig. 2. There, simulations were performed using a double-projector QMC, with a simple “columnar” trial state | V 〉 (alternat- ing nearest-neighbor bonds in 1D). In the following data, four random bond operators were changed in each oper- ator string per Monte Carlo step, and (unless otherwise stated) m = 40 N . Exact results for S$_{2}$ were obtained with a DMRG simulation, converged with 1000 states. One immediately sees that the naive expectation value 〈 Swap$_{A}$ 〉 results in very large statistical error bars when the region A grows large (in 1D, A is the linear region of size i ). One may understand this by considering the ex- pected scaling of S$_{2}$ . Namely, S$_{2}$ scales logarithmically with the size of A , so the expectation value 〈 Swap$_{A}$ 〉 should be polynomially small in A . As a result, there are only very few configurations contributing significantly to Equation (5) is basis independent [13]. In particular it holds in the valence bond basis, where the Swap$_{A}$ oper- ator is not defined by Eq. (3). Rather, it is defined to swap the endpoints of valence bonds between the real and ancillary copies in the region A , as illustrated in Fig. 1. Measuring the Swap operator in QMC. – We now present a procedure for measuring the entropy S$_{2}$ via importance sampling of the Swap$_{A}$ operator in QMC. A projector scheme in the valence bond basis has recently been pioneered by Sandvik, and we only briefly review 2