Renyi Entanglement Entropy of Interacting Fermions Calculated Using Continuous-Time Quantum Monte Carlo Method

Renyi Entanglement Entropy of Interacting Fermions Calculated Using
Continuous-Time Quantum Monte Carlo Method

Lei Wang and Matthias Troyer Theoretische Physik, ETH Zurich, 8093 Zurich, Switzerland

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.

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.

In the replica approach Calabrese and Cardy (2004); Melko et al. (2010), the rank- Renyi entanglement entropy is calculated as


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.

Eq. (12)

Eqs. (10-11)

Eq. (12)

Eqs. (10-11)

Figure 1: Key concepts of the algorithm. (a) A configuration with vertices in the time interval and vertex in . The weights of this configuration are given in Eqs. (8-9). (b) The extended configuration space combines two ensembles and . MC updates Eqs. (10-11) change the vertex configuration and Eq. (12) switch between the two ensembles.

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


In the CTQMC method, the partition function ratios are expanded in terms of the interaction vertices Rubtsov et al. (2005); Gull et al. (2011); Wang et al. (2014b),


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 relative time spent in each ensemble Humeniuk and Roscilde (2012) provides an estimator of the ratio between Eq. (6) and Eq. (7), thus


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.

Figure 2: Renyi EE of an -site open chain with equal partitions. The black solid lines are the exact diagonalization results, which agree perfectly with the results calculated by our algorithm (circles).

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.

Figure 3: The rank-2 Renyi entanglement entropy versus the subregion size of a cylinder embedded in a torus for (a) and (b) .

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.

Figure 4: The scaled mutual information for various system sizes at (a) and (b) . The crossing point agrees with transition point determined in Ref. Gubernatis et al. (1985) using correlation functions.

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.

We next proceed to a quantitative determination of the phase boundaries. We consider the mutual information Wolf et al. (2008); Melko et al. (2010)


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.


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

The Green’s function matrices appeared in Eqs. (8-9) reads


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).

Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description