Quantum critical properties of Bose-Hubbard models

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 $\alpha$ describing the singularity of the specific heat could be determined with exceptional accuracy [4], finally resulting in the value $\alpha_{\rm ex}=-0.0127\pm 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 $\phi^{4}$ lattice model and the dynamically diluted $XY$ model, giving $\alpha_{\rm th}=-0.0151\pm 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 $\mu/U$ and the scaled strength $J/U$ of hopping between adjacent lattice sites. Keeping $\mu/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)_{\rm c}$. For almost all $\mu/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}(\mu/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$ – $\mu/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 $\beta$ [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.

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 $\widehat{b}_{i}^{\phantom{\dagger}}$ and $\widehat{b}_{j}^{\dagger}$ 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 $\langle i,j\rangle$, 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 $\eta$ by means of a Legendre transformation [23]. We expand the intensive ground-state energy ${\mathcal{E}}$, that is, the ground-state energy per lattice site of the extended model (1), as a power series in $\eta$,

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 $n$th order, with $n\geq 2k$, therefore allows one to include $n-2k$ tunneling processes. After evaluating the series (5) for given $\mu/U$ to an achievable order $\nu_{\rm 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)_{\rm c}$ and the exponent $\epsilon_{2k}(\mu/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 $\epsilon_{2}$ and $\epsilon_{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.

We first consider the phase boundary $(J/U)_{\rm 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 $(\mu/U)_{\pm}$ limiting a Mott lobe with integer filling factor $g$ can be parametrized by the scaling ansatz [27]

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 $\nu=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 $\nu$.

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 $9$th 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)_{\rm c}$ as a function of $\mu/U$. We then have taken 100 data points computed in a fairly narrow interval of size $\Delta(\mu/U)=0.002$ around the lobe tip, and have fitted them to a function of the form

obtaining

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 $\nu$. Instead, according to our numerical data the phase boundary $(\mu/U)_{\pm}$ 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 $\Delta(\mu/U)=0.002$ around the lobe tip, and taking 100 equidistant fit points, the ansatz (9) yields

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 $\nu=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 the fit

to the ansatz (9). We emphasize that this fit results from the evaluation of merely 5 data points taken from a comparatively wide interval $\Delta(\mu/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.

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 $\epsilon_{2}$, and the divergence exponent $\epsilon_{4}$ multiplied by $2/7$, for chemical potentials covering the lowest four Mott lobes, $0\leq\mu/U\leq 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 $\mu/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].

In contrast, the 2-dimensional models do become critical at the lobe tips. The corresponding comparison of $c_{2}$ and $2/7\cdot 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 $\alpha_{2k}^{(\nu)}(\mu/U)$ constituting the series (5). To this end, we plot in Fig. 7 the ratios $\alpha_{2}^{(\nu)}/\alpha_{2}^{(\nu-1)}$ for $\mu/U=0.3769$ (upper panel) and $\mu/U=0.991$ (lower panel). In the latter case, for $\mu/U$ close to the upper edge of the lowest lobe, the series representation (5) of $c_{2}$ is close to geometric, allowing for a comparatively simple and reliable determination of the limit $\alpha_{2}^{(\nu)}/\alpha_{2}^{(\nu-1)}\to 1/(J/U)_{\rm c}$ [24, 25]. In the former case however, close to the critical lobe tip, the ratios decrease markedly with increasing hopping order $\nu$, and appear to oscillate slightly around a common mean trend, so that their fit to a Gaussian hypergeometric function ${}_{2}F_{1}$, or to generalized hypergeometric functions ${}_{q+1}F_{q}$, indeed is indispensable for correctly estimating the limit $\nu\to\infty$ [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.

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 divergence-exponent inequality (7) previously surmised in Ref. [14]: In general, the divergence exponent $\epsilon_{2}$ which characterizes the one-particle correlation function $c_{2}$ at the phase transition is smaller than $2/7$ times the divergence exponent $\epsilon_{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 $\beta=(\epsilon_{4}-3\epsilon_{2})/2$ for the critical exponent $\beta$ of the order parameter [14]; together with $\epsilon_{4}\stackrel{{\scriptstyle!}}{{=}}7\epsilon_{2}/2$, this yields the remarkable identity

Indeed, reading off $\epsilon_{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 $\beta=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

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.