Three-dimensional Hubbard model in the thermodynamic limit

# Three-dimensional Hubbard model in the thermodynamic limit

## Abstract

We employ the numerical linked-cluster expansion to study finite-temperature properties of the uniform cubic lattice Hubbard model in the thermodynamic limit for a wide range of interaction strengths and densities. We carry out the expansion to the 9th order and find that the convergence of the series extends to lower temperatures as the strength of the interaction increases, giving us access to regions of the parameter space that are difficult to reach by most other numerical methods. We study the precise trends in the specific heat, the double occupancy, and magnetic correlations at temperatures as low as 0.2 of the hopping amplitude in the strong-coupling regime. We show that in this regime, accurate estimates for transition temperatures to the Néel ordered phase, in agreement with the predicted asymptotic behavior, can be deduced from the low-temperature magnetic structure factor. We also find evidence for possible instability to the magnetically ordered phase away from, but close to, half filling. Our results have important implications for parametrizing fermionic systems in optical lattice experiments and for benchmarking other numerical methods.

###### pacs:
67.85.-d, 05.30.Jp, 05.30.Fk, 05.70.-a

## I Introduction

The Hubbard model Mott (1949); Hubbard (1963) is the archetypal model for studying strongly-correlated systems and their phase transitions in condensed matter physics. Its Hamiltonian is expressed as

 H=−t∑⟨ij⟩σc†iσcjσ+U∑ini↑ni↓−μ∑iσniσ, (1)

where () annihilates (creates) a fermion with spin on site , is the number operator, is the onsite Coulomb interaction, denotes nearest neighbors, and is the corresponding hopping integral. The two-dimensional (2D) versions have become de facto models to describe Mott transition Imada et al. (1998), potentially high-temperature superconductivity Anderson (1987), and a wealth of other low-temperature phases in materials such as cuprates Dagotto (1994), pnictides Stewart (2011), iron-selenides Dagotto (2013), and heavy fermions Gegenwart et al. (2008), to name a few Scalapino (2012).

Other than a variety of numerical techniques that have been developed over the past three decades to solve this often unforgivingly difficult model Hirsch (1987); Scalettar et al. (1989); Pan and Wang (1997); Kent et al. (2005); Kotliar et al. (2006); Fuchs et al. (2011); Kozik et al. (2013), in recent years, the community has witnessed remarkable efforts in simulating, among other correlated theoretical models, the Hubbard model, using ultracold atoms in optical lattices Jaksch et al. (1998); Bloch et al. (2008). The observation of the Mott transition with fermions as well as bosons in different dimensions Greiner et al. (2002); Stöferle et al. (2004); Spielman et al. (2007); Jordens et al. (2008); Schneider et al. (2008) provided the added impetus for attempts to bring down the temperature to a range relevant to superconductivity. Even though the latter has proven extremely difficult, there has been significant progress along the way.

Some experimental groups have focused on the three-dimensional (3D) system Imriška et al. (2014); Hart et al. (2015), where the Néel transition to the long-range antiferromagnetic (AF) phase is expected to occur at temperatures about an order of magnitude larger than that predicted for superconductivity to develop in the two-dimensional system Staudt et al. (2000); Maier et al. (2005). In a groundbreaking study Hart et al. (2015), long-range AF correlations near the critical temperature of the 3D model were observed for the first time. This milestone was reached through characterization of the experimental system via comparisons of observed thermodynamic properties to results from two state-of-the-art numerical techniques, the determinant quantum Monte Carlo (DQMC) Scalettar et al. (1989) and the numerical linked-cluster expansion (NLCE) Rigol et al. (2006, 2007a, 2007b); Tang et al. (2013a). The latter can provide exact results in the thermodynamic limit, for a wide range of temperatures and average densities.

The 3D Hubbard model has long been a playground for numerical methods to study finite-temperature critical behavior in a system where the electronic correlations can be tuned. Its finite-temperature phase transition to the Néel ordered phase has been studied carefully by a variety of numerical techniques Scalettar et al. (1989); Staudt et al. (2000); Kent et al. (2005); Paiva et al. (2011); Kozik et al. (2013); Hirschmeier et al. (2015). However, as these methods are often not well-equipped to handle the strong-coupling regime of the model (), the complete mapping of the ground state phase diagram in the temperature-interaction plane has been assisted by an asymptotic behavior based on the critical temperature of the low-energy theory in the limit of large interaction strengths Staudt et al. (2000); Hirschmeier et al. (2015).

