Quantum critical properties of Bose-Hubbard models

The Mott insulator-to-superfluid transition exhibited by the Bose-Hubbard model on a two-dimensional square lattice occurs for any value of the chemical potential, but becomes critical at the tips of the so-called Mott lobes only. Employing a numerical approach based on a combination of high-order perturbation theory and hypergeometric analytic continuation we investigate how quantum critical properties manifest themselves in computational practice. We consider two-dimensional triangular lattices and three-dimensional cubic lattices for comparison, providing accurate parametrizations of the phase boundaries at the tips of the respective first lobes. In particular, we lend strong support to a recently suggested inequality which bounds the divergence exponent of the one-particle correlation function in terms of that of the two-particle correlation function, and which sharpens to an equality if and only if a system becomes critical.


Introduction
The application of the theory of critical phenomena [1][2][3] to the description of the lambda transition undergone by liquid helium at about 2.18 K is still confronting physicists with a seemingly minor, yet annoying challenge: Owing to measurements performed under zero-gravity conditions in Earth orbit in order to reduce the rounding of the transition caused by gravitationally induced pressure gradients, the critical exponent α describing the singularity of the specific heat could be determined with exceptional accuracy [4], finally resulting in the value α ex = −0.0127 ± 0.0003. On the other hand, theoretical calculations so far have not managed to reproduce this figure to an accuracy compatible with the experimental uncertainty, despite substantial effort [5][6][7][8][9]. To our knowledge, the currently most accurate theoretical value has been reported by Campostrini et al. [8], based on finite-size scaling analyses of high-statistics Monte Carlo simulations and resummations of 22nd-order high-temperature expansions for the φ 4 lattice model and the dynamically diluted XY model, giving α th = −0.0151 ± 0.0003. Since this comparison of measurement with calculation constitutes a core test of the renormalization group theory, the remaining discrepancy needs to be taken seriously.
It is therefore of interest to observe that the quantum phase transition exhibited at zero temperature by the Bose-Hubbard model [10] may contribute to solving this long-standing puzzle. In a grand-canonical setting this model depends on two dimensionless parameters, the scaled chemical potential µ/U and the scaled strength J/U of hopping between adjacent lattice sites. Keeping µ/U fixed at a noninteger positive value, and increasing J/U from zero to positive values, the d-dimensional model features a transition from an incompressible Mott insulator to a superfluid at a certain (J/U ) c . For almost all µ/U that transition is mean field-like, but at particular multicritical points corresponding to the tips of the so-called Mott lobes the transition falls into the universality class of the (d + 1)-dimensional XY model [10]. In particular, the comparatively simple 2-dimensional Bose-Hubbard model then belongs to the 3-dimensional XY universality class, the same class which also covers the lambda transition of liquid helium. Thus, in principle it should be possible to deduce the critical exponents of the lambda transition from a 2-dimensional Bose-Hubbard system [11].
We have previously outlined a strategy to exploit this connection: Utilizing numerically executed strong-coupling perturbation theory to high orders, in the guise of Eckardt's process-chain approach [12], one obtains divergent-series representations of the model's k-particle correlation functions c 2k (µ/U, J/U ); these divergent series can then be analytically continued with the help of generalized hypergeometric functions, taking up a recent suggestion by Mera et al. [13] The locus in the J/U -µ/U parameter plane where the correlation functions diverge then marks the phase boundary; the manner how they diverge at the multicritical points provides information on the critical exponents [14,15]. Applied to the 2-dimensional Bose-Hubbard model on a square lattice, the first hard test of this strategy indeed has yielded a fairly good estimate of the critical exponent β [14]. Thus, the combination of high-order perturbation theory with hypergeometric analytic continuation provides an alternative, self-contained "bootstrap technique" for computing critical exponents that is independent from other approaches; it is hoped that further refinements of this technique will contribute to resolving the current discrepancy between the experimental and theoretical benchmark data [4,8].
Here we report further results of our approach not only for the 2-dimensional square lattice, but also for the triangular lattice [16]. As an essential consistency check we also consider the Bose-Hubbard model on a 3-dimensional cubic lattice, which falls into the 4-dimensional XY class and therefore should not become critical at all.

