Unusual Corrections to Scaling and Convergence
of Universal Renyi Properties at Quantum Critical Points
Abstract
At a quantum critical point, bipartite entanglement entropies have universal quantities which are subleading to the ubiquitous area law. For Renyi entropies, these terms are known to be similar to the von Neumann entropy, while being much more amenable to numerical and even experimental measurement. We show here that when calculating universal properties of Renyi entropies, it is important to account for unusual corrections to scaling that arise from relevant local operators present at the conical singularity in the multi-sheeted Riemann surface. These corrections grow in importance with increasing Renyi index. We present studies of Renyi correlation functions in the transverse-field Ising model (TFIM) using conformal field theory, mapping to free fermions, and series expansions, and the logarithmic entropy singularity at a corner in for both free bosonic field theory and the TFIM, using numerical linked cluster expansions. In all numerical studies, accurate results are only obtained when unusual corrections to scaling are taken into account. In the worst case, an analysis ignoring these corrections can get qualitatively incorrect answers, such as predicting a decrease in critical exponents with the Renyi index, when they are actually increasing. We discuss a two-step extrapolation procedure that can be used to account for the unusual corrections to scaling.
I Introduction
Quantum entanglement in the ground state of a many-body system contains universal signatures of a variety of low-energy and long-lengthscale physical phenomena. These include Goldstone modes due to spontaneous breaking of a continuous symmetry Kallin et al. (2011); Metlitski and Grover ; Luitz et al. (2015); Kulchytskyy et al. (2015); topological order Kitaev and Preskill (2006); Levin and Wen (2006); Isakov et al. (2011); Fermi surfaces Wolf (2006); Gioev and Klich (2006); quantum criticality Metlitski et al. (2009); Casini and Huerta (2012); Grover (2014); Bueno et al. (2015); and more. If a system is divided into two spatial regions and , the entanglement in a pure state can be characterized by the von Neumann entropy associated with the reduced density matrix of either subsystem,
(1) |
However, numerically computing the von Neumann entropy is a significant challenge, requiring an explicit representation of the reduced density matrix or the ground state wavefunction. In continuum field theories, the von Neumann entropy can be computed from the Renyi entropies,
(2) |
by introducing identical replicas of the system, then formally taking the limit Calabrese and Cardy (2004). Renyi entropies for integer values of can be calculated by a variety of methods employing this replica-trick, however, the limit is difficult to obtain from numerical data for finite, discrete , which can often be most accurate for small .
For this reason, integer Renyi entropies are often studied in quantum many-body systems on their own merit. A variety of recent numerical techniques have been developed for this task. One way to calculate them in Quantum Monte Carlo (QMC) simulations is by evaluating the expectation value of a swap operator between the different replicas Hastings et al. (2010). Alternatively, these integer Renyi entropies can be expressed as a ratio of two partition functions Melko et al. (2010) and evaluated by stochastically sampling an extended ensemble where one switches between the configurations defining the two partition functions Humeniuk and Roscilde (2012). Such sampling techniques make the Renyi entropies available in other numerical approximations of the wave function, such as tensor network states Wang et al. (2013). They can also be calculated by series expansions in various coupling constants of the model Singh et al. (2011); Kallin et al. (2011); Singh et al. (2012); Devakul and Singh (2014a, b). Most importantly, there are a number of proposals Cardy (2011); Abanin and Demler (2012); Daley et al. (2012) to measure the Renyi entropy experimentally.
This manuscript considers the universal properties associated with the Renyi entropies and associated Renyi correlation functions in a critical system. While Renyi entropies do not satisfy strong subadditivity a priori, and thus may not provide as rigorous an information measure of entanglement as Eq. (1), explicit calculations have found that their singular behavior is very analogous to the von Neumann entropy in typical many-body systems. Universal terms associated with Goldstone modes, topological order, Fermi surfaces and quantum critical points often depend in a simple parametric way on the Renyi index . Hence, computing universal terms in any one Renyi entropy can be sufficient to determine the universality class of the system. If there is an entanglement property that exhibits monotonicity under renormalization Zamolodchikov (1986), one would expect it to remain valid for the Renyi entropies as well.
However, there is a subtlety in obtaining universal, asymptotic Renyi properties, which is the main focus of this work. Renyi entropies, and associated correlation functions, suffer from “unusual” corrections to scaling that are absent for von Neumann entropies Cardy and Calabrese (2010). These corrections become more and more severe as the Renyi index increases and cause a slow convergence to the asymptotic limit. Rapid convergence in the von Neumann entropy and slower convergence for higher Renyi cases were previously noted for topological entanglement entropy, where they were found to be related to non-universal perturbative terms Jiang et al. (2013). Here, we discuss a two step extrapolation process needed to obtain the asymptotic universal critical Renyi properties.
In light of our work, there may be a need to revisit previous calculations of Renyi entropies and their universal terms in past computational studies. For example, in this paper we re-examine recent numerical studies of models in dimension, where it was found that the coefficients of the logarithmic singularities associated with a corner approximately scales with Kallin et al. (2013, 2014); Stoudenmire et al. (2014); Helmes and Wessel (2014). These studies also found significant differences between the entropy of free fields and a single copy of the interacting theory (differences of order of greater than 10 percent). In light of a recent conjecture that the corner entanglement mimics a central charge Bueno et al. (2015), one might expect a much smaller difference. Here, we find that some of the observed difference could be due to the unusual corrections to scaling. Our work also suggests that one should be especially careful in computing universal properties at large values of the Renyi index .
Following a discussion of how unusual corrections to Renyi properties arise in Section II, we illustrate their importance in Section III by studying Renyi correlation functions in (space+time) dimension , for which there is a well-developed theory. In Section IV we apply these insights in to universal entanglement entropy contributions from sharp corners. Accounting for the unusually strong Renyi finite-size corrections significantly improves our numerical results for the universal corner term for the free boson theory; we apply the same procedure to improve our estimates of the corner term for the critical Ising model.
Ii Unusual corrections to scaling
Corrections to critical scaling are well known for bulk quantities and are due to the presence of irrelevant operators Wegner (1972). When approaching a critical point, a thermodynamic quantity such as the order parameter susceptibility does not behave as a single power-law, but as a sum of power laws. For small non-zero ,
(3) |
or for a finite system of size at the critical point,
(4) |
The subleading terms above lead to violations of simple scaling and data collapse, such as the failure of Binder ratios to cross precisely at a critical point. The subleading terms can be seen to directly affect the calculation of exponents, by taking a logarithm of Eq. (3) and differentiating with respect to for small ,
(5) |
Thus the effective exponent depends on and, by scaling theory, on ,
(6) |
Note that the same is true for the coefficient of a logarithmic singularity.
When studying lattice theories in the bulk, such as the -dimensional classical Ising model, including or ignoring such conventional scaling corrections in the analysis of finite-size data typically makes a difference of less than one percent in estimates of critical exponents such as . For this reason, analyses of Monte Carlo simulation or series expansion results often ignore subleading corrections to finite-size scaling. However, as first discussed by Cardy and Calabrese Cardy and Calabrese (2010) in the context of -dimensional conformal field theories (CFTs), unusual corrections to scaling arise in the calculation of Renyi entropies. This can be understood in the path-integral picture, where the entropy is proportional to the partition function of a -sheeted Riemann surface Holzhey et al. (1994); Calabrese and Cardy (2004) (Fig. 1), which at finite-temperature involves a branch cut located at for between each of the replicas defined on region . A conical singularity exists on the boundary between regions and where this branch cut intersects the -periodic imaginary time structure of the path integral in region . Local operators, which might even be relevant in the bulk, are introduced at this conical singularity in the multi-sheeted Riemann surface. Even though their contribution is only local, these operators yield subleading corrections to scaling with a power that nevertheless diminishes with the Renyi index and ultimately goes to zero as goes to infinity, strongly altering the scaling behavior.
The unusual correction to scaling exponents in are well known Cardy and Calabrese (2010) and have been clearly observed in numerical and analytical lattice calculations of spin chains Fagotti and Calabrese (2010, 2011); Calabrese et al. (2010a); Calabrese and Essler (2010); Calabrese et al. (2010b). For the case of the transverse-field Ising model, which we use as an example in the next section, the most relevant operator that does not break the symmetry is the energy operator with a scaling dimension . For the Renyi correlation functions (defined in the next section) it leads to a finite-size correction to scaling like that of Eq. (4) but with an exponent,
(7) |
which goes to zero for large .
While there is no theory for the unusual corrections to scaling in , it is reasonable to expect that they are again dominated by the presence of the energy operator at the conical singularity. The scaling dimension of the energy operator is related to the correlation length exponent by the relation Cardy (1996),
(8) |
In particular, for the free boson theory, and approximately for the Ising theory Kos et al. (2014). In there is no particular reason why the correction should be given by the simple formula , as in . In Section IV, we will see that while the exponent decreases with , the numerical results suggest that the dependence on is slightly different.
Iii Renyi correlators in
To illustrate a case where failing to account for the unusual corrections to scaling can lead to qualitative changes in finite-size scaling, consider the transverse-field Ising model (TFIM) in spatial dimension . The Hamiltonian is,
(9) |
where are the Pauli operators. We focus on the critical point , or its vicinity. The effect of unusual corrections to scaling has been well studied in this model for the Renyi entropies Calabrese and Essler (2010) and related quantities Alba (2013); Coser et al. (2014). Here, we instead examine the scaling of estimators weighted by the Renyi moments. One advantage of looking at these Renyi correlations is that one can observe more directly how the conical singularity affects correlations, not just the partition function (the entropy). As we will demonstrate, these correlation can easily be computed with conformal field theory techniques.
To define the Renyi correlators, we divide the system into two halves and and calculate the expectation value of operators (or pairs of operators) in the Renyi multi-sheeted geometry . The expectation value of an observable in subsystem is defined as,
(10) |
where is the reduced density matrix for subsystem and is the Renyi index. Note that setting gives us normal bulk expectation values, which can also be obtained without partitioning of the system, providing a non-trivial check on our computational procedures.
iii.1 Renyi correlators in conformal field theory
We consider the simplest setup, which is a subsystem of length in an infinite system. As usual for such computations we use the replica trick: is first assumed to be an integer, but the formulae are expected to hold for any . As explained in e.g. Refs. Holzhey et al. (1994); Calabrese and Cardy (2004), the entropy is – up to some normalization – nothing but the partition function of a -sheeted Riemann surface . The Riemann surface may be mapped onto the complex plane by the conformal transformation,
(11) |
Here, is a complex coordinate on the Riemann surface , while is a complex coordinate in the complex plane . Subsystem corresponds to real in the interval .
To obtain a Renyi correlator in the ground state, what we need is to compute the correlation on . Let us do that for some two point function of a primary operator. We use the two point function of a primary operator with scaling dimension , together with the transformation law of a primary operator Di Francesco et al. (1997). We obtain,
(12) |
Here (resp. ) has to be understood as (resp. ). Using (11) this becomes,
(13) |
Eq. (13) is the main result of this section. The Renyi correlations are translational invariant for , as they should be. This is also the case far from the boundaries for any . Similar to the entanglement entropy Calabrese and Cardy (2004), various generalizations to finite systems and finite temperatures are straightforward ^{1}^{1}1For example the result for an interval of size in a periodic system of size follows from Eq. (12) with the conformal mapping .
In the following we are going to check our results in the Ising chain. The most natural example is the case when one of the points lies at the end of the interval ( is a UV cutoff of the order of a lattice spacing; it ensures that the correlations stay finite), and the other at some distance from it:
(14) | |||||
(15) |
where is some -dependent UV cutoff. The correction terms lead to the correction to scaling exponent in Eq. (7). We are going to test formula (15) in the following subsections. Note also, we recover that the scaling dimension of an operator near the cone is modified to (Eq. (7)), as pointed out in Ref. Cardy and Calabrese (2010).
iii.2 A numerical check using free fermions
It is well known that the Hamiltonian (9) can be mapped onto a system of free fermions through a Jordan-Wigner transformation,
(16) | ||||
(17) |
and then diagonalized by a Bogoliubov transformation. The correlations in the ground state may then be obtained analytically Lieb et al. (1961).
It is also well known how to compute the entanglement entropy for such quadratic fermion systems Chung and Peschel (2001); Peschel (2003). The method relies on the following observation. Since Wick’s theorem holds for any fermionic correlation in subsystem , the reduced density matrix (RDM) itself is the exponential of a quadratic fermion Hamiltonian , . Finding the precise form of can then be done by demanding that it reproduce all correlation functions in subsystem , . To perform the identification one of the easiest way is to use the correlation matrix,
(18) |
where (resp. ) denotes the transpose (resp. conjugate transpose) of any matrix . The matrices and encode the two types of correlators:
(19) | |||||
(20) |
The eigenvalues of are simply related to the single particle eigenenergies of the entanglement Hamiltonian Peschel (2003), and the knowledge of the former allows to determine the latter. Said differently, given a subsystem correlation matrix it is easy to determine the corresponding entanglement Hamiltonian, and then to compute the entropy.
Here we want to compute the Renyi correlations, as defined in Eq. (10). Since the RDM is the exponential of a quadratic form, so is any power of it. The new entanglement Hamiltonian corresponding to the new RDM is then simply . One can then find the modified Renyi correlation matrix that would have as the entanglement Hamiltonian. The details are explained in Appendix A. In compact form the result may be written,
(21) |
Note that in absence of pairing terms , the same formula holds with and replaced by and . Using (21), one can access all Renyi correlations between fermions, , , the spin correlations in the basis being recovered by using the Jordan Wigner transformation (16, 17), and appropriate treatment of the resulting Jordan-Wigner string (see e.g. Ref. Lieb et al., 1961). For the spin-spin correlation at distance we obtain,
(22) |
with and .
We now use (22) to numerically compute the correlation,
(23) |
in the ground state of the chain (9) at the critical point =1. We consider several values of and several distances in a subsystem of length . As can be seen in Fig. 2, the agreement with the prediction (15) with Ising/Onsager exponent is very good. Note however that finite-size effects increase significantly with ; indeed the formula is expected to be valid only in the regime , which for moderately large values of already requires immense system sizes.
Another interesting feature is that the correlation near blows up again at for . One possible interpretation would be that since the scaling dimension near the cone is modified to , the spin configurations near the (fictitious) endpoints at are more constrained, so more likely to match. Since translational invariance is broken, the effect appears similar to that of a boundary with partial reflection.
Of course, the fact that we are able to test our results with good precision relies heavily on the free fermion structure. For more general system the computation of Renyi correlations is more difficult, and the accessible system sizes are much smaller. To illustrate this we look at similar correlations with a different method below.
iii.3 Correlation Function and Structure Factors
Let us calculate correlation functions directly in the high-field (disordered) phase of Eq. (9) by a series expansion in Oitmaa et al. (2006). We define via Eq. (23), where is a boundary site of directly at the interface with and refers to any other site in subsystem . For correlation function at in Eq. (23), the series coefficients will be zero below order . One needs at least powers of the perturbation to get a non-zero correlation between site and site . This means as one goes further away in distance, the series become effectively shorter and shorter. For this reason, analyzing real-space correlations is not very useful. Series extrapolation methods, which are necessary for studying the critical region, simply cannot be applied. Instead, we need to focus on structure factors or correlation sums defined as,
(24) |
It is useful but not necessary to multiply all terms in the sum in Eq. (24) by a factor of 2, as that would be the equal-time structure factor of the 1D system in the absence of a partition, where one has two neighbors at each such distance. This is what we have done in the analysis discussed below.
We describe the series analysis using the method of differential approximants Fisher and Au-Yang (1979); Hunter and Baker (1979) in Appendix B. There, we define as the order of the series used in the estimation of the critical exponent (, , and are individually the orders of three polynomials defining the extrapolation scheme). In other words, is the maximum power used in determining the polynomial coefficients. The order of the series expansion is related to a length scale Affleck and Bonner (1990), thus effective exponents can vary with the order as,
(25) |
Again, theory predicts that the scaling dimension of the bulk spin is given by, . The scaling dimension of the boundary spin with Renyi index is given by,
(26) |
Thus the correlation function should decay for large as,
(27) |
Summing over , the structure factor should diverge on approach to the critical point with an exponent,
(28) |
Since, for the model, the exponent simplifies to,
(29) |
Figure 3(a) shows the estimated structure factor exponent obtained from order series expansions with various choices for the order of one of the polynomials used in the fit. Note that the exponent is converged very well to the theoretical value for (which gives the bulk exponent), but the convergence is poor for larger . In fact, while theory predicts that the exponent should increase with , the series extrapolation estimate appears to systematically decrease.
Figure 3(b) shows plots of the estimated exponent values for different orders of the series. The axis has been scaled differently for different values. The rescaling pushes the data for larger values more and more to the right. Thus for larger one needs to extrapolate farther to reach the asymptotic value. The solid squares are the asymptotic theoretical values and the dashed lines are a guide to the eye showing the results are clearly consistent with theory, only when the unusual corrections to scaling are taken into account.
iii.4 Magnetization in the ordered phase
We now calculate series expansions for the Renyi magnetization at the boundary site of by a series expansion in (a low-field expansion in the ordered phase). This boundary magnetization is defined as,
(30) |
Once again, setting just gives us the usual bulk magnetization that is well known to vanish at the critical point with exponent .
When extrapolating, we do not allow any terms since the magnetization vanishes identically at the critical point with no background value (see Appendix B). Without such terms, these approximants are also known as d-log Pade approximants. The order now equals .
The boundary magnetization vanishes with the scaling dimension of the boundary operator (with )
(31) |
Figure 4(a) and 4(b) show corresponding results calculated up to order for the magnetization. Note, however, that the magnetization series only has even powers, so it is effectively a ten-term series in .
Notice that for the magnetization, the variation of the estimated exponents with order appears to adhere more closely to a dependence. Such terms are always present in addition to the terms, but should usually be weaker. However, we must use even more caution in extrapolating the magnetization exponents to asymptotic values, especially for , as our data is still very far from the final values.
Iv Renyi entropies in corners
Having seen the importance of unusual corrections to scaling for Renyi quantities in space-time dimension , we turn to a case of intense recent interest: universal Renyi entropy corrections for critical systems. In two spatial dimensions, universal quantities appearing in the Renyi entropies have a rich dependence on subregion geometry. For example, various universal terms (sub-leading to the area law) can arise from geometries such as smooth circular bipartitions Casini and Huerta (2009, 2012); sharp corners Casini and Huerta (2007); Fradkin and Moore (2006); Zaletel et al. (2011); Stéphan et al. (2011); bipartitioned infinite cylinders Metlitski et al. (2009) or cylindrically-bifurcated tori Ju et al. (2012); Stéphan et al. (2013); Inglis and Melko (2013); Chen et al. (2015).
Significant computational effort has been dedicated recently toward calculating the subleading contributions to the entanglement entropy arising from corners in quantum models on the square lattice. Here we study the free scalar field theory, as well as the two-dimensional Ising chain in transverse field. A bipartition with a single corner (e.g. a quadrant of the infinite plane) in the entangling boundary contributes a subtractive, divergent logarithmic correction to the Renyi entropy,
(32) |
where is the linear length scale of the bipartition , and is the UV cutoff corresponding to the lattice length-scale. The constant is not polluted by UV details, and constitutes a universal number that can be explored easily for using lattice numerics.
Previous numerical studies reveal that the value of obtained for the model is somewhat close to times the value obtained for a free scalar field theory by Casini and Huerta Casini and Huerta (2007). In this section, we use the numerical linked-cluster expansion (NLCE) to examine the finite-size scaling behavior of for free bosons, and perform a parallel comparison of data obtained for the transverse field Ising model (TFIM) in two spatial dimensions, at its fixed-point. For the latter we use the density matrix renormalization group (DMRG) as our cluster solver Schollwöck (2005); ITe .
In our NLCE procedure, described in detail in Appendix C, the leading-order area-law piece of Eq. (32) cancels, leaving only the logarithmic divergence,
(33) |
In the next two sections, we interpret the order of the NLCE (see Appendix C) as a length scale , a procedure used successfully several times in the past Kallin et al. (2013, 2014); Stoudenmire et al. (2014). This length scale is used to fit the two parameters of Eq. (33) to extract the universal constant .
iv.1 Free scalar field theory
We begin by employing the NLCE method to study the free scalar field theory, described in the Appendices. For the free boson, the NLCE procedure can be used with a correlation-matrix technique as the cluster solver, as detailed in Appendix A. The correlation-matrix technique permits the calculation of Renyi entropies on very large finite-size clusters, allowing the NLCE sum to be carried to extremely high order.
In Fig. 5(a), we extract the universal coefficient from NLCE data over a range of orders, from to , by fitting Eq. (33) as a function of NLCE order with two free parameters ( and ). For reference, Fig. 5(a) includes the exact values for the thermodynamic limit obtained by Casini and Huerta Casini and Huerta (2007). The von Neumann corner coefficient converges quite rapidly to the exact value, while the Renyi coefficients and converge much more slowly, even for very large linear cluster sizes .
As discussed in Section II, one may expect that for a second extrapolation could become necessary because of the unusual corrections to scaling arising from the conical singularities. We perform this second extrapolation in Fig. 5(b). The individual points are obtained by first fitting the NLCE data as a function of order to Eq. (33) from up to a fixed maximum order . This is the same procedure discussed above for producing the curves labeled NLCE in Fig. 5(a), but keeping fixed and combining results for various into a single curve. Next, we identify each , a length scale which depends on the NLCE scheme (see Appendix C). We then fit the finite-size results to the function with and as fitting parameters. From Eq. (8), taking for , one might guess the correction should go as , i.e. , in analogy to the result. However, as we see in Fig. 5(b), the optimized values of which produce the best fit are somewhat different from and yield extrapolations closer to the exact results Casini and Huerta (2007), likely indicating different scaling exponents associated with the higher-dimensional system.
From these extrapolations, we estimate the following values for in the thermodynamic limit: for we obtain which is 1% off of the exact value of Ref. Casini and Huerta (2007); for we obtain which is 3% from the exact value of . Clearly, our NLCE estimates will approach the exact results with increasing accuracy as increases to infinity.
iv.2 Transverse-field Ising model
We turn now to NLCE calculations of the transverse-field Ising model (TFIM) in spatial dimension . The Hamiltonian on a square lattice is
(34) |
and we work at the quantum critical point . The corner contribution to the Renyi entropy has been studied by NLCE several times in the past Kallin et al. (2013, 2014). The main technical difference from the free boson calculation of the preceding section is that each individual cluster must be solved using a method for strongly correlated Hamiltonians. This severely restricts the value of , even if one uses the powerful DMRG method as a cluster solver.
In Fig. 6, the universal coefficient is extracted from NLCE data in a range of orders, from to , using the fitting form Eq. (33) with two free parameters ( and ). Up to the maximum order we reach, the results obtained this way for the TFIM are significantly different from the exact results for the free boson, particularly for . However, if one performs a second finite-size extrapolation for the Renyi coefficients, using the form with values taken to be , where are the optimal exponents for the free boson and is the TFIM energy operator scaling dimension Eq. (8), then our estimates for in the limit are 0.0059 for ; 0.0045 for ; and 0.0040 for . As shown in Fig. 6(a), if we instead assume the exponents , our estimates are 0.0061 for ; 0.0048 for ; and 0.0045 for , higher than the extrapolation based on the boson exponents, especially for larger .
The data indicate that a large part of the observed difference between free bosons and the TFIM in past studies may be the unusual corrections to scaling, which are manifest as large finite-size effects in the NLCE procedure for . When such unusual corrections, arising from the presence of the conical singularity in the Renyi entropies, are taken into account, it becomes much more difficult to distinguish the value of obtained for the TFIM from that for free bosons. Since it will be very difficult to obtain larger values of using the DMRG method, it is hard to judge how much of the remaining discrepancy can be attributed to the very limited cluster sizes of the TFIM NLCE calculation.
iv.3 Discussion
In this paper, we have performed a systematic study of the finite-size scaling of Renyi entropies and Renyi correlators at a quantum critical point. In particular, special attention was paid to the first subleading correction, which can play a crucial role in a proper extraction of universal terms.
In space-time dimension 1+1, a known scaling ansatz has been derived by Calabrese and Cardy and applied previously to numerical studies of Renyi entropies. We have applied this ansatz to Renyi correlators, calculated with conformal field theory, a mapping to free fermions, and series expansions for the transverse-field Ising model (TFIM) at its quantum critical point. For data obtained with finite-size (or finite-series length) calculations, we find that a second extrapolation in the size (or length of the series) using this ansatz is needed to show consistency with the theoretically known values. Failing to take the correction to scaling into account can lead to qualitatively incorrect answers. Thus, it appears that the conical singularities are far more important for Renyi correlators in than irrelevant operators are, e.g. for computing bulk properties of the -dimensional classical Ising model.
We also studied space-time dimensions, where the influence on scaling of the conical singularity is generally not known exactly. We focus on the universal term which arises in Renyi entropies due to a corner, using a recently developed numerical linked cluster expansion (NLCE). In a scalar field-theory corresponding to free bosons, where the exact result is known, very large systems can be studied. As in 1+1, we see that strong corrections to scaling are present for Renyi entropies, which retain a significant discrepancy from exact results even for relatively large cluster sizes. However, if we do an additional extrapolation in the maximum linear dimension as , we find that fitting as a function of produces results that agree well with the exact value.
Finally, following a similar procedure for the corner Renyi entropies of the TFIM at its quantum critical point, we extrapolate existing NLCE data. Using the optimal found for free bosons, but including a non-trivial scaling dimension for the energy operator, we find that the values for are much closer to the free boson values than concluded from previous analyses.
Improving numerical accuracy with such extrapolations is crucially important for the interpretation of many lattice calculations of Renyi entropies. For example, recent calculations on several models tuned to the Wilson-Fisher fixed-point have revealed that is a universal function of that scales to within numerical accuracy as , where is the number of field components, and is some universal function that appears to be the same for all models studied to date. Intriguingly, the central charge also grows linearly with up to small corrections Kos et al. (2014); Bueno et al. (2015). Our results indicate that, when unusual corrections to scaling are taken into account, for can be very close numerically to the value for free bosons, however for other values, the conclusion that remains well justified. Thus, we confirm that the universal subleading corner coefficient is sensitive to number of degrees of freedom in the low-energy field theory describing a strongly-interacting quantum system at its critical point.
Acknowledgments
We would like to acknowledge crucial discussions with J. Cardy, P. Fendley, R. Myers, F. Pollmann, and W. Witczak-Krempa. This research was supported in part by the National Science Foundation Grant No. NSF PHY11-25915 (RRPS, RGM), NSF DMR-1306048 (RRPS) and Grant No. 1006549. The simulations were performed on the computing facilities of SHARCNET and on the Perimeter Institute HPC. Support was provided by NSERC, the Canada Research Chair program, the Ontario Ministry of Research and Innovation, the John Templeton Foundation, and the Perimeter Institute (PI) for Theoretical Physics. Research at PI is supported by the Government of Canada through Industry Canada and by the Province of Ontario through the Ministry of Economic Development & Innovation.
Appendix A Reduced density matrices and entropies for free theories
We recall in this appendix some general results regarding reduced density matrices for free bosonic and fermionic systems. The techniques are well-known; we refer to Ref. Eisler and Peschel (2009) for a review. We first present the method to compute entanglement for free bosons, relevant to Sec. IV.1, before focusing on Renyi correlations (relevant to Sec. III.2).
a.1 Entropies for the free lattice scalar field
It has been long known how to evaluate the entanglement entropy for a system of coupled oscillators, see e.g. Refs. Bombelli et al. (1986); Srednicki (1993). Here we follow the method put forward by Peschel Peschel (2003), that was used by Casini and Huerta Casini and Huerta (2007) to compute exactly the corner contribution mentioned in the text. Numerically, the technique may also be used to compute the entropy on very large clusters.
Consider a finite two-dimensional square lattice, such that a free scalar field and its conjugate momentum exist at each lattice point. A non-interacting system is described by following lattice Hamiltonian,
(35) |
where and are the linear dimensions of the lattice. The Hamiltonian can be rewritten in terms of bosonic operators, but we refrain from doing so for now. We are mainly interested in the massless case , but keep general in all the formulas.
The total number of sites is . The Hamiltonian may be rewritten as
(36) |
The matrix is, up to a diagonal term, nothing but the discrete laplacian. The ground state correlations are given
(37) | ||||
(38) |
and may be obtained explicitly in many geometries, numerically in others.
Let us now cut our system in two parts and , and focus only on the degrees of freedom inside subsystem . Since for any observable in subsystem , the reduced density matrix itself is the exponential of a quadratic boson form
(39) |
where the operators are linear combinations of the original fields . They also satisfy the canonical commutation relations . Here is a normalization constant that ensures . Such an ansatz has by definition the Wick factorization property build in. The next step is then to adjust the energies and the operators so as to reproduce the ground-state correlations, which can always be done. After some algebra, we find that the single particle energies satisfy
(40) |
where the are the positive eigenvalues of the matrix . Here (resp. ) is the matrix (resp. ) restricted to subsystem . We note it , to emphasize the fact that it depends on the choice of subsystem. The Renyi entropy then follows from Eq. (39)
(41) |
This is of course a huge simplification, because what is only needed is the diagonalization of the matrix , whose size is given by the number of sites in . This last step can be done using standard linear algebra routines.
For the geometry we choose the full system to be an open-boundary rectangle of size . The correlators come from diagonalizing the discrete laplacian on the rectangle, and using Eqs. (37, 38). We obtain
(42) | |||||
(43) | |||||
where the quasi-momenta are quantized as
(44) | |||||
(45) |
and
(46) |
Combined with the NLCE procedure described in Appendix C, this produces the data shown in Sec. IV.1. We note that the choice of boundary conditions is not completely innocent. For example, if we were to choose a periodic-boundary torus instead of an open-boundary rectangle, the correlators would be given by similar but translation invariant formulas
(47) | ||||
(48) |
with
(49) |
and the momenta quantized as
(50) | |||||
(51) |
Unfortunately, the correlator (47) is divergent in the massless case , due to the presence of a zero mode . Thus, also diverges, as shown previously in Ref. Metlitski and Grover , significantly complicating the calculation of Renyi entropies on finite-size tori. To remove this divergence, a possible way is keep a very small mass term; this results in an additional term proportional to in the entropy, where is the length of the boundary. This may then be systematically decreased and the results extrapolated to .
This complication is not a fundamental one; however, it motivates us to only consider clusters with open boundary conditions in the procedure of Appendix C, which do not have zero modes.
a.2 Renyi correlations
We consider the following Hamiltonian
(52) |
Such a quadratic fermion form can always be diagonalized by a Bogoliubov transformation. Namely, there is a set of fermion operators
(53) |
which obey the canonical commutation relations , and which diagonalize :
(54) |
The coefficients and the single particle energies may be e.g. obtained by diagonalizing the following Bogoliubov matrix
(55) |
where (resp. ) is the transpose of (resp. conjugate transpose of ). Using this method, computing for example the ground-state correlations is a straightforward task.
We now focus on a subsystem, which we choose to be the first sites for simplicity. For later convenience we introduce the following correlation matrix
(56) |
and encode the two types of correlators:
(57) | |||||
(58) |
Now, any correlation in subsystem may be recovered from the reduced density matrix,
(59) |
Because of Wick’s theorem, the reduced density matrix itself can be written as
(60) |
where is a quadratic fermion form for a system of size , similar to (52). We name the entanglement Hamiltonian, as it is different from the original Hamiltonian . The crucial point is that the Bogoliubov matrix corresponding to and the subsystem correlation matrix are diagonalized by the same transformation Peschel (2003). More precisely, it is possible to find a unitary matrix which diagonalizes
(61) |
and such that
(62) |
with
(63) |
The single particle eigenenergies of the entanglement Hamiltonian may then be obtained from the by requiring consistency with (59). We obtain
(64) |
This relation, together with (63), uniquely determines the entanglement Hamiltonian from the eigenvalues and the eigenvectors of the correlation matrix .
With this correspondence at hand, computing Renyi correlations in subsystem becomes straightforward. Indeed the new (powered up) reduced density is
(65) |
Therefore, the new entanglement Hamiltonian is simply . We now wish to obtain the new correlation matrix
(66) |
with
(67) | |||||
(68) |
corresponding to the Renyi correlators
(69) |
To do so it is sufficient to reverse the logic leading to the determination of . From (61) and (64), we have
(70) |
so the new correlation matrix may be obtained by simply multiplying all single particle energies by . We obtain
(71) |
In terms of the initial correlation matrix , the previous relation reads
(72) |
The result is true for any eigenstate and also at finite temperature, the only requirement being a free fermion Hamiltonian. In fact, (72) even holds for quadratic boson Hamiltonians, provided the correlation matrix be redefined as
(73) |
instead of (56). This can be shown by repeating the steps described above with bosonic operators instead of fermions.
Appendix B Series analysis
Series are analyzed by the method of Differential Approximants Fisher and Au-Yang (1979); Hunter and Baker (1979). These are generalizations of the better known d-log Pade method and allow for a power-law singularity. The function of interest is represented by a solution to a first order ordinary differential equation:
(74) |
Here, , and are polynomials of order , and respectively that are determined by demanding that the series expansion for the function match to some order. It is easy to show that these differential equations have power-law singularities at critical points given by the roots of the polynomial with exponents related to the residue of and the value of at that critical point. The less important term allows for a background function. Since the critical point of Eq. (9) is known to be at , we further bias the analysis by demanding that
(75) |
thus ensuring that critical point is at .
The calculation proceeds by first choosing orders for the polynomials , and . We fix that just allows for a constant background term, and keep to at least . The sum of the orders () can not exceed the number of terms (more correctly number of terms minus one after the biasing) that one knows in the series expansion. In general, different choice of the order of the three polynomials , and will give different values for critical exponents.
Appendix C The Numerical Linked Cluster Expansion
The procedure of the Appendix A allows one to calculate the entanglement entropy for free bosons on clusters of arbitrary shape and size. Thus, it may be used as a “cluster solver” for a NLCE procedure, similar to previous studies where Lanczos or DMRG were used as cluster solvers for interacting models Kallin et al. (2013, 2014); Stoudenmire et al. (2014).
We refer the reader to the relevant literature Rigol et al. (2006, 2007a, 2007b); Tang et al. (2013). The NLCE method is based on the fact that an extensive property of a lattice model can be expressed as a sum over contributions from all distinct clusters which are embeddable in the lattice. Then, this property per site is,
(76) |
where is the “weight”, or the unique contribution to that comes from a cluster ^{2}^{2}2Note, unique embedding of the same cluster shape must be accounted for by the sum.. It is defined recursively,
(77) |
where is any subcluster of . As dictated by the inclusion-exclusion principle, only connected clusters give non-zero weights to the sum. The NLCE procedure builds up the value of starting from the smallest cluster and ending with some maximal cluster size. For this paper, we employ only rectangular clusters Kallin et al. (2013), which makes the calculation of cluster embeddings trivial; then, the computational bottleneck is shifted to the cluster solver.
For interacting models, the computational cost of calculating the given property on the cluster is extremely expensive. However, as described in the last section, a very efficient numerical method can be devised for calculating the Renyi entropies for a free bosons on the square lattice using Equations (42,43). This allows us to perform the NLCE to extremely high orders.
As with any cluster solver on the square lattice, the NLCE is able to isolate the entanglement due to a corner by subtracting off the contributions to from the linear portions of the boundary (the area-law term in Eq. (32)). As described in Ref. Kallin et al. (2013) the property then becomes the isolated contribution for a single corner, given by Eq. (33). This property is made extensive by adding the values for all possible translations of the corner around the cluster .
Finally, for studies of critical systems such as the free boson, the definition of a length scale is important. In this paper, we are interested in benchmarking the success of the NLCE in converging the exact results for in the limit where this length scale approaches infinity. A second advantage of using rectangular clusters is the easy definition of this length scale in relation to the cluster order . In past studies, several definitions have been used, including the arithmetic and geometric means of the rectangle linear dimensions Kallin et al. (2013, 2014); Stoudenmire et al. (2014); Kallin (2014). In the main text, we define the cluster length scale as the maximum of the two linear dimensions. This definition, which also defines the expansion order , is simple to use and is motivated by the intuition that the role of the smaller clusters in the NLCE sum is to cancel boundary effects. Thus long and narrow clusters of dimension should only be included at the same order as square clusters with ; otherwise these narrow clusters are not the boundary of any more isotropic cluster appearing to this order in the sum.
References
- Kallin et al. (2011) A. B. Kallin, M. B. Hastings, R. G. Melko, and R. R. P. Singh, Phys. Rev. B 84, 165134 (2011), ISSN 1098-0121.
- (2) M. A. Metlitski and T. Grover (????), eprint arXiv:1112.5166.
- Luitz et al. (2015) D. J. Luitz, X. Plat, F. Alet, and N. Laflorencie, Phys. Rev. B 91, 155145 (2015).
- Kulchytskyy et al. (2015) B. Kulchytskyy, C. M. Herdman, S. Inglis, and R. G. Melko (2015), eprint arXiv:1502.01722.
- Kitaev and Preskill (2006) A. Kitaev and J. Preskill, Phys. Rev. Lett. 96, 110404 (2006).
- Levin and Wen (2006) M. Levin and X.-G. Wen, Phys. Rev. Lett. 96, 110405 (2006).
- Isakov et al. (2011) S. V. Isakov, M. B. Hastings, and R. G. Melko, Nat. Phys. 7, 772 (2011).
- Wolf (2006) M. M. Wolf, Phys. Rev. Lett. 96, 010404 (2006).
- Gioev and Klich (2006) D. Gioev and I. Klich, Phys. Rev. Lett. 96, 100503 (2006).
- Metlitski et al. (2009) M. A. Metlitski, C. A. Fuertes, and S. Sachdev, Phys. Rev. B 80, 115122 (2009).
- Casini and Huerta (2012) H. Casini and M. Huerta, Phys. Rev. D 85, 125016 (2012).
- Grover (2014) T. Grover, Phys. Rev. Lett. 112, 151601 (2014).
- Bueno et al. (2015) P. Bueno, R. C. Myers, and W. Witczak-Krempa, Phys. Rev. Lett. 115, 021602 (2015).
- Calabrese and Cardy (2004) P. Calabrese and J. Cardy, Journal of Statistical Mechanics: Theory and Experiment 2004, P06002 (2004).
- Hastings et al. (2010) M. B. Hastings, I. González, A. B. Kallin, and R. G. Melko, Phys. Rev. Lett. 104, 157201 (2010).
- Melko et al. (2010) R. G. Melko, A. B. Kallin, and M. B. Hastings, Phys. Rev. B 82, 100409 (2010), ISSN 1098-0121.
- Humeniuk and Roscilde (2012) S. Humeniuk and T. Roscilde, Phys. Rev. B 86, 235116 (2012).
- Wang et al. (2013) L. Wang, D. Poilblanc, Z.-C. Gu, X.-G. Wen, and F. Verstraete, Phys. Rev. Lett. 111, 037202 (2013).
- Singh et al. (2011) R. R. P. Singh, M. B. Hastings, A. B. Kallin, and R. G. Melko, Phys. Rev. Lett. 106, 135701 (2011).
- Singh et al. (2012) R. R. P. Singh, R. G. Melko, and J. Oitmaa, Phys. Rev. B 86, 075106 (2012).
- Devakul and Singh (2014a) T. Devakul and R. R. P. Singh, Phys. Rev. B 90, 064424 (2014a).
- Devakul and Singh (2014b) T. Devakul and R. R. P. Singh, Phys. Rev. B 90, 054415 (2014b).
- Cardy (2011) J. Cardy, Phys. Rev. Lett. 106, 150404 (2011).
- Abanin and Demler (2012) D. A. Abanin and E. Demler, Phys. Rev. Lett. 109, 020504 (2012).
- Daley et al. (2012) A. J. Daley, H. Pichler, J. Schachenmayer, and P. Zoller, Phys. Rev. Lett. 109, 020505 (2012).
- Zamolodchikov (1986) A. Zamolodchikov, JETP Lett. 43, 730 (1986).
- Cardy and Calabrese (2010) J. Cardy and P. Calabrese, Journal of Statistical Mechanics: Theory and Experiment 2010, P04023 (2010).
- Jiang et al. (2013) H.-C. Jiang, R. R. P. Singh, and L. Balents, Phys. Rev. Lett. 111, 107205 (2013).
- Kallin et al. (2013) A. B. Kallin, K. Hyatt, R. R. P. Singh, and R. G. Melko, Phys. Rev. Lett. 110, 135702 (2013).
- Kallin et al. (2014) A. B. Kallin, E. M. Stoudenmire, P. Fendley, R. R. P. Singh, and R. G. Melko, Journal of Statistical Mechanics: Theory and Experiment 2014, P06009 (2014).
- Stoudenmire et al. (2014) E. M. Stoudenmire, P. Gustainis, R. Johal, S. Wessel, and R. G. Melko, Phys. Rev. B 90, 235106 (2014).
- Helmes and Wessel (2014) J. Helmes and S. Wessel, Phys. Rev. B 89, 245120 (2014).
- Wegner (1972) F. J. Wegner, Phys. Rev. B 5, 4529 (1972).
- Holzhey et al. (1994) C. Holzhey, F. Larsen, and F. Wilczek, Nucl. Phys. B 424, 443 (1994).
- Fagotti and Calabrese (2010) M. Fagotti and P. Calabrese, Journal of Statistical Mechanics: Theory and Experiment 2010, P04016 (2010).
- Fagotti and Calabrese (2011) M. Fagotti and P. Calabrese, Journal of Statistical Mechanics: Theory and Experiment 2011, P01017 (2011).
- Calabrese et al. (2010a) P. Calabrese, J. Cardy, and I. Peschel, Journal of Statistical Mechanics: Theory and Experiment 2010, P09003 (2010a).
- Calabrese and Essler (2010) P. Calabrese and F. H. L. Essler, Journal of Statistical Mechanics: Theory and Experiment 2010, P08029 (2010).
- Calabrese et al. (2010b) P. Calabrese, M. Campostrini, F. Essler, and B. Nienhuis, Phys. Rev. Lett. 104, 095701 (2010b).
- Cardy (1996) J. Cardy, Scaling and Renormalization in Statistical Physics (Cambridge University Press, 1996).
- Kos et al. (2014) F. Kos, D. Poland, and D. Simmons-Duffin, Journal of High Energy Physics 2014 (2014).
- Alba (2013) V. Alba, Journal of Statistical Mechanics: Theory and Experiment 2013, P05013 (2013).
- Coser et al. (2014) A. Coser, L. Tagliacozzo, and E. Tonni, Journal of Statistical Mechanics: Theory and Experiment 2014, P01008 (2014).
- Di Francesco et al. (1997) P. Di Francesco, P. Mathieu, and D. Sénéchal, Conformal Field Theory (Springer Verlag, 1997).
- (45) Note1, for example the result for an interval of size in a periodic system of size follows from Eq. (12) with the conformal mapping .
- Lieb et al. (1961) E. Lieb, T. Schultz, and D. Mattis, Ann. Phys. 16, 407 (1961).
- Chung and Peschel (2001) M.-C. Chung and I. Peschel, Phys. Rev. B 64, 064412 (2001).
- Peschel (2003) I. Peschel, J. Phys. A: Math. Gen. 36, L205 (2003).
- Oitmaa et al. (2006) J. Oitmaa, C. Hamer, and W. Zheng, Series Expansion Methods for strongly interacting lattice models (Cambridge University Press, 2006).
- Fisher and Au-Yang (1979) M. E. Fisher and H. Au-Yang, Journal of Physics A: Mathematical and General 12, 1677 (1979).
- Hunter and Baker (1979) D. L. Hunter and G. A. Baker, Phys. Rev. B 19, 3808 (1979).
- Affleck and Bonner (1990) I. Affleck and J. C. Bonner, Phys. Rev. B 42, 954 (1990).
- Casini and Huerta (2009) H. Casini and M. Huerta, Journal of Physics A: Mathematical and Theoretical 42, 504007 (2009).
- Casini and Huerta (2007) H. Casini and M. Huerta, Nuclear Physics B 764, 183 (2007).
- Fradkin and Moore (2006) E. Fradkin and J. E. Moore, Phys. Rev. Lett. 97, 050404 (2006).
- Zaletel et al. (2011) M. P. Zaletel, J. H. Bardarson, and J. E. Moore, Phys. Rev. Lett. 107, 020402 (2011).
- Stéphan et al. (2011) J.-M. Stéphan, G. Misguich, and V. Pasquier, Phys. Rev. B 84, 195128 (2011).
- Ju et al. (2012) H. Ju, P. Kallin, Ann B.Fendley, M. B. Hastings, and R. G. Melko, Phys. Rev. B 85, 165121 (2012).
- Stéphan et al. (2013) J.-M. Stéphan, H. Ju, P. Fendley, and R. G. Melko, New. J. Phys. 15, 015004 (2013).
- Inglis and Melko (2013) S. Inglis and R. G. Melko, New Journal of Physics 15, 073048 (2013).
- Chen et al. (2015) X. Chen, G. Y. Cho, T. Faulkner, and E. Fradkin, Journal of Statistical Mechanics: Theory and Experiment 2015, P02010 (2015).
- Schollwöck (2005) U. Schollwöck, Rev. Mod. Phys. 77, 259 (2005).
- (63) Calculations were performing using the ITensor Library: http://itensor.org/.
- Eisler and Peschel (2009) V. Eisler and I. Peschel, J. Phys. A: Math. Gen. 42, 504003 (2009).
- Bombelli et al. (1986) L. Bombelli, R. K. Koul, J. Lee, and R. D. Sorkin, Phys. Rev. D 34, 373 (1986).
- Srednicki (1993) M. Srednicki, Phys. Rev. Lett. 71, 666 (1993).
- Rigol et al. (2006) M. Rigol, T. Bryant, and R. R. P. Singh, Phys. Rev. Lett. 97, 187202 (2006).
- Rigol et al. (2007a) M. Rigol, T. Bryant, and R. R. P. Singh, Phys. Rev. E 75, 061118 (2007a).
- Rigol et al. (2007b) M. Rigol, T. Bryant, and R. R. P. Singh, Phys. Rev. E 75, 061119 (2007b).
- Tang et al. (2013) B. Tang, E. Khatami, and M. Rigol, Computer Physics Communications 184, 557 (2013).
- (71) Note2, note, unique embedding of the same cluster shape must be accounted for by the sum.
- Kallin (2014) A. B. Kallin, Ph.D. thesis, University of Waterloo (2014), URL https://uwspace.uwaterloo.ca/handle/10012/8539.