# Entanglement, noise, and the cumulant expansion

###### Abstract

We put forward a simpler and improved variation of a recently proposed method to overcome the signal-to-noise problem found in Monte Carlo calculations of the entanglement entropy of interacting fermions. The present method takes advantage of the approximate lognormal distributions that characterize the signal-to-noise properties of other approaches. In addition, we show that a simple rewriting of the formalism allows circumvention of the inversion of the restricted one-body density matrix in the calculation of the -th Rényi entanglement entropy for . We test our technique by implementing it in combination with the hybrid Monte Carlo algorithm and calculating the Rényi entropies of the 1D attractive Hubbard model. We use that data to extrapolate to the von Neumann () and cases.

###### pacs:

03.65.Ud, 05.30.Fk, 03.67.Mn## I Introduction

Recently HMCEE (), we proposed an algorithm to compute the Rényi entanglement entropy of interacting fermions. Many algorithms have been proposed to this effect in the last few years Buividovich (); Melko (); Humeniuk (); McMinis (); Broecker (); WangTroyer (); Luitz (). Our proposal, based on the free-fermion decomposition approach of Ref. Grover (), overcomes the signal-to-noise problem present in that approach and is compatible with the hybrid Monte Carlo (HMC) method HMC () widely used in the context of lattice quantum chromodynamics. The core idea of our method is that, by differentiating with respect to an auxiliary parameter , one may carry out a Monte Carlo (MC) calculation of with a probability measure that includes entanglement properties explicitly. [This was not the case in the approach of Ref. Grover (), where the probability measure factored across auxiliary field replicas; we identified this as the cause of the signal-to-noise problem (see below)]. Once the MC calculation is done, integration with respect to returns the desired entanglement entropy relative to that of a noninteracting system (which is easily computed separately).

In this work, we describe and implement a variation on that Monte Carlo algorithm which, while sharing the properties and core idea mentioned above, differs from it in two important ways; the new method, in fact, is different enough that we advocate its use over our original proposal. First, the new method takes advantage of the approximate lognormal shape of the underlying statistical distributions of the fermion determinants, which we already noted in Ref. HMCEE () and which we explain in detail below. Second, and more importantly, the present method is simpler than our original proposal: whereas in the latter the parameter multiplied the coupling constant (thus generating a rather involved set of terms upon differentiation of the fermion determinant), here is coupled to the number of fermion species . As we show below, this choice not only simplifies the implementation, but also exposes the central role of the logarithm of the fermion determinant in our calculation of , and thus brings to bear the approximate lognormality property mentioned above.

Below, we present the basic formalism, review the evidence for approximate lognormal distributions, and explain our method. Besides the points mentioned above, in our calculations we have found the present method to be more numerically stable than its predecessor. We explain this in detail in our Results section.

In addition to the new method, we show that it is possible to rewrite part of the formalism in order to bypass the calculation of inverses of the restricted density matrix (see e.g. Broecker (); WangTroyer (); HMCEE ()) in the determination of Rényi entropies of order . To test our method, we computed the Rényi entropy of the 1D attractive Hubbard model using the previous as well as the new formalism, and checked that we obtained identical results. Going beyond the case, we present results for the Rényi entropies and find that higher Rényi entropies display lower statistical uncertainty in MC calculations.

## Ii Basic formalism

As in our previous work, we set the stage by briefly presenting the formalism of Ref. Grover (). The -th Rényi entropy of a sub-system of a given system is

(1) |

where is the reduced density matrix of sub-system . For a system with density matrix , the reduced density matrix is defined via a partial trace over the Hilbert space corresponding to the complement of our sub-system:

(2) |

An auxiliary-field path-integral form for , from which can be computed using MC methods for a wide variety of systems, was presented in Ref. Grover (), which we briefly review next.

As is well known from conventional many-body formalism, the full density matrix can be written as a path integral by means of a Hubbard-Stratonovich auxiliary-field transformation:

(3) |

for some normalized probability measure determined by the details of the underlying Hamiltonian (for more detail, see below and also Ref. MCReviews ()). Here, is the partition function, and is the density matrix of noninteracting particles in the external auxiliary field . One of the main contributions of Ref. Grover () was to show that the above decomposition determines not only the full density matrix but also the restricted one. Indeed, Ref. Grover () shows that