Technical background
We start by summarizing the technical details of our procedure: Employing the repulsion energy U of a pair of bosons occupying a common lattice site as reference energy, we write the Bose-Hubbard Hamiltonian in the dimensionless form where the indices i, j label individual sites, and the operators b i and b † j are bosonic annihilation and creation operators, obeying the canonical commutation relation Moreover, denotes the number operator which counts the particles on site i. Thus, the first two terms of the Hamiltonian (1) specify, respectively, the total scaled on-site repulsion energy, and the system's coupling to the chemical potential. The sum in the third term ranges over all pairs of sites i, j that are connected by a tunneling contact, as indicated by the symbol i, j , describing the hopping of Bosons from sites j to i. Taken together, these three terms constitute the conventional Bose-Hubbard model that has been studied in great detail by many authors, using a wide variety of techniques [17][18][19][20][21][22]. The fourth term represents spatially homogeneous sources and drains which explicitly break particle number conservation; this allows one to introduce an effective potential depending on the source strength η by means of a Legendre transformation [23]. We expand the intensive ground-state energy E, that is, the ground-state energy per lattice site of the extended model (1), as a power series in η, and then apply the process-chain approach [12] for estimating the k-particle correlation functions To this end, we start from the ground state of the site-diagonal first two terms of the system (1), and include the other two terms as perturbations [24][25][26]. The computation of c 2k then prompts one to account for k creation and k annihilation processes; invoking perturbation theory to nth order, with n ≥ 2k, therefore allows one to include n − 2k tunneling processes. After evaluating the series (5) for given µ/U to an achievable order ν max in the hopping strength J/U , we fit this approximant by hypergeometric functions q+1 F q , allowing us to determine both the respective transition point (J/U ) c and the exponent 2k (µ/U ) which governs the divergence of c 2k at that point [15], From general arguments invoking neither the dimension d nor the geometry of the underlying lattice we have inferred the inequality which is expected to reduce to an exact equality if and only if the system becomes critical [14]. Thus, a comparison of the divergence exponents 2 and 4 , computed independently for a given variant of the model, enables one to detect quantum criticality. One of the main purposes of the present paper is to lend further support to this relation (7), which indicates quantum criticality through equality of both sides, by reporting the results of large-scale numerical calculations for three variants of the Bose-Hubbard model.

Accurate parametrization of phase boundaries
We first consider the phase boundary (J/U ) c which, in accordance with the asymptotic relation (6), can already be deduced from c 2 alone. This study is motivated by a suggestion due to Freericks et al. [27], claiming that the two values (µ/U ) ± limiting a Mott lobe with integer filling factor g can be parametrized by the scaling ansatz [27] (µ/U ) where x = dJ/U , and A(x), B(x), and the scaling polynomial are certain polynomials which have been stated for a 2-dimensional square lattice and g = 1 as Eq. (104) in Ref. [27], and ν = 0.6717(1) is the critical exponent which governs the correlation length [8]. This is an interesting suggestion, implying that accurate knowledge of the phase boundary suffices to determine ν.
We therefore have made a precise computation of the phase boundary for the first Mott lobe of a 2-dimensional Bose-Hubbard model on a square lattice by combining process-chain calculations to 9th order in J/U with hypergeometric analytic continuation based on 2 F 1 [14,15], and have plotted the very tip of this lobe in Fig. 1; here we consider (J/U ) c as a function of µ/U . We then have taken 100 data points computed in a fairly narrow interval of size ∆(µ/U ) = 0.002 around the lobe tip, and have fitted them to a function of the form The deviation of the exponent y = 2.00018 "measured" here from the exponent 2 of a perfect parabola is fully accounted for by the limits to our numerical accuracy; clearly, the quality of our fit is not compatible with the suggested appearance of the critical exponent ν. Instead, according to our numerical data the phase boundary (µ/U ) ± right at the tip appears to be given by a smooth square root function.
We have also performed such calculations for a triangular lattice: While each site of a square lattice is linked by a tunneling contact to four nearest neighbors, each site is surrounded by even six nearest neighbors on a triangular one. The results are shown in Fig. 2. As expected on the grounds of the respective coordination number, in comparison with the square lattice the lobe tip shifts to lower J/U in the triangular case [16]. Again selecting an interval of width ∆(µ/U ) = 0.002 around the lobe tip, and taking 100 equidistant fit points, the ansatz (9) yields (J/U ) c = 0.03779 − 0.2926 · |µ/U − 0.3834| 2.00003 , likewise suggesting that the lobe tip is well approximated by an exact parabola, without correction involving a critical exponent. For comparison, we also have considered a 3-dimensional cubic lattice, for which the Mott transition does not become critical [10]. From the perspective of perturbation theory, the process chains defining c 2 for the 2-dimensional square lattice and those for the cubic lattice are virtually the same up to and including the 7th order in J/U , with only different weight factors accounting for the different dimensions. Starting with the order ν = 8 additional diagrams begin to appear in three dimensions, such as depicted exemplarily in Fig. 3. The tip of the first Mott lobe is displayed in Fig. 4, together with c 2 = · · · + 1 3! · (2d)! (2d−5)! · + · · · Figure 3. Diagram contributing to the 8th order in J/U of the perturbation series for the one-particle correlation function c 2 of the d-dimensional Bose-Hubbard model on a (hyper-)cubic lattice with d ≥ 3. Arrows denote hopping processes between neighboring lattice sites, the circle symbolizes a creation process, and the cross an annihilation process. The dotted arrow indicates an out-of-plane hopping process which distinguishes the 3-dimensional model from the 2-dimensional one on a square lattice.  (12) to the ansatz (9). We emphasize that this fit results from the evaluation of merely 5 data points taken from a comparatively wide interval ∆(µ/U ) = 0.0319, more than 15 times as large as employed in Fig. 1. Yet, the deviation of the observed exponent y = 2.00559 from the parabola exponent 2 is not significant.