Here, we use the NLCE to explore the thermodynamic properties, such as the heat capacity, double occupancy, and spin correlations, of the 3D Hubbard model in the thermodynamic limit. An outstanding advantage of the NLCE for the Hubbard models over more typical methods is that one will have access to lower temperatures at larger values of the Coulomb repulsion Khatami and Rigol (2011). This, in turn allows us to explore the divergent behavior of the spin structure factor by using different extrapolation techniques to obtain reliable estimates for the Néel temperature () in the strong-coupling regime. We find that our estimates match the calculated using methods based on quantum Monte Carlo (QMC) within the errorbars in the intermediate-coupling regime (), and are in very good agreement with the asymptotic form for . We further confirm our extrapolation schemes in obtaining the Néel transition temperature for the 3D antiferromagnetic Heisenberg model, for which the NLCE can be carried out to significantly higher orders. The critical temperature is well known for this model from a large-scale QMC study Sandvik (1998).

## Ii the numerical method

To study Hamiltonian (1), we have implemented the NLCE for the 3D cubic lattice. In the NLCE, an extensive property of the lattice model in the thermodynamic limit, normalized to the number of sites, is expressed in terms of contributions from all finite clusters, up to a certain size, that can be embedded in the lattice. These contributions are calculated according to the inclusion-exclusion principle. The method can be summarized in the following series:

 P=1L∑cWP(c) (2)

where represents the extensive property per site in the thermodynamic limit, is the symbolic lattice size (), and the contribution of cluster to the property, also known as the weight, is shown by . If the model does not break translational symmetry of the underlying lattice, the right-hand side of Eq. (2) can be simplified to a sum, without the factor, over only those clusters that are not related by translational symmetry. Further, if the model does not break point group symmetries of the underlying lattice, the contribution of all clusters that are related by point group symmetry can be expressed as one term that is, the weight of one of the clusters, times a multiplicity factor, which represents the number of ways one can obtain a distinct cluster by applying point group symmetry operations to that cluster.

Equation (2) is a cluster expansion Sykes et al. (1966) that can be written not only for the infinite lattice, but also for a finite cluster Tang et al. (2013a). We use this fact to find the weights. Consider, for example, the equation for , the property calculated for a finite cluster :

 p(c)=WP(c)+∑s⊂cWP(s), (3)

where we have intentionally separated the weight of itself, with running over all subclusters of (clusters obtained by removing different number of sites from ). Note that is not normalized by the number of sites. By reordering the terms in this equation, we can write the weight of each cluster as its property less the contributions of its subclusters:

 WP(c)=p(c)−∑s⊂cWP(s). (4)

The above equation provides a recursive method for calculating all the weights up to a certain size; We start with the smallest cluster, a single site, which does not have any subclusters, i.e., , and generate larger clusters by adding sites to it one by one in the so-called site expansion NLCE. We carry out this expansion to the 9th order, which means we will work with clusters as large as 9 sites.

NLCEs use the same basis as the high-temperature series expansions (HTSEs) Oitmaa et al. (2006). However, the calculation of the extensive quantities at the level of individual clusters [] is left to an exact numerical method, such as the exact diagonalization, as opposed to a perturbative expansion in terms of inverse temperature, as done in the HTSEs. Despite the lack of an explicit small parameter, having a finite series inevitably leads to the loss of convergence below a certain temperature, where the correlations in the system extend beyond a length of the order of the largest sizes considered. However, the exact treatment of clusters leads to minimum convergence temperatures that are lower than those of HTSE with a comparable order (see Fig. 1).

We study the specific heat, , which can be obtained within NLCE from the knowledge of the energy, density, and their correlation, without any numerical derivatives or numerical integration Khatami and M. Rigol (2012), double occupancy, (site averaged), the nearest-neighbor spin correlations,

 Szz=1L∑⟨ij⟩⟨SziSzj⟩, (5)

where is the component of the spin operator at site , and the spin structure factor,

 S(q)=1L∑jkeiq⋅(rj−rk)⟨SzjSzk⟩, (6)

at , where is the displacement vector for site on the lattice. We denote the latter by . In the disordered phase at high temperatures, it remains finite in the thermodynamic limit, while its divergence at low temperatures is the indication for AF ordering in the system.