(4) |

where is the same probability used in Eq. (3),

(5) |

and

(6) |

Here, is the restricted Green’s function of the noninteracting system in the external field (see below), and , are the fermion creation and annihilation operators. The sums in the exponent of Eq. (5) go over those points in the system that belong to the subsystem .

Using the above formalism for the case of -component fermions, the entanglement entropy (c.f. Eq. 1) takes the form

(7) |

where the field integration measure, given by

(8) |

is over the “replicas” of the Hubbard-Stratonovich field (which result from taking the -th power of the path integral representation of shown above), and the normalization

(9) |

was included in the measure. It is worth noting that, by separating a factor of in the denominator of Eq. (7), an explicit form can be identified in the numerator as in the replica trick CalabreseCardy (), which corresponds to a partition function for copies of the system, “glued” together in the region .

The naive probability measure, namely

(10) |

factorizes across replicas, which makes it insensitive to entanglement. This factorization is the main reason why using as a MC probability leads to signal-to-noise issues (see Ref. Grover ()). In Eq. (10), encodes the dynamics of the -th fermion component, including the kinetic energy and the form of the interaction after a Hubbard-Stratonovich transformation. That matrix also encodes the form of the trial state in ground-state approaches (see e.g. Ref. MCReviews ()), which we use here; we have taken to be a Slater determinant. In finite-temperature approaches, is obtained by evolving a complete set of single-particle states in imaginary time.

The quantity that contains the pivotal contributions to entanglement is

(11) |

which we refer to below as the “entanglement determinant,” and where

(12) | |||||

The product played the role of an observable in Ref. Grover (), which is a natural interpretation given Eq. (7). However, we will interpret this differently below. Other than the field replicas, the new ingredient in the determination of is the restricted Green’s function . This is the same as the noninteracting one-body density matrix of the -th fermion component in the background field , but the arguments are restricted to the region (see Ref. Grover () and also Ref. Peschel (), where expressions were originally derived for the reduced density matrix of noninteracting systems, based on reduced Green’s functions).

## Iii Avoiding inversion of the reduced Green’s function for

As noted in Ref. Assaad (), for , no inversion of is actually required in the calculation of the entanglement determinant , as the equations clearly simplify in that case. However, for higher it is not obvious how to avoid such an inversion. Here, however, we show that this calculation can indeed be accomplished without inversion. We begin by noting that

(13) |

where is a block diagonal matrix (one block per replica ):

(14) |

and

(15) |

with

(16) |

The equivalence of the determinants in Eq. (13) can be shown in a straightforward fashion: the factor is easily understood, as that matrix is block diagonal and therefore its determinant reproduces the first r.h.s. factor in the first line of Eq. (12); the remaining factor relies on the identity

(17) |

which is valid for arbitrary block matrices , is a standard result often used in many-body physics (especially when implementing a Hubbard-Stratonovich transformation), and can be shown using so-called elementary operations on rows and columns.

Within the determinant of Eq. (13), we may of course multiply and :

(18) |

where is a block diagonal matrix defined by

(19) |

and

(20) |

Equation (18) shows our claim, as we may use in our calculations instead of , and the former contains no inverses of .

Summarizing, a class of approaches to calculating for , based on the Hubbard-Stratonovich representation of (also known as free-fermion decomposition), requires computing , which in turn requires inverting per Eq. (12). By arriving at Eq. (18), and given that

(21) |

[Eq. (13) and beyond] we have shown that no inversions are actually required, as contains no inverses. While this is a desirable feature from a numerical point of view, it should be mentioned that, from a computational-cost point of view, the price of not inverting reappears in the fact that , though sparse, scales linearly with in size.

For the remaining of this work, calculations carried out at use the approach, which is based on Eq. (12) and the ‘proposed method’ described below. We reproduced those results by switching to the approach, which uses Eq. (18) (as well as the method described below), and then proceeded to higher with the latter.

## Iv A statistical observation: lognormal distribution of the entanglement determinant

In Ref. HMCEE (), we presented examples of the approximate log-normal distributions obeyed by when sampled according to . One such example is reproduced here for reference in Fig. 1.

The fact that such distributions are approximately log-normal, at least visually, suggests that one may use the cumulant expansion to determine . Indeed, in general,

