0 1 2 3 4 5 6 7 0.4 0.45 0.5 0.55 0.6 0.65 0.7 Swap r$_{max}$=4 r$_{max}$=2 r$_{max}$=1 0 4 8 12 16 0.4 0.5 0.6 0.7 Swap r$_{max}$=8 r$_{max}$=4 r$_{max}$=2 r$_{max}$=1 ℓ L = 16 L = 8 S 2 /ℓ S 2 / ℓ 0 20 40 60 80 100 0.4 0.8 1.2 1.6 DMRG Swap  j$_{max}$=5 0.0001 0.001 0.1 0.15 0.2 0.25 0.3 0.35 0.4 i 1 / m S 2 3 FI G. 3: (color online) The Renyi entropy divided by the linear di m ension ℓ of the entangled region A , for 2D lattices with L = 8 and L = 16. The data labeled “Swap” is from a single simulation calculating Eq. (9) directly. The other data sets are derived from the improved ratio estimator, Eq. (10), with different ranges of r ∈ [1 , r max ] (see text). The inset is a periodic L = 8 lattice with region A consisting of the 16 central (dark) sites labeled by ℓ = 4. FIG. 2: (color online) The Renyi entropy S 2 as a function of site index i ∈ A , for a 100-site Heisenberg chain with open boundaries, calculated with DMRG and QMC. Data labeled “Swap” was calculated with Eq. (9) with one QMC simula- tion, while data labeled j max = 5 was calculated with Eq. (10) using 20 separate QMC simulations with a range of j ∈ [1 , 5]. The inset shows the convergence of S 2 to the exact value (dashed line) for i = 6 with up to m = 4000. j$_{max}$ = 5 are sufficient to converge 100-site chain to the exact results along most of its length. To combat this issue, we propose a refinement to the algorithm, which we called “improved ratio” sampling. Consider first a one-dimensional chain, and define regions A $^{1}$, A$^{2}$, ..., A $^{n}$, such that A i contains i sites, and A i is a subset of A i $^{+1}$. In other words, each region A i +1 is obtained by adding one site to region A $^{i}$, and region A 0 is the empty set (thus, Swap$_{A}$ 0 is equal to the identity operator). In a given simulation, one can calculate the ratio 2D Heisenberg Results. – We now extend these concepts to QMC simulations of the spin 1/2 Heisenberg model on N = L × L lattices. For the simulations using the im- proved ratio estimator, 〈 Swap$_{A}$ ℓ + r 〉 / 〈 Swap$_{A}$ ℓ 〉 , we define A ℓ as a square region of linear size ℓ containing ℓ × ℓ sites, so A ℓ + r contains 2 ℓr + r 2 more sites than A $^{ℓ}$. Figure 3 illustrates S$_{2}$ on two 2D Heisenberg mod- els with L = 8 and L = 16. In the case of L = 8, data for a single simulation calculating 〈 Swap$_{A}$ 〉 directly is essentially identical to the improved ratio estimators for all ℓ . Here, the direct estimator (Eq. (9)) corre- sponds to one simulation, while r$_{max}$ = 1 corresponds to 6 different simulations of the improved ratio estima- tor 〈 Swap$_{A}$ ℓ +1 〉 / 〈 Swap$_{A}$ ℓ 〉 and ℓ ∈ [1 , 6]. In contrast to L = 8, for L = 16 the direct estimator is significantly different than the improved ratio estimator. We can see the convergence of the data as one successively decreases the range of r ∈ [1 , r$_{max}$ ], so that r$_{max}$ = 2 is identical to r$_{max}$ = 1, which is the smallest possible range for the improved ratio estimator in the 2D geometry illustrated. It is important to note that for large r$_{max}$ the correct value of S$_{2}$ does not lie within the error bars of the data. This is an indication that simulation ergodicity may be jeopardized by low sampling of the Swap$_{A}$ operator. We use the improved ratio estimator with r$_{max}$ = 1 to scale the results for the Renyi entropy of the Heisenberg model to larger system sizes in 2D. Results are given in Fig. 4 for L = 4 to L = 32, plotted as S$_{2}$/ℓ . From i.e., a unique QMC simulation must be done for each de- sired i . If every site were used for a different simulation weight, the QMC algorithm would obtain an additional multiplicative factor of N in its scaling. In practice, this can be reduced by noting that good statistical control can be retained by calculating 〈 Swap$_{A}$ i + j 〉 / 〈 Swap$_{A}$ i 〉 for fixed i and a range of j ∈ [1 , j$_{max}$ ] (see also the Discus- sion). We illustrate this in Fig. 2, where simulations with where, for each i = 0 , ..., n − 1, the log of this ratio is equal to minus the difference S$_{2}$ ( ρ$_{A}$ i +1 ) − S$_{2}$ ( ρ$_{A}$ i ). Computing this ratio for each i will let us compute S$_{2}$ for all A $^{i}$. To calculate this ratio, the simulation must be performed with the modified sampling weight, W$_{l}$W$_{r}$ 〈 V ( l ) | Swap$_{A}$ i | V ( r ) 〉 , (11) 〈 Swap$_{A}$ i +1 〉 〈 Swap$_{A}$ i 〉 = ∑ $_{rl}$W$_{l}$ W$_{r}$ 〈 V ( l ) | Swap$_{A}$ i +1 | V ( r ) 〉 $_{rl}$W$_{l}$ W$_{r}$ 〈 V ( l ) | Swap$_{A}$ i | V ( r ) 〉 (10)