Similar to the analytic Padé approximations used extensively in HTSEs, here we take advantage of two numerical resummation techniques to improve the convergence of our series at low temperatures. We use the Euler algorithm Press et al. (1999) to resum the last 6 terms of the series or the Wynn algorithm Guttmann (1989) with 3 and 4 cycles of improvement (details about these techniques and their use for NLCEs can be found in Refs. Rigol et al., 2007a; Tang et al., 2013a). Except for the staggered structure factor, for which we know the Euler algorithm does not perform as well as the Wynn (see the discussion surrounding Fig. 3), we take the average of properties from the last two orders after the Euler, and the last orders of each of the Wynn resummations as our best estimate. To quantify our confidence in the accuracy of the resummed results, we define a “confidence region” around this average where all four values that contribute to the average fall. Thus, the errorbars in our figures simply mark the boundaries of this region and should not be mistaken as statistical errorbars.

## Iii Results

The NLCE makes the study of thermodynamic quantities in the thermodynamic limit efficient and easy. In Fig. 1, we show the specific heat of the half-filled system as a function of temperature for an interaction strength that ranges from in the weak-coupling regime to in the strong-coupling regime. We have set , and as the unit of energy throughout the paper. Since the properties of the clusters in the series are calculated using exact diagonalization, we obtain information at all temperatures and chemical potential values for a given in a single run, as opposed to QMC-based methods in which each temperature and chemical potential has to be treated separately. For this reason, one can choose a fine temperature and chemical potential grid to better capture the details and trends in quantities of interest in different regions of the parameter space (see Fig. 1). This is also of great importance for modeling of systems in optical lattice experiments. For instance, because of the existence of a trapping potential, theoretical modeling of the inhomogeneous system is often done through proper averaging of properties over a set of homogeneous systems whose chemical potentials vary only slightly from one to the next. This is known as the local density approximation. The other advantage of exact diagonalization is that one has full access to the partition function of the clusters, and so, there is no need to employ numerical integration or derivations, which can introduce systematic errors, for the calculation of , or the entropy.

As can be seen in Fig. 1, we capture the exact location of the high-temperature peak in for all values of shown, and reach temperatures as low as for the largest of 20. The high-temperature peak marks the temperature region where local moments form as the system is cooled. Like for the 2D counterpart Khatami and M. Rigol (2012), and as expected from the increasingly dominant Mott physics, this region moves to higher temperatures as the interaction is increased. On the other hand, the minimum convergence temperature of the series decreases as increases. In the strong-coupling regime, this temperature is proportional to the exchange interaction of the effective Heisenberg model, which scales as . One can see that the lowest convergence temperatures achieved with the NLCE before resummations (thin dotted lines in Fig. 1) are generally lower than those achieved using the HTSE up to the 10th order (red thin solid lines), especially in the intermediate-coupling regime, .

A feature of that can be resolved here with high accuracy is the unique crossing point of curves for different around , which persists for . The physical implications of this phenomenon, which is ubiquitous in strongly-correlated systems and has also been observed to persist for up to the bandwidth in the two-dimensional version of the Hubbard model Paiva et al. (2001); Khatami and M. Rigol (2012) are discussed in Ref. Vollhardt, 1997.

In Figs. 1(b) and 1(c), we show the double occupancy at, and away from, half filling as a function of temperature. At half filling, the uncorrelated limit of at high temperatures is regardless of . However, the larger the , the faster drops upon decreasing the temperature. Quantum fluctuations leave the system with a nonzero and -dependent double occupancy even at . Away from half filling, the uncorrelated values of vary as with the electron density, , and the ground state values are expected to remain nonzero. A clear upturn in upon decreasing the temperature is found after the initial drop at all fillings, at least with at the accessible temperatures. Close to half filling, this phenomenon can be attributed to the tendency of the system to order antiferromagnetically, thus, giving way to a larger number of allowed virtual hoppings to neighboring sites that were otherwise forbidden by the Pauli’s exclusion principle in the uncorrelated system Gorelik et al. (2010); Khatami and Rigol (2011). Sufficiently far from half filling, or in the weak-coupling regime Paiva et al. (2011), it is believed that the upturn is associated with Fermi liquid physics.Werner et al. (2005); Mikelsons et al. (2009); Khatami and Rigol (2011) Our results match those obtained using DQMC with clusters as large as  Paiva et al. (2011); however, in the weak-coupling regime, we are limited to .