(22) | |||||

where is the -th cumulant of , and the first two nonzero cumulants are given by

(23) |

and

(24) |

for a functional , and where the expectation value here taken with respect to the produce measure . If the distribution of were truly gaussian, the above series would terminate after the first two terms, which would provide us with an efficient way to bypass signal-to-noise issues in the determination of with stochastic methods NoiseSignProblemStatistics (). Unfortunately, the distribution is not exactly gaussian. Moreover, the cumulants beyond are often extremely sensitive to the details of the distribution (i.e. they can fluctuate wildly), they are hard to determine stochastically (the signal-to-noise problem re-emerges), and there is no easy way (that we know of) to obtain analytic insight into the large- behavior of . However, this approximate log-normality does provide a path forward, as it indicates that we may still evaluate with good precision with MC methods. As we will see in the next sections, this is enough to determine if we are willing to pay the price of a one-dimensional integration on a compact domain.

Although (approximate) lognormality in the entanglement determinant seems very difficult to prove analytically in the present case, evidence of its appearance has been found in systems as different as ultracold atoms and relativistic gauge theories NoiseSignProblemStatistics (); LogNormalDeGrand (). The underlying reason for this distribution appears to be connected to a similarity between the motion of electrons in disordered media and lattice fermions in the external auxiliary (gauge) field in MC calculations.

## V Proposed method

Starting from the right-hand side of Eq. (7), we introduce an auxiliary parameter and define a function via

(25) |

At ,

(26) |

while for , yields the entanglement entropy:

(27) |

Using Eq. (25),

(28) |

where

(29) |

In the presence of an even number of flavors and attractive interactions, and are real and non-negative for all , such that there is no sign problem and above is a well-defined, normalized probability measure.

As in our previously proposed method, we can then calculate by taking the point as a reference and computing using

(30) |

where

(31) |

We thus obtain an integral form of the interacting Rényi entropy that can be computed using any MC method (see e.g. MCReviews ()), in particular HMC HMC ().

As in our previous work, we note that the above expectation values are determined with respect to the probability measure , which communicates correlations responsible for entanglement. In contrast to the canonical MC probability , which corresponds to statistically independent copies of the Hubbard-Stratonovich field, this admittedly more complicated distribution does not exhibit the factorization to blame for the signal-to-noise problems present in the approach as originally formulated.

Using Eq. (30) requires Monte Carlo methods to evaluate as a function of , followed by integration over . As in our previous method, we find here that is a smooth function of , which is essentially linear in the present case. It is therefore sufficient to perform the numerical integration using a uniform grid. The stochastic evaluation of , for fixed subregion , can be expected to feature roughly symmetric fluctuations about the mean. As a consequence, the statistical effects on the entropy are reduced after integrating over .

Finally, we note an interesting application of Jensen’s inequality at . At that point

which must be satisfied by our calculations. Our Monte Carlo results at indeed satisfy this bound.

## Vi Results

### vi.1 Second Rényi entropy

As a first test of our algorithm and in efforts to make contact with previous work HMCEE (); Grover (), we begin by showing results for the second Rényi entropy for the one-dimensional Hubbard chain with periodic boundary conditions at half filling, whose Hamiltonian is

(33) |

where the first sum includes and pairs of adjacent sites. We implemented a symmetric Trotter-Suzuki decomposition of the Boltzmann weight, with an imaginary-time discretization of (in lattice units). As mentioned earlier, the many-body factor in the Trotter-Suzuki approximation was treated by introducing a replica auxiliary field for each power of the reduced density matrix. As in our previous work, we implemented a Hubbard-Stratonovich transformation of a compact continuous form MCReviews ().

We present plots for with in Fig. 2. In contrast to the results obtained in Ref. HMCEE () and as mentioned above, the resulting expectation demonstrates surprisingly little curvature as the region size is varied and is stunningly linear as a function of the auxiliary parameter . Even after twice doubling the strength of the interaction, the curvature of constant-subsystem-size slices is increased only marginally. We note that if one assumes such benign curvature is a somewhat universal feature, at least for weakly-coupled systems, our method provides a means by which to rapidly estimate the entanglement entropy for a large portion of parameter space at the very least yielding a qualitative picture of its behavior as a function of the physically relevant input parameters.

