Hypergeometric continuation of divergent perturbation series. I. Critical exponents of the Bose-Hubbard model

We study the connection between the exponent of the order parameter of the Mott insulator-to-superfluid transition occurring in the two-dimensional Bose-Hubbard model, and the divergence exponents of its one- and two-particle correlation functions. We find that at the multicritical points all divergence exponents are related to each other, allowing us to express the critical exponent in terms of one single divergence exponent. This approach correctly reproduces the critical exponent of the three-dimensional $XY$ universality class. Because divergence exponents can be computed in an efficient manner by hypergeometric analytic continuation, our strategy is applicable to a wide class of systems.


Introduction
Continuous phase transitions are often described by Landau's approach [1][2][3][4][5]: assume that the thermodynamical potentialΓ of a given system possesses the form y y G = + + ( ) a a a , 1 0 2 2 4 4 where the coefficients a 0 , a 2 , a 4 depend on a control parameterj, and the system adopts, for each fixed value of j, that value y min ofψ for which the potential(1) takes on its minimum. If then a 4 is positive and thus guarantees stability, and if one may further neglect the dependence of a 4 onj, while a 2 crosses zero at some value j c , being positive for < j j c and negative for > j j c , one finds y = < ( ) for a two-dimensional square lattice in the J U -m U -plane is shown in figure 1; the corresponding diagrams for triangular or hexagonal lattices are available in the literature [8]. Within the so-called Mott lobes confined at low J U between successive integer valuesg 1 and g of the scaled chemical potential m U the system is in an incompressible Mott state with gparticles per site; when increasing J U at fixed m U it enters the superfluid phase at the phase boundary ( ) J U c . This quantum phase transition has been studied in quantitative detail by quite a number of authors with various methods [9][10][11][12][13][14][15]; it reflects the competition between the lowering of the kinetic energy with increasing delocalization, and the lowering of the repulsion energy with increasing localization of the particles.
It is known from the scaling theory of Fisher et al [7] that generally this Mott insulator-to-superfluid transition is mean field-like in character, but with the exception of the multicritical points with particle-hole symmetry at the tips of the Mott lobes, where the transition takes place at fixed density corresponding to an integer filling factorg. At these special points the transition shown by the d-dimensional Bose-Hubbard model falls into the universality class of the + ( ) d 1 -dimensional XYmodel. Thus, the case d=2 is of particular interest, since it leads to the three-dimensional XY universality class, which also covers the intensely studied lambda-transition undergone by liquid helium [16].
When trying to reconcile this existing knowledge with an approach based on an effective potential(1), one faces several questions: how do the Landau coefficients a 2k , which now also depend, besides the control parameter J U , on the scaled chemical potential m U , manage to switch from 'mean field-like' to 'multicritical' upon variation of m U ? How does one obtain nontrivial critical exponents from this approach, as opposed to the trivial exponent b = 1 2 showing up in equation (5)? What effort would be required to compute these critical exponents along this line with sizeable accuracy? These are the questions we address in the present work, which constitutes a clarification and significant extension of our previous brief communication [17].
Traditionally, the calculation of critical exponents is performed within the framework of the renormalization group (RG) theory [4,5,18,19], having produced fairly precise data. Thus, we do not primarily aim at improving the numerical accuracy of known critical exponents. Instead, we intend to establish a novel bootstrap procedure for computing critical exponents which does not make use of RG theory, and therefore might lead to additional insight. The key input into our analysis are the correlation functions which have been obtained in the accompanying work [20], referred to as paperII in the following. In that more technical paperII we have investigated the analytic continuation of divergent strong-coupling perturbation series for the Bose-Hubbard model by means of generalized hypergeometric functions + F q q 1 , in comparison with the more familiar Shanks transformation and Padé approximation methods, and have found hypergeometric analytic continuation to be particularly well-suited for characterizing the divergence of the correlation functions at the transition points. Nonetheless, the present paper can be read independently from paperII, since here we require only certain results obtained therein, while detailed working knowledge of the hypergeometric continuation technique as such is not necessary. Phase diagram for the Bose-Hubbard model on a two-dimensional square lattice at zero temperature. Within the lobes located at low J U the system is in an incompressible Mott state with gparticles per site, outside these lobes in a superfluid state. The tips of the lobes represent multicritical points; here the system falls into the universality class of the three-dimensional XYmodel. This diagram has been computed according to the hypergeometric scheme developed in [20].
We proceed as follows: in section 2 we briefly recapitulate the formal derivation of the appropriate effective potential for the Bose-Hubbard model [12,13], and state the required relations between the Landau coefficients and the correlation functions. In the central section 3 we then show how to evaluate the critical exponent of the order parameter. In view of the existing accurate reference data, this puts hypergeometric continuation to a truly hard, meaningful test. In section 4 we discuss a property that characterizes the Landau coefficients at the multicritical points. Finally, the discussion led in section 5 concludes our investigation. As an interesting conceptual insight gained from our analysis, we find that it may not always suffice to terminate the effective potential after the fourth-order term, as done in the time-honored paradigm(1); rather, for extracting the order-parameter exponent describing the Mott insulator-to-superfluid transition one also has to resort to the Landau coefficienta 6 . While we do deliberately restrict ourselves to the two-dimensional Bose-Hubbard model for the sake of definiteness, it stands to reason that our methods are also applicable to further systems.