The change in the double occupancy as we vary can be better seen in Fig. 2(a), where we show at half filling as a function of in a low-temperature window. The large fluctuations of data at point to the lack of convergence in the series at the lowest temperatures. Regardless, the double occupancy clearly decreases rapidly by increasing with no sharp features or outstanding variation in its behavior as the temperature is decreased from to . In the same time, short-range AF correlations display a non-monotonic behavior. The absolute value of , shown in Fig. 2(b), initially increases with increasing , and then slowly decreases as further increases in the strong-coupling regime. This behavior can be understood by considering the interplay between moment formation, which takes place at higher temperatures for larger values of , and the strength of the exchange interaction between moments, which decreases with increasing . The overall behavior is reminiscent of that for the 2D system Khatami and Rigol (2011). We find that the peak at becomes sharper as the temperature is decreased, hinting that short-range correlations may eventually be strongest at a slightly larger interaction strength than the one corresponding to the largest Néel transition temperature to the long-range order ()  Staudt et al. (2000); Kent et al. (2005); Paiva et al. (2011); Rohringer et al. (2011); Kozik et al. (2013); Hirschmeier et al. (2015).

In Figs. 3 and 4 we explore the AF structure factor, , which contains information about correlations at all length scales. In Fig. 3, we show the temperature dependence of for various interaction strengths. The Euler resummation technique behaves poorly for fast growing properties such as whose terms in the series at a given do not necessarily alternate in sign. For this reason, we perform only Wynn resummations with 3 and 4 cycles of improvement for the structure factor and show the latter at temperatures where the two are within roughly 10% of each other. The reliability of the Wynn resummations for have been confirmed also through previous comparisons to DQMC results Hart et al. (2015). It is worth noting that the lowest convergence temperature decreases by increasing , something we already saw for in Fig. 1. Of course, this does not coincide with a larger at low temperatures for a larger in the strong-coupling regime. In fact, the extent of the correlations in the system, which essentially controls the convergence of the NLCE, is becoming smaller as the interaction is getting stronger.

The model has a finite-temperature phase transition to the long-range Néel ordered state. Therefore, the structure factor for any is expected to diverge at a nonzero temperature. The transition temperature, , is expected to be largest around . So, it is not surprising to find that is also largest for at the lowest accessible temperatures. In the strong-coupling regime, is expected to be around () Sandvik (1998), the strength of the exchange interaction in the effective low-energy Heisenberg model. The inset of Fig. 3 shows the inverse of for , 9, and 16. A simple fit to (constant and ) for suggests that for , which is consistent with current best estimates Paiva et al. (2011); Rohringer et al. (2011); Kozik et al. (2013); Hirschmeier et al. (2015). It is also evident that will be smaller for the other two values of in the inset. The critical temperature as a function of and its different estimates within the NLCE will be discussed below.

Similarly to the nearest-neighbor correlations, the structure factor as a function of , plotted in Fig.  4, exhibits a peak at , which develops faster than that for the former as the temperature is lowered. This is an indication of the fast growing long-range correlations in the system as one approaches the critical temperature. The dashed line in Fig. 4 is a fit proportional to () using the structure factor for the largest three values at . It makes clear the asymptotic behavior of the magnetic correlations in the strong-coupling regime.

So far, we showed results for the structure factor only at half filling. But, what happens to the divergent AF correlations in a system with a ? To answer this question, we plot in Fig. 5 as a function of density for at different temperatures. A very sharp peak develops at as the temperature is lowered, indicating that the correlations in the system remain large only in the close proximity of half filling. These results are of special importance for the simulation of the model using ultracold fermionic atoms in optical lattices as a range of densities are present simultaneously at different radii from the center of the trap Hart et al. (2015).

The system can in principle make a transition to the long-range Néel phase even away from half filling. To explore this possibility, we plot as a function of temperature for different in the inset of Fig. 5. The structure factor at shows strong indication of a nonzero critical temperature that is nonetheless, smaller than that for the half-filled system. At a smaller density of , we do not have enough low-temperature data from the NLCE to draw any conclusion about the critical temperature. We point out that despite the divergent behavior of close to half filing, an instability to a different type of order may be dominant in this region. We have not studied such a scenario here.