We observe that this surface displays almost no torsion, its dominant features being those present in the noninteracting case i.e. an alternating shell-like structure. Toward larger region sizes, we observe a combination of twisting and translation culminating in the required, and somewhat delicate, cancellation upon reaching the full system size. Presented with this relatively forgiving geometry, we performed the required integration via cubic-spline interpolation. Using a uniformly spaced lattice of size points, we determine the desired entropy to a precision limited by statistical rather than systematic considerations.

### vi.2 Comparison to exact diagonalization

Shown in Fig. 3 are results for a system of size with a number of sites . For couplings and and region sizes , we find solid agreement with previous calculations in Refs. HMCEE (); Grover (), and as in the former, we observe convergence rather quickly with only decorrelated samples as can be seen in Fig. 4. Further, for large sample sizes , we observe that the standard error in the entropy , computed from the envelope defined by the MC uncertainty in the source for at each value in -space, scales asymptotically as up to minute corrections.

### vi.3 Results for

In this section, we extend the results presented above to . In order to highlight the differences between and , we show in Fig. 5 the Rényi entropies for (top to bottom) of the 1D attractive Hubbard model, as obtained with our method and the reformulation of the fermion determinant shown in Eq. (18).

As evident from the figure, increasing leads to lower values of at fixed subsystem size consistent with knowledge that the Rényi entropy is a nonincreasing function of its order. However, increasing also amplifies the fluctuations as a function of . Interestingly, the approach of our system to the large- regime is quite rapid, and after only the first few orders, the difference between consecutive entropies is only marginal, most obviously so at weak coupling. We also observe that, as is increased, the statistical fluctuations that define the error bars appear to be progressively more suppressed, which is particularly evident for the strongest coupling we studied, namely .

At the level of the auxiliary function , we again see very predictable changes in the geometry of this surface as a function both arguments as shown in Fig. 6. With fixed coupling and particle content, increasing the Rényi order results in a tilting effect reminiscent of that seen previously with increasing coupling, but rather than being localized away from vanishing subsystem size, the change is much more global, affecting all subsystems in a qualitatively similar fashion and leaving each surface’s characteristic quasi-linearity in intact. Although the shell-like structure present in this function’s dependence is amplified, this increased fluctuation affects the quality of the results negligibly at most, as again, the geometry remains amenable to fairly naive quadratures.

With the data presented above, we would be remiss if we did not attempt an extrapolation not only to the limit of infinite Rényi order , but also to the von Neumann entropy, despite knowledge of the formidable challenges presented by these extrapolations, particularly in the case of the latter. The former limit provides a lower bound on all finite-order entropies, whereas the latter is of interest to a variety of disciplines and has proven difficult to study. At fixed coupling and with the knowledge that the Rényi entropy is nonincreasing in the order, we found that our results at each fixed region size and at every studied coupling were well-characterized by exponential decays.

Interestingly, the relative speed of this decay oscillates as a function of the region size as can be seen in Fig. 7. Regions corresponding to an even number of lattice sites demonstrate a much more sudden initial decay than do those regions comprised of an odd number of sites. This peculiar oscillation results in an inverted shell structure for the extrapolation to , in contrast to the case where in which this feature is preserved. A representative example of this procedure is shown in Fig. 8.

## Vii Summary and Conclusions

We have presented a method to compute the entanglement entropy of interacting fermions which takes advantage of an approximate log-normality property of the distribution of fermion determinants. The resulting approach overcomes the signal-to-noise problem of naive methods, and is very close in its core idea to another method we proposed recently: both methods involve defining an auxiliary parameter , differentiating, and then integrating to recover after a MC calculation. The order of the steps is important, as the differentiation with respect to induces the appearance of entanglement-sensitive contributions in the MC probability measure. Beyond those similarities, the present method has the distinct advantages of being simultaneously simpler to formulate (algebraically as well as computationally) and of explicitly using the approximate log-normality property. Moreover, we have found that the integration step displays clearly more stable numerical behavior in the present approach than in its predecessor: it is approximately linear in the present case and markedly not so in the original incarnation. We therefore strongly advocate using the present algorithm over the former.