The effective potential for the Bose-Hubbard model
The Bose-Hubbard model is formulated in terms of operators  † b i and  b i which create and annihilate, respectively, a Bose particle at the ith lattice site [6,7]. Thus, they obey the usual commutation relations and the local number operators yield the number of particles occupying the ith site. Employing the pair repulsion energyU as the energy scale of reference, the model Hamiltonian is written in dimensionless form as where the first two sums extend over the entire lattice, and the symbol á ñ i j , is meant to indicate that the third sum ranges over all pairs of neighboring sites i and j. Hence, the first term on the right-hand side gives the total repulsion energy, the second specifies the interaction with the given chemical potential, and the third corresponds to the kinetic energy of the particles. As usual in field theory, we couple this system(8) to external sources and drains which we choose to be spatially uniform with strength η, giving the extended system Without loss of generality we have taken η to be real, since any phase could be removed by an appropriate redefinition of  † b i and  b i . The key quantity of interest for the theoretical description of this model at zero temperature now is the intensive energy landscape , where M denotes the total number of lattice sites, which is assumed to be so large that finite-size effects do not matter, and the expectation value is taken with respect to the ground state of the extended system(9) which, in contrast to the basic system(8), does not conserve the number of particles. Since this ground state energy is an even function of the source strengthη, we expand it in the form represent k-particle correlation functions. In the accompanying paperII we have shown how to evaluate these correlation functions by means of a combination of high-order perturbation theory and hypergeometric analytic continuation. In particular, we have studied the one-particle correlation function c 2 and the two-particle correlation function c 4 : when approaching the phase boundary from within a Mott lobe by varying the control parameter J U at fixed chemical potential m U , these functions diverge as numerically for chemical potentials pertaining to the lowest lobes [20]. Following standard procedures of field theory, the connection between these correlation functions and Landau's description of phase transitions is made by means of a Legendre transformation [4,5,21]: an effective potential Γ is obtained as the Legendre transform of , with respect to the source strengthη [12]. Thus, we introduce a variable ψ according to the Hellmann-Feynman theorem, applied to the extended, particle number non-conserving Hamiltonian (9) then immediately gives the relation The series(10) now yields the representation which, upon inversion, allows one to express the source strengthη in terms of its conjugate variableψ: Performing the Legendre transformation according to the prescription we then obtain the desired effective potential is the intensive ground-state energy of the basic system (8), and the Landau coefficients emerge as certain combinations of the correlation functions. In field-theoretic jargon, these Landau coefficients represent one-particle irreducible (1PI) vertex functions [4,5]. Now we are in a situation analogous to the one considered in the introduction: since η and ψ are Legendreconjugated variables, we have the identity [21] y h ¶G ¶ = -( ) 2 ; 18 since the actual system of interest (8) is recovered by setting h = 0, it corresponds to the stable stationary points y min of Γ. In accordance with equation (13) the Mott insulating phase is characterized by y = 0 min , whereas a nonvanishing value y ¹ 0 min indicates the presence of a superfluid phase, so that y min constitutes a bona fide order parameter of the Mott insulator-to-superfluid transition. However, we are not free to make convenient assumptions concerning the dependence of the Landau coefficients on the control parameter J U and the scaled chemical potential m U , but rather have to respect the above connections between the Landau coefficients and the correlation functions determined in paperII.