As we saw in Fig. 3, the critical temperatures can be estimated from the extrapolations of the structure factor in the intermediate- to strong-coupling regime, where enough information at low temperatures are available. In Fig. 6(a), we plot the Néel temperatures deduced in this way as a function of as filled circles. We also plot for from the dynamical cluster approximation (DCA) Hettler et al. (1998); Kent et al. (2005) and DQMC, which match our results within the errorbars. The NLCE results are in very good agreement with the theoretical prediction for the large- Heisenberg limit for as well Affleck et al. (1988); Sandvik (1998).

In the ordered phase, we expect the maximum to scale linearly with the cluster size, , for finite clusters since the correlations extend to all sites. On the other hand, in the disordered phase above , increases with linearly so long as the linear size of the cluster is smaller in order than the correlation length. For larger systems, as a function of saturates to a temperature-dependent value. Within our NLCE, the order refers to the size of the largest clusters in the expansion. Thus, the order of the expansion does not exactly represent as in a finite-size calculation due to the existence of a large number of smaller clusters in the series. However, the role of the latter is to eliminate boundary effects Sahoo et al. (2016), and so, one could expect that the expansion order would approximately play the role of . In fact, the NLCE order has been successfully used as a length scale to study the scaling of Réyni entropies at quantum critical points Kallin et al. (2013, 2014); Stoudenmire et al. (2014); Sahoo et al. (2016); Devakul and Singh (2015).

With this assumption, we plot in Fig. 6(b) the bare sums (partial sums without resummations) of for vs the NLCE order at five different temperatures, ranging from to . The fluctuations from one order to the other are smaller at higher , where the convergence is achieved faster. For , we expect the function to be weaker than linear, and eventually saturate to a finite value in the thermodynamic limit. On the other hand, for , it is expected to show at least a linear overall increase by increasing the order and to diverge in the thermodynamic limit. In other words, the critical temperature can be located by monitoring the overall curvature of as a function of NLCE order as it changes sign from negative above to positive below . Hence, we fit the data to a second-degree polynomial and locate as the temperature where the quadratic coefficient of the polynomial vanishes. For , the partial sums fluctuate wildly between even and odd orders at temperatures near the expected , and the fits largely underestimate as compared to available estimates. At , we obtain a that is about 10% less than the average of DCA and DQMC estimates. We find that if instead we perform the fits to data only using even orders [shown in Fig. 6(b) as dashed lines], the agreement between the resulting and those from the extrapolations of the structure factor to low temperatures, and the DCA and DQMC estimates, improves in the intermediate-coupling region. The former are plotted in Fig. 6(a) as open circles. Our results are within the errorbars of the results of the DQMC and the DCA, as shown in that figure as triangles and diamonds. The difference between ’s obtained from fits to all orders and fits to only even orders of the NLCE decreases as the interaction increases in the strong-coupling region. For , this difference is only about 3%.