In addition to presenting an improved method, we have put forward a straightforward algebraic reformulation of the equations which, while exactly equivalent to the original formalism, avoids the numerical burden of computing inverses of restricted Green’s functions in the calculation of -th order Rényi entropies for . This issue had been pointed out by us and others (see e.g. Ref. Assaad ()) as an inconvenience, as it is perfectly possible for those matrices to be singular.

As a test of our algorithm, we have presented results for the Rényi entropy of the half-filled 1D Hubbard model with periodic boundary conditions. The present and old formalisms were used for calculations at , which matched exactly. The rewritten form based on Eq. (18) was then used to extend our computations to , allowing us to attempt extrapolations in the Rényi order in both directions.

Our results show that, with increasing Rényi order , the value of decreases for all , and the fluctuations as a function of become more pronounced. Remarkably, the statistical MC fluctuations decrease as is increased. Since the problem we set out to solve was in fact statistical in nature, our observations indicate that calculations for large systems and in higher dimensions will benefit from pursuing orders .

###### Acknowledgements.

This material is based upon work supported by the National Science Foundation under Grants No. PHY1306520 (Nuclear Theory program) and No. PHY1452635 (Computational Physics program).## References

- (1) J. E. Drut and W. J. Porter, Phys. Rev. B 92, 125126 (2015).
- (2) P. V. Buividovich, M. I. Polikarpov, Nucl. Phys. B 802, 458 (2008).
- (3) R. G. Melko, A. B. Kallin, and M. B. Hastings, Phys. Rev. B 82, 100409 (2010); M. B. Hastings, I González, A. B. Kallin, and R. G. Melko, Phys. Rev. Lett. 104, 157201 (2010); S. V. Isakov, M. B. Hastings, and R. G. Melko, Nature Phys. 7, 772 (2011); R. R. P. Singh, M. B. Hastings, A. B. Kallin, and R. G. Melko, Phys. Rev. Lett. 106, 135701 (2011); S. Inglis and R. G. Melko, Phys. Rev. E 87, 013306 (2013).
- (4) S. Humeniuk and T. Roscilde, Phys. Rev. B 86, 235116 (2012).
- (5) J. McMinis and N. M. Tubman, Phys. Rev. B 87, 081108(R) (2013).
- (6) P. Broecker and S. Trebst, J. Stat. Mech. (2014) P08015.
- (7) L. Wang, M. Troyer, Phys. Rev. Lett. 113, 110401 (2014).
- (8) D. J. Luitz, X. Plat, N. Laflorencie, and F. Alet, Phys. Rev. B 90, 125105 (2014).
- (9) S. Duane, A. D. Kennedy, B. J. Pendleton, D. Roweth, Phys. Lett. B 195, 216 (1987); S. A. Gottlieb, W. Liu, D. Toussaint, R. L. Renken, Phys. Rev. D 35, 2531 (1987).
- (10) T. Grover, Phys. Rev. Lett. 111, 130402 (2013).
- (11) P. Calabrese, J. L. Cardy, J. Stat. Mech. 0406 (2004) P06002.
- (12) D. Lee, Phys. Rev. C 78, 024001 (2008); D. Lee, Prog. Part. Nucl. Phys. 63, 117 (2009); F. F. Assaad and H. G. Evertz, Worldline and Determinantal Quantum Monte Carlo Methods for Spins, Phonons and Electrons, in Computational Many-Particle Physics, H. Fehske, R. Shnieider, and A. Weise Eds., Springer, Berlin (2008); J. E. Drut and A. N. Nicholson, J. Phys. G: Nucl. Part. Phys. 40, 043101 (2013);
- (13) M.-C. Chung and I. Peschel, Phys. Rev. B 64, 064412 (2001); S.-A. Cheong and C. L. Henley, Phys. Rev. B 69, 075111 (2004); I. Peschel, J. Phys. A 36, L205 (2003).
- (14) F. F. Assaad, T. C. Lang, and F. P. Toldin, Phys. Rev. B 89, 125121 (2014); F. F. Assaad, Phys. Rev. B 91, 125146 (2015).
- (15) M. G. Endres, D. B. Kaplan, J.-W. Lee, A. N. Nicholson, Phys. Rev. Lett. 107, 201601 (2011).
- (16) T. DeGrand, Phys. Rev. D 86, 014512 (2012).