Evaluation of the critical exponent
The purpose of this section is to investigate the exponentsb b m = ( ) U which govern the emergence of the order parameter according to when J U is increased at fixed m U beyond the respective transition point ( ) J U c . In particular, we will evaluate the exponent b crit which belongs to the multicritical points at the tips of the Mott lobes shown in figure 1, where we do expect numerical agreement with the critical exponent b = ( ) 0.3485 2 crit characterizing the three-dimensional XY universality class [22].
In paperII the correlation functions m ( ) c U J U , 2 and m ( ) c U J U , 4 have been obtained by fitting their strong-coupling perturbation series in the Mott-insulator regime, that is, for < ( ) J U J U c to hypergeometric functions + F q q 1 , thereby determining their divergence exponents [20]. Here we do not utilize the analytically continued hypergeometric functions for > ( ) J U J U c . Instead, the following analysis relies on the proposition that the asymptotic relations(11), namely on both sides of the pole, thus allowing us to make the decisive step into the superfluid regime.
For the sake of the argument, let us for the moment assume that for certain m U we may neglect terms of order  y ( ) 6 in the full effective potential (17). This assumption reduces the effective potential to the archetypal form(1) reviewed in the introduction. Its minimum y min then is given by According to equation (21) the exponentβ would then be given by We observe that, in contrast to the example reviewed in the introduction, the relation(23) allows a 4 to vanish at the phase boundary. This is indeed what happens: figure 2, which displays the Landau coefficients a 2 and a 4 for the arbitrary value m = U 0.2652 of the scaled chemical potential, shows that both a 2 and a 4 , when considered as a function of J U , become zero at ( ) J U ; c the same behavior is found for all chemical potentials. Therefore, at the phase boundary we are not entitled to neglect terms of order  y ( ) 6 , and have at least to consider the effective potential in the form  This means that c c c 2  In order to deduce the exponent β from equation (27), we distinguish two cases: (i) In case we have a strict inequality    > -2 6 4 2 , the combination c c c Based on the particular shape of the curves drawn in figure 3, we surmise that this actually is a strict inequality aside from the tips, whereas at the tips of the lobes. 39 4 2 Under the assumption of the validity of this equation, the divergence of c 4 2 and that of c 2 7 cancel each other at the tips, and c c 4 2 2 7 has a finite, non-zero limit at the phase boundary. If we further assume that the ratio c c 6 2 6 , and hence a 6 , shares the same behavior, we deduce     = = -6 7 2 . This is precisely the second case(ii) in the above distinction, for which we have derived the equality(36). This leads to a decisive conclusion: inserting the relation(39) between the divergence exponents into this formula(36) for the exponentβ of the order parameter, we obtain the identity for the critical exponent b crit at the tips of the lobes. While the inequality(37) is a general result, this equality(40) hinges on the observations made in figure 3, and applies to the multicritical points only. We thus arrive at an interesting characterization of the multicritical points: the Landau coefficient a 6 diverges for all chemical potentials when  ( ) J U J U c , with the exception of these points. The key result(40), stating that the critical exponent b crit at the tips of the Mott lobes is given by one fourth of the divergence exponent 2 of the two-particle correlation functionc 2 , is amenable to quantitative verification: in table 1 we list the values of b crit as obtained from equation (40) at the tip of the lowest Mott lobe g=1 by hypergeometric continuation based on + F q q 1 withq ranging from 0 to 4; the third column states the relative deviation of the respective result from the reference valueb = ( ) 0.3485 2 crit which has been derived from the f 4 lattice model and the dynamically diluted XY-model by combining Monte Carlo simulations based on finite-size scaling methods, and high-temperature expansions [22]. Evidently, the agreement is quite good.  Speculating further that F 2 1 might yield the most accurate numerical estimates, representing a good compromise between flexibility, as provided by the number of fitting parameters, and the number of available input data still accessible to high-order perturbation theory, we also present estimates for b crit computed with

Characterization of the critical effective potential
The previous observation that the Landau coefficient a 6 diverges for all chemical potentials at the phase boundary, except at the multicritical points, necessitates further investigations. For motivation, let us once again consider the truncated potential(1), which yields the necessary condition  2 . This is tantamount to the observation that all terms in the effective potential(17) display the same asymptotic behavior, and allows us to express the critical exponent b crit in terms of  2 alone, see equation (39).

Discussion
In this paper we have established a connection between the divergence exponents  k 2 of the k-particle correlation functions c 2k pertaining to the two-dimensional Bose-Hubbard model, as defined by equations (10) and (11), and the critical exponent of the order parameter of the Mott insulator-to-superfluid transition. This allows us to take advantage of the fact that the divergence exponents  2 and  4 can be computed numerically with tolerable effort for any value of the scaled chemical potential m U . This is achieved by means of hypergeometric analytic continuation of the strong-coupling perturbation series of c 2 and c 4 , respectively, as demonstrated in detail in paperII [20].
Under the assumption that the effective potential (17) can be truncated after the sixth order term, requiring the Landau coefficient a 6 to be positive, we have derived the lower bound for the exponent b m ( ) U with which the order parameter emerges at the Mott insulator-to-superfluid transition.
For all chemical potentials except those marking the multicritical points the transition is expected to be mean field-like [7]; the bound then is well compatible with the mean-field exponent b = 1 2 mf . At multicriticality, that is, at the lobe tips the bound becomes sharp, and actually yields, to within the numerical accuracy achieved here, the critical exponent b crit of the three-dimensional XY class. Moreover, at multicriticality the divergence exponents are no longer independent of each other, but can all be expressed in terms of  2 . Utilizing the 'multicritical' equality   = · 7 2 Our approach to critical exponents is essentially self-contained, and comparatively straightforward. Along the lines pioneered in this paper, hypergeometric continuation for evaluating divergence exponents may provide critical exponents for wide classes of models, thus shedding further light on the universality hypothesis of statistical physics.