Neglecting odd orders in the fits is an arbitrary choice, which is likely due to the fact that we have a small number of terms in the series. Obviously, we cannot obtain results at higher orders for the Hubbard model. However, we are mostly focused on the strong-coupling regime of this model whose approximate low-energy theory is the Heisenberg model Affleck et al. (1988). The Heisenberg model has the advantage of having a much smaller Hilbert space, which makes it possible for us to go to significantly higher orders in the NLCE. We have carried out the NLCE for the 3D Heisenberg model up to the 14th order. The resulting inverse is plotted in Fig. 6(c), where both the last two orders of the bare sums and the Wynn resummation after 6 cycles of improvement can be seen. The latter matches with the resummation after 5 cycles all the way to the transition temperature for this model, which is found to be . This value is close to that obtained by finite-size scalings in a QMC study (Sandvik (1998). Now, the question is: can we arrive at a similar with fits to the vs the NLCE order? We find that here, with more terms in the expansion, regardless of whether we fit the polynomial using even or all orders of the NLCE, we find . This value is about 6% more than that obtained from QMC. However, this is also roughly the error between the different fitting schemes and between the NLCE and DCA/DQMC estimates for in the Hubbard model when .

In summary, we have implemented a NLCE in three dimensions, up to the 9th order in the site expansion, to study the exact thermodynamic and critical behavior of the Hubbard model in the thermodynamic limit. We study trends in the specific heat, double occupancy and magnetic correlations in the model as we tune the strength of the interaction. We find that both short-range and long-range AF correlations are largest around at the lowest temperatures available. We further extract the Néel temperature by extrapolating the AF structure factor to lower temperatures and find strong evidence that the instability to the AF phase persists at densities close to half filling. We also explore a different scheme in which polynomial fits to bare partial sums of the series for the structure factor can provide accurate estimates of the transition temperature in the strong-coupling region of the model. We confirm this method by extracting for the 3D Heisenberg model, where we can obtain a larger number of terms. This scheme can be exploited in future to study critical phenomena in other models.

We thank Vito W. Scarola for providing the high-temperature series expansion results and are grateful to Rajiv. R. P. Singh, Richard T. Scalettar and Marcos Rigol for their insightful comments on the manuscript. This work used the Extreme Science and Engineering Discovery Environment (XSEDE) under Project No. TG-DMR130143, which is supported by NSF Grant No. ACI-1053575. This work was partly supported by Research Scholarship and Creative Activity grants at San José State University.

### References

1. M. Imada, A. Fujimori,  and Y. Tokura, Rev. Mod. Phys. 70, 1039 (1998).
2. P. W. Anderson, Science 235, 1196 (1987).
3. E. Dagotto, Rev. Mod. Phys. 66, 763 (1994).
4. G. R. Stewart, Rev. Mod. Phys. 83, 1589 (2011).
5. E. Dagotto, Rev. Mod. Phys. 85, 849 (2013).
6. P. Gegenwart, Q. Si,  and F. Steglich, Nat Phys 4, 186 (2008).
7. D. J. Scalapino, Rev. Mod. Phys. 84, 1383 (2012).
8. J. E. Hirsch, Phys. Rev. B 35, 1851 (1987).
9. R. T. Scalettar, D. J. Scalapino, R. L. Sugar,  and D. Toussaint, Phys. Rev. B 39, 4711 (1989).
10. K.-K. Pan and Y.-L. Wang, Phys. Rev. B 55, 2981 (1997).
11. P. R. C. Kent, M. Jarrell, T. A. Maier,  and T. Pruschke, Phys. Rev. B 72, 060411 (2005).
12. G. Kotliar, S. Y. Savrasov, K. Haule, V. S. Oudovenko, O. Parcollet,  and C. A. Marianetti, Rev. Mod. Phys. 78, 865 (2006).
13. S. Fuchs, E. Gull, L. Pollet, E. Burovski, E. Kozik, T. Pruschke,  and M. Troyer, Phys. Rev. Lett. 106, 030401 (2011).
14. E. Kozik, E. Burovski, V. W. Scarola,  and M. Troyer, Phys. Rev. B 87, 205102 (2013).
15. D. Jaksch, C. Bruder, J. I. Cirac, C. W. Gardiner,  and P. Zoller, Phys. Rev. Lett. 81, 3108 (1998).
16. I. Bloch, J. Dalibard,  and W. Zwerger, Rev. Mod. Phys. 80, 885 (2008).
17. M. Greiner, O. Mandel, T. Esslinger, T. Hansch,  and I. Bloch, Nature 415, 39 (2002).
18. T. Stöferle, H. Moritz, C. Schori, M. Köhl,  and T. Esslinger, Phys. Rev. Lett. 92, 130403 (2004).
19. I. B. Spielman, W. D. Phillips,  and J. V. Porto, Phys. Rev. Lett. 98, 080404 (2007).
20. R. Jordens, N. Strohmaier, K. Gunter, H. Moritz,  and T. Esslinger, Nature 455, 204 (2008).
21. U. Schneider, L. HackermÃ¼ller, S. Will, T. Best, I. Bloch, T. A. Costi, R. W. Helmes, D. Rasch,  and A. Rosch, Science 322, 1520 (2008).
22. J. Imriška, M. Iazzi, L. Wang, E. Gull, D. Greif, T. Uehlinger, G. Jotzu, L. Tarruell, T. Esslinger,  and M. Troyer, Phys. Rev. Lett. 112, 115301 (2014).
23. R. A. Hart, P. M. Duarte, T. L. Yang, X. Liu, T. Paiva, E. Khatami, R. T. Scalettar, N. Trivedi, D. A. Huse,  and R. G. Hulet, Nature 519, 211 (2015).
24. R. Staudt, M. Dzierzawa,  and A. Muramatsu, The European Physical Journal B - Condensed Matter and Complex Systems 17, 411 (2000).
25. T. A. Maier, M. Jarrell, T. C. Schulthess, P. R. C. Kent,  and J. B. White, Phys. Rev. Lett. 95, 237001 (2005).
26. M. Rigol, T. Bryant,  and R. R. P. Singh, Phys. Rev. Lett. 97, 187202 (2006).
27. M. Rigol, T. Bryant,  and R. R. P. Singh, Phys. Rev. E 75, 061118 (2007a).
28. M. Rigol, T. Bryant,  and R. R. P. Singh, Phys. Rev. E 75, 061119 (2007b).
29. B. Tang, E. Khatami,  and M. Rigol, Computer Physics Communications 184, 557 (2013a).
30. T. Paiva, Y. L. Loh, M. Randeria, R. T. Scalettar,  and N. Trivedi, Phys. Rev. Lett. 107, 086401 (2011).
31. D. Hirschmeier, H. Hafermann, E. Gull, A. I. Lichtenstein,  and A. E. Antipov, Phys. Rev. B 92, 144409 (2015).
32. E. Khatami and M. Rigol, Phys. Rev. A 84, 053611 (2011).
33. A. W. Sandvik, Phys. Rev. Lett. 80, 5196 (1998).
34. M. F. Sykes, J. W. Essam, B. R. Heap,  and B. J. Hiley, Journal of Mathematical Physics 7, 1557 (1966).
35. J. Oitmaa, C. Hamer,  and W. Zhang, Series Expansion Methods for Strongly Interacting Lattice Models (Cambridge University Press, Cambridge, England, 2006).
36. E. Khatami and M. Rigol, Phys. Rev. A 86, 023633 (2012).
37. H. W. Press, B. P. Flannery, S. A. Teukolsky,  and W. T. Vetterling, Numerical Recipes in Fortran (Cambridge University Press, Cambridge, England, 1999) Chap. 5.1.
38. A. J. Guttmann, Phase Transitions and Critical Phenomena, edited by C. Domb and J. Lebowitz, Vol. 13 (Academic Press, London, 1989).
39. V. W. Scarola, L. Pollet, J. Oitmaa,  and M. Troyer, Phys. Rev. Lett. 102, 135302 (2009).
40. R. Jördens, L. Tarruell, D. Greif, T. Uehlinger, N. Strohmaier, H. Moritz, T. Esslinger, L. De Leo, C. Kollath, A. Georges, V. Scarola, L. Pollet, E. Burovski, E. Kozik,  and M. Troyer, Phys. Rev. Lett. 104, 180401 (2010).
41. L. De Leo, J.-S. Bernier, C. Kollath, A. Georges,  and V. W. Scarola, Phys. Rev. A 83, 023606 (2011).
42. T. Paiva, R. T. Scalettar, C. Huscroft,  and A. K. McMahan, Phys. Rev. B 63, 125116 (2001).
43. D. Vollhardt, Phys. Rev. Lett. 78, 1307 (1997).
44. E. V. Gorelik, I. Titvinidze, W. Hofstetter, M. Snoek,  and N. Blümer, Phys. Rev. Lett. 105, 065301 (2010).
45. F. Werner, O. Parcollet, A. Georges,  and S. R. Hassan, Phys. Rev. Lett. 95, 056401 (2005).
46. K. Mikelsons, E. Khatami, D. Galanakis, A. Macridin, J. Moreno,  and M. Jarrell, Phys. Rev. B 80, 140505 (2009).
47. M. H. Hettler, A. N. Tahvildar-Zadeh, M. Jarrell, T. Pruschke,  and H. R. Krishnamurthy, Phys. Rev. B 58, R7475 (1998).
48. G. Rohringer, A. Toschi, A. Katanin,  and K. Held, Phys. Rev. Lett. 107, 256402 (2011).
49. I. Affleck, Z. Zou, T. Hsu,  and P. W. Anderson, Phys. Rev. B 38, 745 (1988).
50. S. Sahoo, E. M. Stoudenmire, J.-M. Stéphan, T. Devakul, R. R. P. Singh,  and R. G. Melko, Phys. Rev. B 93, 085120 (2016).
51. A. B. Kallin, K. Hyatt, R. R. P. Singh,  and R. G. Melko, Phys. Rev. Lett. 110, 135702 (2013).
52. 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).
53. E. M. Stoudenmire, P. Gustainis, R. Johal, S. Wessel,  and R. G. Melko, Phys. Rev. B 90, 235106 (2014).
54. T. Devakul and R. R. P. Singh, Phys. Rev. Lett. 115, 187201 (2015).
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 minumum 40 characters