Renyi Entanglement Entropy of Interacting Fermions Calculated Using
Continuous-Time Quantum Monte Carlo Method
We present a new algorithm for calculating the Renyi entanglement entropy of interacting fermions using the continuous-time quantum Monte Carlo method. The algorithm only samples interaction correction of the entanglement entropy, which by design ensures efficient calculation of weakly interacting systems. Combined with Monte Carlo reweighting, the algorithm also performs well for systems with strong interactions. We demonstrate the potential of this method by studying the quantum entanglement signatures of the charge-density-wave transition of interacting fermions on a square lattice.
pacs:03.65.Ud, 02.70.Ss, 71.10.Fd
We have entered an information age of condensed matter physics. Many information theoretical tools have been used to identify exotic phases and phase transitions Feiguin et al. (2007); Bonesteel and Yang (2007); Isakov et al. (2011); Jiang et al. (2012); Depenbrock et al. (2012). These tools are especially useful when a conventional local order parameter characterization fails, e.g. in the case of topological order Wen (1990). The reduced density matrix plays a central role in this quantum information perspective. It contains all the information about a subregion A when viewing the remaining part of the system as environment. The entanglement entropy (EE) quantitatively measures the entanglement between the subregion A and its environment. The rank- Renyi EE, in particular, reads
Compared to the more familiar von Neumann entropy the Renyi EE is easier to evaluate both analytically and numerically. Quantum Monte Carlo methods have been used to calculate entanglement properties of many bosonic Buividovich and Polikarpov (2008); Hastings et al. (2010); Herdman et al. (2014) and fermionic Zhang et al. (2011); McMinis and Tubman (2013); Grover (2013); Broecker and Trebst (2014) systems. Some interesting applications include the calculation of the topological entanglement entropy Levin and Wen (2006); Kitaev and Preskill (2006); Isakov et al. (2011) and the identification of interacting topological insulators Assaad et al. (2014); Wang et al. (2014a).
Since the Monte Carlo method samples but not directly, it will suffer from a severe problem when the entanglement entropy is large Hastings et al. (2010); Zhang et al. (2011); Grover (2013). For example, the Renyi EE of a generic gapped many-body system follows an area law Eisert et al. (2010) and as increases linearly with the boundary size, the Monte Carlo simulation has to sample rare events whose probability vanishes exponentially. The problem is even more severe for highly entangled systems which violate the area law Wolf (2006); Gioev and Klich (2006); Swingle (2010); Ding et al. (2012); Lai et al. (2013); Leschke et al. (2014), including the systems with Fermi surface or at nonzero temperature. One way to alleviate the problem is to employ the ratio method Hastings et al. (2010); Humeniuk and Roscilde (2012), which splits the estimator of into products of intermediate ratios and samples each term individually. However, too many intermediate ratios will lead to an accumulation of errors in the final result Humeniuk and Roscilde (2012).
In this paper we introduce an algorithm for computing the Renyi EE of interacting fermions using the continuous-time quantum Monte Carlo method (CTQMC) Rubtsov et al. (2005); Gull et al. (2011). The main advantage of the algorithm is that it samples the interaction correction of the Renyi EE. By design the method allows to easily calculate large Renyi EE of weakly interacting systems. Our benchmark shows that it is also applicable to strongly interacting systems when combining with the ratio method Hastings et al. (2010); Humeniuk and Roscilde (2012). As an application we calculate the Renyi EE of spinless fermions across a charge-density-wave transition and show it correctly predicts the critical point.
where is the partition function at inverse temperature and is the partition function defined on a -sheeted Riemann surface Melko et al. (2010). Ref.Broecker and Trebst (2014) introduced an artificial system with imaginary-time dependent Hamiltonian whose partition function at inverse temperature equals to . We here formulate an efficient algorithm to calculate Eq. (2) using CTQMC.
The first and crucial step of our method is to split Eq. (2) into , in which
and are the noninteracting counterparts of and 111In general one can choose and as partition functions of any quadratic Hamiltonian. It amounts to redefine the reference Hamiltonian around which one performs the interaction expansion. A particular useful example is to expand around a mean-field Hamiltonian and the resulting will be the entanglement entropy difference to a mean-field state. Employing this trick will further improve efficiency of our algorithm, especially in the long-range-ordered phase. We thank Fakher Assaad for pointing this out to us.. is the Renyi EE of noninteracting fermions and can be calculated easily using the correlation matrix method Peschel (2003); SM (). Eq. (4) describes the interaction corrections to the Renyi EE which is presumably much smaller than for weakly interacting systems. We use Monte Carlo (MC) sampling to calculate the quantity in the square bracket of Eq. (4), where we have introduced a free parameter to further control the MC dynamics. For an optimal choice of the sampled quantity in the square bracket is close to one and will contribute mostly to . The fact that MC sampling corrects an educated guess is another appealing feature of the present algorithm. In practice, can be determined from a rough estimate of (either based on existing theory of EE scaling laws or in the MC equilibration steps). MC sampling will correct the estimate and restore the exact result no matter what the initial choice of was.
Our method is general and applicable to any fermionic system which is accessible to Monte Carlo simulations. For illustration purposes we here focus on calculating the rank- Renyi EE of an interacting spinless fermions model
We can treat Eq. (6) and Eq. (7) on equal footing and use the same configuration for each term in the sampling. Fig. 1(a) shows an example configuration with vertices, where the imaginary times satisfy . Any of the configurations is a valid configuration in both ensembles, but with different weights
where and are and matrices whose matrix elements only depend on the noninteracting Green’s functions SM (). On a bipartite lattice with repulsive interaction the weights are positive Huffman and Chandrasekharan (2014); Wang et al. (2014b) and there is no sign problem in the Monte Carlo simulation.
We introduce an ensemble flag and perform MC simulation in an extended ensemble Burovski et al. (2006); Humeniuk and Roscilde (2012) with the partition function . Now a MC configuration corresponds to a given set of variables: the ensemble flag , the perturbation order and the vertex configurations. Two kinds of MC updates are necessary to ensure ergodicity of the sampling, shown in Fig.1(b). First, we keep the ensemble flag unchanged and update the vertex configuration by either adding or removing one vertex. This is done either by proposing a candidate vertex at a random time (in the interval ) and a random bond (out of possible ones) or randomly choosing an existing vertex (out of possible ones) to be removed. The acceptance ratios are
The second class of MC updates switch the ensemble to while keeping the configuration fixed. The acceptance probability is
A direct evaluation of is numerically expensive and unstable, since it requires the calculation of two determinants Eq. (8) and Eq. (9) and their ratio. We utilize the fact that the determinant ratio of the zeroth expansion order is one and keep updating the ratio Eq. (12) during the MC updates. Therefore the ensemble switch can then be implemented without any matrix operations and is very cheap.
The fact that MC simulation only samples the interaction corrections of the Renyi EE is the major advantage of the present approach. On top of that, the parameter provides additional control over the MC dynamics. An ideal choice of will balance the probability of each ensemble and the observable in the square bracket of Eq. (13) will be of order one.
Fig. 2 shows benchmark results of Renyi EE of an -site chain with open boundaries. Our CTQMC results perfectly reproduce the exact diagonalization results (black solid lines) for all interaction strengths and temperatures. In the strong coupling limit the system has two fold quasi-degenerated ground states corresponding to the two staggered charge-density-wave (CDW) configurations. They are separated from other states by an excitation gap . The Renyi EE indeed approaches to when the excitation gap is much larger than the temperature.
Although the above algorithm can already handle a large class of interesting problems, it will suffer from very low transition probability when the two ensembles and have vanishing small overlap in configurational space (for example if their average perturbation order differs substantially). To achieve better performance, we additionally employ the ratio method Hastings et al. (2010); Humeniuk and Roscilde (2012) and split the partition function ratio into products of a series of intermediate ratios, each one corresponds to two ensembles where the size of region A only differs by a small amount of sites 222In the case the weight Eq. (9) is replaced by Eq. (8) with a different subregion A. Because can be much smaller than the total entanglement entropy , we need far less intermediate ratios compare to the other approach Broecker and Trebst (2014) to achieve accurate results.
Our method is efficient for simulating highly entangled systems because it samples the difference of entanglement entropy to that of the free fermions. For obvious application like weakly interacting fermions, the algorithm enjoys additional advantage of dealing with small average expansion order, which is proportional to the interaction energy Rubtsov et al. (2005); Gull et al. (2011).
We will demonstrate the power of our method with a more challenging example, where interaction turns the free fermions state into another phase which follows different entanglement entropy scaling. We consider the model (5) on a square lattice at half-filling Scalapino et al. (1984); Gubernatis et al. (1985). At zero temperate it has CDW ground state for arbitrary weak repulsive interaction because of Fermi surface nesting. The CDW ground state in the strong coupling limit can be interpreted as an “antiferromagnetic Ising” ground state of a classical lattice gas. Upon increasing the temperature, the CDW state undergoes an Ising phase transition. The transition temperature is exponentially small in the weak coupling limit and is proportional to the interaction strength in the strong coupling limit. Figure 20 of Ref.Gubernatis et al. (1985) reports the phase diagram of this model.
We here revisit this problem from a quantum information perspective by calculating the Renyi EE across the phase transition. We consider a cylindrical subregion A embedded in a torus, Fig. 3(a) inset. The boundary length is chosen to be and is independent of the subregion size . Figure 3 shows the rank- Renyi EE as a function at various temperature and interaction strengths, which clearly exhibits two different characteristic behaviors. At high temperature or weak interaction the Renyi EE increases linearly with , indicating a highly entangled state with a volume scaling law for the EE. At low temperature or for strong interaction the Renyi EE is strongly suppressed and is essentially flat as increases, which is a characteristic behavior of a gapful CDW state. This state is approximately a linear superposition of two simple product states. Figure 3(a) also suggests that simulations deep inside the CDW state may not be ideal for the present method because the EE is much smaller than the free fermion state. The fact that we can nevertheless easily obtain correct results in this region provides a stringent test of the proposed method.
which cancels the bulk contribution in and exhibits boundary law even at nonzero temperature. Studies of classical and quantum spin models Melko et al. (2010); Singh et al. (2011); Iaconis et al. (2013); Storms and Singh (2014) show that the scaled mutual information crosses around the transition point. We consider a region A with and the mutual information can be directly calculated using two data points of in the Fig.3.
Figure 4 shows the scaled mutual information versus interaction strength and temperature for several system sizes. The crossing points of provide an estimation of transition point at or at , which is consistent with the phase diagram reported in Gubernatis et al. (1985). Similar to the case of Ref. Melko et al. (2010); Iaconis et al. (2013), the curves of mutual information also cross around (not shown). These results indicate that our approach can be a versatile tool to calculate the Renyi EE and detect phase transitions in the interacting fermionic models.
The proposed method allows efficient calculation of the Renyi entanglement entropy of interacting fermions by reallocating computational resources to interaction corrections. Compare to the direct approaches Grover (2013); Broecker and Trebst (2014), it can avoid sampling exponentially rare events caused by fast increase of free fermion entanglement entropy. The method is applicable to any fermionic system which could be simulated with CTQMC method and may provide further insights to the unconventional quantum critical point Assaad and Herbut (2013); Wang et al. (2014b) and the topological phase transitions Hohenadler et al. (2011); Wang et al. (2014c, a). Our method is an ideal tool to offer an entanglement perspective on the Kondo problem Sørensen et al. (2007a, b); Eriksson and Johannesson (2011); Bayat et al. (2012), which is intimately related to the topological entanglement entropy Levin and Wen (2006); Kitaev and Preskill (2006) of two dimensional gapful states Fendley et al. (2007). Most interestingly, the ability to calculate entanglement entropy in CTQMC offers a portal to study entanglement in the framework of dynamical-mean-field-theory Georges et al. (1996); Gull et al. (2011) and will shed light on entanglement properties of realistic correlated materials.
We thank Nikolay Prokof’ev for insightful suggestions. L.W. acknowledges KITPC Beijing for hospitality during the workshop “Precision Many-body Physics of Strongly Correlated Quantum Matter”, where part of this work was done. Simulations were performed on the Mönch cluster of Platform for Advanced Scientific Computing (PASC), the Brutus cluster at ETH Zurich and the “Monte Rosa” Cray XE6 at the Swiss National Supercomputing Centre (CSCS). We have used ALPS libraries Bauer et al. (2011) for Monte Carlo simulations and data analysis. This work was supported by ERC Advanced Grant SIMCOFE.
- Feiguin et al. (2007) A. Feiguin, S. Trebst, A. Ludwig, M. Troyer, A. Kitaev, Z. Wang, and M. Freedman, Phys. Rev. Lett. 98, 160409 (2007).
- Bonesteel and Yang (2007) N. E. Bonesteel and K. Yang, Phys. Rev. Lett. 99, 140405 (2007).
- Isakov et al. (2011) S. V. Isakov, M. B. Hastings, and R. G. Melko, Nature Physics 7, 772 (2011).
- Jiang et al. (2012) H.-C. Jiang, Z. Wang, and L. Balents, Nature Physics 8, 902 (2012).
- Depenbrock et al. (2012) S. Depenbrock, I. P. McCulloch, and U. Schollwöck, Phys. Rev. Lett. 109, 067201 (2012).
- Wen (1990) X.-G. Wen, Int. J. Mod. Phys. B 4, 239 (1990).
- Buividovich and Polikarpov (2008) P. V. Buividovich and M. I. Polikarpov, Nuclear Physics B 802, 458 (2008).
- Hastings et al. (2010) M. B. Hastings, I. González, A. B. Kallin, and R. G. Melko, Phys. Rev. Lett. 104, 157201 (2010).
- Herdman et al. (2014) C. M. Herdman, P. N. Roy, R. G. Melko, and A. Del Maestro, Phys. Rev. B 89, 140501 (2014).
- Zhang et al. (2011) Y. Zhang, T. Grover, and A. Vishwanath, Phys. Rev. Lett. 107, 067202 (2011).
- McMinis and Tubman (2013) J. McMinis and N. M. Tubman, Phys. Rev. B 87, 081108 (2013).
- Grover (2013) T. Grover, Phys. Rev. Lett. 111, 130402 (2013).
- Broecker and Trebst (2014) P. Broecker and S. Trebst, arXiv:1404.3027 (2014).
- Levin and Wen (2006) M. Levin and X.-G. Wen, Phys. Rev. Lett. 96, 110405 (2006).
- Kitaev and Preskill (2006) A. Kitaev and J. Preskill, Phys. Rev. Lett. 96, 110404 (2006).
- Assaad et al. (2014) F. F. Assaad, T. C. Lang, and F. Parisen Toldin, Phys. Rev. B 89, 125121 (2014).
- Wang et al. (2014a) D. Wang, S. Xu, Y. Wang, and C. Wu, arXiv:1405.2043 (2014a).
- Eisert et al. (2010) J. Eisert, M. Cramer, and M. B. Plenio, Rev. Mod. Phys. 82, 277 (2010).
- 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).
- Swingle (2010) B. Swingle, Phys. Rev. Lett. 105, 050502 (2010).
- Ding et al. (2012) W. Ding, A. Seidel, and K. Yang, Phys. Rev. X 2, 011012 (2012).
- Lai et al. (2013) H.-H. Lai, K. Yang, and N. E. Bonesteel, Phys. Rev. Lett. 111, 210402 (2013).
- Leschke et al. (2014) H. Leschke, A. V. Sobolev, and W. Spitzer, Phys. Rev. Lett. 112, 160403 (2014).
- Humeniuk and Roscilde (2012) S. Humeniuk and T. Roscilde, Phys. Rev. B 86, 235116 (2012).
- Rubtsov et al. (2005) A. Rubtsov, V. Savkin, and A. Lichtenstein, Phys. Rev. B 72, 035122 (2005).
- Gull et al. (2011) E. Gull, A. J. Millis, A. I. Lichtenstein, A. N. Rubtsov, M. Troyer, and P. Werner, Rev. Mod. Phys. 83, 349 (2011).
- Calabrese and Cardy (2004) P. Calabrese and J. Cardy, J. Stat. Mech.: Theor. Exp. 2004, P06002 (2004).
- Melko et al. (2010) R. G. Melko, A. B. Kallin, and M. B. Hastings, Phys. Rev. B 82, 100409 (2010).
- Peschel (2003) I. Peschel, J. Phys. A: Math. Gen. 36, L205 (2003).
- (31) See Supplemental Material for details of calculating noninteracting Renyi entanglement entropy and noninteracting Green’s functions .
- Wang et al. (2014b) L. Wang, P. Corboz, and M. Troyer, arXiv:1407.0029 (2014b).
- Huffman and Chandrasekharan (2014) E. F. Huffman and S. Chandrasekharan, Phys. Rev. B 89, 111101 (2014).
- Burovski et al. (2006) E. Burovski, N. Prokof’ev, B. Svistunov, and M. Troyer, New J. Phys. 8, 153 (2006).
- Scalapino et al. (1984) D. J. Scalapino, R. L. Sugar, and W. D. Toussaint, Phys. Rev. B 29, 5253 (1984).
- Gubernatis et al. (1985) J. E. Gubernatis, D. J. Scalapino, R. L. Sugar, and W. D. Toussaint, Phys. Rev. B 32, 103 (1985).
- Wolf et al. (2008) M. Wolf, F. Verstraete, M. Hastings, and J. Cirac, Phys. Rev. Lett. 100, 070502 (2008).
- Singh et al. (2011) R. R. P. Singh, M. B. Hastings, A. B. Kallin, and R. G. Melko, Phys. Rev. Lett. 106, 135701 (2011).
- Iaconis et al. (2013) J. Iaconis, S. Inglis, A. B. Kallin, and R. G. Melko, Phys. Rev. B 87, 195134 (2013).
- Storms and Singh (2014) M. Storms and R. R. P. Singh, Phys. Rev. E 89, 012125 (2014).
- Assaad and Herbut (2013) F. F. Assaad and I. F. Herbut, Phys. Rev. X 3, 031010 (2013).
- Hohenadler et al. (2011) M. Hohenadler, T. C. Lang, and F. F. Assaad, Phys. Rev. Lett. 106, 100403 (2011).
- Wang et al. (2014c) L. Wang, H.-H. Hung, and M. Troyer, arXiv:1402.6704 (2014c).
- Sørensen et al. (2007a) E. S. Sørensen, M.-S. Chang, N. Laflorencie, and I. Affleck, J. Stat. Mech.: Theor. Exp. 2007, L01001 (2007a).
- Sørensen et al. (2007b) E. S. Sørensen, M.-S. Chang, N. Laflorencie, and I. Affleck, J. Stat. Mech.: Theor. Exp. 2007, P08003 (2007b).
- Eriksson and Johannesson (2011) E. Eriksson and H. Johannesson, Phys. Rev. B 84, 041107 (2011).
- Bayat et al. (2012) A. Bayat, S. Bose, P. Sodano, and H. Johannesson, Phys. Rev. Lett. 109, 066403 (2012).
- Fendley et al. (2007) P. Fendley, M. P. A. Fisher, and C. Nayak, J Stat Phys 126, 1111 (2007).
- Georges et al. (1996) A. Georges, G. Kotliar, W. Krauth, and M. J. Rozenberg, Rev. Mod. Phys. 68, 13 (1996).
- Bauer et al. (2011) B. Bauer, L. D. Carr, H. G. Evertz, A. Feiguin, J. Freire, S. Fuchs, L. Gamper, J. Gukelberger, E. Gull, S. Guertler, A. Hehn, R. Igarashi, S. V. Isakov, D. Koop, P. N. Ma, P. Mates, H. Matsuo, O. Parcollet, G. Pawlowski, J. D. Picon, L. Pollet, E. Santos, V. W. Scarola, U. Schollwock, C. Silva, B. Surer, S. Todo, S. Trebst, M. Troyer, M. L. Wall, P. Werner, and S. Wessel, J. Stat. Mech.: Theor. Exp. 2011, P05001 (2011).
- Assaad and Evertz (2008) F. F. Assaad and H. G. Evertz, Computational Many Particle Physics, Lecture Notes in Physics 739, 277 (2008).
- Hirsch (1988) J. Hirsch, Phys. Rev., B Condens. Matter 38, 12023 (1988).
Appendix A Supplementary Materials
a.1 Renyi Entanglement Entropy of Noninteracting Fermions
The reduced density matrix of a free fermion system is fully determined by the correlation matrix Peschel (2003)
where denotes thermal average with respect the noninteracting Hamiltonian. To calculate the Renyi entanglement entropy, we restrict to a subregion A and diagonalize it to get eigenvalues . Introducing , the rank- Renyi entanglement entropy of free fermions is calculated as
a.2 Noninteracting Green’s Functions
where , are the vertex indices and are the site indices. The term arises from the constant shift defined in the interaction term. For and , the corresponding site index is interpreted as a site in the artificially enlarged part of the system Broecker and Trebst (2014). and are the noninteracting Green’s functions for the and ensembles. They have the following general form Assaad and Evertz (2008)
where is the noninteracting time-ordered propagator, with the matrix being the quadratic part of the Hamiltonian.
For one has and the fact that is imaginary-time independent allows simplification of Eq. (16),
Since only depends the imaginary-time difference we precompute it on a fine -mesh and use linear interpolations to get the value of for arbitrary imaginary-times.
For one has and is an artificial imaginary-time dependent Hamiltonian Broecker and Trebst (2014). In general depends both on and . We can also precompute on a fine mesh and perform bilinear interpolations. However, storing is not possible for large system size, and as a compromise we only pre-calculate on a coarse imaginary-time grid using Hirsch’s matrix inversion method Hirsch (1988). In the CTQMC calculation, whenever a new vertex is proposed, we find the nearest imaginary time point such that and and compute the required Green’s function using
In this way we avoid CPU expensive matrix inversion in Eq. (16).