Manifestation of quantum criticality
Since it does not appear to be feasible to recognize criticality by inspection of the phase boundary alone, we now turn to the more demanding criterion (7), prompting us to evaluate c 4 besides c 2 . We start with the 3-dimensional cubic lattice: Figure 5 shows both the divergence exponent 2 , and the divergence exponent 4 multiplied by 2/7, for chemical potentials covering the lowest four Mott lobes, 0 ≤ µ/U ≤ 4. To produce this figure, the one-particle correlation function c 2 has been computed to 9th order in J/U , and the two-particle correlation function c 4 to 7th order, amounting to the evaluation of perturbation theory to 11th order in both cases. Evidently the inequality (7) is well satisfied, equality being excluded for all µ/U . This is a fairly important consistency check, vindicating the recognition that the 3-dimensional Bose-Hubbard model does not become critical even at the tips of its Mott lobes [10,11].  Figure 5. Comparison of the divergence exponent 2 of the one-particle correlation function c 2 (full line) to 2/7 times the divergence exponent 4 of the two-particle correlation function c 4 (dotted), for the 3-dimensional Bose-Hubbard model on a cubic lattice. Here c 2 has been evaluated to 9th order in J/U , and c 4 to 7th order. The inequality (7) is well satisfied for all values of the scaled chemical potential µ/U , confirming that this model does not become critical. Here c 2 has been evaluated to 10th order in J/U , and c 4 to 7th order. To numerical accuracy the inequality (7) actually becomes an equality at the tips of the four Mott lobes considered here, indicating quantum criticality at these points.  In contrast, the 2-dimensional models do become critical at the lobe tips. The corresponding comparison of c 2 and 2/7 · c 4 for the square lattice is shown in Fig. 6. In marked contrast to Fig. 5, here the two lines touch each other at the lobe tips to within numerical accuracy [14], so that the inequality (7) reduces to an equality at these points, thereby heralding criticality.
For gaining insight into the working principles of hypergeometric analytic continuation it is instructive to inspect the behavior of the coefficients α (ν) 2k (µ/U ) constituting the series (5). To this end, we plot in Fig. 7 the ratios α  [24,25]. In the former case however, close to the critical lobe tip, the ratios decrease markedly with increasing hopping order ν, and appear to oscillate slightly around a common mean trend, so that their fit to a Gaussian hypergeometric Here c 2 has been evaluated to 9th order in J/U , and c 4 to 6th order. To numerical accuracy the inequality (7) actually becomes an equality at the tips of the four Mott lobes considered here, indicating quantum criticality at these points. function 2 F 1 , or to generalized hypergeometric functions q+1 F q , indeed is indispensable for correctly estimating the limit ν → ∞ [15]. It is now of major interest to check the crucial inequality (7) for the triangular lattice. The numerical data shown in Fig. 8, where once again we have plotted both sides of this inequality, are fairly convincing indeed: To within numerical accuracy, the inequality (7) is confirmed, with exact equality being reached, again to within numerical accuracy, at the lobe tips. Thus, an approach based on high-order perturbation theory combined with hypergeometric analytic continuation definitively is sensitive to quantum criticality.

Discussion
The most important findings of the present numerical study are encoded in Figs. 5, 6, and 8, representing a computationally rather demanding verification of the divergenceexponent inequality (7) previously surmised in Ref. [14]: In general, the divergence exponent 2 which characterizes the one-particle correlation function c 2 at the phase transition is smaller than 2/7 times the divergence exponent 4 of c 4 , but both sides become equal if and only if the system becomes critical.
Yet, there is one more issue that is at stake here. At criticality one expects the relation β = ( 4 −3 2 )/2 for the critical exponent β of the order parameter [14]; together with 4 ! = 7 2 /2, this yields the remarkable identity Indeed, reading off 2 from the lobe tips seen in Fig. 6 for the square lattice with its coordination number Z = 4, and dividing by four, one obtains fairly good agreement with the expected value β = 0.3486(1) reported by Campostrini et al. [8]. While we do no yet reach the numerical accuracy level of this benchmark value, this may become feasible after honing our method still further, possibly allowing for a rather stringent test of the universality hypothesis of statistical physics. But the same procedure applied to the data depicted in Fig. 8 for the triangular lattice with coordination number Z = 6 yields the markedly different value β = 0.332(1) .
This is an interesting observation, suggesting that critical exponents of the Mott insulator-to-superfluid quantum phase transition in the Bose-Hubbard model depend not only on the dimensionality, but also on the coordination number of the respective lattice. The experimental verification of this prediction, possibly with ultracold atoms in optical lattices [28], could provide important guidance for the further development of the theory.