Improving entanglement and thermodynamic Rényi entropy measurements in quantum Monte Carlo
Abstract
We present a method for improving measurements of the entanglement Rényi entropies in quantum Monte Carlo simulations by relating them with measurements of participation Rényi entropies. Exploiting the capability of building improved estimators for the latter allows to obtain very good estimates for entanglement Rényi entropies. When considering a full system instead of a bipartition, the method can be further ameliorated providing access to the thermodynamic Rényi entropies with high accuracy. We also explore a recentlyproposed method for the reconstruction of the entanglement spectrum from entanglement Rényi entropies and finally show how potential entanglement Hamiltonians may be tested for their validity using a comparison with thermal Rényi entropies.
pacs:
02.70.Ss,03.67.Mn,75.10.Jm,05.10.LnI Introduction
Quantum entanglement has been of interest since the early days of quantum mechanics Einstein et al. (1935). The quantification of the entanglement in interacting many body quantum systems has attracted a lot of attention during the last decade for several fundamental and practical reasons Calabrese et al. (2009). Entanglement properties of onedimensional quantum problems can be treated fully analytically only in a limited number of cases (see e.g. Refs. Jin and Korepin, 2004; Fan et al., 2004), or asymptotically for conformally invariant Calabrese and Cardy (2004) or some disordered systems Refael and Moore (2004). They remain generically accessible to numerical Density Matrix Renormalization Group (DMRG) calculations, provided the entanglement between the subsystems is not too large Schollwöck (2005). For higher dimensional systems however, exact methods are much more difficult to implement. Nevertheless, remarkable progresses have been made recently, e.g. using series expansions Kallin et al. (2011); Singh et al. (2012), numerical linked cluster expansion Kallin et al. (2013, 2014), or using quantum Monte Carlo (QMC) simulations Alet et al. (2007); Kallin et al. (2009); Hastings et al. (2010); Melko et al. (2010); Humeniuk and Roscilde (2012); Grover (2013); Assaad et al. (2014); Chung et al. (2013); Helmes and Wessel (2014); Herdman et al. (2014), which is precisely the topic of the present work.
In nonfrustrated quantum spin systems, standard thermodynamic observables can be obtained to very high accuracy within QMC simulations Sandvik (2010). Here we are interested in the Rényi entanglement entropy (EE)
(1) 
where is the reduced density matrix, assuming that is a subsystem imbedded in a larger system. Clearly cannot be related to a simple thermodynamic observable, e.g. a correlation function (except for noninteracting systems Chung and Peschel (2001); Audenaert et al. (2002)). At zero temperature, Hastings et al. Hastings et al. (2010) developed a technique based on the introduction of a “swap”operator in a projector Monte Carlo approach to tackle this issue. At finite temperature, several techniques have been explored, including temperature integration Melko et al. (2010) and Wang Landau sampling Inglis and Melko (2013).
Perhaps the most elegant method was brought forward by Humeniuk and Roscilde Humeniuk and Roscilde (2012). Their method for the calculation of entanglement Rényi entropies of order for a subsystem in path integral QMC methods is based on the observation Calabrese and Cardy (2004) that they are related to the ratio of partition functions . Here, is the partition function of replicas glued together at one imaginary time slice on the subsystem only. is the partition function of independent replicas. Here and from now on, is an integer .
In a simulation which samples both partition functions in a generalized ensemble, proposing moves between the two ensembles, the estimator for the entanglement Rényi entropy is given byHumeniuk and Roscilde (2012):
(2) 
where is the number of QMC configurations observed in the glued ensemble, while is the number of QMC configurations seen in the independent ensemble. This method becomes inefficient for too large entropies, which lead to very small and large . This is a problem of rare events which is also known in the related context of participation Rényi (PR) entropiesLuitz et al. (2014a, b, c) and prohibits the estimation of large entropies in finite simulation time. The problem can however be cured by the application of the “ratio trick” de Forcrand et al. (2001); Hastings et al. (2010), calculating the entanglement entropy by a stepwise increase of the subsystem .
Let us give a description of the results presented in this article, along with its organization. We will first start from equation (2) to show how the measurement of the entanglement entropy is related to the basis dependent participation Rényi entropies Stéphan et al. (2009); Luitz et al. (2014c) (also called ShannonRényi entropies in the litterature). The basic idea is to split the extended ensemble in its two parts and to simulate the ensemble of independent replicas and the ensemble of replicas that are glued together on subsystem (see Fig. 1) separately. In Sec. II, we derive the following relation:
(3) 
which relates the entanglement entropy to the difference between participation Rényi entropy and the replica correlation , which is introduced in Sec. II and is defined in the glued ensemble.
Remarkably, this combination of two basisdependent quantities will hint towards a more efficient calculation of the entanglement entropy, for two different reasons. First, large entanglement entropies can be obtained (both in our setup and the one used in Ref. Humeniuk and Roscilde, 2012) when using a QMC computational basis where the participation entropies are small. Second, we will introduce in Sec. III several improved Monte Carlo estimators which will greatly increase the precision on and consequently on .
As an interesting byproduct, our scheme allows to compute the thermodynamic Rényi entropy in the specific case where the full system and the subsystem are identical, as developed in Sec. IV. We will show there how the translation invariance in imaginary time can also be used to construct an improved estimator for the replica correlation , leading to a very accurate result for the thermodynamic Rényi entropy. Let us emphasize that the thermodynamic Rényi entropy (for integer ) can be calculated from two standard QMC simulations (using independent replicas) without the need of implementing a different stochastic process or WangLandau sampling.
The new methods are extensively tested and their efficiency discussed in Sec. V where we provide several results on quantum spin chains and ladders. We also explore in this section the possibility to reconstruct the entanglement spectrum from Rényi entanglement entropies Song et al. (2012); Chung et al. (2013), discussing the limitations of this approach. There, we also propose an alternative method to test a putative entanglement Hamiltonian, based on the comparison between Rényi entanglement entropies and the Rényi thermodynamic entropies of the entanglement Hamiltonian.
Finally, Sec. VI draws conclusions on our work while the appendices contain details on the improved estimator derivation as well as on using symmetry sectors when measuring entanglement entropies.
Ii Method
The method proposed by Humeniuk and Roscilde Humeniuk and Roscilde (2012) uses an extended ensemble simulation (see Fig. 1), which dynamically moves between the glued () ensemble and the independent ensemble () and records the ratio of Monte Carlo steps performed in the glued vs. independent ensembles. For the equilibrium (Monte Carlo) time evolution of the probability to be in ensemble , the Master equation
(4) 
holds, where is the probability of moving from the independent ensemble to the glued ensemble and is the probability of the inverse move. In an equilibrated Markov chain of QMC configurations, the probability of finding the glued ensemble is timeindependent, and we therefore obtain:
(5) 
and with equation (2)
(6) 
Instead of calculating this ratio of probabilities in a single QMC calculation, let us now concentrate on an estimation of and in separate calculations.
ii.1 Probability of leaving the independent ensemble
If the simulation is in the independent ensemble, the condition for moving to the glued ensemble is given by finding identical states on the subsystem in all replicas. This corresponds to identical states for all replicas following the convention of Fig. 1. It is important to note that in principle due to the cyclicity of the trace the time slice where the replicas would be sewed together does not matter. However, we will not actually perform the step of moving to the glued ensemble here, but will just think about how probable it is. Clearly, we then have:
(7) 
But the probability of finding identical states on the subsystem in all replicas is just given by the participation Rényi entropiesLuitz et al. (2014a) of subsystem
(8) 
We will present in Sec. III improved estimators for estimating efficiently the participation Rényi entropies.
ii.2 Probability of leaving the glued ensemble
When the simulation explores the glued ensemble , the condition of moving to the independent ensemble is given by having identical states on top and bottom of each replica individually, such as to meet the condition of the trace. Therefore, can be estimated by performing a simulation in the glued ensemble and recording how often this condition is met, relative to the total number of QMC steps. We have
(9) 
This step turns out to be the bottleneck of the method in Ref. Humeniuk and Roscilde, 2012, as this probability decays exponentially with the number of degrees of freedom in the subsystem and thus the replica correlation
(10) 
exhibits a volume law. Note that the participation Rényi part suffers similar exponentially small probabilities for which we can however improve the estimate (see Sec. III). This volume law is directly related to the problem of low acceptance rates in the standard method Humeniuk and Roscilde (2012). It should be noted that is not an entropy in the sense of equation (2) as in general no density matrix can be found that provides for all . In particular, can grow with , which is not possible for Rényi entropies.
ii.3 Entanglement entropy
Combining the two results, we therefore arrive at the previously announced Eq. (3) for measuring entanglement entropies, which is restated for clarity using the spatial region index:
(11) 
This relation between the Rényi form of the participation and entanglement entropies provides very interesting insights in the performance of any method based on the idea by Humeniuk and Roscilde Humeniuk and Roscilde (2012). It has been established that PR entropies show a volume law in local bases, Stéphan et al. (2009, 2010, 2011); Luitz et al. (2014a, b, c). However the coefficient of the volume term is nonuniversal and depends on the basis. As entanglement entropies usually display an area law for condensedmatter groundstates Eisert et al. (2010), the replica correlation necessarily has to exhibit the same volume law (and correspondingly, the probability to leave the glued ensemble decreases exponentially with the number of degrees of freedom in ).
The same reasoning goes for the basis dependence: depends on the basis, while does not. Thus, needs to be a basis dependent quantity, too. As the probabilities to be observed in the QMC calculations are potentially very small, it is beneficial to choose the basis in which they assume larger values. Therefore, one should try to choose a computational basis in which the participation entropy is the smallest. One notable example is a simulation of the XX model in the basis in which is diagonal instead of the usual basis.
Iii Improved estimators
Too large entropies lead in general to statistical issues in the QMC simulations. In order to tackle problems connected to the corresponding rare events, it is useful to increase the number of Monte Carlo measurements as much as possible. Here, we present improved estimators that greatly enhance the precision of participation entropies using all possible symmetries in imaginary time and real space.
The starting point is the replica method introduced in Ref. Luitz et al., 2014a. The basic idea is that in order to measure the PR entropy of a subsystem (which may coincide with the full system), it is sufficient to estimate the probability of finding the same state on the part corresponding to subsystem of the state at operator string slice . According to the convention given in Fig. 1, this corresponds to the estimator
(12) 
where the first sum runs over the Markov chain of length and the second sum over all slices of the operator strings in the replicas. For simplicity, we enforce the same cutoff for all replicas.
This method can be greatly enhanced by the observation that all operator strings in the ensemble are independent. This leads to two possible improvements:

Due to the cyclicity of the trace, each operator string has a cylinder topology and can thus be independently translated cyclically by any number of states in the “imaginary time” direction. Each of the such transformed Monte Carlo configurations has exactly the same weight.

If the system is invariant under (spatial) symmetry transformations, we can transform the whole operator string of each replica with independent transformations without changing the weight of the configuration.
These two recipes can be used to greatly improve the quality of the estimation of . However, a naive application of these ideas is not possible as it is far too expensive to try all combinations of shifted and transformed operator strings. As discussed below, it is possible to exploit all symmetries by only one pass on each operator string.
iii.1 PR entropy of the full system
Let us first start with the full system as the symmetries are simpler to apply and the resulting improved estimator formulae are clearer.
For each state in the operator string of replica , we calculate the parent state by applying all model symmetries. The application of all symmetries classes all the basis states into nonoverlapping families of states, each of which is represented by the uniquely defined parent state – here, we will take the state in the family with the smallest binary representation. It is important to record the multiplicity (i.e. the number of states belonging to the family) of each state family for the purpose of correct normalization. Note that most state families have the maximal multiplicity given by the number of symmetries , with the exception of high symmetry states for which .
While transversing all the operator strings, we record the histogram of the number of occurrences of states with parent in operator string .
In the next step, for each parent state that has been observed in one of the replicas we have to count the number of configurations (of symmetrytransformed operator strings) in which we can find identical states in all replicas. It is clearly given by
(13) 
Note that the multiplicity of the parent state accounts for the fact that we can have any of the states in the family of parent as the identical state in all replicas.
Finally, we have to normalize equation (13) by the total number of symmetry equivalent configurations of all replicas. As the family of parent consists of states, the correct normalization is . Together with the possible cyclical shifts of the operator strings, this yields the improved estimator in Monte Carlo configuration for the probability of moving from the independent ensemble to the glued ensemble
(14) 
Here, is the number of states in the operator string (forced to be equal in all replicas for simplicity). The sum runs over all parent states recorded in the histogram . In equation (14) it is immediately clear why this method is extremely beneficial for the observation of small probabilities: the normalization factor can become a very small number, as typically the number of applied symmetries for a translationally invariant system with sites and in our calculations. For large values of , very small numbers can be obtained by one Monte Carlo measurement and the variance of the estimator is greatly reduced.
Because of the tremendous number of possible combinations of symmetry transformed replica operator strings, it should be noted that the evaluation of Eq. (14) may cause numerical problems. The products of the numbers for all replicas can easily become too large to be stored as 64bit integers. One solution for this problem is to use extended precision floating point numbers. We have found, however, that it is sufficient to perform the products using double precision floating point numbers and performing the sum using Kahan’s summation algorithmKahan (1965) to avoid precision loss and cancellation effects.
An interesting further improvement of the method stems from the fact that a single simulation of replicas can be used for the calculations of PR entropies (or equivalently the probabilities ) with ranging from to . Indeed, at the time of performing a Monte Carlo measurement, we can select any combination (without repetition) of replicas out of the copies and apply Eq. (14) without the need of creating a new histogram. This can be repeated for all possibilities, further improving the precision of the estimate:
(15) 
where the additional sum runs over all combinations of the replicas.
iii.2 PR entropy of subsystem
If the subsystem and the full system are not identical, we have to slightly modify the procedure described above. The reason for this is the fact that if we cut out the part of subsystem from every state in the family corresponding to a parent and for a different family corresponding to , we will find that the families can now have an overlap if the full system has more or different symmetries than the subsystem, which is generally the case. In some cases, where the exploited symmetries of the subsystem are identical with the symmetries of the full system, the same algorithm as for the full system may be used. This is generally the case if any symmetry transformation maps the subsystem on itself, e.g. in the case of the periodic ladder (see below, Sec. V.2) when solely using translation symmetries along the ladder.
In the general case, accounting for these overlaps is expensive as all pairs of parent states in the histogram have to be treated. We therefore choose to go through all parents in the histogram and create the family of states from which we deduce the subsystem states by cutting out the corresponding part. In doing so, we generate a new histogram filled with the cut states where we accumulate the corresponding from histogram . Note that the histogram may become very large if the number of symmetries and the number of lattice sites in subsystem is large. The size of the histogram is always smaller or at most equal to (typically, its size is reduced by the number of symmetries to ). The maximal size of the histogram is, however, given by , where is the dimension of the Hilbert space of states on subsystem . Note that the relevant number of symmetries here is the number of applied symmetries of the full system as we discuss here the general case, in which the subsystem has no (or less) symmetries.
The equation for the estimator in the previous section only have to be slightly modified:
(16) 
where is then the number of times the subsystem state has been observed in all symmetry realizations of replica . This yields the final estimator
(17) 
iii.3 Autocorrelation problem for large values of
While the improved estimator performs remarkably well for a wide range of and compares perfectly with exact results for small systems (see Sec. V), we have found that it may yield wrong results for large values of in some extreme cases. A detailed investigation shows that this behavior stems from an increasing variance of the improved estimator with together with an increasing autocorrelation time which may exceed the simulation time and thus yield systematic errors^{1}^{1}1It should be noted that this problem only occurs in the regime of very large PR entropies, far beyond the reach of the simple estimator given in Eq. (12)..
The reason for this behavior is identified by studying the time series of the estimator for different values of (see Fig.2). Clearly, for larger values of , the time series shows more and more pronounced “spikes” that at the same time become rare and large for large . We have found that the most severe spikes are created by the occurrence of the most probable states (for antiferromagnetic systems, these are the Néel states and , see Refs. Luitz et al., 2014a, b, c). These states occur with the highest probability in any operator string and have an enhanced symmetry (for instance the symmetry family of the Néel states only has a multiplicity). Consequently, the factor in Eq. (15) is much larger than for less symmetric states (for which it is usually ), creating extremely large spikes for large if the Néel states are observed in all replicas simultaneously. The frequency of the spikes becomes however rare with growing , as the simultaneous observation probability of the Néel states in all replicas decreases rapidly as .
However, we can remedy this situation by calculating the probability separately, recording the frequency of the Néel states across all replicas. We then deliberately exclude the most probable states with common parent from the estimator in Eq. (15) (similarly for Eq. (17)) and obtain
(18) 
In order to pull out the term out of the sum, we have used the fact that on average the frequency of the Néel states in replica of length is given by .
In Fig. 2, we compare the time series of (minus log of) the estimator given by Eq. (15) and the improved estimator (excluding the additive term ) given by Eq. (18) in a test example discussed deeper later. Clearly, for the two estimators behave nearly exactly equally. However, as grows, the estimator not excluding the most probable states (shown in red) has a much larger variance than the estimator that excludes the most probable state. In addition to that, the spikes created by the Néel state become rare with growing and eventually may not be even recorded once in a simulation, leading to incorrect results (especially since will dominate the PR entropy for large values of ).
Note that usually can be calculated with a much higher precision than any other state, leaving us with the possibility to correct efficiently this systematic error. Any other state can be taken out of the estimator as in Eq. (18) if its probability can be measured to a sufficient accuracy directly. This way, a hybrid method capturing the histogram of the most probable states with high precision and calculating the correction caused by the less probable states by the replica trick can be easily constructed.
iii.4 Measurement of the replica correlation entropy
The final element to compute the entanglement entropy is the measurement of the replica correlation entropy , for which the probability has to be estimated efficiently.
In the method of Ref. Humeniuk and Roscilde, 2012, the condition for moving from the glued ensemble to the independent ensemble is given by . In this case, the glue can be cut and rewired thus moving from the glued to the independent ensemble without changing the weight as illustrated in Fig. 1.
In our scheme, we can simply simulate the glued ensemble and measure the number of times we observe relative to the total number of QMC steps. Unfortunately, the glued replicas can not be changed by separate symmetry transformations and the topology of the glued operator string is rigid in imaginary time, i.e. it can not be translated. Generally, it is therefore not possible to construct an improved estimators for in the same way as for .
Iv Thermodynamic Rényi entropies
Let us discuss in more detail the special case in which the subsystem and the full system are identical. In that case, the entanglement Rényi entropies reduce to the thermodynamic Rényi entropies at inverse temperature , which are defined by:
(19)  
Here, denotes the free energy of the system governed by the Hamiltonian at inverse temperature . Noting that the thermal density matrix is nothing else but the reduced density matrix of the subsystem , we can use Eq. (3) to calculate the thermodynamic Rényi entropy.
As in the case of the entanglement Rényi entropy, we will decompose the thermodynamic Rényi entropy in a difference of the participation PR entropy and the replica correlation. The discussion of Sec. III on how to build improved estimators exploiting spatial and imaginary time symmetries carries on for the calculation of the PR entropy.
The calculation of the replica correlation for can be much further improved than in the generic case. Indeed this is the case where the replicas are glued on the full system and the periodicity in (or in the cutoff for SSE) is replaced by a periodicity for a larger operator string. Clearly, this situation can be achieved by simply performing a standard SSE simulation with a single replica at inverse temperature . The estimator for the replica correlation is then given by the probability of cutting the enlarged configurations in valid SSE replicas at inverse temperatures , with the condition that these replicas have to be periodic in .
It should be observed that there are multiple valid ways of slicing the large configuration in parts and in fact any partition is valid as long as every cutoff is large enough such as to correctly sample the inverse temperature . The probability of slicing the large replica in valid parts can be estimated (see Appendix A for a detailed derivation) by simply measuring the observable
(20) 
for any Monte Carlo configuration in the replica simulated at inverse temperature . Here, the Kronecker delta yields if the states at the beginning of each replica of length are identical, which is precisely the condition for obtaining valid, periodic replicas after performing the cut. Note that one should of course average over several possible partitions with and all translations of the partitions in imaginary time.
We obtain
(21) 
where the average is performed over the Markov chain and different partitions in order to improve the statistics. The thermodynamic Rényi entropy is thus finally given by
(22) 
where the absence of the subsystem index indicates that the full system is to be considered. We emphasize that both and are obtained within a standard SSE independent ensemble simulation.
V Results
In this section we present various results obtained on simple model Hamiltonians, such as Heisenberg chains and ladders, in order to carefully test the method. We compare, when possible, our QMC estimates with exact diagonalization (ED) or DMRG results. The quantities of interest we discuss in the rest are the Rényi entanglement (zero temperature) and thermodynamic (finite temperature) entropies. We also compare the efficiency of our new method for calculating entanglement entropies with the method of Ref. Humeniuk and Roscilde, 2012.
The last part of this section deals with a careful analysis of the reconstruction of the entanglement spectrum given entanglement entropies. We finally provide a quantitative analysis of how to probe an ansatz entanglement Hamiltonian in the case of Heisenberg ladders.
v.1 Heisenberg chain
As a first application of the different aspects of the method introduced in Sec. II we perform calculations for the wellstudied antiferromagnetic Heisenberg chain of spins , described by the Hamiltonian
(23) 
using periodic boundary conditions .
v.1.1 Entanglement entropies
We test our implementation on the example of a chain of spins, which can be easily solved using ED. As a way to show the different elements in our method, we first display in the main panel of Fig. 3 the participation entropies as a function of the subsystem size and Rényi index . The correspondence between the QMC and exact result is perfect and the precision is such that the error bars are not even visible in the graph. The inset displays the probabilities of finding the same basis state in replicas for ranging from 2 to 10, for the full system . The exponential decay of is perfectly reproduced by the Monte Carlo result and even the smallest probabilities of the order of are estimated with extremely high accuracy (errorbar of the order of ) in a calculation containing Monte Carlo measures. In order to appreciate this result, let us mention that the simple estimator given in Eq. (12) can hardly see a single event occurring with such low probabilities. The reason for this is that the total denominator for the present case is given by and the expansion order (for the case of ), thus it is of the order of . Consequently, events of probabilities will typically never be seen in a Markov chain of length .
Let us now combine the result in Fig. 3 with the result for the calculation of the replica correlation entropy to obtain the entanglement Rényi entropy according to equation (3). Fig. 4 shows the comparison of this result with the one obtained by ED. The correspondence is again perfect and shows that the method works very well: for instance, the evenodd oscillations for are perfectly reproduced Calabrese et al. (2010). Note that the errorbars of the two Monte Carlo results yield the final error of by
(24) 
We usually find that the error bar of the replica correlation entropy is larger than the error of the PR entropy due to the lack of an improved estimator of and therefore dominates the total error.
v.1.2 Thermodynamic Rényi entropies
We here now consider the fullsystem made of the chain to obtain the thermodynamic Rényi entropy. We have performed calculations with replicas for the periodic Heisenberg chain at different finite temperature for different system sizes and calculated the participation entropies . A second set of simulations at inverse temperatures has then been carried out in order to obtain the replica correlation using Eq. (21). We extracted the thermodynamic Rényi entropies and from the result and compare to ED in Figure 5. Clearly the QMC result for matches the exact values perfectly. Furthermore, one can access much larger chain sizes as compared to ED techniques, limited to for the full diagonalization required to access finite temperature behaviour.
We also show the temperature dependence of the individual terms (see inset of Fig. 5) from which the thermodynamic Rényi entropy is obtained. While the participation entropy decreases with inverse temperature (it assumes its maximal value of at ), the replica correlation increases to eventually match the value of the PR entropy at zero temperature.
One can also test with chains the conformal field theory prediction for . In the regime (and ignoring logarithmic corrections due to marginal operators Cardy (1986)), the free energy obeys the following scaling Affleck (1986)
(25) 
where is the groundstate energy, the central charge, and the velocity of excitations. Using Eq. (19), one arrives for the low temperature scaling to:
(26) 
This behavior is checked with for XXX chains of various lengths in Fig. 6 (left), where the low temperature linear form is well reproduced using and . Finite size convergence effects are due to the finite length gap such that the asymptotic low T behavior Eq. (26) is expected to be valid for . Below this gap, displays an activated shape, controlled by . We checked this finitesize effect using ED at the freefermion point (XX chain) with open boundary conditions^{2}^{2}2The use of open chains simplifies the problem of computing thermal averages in the free fermion representation of the XX spin chain. where the asymptotic linear scaling is perfectly well reproduced for large enough sizes , as displayed in Fig. 6 (right).
v.2 Heisenberg ladders
Let us now consider Heisenberg ladders consisting of two neighboring one dimensional periodic Heisenberg chains (the “legs”) with an additional “rung” coupling between the chains:
(27) 
where is the spin operator on site of chain , corresponding to the left and right leg respectively (see Fig. 7). We use periodic boundary conditions along the legs.
For the calculation of entanglement properties, we consider the cut where is one leg of the ladder and perform calculations in the strongly gapped rungsinglet regime , where entanglement entropies are known to be quite large from ED studies Poilblanc (2010); Läuchli and Schliemann (2012). The motivation for this regime is to test our method in a difficult, largeentanglement, regime. Such a cut has also been used in several other works on ladder systems Poilblanc (2010); Peschel and Chung (2011); Cirac et al. (2011); Läuchli and Schliemann (2012); Schliemann and Läuchli (2012); Lundgren et al. (2012); Chen and Fradkin (2013); Lundgren et al. (2013).
v.2.1 Entanglement entropies
Fig. 8 displays our QMC result for various values of , system sizes ranging from to and . For comparison, we also display the numerically exact DMRG result for . We are still able to perform the calculation for up to as large as for and begin to see limitations at for as the errorbar becomes larger. Clearly, the situation becomes worse for , while the result for smaller values of remains extremely good. For comparison, ED (due to the Hilbert space size) or DMRG (due to the large entanglement in this regime) cannot reach systems larger than .
Interestingly, the finite size effects on strongly depend on the Rényi index . For , no difference between the result for and the one for is visible, however, for , displays a sizeable finite length dependence. This can be easily understood if one realizes that the Rényi index plays the role of an inverse temperature in the entanglement spectrum. This behavior points to a stronger finite size dependence of the lowest lying level of the entanglement spectrum (i.e. the groundstate energy of the entanglement Hamiltonian — see discussion later) than for high temperature quantities, which are averaged over the whole spectrum.
v.2.2 Comparison with the mixed ensemble method
In order to get an estimate of the efficiency of the method discussed in this article, we performed calculations for the , Heisenberg ladder, where subsystem corresponds to one leg of the ladder (Fig. 7), and compare to results obtained using the method of Humeniuk and Roscilde Humeniuk and Roscilde (2012) where for every , we optimized the subsystem increment used for the ratio trick Hastings et al. (2010).
Using the same amount of CPU time ( for and for for our method), we compare the values of the errorbars between the two different methods. The results are shown in Table 1. Clearly, for small values of , the error bars obtained from our method are reduced by one order of magnitude. For very large values of , the situation changes and the ratio trick provides a better accuracy, leading to an errorbar that is roughly times smaller. This comes from the fact that computing the replica correlation part scales exponentially with , and, as it was explained previously, no improved estimator is available. For practical purposes, it means that we are limited to values of such that (corresponding to in this particular case) and the associated error will dominate the total error. On the other hand, in the mixed ensemble calculation the ratio trick offers some flexibility in a certain range of . For , the best error bar is generally obtained for an increment larger that one, whereas for larger it becomes quickly much more efficient to set it to unity. Thus, it becomes computationally more interesting to make several simulations for which the probabilities are larger, yielding an exponential gain, while the cpu time required by the increasing number of simulations in order to maintain a constant error bar increases as ( the number of increments). This explains why the mixed ensemble method Humeniuk and Roscilde (2012) becomes more efficient when grows.
It should be noted that this comparison is rather rough, as all three implementations have slightly different optimization goals. Additionally, we did not optimize the CPU time ratio between the calculation of and , which can certainly lead to some improvement. For the calculation of we used replicas and obtained the result for in one single simulation, while in the other two simulations, every has to be done separately. Therefore the comparison gives a slight advantage to the method of Humeniuk and Roscilde as the calculation provides more information (on and ) than needed.
Let us finally mention that the biggest advantage of the method proposed in this article is found in situations of very large entropies, such as the example of the ladders presented here. For the case of weak entanglement entropies such as the one dimensional Heisenberg chain, we obtained roughly the same errorbars in both methods for , indicating that the mixed ensemble method performs very efficiently here.
v.3 Entanglement spectrum reconstruction
Having access to the entanglement entropies for various values of the Rényi index, we are able to explore a recently proposed method for the reconstruction of the entanglement spectrum from Rényi entanglement entropies measured in QMCSong et al. (2012); Chung et al. (2013). The method relies on the NewtonGirard identities, linking the coefficients of a polynomial to the power sums of its roots. This means that a polynomial with roots at the corresponding to the entanglement spectrum can be constructed from the knowledge of Rényi entanglement entropies ( are the eigenvalues of the reduced density matrix).
However, in a practical QMC calculation, the Rényi entanglement entropies are neither known for arbitrarily many values of nor to unlimited precision. Therefore, the polynomial has to be truncated and its order is limited to the maximal Rényi index . For the onedimensional extended BoseHubbard model, Chung et al. Chung et al. (2013) obtained interesting results for the low lying entanglement spectrum using . Here, we perform a similar calculation for the entanglement spectrum of a Heisenberg ladder (see Fig. 7, using the same bipartition as previously) for and focus on the role of systematic and statistical errors. As in Ref. Chung et al., 2013, we split the reduced density matrix in its symmetry sectors and perform the calculation in each sector (see Appendix B).
We have calculated the Rényi entanglement entropies up to with the QMC method described above for the sectors , and ( decreases with due to too large entropies) and reconstructed the entanglement spectrum, systematically varying in order to demonstrate the rate of convergence towards the exact entanglement spectrum obtained from DMRG. In Fig. 9, we also show the reconstructed spectrum from the exact DMRG entanglement entropies. The deviation of the reconstructed DMRG spectrum from the exact spectrum gives an impression of the systematic error due to the truncation of the polynomial, while the reconstructed QMC spectrum carries additional errors due to the statistical uncertainty of the entanglement entropies.
The convergence of the lowest level (largest eigenvalue of ) in the sector is very good and the QMC result is trustworthy. This result corresponds to the single copy entanglement entropy for which no direct QMC estimate is available. The lowest level in the sector also seems to be converged (within errorbars), however, judging from the QMC data only, it is not possible to decide whether the result is trustworthy or not. All other levels can not be trusted due to the statistical uncertainty. Therefore, one should bear in mind that a careful convergence analysis with has to be carried out and that in general only the lowest part of the entanglement spectrum can be extracted from QMC data bearing statistical uncertainties. The reconstructed spectrum from the exact DMRG entanglement entropies shows in particular, how slowly higher levels of the entanglement spectrum converge with . This is the reason why it is beneficial to split the calculation in symmetry sectors.
A global view on the full entanglement spectrum, resolved in spin sectors, is provided in Fig. 10 (left) where the exact entanglement levels from DMRG are compared to the reconstruction from QMC data. A parabolic enveloppe is shown, as expected from the low energy spectrum of XXZ chains Alcaraz et al. (1987). The right panel of Fig. 10 shows that besides the reconstruction method, it is also possible to accept the single copy entanglement using an extrapolation in of the Rényi entanglement entropies .
v.4 Entanglement Hamiltonian
Given the difficulty of the extraction of the entanglement spectrum from QMC data, different methods of the verification of effective entanglement Hamiltonians and the extraction of the inverse entanglement temperature should be explored. We have proposed the usage of the participation spectrum for this purpose in a previous workLuitz et al. (2014c), which can only provide partial proof that the effective model is correct. Here, we propose a different and in fact complementary method that relies on the comparison of the Rényi entanglement entropy and the finite temperature thermodynamic Rényi entropy of a putative entanglement Hamiltonian.
The reduced density matrix can be expressed as a thermal mixed state of an effective entanglement Hamiltonian at inverse temperature by
(28) 
with . Therefore, the thermodynamic entropies of the effective model at inverse temperature have to be equal to the entanglement entropies for all , which is obvious from the above definition Eq. (28) of .
In the case of the strongly entangled Heisenberg ladder (), it was shown within first order perturbation theory Peschel and Chung (2011); Läuchli and Schliemann (2012); Schliemann and Läuchli (2012); Chen and Fradkin (2013) that the entanglement Hamiltonian is given by the simple Heisenberg chain at an effective inverse temperature of . The second order correction for the effective temperature (cf. Eq. (21) in Ref. Schliemann and Läuchli, 2012, ignoring next nearest neighbor interactions) yields
(29) 
In Figure 11, we explore the regime of validity of the firstorder entanglement Hamiltonian. We display of the Heisenberg chain for different Rényi indices as a function of inverse temperature . Then, we calculate the entanglement entropy of the Heisenberg ladder for different values of with our QMC method and extract the effective inverse temperature for which the two quantities match. If the entanglement Hamiltonian is correct, the resulting effective inverse temperature has to be independent of . This is clearly the case if becomes large, as visible both in main panel and inset of Fig. 11, where the deviation from Eq. 29 is getting smaller and flatter (as a function of ) when increases.
Vi Conclusion
We have shown that the calculation of entanglement Rényi entropies may be split into two independent Monte Carlo simulations, one of which boils down to a standard calculation of the participation Rényi entropy obtained from a simulation of independent replicas, while the other part is a “replica correlation” entropy obtained in a simulation of replicas glued together on subsystem . As the PR entropy is a basis dependent quantity, has to be basis dependent, too. Both quantities have to be calculated in the same basis to obtain the correct entanglement entropy.
In a second step, we have developed an improved estimator for the PR entropy, exploiting the fact that the independent replicas can be transformed independently under imaginary time and space symmetry transformations leaving the weight of the Monte Carlo configuration invariant. The number of configurations that are averaged over is therefore multiplied by a number growing exponentially with the number of replicas and counteracts the exponential decay of the probability of finding identical states in all replicas with . This improved estimator allows us to measure extremely low probabilities, crucial for the calculations of for large values of the Rényi index .
Given that the two terms and exhibit a volume law, the realm of applicability of the method is limited to cases, where is not too large. This is precisely the case in situations where the circumference of subsystem is identical to its volume such as for the example of the ladders studied in this article. Here, the largest contribution to the entanglement entropy stems from the PR entropy, which can be calculated with very good precision due to the improved estimator.
For situations where the volume of the subsystem becomes large, such as the half system of a two dimensional lattice, it is possible to combine the two methods to calculate the entanglement entropy. This would result for example to perform the first increment (i.e. a lineshaped subsystem) with the method presented here with very good precision and start from there exploiting a ratio trick Hastings et al. (2010); de Forcrand et al. (2001) using the method introduced in Ref. Humeniuk and Roscilde, 2012. This way, the largest growth of entanglement entropy is dealt with by the improved estimator and the addition of further lattice sites to the subsystem does not increase the entanglement dramatically, therefore the ratio trick method is supposed to work very well without accumulating larger errors.
Let us finally note that the lessons from the improved estimator can certainly be implemented in the method introduced by Humeniuk and RoscildeHumeniuk and Roscilde (2012). If the QMC configuration is in the independent ensemble, one can check if any of the symmetry transformed replica configurations introduced here matches the gluing condition. If so, one has to actually perform the transformation of the whole operator string: in practice, this may turn computationally expensive and one would need to check in which situations the improvement in statistics will be worth the additional computational extracost.
Vii Acknowledgements
We wish to thank S. Capponi for discussions and for providing test data, D. Poilblanc, S. Pujari and S. Wessel for useful related discussions and I. McCulloch for providing access to his code ^{3}^{3}3See http://physics.uq.edu.au/people/ianmcc/mptoolkit/ used to perform the test DMRG calculations. Our QMC codes are partly based on the ALPS libraries Albuquerque et al. (2007); Bauer et al. (2011). This work was performed using HPC resources from GENCI (grant x2014050225) and CALMIP (grant 2014P0677) and is supported by the French ANR program ANR11IS0400501. We also acknowledge the technical assistance by the IDRIS supercomputer center.
Appendix A Improved estimator of thermal Rényi entropies
In order to be explicit, let us start from Eq. (19) and concentrate on the ratio of partition functions that have to be calculated. Here, we will use equation (263) from Ref. Sandvik, 2010 for the stochastic series expansion of the partition function in terms of an operator string composed of bond Hamiltonians linking the lattice sites and .
(30) 
Let us now introduce the observable to compute the PR entropy in the ensemble of independent replicas. It is given by the Kronecker delta .
(31) 
where in the last step, we reexpressed the independent replicas in terms of a unique system at inverse temperature . This is done by introducing a partition of the operator string into parts such that the cutoffs sum up to : . Clearly, in all the expressions above, the numbers of (offdiagonal) operators () in slice of the operator string also have to sum up to the complete number of (offdiagonal) operators (). The introduction of the Kronecker delta expresses the fact that only operator strings in which the states at the end of each slice are identical will contribute to our result.
Alltogether, we see that we have to perform an SSE calculation at inverse temperature and measure the observable:
(32) 
Appendix B Rényi entanglement entropies by sector
Here, we provide some additional details on the reconstruction of the entanglement spectrum by symmetry sectors (here sectors) of the reduced density matrix using the method described in Refs. Song et al., 2012; Chung et al., 2013. Since the reduced density matrix is blockdiagonal with of subsystem , we can split the calculation into the blocks .
The first observation is that the normalization of each block is no longer given by , but by the probability of finding a state with subsystem magnetization equal to :
(33) 
which is easily measured in the independent ensemble. This is of course the first power sum of the eigenvalues of .
We can calculate the th power sum of by using the same method as described in the main text by just ignoring in the measurement states that are not in the correct sector. However, this would corresponds to a reduced density matrix which is normalized to . In order to obtain the correct normalization, we calculate
(34) 
where and correspond to the transition probabilities estimated from measurements in the corresponding sector only. Note that these probabilities can be obtained by binning the measurements of the transition probabilities by their sectors without changing their normalization and then by dividing by the probability of being in the correct sector in the corresponding ensemble (for the independent ensemble, this cancels exactly the factor ).
From the knowledge of the power sums of the eigenvalues of we can now use the method described by Song et al. Song et al. (2012) to reconstruct the entanglement spectrum. Note however, that in equation (2.30) of Ref. Song et al., 2012 the on the diagonal has to be replaced by the first power sum of the eigenvalues of , i.e. by .
References
 Einstein et al. (1935) A. Einstein, B. Podolsky, and N. Rosen, Phys. Rev. 47, 777 (1935), URL http://link.aps.org/doi/10.1103/PhysRev.47.777.
 Calabrese et al. (2009) P. Calabrese, J. Cardy, and B. Doyon, Journal of Physics A: Mathematical and Theoretical 42, 500301 (2009), URL http://stacks.iop.org/17518121/42/i=50/a=500301.
 Jin and Korepin (2004) B.Q. Jin and V. Korepin, Journal of Statistical Physics 116, 79 (2004), ISSN 00224715, URL http://dx.doi.org/10.1023/B%3AJOSS.0000037230.37166.42.
 Fan et al. (2004) H. Fan, V. Korepin, and V. Roychowdhury, Phys. Rev. Lett. 93, 227203 (2004), URL http://link.aps.org/doi/10.1103/PhysRevLett.93.227203.
 Calabrese and Cardy (2004) P. Calabrese and J. Cardy, J. Stat. Mech. 2004, P06002 (2004), ISSN 17425468, URL http://iopscience.iop.org/17425468/2004/06/P06002.
 Refael and Moore (2004) G. Refael and J. E. Moore, Phys. Rev. Lett. 93, 260602 (2004), URL http://link.aps.org/doi/10.1103/PhysRevLett.93.260602.
 Schollwöck (2005) U. Schollwöck, Rev. Mod. Phys. 77, 259 (2005), URL http://link.aps.org/doi/10.1103/RevModPhys.77.259.
 Kallin et al. (2011) A. B. Kallin, M. B. Hastings, R. G. Melko, and R. R. P. Singh, Phys. Rev. B 84, 165134 (2011), URL http://link.aps.org/doi/10.1103/PhysRevB.84.165134.
 Singh et al. (2012) R. R. P. Singh, R. G. Melko, and J. Oitmaa, Phys. Rev. B 86, 075106 (2012), URL http://link.aps.org/doi/10.1103/PhysRevB.86.075106.
 Kallin et al. (2013) A. B. Kallin, K. Hyatt, R. R. P. Singh, and R. G. Melko, Phys. Rev. Lett. 110, 135702 (2013), URL http://link.aps.org/doi/10.1103/PhysRevLett.110.135702.
 Kallin et al. (2014) A. B. Kallin, E. M. Stoudenmire, P. Fendley, R. R. P. Singh, and R. G. Melko, arXiv eprint 1401.3504 (2014), URL http://arxiv.org/abs/1401.3504.
 Alet et al. (2007) F. Alet, S. Capponi, N. Laflorencie, and M. Mambrini, Phys. Rev. Lett. 99, 117204 (2007), URL http://link.aps.org/doi/10.1103/PhysRevLett.99.117204.
 Kallin et al. (2009) A. B. Kallin, I. González, M. B. Hastings, and R. G. Melko, Phys. Rev. Lett. 103, 117203 (2009), URL http://link.aps.org/doi/10.1103/PhysRevLett.103.117203.
 Hastings et al. (2010) M. B. Hastings, I. González, A. B. Kallin, and R. G. Melko, Phys. Rev. Lett. 104, 157201 (2010), URL http://link.aps.org/doi/10.1103/PhysRevLett.104.157201.
 Melko et al. (2010) R. G. Melko, A. B. Kallin, and M. B. Hastings, Phys. Rev. B 82, 100409 (2010), URL http://link.aps.org/doi/10.1103/PhysRevB.82.100409.
 Humeniuk and Roscilde (2012) S. Humeniuk and T. Roscilde, Phys. Rev. B 86, 235116 (2012), URL http://link.aps.org/doi/10.1103/PhysRevB.86.235116.
 Grover (2013) T. Grover, Phys. Rev. Lett. 111, 130402 (2013), URL http://link.aps.org/doi/10.1103/PhysRevLett.111.130402.
 Assaad et al. (2014) F. F. Assaad, T. C. Lang, and F. Parisen Toldin, Phys. Rev. B 89, 125121 (2014), URL http://link.aps.org/doi/10.1103/PhysRevB.89.125121.
 Chung et al. (2013) C.M. Chung, L. Bonnes, P. Chen, and A. M. Läuchli, arXiv:1305.6536 [condmat] (2013), URL http://arxiv.org/abs/1305.6536.
 Helmes and Wessel (2014) J. Helmes and S. Wessel, ArXiv eprints (2014), eprint 1403.7395.
 Herdman et al. (2014) C. M. Herdman, P.N. Roy, R. G. Melko, and A. Del Maestro, Phys. Rev. B 89, 140501 (2014), URL http://link.aps.org/doi/10.1103/PhysRevB.89.140501.
 Sandvik (2010) A. W. Sandvik, AIP Conference Proceedings 1297, 135 (2010), ISSN 0094243X, URL http://proceedings.aip.org/resource/2/apcpcs/1297/1/135_1.
 Chung and Peschel (2001) M.C. Chung and I. Peschel, Phys. Rev. B 64, 064412 (2001), URL http://link.aps.org/doi/10.1103/PhysRevB.64.064412.
 Audenaert et al. (2002) K. Audenaert, J. Eisert, M. B. Plenio, and R. F. Werner, Phys. Rev. A 66, 042327 (2002), URL http://link.aps.org/doi/10.1103/PhysRevA.66.042327.
 Inglis and Melko (2013) S. Inglis and R. G. Melko, Phys. Rev. E 87, 013306 (2013), URL http://link.aps.org/doi/10.1103/PhysRevE.87.013306.
 Luitz et al. (2014a) D. J. Luitz, F. Alet, and N. Laflorencie, Phys. Rev. Lett. 112, 057203 (2014a), URL http://link.aps.org/doi/10.1103/PhysRevLett.112.057203.
 Luitz et al. (2014b) D. J. Luitz, F. Alet, and N. Laflorencie, Phys. Rev. B 89, 165106 (2014b), URL http://link.aps.org/doi/10.1103/PhysRevB.89.165106.
 Luitz et al. (2014c) D. J. Luitz, N. Laflorencie, and F. Alet, arXiv eprint 1403.3717 (2014c), URL http://arxiv.org/abs/1404.3713.
 de Forcrand et al. (2001) P. de Forcrand, M. D’Elia, and M. Pepe, Phys. Rev. Lett. 86, 1438 (2001), URL http://link.aps.org/doi/10.1103/PhysRevLett.86.1438.
 Stéphan et al. (2009) J.M. Stéphan, S. Furukawa, G. Misguich, and V. Pasquier, Physical Review B 80, 184421 (2009), ISSN 10980121, 1550235X, URL http://link.aps.org/doi/10.1103/PhysRevB.80.184421.
 Song et al. (2012) H. F. Song, S. Rachel, C. Flindt, I. Klich, N. Laflorencie, and K. Le Hur, Phys. Rev. B 85, 035409 (2012), URL http://link.aps.org/doi/10.1103/PhysRevB.85.035409.
 Stéphan et al. (2010) J.M. Stéphan, G. Misguich, and V. Pasquier, Physical Review B 82, 125455 (2010), ISSN 10980121, 1550235X, URL http://li